On Fri, 20 Jul 2007 13:03:09 -0400
"Kevin Jacobs <[EMAIL PROTECTED]>"
<[EMAIL PROTECTED]> wrote:
On 7/20/07, Anne Archibald <[EMAIL PROTECTED]>
wrote:
On 20/07/07, Nils Wagner <[EMAIL PROTECTED]>
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
from scipy import *
def sqrtm_eig(x):
d,e = linalg.eig(x)
d = (d**0.5).astype(float)
return dot(e,transpose(d*e))
def sqrtm_svd(x):
u,s,vt = linalg.svd(x)
return dot(u,transpose((s**0.5)*transpose(vt)))
#
# A is defective
#
A = array((
[4.,1.,1.],
[2.,4.,1.],
[0.,1.,4.]))
X = linalg.sqrtm(A)
Y = sqrtm_eig(A)
Z = sqrtm_svd(A)
res_1 = linalg.norm(A-dot(X,X))
res_2 = linalg.norm(A-dot(Y,Y))
res_3 = linalg.norm(A-dot(Z,Z))
print res_1, res_2, res_3
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
