# Speeding Up Sparse Matrix Generalized Eigenvalues

I have a generalized eigenvalue problem $${\mathbf A}\vec v = \lambda {\mathbf B} \vec v$$ for which I am trying to find the smallest eigenvalue $$\lambda$$ and the associated eigenvector $$\vec v$$. Note that $${\mathbf A}$$ and $${\mathbf B}$$ are both sparse (they are $$N\times N$$ but contain $${\mathcal O}(N)$$ non-zero entries). I then use the following function to find the smallest magnitude eigenvalue and associated eigenvectors $$\{{\lambda,v}\} = {\texttt{Eigensystem[{A, B}, -1]}}\,.$$ This works great when $${\mathbf A}$$ and $${\mathbf B}$$, which are both defined as $$\texttt{SparseArray}$$ objects, are smaller than about $$600\times 600$$, but becomes extremely slow when the matrices become marginally bigger. Clearly I am doing something wrong, because the same code runs in matlab in a fraction of a second for much bigger matrices.

How can I speed up this computation? Are there options I have neglected to set correctly? Is there another function I should be using?

Here is an example problem that exhibits the same behavior

Clear["Global*"]
GridSize = 500;
dx = 0.1;
ones = Table[1, {ii, 1, GridSize}];
DL = 1/dx^2 ones;
DL = Drop[DL, -1];
DM = -(2/dx^2) ones;
DR = DL;
RHS = ones;
RHS[[-1]] = 0;
A = SparseArray[{Band[{1, 2}] -> DL, Band[{1, 1}] -> DM,
Band[{2, 1}] -> DR}, {GridSize, GridSize}];
B = SparseArray[{Band[{1, 1}] -> RHS}, {GridSize, GridSize}];

EV = Eigensystem[{A, B}, -1];

• How are A and B defined? Are they actually SparseArray objects? Can you give an example? Sep 4, 2022 at 22:18
• I am trying to think of a good example that isn't the actual problem I'm working with, but yes they are defined as sparse-array objects.
– Guy
Sep 4, 2022 at 22:20
• @LukasLang I've included an example. What makes this problem slow is the non-invertability of B.
– Guy
Sep 4, 2022 at 22:38
• “Clearly I am doing something wrong” Not necessarily. Unlike the usual eigenvalue problem, there doesn't exist an efficient algorithm for generalized eigenvalue problem AFAIK. I sincerely hope I'm wrong. Sep 5, 2022 at 2:36
• @xzczd "For linear algebra, MATLAB and Mathematica uses the same library under the hood AFAIK." This is very true in general. They both use the Intel MKL on x86 processors. But I think to remember that some guys from WR once told me on this site that the eigensolvers where home baked. And btw., working now with an Apple M1, I can say that Mathematica does not use MKL there (because it runs "natively"). I know that very well, because I am deeply missing Intel Pardiso. (sniff) But Matlab still requires MKL (and thus is run through the x86 compatibility mode Rosetta 2). Jan 5 at 23:04

Since you seek a numerical result and your matrices are machine-precision, you can use the Arnoldi method, which is optimized to find only some eigenvalues. You also ask specifically for eigenvalues, so you could use Eigenvalues rather than Eigensystem:

GridSize = 100000;

dx = 0.1;
ones = Table[1, {ii, 1, GridSize}];
DL = 1/dx^2 ones;
DL = Drop[DL, -1];
DM = -(2/dx^2) ones;
DR = DL;
RHS = ones;
RHS[[-1]] = 0;
A = SparseArray[{Band[{1, 2}] -> DL, Band[{1, 1}] -> DM, Band[{2, 1}] -> DR}, {GridSize, GridSize}];
B = SparseArray[{Band[{1, 1}] -> RHS}, {GridSize, GridSize}];

EV = Eigenvalues[{A, B}, -1, Method -> "Arnoldi"]

(* Out: {-9.86941*10^-8} *)


Timing indicates that the calculation does indeed take a fraction of a second even for a grid size of $$100 000\times100000$$ in MMA as well:

RepeatedTiming[Eigenvalues[{A, B}, -1, Method -> "Arnoldi"];]

(* Out: {0.621861, Null} *)


These calculations are done at machine precision, which would be the default in Matlab. If you think you need more precision, you can change the definition of dx to include an arbitrary precision specification, e.g. dx = 0.110 to carry out calculations with 10 digits of precision. This will inevitably be much slower and, in this case, it did not change the result when I tried it. It would be helpful if you could add to your question what you expect the result to be.

• Thanks for this answer! I suppose my simplified example failed to capture some of the key features of my problem. In particular, my matrix is non-hermitian, so arnoldi doesn't work. I'll work out a better non-hermitian example problem and put it in the main question.
– Guy
Sep 5, 2022 at 20:24
• Perhaps there is a way to embed a non-hermitian operator in a larger space where it is hermitian? Sounds unlikely but maybe that's what matlab is doing?
– Guy
Sep 5, 2022 at 20:33
• @Guy I don't think the Arnoldi method requires its input to be Hermitian. Consider: a = {{1, 2, 3}, {2, 3, 3}, {2, 1, 0}}; HermitianMatrixQ[a] returns False, and yet Eigenvalues[N@a, 1, Method -> "Arnoldi"] works fine. The input must be either machine-precision, or of finite arbitrary precision though (which you could achieve with N or SetPrecision). Sep 5, 2022 at 21:57
• Interesting - it doesn't seem to work for a generalized eigenvalue problem though. For instance a = {{1, 2, 3}, {2, 3, 3}, {2, 1, 0}}; b = {{1, 1, 3}, {2, 5, 3}, {4, 1, 0}}; Eigenvalues[{N@a, N@b}, 1, Method -> "Arnoldi"] leads to the output Eigenvalues::herm: The matrix {{1.,1.,3.},{2.,5.,3.},{4.,1.,0.}} is not Hermitian or real and symmetric.
– Guy
Sep 5, 2022 at 23:38