[R] Matrix question: obtaining the square root of a positive definite matrix?

[Peter Dalgaard]

Prof Brian Ripley wrote:
> On Wed, 24 Jan 2007, gallon li wrote:
>
>   
>> I want to compute B=A^{1/2} such that B*B=A.
>>     
>
> According to your subject line A is positive definite and hence 
> symmetric?  The usual definition of a matrix square root involves a 
> transpose, e.g. B'B = A.  There are many square roots: were you looking 
> for a symmetric one?
>
>   
If not, Choleski decomposition by chol() is often the expedient way.

> For such an A,
>
>   
>> e <- eigen(A)
>> V <- e$vectors
>> V %*% diag(e$values) %*% t(V)
>>     
>
> recovers A (up to rounding errors), and
>
>   
>> B <- V %*% diag(sqrt(e$values)) %*% t(V)
>>     
>
> is such that B %*% B = A.  Even that is not unique, e.g. -B is an equally 
> good answer.
>
>
>   
and you can flip the sign of the individual root eigenvalues too, and if
the eigenvalues are not unique, you can rotate the eigenspace coordinate
systems at will and then flip signs.

-- 
   O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907