Methods for QR Factorization

Computes the pivoted QR factorization of an \(m \times n\) real matrix \(A\), which has the general form $$P_{1} A P_{2} = Q R$$ or (equivalently) $$A = P_{1}' Q R 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 equal to the product of \(n\) Householder matrices \(H_{j}\), and \(R\) is an \(m \times n\) upper trapezoidal matrix.

denseMatrix use the default method implemented in base, namely qr.default. It is built on LINPACK routine dqrdc and LAPACK routine dgeqp3, which do not pivot rows, so that \(P_{1}\) is an identity matrix.

Methods for sparseMatrix are built on CXSparse routines cs_sqr and cs_qr, which require \(m \ge n\).

qr(x, ...)
# S4 method for class 'dgCMatrix'
qr(x, order = 3L, ...)

Arguments

x

a finite matrix or Matrix to be factorized, satisfying nrow(x) >= ncol(x) if sparse.

order

an integer in 0:3 passed to CXSparse routine cs_sqr, indicating a strategy for choosing the column permutation \(P_{2}\). 0 means no column permutation. 1, 2, and 3 indicate a fill-reducing ordering of \(A + A'\), \(\tilde{A}' \tilde{A}\), and \(A' A\), where \(\tilde{A}\) is \(A\) with “dense” rows removed. Do not set to 0 unless you know that the column order of \(A\) is already sensible.

...

further arguments passed to or from methods.

Value

An object representing the factorization, inheriting from virtual S4 class QR or S3 class qr. The specific class is qr unless x inherits from virtual class sparseMatrix, in which case it is sparseQR.

Details

If x is sparse and structurally rank deficient, having structural rank \(r < n\), then x is augmented with \((n-r)\) rows of (partly non-structural) zeros, such that the augmented matrix has structural rank \(n\). This augmented matrix is factorized as described above: $$P_1 A P_2 = P_1 \begin{bmatrix} A_{0} \\ 0 \end{bmatrix} P_2 = Q R$$ where \(A_0\) denotes the original, user-supplied \((m-(n-r)) \times n\) matrix.

See also

Class sparseQR and its methods.

Class dgCMatrix.

Generic function qr from base, whose default method qr.default “defines” the S3 class qr of dense QR factorizations.

Generic functions expand1 and expand2, for constructing matrix factors from the result.

Generic functions Cholesky, BunchKaufman, Schur, and lu, for computing other factorizations.

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

showMethods("qr", inherited = FALSE)
#> Function: qr (package base)
#> x="ANY"
#> x="dgCMatrix"
#> x="sparseMatrix"
#> 

## Rank deficient: columns 3 {b2} and 6 {c3} are "extra"
M <- as(cbind(a1 = 1,
              b1 = rep(c(1, 0), each = 3L),
              b2 = rep(c(0, 1), each = 3L),
              c1 = rep(c(1, 0, 0), 2L),
              c2 = rep(c(0, 1, 0), 2L),
              c3 = rep(c(0, 0, 1), 2L)),
        "CsparseMatrix")
rownames(M) <- paste0("r", seq_len(nrow(M)))
b <- 1:6
eps <- .Machine$double.eps

## .... [1] full rank ..................................................
## ===> a least squares solution of A x = b exists
##      and is unique _in exact arithmetic_

(A1 <- M[, -c(3L, 6L)])
#> 6 x 4 sparse Matrix of class "dgCMatrix"
#>    a1 b1 c1 c2
#> r1  1  1  1  .
#> r2  1  1  .  1
#> r3  1  1  .  .
#> r4  1  .  1  .
#> r5  1  .  .  1
#> r6  1  .  .  .
(qr.A1 <- qr(A1))
#> QR factorization of Formal class 'sparseQR' [package "Matrix"] with 7 slots
#>   ..@ beta    : num [1:4] 0.293 0.293 0.5 0.423
#>   ..@ V       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:9] 0 3 1 2 2 4 5 3 4
#>   .. .. ..@ p       : int [1:5] 0 2 4 7 9
#>   .. .. ..@ Dim     : int [1:2] 6 4
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:9] 2.41 1 2.41 1 1.41 ...
#>   .. .. ..@ factors : list()
#>   ..@ R       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:9] 0 1 0 1 2 0 1 2 3
#>   .. .. ..@ p       : int [1:5] 0 1 2 5 9
#>   .. .. ..@ Dim     : int [1:2] 6 4
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:9] -1.41 -1.41 -1.41 -1.41 -1.41 ...
#>   .. .. ..@ factors : list()
#>   ..@ p       : int [1:6] 0 1 4 3 2 5
#>   ..@ q       : int [1:4] 2 3 0 1
#>   ..@ Dim     : int [1:2] 6 4
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:6] "r1" "r2" "r3" "r4" ...
#>   .. ..$ : chr [1:4] "a1" "b1" "c1" "c2"

