R qr-methods QR 分解方法

说明

计算 m \times n 实数矩阵 A 的旋转 QR 分解，其具有一般形式

或(等效地)

其中 P_{1} 和 P_{2} 是置换矩阵，Q = \prod_{j = 1}^{n} H_{j} 是 m \times m 正交矩阵，等于 n Householder 矩阵 H_{j} 的乘积，R 是 m \times n 上陷阱矩状矩阵。

denseMatrix 使用 base 中实现的默认方法，即 qr.default 。它基于 LINPACK 例程 dqrdc 和 LAPACK 例程 dgeqp3 构建，它们不进行行旋转，因此 P_{1} 是单位矩阵。

sparseMatrix 的方法基于 CSparse 例程 cs_sqr 和 cs_qr 构建，需要 m \ge n 。

用法

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

参数

细节

如果 x 稀疏且结构秩不足，具有结构秩 r < n ，则用 (n-r) 行(部分非结构性)零来增强 x ，使得增强矩阵具有结构秩 n 。该增广矩阵按如上所述进行因式分解：

其中A_0 表示用户提供的原始(m-(n-r)) \times n 矩阵。

值

表示因式分解的对象，继承自虚拟 S4 类 QR 或 S3 类 qr 。特定类是 qr ，除非 x 继承自虚拟类 sparseMatrix ，在这种情况下它是 sparseQR 。

例子

showMethods("qr", inherited = FALSE)

## 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)])
(qr.A1 <- qr(A1))

stopifnot(exprs = {
    rankMatrix(A1) == ncol(A1)
    { d1 <- 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)
})

## .... [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)
(qr.A2 <- qr(A2))

stopifnot(exprs = {
    rankMatrix(A2) == ncol(A2) - 2L
    { d2 <- 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)
})

## .... [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
(qr.A3 <- qr(A3)) # with a warning ... "additional 2 row(s) of zeros"

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 <- diag(qr.A3@R); sum(d3 < max(d3) * eps) == 2L }
    rcond(crossprod(A3)) < eps
})

## 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)
qr.fitted(qr.A3, b)
qr.resid (qr.A3, b)

## .... [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

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

参考

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

也可以看看

类sparseQR 及其方法。

类dgCMatrix。

来自 base 的通用函数 qr ，其默认方法 qr.default “defines” 是密集 QR 分解的 S3 类 qr。

通用函数 expand1 和 expand2 ，用于根据结果构造矩阵因子。

通用函数 Cholesky 、 BunchKaufman 、 Schur 和 lu 用于计算其他因式分解。