[R] R: Re: R: Re: chol( neg.def.matrix ) WAS: Re: Choleski and Choleski with pivoting of matrix fails

[Ravi Varadhan]

I do not understand what the problem is, as it works just fine for me:

A <- matrix(c(0.5401984,-0.3998675,-1.3785897,-0.3998675,1.0561872,  
0.8158639,-1.3785897, 0.8158639, 1.6073119), 3, 3, byrow=TRUE)

eA <- eigen(A)

chA <-  eA$vec %*% diag(sqrt(eA$val+0i)) %*% t(eA$vec)

all.equal(A, Re(chA %*% t(chA)))

Y <- diag(c(1,2,3))

solve(chA %*% Y)

Ravi.

-----------------------------------------------------------------------------------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.html



------------------------------------------------------------------------------------

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of simona.racioppi at libero.it
Sent: Wednesday, November 25, 2009 9:59 AM
To: P.Dalgaard at biostat.ku.dk
Cc: r-help at r-project.org
Subject: [R] R: Re: R: Re: chol( neg.def.matrix ) WAS: Re: Choleski and Choleski with pivoting of matrix fails

Dear Peter,
thank you very much for your answer.

My problem is that I need to calculate the following quantity:

solve(chol(A)%*%Y)

Y is a 3*3 diagonal matrix and A is a 3*3 matrix. Unfortunately one 
eigenvalue of A is negative. I can anyway take the square root of A but when I 
multiply it by Y, the imaginary part of the square root of A is dropped, and I 
do not get the right answer.

I tried to exploit the diagonal structure of Y by using 2*2 matrices for A 
and Y. In this way the problem mentioned above disappears (since all 
eigenvalues of A are positive) and when I perform the calculation above I get 
approximately the right answer. The approximation is quite good. However it is 
an approximation.

Any suggestion?
Thank you very much!
Simon




>----Messaggio originale----
>Da: P.Dalgaard at biostat.ku.dk
>Data: 23-nov-2009 14.09
>A: "simona.racioppi at libero.it"<simona.racioppi at libero.it>
>Cc: "Charles C. Berry"<cberry at tajo.ucsd.edu>, <r-help at r-project.org>
>Ogg: Re: R: Re: [R] chol( neg.def.matrix ) WAS: Re: Choleski and Choleski 
with pivoting of matrix fails
>
>simona.racioppi at libero.it wrote:
>> It works! But Once I have the square root of this matrix, how do I convert 
it 
>> to a real (not imaginary) matrix which has the same property? Is that 
>> possible?
>
>No. That is theoretically impossible.
>
>If A = B'B, then x'Ax = ||Bx||^2 >= 0
>
>for any x, which implies in particular that all eigenvalues of A should
>be nonnegative.
>
>> 
>> Best,
>> Simon
>> 
>>> ----Messaggio originale----
>>> Da: p.dalgaard at biostat.ku.dk
>>> Data: 21-nov-2009 18.56
>>> A: "Charles C. Berry"<cberry at tajo.ucsd.edu>
>>> Cc: "simona.racioppi at libero.it"<simona.racioppi at libero.it>, <r-help at r-
>> project.org>
>>> Ogg: Re: [R] chol( neg.def.matrix ) WAS: Re: Choleski and Choleski with 
>> pivoting of matrix fails
>>> Charles C. Berry wrote:
>>>> On Sat, 21 Nov 2009, simona.racioppi at libero.it wrote:
>>>>
>>>>> Hi Everyone,
>>>>>
>>>>> I need to take the square root of the following matrix:
>>>>>
>>>>>                [,1]               [,2]                [,3]
>>>>> [1,]  0.5401984 -0.3998675 -1.3785897
>>>>> [2,] -0.3998675  1.0561872  0.8158639
>>>>> [3,] -1.3785897  0.8158639  1.6073119
>>>>>
>>>>> I tried Choleski which fails. I then tried Choleski with pivoting, but
>>>>> unfortunately the square root I get is not valid. I also tried eigen
>>>>> decomposition but i did no get far.
>>>>>
>>>>> Any clue on how to do it?!
>>>>
>>>> If you want to take the square root of a negative definite matrix, you 
>>>> could use
>>>>
>>>>     sqrtm( neg.def.mat )
>>>>
>>>> from the expm package on rforge:
>>>>
>>>>     http://r-forge.r-project.org/projects/expm/
>>> But that matrix is not negative definite! It has 2 positive and one 
>>> negative eigenvalue. It is non-positive definite.
>>>
>>> It is fairly easy in any case to get a matrix square root from the 
eigen 
>>> decomposition:
>>>
>>>> v%*%diag(sqrt(d+0i))%*%t(v)
>>>                       [,1]                  [,2]                  [,3]
>>> [1,]  0.5164499+0.4152591i -0.1247682-0.0562317i -0.7257079+0.3051868i
>>> [2,] -0.1247682-0.0562317i  0.9618445+0.0076145i  0.3469916-0.0413264i
>>> [3,] -0.7257079+0.3051868i  0.3469916-0.0413264i  1.0513849+0.2242912i
>>>> ch <- v%*%diag(sqrt(d+0i))%*%t(v)
>>>> t(ch)%*% ch
>>>               [,1]          [,2]          [,3]
>>> [1,]  0.5401984+0i -0.3998675-0i -1.3785897-0i
>>> [2,] -0.3998675-0i  1.0561872+0i  0.8158639-0i
>>> [3,] -1.3785897-0i  0.8158639-0i  1.6073119-0i
>>>
>>> A triangular square root is, er, more difficult, but hardly impossible.
>>>
>>> -- 
>>>    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
>>>
>> 
>> 
>
>
>-- 
>   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
>
>

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.