[Numpy-discussion] NumPy/SciPy LAPACK version

[Charles R Harris]

On 7/16/07, Kevin Jacobs <jacobs at bioinformed.com> <bioinformed at gmail.com>
wrote:
>
> 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.


Hmm,

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.

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