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

[Spencer Graves]

There are two square roots for any real number.  With matrices, there 
are infinite.  There are several standard "square root" matrices, one of 
which is the Cholesky decomposition.  I just did "R <- solve(chol(W))", 
and got the transpose of the matrix you reported from Burnaby (1966). 
Since the Cholesky decomposition is not symmetric, it is essentially as 
correct to specify either R or its transpose.

I think it's interesting that solve(chol(W)) is not the same as 
chol(solve(W)), but I'll let someone else comment on that.

Spencer Graves


Ralf Engelhorn wrote:
> 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 symmetric solution (A=A'=KD^{1/2}K', where K has eigenvectors and D
> has eigenvalues) suggested by Jan de Leeuw yields results (discriminant
> function coefficients, D^2) as expected, though a matrix calculated in an
> intermediate step differs from its counterpart in the example.
> 
> To me it looks like having to get deeper into matrix algebra to solve
> this.
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> http://www.stat.math.ethz.ch/mailman/listinfo/r-help