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.
-Kevin
_______________________________________________
Numpy-discussion mailing list
[email protected]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
