0

Just puzzled by that : I'm computing the QR decomposition of the following matrix X and I expected to get back to X by calculating the product QR. Well, the product gives "almost" X, but with columns 1 and 2 inverted... How can it be ?

> X
     [,1]  [,2] [,3]
[1,] 1-0i -1+1i 1-1i
[2,] 1-1i -1+2i 1-1i
[3,] 1-1i -1+1i 1+0i
> 
> QRX <- qr(X)
> X <- qr.X(QRX)
> Q <- qr.Q(QRX)
> R <- qr.R(QRX)
> 
> Q%*%R
      [,1] [,2] [,3]
[1,] -1+1i 1-0i 1-1i
[2,] -1+2i 1-1i 1-1i
[3,] -1+1i 1-1i 1+0i
Andrew
  • 926
  • 2
  • 17
  • 24
  • See `?qr.R` which explains that the output 'may be pivoted'. – Andrew Gustar Nov 28 '17 at 11:27
  • You can use `qr(X)$pivot` to access what pivot. So to get back your original try something like `library(binhf);(qr.Q(qr(X)) %*%qr.R(qr(X)))[,shift(qr(X)$pivot,1)]` – Michael Bird Nov 28 '17 at 11:35
  • OK, thanks. But having to use `Q%*%R[, QRX$pivot]` instead of simply `Q%*%R` is a bit strange to me. I wonder why we should have to do so. Yes, it's how the QR program code was made, I suppose ; but why was it made in this way ...? – Andrew Nov 28 '17 at 12:08
  • 1
    R's `qr()` is just a wrapper for a linear algebra library written in Fortran called LAPACK which was written before I was born. Pivots are used when `X` is nearly rank deficient see [here](https://en.wikipedia.org/wiki/QR_decomposition#Column_pivoting). So when you tell R to do just a standard matrix multiplication, it has no indication that a pivot was used. – Michael Bird Nov 28 '17 at 13:16

0 Answers0