stopifnot(exprs = {
    rankMatrix(A1) == ncol(A1)
    { d1 <- abs(diag(qr.A1@R)); sum(d1 < max(d1) * eps) == 0L }
    rcond(crossprod(A1)) >= eps
    all.equal(qr.coef(qr.A1, b), drop(solve(crossprod(A1), crossprod(A1, b))))
    all.equal(qr.fitted(qr.A1, b) + qr.resid(qr.A1, b), b)
})
#> Warning: 'rcond' via sparse -> dense coercion

## .... [2] numerically rank deficient with full structural rank .......
## ===> a least squares solution of A x = b does not
##      exist or is not unique _in exact arithmetic_

(A2 <- M)
#> 6 x 6 sparse Matrix of class "dgCMatrix"
#>    a1 b1 b2 c1 c2 c3
#> r1  1  1  .  1  .  .
#> r2  1  1  .  .  1  .
#> r3  1  1  .  .  .  1
#> r4  1  .  1  1  .  .
#> r5  1  .  1  .  1  .
#> r6  1  .  1  .  .  1
(qr.A2 <- qr(A2))
#> QR factorization of Formal class 'sparseQR' [package "Matrix"] with 7 slots
#>   ..@ beta    : num [1:6] 2.93e-01 2.76e-01 3.98e-01 4.38e-01 1.54e+31 ...
#>   ..@ V       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:14] 0 1 1 2 3 2 3 4 3 4 ...
#>   .. .. ..@ p       : int [1:7] 0 2 5 8 11 13 14
#>   .. .. ..@ Dim     : int [1:2] 6 6
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:14] 2.41 1 -2.29 1 1 ...
#>   .. .. ..@ factors : list()
#>   ..@ R       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:18] 0 0 1 1 2 1 2 3 0 1 ...
#>   .. .. ..@ p       : int [1:7] 0 1 3 5 8 13 18
#>   .. .. ..@ Dim     : int [1:2] 6 6
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:18] -1.414 -0.707 1.581 0.632 -1.265 ...
#>   .. .. ..@ factors : list()
#>   ..@ p       : int [1:6] 2 5 0 1 3 4
#>   ..@ q       : int [1:6] 5 1 3 4 2 0
#>   ..@ Dim     : int [1:2] 6 6
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:6] "r1" "r2" "r3" "r4" ...
#>   .. ..$ : chr [1:6] "a1" "b1" "b2" "c1" ...

stopifnot(exprs = {
    rankMatrix(A2) == ncol(A2) - 2L
    { d2 <- abs(diag(qr.A2@R)); sum(d2 < max(d2) * eps) == 2L }
    rcond(crossprod(A2)) < eps

    ## 'qr.coef' computes unique least squares solution of "nearby" problem
    ## Z x = b for some full rank Z ~ A, currently without warning {FIXME} !
    tryCatch({ qr.coef(qr.A2, b); TRUE }, condition = function(x) FALSE)

    all.equal(qr.fitted(qr.A2, b) + qr.resid(qr.A2, b), b)
})
#> Warning: 'rcond' via sparse -> dense coercion

## .... [3] numerically and structurally rank deficient ................
## ===> factorization of _augmented_ matrix with
##      full structural rank proceeds as in [2]

##  NB: implementation details are subject to change; see (*) below

