[Rd] Column names in qr() and chol() (PR#8258)
david.clayton@cimr.cam.ac.uk
david.clayton at cimr.cam.ac.uk
Thu Oct 27 16:50:51 CEST 2005
I am using 2.2.0
If the QR decomposition of an N*M matrix is such that the pivoting order
is not 1:M, Q%*%R does not result in the original matrix but in a
matrix with the columns permuted. This is clearly intentional, and
probably to be expected if pivoting is used --- chol() behaves in the
same manner (it would perhaps be nice if the qr help page made that
clear in the same way that the chol() help page does).
The small bug is that column names are not permuted in the same way, so
that, after permuting the columns of Q%*%R so that it is equal to the
input matrix, it is mislabelled. The same small bug is present in
chol() when pivoting is switched on.
It's a minor matter but, perhaps, easily fixed. I think one line which
reordered the column names of the $qr list element of the object
returned by qr() would fix things (and similarly for chol)
Example:
> A <- 1:10
> B <- A^2
> C <- B-A
> D <- rep(1, 10)
> X <- cbind(A,B,C,D)
> qrX <- qr(X)
> qrX$pivot
[1] 1 2 4 3
> oo <- order(qrX$pivot)
> Q <- qr.Q(qrX)
> R <- qr.R(qrX)
> Q%*%R
A B C D
[1,] 1 1 1 2.875732e-14
[2,] 2 4 1 2.000000e+00
[3,] 3 9 1 6.000000e+00
[4,] 4 16 1 1.200000e+01
[5,] 5 25 1 2.000000e+01
[6,] 6 36 1 3.000000e+01
[7,] 7 49 1 4.200000e+01
[8,] 8 64 1 5.600000e+01
[9,] 9 81 1 7.200000e+01
[10,] 10 100 1 9.000000e+01
> (Q%*%R)[,oo]
A B D C
[1,] 1 1 2.875732e-14 1
[2,] 2 4 2.000000e+00 1
[3,] 3 9 6.000000e+00 1
[4,] 4 16 1.200000e+01 1
[5,] 5 25 2.000000e+01 1
[6,] 6 36 3.000000e+01 1
[7,] 7 49 4.200000e+01 1
[8,] 8 64 5.600000e+01 1
[9,] 9 81 7.200000e+01 1
[10,] 10 100 9.000000e+01 1
# Note: column names of last two columns have been interchanged
> XtX <- t(X)%*%X
> U <- chol(XtX, pivot=T)
> oo <- order(attr(U, "pivot"))
> XtX
A B C D
A 385 3025 2640 55
B 3025 25333 22308 385
C 2640 22308 19668 330
D 55 385 330 10
> t(U[,oo])%*% U[,oo]
D A B C
D 385 3025 2640 55
A 3025 25333 22308 385
B 2640 22308 19668 330
C 55 385 330 10
# Matrix has been successfully reproduced but mislabelled
David Clayton
More information about the R-devel
mailing list