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 <[EMAIL PROTECTED]>" <[EMAIL PROTECTED]> 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'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.<br><br>The routines of most interest to me are:
> <br>DGELSD<br>DGGGLM<br>DGGLSE<br><br>I've also found that
SciPy's sqrtm is broken:<br><br>>>> a=matrix([[59, 64,
69],<br> [64, 72,
80],<br> [69, 80,
91]])<br>>>> sqrtm(a)<br>array([[
> 5.0084801 +1.03420519e-08j, 4.40747825
-2.06841037e-08j,<br>
3.8064764
+1.03420519e-08j],<br> [
4.40747825 -2.06841037e-08j, 4.88353492
+4.13682075e-08j,<br>
5.3595916
-2.06841037e-08j],<br> [
> 3.8064764 +1.03420519e-08j, 5.3595916
-2.06841037e-08j,<br>
6.9127068 +1.03420519e-08j]])<br>>>>
sqrtm(a)*sqrtm(a)<br>array([[ 25.08487289 +1.03595922e-07j,
19.42586452
-1.82329475e-07j,<br>
> 14.48926259
+7.87335527e-08j],<br> [
19.42586452 -1.82329475e-07j, 23.84891336
+4.04046172e-07j,<br>
28.72522212 -2.21716697e-07j],<br>
[ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j,<br>
> 47.78551529
+1.42983144e-07j]])<br><br>So not even close... (and yes, it
does deserve a bug report if one doesn't already
exist)<br><br>-Kevin<br><br>
>
> ------=_Part_59405_32758974.1184593945795--
>
--
Hanno Klemm
[EMAIL PROTECTED]
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
