[Numpy-discussion] NumPy/SciPy LAPACK version

Kevin Jacobs <jacobs@bioinformed.com> bioinformed@gmail....
Mon Jul 16 11:24:52 CDT 2007


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.

-Kevin


On 7/16/07, Hanno Klemm <klemm@phys.ethz.ch> wrote:
>
>
> Kevin,
>
> the problem appears to be that sqrtm() gives back an array, rather
> than a matrix:
>
>
> >>> import scipy.linalg as sl
> >>> a = s.matrix([[59, 64, 69],[64, 72, 80],[69, 80, 91]])
> >>> type(a)
> <class 'numpy.core.defmatrix.matrix'>
> >>> a
> matrix([[59, 64, 69],
>         [64, 72, 80],
>         [69, 80, 91]])
> >>> a*a - N.dot(a,a)
> matrix([[0, 0, 0],
>         [0, 0, 0],
>         [0, 0, 0]])
> >>> b = sl.sqrtm(a)
> >>> type(b)
> <type 'numpy.ndarray'>
> >>> N.dot(b,b)
> array([[ 59. +1.85288457e-22j,  64. -6.61744490e-23j,  69.
> +1.85288457e-22j],
>        [ 64. -2.64697796e-23j,  72. -3.70576914e-22j,  80.
> -5.29395592e-23j],
>        [ 69. +1.85288457e-22j,  80. -1.32348898e-22j,  91.
> +2.38228016e-22j]])
> >>> b*b - N.dot(b,b)
> array([[-33.91512711 +1.03595922e-07j, -44.57413548 -1.82329475e-07j,
>         -54.51073741 +7.87335527e-08j],
>        [-44.57413548 -1.82329475e-07j, -48.15108664 +4.04046172e-07j,
>         -51.27477788 -2.21716697e-07j],
>        [-54.51073741 +7.87335527e-08j, -51.27477788 -2.21716697e-07j,
>         -43.21448471 +1.42983144e-07j]])
> >>>
>
> This certainly is a slightly unexpected behaviour.
>
> Hanno
>
>
> "Kevin Jacobs <jacobs@bioinformed.com>" <bioinformed@gmail.com> said:
>
> > ------=_Part_59405_32758974.1184593945795
> > Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> > Content-Transfer-Encoding: 7bit
> > Content-Disposition: inline
> >
> > Hi all,
> >
> > This is a bit of a SciPy question, but I thought I'd ask here since I'm
> > already subscribed.  I'd like to add some new LAPACK bindings to
> SciPy and
> > was wondering if there was a minimum version requirement for LAPACK,
> since
> > it would be ideal if I could use some of the newer 3.0 features.  In
> > addition to using some block methods only added in 3.0, it is very
> > convenient to use the WORK=-1 for space queries instead of
> reimplementing
> > the underlying logic in the calc_work module.
> >
> > The routines of most interest to me are:
> > DGELSD
> > DGGGLM
> > DGGLSE
> >
> > I've also found that SciPy's sqrtm is broken:
> >
> > >>> a=matrix([[59, 64, 69],
> >         [64, 72, 80],
> >         [69, 80, 91]])
> > >>> sqrtm(a)
> > array([[ 5.0084801  +1.03420519e-08j,  4.40747825 -2.06841037e-08j,
> >          3.8064764  +1.03420519e-08j],
> >        [ 4.40747825 -2.06841037e-08j,  4.88353492 +4.13682075e-08j,
> >          5.3595916  -2.06841037e-08j],
> >        [ 3.8064764  +1.03420519e-08j,  5.3595916  -2.06841037e-08j,
> >          6.9127068  +1.03420519e-08j]])
> > >>> sqrtm(a)*sqrtm(a)
> > array([[ 25.08487289 +1.03595922e-07j,  19.42586452 -1.82329475e-07j,
> >          14.48926259 +7.87335527e-08j],
> >        [ 19.42586452 -1.82329475e-07j,  23.84891336 +4.04046172e-07j,
> >          28.72522212 -2.21716697e-07j],
> >        [ 14.48926259 +7.87335527e-08j,  28.72522212 -2.21716697e-07j,
> >          47.78551529 +1.42983144e-07j]])
> >
> > So not even close...  (and yes, it does deserve a bug report if one
> doesn't
> > already exist)
> >
> > -Kevin
> >
> > ------=_Part_59405_32758974.1184593945795
> > Content-Type: text/html; charset=ISO-8859-1
> > Content-Transfer-Encoding: 7bit
> > Content-Disposition: inline
> >
> > Hi all,<br><br>This is a bit of a SciPy question, but I thought
> I&#39;d ask here since I&#39;m already subscribed.&nbsp; I&#39;d like
> to add some new LAPACK bindings to SciPy and was wondering if there
> was a minimum version requirement for LAPACK, since it would be ideal
> if I could use some of the newer
> > 3.0 features.&nbsp; In addition to using some block methods only
> added in 3.0, it is very convenient to use the WORK=-1 for space
> queries instead of reimplementing the underlying logic in the
> calc_work module.<br><br>The routines of most interest to me are:
> > <br>DGELSD<br>DGGGLM<br>DGGLSE<br><br>I&#39;ve also found that
> SciPy&#39;s sqrtm is broken:<br><br>&gt;&gt;&gt; a=matrix([[59, 64,
> 69],<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [64, 72,
> 80],<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [69, 80,
> 91]])<br>&gt;&gt;&gt; sqrtm(a)<br>array([[
> > 5.0084801&nbsp; +1.03420519e-08j,&nbsp; 4.40747825
> -2.06841037e-08j,<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
> 3.8064764&nbsp;
> +1.03420519e-08j],<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [
> 4.40747825 -2.06841037e-08j,&nbsp; 4.88353492
> +4.13682075e-08j,<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
> 5.3595916&nbsp;
> -2.06841037e-08j],<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [
> > 3.8064764&nbsp; +1.03420519e-08j,&nbsp; 5.3595916&nbsp;
> -2.06841037e-08j,<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
> 6.9127068&nbsp; +1.03420519e-08j]])<br>&gt;&gt;&gt;
> sqrtm(a)*sqrtm(a)<br>array([[ 25.08487289 +1.03595922e-07j,&nbsp;
> 19.42586452
> -1.82329475e-07j,<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
> > 14.48926259
> +7.87335527e-08j],<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [
> 19.42586452 -1.82329475e-07j,&nbsp; 23.84891336
> +4.04046172e-07j,<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
> 28.72522212 -2.21716697e-07j],<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
> [ 14.48926259 +7.87335527e-08j,&nbsp; 28.72522212 -2.21716697e-07j,<br>
> > &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 47.78551529
> +1.42983144e-07j]])<br><br>So not even close...&nbsp; (and yes, it
> does deserve a bug report if one doesn&#39;t already
> exist)<br><br>-Kevin<br><br>
> >
> > ------=_Part_59405_32758974.1184593945795--
> >
>
>
>
> --
> Hanno Klemm
> klemm@phys.ethz.ch
>
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/numpy-discussion/attachments/20070716/5ff3954a/attachment-0001.html 


More information about the Numpy-discussion mailing list