CODEWITHSUNDEEP
rmatrix
power
power

Reputation: 1760

How do I do a QR decomposition on an object of class "sparseMatrix" in the Matrix package?

I want to do a QR decomposition with the Matrix:::qr() function on a Matrix that I created with B<-as(A, "sparseMatrix"). I know that I can get the R matrix with Matrix:::qr.R(). However, I also need to the Q Matrix. There seems to be no qr.Q() function in the Matrix package. How do I get the Q matrix?

Upvotes: 2

Views: 372

Answers (1)

user1539634
user1539634

Reputation:

The Q matrix is actually stored in the V slot. It seems that the current R Matrix version contains a bug --- it just mysteriously adds zero rows into the matrix a before doing the qr decomposition. I wish the developers could come and explain it. Therefore the following codes 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 have tested by randomly set the matrix size and sparsity and they work smoothly. However this is of the style "just make it work without knowing why", and is posted here only for discussion. Because I have not hacked into the implementation details of the package "Matrix".

Upvotes: 1

Related Questions