[Numpy-discussion] NumPy/SciPy LAPACK version

Hanno Klemm klemm@phys.ethz...
Mon Jul 16 09:48:06 CDT 2007


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




More information about the Numpy-discussion mailing list