Sparse QR Factorizations

sparseQR is the class of sparse, row- and column-pivoted QR factorizations of \(m \times n\) (\(m \ge n\)) real matrices, having the general form P_1 A P_2 = Q R = \\\\beginbmatrix Q_1 & Q_2 \\\\endbmatrix \\\\beginbmatrix R_1 \\\\\\\\ 0 \\\\endbmatrix = Q_1 R_1P1 * A * P2 = Q * R = [Q1, Q2] * [R1; 0] = Q1 * R1 or (equivalently) A = P_1' Q R P_2' = P_1' \\\\beginbmatrix Q_1 & Q_2 \\\\endbmatrix \\\\beginbmatrix R_1 \\\\\\\\ 0 \\\\endbmatrix P_2' = P_1' Q_1 R_1 P_2'A = P1' * Q * R * P2' = P1' * [Q1, Q2] * [R1; 0] * P2' = P1' * Q1 * R1 * P2' where \(P_1\) and \(P_2\) are permutation matrices, \(Q = \prod_{j = 1}^{n} H_j\) is an \(m \times m\) orthogonal matrix (\(Q_1\) contains the first \(n\) column vectors) equal to the product of \(n\) Householder matrices \(H_j\), and \(R\) is an \(m \times n\) upper trapezoidal matrix (\(R_1\) contains the first \(n\) row vectors and is upper triangular).

qrR(qr, complete = FALSE, backPermute = TRUE, row.names = TRUE)

Arguments

