title: Error when computing eigenvalues of a scipy LinearOperator: "gmres did not converge"
body: <p>I'm trying to solve a large eigenvalue problem with Scipy where the matrix <code>A</code> is dense but I can compute its action on a vector without having to assemble <code>A</code> explicitly. So in order to avoid memory issues when the matrix A gets big I'd like to use the sparse solver <code>scipy.sparse.linalg.eigs</code> with a <code>LinearOperator</code> that implemements this action.</p>
  
  <p>Applying <code>eigs</code> to an explicit numpy array <code>A</code> works fine. However, if I apply <code>eigs</code> to a <code>LinearOperator</code> instead then the iterative solver fails to converge. This is true even if the <code>matvec</code> method of the <code>LinearOperator</code> is simply matrix-vector multiplication with the given matrix <code>A</code>.</p>
  
  <p>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 <code>A</code> just fine, but fails when applied to a <code>LinearOperator</code> that is directly converted from <code>A</code>. I tried to fiddle with the parameters for the iterative solver (<code>v0</code>, <code>ncv</code>, <code>maxiter</code>) but to no avail.</p>
  
  <p>Am I missing something obvious? Is there a way to make this work? Any suggestions would be highly appreciated. Many thanks!</p>
  
  <p><strong>Edit:</strong> 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 <code>1j</code>). 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 <code>which='SM'</code>. Therefore I'm using shift-invert mode by passing the parameters <code>sigma=0.0, which='LM'</code>. I'm happy to try a different approach so long as it allows me to compute a bunch of smallest-magnitude eigenvalues. </p>
  
  <pre><code>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)
  </code></pre>