From david.clayton at cimr.cam.ac.uk  Thu Oct 27 16:50:51 2005
From: david.clayton at cimr.cam.ac.uk (david.clayton@cimr.cam.ac.uk)
Date: Thu, 27 Oct 2005 16:50:51 +0200 (CEST)
Subject: [Rd] Column names in qr() and chol() (PR#8258)
Message-ID: <20051027145051.7B03820BF3@slim.kubism.ku.dk>

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



[Rd] Column names in qr() and chol() (PR#8258)
david.clayton@cimr.cam.ac.uk