How to do QR decomposition on the "sparseMatrix" class object in the Matrix package?
I want to do QR decomposition using a function Matrix:::qr()
on a matrix that I created with B<-as(A, "sparseMatrix")
. I know that I can get an R matrix with Matrix:::qr.R()
. However, I also need a Q-matrix. There seems to be no qr.Q () function in the Matrix package. How to get the Q-matrix?
source to share
The matrix is Q
actually stored in the slot V
. The current version of R Matrix seems to have a bug - it just mysteriously adds zero rows to the matrix a
before performing qr decomposition. I would like the developers to come and explain this. Therefore, the following codes will help you recover both R and Q:
gx.qr.Q <- function(a){
if(is(a, "qr")){
return(qr.Q(a, complete = T))
}else if(is(a, "sparseQR")){
Q <- diag(nrow = a@V@Dim[1], ncol = a@V@Dim[1])
for(i in rev(1:a@V@Dim[2])){
Q <- Q - (a@V[ ,i] * a@beta[i]) %*% (a@V[ ,i] %*% Q)
}
return(Q[order(a@p), ][1:a@Dim[1], 1:a@Dim[1]])
}else{
stop(gettextf("gx.qr.Q() fails on class '%s'", class(a)[1]))
}
}
gx.qr.R <- function(a){
if(is(a, "qr")){
return(qr.R(a, complete = T)[ ,order(a$pivot)])
}else if(is(a, "sparseQR")){
if(length(a@q) == 0){
return(a@R[1:a@Dim[1], ])
}else{
return(a@R[1:a@Dim[1] ,order(a@q)])
}
}else{
stop(gettextf("gx.qr.R() fails on class '%s'", class(a)[1]))
}
}
I've randomly tested matrix size and sparsity and they work smoothly. However, this has a style of "just make it work without knowing why" and is posted here for discussion only. Because I haven't cracked the implementation details of the Matrix package.
source to share