sparseQR-class function

Sparse QR Factorizations

sparseQR is the class of sparse, row- and column-pivoted QR factorizations of $m-by-n$ ($m >= n$) real matrices, having the general form [REMOVE_ME]$P_1 A P_2 = Q R = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1 R_1P1 * A * P2 = Q * R = [Q1, Q2] * [R1; 0] = Q1 * R1 [REMOVE_ME_2]$

or (equivalently) [REMOVE_ME]$A = P_1' Q R P_2' = P_1' \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} 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' [REMOVE_ME_2]$

where $P1$ and $P2$ are permutation matrices, $Q = prod(Hj : j = 1,...,n)$

is an $m-by-m$ orthogonal matrix ($Q1$ contains the first $n$ column vectors) equal to the product of $n$ Householder matrices $Hj$, and $R$ is an $m-by-n$ upper trapezoidal matrix ($R1$ contains the first $n$ row vectors and is upper triangular).

class

Description

sparseQR is the class of sparse, row- and column-pivoted QR factorizations of $m-by-n$ ($m >= n$) real matrices, having the general form

$$P_1 A P_2 = Q R = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = 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' \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} 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 $P1$ and $P2$ are permutation matrices, $Q = prod(Hj : j = 1,...,n)$

is an $m-by-m$ orthogonal matrix ($Q1$ contains the first $n$ column vectors) equal to the product of $n$ Householder matrices $Hj$, and $R$ is an $m-by-n$ upper trapezoidal matrix ($R1$ 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 $R1$.

• backPermute: a logical indicating if $R$ or $R1$

  should be multiplied on the right by $P2'$.

• 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 $Hj$ that compose the orthogonal $Q$ factor. Specifically, $Hj$ 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 syntax, the matrix $P1 * A * P2$

   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 $P1' * Q$ or $P1' * Q1$, depending on optional argument complete. The default is FALSE, indicating $P1' * Q1$.

• qr.R: signature(qr = "sparseQR"): qrR returns $R$, $R1$, $R * P2'$, or $R1 * P2'$, depending on optional arguments complete and backPermute. The default in both cases is FALSE, indicating $R1$, 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 $P1 * Q * J$, where $J$ is composed of columns $(n+1)$ through ncol of the $m-by-m$ identity matrix.

• qr.coef: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by $P2 * R1^{-1} * Q1' * P1$.

• qr.fitted: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by $P1' * Q1 * Q1' * P1$.

• qr.resid: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by $P1' * Q2 * Q2' * P1$.

• qr.qty: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by $Q' * P1$.

• qr.qy: signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by $P1' * 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 $P1' * 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 $P1$ to be an identity matrix. It follows that qr.Q(qr.default(x)) also returns $P1' * Q$.

Similarly, the methods for qr.qy and qr.qty multiply on the left by $P1' * Q$ and $Q' * P1$

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 $P1$ and $P2$.

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.

See Also

Class dgCMatrix.

Generic function qr from base, whose default method qr.default defines

the S3 class qr of dense QR factorizations.

qr-methods for methods defined in Matrix.

Generic functions expand1 and expand2.

The many auxiliary functions for QR factorizations: qr.Q, qr.R, qr.X, qr.coef, qr.fitted, qr.resid, qr.qty, qr.qy, and qr.solve.

References

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

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

Examples

showClass("sparseQR")
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))
str(e.qr.A <- expand2(qr.A, complete = FALSE), max.level = 2L)
str(E.qr.A <- expand2(qr.A, complete =  TRUE), max.level = 2L)

t(sapply(e.qr.A, dim))
t(sapply(E.qr.A, dim))

## 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)
})

## 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))
})