1

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?

yrx1702
  • 1,619
  • 15
  • 27

1 Answers1

2

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
Mikael Jagan
  • 9,012
  • 2
  • 17
  • 48