Why does the resulting matrix from Cholesky decomposition of a covariance matrix when multiplied by its transpose not give back the covariance matrix?

Question

I have a covariance matrix, S, which I use Cholesky decomposition to find A. It is stated that AA'=S, however, I am not recovering S when I do AA'. The example code in R is as follows.

S <- matrix(c(1.091385, 1.949606, 1.949606, 4.520746), 2, 2)
A <- chol(S)
T <- t(A)
R <- A %*% T

As can be seen, S is as follows.

1.091385, 1.949606
1.949606, 4.520746

But R is as follows.

4.574082, 1.901370
1.901370, 1.038049

And for completeness, A is as follows.

1.044694, 1.866199
0.000000, 1.018847

However, when I do A'A with the code R <- T %*% A, I do recover S.

Any ideas on what I'm doing wrong? Or is the link on Wikipedia wrong? The R examples on Cholesky decomposition seems to suggest A'A=S.

Cholesky decomposition is such that the positive definite matrix $K = LL^T$. $L$ is *lower*-triangular. — usεr11852, Mar 22 '18 at 22:29
`crossprod(A)` returns `S` (the expected answer). So, insofar as the link on Wikipedia, should I use `A` or `A'` as returned by `chol`? — Jane Wayne, Mar 22 '18 at 22:35
The inconvenient truth is that the Cholesky decomposition while usually defined as $K=LL^T$ where $L$ is lower triangular is equally valid as $K=U^TU$ where $U$ is upper triangular. The implementation of Cholesky decomposition in LAPACK (the libraries our computer use to compute Linear Algebra tasks) allow both expressions. R unfortunately has [hard-coded](https://github.com/wch/r-source/blob/3293d437efa2843ef3a74ad53ca4087d5d4f48b0/src/modules/lapack/Lapack.c) the upper one. (There is a `U` in the call of the routine `dpstrf` that actually compute the Cholesky.) — usεr11852, Mar 22 '18 at 22:49

score 11 · Accepted Answer · answered Mar 22 '18 at 23:31

As explained in my comment, the inconvenient truth is that the Cholesky decomposition while usually defined as $K=LL^T$ where $L$ is lower triangular, is equally valid as $K=U^TU$ where $U$ is upper triangular. The implementation of Cholesky decomposition in LAPACK (the libraries our computer use to compute Linear Algebra tasks) allow both expressions. R unfortunately has hard-coded the upper one. (There is a U in the call of the routine dpstrf that actually compute the Cholesky.)

This means one has to transpose the results for chol in order to get a lower triangular matrix. After that is done, as you have already discovered yourself, the result follow-up directly. So for example, given the matrix S of the original post:

U <- chol(S);
L <- t(chol(S));
S -  crossprod(U) # This is equivalent to S- U^T*U and should be approx. 0
S - tcrossprod(L) # This is equivalent to S- L*L^T and should be approx. 0

I hope it is therefore clear that the Wikipedia page is not wrong. Being somewhat critical, Wikipedia's Cholesky decomposition article should probably mention that the Cholesky decomposition $K=LL^T$ is equivalent to $K=U^TU$ where $U=L^T$. It was probably omitted for consistency of notation.

Can the person who downvoted this post please give some feedback about why they considered this a bad answer? — usεr11852, Feb 03 '20 at 21:06

Aksakal · Answer 2 · 2018-03-22T22:33:43.657

1

So, why do you think that chol(S) returns your A and not A'?

In fact it does return A' if you look at the values or read the documentation. It returns upper triangular, which corresponds to A' in your Wiki reference

edited Mar 22 '18 at 22:33

answered Mar 22 '18 at 22:30

Aksakal

55,939
5
90
176

I would say that the Wikipedia article is not clear on whether `A` is upper or lower matrix. Additionally, if `chol` returns `A'` and that matrix is the same as `A'` on the Wikipedia article, then `AA'` should have recovered `S`. The way I'm reading all of this information is that `A'` returned from `chol` is the upper matrix which is `A` in the Wikipedia article. – Jane Wayne Mar 22 '18 at 22:45
Your code `T %*% A` demonstrates my point. You got A' from chol() function, so you need to get A=t(A'), which you called T. That's why you're able to recover R with TA. I understand that it's a bit confusing but when you printed A you could see immediately that it's upper triangle. It would be a little less confusing if you called `U=chol(R)`, then you wouldn't have a weird situation where your A is really A' in the wiki page that you refer to – Aksakal Mar 23 '18 at 01:03

Why does the resulting matrix from Cholesky decomposition of a covariance matrix when multiplied by its transpose not give back the covariance matrix?

2 Answers2