rsparse-matrixqr-decomposition

How to get the pivot and rank from Matrix::qr() like that of base::qr()?


When applying Matrix::qr() on the sparse matrix in R, the output is quite different from that of base::qr. There are V, beta, p, R, q but not rank and pivot. Below is a small example code. I want to detect linear dependent columns of the A sparse matrix, which requires the pivot and rank. How should I get these information?

library(Matrix)
A <- matrix(c(0, -2, 1, 0, 
              0, -4, 2, 0,
              1, -2, 1, 2,
              1, -2, 1, 2,
              1, -2, 1, 2), nrow = 5, byrow = T)
A.s <- as(A, "dgCMatrix")

qrA.M <- Matrix::qr(A.s)
qrA.R <- base::qr(A)

There is another related but not answered question, Get base::qr pivoting in Matrix::qr method


Solution

  • I would reconstruct your example matrix A a little bit:

    A <- A[, c(1,4,3,2)]
    #     [,1] [,2] [,3] [,4]
    #[1,]    0    0    1   -2
    #[2,]    0    0    2   -4
    #[3,]    1    2    1   -2
    #[4,]    1    2    1   -2
    #[5,]    1    2    1   -2
    

    You did not mention in your question why rank and pivot returned by a dense QR factorization are useful. But I think this is what you are looking for:

    dQR <- base::qr(A)
    with(dQR, pivot[1:rank])
    #[1] 1 3
    

    So columns 1 and 3 of A gives a basis for A's column space.

    I don't really understand the logic of a sparse QR factorization. The 2nd column of A is perfectly linearly dependent on the 1st column, so I expect column pivoting to take place during the factorization. But very much to my surprise, it doesn't!

    library(Matrix)
    sA <- Matrix(A, sparse = TRUE)
    sQR <- Matrix::qr(sA)
    sQR@q + 1L
    #[1] 1 2 3 4
    

    No column pivoting is done! As a result, there isn't an obvious way to determine the rank of A.

    At this moment, I could only think of performing a dense QR factorization on the R factor to get what you are looking for.

    R <- as.matrix(Matrix::qrR(sQR))
    QRR <- base::qr(R)
    with(QRR, pivot[1:rank])
    #[1] 1 3
    

    Why does this work? Well, the Q factor has orthogonal hence linearly independent columns, thus columns of R inherit linear dependence or independence of A. For a matrix with much more rows than columns, the computational costs of this 2nd QR factorization is negligible.

    I need to figure out the algorithm behind a sparse QR factorization before coming up with a better idea.