Generalized eigenvalue problem for large, potentially ill-conditioned systems

Question

Say that I have a generalized eigenvalue problem of the form $$Ax=\lambda Bx.$$ Using MATLAB, some naive ways that one may solve this is by either

1. directly inverting $B$ then applying the eig function in order to compute the eigenvalues and eigenvectors, or

2. simply applying eig to both matrices i.e. eig(A,B).

The issue with these approaches is the following. Imagine that $A, B$ are large matrices with relatively large condition numbers such as a matrix representing finite difference approximations. Then, the results from method 1) are less trustworthy and prone to error, where results from method 2) take a long time to be generated. For example, on a test problem that I ran, it took over 2 times as long to do method 2), around $65$ vs. $140$ seconds, and the grid resolution was far too coarse for me to realistically use.

What types of methods are best to apply to this problem which are both speedy and accurate — obviously a tall order where some degree of tradeoff must occur? I thought about trying to use an iterative method like multigrid, this isn't in the standard $Ax=b$ form that I'm used to seeing.

For context, this problem actually comes from a quadratic eigenvalue problem $$(\lambda^2 A_2+\lambda A_1+A_0)x=0,$$ but MATLAB's polyeig function takes a very long time due to calculating the QZ factorization (just like how method 2) takes a long time). In this case, the matrix which I called $A$ originally is a block matrix of the form $$\begin{bmatrix} U & V\\ 0 & I\end{bmatrix},$$ and $B$ is of the form $$\begin{bmatrix} 0 & W\\ I & 0\end{bmatrix},$$ where $U,V$, and $W$ involve finite difference matrices.

Answer 1

For large systems, any direct solver methods tend to be a dead end as what often starts as a sparse system ends up becoming dense. In fact, just storing all eigenvectors is itself typically impossibly big.

Some common algorithms for a generalized eigenvalue problem are

1. Power Algorithm
2. Inverse Iteration
3. Subspace Iteration
4. Lanczos Method

all of which come with a lots of options and their own sets of pros and cons (and typically even more options on what to use for the linear solver itself).

In matlab I've fortunately always been well served by the built in eigs which offers a pretty good selection of selection criteria. I'm pretty sure eigs uses the Lanzcos method.

If this isn't enough, I would consider reaching for SLEPc. It seems to no longer offer ready made bindings for MATLAB but writing a MEX function is an option if you are desperate enough.

Unfortunately, it's not trivial to use SLEPc, but it can solve quadratic eigenvalue problems in it's PEP module see example 16

Easiest way to use it is probably to call it from Python https://gitlab.com/slepc/slepc#supported-problem-classes which does cover the quadratic eigenvalue problem, which might work well for you.