I assumed that chol2inv(chol(x))
and solve(x)
are two different methods that arrive at the same conclusion in all cases. Consider for instance a matrix S
A <- matrix(rnorm(3*3), 3, 3)
S <- t(A) %*% A
where the following two commands will give equivalent results:
solve(S)
chol2inv(chol(S))
Now consider the transpose of the Cholesky decomposition of S:
L <- t(chol(S))
where now the results of the following two commands do not give equivalent results anymore:
solve(L)
chol2inv(chol(L))
This surprised me a bit. Is this expected behavior?
CodePudding user response:
chol
expects (without checking) that its first argument x
is a symmetric positive definite matrix, and it operates only on the upper triangular part of x
. Thus, if L
is a lower triangular matrix and D = diag(diag(L))
is its diagonal part, then chol(L)
is actually equivalent to chol(D)
, and chol2inv(chol(L))
is actually equivalent to solve(D)
.
set.seed(141339L)
n <- 3L
S <- crossprod(matrix(rnorm(n * n), n, n))
L <- t(chol(S))
D <- diag(diag(L))
all.equal(chol(L), chol(D)) # TRUE
all.equal(chol2inv(chol(L)), solve(D)) # TRUE