On 7/16/07, Kevin Jacobs <[EMAIL PROTECTED]> <[EMAIL PROTECTED]>
wrote:
On 7/16/07, Charles R Harris <[EMAIL PROTECTED]> wrote:
>
>
>
> On 7/16/07, Robert Kern < [EMAIL PROTECTED]> wrote:
> >
> > Kevin Jacobs <[EMAIL PROTECTED]> 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
_______________________________________________
Numpy-discussion mailing list
[email protected]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
