I need to perform a QR-decomposition of a sparse matrix, so I've been
trying to use the Matrix package. Unfortunately I don't seem to be getting
exactly the same results as if I had used the qr() command from the base
package. Here is an example of what I'm doing.
> spdata <-rpois(50,1)
> y <- rnorm(10,0,1)
> S <- matrix(round(ifelse(spdata==0,0,1/spdata),2),ncol=5,nrow=10)
> Ssp <- as(S,"sparseMatrix")
>
> decomp <- qr(S)
> Y <- qr.qty(decomp,y)
> R <- qr.R(decomp)
>
> decompSp <- qr(Ssp)
> Ysp <- qr.qty(decompSp,y)
> Rsp <- qr.R(decompSp)
Warning message:
In qr.R(decompSp) :
qr.R(<sparse>) may differ from qr.R(<dense>) because of permutations
>
>
> R
[,1] [,2] [,3] [,4] [,5]
[1,] -2.54951 -1.372813 -1.7650452 -1.53362818 -1.2551433
[2,] 0.00000 -1.270978 0.3328751 -0.04297115 0.4902343
[3,] 0.00000 0.000000 -1.3318444 -0.47857013 0.1340783
[4,] 0.00000 0.000000 0.0000000 -1.06532095 0.2109142
[5,] 0.00000 0.000000 0.0000000 0.00000000 -0.9548947
> Rsp
5 x 5 sparse Matrix of class "dtCMatrix"
[1,] 2.54951 1.372813 1.7650452 1.53362818 1.2551433
[2,] . 1.270978 -0.3328751 0.04297115 -0.4902343
[3,] . . 1.3318444 0.47857013 -0.1340783
[4,] . . . 1.06532095 -0.2109142
[5,] . . . . 0.9548947
>
>
> Y
[1] -1.3021860 2.4266858 1.6970921 -1.8412319 -0.4492685 1.2129510
[7] 0.5840221 0.4108283 0.3204408 -0.2110558
> Ysp
10 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 1.30218603
[2,] -2.42668582
[3,] -1.69709211
[4,] 1.84123190
[5,] 0.44926846
[6,] -0.86991028
[7,] 0.95121324
[8,] -0.09605988
[9,] 0.62723267
[10,] -0.25325849
It gives a warning saying that R may be different because of permutations,
which I believe has something to do with the permutation of columns of the
original matrix that qr() uses in its calculation. If anyone could provide
some insight on how this QR method works for sparse matrices, it would be
greatly appreciated.
Karl Pazdernik
Iowa State University
Department of Statistics & Statistical Laboratory
Snedecor Hall
Ames, IA 50011
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
