sparseQR-class {Matrix}R Documentation

Sparse QR Factorizations

Description

sparseQR is the class of sparse, row- and column-pivoted QR factorizations of m×nm \times n (mnm \ge n) real matrices, having the general form

P1AP2=QR=[Q1Q2][R10]=Q1R1P_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=P1QRP2=P1[Q1Q2][R10]P2=P1Q1R1P2A = 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 P1P_1 and P2P_2 are permutation matrices, Q=j=1nHjQ = \prod_{j = 1}^{n} H_j is an m×mm \times m orthogonal matrix (Q1Q_1 contains the first nn column vectors) equal to the product of nn Householder matrices HjH_j, and RR is an m×nm \times n upper trapezoidal matrix (R1R_1 contains the first nn row vectors and is upper triangular).

Usage

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 RR should be returned instead of R1R_1.

backPermute

a logical indicating if RR or R1R_1 should be multiplied on the right by P2P_2'.

row.names

a logical indicating if dimnames(qr)[1] should be propagated unpermuted to the result. If complete = FALSE, then only the first nn names are kept.

Details

The method for qr.Q does not return QQ but rather the (also orthogonal) product P1QP_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 P1P_1 to be an identity matrix. It follows that qr.Q(qr.default(x)) also returns P1QP_1' Q.

Similarly, the methods for qr.qy and qr.qty multiply on the left by P1QP_1' Q and QP1Q' P_1 rather than QQ and QQ'.

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 P1P_1 and P2P_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 HjH_j that compose the orthogonal QQ factor. Specifically, HjH_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 RR 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 P1AP2P_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 AA 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 P1QP_1' Q or P1Q1P_1' Q_1, depending on optional argument complete. The default is FALSE, indicating P1Q1P_1' Q_1.

qr.R

signature(qr = "sparseQR"): qrR returns RR, R1R_1, RP2R P2', or R1P2R_1 P2', depending on optional arguments complete and backPermute. The default in both cases is FALSE, indicating R1R_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 AA as a dgeMatrix, by default. If m>nm > n and optional argument ncol is greater than nn, then the result is augmented with P1QJP_1' Q J, where JJ is composed of columns (n+1)(n+1) through ncol of the m×mm \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 P2R11Q1P1P_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 P1Q1Q1P1P_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 P1Q2Q2P1P_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 QP1Q' P_1.

qr.qy

signature(qr = "sparseQR", y = .): returns as a dgeMatrix or vector the result of multiplying y on the left by P1QP_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))
})

[Package Matrix version 1.7-0 Index]