How to find the smallest positive eigenvalue of a large general system if they are all in +/- pairs of real eigenvalues

Question

I used MKL Lapack SBGVX because it can solve for "selected" eigenvalues/modes (both positive and negative), thinking it would be efficient, but it is extremely slow when compared with Bathe's subspace iteration eigensolver, which unfortunately can solve only for positive eigenvalues. I need only the smallest positive eigenvalue and its mode.

I am trying to solve the problem $$ (\mathbf{K} - e\,\mathbf{M})\,\mathbf{u} = \mathbf{0} $$ with $\mathbf{M}$ not +definite, and $\mathbf{K}$-- A-- that is +definite.

Here, $\mathbf{K}$ is has the form of a stiffness matrix and is invertible while $\mathbf{M}$ is like a mass matrix, some of the components of which can be negative (for instance, for negative-mass metamaterials). It is possible that some of the diagonal terms are zero, but let us assume that the matrix is invertible.

Since SBGVX needs the second matrix to be +definite, I multiply by $-1/e$ and rearrange it as $$ (\mathbf{M} - e'\,\mathbf{K})\,\mathbf{u} = \mathbf{0} $$ then look for the highest $e'$, which is positive, to then compute the lowest positive eigenvalue $e=1/e'$ . It works, but it is very slow.

Does anyone know how to modify Bathe's subspace iteration eigensolver to compute the highest eigenvalue instead of the smaller?

Or, how can I transform $$ (\mathbf{K} - e\,\mathbf{M})\,\mathbf{u} = \mathbf{0} $$ to solve for eigenvalues $e^2$ instead of $e$, so I can use Bathe's solver to find the smallest $e$? Bathe's solver crashes trying to solve for the smallest $e$ when the problem has negative eigenvalues.

Or, does anyone know another solver that I can try? A typical system has 5000 equations.

Answer 1

I am not an expert in this area, so there may be better methods for handling this problem, but what I suggest should work in principle; I don't know if it will converge faster than the methods you are currently using.

If $M$ is invertible, then you can transform your generalized eigenvalue problem

\begin{align} Ku = \lambda Mu \end{align}

into an ordinary eigenvalue problem

\begin{align} M^{-1}Ku = \lambda u \end{align}

and then solve it. Furthermore, if $K$ is positive definite, then $M^{-1}K$ should be invertible, and all (generalized) eigenvalues should be nonzero. If all (generalized) eigenvalues occur in pairs with opposite signs, then your problem reduces to finding the (generalized) eigenvalue with smallest magnitude, which you can do with something like an inverse power method. Once you find the smallest magnitude (generalized) eigenvalue, it should be straightforward to calculate the corresponding (generalized) eigenvector.

For a problem with 5000 degrees of freedom, it should be possible to invert $M$ with direct methods, calculate $M^{-1}K$ explicitly, and then use standard methods for calculating eigenvalues and eigenvectors. If $M$ is not invertible, you will need to use methods specialized for generalized eigenvalue problems. One software package that can handle these problems is SLEPc; it interfaces to a number of iterative eigenvalue solvers (ARPACK, BLOPEX, TRLAN, etc.).

Answer 2

$\def\m#1{\mathbf{#1}}$ I know my answer below is a bit off. Somehow, I misread the question that the eigenvalues would come in pairs $\pm\lambda$. On second reading you did not say so.

But, I do not delete the answer yet since it shows the principle how to do a double-shift. You can apply the same principle to calculate $\lambda^2$ instead of $\lambda$. Just solve the system twice in each step. To get $\lambda^2$ is even simplier than what I wrote below. You can choose the identity matrix $\m1$ for $f$ and shift two times with $\lambda$, i.e., $\bar K = K-\lambda\m1$.

---

If you apply an iterative method with zero shift to a matrix with eigenvalues of the same absolute value you can run into problems. A good example is the matrix \begin{align} K:= \begin{pmatrix} 0&1\\ 1&0 \end{pmatrix} \end{align} with eigenvalues $\lambda_{12}=\pm 1$ and eigen directions $e^{(1)}:=\frac1{\sqrt{2}}\begin{pmatrix}1\\-1\end{pmatrix}$, $e^{(2)}:=\frac1{\sqrt{2}}\begin{pmatrix}1\\1\end{pmatrix}$. If you iterate \begin{align} v^{(k+1)} := K v^{(k)} \end{align} with initial guess $v^{(0)}=\begin{pmatrix}1\\0\end{pmatrix}$ you get the oscillating sequence $v^{(2k)}=v^{(0)}, v^{(2k+1)}=\begin{pmatrix}0\\1\end{pmatrix}$ for $(k=0,1,\ldots)$.

The iterated for the eigenvalues remain zero: $$ \lambda^{(k+1)} = \frac{(v^{(k)})^T K v^{(k)}}{|v^{(k)}|^2}=\frac{(v^{(k)})^T v^{(k)}}{|v^{(k)}|^2}=0 $$

In such cases it is important to have a good initial guess for the eigenvalue and apply the corresponding shift to the matrix $K$.

Since you are looking for the smallest eigenvalue you should apply inverse iteration.

---

You told us that the eigenvalues come in pairs $\pm\lambda$. This looks like there is structure in the problem. If you have some (linear) mapping $f:\mathbb{R}^n\rightarrow\mathbb{R}^n$ which maps the eigen vector $e$ for $\lambda$ to the eigen vector $\bar e$ for $-\lambda$, you could apply a double shift strategy:

$\lambda^{(k)}$ is the current approximation for the eigenvalue, $x^{(k)}$ the current guess for the corresponding eigen vector.

Shift: $\bar K^{(k)} := K+\lambda^{(k)}\m1$

Solve $$ \bar K^{(k)} y^{(k)} = M f(x^{(k)}) $$ for $y^{(k)}$.

With $\tilde K^{(k)} := K-\lambda^{(k)}\m1$ solve $$ \tilde K^{(k)} z^{(k)} = M f(y^{(k)}) $$ for $z^{(k)}$ and set $$ x^{(k+1)}:=\frac{z^{(k)}}{|z^{(k)}|} $$ The new guess for the eigen value should be $$ \lambda^{(k+1)} = \bigl((x^{(k)})^T z^{(k)} + (\lambda^{(k)})^{-2}\bigr)^{-\frac12}. $$

The first guess $x^{(0)}$ should be normalized $|x^{(0)}|=1$.

You can try it with $M=\m1$, $K$ from above and $f(v)=\begin{pmatrix}v_1\\-v_2\end{pmatrix}$.

Answer 3

This is an engineer's solution :) First, transform the problem like I explained in my question so that we (a) take advantage of $\mathbf{K}$ being invertible (while $\mathbf{M}$ is not), and BTW we look for the largest eigenvalue $e'$, which of course will be positive. Then, the smallest positive eigenvalue is $e=1/e'$, and that is what we want! $$ \mathbf{A}\,\mathbf{u} = e'\,\mathbf{u} $$ with $$ \mathbf{A}=\mathbf{K}^{-1}\,\mathbf{M} $$ Now, the power method finds the largest eigenvalue (and its eigenvector) by iterating on: $$ \mathbf{x}^k = \mathbf{A}\,\mathbf{x}^{k-1}\quad\text{for}\quad k=1,2,... $$ But you do not need to invert $\mathbf{K}$. Instead solve $$ \mathbf{LU}\, \mathbf{x}^k = \mathbf{M}\,\mathbf{x}^{k-1} $$ where $\mathbf{LU}$ is the Cholesky factorization of $\mathbf{K}$. Super fast and easy!