[R] How to solve A'A=S for A
Douglas Bates
bates at stat.wisc.edu
Tue Feb 18 20:25:04 CET 2003
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/
More information about the R-help
mailing list