typo or stale info in qr man

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

typo or stale info in qr man

Wojciech Musial (Voitek)
man for `qr` says that the function uses LINPACK's DQRDC, while it in
fact uses DQRDC2.

```
The QR decomposition of the matrix as computed by LINPACK or LAPACK.
The components in the returned value correspond directly to the values
returned by DQRDC/DGEQP3/ZGEQP3
```

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: typo or stale info in qr man

Martin Maechler
>>>>> Wojciech Musial (Voitek) <[hidden email]>
>>>>>     on Mon, 24 Oct 2016 15:07:55 -0700 writes:

    > man for `qr` says that the function uses LINPACK's DQRDC, while it in
    > fact uses DQRDC2.

which is a modification of LINPACK's DQRDC.

But you are right, and I have added to the help file (and a tiny
bit to the comments in the Fortran source).

When this change was done > 20 years ago, it was still hoped
that the numerical linear algebra community or more specifically
those behind LAPACK would eventually provide this functionality
with LAPACK (and we would then use that),
but that has never happened according to my knowledge.

Thank you for the 'heads up'.

Martin Maechler
ETH Zurich

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: typo or stale info in qr man

Jari Oksanen
And that missing functionality is that Linpack/Lapack routines do not return rank and have a different style of pivoting? For other aspects, the user-interface is very similar in dqrdc2 in R and in dqrdc in Linpack. Another difference seems to be that the final pivoting reported to the user is different: R keeps the original order except for aliased variables, but Linpack makes either wild shuffling or no pivoting at all. I haven't looked at dqpq3 in Lapack, but it appears to return no rank either (don't know about shuffling the columns). It seems that using Linpack dqrdc directly is not always compatible with dqrdc2 of R although it returns similar objects. That is, when packing up the Linpack function to produce an object with same items as qr.default (qr, rank, qraux, pivot, class "qr"), the result object may not yield similar results in base::qr.fitted, base::qr.resid etc as base::qr.default result (but I haven't had time for thorough testing).

This is how I tried to do the packing (apologies for clumsy coding):

SEXP do_QR(SEXP x, SEXP dopivot)
{
    /* set up */
    int i;
    int nr = nrows(x), nx = ncols(x);
    int pivoting = asInteger(dopivot);
    SEXP qraux = PROTECT(allocVector(REALSXP, nx));
    SEXP pivot = PROTECT(allocVector(INTSXP, nx));
    /* do pivoting or keep the order of columns? */
    if (pivoting)
        memset(INTEGER(pivot), 0, nx * sizeof(int));
    else
        for(i = 0; i < nx; i++)
            INTEGER(pivot)[i] = i+1;
    double *work = (double *) R_alloc(nx, sizeof(double));
    int job = 1;
    x = PROTECT(duplicate(x));

    /* QR decomposition with Linpack */
    F77_CALL(dqrdc)(REAL(x), &nr, &nr, &nx, REAL(qraux),
                    INTEGER(pivot), work, &job);

    /* pack up */
    SEXP qr = PROTECT(allocVector(VECSXP, 4));
    SEXP labs = PROTECT(allocVector(STRSXP, 4));
    SET_STRING_ELT(labs, 0, mkChar("qr"));
    SET_STRING_ELT(labs, 1, mkChar("rank"));
    SET_STRING_ELT(labs, 2, mkChar("qraux"));
    SET_STRING_ELT(labs, 3, mkChar("pivot"));
    setAttrib(qr, R_NamesSymbol, labs);
    SEXP cl = PROTECT(allocVector(STRSXP, 1));
    SET_STRING_ELT(cl, 0, mkChar("qr"));
    classgets(qr, cl);
    UNPROTECT(2); /* cl, labs */
    SET_VECTOR_ELT(qr, 0, x);
    SET_VECTOR_ELT(qr, 1, ScalarInteger(nx)); /* not really the rank,
                                                 but no. of columns */
    SET_VECTOR_ELT(qr, 2, qraux);
    SET_VECTOR_ELT(qr, 3, pivot);
    UNPROTECT(4); /* qr, x, pivot, qraux */
    return qr;
}


cheers, Jari Oksanen
________________________________________
From: R-devel <[hidden email]> on behalf of Martin Maechler <[hidden email]>
Sent: 25 October 2016 11:08
To: Wojciech Musial (Voitek)
Cc: [hidden email]
Subject: Re: [Rd] typo or stale info in qr man

>>>>> Wojciech Musial (Voitek) <[hidden email]>
>>>>>     on Mon, 24 Oct 2016 15:07:55 -0700 writes:

    > man for `qr` says that the function uses LINPACK's DQRDC, while it in
    > fact uses DQRDC2.

which is a modification of LINPACK's DQRDC.

But you are right, and I have added to the help file (and a tiny
bit to the comments in the Fortran source).

When this change was done > 20 years ago, it was still hoped
that the numerical linear algebra community or more specifically
those behind LAPACK would eventually provide this functionality
with LAPACK (and we would then use that),
but that has never happened according to my knowledge.

Thank you for the 'heads up'.

Martin Maechler
ETH Zurich

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: typo or stale info in qr man

Peter Dalgaard-2
In reply to this post by Martin Maechler

On 25 Oct 2016, at 10:08 , Martin Maechler <[hidden email]> wrote:

>>>>>> Wojciech Musial (Voitek) <[hidden email]>
>>>>>>    on Mon, 24 Oct 2016 15:07:55 -0700 writes:
>
>> man for `qr` says that the function uses LINPACK's DQRDC, while it in
>> fact uses DQRDC2.
>
> which is a modification of LINPACK's DQRDC.
>
> But you are right, and I have added to the help file (and a tiny
> bit to the comments in the Fortran source).
>
> When this change was done > 20 years ago, it was still hoped
> that the numerical linear algebra community or more specifically
> those behind LAPACK would eventually provide this functionality
> with LAPACK (and we would then use that),
> but that has never happened according to my knowledge.
>

I had some thoughts on this recently and resolved that the base issue is that R wants successive (Gram/Schmidt-type) orthogonalization of the design matrix, not really QR as such.

The LINPACK QR routine happens to work by orthogonalization, but it is far from the only way of doing QR, and most likely not the "best" one (speedwise/precisionwise) if a QR decompositiion as such is the target. (Pivoting is only part of the story)

lm() and associates (notably anova()) relies so much on successive terms being orthogonalized that method="qr" really is a misnomer. For much the same reason, it really is too much to expect that numerical analysts would enforce orthogonality features on a general QR-decomposer.

I suppose that if we want to be free of LINPACK, we may need to step back and write our own orthogonalization routines based on other routines in LAPACK or on the BLAS directly.

-pd

> Thank you for the 'heads up'.
>
> Martin Maechler
> ETH Zurich
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel