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.

  • 1
    $\begingroup$ 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"? $\endgroup$ – b3m2a1 Jan 16 at 19:53
  • $\begingroup$ 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. $\endgroup$ – Sjoerd Smit Jan 16 at 20:05
  • 1
    $\begingroup$ 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? $\endgroup$ – Henrik Schumacher Jan 16 at 21:32
  • 1
    $\begingroup$ @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. $\endgroup$ – 56th Jan 17 at 1:12
  • 1
    $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Jan 17 at 7:23

Your Answer

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

Browse other questions tagged or ask your own question.