# Matrix-free Arnoldi method for eigensystems

I am solving a generalized eigenvalue problem

$$\mathbf S\,\mathbf x = \lambda\,\mathbf M\,\mathbf x$$

w/ $$\mathbf S := \mathbf B\,\mathbf A^{-1}\,\mathbf B^{T}$$, and $$\mathbf A$$ is a sparse matrix (SparseArray[]):

S = B.Inverse@A.Transpose@B;
Eigenvalues[{S, M}, 1, Method -> "Arnoldi"]


The problem is, computing Inverse@A is not feasible already for the matrix size 32 736 ($$\mathbf A^{-1}$$ is dense and takes ~8.5 gigabytes). I do not really need $$\mathbf A^{-1}$$ explicitly for the iterative Arnoldi method, and in principle I can approximate its action e.g. via a preconditioned Krylov solver.

Is there a way to do this? I.e. instead of passing $$\mathbf S$$ to Eigenvalues[] pass something that mimics multiplication of $$\mathbf S$$ by a vector.

• I think the title on this could be a bit more descriptive of what your specific issue is. Maybe something like "Solving generalized Eigenvalue problems with SparseArray"? – b3m2a1 Jan 16 at 19:53
• Well, there is LinearSolve[A] which is the practical alternative to Inverse[A] for large matrices. I don't know if you can throw that into Eigenvalues`, but it's a start. – Sjoerd Smit Jan 16 at 20:05
• Unfortunately, I don't know an Anoldi imeplementation in Mathematica that allows a function that computes a matrix-vector as instead of a matrix itself as argument. It would be great to have one (therefore +1 for bringing this up). If $\mathbf{A}$ is a dense matrix of size $32 736 \times 32 736$, that's really bad anyway. Do you have a very good preconditioner for it? In general, we need much more context. Where does $A$ come from? Maybe there is some structure that can be exploited. What are $\mathbf{B}$ and $\mathbf{M}$? Maybe their inverses are easier to compute? – Henrik Schumacher Jan 16 at 21:32
• @henrik-schumacher $\mathbf A$ itself is sparse; $\mathbf A^{-1}$ is not. It comes from Trace/Cut FEM discretization of the surface Stokes problem, so $\mathbf A$ is a surface vector Laplace block, $-\mathbf B$/$\mathbf B^{T}$ are divergence/gradient blocks, and $\mathbf M$ is a pressure mass matrix. All of them are sparse and a good preconditioner for $\mathbf A$ is available. – 56th Jan 17 at 1:12
• Hm. I see. So you really need only a matrix-free Arnoldi routine. Unfortunately, Mathematica is shipped without one (to my knowledge). =( You would have to use Library Link to a library that is capable of doing this, for example the "Extended Eigensolver Interface" provided by Intel MKL. – Henrik Schumacher Jan 17 at 7:23