R sparseQR-class 稀疏 QR 分解

说明

sparseQR 是 m \times n ( m \ge n ) 实数矩阵的稀疏、行和 column-pivoted QR 分解类，具有一般形式

或(等效地)

其中 P_1 和 P_2 是置换矩阵，Q = \prod_{j = 1}^{n} H_j 是 m \times m 正交矩阵(Q_1 包含第一个 n 列向量)，等于 n Householder 矩阵 H_j 的乘积, 和R 是具有非负对角线元素的 m \times n 上梯形矩阵(R_1 包含第一个 n 行向量，并且是上三角矩阵)。

用法

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

参数

细节

qr.Q 的方法不会返回 Q ，而是返回(也是正交的)乘积 P_1' Q 。此行为在代数上与 base 实现一致(请参阅 qr )，可以通过注意到 base 中的 qr.default 不旋转行来看出这一点，从而将 P_1 限制为单位矩阵。因此 qr.Q(qr.default(x)) 也返回 P_1' Q 。

同样， qr.qy 和 qr.qty 的方法在左侧乘以 P_1' Q 和 Q' P_1 而不是 Q 和 Q' 。

预期 qr.Q (或 qr.R 、 qr.qy 、 qr.qty )的值是根据 “equivalent” 稀疏和密集分解(例如，qr(x) 和 qr(as(x, "matrix")) 对于类的 x 计算得出的)是错误的dgCMatrix ) 比较相等。底层分解算法有很大不同，特别是它们采用不同的旋转策略，并且一般来说，即使对于固定的 P_1 和 P_2 ，分解也不是唯一的。

另一方面， qr.X 、 qr.coef 、 qr.fitted 和 qr.resid 的值是明确定义的，在这些情况下，稀疏和密集计算应该比较相等(在一定的容差范围内)。

qr.R 的方法是 qrR 的简单包装，但默认情况下不是 back-permuting 并且从不给出行名称。它不支持 backPermute = TRUE 直到 Matrix 1.6-0 ，因此如果 Matrix >= 1.6-0 未知，则需要 back-permuted 结果的代码应调用 qrR。

插槽

Dim , Dimnames

从虚拟类MatrixFactorization继承。

beta

长度为 Dim[2] 的数值向量，用于构造 Householder 矩阵；请参阅下面的V。

V

具有 Dim[2] 列的类 dgCMatrix 的对象。行数 nrow(V) 至少为 Dim[1] 且最多为 Dim[1]+Dim[2] 。 V 是下梯形，其列向量生成组成正交 Q 因子的 Householder 矩阵 H_j。具体来说， H_j 被构造为 diag(Dim[1]) - beta[j] * tcrossprod(V[, j]) 。

R

类dgCMatrix 的对象，具有nrow(V) 行和Dim[2] 列。 R 是具有非负对角线元素的上梯形 R 因子。

p , q

长度从 0 开始的整数向量nrow(V)和Dim[2]，分别指定应用于因式分解矩阵的行和列的排列。q长度为 0 的值是有效的并且等价于恒等排列，意味着没有列旋转。使用R语法，矩阵P_1 A P_2正是A[p+1, q+1](A[p+1, ]当q长度为 0)。

扩展

直接类 QR 。类 MatrixFactorization ，按类 QR ，距离 2。

实例化

对象可以通过 new("sparseQR", ...) 形式的调用直接生成，但它们更通常作为继承自 sparseMatrix (通常是 dgCMatrix )的 x 的 qr(x) 值获得。

方法

determinant

signature(from = "sparseQR", logarithm = "logical") ：计算因式分解矩阵A 的行列式或其对数。

expand1

signature(x = "sparseQR") ：参见expand1-methods 。

expand2

signature(x = "sparseQR") ：参见expand2-methods 。

qr.Q

signature(qr = "sparseQR") ：返回 dgeMatrix P_1' Q 或 P_1' Q_1 ，具体取决于可选参数 complete 。默认为 FALSE ，表示 P_1' Q_1 。

qr.R

signature(qr = "sparseQR") : qrR 返回 R 、 R_1 、 R P2' 或 R_1 P2' ，具体取决于可选参数 complete 和 backPermute 。两种情况下的默认值都是 FALSE ，表示 R_1 ，以与 base 兼容。这种情况下结果的类是 dtCMatrix 。在其他三种情况下，它是 dgCMatrix 。

qr.X

signature(qr = "sparseQR") ：默认情况下，将 A 作为 dgeMatrix 返回。如果 m > n 和可选参数 ncol 大于 n ，则结果将用 P_1' Q J 进行扩充，其中 J 由 m \times m 单位矩阵的 (n+1) 到 ncol 列组成。

qr.coef

signature(qr = "sparseQR", y = .) ：以 dgeMatrix 或向量形式返回左侧 y 乘以 P_2 R_1^{-1} Q_1' P_1 的结果。

qr.fitted

signature(qr = "sparseQR", y = .) ：以 dgeMatrix 或向量形式返回左侧 y 乘以 P_1' Q_1 Q_1' P_1 的结果。

qr.resid

signature(qr = "sparseQR", y = .) ：以 dgeMatrix 或向量形式返回左侧 y 乘以 P_1' Q_2 Q_2' P_1 的结果。

qr.qty

signature(qr = "sparseQR", y = .) ：以 dgeMatrix 或向量形式返回左侧 y 乘以 Q' P_1 的结果。

qr.qy

signature(qr = "sparseQR", y = .) ：以 dgeMatrix 或向量形式返回左侧 y 乘以 P_1' Q 的结果。

solve

signature(a = "sparseQR", b = .) ：参见solve-methods 。

例子

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

参考

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

也可以看看

类dgCMatrix。

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

qr-methods 用于 Matrix 中定义的方法。

通用函数 expand1 和 expand2 。

QR 分解的许多辅助函数： qr.Q , qr.R , qr.X , qr.coef , qr.fitted , qr.resid , qr.qty , qr.qy 和 qr.solve 。