[SciPy-dev] Scipy cannot handle singular eigenvalue problems

[Nils Wagner]

Pearu Peterson schrieb:
> 
> On Wed, 10 Dec 2003, Nils Wagner wrote:
> 
> > Hi all,
> >
> > Scipy cannot handle singular eigenvalue problems. Note, that both
> > matrices are singular.
> >
> > from scipy import *
> > A = array(([1.,-1.],[-1.,1.]))
> > B = array(([1.0,0.0],[0.0,0.0]))
> > w,vr = linalg.eig(A,B)
> > print w
> > print vr
> >
> > Traceback (most recent call last):
> >   File "sevp.py", line 4, in ?
> >     w,vr = linalg.eig(A,B)
> >   File "/usr/lib/python2.3/site-packages/scipy/linalg/decomp.py", line
> > 109, in eig
> >     return _geneig(a1,b,left,right,overwrite_a,overwrite_b)
> >   File "/usr/lib/python2.3/site-packages/scipy/linalg/decomp.py", line
> > 68, in _geneig
> >     vr = _make_complex_eigvecs(w, vr, t)
> >   File "/usr/lib/python2.3/site-packages/scipy/linalg/decomp.py", line
> > 29, in _make_complex_eigvecs
> >     vnew.real = scipy_base.take(vin,ind[::2],1)
> > ValueError: matrices are not aligned for copy
> 
> This error is caused by the fact that in Numeric (1+0j)/0.0
> returns Inf+Nan*1j while Matlab returns Inf. It is a matter
> of taste which result to prefer. I (or Travis?) will see if scipy.eig
> can be made robust in this case. The code to be fixed is given in line
>    w = (alphar+_I*alphai)/beta
> of linalg/decomp.py where for your example
>   alphar = [0,1.41421356]
>   alphai = [0,0]
>   beta = [0.70710678,0]
> This results in
>   w = [0+0j,inf+nanj]
> which should be
>   w = [0+0j,inf+0j]
> in order to get Matlab's result.
> 
> > Any suggestion ?
> 
> Use small perturbation for now:
> 
> In [1]: from scipy import *
> 
> In [2]: w,vr=linalg.eig([[1,-1],[-1,1]],[[1,0],[0,1e-18j]])
> 
> In [3]: print w
> [  0.00000000e+000 +0.00000000e+000j               nan             +infj]
> 
> In [4]: print vr
> [[ 1.+0.j  0.+0.j]
>  [ 1.+0.j  1.+0.j]]
> 
> Pearu
> 
Pearu,

The following reference might be of interest in this context
http://www.cs.umu.se/research/nla/singular_pairs/guptri/

Nils

> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-dev