[Numpy-discussion] expm

[Nils Wagner]

On Fri, 20 Jul 2007 13:03:09 -0400
  "Kevin Jacobs <jacobs at bioinformed.com>" 
<bioinformed at gmail.com> wrote:
> On 7/20/07, Anne Archibald <peridot.faceted at gmail.com> 
>wrote:
>>
>> On 20/07/07, Nils Wagner <nwagner at iam.uni-stuttgart.de> 
>>wrote:
>> > lorenzo bolla wrote:
>> > > hi all.
>> > > is there a function in numpy to compute the exp of a 
>>matrix, similar
>> > > to expm in matlab?
>> > > for example:
>> > > expm([[0,0],[0,0]]) = eye(2)
>> > Numpy doesn't provide expm but scipy does.
>> > >>> from scipy.linalg import expm, expm2, expm3
>>
>> Just as a warning, numpy does provide expm1, but it does 
>>something
>> different (exp(x)-1, computed directly).
>>
> 
> On a separate note, I'm working to provide faster and 
>more accurate versions
> of sqrtm and expm.  The current versions do not take 
>full advantage of
> LAPACK.  Here are some preliminary benchmarks:
> 
> Ill-conditioned
> ----------------
> linalg.sqrtm   : error=9.37e-27, 573.38 usec/pass
> sqrtm_svd      : error=2.16e-28, 142.38 usec/pass
> sqrtm_eig      : error=4.79e-27, 270.38 usec/pass
> sqrtm_symmetric: error=1.04e-27, 239.30 usec/pass
> sqrtm_symmetric2: error=2.73e-27, 190.03 usec/pass
> 
> Well-conditioned
> ----------------
> linalg.sqrtm   : error=1.83e-29, 478.67 usec/pass
> sqrtm_svd      : error=8.11e-30, 130.57 usec/pass
> sqrtm_eig      : error=4.50e-30, 255.56 usec/pass
> sqrtm_symmetric: error=2.78e-30, 237.61 usec/pass
> sqrtm_symmetric2: error=3.35e-30, 167.27 usec/pass
> 
> Large
> ----------------
> linalg.sqrtm   : error=5.95e-25, 8450081.68 usec/pass
> sqrtm_svd      : error=1.64e-24, 151206.61 usec/pass
> sqrtm_eig      : error=6.31e-24, 549837.40 usec/pass
> sqrtm_symmetric: error=8.55e-25, 177422.29 usec/pass
> 
> where:
> 
> def sqrtm_svd(x):
>  u,s,vt = linalg.svd(x)
>  return dot(u,transpose((s**0.5)*transpose(vt)))
> 
> def sqrtm_eig(x):
>  d,e = linalg.eig(x)
>  d = (d**0.5).astype(float)
>  return dot(e,transpose(d*e))
> 
> def sqrtm_symmetric(x,cond=1e-7):
>  d,e = linalg.eigh(x)
>  d[d<cond] = 0
>  return dot(e,transpose((d**0.5)*e)).astype(float)
> 
> def sqrtm_symmetric2(x):
>  # Not as robust due to initial Cholesky step
>  l=linalg.cholesky(x,lower=1)
>  u,s,vt = linalg.svd(l)
>  return dot(u,transpose(s*u))
> 
> with SciPy linked against ACML.
> 
> -Kevin

  
Kevin,

Your sqrtm_eig(x) function won't work if x is defective.
See test_defective.py for details.

Have you considered the algorithms proposed by
Nick Higham for various matrix functions ?

http://www.maths.manchester.ac.uk/~higham/pap-mf.html

Nils
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_defective.py
Type: text/x-python
Size: 467 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070720/f951a4fd/attachment.py>