[SciPy-User] Error when computing eigenvalues of a LinearOperator: "gmres did not converge"

[Maximilian Albert]

Hi,

This is my first post to this mailing list, so let me take the opportunity
to thank everyone involved for their great work! I use scipy on a daily
basis and couldn't imagine what I would do without it.

I'm trying to solve large eigenvalue problems using scipy's sparse linear
algebra module. It works fine if I apply the solver
scipy.sparse.linalg.eigs to a regular numpy array A. However, my matrices
get quite big and I'm having problems with memory consumption. So ideally I
would like to use a LinearOperator instead. (The matrix A is dense, but I
can compute the action of A on a vector without having to assemble A
explicitly, so I'm hoping this is still a valid use case for
LinearOperator.)

The problem is that if I try this then the iterative solver doesn't
converge. This is true even if the 'matvec' method of the LinearOperator is
simply matrix-vector multiplication with the (dense) matrix A.

A minimal example illustrating the failure is attached below. This computes
the eigenvalues of a random matrix A just fine, but fails when applied to a
LinearOperator that is directly converted from A. I tried to fiddle with
the parameters for the iterative solver (v0, ncv, maxiter) but to no avail.

Am I missing something obvious? Is there a way to make this work? Any
suggestions would be highly appreciated.

Many thanks and best regards,
Max

==>
from scipy.sparse.linalg import eigs, LinearOperator, aslinearoperator
import numpy as np

# Set a seed for reproducibility
np.random.seed(0)

# Size of the matrix
N = 100

# Generate a random matrix of size N x N
# and compute its eigenvalues
A = np.random.random_sample((N, N))
eigvals = eigs(A, sigma=0.0, which='LM', return_eigenvectors=False)
print eigvals

# Convert the matrix to a LinearOperator and
# try to solve the same eigenproblem again.
# The call to 'eigs' will produce an error:
#
# ValueError: Error in inverting M: function gmres did not converge (info =
1000).
A_op = aslinearoperator(A)
eigvals2 = eigs(A_op, sigma=0.0, which='LM', return_eigenvectors=False)
<==
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20140612/efdb4342/attachment.html>