A3 <- M
A3[, c(3L, 6L)] <- 0
A3
#> 6 x 6 sparse Matrix of class "dgCMatrix"
#>    a1 b1 b2 c1 c2 c3
#> r1  1  1  .  1  .  .
#> r2  1  1  .  .  1  .
#> r3  1  1  .  .  .  .
#> r4  1  .  .  1  .  .
#> r5  1  .  .  .  1  .
#> r6  1  .  .  .  .  .
(qr.A3 <- qr(A3)) # with a warning ... "additional 2 row(s) of zeros"
#> Warning: matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
#> QR factorization of Formal class 'sparseQR' [package "Matrix"] with 7 slots
#>   ..@ beta    : num [1:6] 0.293 0.25 0.667 0.346 0 ...
#>   ..@ V       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:16] 0 1 1 2 3 6 7 2 3 6 ...
#>   .. .. ..@ p       : int [1:7] 0 2 7 11 14 15 16
#>   .. .. ..@ Dim     : int [1:2] 8 6
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:16] 2.41 1 2 1 1 ...
#>   .. .. ..@ factors : list()
#>   ..@ R       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:8] 0 0 1 1 2 0 1 3
#>   .. .. ..@ p       : int [1:7] 0 1 3 5 8 8 8
#>   .. .. ..@ Dim     : int [1:2] 8 6
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:8] -1.41 -1.41 -2 -1 -1 ...
#>   .. .. ..@ factors : list()
#>   ..@ p       : int [1:8] 1 4 0 2 6 7 3 5
#>   ..@ q       : int [1:6] 4 0 3 1 2 5
#>   ..@ Dim     : int [1:2] 6 6
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:6] "r1" "r2" "r3" "r4" ...
#>   .. ..$ : chr [1:6] "a1" "b1" "b2" "c1" ...

stopifnot(exprs = {
    ## sparseQR object preserves the unaugmented dimensions (*)
    dim(qr.A3  ) == dim(A3)
    dim(qr.A3@V) == dim(A3) + c(2L, 0L)
    dim(qr.A3@R) == dim(A3) + c(2L, 0L)

    ## The augmented matrix remains numerically rank deficient
    rankMatrix(A3) == ncol(A3) - 2L
    { d3 <- abs(diag(qr.A3@R)); sum(d3 < max(d3) * eps) == 2L }
    rcond(crossprod(A3)) < eps
})
#> Warning: 'rcond' via sparse -> dense coercion

## Auxiliary functions accept and return a vector or matrix
## with dimensions corresponding to the unaugmented matrix (*),
## in all cases with a warning
qr.coef  (qr.A3, b)
#> Warning: matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
#> a1 b1 b2 c1 c2 c3 
#>  6 -3  0 -2 -1  0 
qr.fitted(qr.A3, b)
#> Warning: matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
#> [1] 1 2 3 4 5 6
qr.resid (qr.A3, b)
#> Warning: matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
#> [1]  1.912284e-16  1.218711e-16 -3.130995e-16 -1.912284e-16 -1.218711e-16
#> [6]  3.130995e-16

## .... [4] yet more examples ..........................................

## By disabling column pivoting, one gets the "vanilla" factorization
## A = Q~ R, where Q~ := P1' Q is orthogonal because P1 and Q are

(qr.A1.pp <- qr(A1, order = 0L)) # partial pivoting
#> QR factorization of Formal class 'sparseQR' [package "Matrix"] with 7 slots
#>   ..@ beta    : num [1:4] 0.118 0.517 0.554 0.907
#>   ..@ V       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:18] 0 1 2 3 4 5 1 2 3 4 ...
#>   .. .. ..@ p       : int [1:5] 0 6 11 15 18
#>   .. .. ..@ Dim     : int [1:2] 6 4
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:18] 3.45 1 1 1 1 ...
#>   .. .. ..@ factors : list()
#>   ..@ R       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. ..@ i       : int [1:9] 0 0 1 0 1 2 0 2 3
#>   .. .. ..@ p       : int [1:5] 0 1 3 6 9
#>   .. .. ..@ Dim     : int [1:2] 6 4
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:9] -2.45 -1.22 -1.22 -8.16e-01 -2.22e-16 ...
#>   .. .. ..@ factors : list()
#>   ..@ p       : int [1:6] 0 1 2 3 4 5
#>   ..@ q       : int(0) 
#>   ..@ Dim     : int [1:2] 6 4
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:6] "r1" "r2" "r3" "r4" ...
#>   .. ..$ : chr [1:4] "a1" "b1" "c1" "c2"

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

stopifnot(exprs = {
    length(qr.A1   @q) == ncol(A1)
    length(qr.A1.pp@q) == 0L # indicating no column pivoting
    ae2(A1[, qr.A1@q + 1L], qr.Q(qr.A1   ) %*% qr.R(qr.A1   ))
    ae2(A1                , qr.Q(qr.A1.pp) %*% qr.R(qr.A1.pp))
})