[R] How to solve A'A=S for A

[Douglas Bates]

Ralf Engelhorn <rengelho at ix.urz.uni-heidelberg.de> writes:

> Dear R-helpers,
> thank you very much for your immediate and numerous reactions
> (also via E-mail).
> I see that the problem I encountered lies in the field of
> matrix algebra and not R. Beeing a biologist, I think I should
> direct further questions to a local statistician.
> 
> To give you some feedback on your engagement, I'd like to report,
> what I was struggling with: Following a numerical example
> from Burnaby 1966 (Growth-invariant discriminant functions and
> generalized distances. Biometrics 22:96-110) one of the first
> steps is to solve R'R=W^{-1} for R (where all matrices have p x p
> order and W^{-1} is the inverse of the within-group covariance
> matrix).
> 
> When I tried
> 
> W <- matrix(c(1,0,2,0,4,2,2,2,6),3,3)
> W.inv <- solve(W)
> A.chol <- chol(W.inv);A.chol
>          [,1]      [,2]       [,3]
> [1,] 2.236068 0.4472136 -0.8944272
> [2,] 0.000000 0.5477226 -0.1825742
> [3,] 0.000000 0.0000000  0.4082483
> 
> I received a matrix A.chol which was different from R given
> in Burnaby's example
> 
> R <- matrix(c(1,0,-2,0,0.5,-0.5,0,0,1),3,3);R
>      [,1] [,2] [,3]
> [1,]    1  0.0    0
> [2,]    0  0.5    0
> [3,]   -2 -0.5    1
> 
> The guidance the author provides for finding R is:
> "The matrix R may be computed during the process of inverting W (it is not
> really necessary to complete the inversion)."
> The only way that I know to compute the inverse of a matrix is via the
> matrix of cofactors (Searle 1982). But this matrix
> did not resemble the R that I was looking for.

The author is referring to the fact that you can determine the inverse
of a positive-definite, symmetric matrix by computing its Cholesky
decomposition W = R'R and inverting only the triangular matrix R.
This is because W^{-1} = R^{-1} R^{-1}' so a lower triangular factor
for W^{-1}, in the sense of L where W^{-1}=L'L, is R^{-1}'.  You can
reproduce this in R as

> W <- matrix(c(1,0,2,0,4,2,2,2,6),3,3)
> R = La.chol(W)
> R
     [,1] [,2] [,3]
[1,]    1    0    2
[2,]    0    2    1
[3,]    0    0    1
> L = t(solve(R))
> t(L) %*% L
     [,1] [,2] [,3]
[1,]    5  1.0 -2.0
[2,]    1  0.5 -0.5
[3,]   -2 -0.5  1.0
> crossprod(L)    # same as above
     [,1] [,2] [,3]
[1,]    5  1.0 -2.0
[2,]    1  0.5 -0.5
[3,]   -2 -0.5  1.0
> round(W %*% crossprod(L), 6) # shows t(L) %*% L is W^{-1}
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/