5

I'm trying to solve a large eigenvalue problem with Scipy where the matrix A is dense but I can compute its action on a vector without having to assemble A explicitly. So in order to avoid memory issues when the matrix A gets big I'd like to use the sparse solver scipy.sparse.linalg.eigs with a LinearOperator that implemements this action.

Applying eigs to an explicit numpy array A works fine. However, if I apply eigs to a LinearOperator instead then the iterative solver fails to converge. This is true even if the matvec method of the LinearOperator is simply matrix-vector multiplication with the given matrix A.

A minimal example illustrating the failure is attached below (I'm using shift-invert mode because I am interested in the smallest few eigenvalues). 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!

Edit: I should clarify what I mean by "make this work" (thanks, Dietrich). The example below uses a random matrix for illustration. However, in my application I know that the eigenvalues are almost purely imaginary (or almost purely real if I multiply the matrix by 1j). I'm interested in the 10-20 smallest-magnitude eigenvalues, but the algorithm doesn't behave well (i.e., never stops even for small-ish matrix sizes) if I specify which='SM'. Therefore I'm using shift-invert mode by passing the parameters sigma=0.0, which='LM'. I'm happy to try a different approach so long as it allows me to compute a bunch of smallest-magnitude eigenvalues. 

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
A_op = aslinearoperator(A)

# Try to solve the same eigenproblem again.
# This time it produces an error:
#
# ValueError: Error in inverting M: function gmres did not converge (info = 1000).
eigvals2 = eigs(A_op, sigma=0.0, which='LM', return_eigenvectors=False)
cilix
  • 1,262
  • 1
  • 14
  • 14
  • Did you ever find an answer? This is pretty much the only result in Google for that specific error. – Jared Feb 22 '22 at 17:03

2 Answers2

0

I tried running your code, but not passing the sigma parameter to eigs() and it ran without problems (read eigs() docs for its meaning). I didn't see the benefit of it in your example.

Dietrich
  • 5,241
  • 3
  • 24
  • 36
  • Dietrich, many thanks for the reply! I agree that from the way my question is phrased the need for the 'sigma' parameter isn't clear. Essentially, I need it to compute the first few smallest-magnitude eigenvalues using shift-invert mode. Unfortunately, this means I can't just leave out the 'sigma' parameter because it would then compute the largest-magnitude eigenvalues. Also, simply passing "which='SM'" doesn't work because the algorithm never finishes. I've edited the question to include this information. Does it make more sense now? – cilix Jun 15 '14 at 12:48
  • I am not an expert on eigenvalues, but I'm not sure if you're approach is feasible: All the algorithms I know of do not return the eigenvalues in a specific order, so when you have one, you don't know if it is one of the smaller ones. Furthermore, what I know of gmres is that it doesn't have a guarantee to converge and small eigenvalues, especially in large matrices, are quite susceptible to numeric errors. Perhaps you have some other structure in your matrix (like arranging it as a banded matrix), which can help you out. – Dietrich Jun 15 '14 at 15:15
0

Eigs can already find the smallest eigenvalues first. Set which = 'SM'

Phooey
  • 1
  • Thanks for the suggestion. I did try which='SM' (see the edit to my question) but unfortunately for my application it didn't behave well and the computations never finished. – cilix Apr 08 '16 at 13:09