1

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

Chris Lotus
  • 33
  • 4
  • Have you looked at each of these three help topics, including the examples: `qr`, `sparseQR-class`, `rankMatrix`? – Mikael Jagan Aug 09 '22 at 23:17
  • @MikaelJagan Yes, I read them but it seems that none of them stated how to obtain the pivot and rank as in the dense case. – Chris Lotus Aug 12 '22 at 19:18

2 Answers2

2

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.

Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • `Matrix` uses a somewhat modified version of Timothy Davis' algorithm for sparse QR factorization, which is laid out in his text _Direct Methods for Sparse Linear Systems_ (2006). Column pivoting serves a different purpose here: to make Q and R as sparse as possible, rather than to facilitate numerical rank approximation (as in the dense case). – Mikael Jagan Aug 23 '22 at 03:13
0

I've been looking at a similar problem and I ended up not relying on Matrix::qr() to calculate rank and to detect linear dependency. Instead I programmed the function GaussIndependent in the package SSBtools.

In the package examples I included an example that demonstrates wrong conclusion from rankMatrix(x, method = "qr"). Input x is a 44*20 dummy matrix.

Starting with your example matrix, A.s:

library(SSBtools)
GaussIndependent(A.s) # List of logical vectors specifying independent rows and columns
# $rows
# [1]  TRUE FALSE  TRUE FALSE FALSE
#
# $columns
# [1]  TRUE  TRUE FALSE FALSE


GaussRank(A.s) # the rank
# [1] 2
Øyvind Langsrud
  • 181
  • 6
  • Dear Oyvind, thanks for your answer. I tested the function. it looks like it can return different results from `base::qr` when the number of columns is large. I have opened an issue in the pkg GitHub. Please check it a bit when you got time. Thanks! – Chris Lotus Sep 09 '22 at 01:54
  • Dear Chris, thanks for the example on GitHub. As I comment “I will change the documentation and warn that the function is written primarily for integer/dummy matrices”. – Øyvind Langsrud Sep 09 '22 at 10:37