qr: an object of class sparseQR, almost always the result of a call to generic function qr with sparse x.
complete: a logical indicating if \(R\) should be returned instead of \(R_1\).
backPermute: a logical indicating if \(R\) or \(R_1\) should be multiplied on the right by \(P_2'\).
row.names: a logical indicating if dimnames(qr)[1] should be propagated unpermuted to the result. If complete = FALSE, then only the first \(n\) names are kept.

Slots

Dim, Dimnames: inherited from virtual class MatrixFactorization.
beta: a numeric vector of length Dim[2], used to construct Householder matrices; see V below.
V: an object of class dgCMatrix with Dim[2] columns. The number of rows nrow(V) is at least Dim[1] and at most Dim[1]+Dim[2]. V is lower trapezoidal, and its column vectors generate the Householder matrices \(H_j\) that compose the orthogonal \(Q\) factor. Specifically, \(H_j\) is constructed as diag(Dim[1]) - beta[j] * tcrossprod(V[, j]).
R: an object of class dgCMatrix with nrow(V) rows and Dim[2] columns. R is the upper trapezoidal \(R\) factor.
p, q: 0-based integer vectors of length nrow(V) and Dim[2], respectively, specifying the permutations applied to the rows and columns of the factorized matrix. q of length 0 is valid and equivalent to the identity permutation, implying no column pivoting. Using R syntax, the matrix \(P_1 A P_2\) is precisely A[p+1, q+1] (A[p+1, ] when q has length 0).

Extends

Class QR, directly. Class MatrixFactorization, by class QR, distance 2.

Instantiation

Objects can be generated directly by calls of the form new("sparseQR", ...), but they are more typically obtained as the value of qr(x) for x inheriting from sparseMatrix (often dgCMatrix).

Methods

determinant: signature(from = "sparseQR", logarithm = "logical"): computes the determinant of the factorized matrix \(A\) or its logarithm.
expand1: signature(x = "sparseQR"): see expand1-methods.
expand2: signature(x = "sparseQR"): see expand2-methods.
qr.Q: signature(qr = "sparseQR"): returns as a dgeMatrix either \(P_1' Q\) or \(P_1' Q_1\), depending on optional argument complete. The default is FALSE, indicating \(P_1' Q_1\).
qr.R: signature(qr = "sparseQR"): qrR returns \(R\), \(R_1\), \(R P2'\), or \(R_1 P2'\), depending on optional arguments complete and backPermute. The default in both cases is FALSE, indicating \(R_1\), for compatibility with base. The class of the result in that case is dtCMatrix. In the other three cases, it is dgCMatrix.
qr.X: signature(qr = "sparseQR"): returns \(A\) as a dgeMatrix, by default. If \(m > n\) and optional argument ncol is greater than \(n\), then the result is augmented with \(P_1' Q J\), where \(J\) is composed of columns \((n+1)\) through ncol of the \(m \times m\) identity matrix.
qr.coef: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by \(P_2 R_1^{-1} Q_1' P_1\).
qr.fitted: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by \(P_1' Q_1 Q_1' P_1\).
qr.resid: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by \(P_1' Q_2 Q_2' P_1\).
qr.qty: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by \(Q' P_1\).
qr.qy: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by \(P_1' Q\).
solve: signature(a = "sparseQR", b = .): see solve-methods.

Details

The method for qr.Q does not return \(Q\) but rather the (also orthogonal) product \(P_1' Q\). This behaviour is algebraically consistent with the base implementation (see qr), which can be seen by noting that qr.default in base does not pivot rows, constraining \(P_1\) to be an identity matrix. It follows that qr.Q(qr.default(x)) also returns \(P_1' Q\).

Similarly, the methods for qr.qy and qr.qty multiply on the left by \(P_1' Q\) and \(Q' P_1\) rather than \(Q\) and \(Q'\).

It is wrong to expect the values of qr.Q (or qr.R, qr.qy, qr.qty) computed from “equivalent” sparse and dense factorizations (say, qr(x) and qr(as(x, "matrix")) for x of class dgCMatrix) to compare equal. The underlying factorization algorithms are quite different, notably as they employ different pivoting strategies, and in general the factorization is not unique even for fixed \(P_1\) and \(P_2\).

On the other hand, the values of qr.X, qr.coef, qr.fitted, and qr.resid are well-defined, and in those cases the sparse and dense computations should compare equal (within some tolerance).

The method for qr.R is a simple wrapper around qrR, but not back-permuting by default and never giving row names. It did not support backPermute = TRUE until Matrix 1.6-0, hence code needing the back-permuted result should call qrR if Matrix >= 1.6-0 is not known.

References

Davis, T. A. (2006). Direct methods for sparse linear systems. Society for Industrial and Applied Mathematics. doi:10.1137/1.9780898718881

Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944

Examples

showClass("sparseQR")
#> Class "sparseQR" [package "Matrix"]
#> 
#> Slots:
#>                                                                             
#> Name:       beta         V         R         p         q       Dim  Dimnames
#> Class:   numeric dgCMatrix dgCMatrix   integer   integer   integer      list
#> 
#> Extends: 
#> Class "QR", directly
#> Class "MatrixFactorization", by class "QR", distance 2
set.seed(2)

m <- 300L
n <- 60L
A <- rsparsematrix(m, n, 0.05)

## With dimnames, to see that they are propagated :
dimnames(A) <- dn <- list(paste0("r", seq_len(m)),
                          paste0("c", seq_len(n)))

(qr.A <- qr(A))
#> QR factorization of Formal class 'sparseQR' [package "Matrix"] with 7 slots
#>   ..@ beta    : num [1:60] 0.092 0.1458 0.0993 0.066 0.1095 ...
#>   ..@ V       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:10202] 0 56 57 58 59 174 181 217 269 272 ...
#>   .. .. ..@ p       : int [1:61] 0 12 23 32 44 70 80 90 115 128 ...
#>   .. .. ..@ Dim     : int [1:2] 300 60
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:10202] -3.22 -0.13 1 -0.56 0.83 ...
#>   .. .. ..@ factors : list()
#>   ..@ R       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:1595] 0 1 2 3 3 4 5 6 5 6 ...
#>   .. .. ..@ p       : int [1:61] 0 1 2 3 4 6 7 8 11 12 ...
#>   .. .. ..@ Dim     : int [1:2] 300 60
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:1595] 3.3731 3.1545 2.9568 4.7129 0.0271 ...
#>   .. .. ..@ factors : list()
#>   ..@ p       : int [1:300] 4 23 24 15 99 42 96 101 21 30 ...
#>   ..@ q       : int [1:60] 8 23 26 46 47 2 21 50 55 0 ...
#>   ..@ Dim     : int [1:2] 300 60
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:300] "r1" "r2" "r3" "r4" ...
#>   .. ..$ : chr [1:60] "c1" "c2" "c3" "c4" ...
str(e.qr.A <- expand2(qr.A, complete = FALSE), max.level = 2L)
#> List of 4
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ Q1 :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
#>  $ R1 :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ P2.:Formal class 'pMatrix' [package "Matrix"] with 5 slots
str(E.qr.A <- expand2(qr.A, complete =  TRUE), max.level = 2L)
#> List of 4
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ Q  :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
#>  $ R  :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>  $ P2.:Formal class 'pMatrix' [package "Matrix"] with 5 slots

t(sapply(e.qr.A, dim))
#>     [,1] [,2]
#> P1.  300  300
#> Q1   300   60
#> R1    60   60
#> P2.   60   60
t(sapply(E.qr.A, dim))
#>     [,1] [,2]
#> P1.  300  300
#> Q    300  300
#> R    300   60
#> P2.   60   60

## Horribly inefficient, but instructive :
slowQ <- function(V, beta) {
    d <- dim(V)
    Q <- diag(d[1L])
    if(d[2L] > 0L) {
        for(j in d[2L]:1L) {
            cat(j, "\n", sep = "")
            Q <- Q - (beta[j] * tcrossprod(V[, j])) %*% Q
        }
    }
    Q
}

ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...)
ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...)

## A ~ P1' Q R P2' ~ P1' Q1 R1 P2' in floating point
stopifnot(exprs = {
    identical(names(e.qr.A), c("P1.", "Q1", "R1", "P2."))
    identical(names(E.qr.A), c("P1.", "Q" , "R" , "P2."))
    identical(e.qr.A[["P1."]],
              new("pMatrix", Dim = c(m, m), Dimnames = c(dn[1L], list(NULL)),
                  margin = 1L, perm = invertPerm(qr.A@p, 0L, 1L)))
    identical(e.qr.A[["P2."]],
              new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]),
                  margin = 2L, perm = invertPerm(qr.A@q, 0L, 1L)))
    identical(e.qr.A[["R1"]], triu(E.qr.A[["R"]][seq_len(n), ]))
    identical(e.qr.A[["Q1"]],      E.qr.A[["Q"]][, seq_len(n)] )
    identical(E.qr.A[["R"]], qr.A@R)
 ## ae1(E.qr.A[["Q"]], slowQ(qr.A@V, qr.A@beta))
    ae1(crossprod(E.qr.A[["Q"]]), diag(m))
    ae1(A, with(e.qr.A, P1. %*% Q1 %*% R1 %*% P2.))
    ae1(A, with(E.qr.A, P1. %*% Q  %*% R  %*% P2.))
    ae2(A.perm <- A[qr.A@p + 1L, qr.A@q + 1L], with(e.qr.A, Q1 %*% R1))
    ae2(A.perm                               , with(E.qr.A, Q  %*% R ))
})

