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.

Answer 1

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.

Answer 2

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