A misunderstanding or a bug in LAPACK's solver for generalized eigenvalue problems?

Question

In my application, I have two general real matrices $A$,$B$ defined as follows, $$ A=\begin{bmatrix} -s I_3 & A_0 & 0 & 0 \\ A_0^T & -s I_3 & 0 & 0 \\ 0 & A_1 & -s I_3 & A_0 \\ A_1^T & 0 & A_0^T & -s I_3 \end{bmatrix}, B=\begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ I_3 & 0 & 0 & 0 \\ 0 & I_3& 0 & 0 \end{bmatrix}, $$ where $A_0,A_1 \in R^{3\times3}$ are general matrices, $s$ is a real scalar, and $I_3 \in R^{3\times3}$ is the identity matrix. I would like to solve the eigenvalue problem $AX= B X D$ numerically, where $X$ is the eigenvector matrix and $D$ is the eigenvalue matrix. I only want the right (not left) eigenvalues and eigenvectors. Consider the following as a concrete example, where $$ A_0=\begin{bmatrix} 1 & & \\ & 2 & \\ & & 3 \end{bmatrix}, A_1=\begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}, s=1, $$ From my application, I know there exists at least one eigenvalue $d_1$ and one eigenvector $x_1$, which are $$ x_1=\left(0.625543,0,0,0.625543, 0, 0, 0, -0.417029, 0, 0, -0.208514, 0\right), d_1=0, $$ because $ A x_1 -d_1 B x_1 = \textbf{0},\left\|x_1\right\|=1$. However, no matter which software I used, I couldn't get the correct eigenvalues, even for $d_1$. The software I have tried are Matlab, Eigen(c++ library) and LAPACK.

1. The eigenvalues from Matlab are infinity.
2. The right eigenvalues from Eigen and LAPACK are represented by $d_i=\alpha_i / \beta_i$, where $\alpha_i$ is a complex number and $\beta_i$ is a real number. The output $\beta_i$ is zero for all eigenvalues. In other words, all right eigenvalues are infinity.

I'm not sure which part is wrong. Is it that the software cannot compute a correct eigenvalue or that $d_1$ cannot be considered as an eigenvalue?

Crosspost my own question: A misunderstanding or a bug in LAPACK's solver for generalized eigenvalue problems?, as it may be a better fit here.

Answer 1

I suspect the root of your trouble is what has been detected in the comments by Vibe: For any number $\omega\in \mathbb{K}$ (with $\mathbb{K}= \mathbb{R}$ or $\mathbb{C}$) you can find $\boldsymbol{X}$ such that $AX = \omega BX$ (with $A$ and $B$ taken in your concrete example).

You have already decomposed the problem in 4 blocks of 3 variables. Then let us define $X = (X_1, X_2, X_3, X_4)^T$ where $X_i\in \mathbb{K}^3$.

The first 6 lines of your system yield $X_1 = X_2 = (a,0,0)$, with $a\in\mathbb{K}$ an arbitrary number.

With this in mind, it becomes obvious upon inspection that the last 6 lines of the system are of the form $C (X_3, X_4)^T = R(a,\omega)$ where $C$ is a $6\times6$ non singular matrix (essentially the bottom right blocks in $A$), and $R\in\mathbb{K}^6$ is a right-hand side that depends on $a,\omega\in\mathbb{K^2}$. Since $C$ is invertible, you can always find solutions to this non homogeneous problem.

In conclusion, for all $\omega$, you can invert this system and compute $X\ne 0$ such that $AX=\omega BX$

Answer 2

The other answers already tell you what went wrong, but I will add a terminology note: the term for what is happening is that the pencil $A - \lambda B$ is a singular matrix pencil, i.e., $\det (A - \lambda B)$ is identically equal to zero. So there are no generalized eigenvalues (or, at least, they cannot be defined as usual as the roots of the generalized characteristic polynomial).

An useful tool to analyze singular pencils is the Kronecker canonical form, which is a generalization of the Jordan canonical form to pencils and can include singular (rectangular) blocks.

Note that small perturbations of a singular pencil can give regular (=non-singular) pencils with any eigenvalue in $\mathbb{C} \cup \{\infty\}$, so numerically computing those eigenvalues is an ill-posed problem, and in practice you can expect just about any value to show up among the computed eigenvalues.