## More identities
b <- rnorm(m)
stopifnot(exprs = {
    ae1(qrX <- qr.X     (qr.A   ), A)
    ae2(qrQ <- qr.Q     (qr.A   ), with(e.qr.A, P1. %*% Q1))
    ae2(       qr.R     (qr.A   ), with(e.qr.A, R1))
    ae2(qrc <- qr.coef  (qr.A, b), with(e.qr.A, solve(R1 %*% P2., t(qrQ)) %*% b))
    ae2(qrf <- qr.fitted(qr.A, b), with(e.qr.A, tcrossprod(qrQ) %*% b))
    ae2(qrr <- qr.resid (qr.A, b), b - qrf)
    ae2(qrq <- qr.qy    (qr.A, b), with(E.qr.A, P1. %*% Q %*% b))
    ae2(qr.qty(qr.A, qrq), b)
})
#> Error: ae2(qrf <- qr.fitted(qr.A, b), with(e.qr.A, tcrossprod(qrQ) %*%  .... is not TRUE

## Sparse and dense computations should agree here
qr.Am <- qr(as(A, "matrix")) # <=> qr.default(A)
stopifnot(exprs = {
    ae2(qrX, qr.X     (qr.Am   ))
    ae2(qrc, qr.coef  (qr.Am, b))
    ae2(qrf, qr.fitted(qr.Am, b))
    ae2(qrr, qr.resid (qr.Am, b))
})
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'target' in selecting a method for function 'all.equal': error in evaluating the argument 'obj' in selecting a method for function 'unname': object 'qrr' not found