8

I am trying to find the first three lowest eigenvalues of large sparse matrices of size range $10^3 - 10^5$. The matrices depend on some parameter $x$, so I first construct the matrices and then use Eigenvalues[Mat(x), 3]. Mathematica however orders the eigenvalues by absolute value so there is no guarantee in finding the lowest eigenvalue by the method above. See figure below

enter image description here

I can find all the eigenvalues, order them, and then find the minimum using the instructions found here but that defeats the point of using sparse matrices or Mathematica altogether. Moreover, the code slows down even more when cycling through all values of $x$. So, can this problem be done in Mathematica?

WikawTirso
  • 83
  • 4

2 Answers2

12

Use the Arnoldi method with shift-inversion:

Eigenvalues[A, 3, Method -> {"Arnoldi", "Criteria" -> "Magnitude", "Shift" -> 0}]

gives you the three smallest eigenvalues by absolute value (by magnitude). See here: Efficiently find all values of parameter such that any of the eigenvalues of a matrix is equal to 1

After comments by @HenrikSchumacher it appears that the same can be achieved with

Eigenvalues[A, -3, Method -> {"Arnoldi", "Criteria" -> "Magnitude"}]

and there is no need for explicit shift-inversion.

And as @CarlWoll points out, this method is the default method for sparse matrices, so even

Eigenvalues[A, -3]

achieves the same effect. For non-sparse matrices, however, specifying the method can give a large speedup.

update

It now looks like you want the eigenvalues with smallest real part, not those with smallest magnitude. This you can also achieve with the Arnoldi method:

Eigenvalues[M, 3, Method -> {"Arnoldi", "Criteria" -> "RealPart"}]

gives the three eigenvalues with largest real part. To get the smallest ones, do

-Eigenvalues[-M, 3, Method -> {"Arnoldi", "Criteria" -> "RealPart"}]

The trick is that M and -M have the same eigenvalues (up to sign) and eigenvectors.

Of course this trick also works with the other methods presented above.

Roman
  • 47,322
  • 2
  • 55
  • 121
  • Reading this, I realized once more that Mathematica's implementation of Arnoldi's method is really inconsistent: Both Eigenvalues[A, 1, Method -> {"Arnoldi", "Criteria" -> "Magnitude", "Shift" -> 0}] and Eigenvalues[A, -1, Method -> {"Arnoldi", "Criteria" -> "Magnitude", "Shift" -> 0}] return the same. – Henrik Schumacher Jun 10 '19 at 21:30
  • @HenrikSchumacher I would report this to Wolfram support. This problem only happens with the Shift option. – Roman Jun 10 '19 at 21:43
  • I have reported it already ([CASE:4270162]), but it cannot harm to build up some peer pressure. =) – Henrik Schumacher Jun 10 '19 at 21:49
  • 1
    I believe Eigenvalues with the default method already uses the "Arnoldi" method when the input matrix is sparse and only a limited set of eigenvalues is requested. So, for sparse matrices, Eigenvalues[A, -1], Eigenvalues[A, -1, Method->"Arnoldi"] and Eigenvalues[A, -1, Method->{"Arnoldi", "Criteria"->"Magnitude"}] should all be equivalent. – Carl Woll Jun 10 '19 at 22:08
  • @CarlWoll I think you're right. – Roman Jun 10 '19 at 22:21
  • Can we use Eigenvalues[A, -3, Method -> {"Arnoldi", "Criteria" -> "RealPart"}] to get the eigenvalues with the smallest real part instead of using the "trick" with the minus sign? – noir1993 Mar 26 '23 at 01:54
  • @noir1993 No, this does not seem to work, unfortunately. Given the name of the criterion, I think it should work; but it doesn't. What you suggest gives a very strange answer, which is neither the smallest by real part nor the smallest positive by real part. – Roman Mar 26 '23 at 07:07
  • @Roman Can you give the matrix you are using to test these? I assume it's a small matrix to quickly check the output. – noir1993 Mar 26 '23 at 08:30
  • 1
    For a minimal example, try A = SparseArray[{{1, 1} -> -1., {3, 3} -> 2.}]. Eigenvalues[A, 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}] gives {2.} and -Eigenvalues[-A, 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}] gives {-1.}; but Eigenvalues[A, -1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}] gives the nonsensical {2.00108}. – Roman Mar 26 '23 at 15:26
  • @Roman This is extremely intriguing! I think the problem occurs if the matrix has a zero or negative eigenvalues. It seems that if all the eigenvalues are positive and greater than zero, then there is likely not going to be a problem. Would you agree with this assessment? Furthermore, for your example, Eigenvalues[Mat, -1, Method -> {"Arnoldi", "Criteria" -> "RealPart", "BasisSize" -> 3, "Shift" -> 10}] cures the problem, the shift has to be larger than the largest eigenvalue. – noir1993 Mar 27 '23 at 13:17
  • @noir1993 if you want to learn more about this topic, I suggest you look for references directly on ARPACK and its known problems instead of looking on this Mathematica forum. – Roman Mar 27 '23 at 13:53
4

To calculate the lowest eigenvalues using Mathematica, I always introduce a "shift" in the following way:

mat1 = mat - IdentityMatrix[Length[mat]]*large

and then add large to the result of Eigenvalues[mat1]. This operation leaves the eigenvectors unchanged and I do not need to use Arnoldi specifically. Of course, Arnoldi does a similar shift inside, however I am not sure if that is used correctly in the Mathematica implementation.

Example: 

mat={{1.,2.},{3.,4.}};

The eigenvalues are in the order given by Eigenvalues

{5.37228, -0.372281}.

The procedure above produces for any large enough large:

{-0.372281, 5.37228}

with the smallest first. Or Eigenvalues[mat1,1]+large just the lowest eigenvalue.

Michael Weyrauch
  • 821
  • 5
  • 10
  • It worked! I also realize original question was needlessly ambiguous. I get the general idea of the shift but is there a way to systematically guess a value for large? Or alternatively, 'large' compared to what (|large| > 1)? Thanks. – WikawTirso Jun 12 '19 at 13:16
  • Well, the idea is to introduce a shift that all eigenvalues of the "shifted" matrix are negative. So you could calculate the largest eigenvalue of the shifted matrix, and if is negative, your shift is fine. if not, shift more. – Michael Weyrauch Jun 12 '19 at 14:16