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
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion