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' 時具有更大的靈活性。