Sparse QR Factorizations

Description

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 = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1 R_1

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'

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).

Usage

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

Arguments

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.

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.

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

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.

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