[R] How to solve A'A=S for A
ripley@stats.ox.ac.uk
ripley at stats.ox.ac.uk
Tue Feb 18 20:21:02 CET 2003
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
More information about the R-help
mailing list