[SciPy-User] scipy.sparse.linalg.eigen error code?

[Pauli Virtanen]

Wed, 13 Oct 2010 15:55:42 +0200, Nico Schlömer wrote:
> I'm using scipy.sparse.linalg.eigen to compute the lowest magnitude
> eigenvalue of a matrix, and I noticed that sometimes the code would
> return 0+0j where I really didn't expect it. Turns out that tweaking the
> number of iterations results in some more meaningful value here,
> suggesting that the underlying Arnoldi (ARPACK?) iteration failed with
> the maxiter value previously given. As far as I can see, there's no way
> to tell that the iteration actually failed.
> 
> Is that correct?

Yep. It seems to try to raise a warning, though,

    elif self.info == -1:
        warnings.warn("Maximum number of iterations taken: %s" % 
self.iparam[2])

But this is bad practice. You'll only see the warning *once*, and there 
is indeed no way to catch it.

The Arpack interface seems to need some fixing, too. Exceeding the number 
of iterations is really an error condition. The interface also apparently
doesn't check that ncv is sensible, so in some cases (small matrices) it
can give "Error -3".

As a work-around, you can probably use

def eigen_fixed(A, k=6, M=None, sigma=None, which='LM', v0=None,
          ncv=None, maxiter=None, tol=0,
          return_eigenvectors=True):
    import numpy as np
    from scipy.sparse.linalg.interface import aslinearoperator
    from scipy.sparse.linalg.eigen.arpack.arpack import _UnsymmetricArpackParams

    A = aslinearoperator(A)
    if A.shape[0] != A.shape[1]:
        raise ValueError('expected square matrix (shape=%s)' % (A.shape,))
    n = A.shape[0]

    matvec = lambda x : A.matvec(x)
    params = _UnsymmetricArpackParams(n, k, A.dtype.char, matvec, sigma,
                           ncv, v0, maxiter, which, tol)

    if M is not None:
        raise NotImplementedError("generalized eigenproblem not supported yet")

    while not params.converged:
        params.iterate()

    if params.info == -1:
        raise RuntimeError("Did not converge")

    return params.extract(return_eigenvectors)

-- 
Pauli Virtanen