Given an $n\times n$ matrix $Q$ (with e.g. $n\approx10^4$) I am only interested in the 3rd smallest eigenvalue of $Q,$ and not the entire spectrum (assume all eigenvalues are real, e.g. a Hermitian matrix). So one direct way would be to solve the Eigensystem, order the eigenvalues and pick the 3rd one:

{eigv, U} = Eigensystem[Q];
order = Ordering[eigv];
eigv = eigv[[order]];
res = eigv[[3]];
  1. (Main question): But is there a way to go about this in a more targeted and efficient way considering that we know which eigenvalue we are looking for?

  2. In the above I order the eigenvalues, but if we don't, does Eigensystem arrange them in a specific order already?


The matrix can has only nonnegative eigenvalues (M-matrix), and its first eigenvalue, $\lambda_0,$ is typically $0.$ Therefore, the idea is to find the 3rd smallest eigenvalue $\lambda$ in magnitude, knowing that the real parts of all $\lambda$ are nonnegative. These properties are e.g. satisfied for Laplacian matrices of graphs, in this case graphs with $10^4$ nodes.

  • 1
    $\begingroup$ Are your eigenvalues guaranteed to be real? If not, what do you mean by smallest? $\endgroup$ – mikado Oct 26 '19 at 16:49
  • $\begingroup$ @mikado good question, yes let us only consider the case where all eigenvalues are indeed real. $\endgroup$ – user52181 Oct 26 '19 at 16:51
  • 3
    $\begingroup$ Eigensystem[Q, 3, Method -> {"Arnoldi", "Criteria" -> "RealPart"}]could be one choice. $\endgroup$ – yarchik Oct 26 '19 at 16:58
  • $\begingroup$ Is you matrix numeric and much larger then 3 by 3, or is it a small symbolic matrix? $\endgroup$ – yarchik Oct 26 '19 at 17:11
  • 1
    $\begingroup$ In that case I'd suggest finding the largest eigenvalues (call it eigbig) using Eigenvalues[mat,1], then find the three largest in magnitude of mat-eigbig*identity and add eigbig to them. Can use Eigensystem for this second step to get the corresponding eigenvectors. $\endgroup$ – Daniel Lichtblau Oct 28 '19 at 16:00

Here is an example of what I had suggested in a comment. We will work with a 100x100 symmetric sparse matrix, and make it diagonally dominant (or close enough to that) so that all eigenvalues are positive.

rng = 100;
nonzerocount = 3*rng;

mat0 = SparseArray[
   Thread[RandomInteger[{1, rng}, {nonzerocount, 2}] -> 
     RandomReal[{0, 1}, nonzerocount]], {rng, rng}];
mat = mat0 + Transpose[mat0] + SparseArray[{a_, a_} -> 5, {rng, rng}];

First check the smallest eigenvalues using the method we want to avoid:

Eigenvalues[mat, -5]

(* Out[337]= {2.6747750632, 2.43141432273, 2.40432837499, 2.30813913679, \
1.93291151486} *)

So these are values we hope to attain without requiring an linear equation solving under the hood.

First get the biggest eigenvalue.

eigbig = First[Eigenvalues[mat, 1]]

(* Out[338]= 8.94206162291 *)

Now shift by the negative of this largest, and obtain the three most negative eigenvalues, along with their corresponding eigenvectors. Shift back to get the correct eigenvalues, that is, the smallest positive ones from the original matrix.

matShifted = mat - SparseArray[{a_, a_} -> eigbig, {rng, rng}];
smallEig = Eigensystem[matShifted, 3];
smallEig[[1]] + eigbig

(* Out[344]= {1.93291151486, 2.30813913679, 2.40432837499} *)

We have recovered the ones we wanted.

When I instead use a dimension of 10000 it takes around 0.1 seconds on my desktop to find the largest eigenvalue, and around 0.5 seconds to get those eigenvalues and corresponding eigenvectors. Finding them using Eigensystem[mat,-3], in contrast, requires 6.8 seconds.

| improve this answer | |
  • $\begingroup$ If I may ask: in using this method for my matrices (which are KirchhoffMatrix[] of graphs, I keep getting the warning: Eigenvalues::maxit2: Warning: maximum number of iterations, 5000, has been reached by the Arnoldi algorithm without convergence to the specified tolerance, but the current best computed value has been returned. You can use method options with Method -> {Arnoldi, opts} to increase the size of basis vectors, the maximum number of iterations, .... Should I be worried? from your experience, how should I workaround the convergence problem? Many thanks for any hints! $\endgroup$ – user52181 Nov 13 '19 at 9:57
  • $\begingroup$ I assume this badness happens when you do Eigenvalues[mat, 1]. It is not a familiar problem for me. My advice would be to increase the iterations, and maybe also extract the two or three largest eigenvalues to see if they are close to one another (multiplicity or near-multiplicity might account for a convergence issue, I think). Since the eigenvalues are apparently all positive, there might also be tricks to play using SingularValueList, again extracting only the top few. $\endgroup$ – Daniel Lichtblau Nov 13 '19 at 15:00
  • $\begingroup$ Good point about closeness of values, maybe it is the case, I'm not sure to be honest. Here are 2 examples of 3 largest eigens: {3.1433639275971204, 3.0693111881309645, 3.0135000156020246} and {3.110447056413233, 3.080469746360026, 2.9917100715725695}. $\endgroup$ – user52181 Nov 14 '19 at 12:36
  • 1
    $\begingroup$ Those values do not strike me as being terribly close. Another possibility occurred to me. You could try working in extended precision. This would slow things substantially but might do a better job of handling the matrix. $\endgroup$ – Daniel Lichtblau Nov 14 '19 at 17:08
  • $\begingroup$ Indeed! Does that mean simply wrapping Eigenvalues[...] in N[] to enforce machine precision? $\endgroup$ – user52181 Nov 14 '19 at 17:40

As yarchik says in the comments, you can use the "Arnoldi" method. However, the "Arnoldi" method doesn't work well for finding the smallest real eigenvalues (perhaps a bug), but you can instead find the largest real eigenvalues after changing the sign of the matrix. Here is an example Hermitian matrix (so eigenvalues are real):

m = (ConjugateTranspose[#] + #)& @ RandomComplex[10+10I, {10, 10}];

Finding the full eigenvalues (I use Eigenvalues, but Eigensystem works similarly):

Eigenvalues[m] //ReverseSort

{98.0228, 23.1452, 18.0616, 12.4434, 8.72238, -0.890917, -9.86992, -13.0016, \ -20.4409, -26.7991}

Now, using the "Arnoldi" method:

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


| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy