[SciPy-dev] Scipy cannot handle singular eigenvalue problems

[Pearu Peterson]

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