3
$\begingroup$

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.

$\endgroup$
10
  • 3
    $\begingroup$ In problems like your example the matrices are very sparse and it is usually necessary to calculate only a few of the largest or smallest eigenvalues. In this case, the matlab eigs function is a good choice. $\endgroup$ Jun 28, 2023 at 17:51
  • $\begingroup$ @BillGreene Indeed, the matrices are relatively sparse. With that being said, I do typically want to extract some of the eigenvalues in different manners depending on the context. For example, I have want the ones with largest imaginary part or largest real part (or both). $\endgroup$
    – user45844
    Jun 28, 2023 at 18:10
  • $\begingroup$ @BillGreene when you said to use eig, do you mean the first or second method listed above? The issue with the first method is poor conditioning/large dimensionality leads to me not particularly wanting me to directly invert, but method two takes too long to do the number of runs that I may need. $\endgroup$
    – user45844
    Jun 28, 2023 at 18:12
  • 3
    $\begingroup$ The function I am referring to is eigs not eig. $\endgroup$ Jun 28, 2023 at 18:30
  • 3
    $\begingroup$ Maybe some of you could write an answer rather than a comment, so that this question does not pop up in the unanswered list? $\endgroup$ Jun 30, 2023 at 16:00

1 Answer 1

4
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ I've been playing with the eigs function, but I end up needing to set the tolerance higher than I'd like (~1e-2) in order to get results in a reasonable amount of time, even after making all of the matrices in my original question into sparse ones. How well does eigs perform for poorly-conditioned problems? $\endgroup$
    – user45844
    Jul 3, 2023 at 17:28
  • $\begingroup$ Does it have to be poorly conditioned? I would not expect any iterative methods to be to happy about that. Take the time to rescale the equations to keep things dimensionless and $U \sim I$ if you have not done so. I'll edit my question to mention a bit more about SLEPc. $\endgroup$ Jul 3, 2023 at 20:22
  • $\begingroup$ Unfortunately, the conditioning is poor, as they're very large matrices which come from finite difference approximations. All of the matrices come from discretizing a differential equation, and all of the terms are dimensionless. $U$ is very much not like $I$, though. $\endgroup$
    – user45844
    Jul 3, 2023 at 20:41
  • $\begingroup$ The unknowns can always be rescaled so that $U$'s components are the order of magnitude 1 $\endgroup$ Jul 3, 2023 at 20:52
  • $\begingroup$ Oh, I thought you meant to look like the identity itself. Sure, I can give that a try. Generally, I think that this method should prove useful, thanks for the comprehensive answer! $\endgroup$
    – user45844
    Jul 5, 2023 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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