2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Cholesky decomposition is such that the positive definite matrix $K = LL^T$. $L$ is lower-triangular. $\endgroup$ – usεr11852 Mar 22 '18 at 22:29
  • $\begingroup$ Try what crossprod(A) returns. $\endgroup$ – usεr11852 Mar 22 '18 at 22:34
  • $\begingroup$ crossprod(A) returns S (the expected answer). So, insofar as the link on Wikipedia, should I use A or A' as returned by chol? $\endgroup$ – Jane Wayne Mar 22 '18 at 22:35
  • $\begingroup$ Cool. Now check what computation crossprod is doing. $\endgroup$ – usεr11852 Mar 22 '18 at 22:40
  • 3
    $\begingroup$ 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.) $\endgroup$ – usεr11852 Mar 22 '18 at 22:49
2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Jane Wayne Mar 22 '18 at 22:45
  • $\begingroup$ 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 $\endgroup$ – Aksakal Mar 23 '18 at 1:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.