[Numpy-discussion] NumPy/SciPy LAPACK version

[Kevin Jacobs <jacobs@bioinformed.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 [2]: sqrtm(eye(2))
> Out[2]:
> 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,
         -5.55111512e-16+0.j],
       [ -1.05471187e-15+0.j,  -5.55111512e-17+0.j,   5.55111512e-17+0.j,
          0.00000000e+00+0.j],
       [ -6.66133815e-16+0.j,   1.11022302e-16+0.j,   1.66533454e-16+0.j,
          1.11022302e-16+0.j],
       [ -5.55111512e-16+0.j,   1.11022302e-16+0.j,   1.38777878e-16+0.j,
          8.32667268e-17+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.

-Kevin
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070716/b9513f2d/attachment.html>