Eigenvalues of a linear matrix pencil (Ax = λ Bx) using the subspace iteration method

Question

I am trying to find the rightmost eigenvalues of the generalized eigenvalue problem ($Ax = \lambda Bx$) using the subspace iteration method. This formulation arises from flow stability analysis where the matrix 'B' is singular and both A and B are sparse matrices of the order of 50,000. I am basically trying to check if the matrix pencil has eigenvalues with a positive real part (which means that the flow is unstable to small perturbations), and when it does, I need to know the imaginary part of these eigenvalues as it gives the frequency of the oscillations.

The idea is to use power iterations (with orthonormalization) to approximate a subset of the eigenvectors (T). Using these vectors as the basis for an invariant subspace, reduced subspace matrices (A' = T^TA T and B' = T^TB T) are formed. Then, this reduced dimension eigenvalue problem is solved directly using the QZ algorithm.

The algorithm goes like this -
Step 1: Generate a random set of $m$ orthonormal vectors ($R^{(0)}$), where m is greater than $p$, the number of eigenvalues sought and much less than $n$, the order of the matrices
Step 2: Calculate iteratively $R^{(i+1)} = (A^{-1} B) R^{(i)}$ by solving $AR^{(i+1)} = BR^{(i)}$ and normalize $R$ after every iteration
Step 3: When $\kappa(R)$ increases above a threshold, perform a QR decomposition with $T = R^{new} = Q$
Step 4: Calculate the subspace matrices $A' = T^TA T$ and $B' = T^TB T$
Step 5: Find the eigenvalues and eigenvectors of $A'y = \lambda B'y$ using QZ algorithm
Step 6: Recalculate $R$ using $R = Ty$ and go back to Step 2.

This process shall return the smallest eigenvalues of the matrix pencil. In order to get the rightmost eigenvalues, I pre-condition the matrix using shift-invert transformation.

The entire algorithm seems to be perfect in theory, but some reflection upon this process made me wonder that the subspace spanned by the eigenvectors of the matrix pencil will not actually be invariant under A and B independently! I consider this to be the reason why the eigenvalues don't converge when I try a sample problem using my code. However, when I choose B to be the identity matrix, I do get the correct eigenvalues.

Is my reasoning correct, or is there some other source of error in this algorithm?
Upon realizing this issue, I had modified the algorithm with $C = A^{-1}B$ calculated in advance and substituted at all places (It makes $C$ singular but does provide the eigenvalues). However, I am still curious to know if I have identified the issue correctly, and why at all does the algorithm work (it's mentioned in literature), if it is flawed.

Thanks in advance!

Note: I am an Engineering undergrad and I studied about invariant subspaces just for the project that I am working on.

Answer 1

i) $\det(A-\lambda B)=0$ can be rewritten $\det(\dfrac{1}{\lambda}I-A^{-1}B)=0$. You calculate the largest eigenvalue of $A^{-1}B$ which is the smallest generalized eigenvalue of the pencil $A-\lambda B$.

ii) You calculate $A^{-1}$ (what is the time of calculation ?).

When a matrix is large, it is often ill-condditionned; what is $CN(A)$, the condition number of $A$? If $CN(A)\approx 10^k$, then you must consider $k$ supplementary significant digits along the calculations. Then, calculate $CN(A)$ (there are approximation methods of $CN$ that have complexity $O(n^2)$) before calculating $A^{-1}$.

iii) When you calculate $A^{-1}Bx$, dont calculate $(A^{-1}B)x$, but calculate $A^{-1}(Bx)$, or better, solve in $y$ the equation $Ay=Bx$ (here, you don't need $A^{-1}$).

EDIT. Answer to the OP.

i) Perhaps you assume that $A,B$ are real symmetric (??); that does not imply that the eigenvalues of $A^{-1}B$ are real! -except if $A>0$-

ii) If I understand correctly, you obtain (after some iterations) an orthonormal basis of a limit subspace $E$ that is invariant for $A^{-1}B$. Firstly, $E$ is not invariant by $A,B$. Secondly, how do you choose the dimension of $E$ ? Thirdly, I don't understand why you would get such a limit (especially if there are very close eigenvalues); but let's admit that.

You must write the matrix of $A^{-1}B$ in $E$ -using your orthonormal basis- and diagonalize this matrix. In particular you must calculate explicitly $A^{-1}B$ and you don't work about $A,B$.

iii) Normally, you should only work on $\mathbb{R}$ -the imaginary parts are parasitic; remove them-

About the digits -if you work on the real numbers-

complex*16 works with $16.8=128$ bits, that is, with $38$ significand digits. Because of $CN(A)$, it remains $38-19=19$ significand digits, that suffices -except if you use complex numbers-