R help - Jul 2012 - about changing order of Choleski factorization and inverse operation of a matrix

Dear All,

My question is simple but I need someone to help me out.
Suppose I have a positive definite matrix A.
The funtion chol()  gives matrix L, such that A = L'L.
The inverse of A, say A.inv,  is also positive definite and can be
factorized as A.inv = M'M.
Then
        A = inverse of (A.inv) = inverse of (M'M) = (inverse of M) %*%
(inverse of M)'
           = ((inverse of M)')'%*% (inverse of M)',
i.e. if define B = transpose of (inverse of M), then A = B' %*% B.
Therefore L = B = transpose of (inverse of M) = transpose of (inverse of
chol(A.inv))
But when I try it in R, the answer is not as expected.

code as below:

 A <- matrix(1:9,3,3)
A <- A + t(A)
diag(A) <- 50
print(A)
L <- chol(A)
B <- t(solve(chol(solve(A))))
print(L)
print(B)

Thank you in advance,

Best wishes,
Jie

	[[alternative HTML version deleted]]

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Jie
> Sent: Tuesday, July 31, 2012 2:23 PM
> To: r-help at r-project.org
> Subject: [R] about changing order of Choleski factorization and inverse
> operation of a matrix
> 
> Dear All,
> 
> My question is simple but I need someone to help me out.
> Suppose I have a positive definite matrix A.
> The funtion chol()  gives matrix L, such that A = L'L.
> The inverse of A, say A.inv,  is also positive definite and can be
> factorized as A.inv = M'M.
> Then
>         A = inverse of (A.inv) = inverse of (M'M) = (inverse of M) %*%
> (inverse of M)'
>            = ((inverse of M)')'%*% (inverse of M)',
> i.e. if define B = transpose of (inverse of M), then A = B' %*% B.
> Therefore L = B = transpose of (inverse of M) = transpose of (inverse
> of
> chol(A.inv))
> But when I try it in R, the answer is not as expected.
> 
> code as below:
> 
>  A <- matrix(1:9,3,3)
> A <- A + t(A)
> diag(A) <- 50
> print(A)
> L <- chol(A)
> B <- t(solve(chol(solve(A))))
> print(L)
> print(B)
> 
> Thank you in advance,
> 
> Best wishes,
> Jie
> 
This is not an R question, it is a question about your understanding of matrix
algebra.  Just because

   L'L == B'B

does not mean that L must equal B.  You have your own counter-example. Just take
your calculated L and B matrices and compute

   t(L) %*% L

and

  t(B) %*% B

or

all.equal((t(L) %*% L ) , (t(B) %*% B))


Hope this is helpful,

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204