$\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.

  • $\begingroup$ Are you interested in the smallest eigenvalues by magnitude (i.e., those closest to zero, even if there are negative ones) or the smallest eigenvalues algebraically? This makes a huge difference. $\endgroup$ Jul 3 '16 at 8:55
  • $\begingroup$ By magnitude, those have physical meaning $\endgroup$
    – brubeck
    Jul 3 '16 at 13:53

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.

  • 2
    $\begingroup$ Thanks! I also found out that Matlab's eigs implements ARPACK and that it can accept the linear operator as a function. For the case of the smallest eigenvalues, you must provide the way to compute the inverse operator instead. For my particular case this means to solve Sylvester's equation. $\endgroup$
    – brubeck
    Jul 3 '16 at 4:46
  • $\begingroup$ I believe if you implement ARPACK directly, say through C or Fortran, you can compute this directly without providing the inverse operation. However I could understand if you don't want to deal with that! $\endgroup$ Jul 3 '16 at 5:45
  • $\begingroup$ As far as I know, iterative solvers converge to the largest-magnitude eigenvalues. Since the smallest ones of A are the largest ones of inv(A), the algorithm requires the inverse operator. How could this be possible the other way? $\endgroup$
    – brubeck
    Jul 3 '16 at 13:59
  • $\begingroup$ Still, I managed to optimize the computation of inverse operator in O(n^3), so this is not a huge bottleneck anymore. $\endgroup$
    – brubeck
    Jul 3 '16 at 14:02
  • $\begingroup$ Yes this is my mistake, you do need the inverse operator for those smallest values (in my field we never compute the smallest ones, whoops). $\endgroup$ Jul 3 '16 at 17:47

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.)


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.


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.