[Numpy-discussion] NumPy/SciPy LAPACK version

Charles R Harris charlesr.harris@gmail....
Mon Jul 16 12:00:45 CDT 2007


On 7/16/07, Robert Kern <robert.kern@gmail.com> wrote:
>
> Kevin Jacobs <jacobs@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 [2]: sqrtm(eye(2))
Out[2]:
array([[ 1.,  0.],
       [ 0.,  1.]])


Perhaps we aren't using the best method here.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/numpy-discussion/attachments/20070716/6905d594/attachment.html 


More information about the Numpy-discussion mailing list