[R] matrix^(-1/2)

David Winsemius dwinsemius at comcast.net
Mon Nov 2 00:52:53 CET 2009 

On Nov 1, 2009, at 3:12 PM, Charles C. Berry wrote:

> On Sun, 1 Nov 2009, David Winsemius wrote:
>
>>
>> On Nov 1, 2009, at 1:46 PM, spencerg wrote:
>>
>>> Hi, Chuck:
>>>
>>>    Thanks very much, but why do I get "package 'expm' is not  
>>> available"
>>>    from install.packages("expm",repos="http://R-Forge.R- 
>>> project.org")?
>>
>> In my case I think it was it is because there is no 2.10 branch to  
>> either the:

snipped details

>>
>> So I wonder if the package installers' expectations for the r-forge  
>> repository are matching up with the tree structures.
>
> Right. FWIW, the source install works OK on my linux box:
>
>> sessionInfo()
> R version 2.10.0 (2009-10-26)
> x86_64-pc-linux-gnu

Yes, as I said, I had already succeeded in a source compilation using:
install.packages("expm",repos="http://R-Forge.R-project.org",  
type="source")

(I am not quite sure why R was able to navigate the tree, but suspect  
my issues may stem from generally using the 64-bit Mac version of R.)

>
>>
>> I should also note that the matpow or "%^%" functions in expm would  
>> not address the OP's question since they require that the exponent  
>> be positive.
>>
>
> Roger that.
>
> If solve(chol(A)) isn't good enough a symmetric inverse square root  
> is available from expm as 'solve( sqrtm( A ) )'

I (as a noobisher matrix mechanic)  have discovered that the symmetric  
part is essential. Experiments with non-symmetric examples have come  
to a "singularly bad end". The OP did offer a symmetric example, so he  
probably was more matrix-savvy than I.  I am not sufficiently  
knowledgeable to design the error checking. I think a useful addition  
to %^% would be properly designed negative fractional matrix powers  
with more informative error messages for the matrix klutzes among us.  
I hope I am correct in thinking that M %^% (-1/integer) has meaning  
for integers other than 2, correct?  And noticing that the current  
version of %^% rounds the exponent, I think that raises the question  
(given that the matrix argument were symmetric), would one round the 1/ 
<negative-exponent>? And if so, Up or down?

-- 
David.

>
> Chuck
>
>
>> -- 
>> David.
>>
>>>
>>>    Best Wishes,
>>>    Spencer Graves
>>> Charles C. Berry wrote:
>>> > On Sun, 1 Nov 2009, spencerg wrote:
>>> > > > A question, a comment, and an alternative answer to  
>>> matrix^(-1/2):
>>> > > > > QUESTION:
>>> > > > > > > What's the status of the "expm" package, mentioned in  
>>> the email you > > cited from Martin Maechler, dated Apr 5 19:52:09  
>>> CEST 2008? I tried > > both install.packages('expm') and > >  
>>> install.packages("expm",repos="http://R-Forge.R-project.org"), and  
>>> got > > "package 'expm' is not available" in both cases.
>>> > > > > > Try
>>> > >    http://r-forge.r-project.org/projects/expm/
>>> > > HTH,
>>> > > Chuck
>>> > > > > > COMMENT:
>>> > > > > > > The solution proposed by Venables rests on Sylvester's  
>>> matrix theorem, > > which essentially says that if a matrix A is  
>>> diagonalizable with > > eigenvalue decomposition eigA <- eigen(A)  
>>> and f: D → C with D ⊂ C > > be a function for which f(A) is  
>>> well defined > > (http://en.wikipedia.org/wiki/Sylvester%27s_matrix_theorem 
>>> ), then f(A) > > = with(eigA, vectors %*% diag(f(values)) %*%  
>>> solve(vectors)). Maechler > > and others have noted that this can  
>>> be one of the least accurate and > > most computationally  
>>> expensive ways to compute f(A).
>>> > > > > > > ALTERNATIVE ANSWER:
>>> > > > > > > For A^(-1/2), if A is symmetric and nonnegative  
>>> definite, then > > solve(chol(A)) would be a very good way to  
>>> compute it.
>>> > > > > > > Hope this helps,
>>> > > Spencer
>>> > > > > > > David Winsemius wrote:
>>> > > > > > > On Oct 31, 2009, at 9:33 PM, David Winsemius wrote:
>>> > > > > > > > >   On Oct 31, 2009, at 4:39 PM, Kajan Saied wrote:
>>> > > > > > >   Dear R-Help Team,
>>> > > > > > > >   as a R novice I have a (maybe for you very simple  
>>> > > > > > > >   question), how do I > >  get
>>> > > > > >   the following solved in R:
>>> > > > > > > >   Let R be a n x n matrix:
>>> > > > > > > >   \mid R\mid^{-\frac{1}{2}}
>>> > > > > > > >   solve(A) gives me the inverse of the matrix R,  
>>> however not > > > > > > >   the ^(-1/2) > >  of
>>> > > > > >   the matrix...
>>> > > > > >   GIYF: (and Bill Venables if friendly, too.)
>>> > > > > >   http://www.lmgtfy.com/?q=powers+of+matrix+r-project
>>> > > > > > > I had assumed that the first hit I got:
>>> > > > > > > https://stat.ethz.ch/pipermail/r-help/2008-April/160662.html
>>> > > > > > > ... would be the first hit anybody got, but that's not  
>>> necessarily > > > true
>>> > > > now and especially for the future. And further searching  
>>> within the
>>> > > > results produced this more recent Maechler posting:
>>> > > > > > > https://stat.ethz.ch/pipermail/r-devel/2008-April/048969.html
>>> > > > > > > For the Mac users, there appears to be no binary, but  
>>> the source > > > compiles
>>> > > > without error on a 64-bit version of R 2.10.0:
>>> > > > > > > install.packages("expm",repos="http://R-Forge.R-project.org 
>>> ",
>>> > > > type="source")
>>> > > > > > > #The suggested code throws an error, so my very minor  
>>> revision would > > > be:
>>> > > > > > > library(expm)
>>> > > > ?"%^%"
>>> > > > > ______________________________________________
>>> > > 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.
>>> > > > > > > Charles C. Berry                            (858)  
>>> 534-2098
>>> >                                            Dept of Family/ 
>>> Preventive > Medicine
>>> > E mailto:cberry at tajo.ucsd.edu                UC San Diego
>>> > http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego  
>>> > 92093-0901
>>
>> David Winsemius, MD
>> Heritage Laboratories
>> West Hartford, CT
>>
>>
>
> Charles C. Berry                            (858) 534-2098
>                                            Dept of Family/Preventive  
> Medicine
> E mailto:cberry at tajo.ucsd.edu	            UC San Diego
> http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego  
> 92093-0901
>

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

More information about the R-help mailing list