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

[ripley@stats.ox.ac.uk]

It's usual to use R for an upper-triangular matrix, not a lower-triangular 
one. (As in the QR decomposition.)

If U'U = W then Winv = Uinv Uinv' and so you can take R = t(Uinv), that is
> t(solve(chol(W)))
     [,1] [,2] [,3]
[1,]    1  0.0    0
[2,]    0  0.5    0
[3,]   -2 -0.5    1

Now, you can do that a bit more efficiently by

t(backsolve(chol(W), diag(ncol(W))))

since triangular matrices are relatively simple to invert.  But R is not 
going to notice unless W is enormous.


On Tue, 18 Feb 2003, 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
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595