Speeding Up Sparse Matrix Generalized Eigenvalues

Question

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];

Answer 1

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.1`10 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.