How should I solve generalized eigenvalue problems in Python? (Orr-Sommerfeld equation)

Question

I am trying to solve the Orr-Sommerfeld equation numerically, using the techniques given in this article. This leads to solving a generalized eigenvalue problem, that is, given two matrices $\mathbf A,\mathbf B$ with possibly complex entries, we would like to find a pair $(\lambda,\boldsymbol v)$ where both $\lambda $ and $\boldsymbol v$ are possibly complex, satisfying $$\mathbf A\boldsymbol v=\lambda \mathbf B\boldsymbol v\tag 1$$

Note that $\mathbf B$ is not necessarily invertible, i.e, we cannot (in general) simply multiply both sides by $\mathbf B^{-1}$ to reduce it to an ordinary eigenvalue problem. And, even if we could, these matrices $\mathbf A,\mathbf B$ are fairly large, say about $100\times 100$ (and $\mathbf B$ is, in practice, probably quite sparse), so this would probably be a waste of computational resources anyway. (I would also like to add that in my application I am more interested in the eigenvalue $\lambda$ than I am in the eigenvector $\boldsymbol v$.)

A common way of solving this kind of generalized eigenvalue problem is with generalized Schur decomposition, also known as "QZ decomposition" in which we write

$$\mathbf A=\mathbf {QS}\mathbf Z^* \\ \mathbf B=\mathbf {QT}\mathbf Z^*\tag {2}$$ Where $\mathbf Q,\mathbf Z$ are unitary and $\mathbf S,\mathbf T$ are upper triangular. (The asterisk denotes the conjugate transpose.) The eigenvalues of $(1)$ then can be given as the ratio of the diagonal entries of these two upper triangular matrices, $$\lambda_i=\frac{T_{ii}}{S_{ii}}$$

---

Solving the problem numerically in python.

Python's scipy package has the module linalg.eig , which, according to the documentation, is able to solve generalized eigenvalue problems of the form of $(1)$. The right-hand-side matrix can be given in the optional argument b. The default value of this argument b is None, in which the module will solve the standard eigenvector problem, i.e $\mathbf B=\mathbf I$ in $(1)$.

However, the scipy package also includes the linalg.qz module, which, according to the documentation is able to perform a decomposition of the form of $(2)$.

My question, is, would I be better off simply inputting my matrices $\mathbf A,\mathbf B$ as the a and b arguments into the scipy.linalg.eig module, or would I be better off first finding their Schur decompositions with scipy.linalg.qz and then inputting the matrices $\mathbf S,\mathbf T$ into scipy.linalg.eig separately, and setting the argument b=None? Which one is faster?

Thanks.

Answer 1

The performance of the two methods should be very close, because scipy.linalg.eig internally calls Lapack's ggev, which itself calls hgeqz to compute a QZ decomposition. So, typically, most of the CPU time will be spent inside the same hgeqz function anyway.

There are two main differences:

1. ggev uses a form of deflation that might increase the accuracy of the computed quantities, so it is better to use it rather than unscaled QZ. [point EDITed after a correction by thijssteel]

2. Some additional work is required to extract the eigenvectors from the QZ decomposition. If you only need the eigenvalues, you can skip that work by calling scipy.linalg.eigvals instead (or eig with right=0, which is the same thing).

Answer 2

You seem to have a number of misconceptions about numerical linear algebra. The biggest is 100x100 isn't large. Large matrices start somewhere in the 10000x10000 range. The second thing is that even though the inverse doesn't exist, you can still do a linear solve (and preserve sparsity). Specifically, you can just compute the eigenvalues of B\A.