[Numpy-discussion] NumPy/SciPy LAPACK version
Charles R Harris
charlesr.harris at gmail.com
Mon Jul 16 13:37:54 EDT 2007
On 7/16/07, Kevin Jacobs <jacobs at bioinformed.com> <bioinformed at gmail.com>
> On 7/16/07, Charles R Harris <charlesr.harris at gmail.com> wrote:
> > On 7/16/07, Robert Kern < robert.kern at gmail.com> wrote:
> > >
> > > Kevin Jacobs <jacobs at bioinformed.com> wrote:
> > > > Mea culpa on the msqrt example, however I still think it is wrong to
> > > get
> > > > a complex square-root back when a real valued result is expected and
> > > exists.
> > >
> > > No, in floating point you accumulate error. Those 1e-22j's are part of
> > > the
> > > actual result. Some systems like MATLAB implicitly silent such small
> > > imaginary
> > > components; we don't.
> > The problem is that the given matrix has a conditon number of about
> > 10**17 and is almost singular. A singular value decomposition works fine,
> > but apparently the sqrtm call suffers from roundoff and takes the sqrt of a
> > negative number. Sqrtm returns real results in better conditioned cases.
> > In : sqrtm(eye(2))
> > Out:
> > array([[ 1., 0.],
> > [ 0., 1.]])
> > Perhaps we aren't using the best method here.
> Here is a better conditioned example:
> >>> a
> array([[ 1. , 0.5 , 0.3333, 0.25 ],
> [ 0.5 , 0.3333, 0.25 , 0.2 ],
> [ 0.3333, 0.25 , 0.2 , 0.1667],
> [ 0.25 , 0.2 , 0.1667, 0.1429]])
> >>> b=sqrtm(a)
> >>> dot(b,b)
> array([[ 1. +0.j, 0.5 +0.j, 0.3333+0.j, 0.25 +0.j],
> [ 0.5 +0.j, 0.3333+0.j, 0.25 +0.j, 0.2 +0.j],
> [ 0.3333+0.j, 0.25 +0.j, 0.2 +0.j, 0.1667+0.j],
> [ 0.25 +0.j, 0.2 +0.j, 0.1667+0.j, 0.1429+0.j]])
> >>> dot(b,b)-a
> array([[ -1.99840144e-15+0.j, -9.43689571e-16+0.j, -5.55111512e-16+0.j,
> [ -1.05471187e-15+0.j, -5.55111512e-17+0.j , 5.55111512e-17+0.j,
> [ -6.66133815e-16+0.j, 1.11022302e-16+0.j, 1.66533454e-16+0.j,
> [ -5.55111512e-16+0.j, 1.11022302e-16+0.j , 1.38777878e-16+0.j,
> Also verified the results against svd and eigenvalue methods for computing
> msqrt. I suppose I'm just used to seeing msqrt() implemented completely
> using real valued arithmetic.
I get a real result for this, although the result is wildly incorrect. Sqrtm
isn't part of numpy, where are you getting it from? Mine is coming from
pylab and looks remarkably buggy.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion