R chol-methods 计算矩阵的 Cholesky 因子

说明

计算 n \times n 实数、对称、半正定矩阵 A 的上三角 Cholesky 因子(可选)在旋转之后。这就是L' 中的因子

或(等效地)

其中P_{1} 是置换矩阵。

denseMatrix 的方法基于 LAPACK 例程 dpstrf 、 dpotrf 和 dpptrf 构建，后两者不会置换行或列，因此 P_{1} 是单位矩阵。

sparseMatrix 的方法基于 CHOLMOD 例程 cholmod_analyze 和 cholmod_factorize_p 构建。

用法

chol(x, ...)
## S4 method for signature 'dsyMatrix'
chol(x, pivot = FALSE, tol = -1, ...)
## S4 method for signature 'dspMatrix'
chol(x, ...)
## S4 method for signature 'dsCMatrix'
chol(x, pivot = FALSE, ...)
## S4 method for signature 'ddiMatrix'
chol(x, ...)
## S4 method for signature 'generalMatrix'
chol(x, uplo = "U", ...)
## S4 method for signature 'triangularMatrix'
chol(x, uplo = "U", ...)

参数

x	要因式分解的 finite 、对称正半定矩阵或 Matrix 。如果x是正方形但不对称，则将其视为对称；请参阅uplo。当 pivot = FALSE 时，密集 x 的方法需要正定性。稀疏(但不是对角线)x 的方法无条件地需要正定性。
pivot	指示 x 的行和列是否应旋转的逻辑。稀疏 x 的方法采用近似最小度 (AMD) 算法来减少 fill-in，即不考虑数值稳定性。
tol	finite 数字容差，仅在 pivot = TRUE 时使用。如果主元小于或等于 tol ，则因式分解算法停止。负数 tol 相当于 nrow(x) * .Machine$double.eps * max(diag(x)) 。
uplo	一个字符串，"U" 或 "L" ，指示应使用 x 的哪个三角形来计算因式分解。默认值为 "U" ，即使对于下三角 x ，也与 base 中的 chol 一致。
...	传入或传出方法的更多参数。

细节

对于继承自 diagonalMatrix 的 x ，直接计算对角线结果，无需旋转，即绕过 CHOLMOD。

对于所有其他 x ， chol(x, pivot = value) 在底层调用 Cholesky(x, perm = value, ...) 。如果除了 Cholesky 因子 L' 之外，您还必须知道排列 P_{1} ，则直接调用 Cholesky ，因为 chol(x, pivot = TRUE) 的结果指定 L' 但不指定 P_{1} 。

值

表示上三角乔列斯基因子 L' 的矩阵 triangularMatrix 或 diagonalMatrix 。如果 x 是传统矩阵，则结果是传统矩阵；如果 x 是稠密的，结果是稠密的；如果 x 是稀疏的，结果是稀疏的。

例子

showMethods("chol", inherited = FALSE)
set.seed(0)

## ---- Dense ----------------------------------------------------------

## chol(x, pivot = value) wrapping Cholesky(x, perm = value)
selectMethod("chol", "dsyMatrix")

## Except in packed cases where pivoting is not yet available
selectMethod("chol", "dspMatrix")

## .... Positive definite ..............................................

(A1 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 5)))
(R1.nopivot <- chol(A1))
(R1 <- chol(A1, pivot = TRUE))

## In 2-by-2 cases, we know that the permutation is 1:2 or 2:1,
## even if in general 'chol' does not say ...

stopifnot(exprs = {
   all.equal(  A1           , as(crossprod(R1.nopivot), "dsyMatrix"))
   all.equal(t(A1[2:1, 2:1]), as(crossprod(R1        ), "dsyMatrix"))
   identical(Cholesky(A1)@perm, 2:1) # because 5 > 1
})

## .... Positive semidefinite but not positive definite ................

(A2 <- new("dpoMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 4)))
try(R2.nopivot <- chol(A2)) # fails as not positive definite
(R2 <- chol(A2, pivot = TRUE)) # returns, with a warning and ...

stopifnot(exprs = {
   all.equal(t(A2[2:1, 2:1]), as(crossprod(R2), "dsyMatrix"))
   identical(Cholesky(A2)@perm, 2:1) # because 4 > 1
})

## .... Not positive semidefinite ......................................

(A3 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 3)))
try(R3.nopivot <- chol(A3)) # fails as not positive definite
(R3 <- chol(A3, pivot = TRUE)) # returns, with a warning and ...

## _Not_ equal: see details and examples in help("Cholesky")
all.equal(t(A3[2:1, 2:1]), as(crossprod(R3), "dsyMatrix"))

## ---- Sparse ---------------------------------------------------------

## chol(x, pivot = value) wrapping
## Cholesky(x, perm = value, LDL = FALSE, super = FALSE)
selectMethod("chol", "dsCMatrix")

## Except in diagonal cases which are handled "directly"
selectMethod("chol", "ddiMatrix")

(A4 <- toeplitz(as(c(10, 0, 1, 0, 3), "sparseVector")))
(ch.A4.nopivot <- Cholesky(A4, perm = FALSE, LDL = FALSE, super = FALSE))
(ch.A4 <- Cholesky(A4, perm = TRUE, LDL = FALSE, super = FALSE))
(R4.nopivot <- chol(A4))
(R4 <- chol(A4, pivot = TRUE))

det4 <- det(A4)
b4 <- rnorm(5L)
x4 <- solve(A4, b4)

stopifnot(exprs = {
    identical(R4.nopivot, expand1(ch.A4.nopivot, "L."))
    identical(R4, expand1(ch.A4, "L."))
    all.equal(A4, crossprod(R4.nopivot))
    all.equal(A4[ch.A4@perm + 1L, ch.A4@perm + 1L], crossprod(R4))
    all.equal(diag(R4.nopivot), sqrt(diag(ch.A4.nopivot)))
    all.equal(diag(R4), sqrt(diag(ch.A4)))
    all.equal(sqrt(det4), det(R4.nopivot))
    all.equal(sqrt(det4), det(R4))
    all.equal(det4, det(ch.A4.nopivot, sqrt = FALSE))
    all.equal(det4, det(ch.A4, sqrt = FALSE))
    all.equal(x4, solve(R4.nopivot, solve(t(R4.nopivot), b4)))
    all.equal(x4, solve(ch.A4.nopivot, b4))
    all.equal(x4, solve(ch.A4, b4))
})

参考

The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dpstrf.f, https://netlib.org/lapack/double/dpotrf.f, and https://netlib.org/lapack/double/dpptrf.f.

The CHOLMOD source code; see https://github.com/DrTimothyAldenDavis/SuiteSparse, notably the header file ‘CHOLMOD/Include/cholmod.h’ defining cholmod_factor_struct.

Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Transactions on Mathematical Software, 35(3), Article 22, 1-14. doi:10.1145/1391989.1391995

Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004). Algorithm 837: AMD, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software, 17(4), 886-905. doi:10.1145/1024074.1024081

Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944

也可以看看

base 和 chol 中的默认方法调用传统矩阵 x 。

通用函数 Cholesky ，尤其是在计算 Cholesky 分解而不仅仅是因子 L' 时具有更大的灵活性。