# [Numpy-discussion] NumPy/SciPy LAPACK version

Charles R Harris charlesr.harris@gmail....
Mon Jul 16 12:37:54 CDT 2007

```On 7/16/07, Kevin Jacobs <jacobs@bioinformed.com> <bioinformed@gmail.com>
wrote:
>
> On 7/16/07, Charles R Harris <charlesr.harris@gmail.com> wrote:
> >
> >
> >
> > 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.
> >
>
>
> 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://projects.scipy.org/pipermail/numpy-discussion/attachments/20070716/75e384e9/attachment.html
```