Eigenvectors of Black-box matrix

Question

$\DeclareMathOperator{\diag}{diag}$ Consider the generalized eigenproblem $A\mathbf{x}=\lambda B\mathbf{x}$. When solving PDEs numerically (specifically, I am interested on finding the Dirichlet eigenmodes of the 2D laplacian in orthogonal coordinates, i.e. $\nabla^2 u=\lambda u$) this problem takes the form

$$ A_1X + XA_2^\text{T} = \lambda B \circ X$$

where $X\in\mathbb{R}^{m\times n}$ is our eigenvector (the eigenfunction $u$ sampled at the meshpoints), $A_1\in\mathbb{R}^{m\times m}$ and $A_2\in\mathbb{R}^{n\times n}$ are differentiation matrices, $B\in\mathbb{R}^{m\times n}$ is the Jacobian determinant (sampled at the meshpoints) and the symbol $\circ$ denotes element-wise multiplication. I am only interested on the eigenmodes corresponding to the smallest eigenvalues.

With the help of the Kronecker product, the problem can be vectorized as follows:

$$ (I \otimes A_1 + A_2 \otimes I)\mathbf{x} = \lambda \diag(B)\mathbf{x}$$

where $\mathbf{x}\in\mathbb{R}^{mn}$ is a vector formed by stacking the columns of $X$ and $\diag{(B)}$ is a diagonal matrix formed with the columns of $B$.

The problem with this is that the resulting matrices are unnecessarily sparse and have size of $mn\times mn$, and when I try to solve it for modestly large sizes of $m$ and $n$, say 200, my computer runs out of memory. Clearly that is not the way.

I understand that modern eigenvalue algorithms are highly sophisticated, and that coding one by myself is not the best option. So I am looking for an iterative solver that only depends on the computation of the matrix-vector products with $A$ and $B$ and that allows the user to provide them as a black-box.

Answer 1

I recommend ARPACK which provides matrix-free routines for generalized eigenvalue problems with examples in their documentation for this purpose. This method is one of the most widely used, especially for relatively large problems where you are only searching for a few eigenvalues/vectors (such as the smallest, in your case).

Further, eigenvalue problems of the size $10^4$ are not especially large and are likely solvable on any modern workstation. However, if this problem gets bigger the gold standard is SLEPc, which is a library containing many eigenvalue solvers and can distribute your problem across many compute nodes.

Answer 2

You are right that writing an efficient and robust eigenvalue solvers for large sparse systems is a difficult task that should only be done yourself as the very last resort. Spencer already gives the major players in his answer (ARPACK for single-node multithreaded computing, which is part of MATLAB and SciPy, and SLEPc for distributed computing, which can be accessed in Python via slepc4py). Both solve generalized eigenvalue problems, and both can take a procedure that computes matrix-vector products instead of a matrix.

Since you are interested in the smallest magnitude eigenvalues, you probably need to make use of the shift-and-invert strategy (i.e., compute the smallest magnitude eigenvalues of $A^{-1}$. Of course, you never actually compute inverse matrices; instead, whenever you need to compute $A^{-1}v$ for some vector $v$, you solve $Aw=v$ and use $w$ for $A^{-1}v$. There are several strategies to accelerate this, based on the fact that ARPACK/SLEPc also only require a procedure that provides the application of the inverse (the obvious ones should already be implemented in ARPACK/SLEPc):

1. Since you need to solve linear systems repeatedly, you compute a LU- or Cholesky factorization of the matrix $A$ ahead of time and use the (much faster) backward/forward substitution part of the linear solve during the iteration.

2. In your case, you probably don't want to build the big tensor product matrix in the first place; in this case it makes sense to use an iterative solver such as CG for the solution of $Aw=v$ as well (especially if you can provide a good preconditioner). In this way, the whole eigenvalue solver only requires a procedure to perform the matrix-vector multiplication. (You probably have to play with the tolerances a bit to get good results.)

Answer 3

I completely agree with the answers posted so far. I would like to point out that the matrix

$$ A = (I \otimes A_1 + A_2 \otimes I)$$

never needs to be formed explicitly. As other answers correctly point out that one only needs to form matrix-vector products (matvecs) with A, or its inverse. Matvecs of the form $Ax$ can be done using the identity

$$ (B^T \otimes C)x = \text{vec}(CXB), $$

where $\text{vec}(X) = x$. So, in your case

$$ Ax = \text{vec}(A_1X) + \text{vec}(XA_2^T).$$

For applying the inverse operator to a vector (as needed for the smallest generalized eigenvalues), you could use a Krylov subspace method. You could use the same trick for the matvecs.