<br><br><div><span class="gmail_quote">On 7/16/07, <b class="gmail_sendername">Kevin Jacobs <<a href="mailto:jacobs@bioinformed.com">jacobs@bioinformed.com</a>></b> <<a href="mailto:bioinformed@gmail.com">bioinformed@gmail.com

</a>> wrote:</span><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div><span class="e" id="q_113d002d9c235707_0"><div><span class="gmail_quote">

On 7/16/07, <b class="gmail_sendername">Charles R Harris</b> <<a href="mailto:charlesr.harris@gmail.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">charlesr.harris@gmail.com</a>> wrote:</span>

<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">

<br><br><div><span><span class="gmail_quote">On 7/16/07, <b class="gmail_sendername">Robert Kern</b> <<a href="mailto:robert.kern@gmail.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">

robert.kern@gmail.com</a>> wrote:</span><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">

Kevin Jacobs <<a href="mailto:jacobs@bioinformed.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">jacobs@bioinformed.com</a>> wrote:<br>> 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.

<br><br>No, in floating point you accumulate error. Those 1e-22j's are part of the<br>actual result. Some systems like MATLAB implicitly silent such small imaginary<br>components; we don't.</blockquote></span><div>

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.

<br><br><span style="font-family: courier new,monospace;">In [2]: sqrtm(eye(2))<br>Out[2]: <br>array([[ 1.,  0.],<br>       [ 0.,  1.]])<br></span><span style="font-family: courier new,monospace;"></span><br style="font-family: courier new,monospace;">

<br style="font-family: courier new,monospace;"></div>Perhaps we aren't using the best method here.<br></div></blockquote></div><br><br></span></div>Here is a better conditioned example:<br><br>>>> a<br>array([[ 1.    ,  

0.5   ,  0.3333,  0.25  ],<br>       [ 0.5   ,  0.3333,  0.25  ,  0.2   ],<br>       [ 0.3333,  0.25  ,  0.2   ,  0.1667],<br>       [ 0.25  ,  0.2   ,  0.1667,  0.1429]])<br>>>> b=sqrtm(a)<br>>>> dot(b,b)

<br>array([[ 1.    +0.j,  0.5   +0.j,  0.3333+0.j,  0.25  +0.j],<br>       [ 0.5   +0.j,  0.3333+0.j,  0.25  +0.j,  0.2   +0.j],<br>       [ 0.3333+0.j,  0.25  +0.j,  0.2   +0.j,  0.1667+0.j],<br>       [ 0.25  +0.j,  0.2

   +0.j,  0.1667+0.j,  0.1429+0.j]])<br>>>> dot(b,b)-a<br>array([[ -1.99840144e-15+0.j,  -9.43689571e-16+0.j,  -5.55111512e-16+0.j,<br>         -5.55111512e-16+0.j],<br>       [ -1.05471187e-15+0.j,  -5.55111512e-17+0.j

,   5.55111512e-17+0.j,<br>          0.00000000e+00+0.j],<br>       [ -6.66133815e-16+0.j,   1.11022302e-16+0.j,   1.66533454e-16+0.j,<br>          1.11022302e-16+0.j],<br>       [ -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.

</blockquote><div> 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.

<br><br>Chuck<br></div><br></div><br>