Quadratic eigenvalue problem - biorthogonal system

I need to solve the following quadratic eigenvalue problem: $$\left(-\lambda^2\lbrack \mathbf{X}\rbrack+\text{i}\lambda\lbrack \mathbf{Y}\rbrack+\lbrack \mathbf{Z}\rbrack\right)\left(\mathbf{u}\right)=\left(\mathbf{0}\right)$$ To solve this I used a linearization to transform this quadratic eigenvalue problem into the following generalized eigenvalue problem: $$\lbrack \mathbf{A}\rbrack(\mathbf{u})=\lambda\lbrack \mathbf{B}\rbrack(\mathbf{u})$$ with $$\lbrack \mathbf{A}\rbrack$$ and $$\lbrack \mathbf{B}\rbrack$$ respectively: $$\lbrack \mathbf{A}\rbrack=\begin{bmatrix}\text{i}\lbrack \mathbf{Y}\rbrack&\lbrack \mathbf{Z}\rbrack\\-\lbrack \mathbf{I}\rbrack&\lbrack \mathbf{0}\rbrack\end{bmatrix}~~~~~~~~~\lbrack \mathbf{B}\rbrack=\begin{bmatrix}-\lbrack \mathbf{X}\rbrack&\lbrack \mathbf{0}\rbrack\\\lbrack \mathbf{0}\rbrack&\lbrack \mathbf{I}\rbrack\end{bmatrix}$$ Due to the circumstance that $$\lbrack \mathbf{A}\rbrack$$ and $$\lbrack \mathbf{B}\rbrack$$ aren't hermitian in general in need to solve the adjoint eigenvalue problem $$(\lbrack \mathbf{A}\rbrack^*)^\text{T}(\mathbf{v})=\lambda^*(\lbrack \mathbf{B}\rbrack^*)^\text{T}(\mathbf{v})$$ aswell to receive a biorthogonal system.

• $$\lbrack \mathbf{X}\rbrack, \lbrack \mathbf{Y}\rbrack, \lbrack \mathbf{Z}\rbrack$$ - random $$n\times n$$ matrices
• $$\lbrack \mathbf{A}\rbrack, \lbrack \mathbf{B}\rbrack$$ - corresponding matrices based on the linearization ($$2n\times 2n$$)
• $$\lbrack \mathbf{I}\rbrack$$ - Identity Matrix ($$n\times n$$)
• $$\lambda$$ - eigenvalue
• $$(\mathbf{u}), (\mathbf{v})$$ - right- and left-eigenvectors

This is my implementation in Mathematica:

SeedRandom[123];

n = 15;
XMat = RandomReal[{1*10^-7, 1*10^-6}, {n, n}];
YMat = RandomReal[{1*10^-12, 1*10^-11}, {n, n}];
ZMat = RandomReal[{1*10^-12, 1*10^-11}, {n, n}];

AMat = Join[Join[I*YMat, ZMat, 2], Join[-IdentityMatrix[n], ConstantArray[0, {n, n}], 2]];
BMat = Join[Join[-XMat, ConstantArray[0, {n, n}], 2], Join[ConstantArray[0, {n, n}], IdentityMatrix[n], 2]];

{EWR, EVR} = Eigensystem[{AMat, BMat}];
{EWL, EVL} = Eigensystem[{ConjugateTranspose[AMat], ConjugateTranspose[BMat]}];

perm = Flatten[Nearest[EWL -> "Index", EWR]];
EWR = EWR[[perm]];
EVR = EVR[[perm]];


My main problem is that the corresponding eigenvalues of the mentioned eigenvalue problems are the same and they should be complex conjugated. Did I make a mistake in my Mathematica code or is there a mistake in my theoretical approach?

Update:

I corrected the Mathematica code and changed Transpose to ConjugateTranspose but this didn't solve the problem for all eigenvalues. As shown below most of the corresponding eigenvalues are still identical.

0.00456861 + 3.28732*10^-6 I    0.00456861 - 3.2873*10^-6 I
-0.00456861 + 3.28731*10^-6 I   -0.00456861 - 3.28733*10^-6 I
-0.00257144 + 0.00316279 I      -0.00257543 + 0.00316029 I
0.00257144 + 0.00316279 I       0.00257543 + 0.00316029 I
0.00257543 - 0.00316029 I       0.00257144 - 0.00316279 I
-0.00257543 - 0.00316029 I      -0.00257144 - 0.00316279 I
-0.00372817 + 0.00101919 I      -0.00374614 + 0.00102208 I
0.00372817 + 0.00101919 I       0.00374614 + 0.00102208 I
0.00374614 - 0.00102208 I       0.00372817 - 0.00101919 I
-0.00374614 - 0.00102208 I      -0.00372817 - 0.00101919 I
0.002509 + 0.00236923 I         0.00250875 + 0.00242257 I
-0.002509 + 0.00236923 I        -0.00250875 + 0.00242257 I
0.00250875 - 0.00242257 I       0.002509 - 0.00236923 I
-0.00250875 - 0.00242257 I      -0.002509 - 0.00236923 I
-0.00323767 - 4.43767*10^-6 I   -0.00323767 + 4.43769*10^-6 I
0.00323767 - 4.43767*10^-6 I    0.00323767 + 4.43767*10^-6 I
-0.00237108 - 0.00201776 I      -0.00236452 - 0.00204887 I
0.00237108 - 0.00201776 I       0.00236452 - 0.00204887 I
0.00236452 + 0.00204887 I       0.00237108 + 0.00201776 I
-0.00236452 + 0.00204887 I      -0.00237108 + 0.00201776 I
1.82496*10^-11 - 0.00262931 I   -5.2503*10^-12 - 0.00262943 I
-2.53022*10^-14 + 0.00262943 I  3.85978*10^-17 + 0.00262931 I
0.00196643 + 0.00104558 I       0.00196388 + 0.00105333 I
-0.00196643 + 0.00104558 I      -0.00196388 + 0.00105333 I
0.00196388 - 0.00105333 I       0.00196643 - 0.00104558 I
-0.00196388 - 0.00105333 I      -0.00196643 - 0.00104558 I
-0.000751694 - 0.00133752 I     -0.000759215 - 0.00134487 I
0.000751694 - 0.00133752 I      0.000759215 - 0.00134487 I
0.000759215 + 0.00134487 I      0.000751694 + 0.00133752 I
-0.000759215 + 0.00134487 I     -0.000751694 + 0.00133752 I

• Should Transpose[AMat] be ConjugateTranspose[AMat] perhaps? Commented Oct 16, 2022 at 13:20
• The eigenvalues of a transposed matrix are the same as those from the original matrix. Commented Oct 16, 2022 at 15:00
• Ya this was a mistake but this didn't fix the problem unfortunatly. Most of the eigenvalues are still identical as shown in the updated post. But thank you very much for pointing out my mistake. Commented Oct 16, 2022 at 19:03
• Re your update: I think for each value in the list on the left, its complex conjugate appears in the list on the right. Example: 4th value on the left, its complex conjugate is the 5th value on the right. Commented Oct 16, 2022 at 19:04
• @user293787: Thank you so much for pointing this out to me. I was able to find my mistake. The way I adjusted the ordering between the left and right eigensystem was based on real valued eigenvalues. I was able to solve my problem by changing perm = Flatten[Nearest[EWL -> "Index", EWR]]; to perm = Flatten[Nearest[Conjugate[EWL] -> "Index", EWR]];. You helped me alot! Commented Oct 16, 2022 at 19:18

This answer is similar to the solution you have found but has some advantages. Your equation often occurs with second order differential equations such as the equations for a vibrating system

$$X u''+Y u'+ Z u=0$$

Here the matrices X, Y and Z are the mass, damping and stiffness matrices. In many dynamic cases they are symmetric matrices which yours may or may not be. A good way of organising the eigenvalue problem is to write as a system, as you have done, but in an alternative manner as follows.

$$\left( \begin{array}{cc} -Z & 0 \\ 0 & X \\ \end{array} \right) \left( \begin{array}{c} U' \\ U'' \\ \end{array} \right)+\left( \begin{array}{cc} 0 & Z \\ Z & Y \\ \end{array} \right) \left( \begin{array}{c} U \\ U' \\ \end{array} \right)=\left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right)$$

This form has the advantage of keeping the symmetry if the individual matrices are symmetric. This is important since it avoids left and right eigenvectors.

Here is an example. First make the three individual matrices, here symmetric.

n = 3; (* Number of rows*)
a = RandomReal[{-1, 1}, {n, n}];
x = a + Transpose[a];
a = RandomReal[{-1, 1}, {n, n}];
y = a + Transpose[a];
a = RandomReal[{-1, 1}, {n, n}];
z = a + Transpose[a];
zero = 0 x;


Then we form the symmetric matrices and solve the eigen problem

a = ArrayFlatten[{{-z, zero}, {zero, x}}];
b = ArrayFlatten[{{zero, z}, {z, y}}];
{vals, vecs} = Eigensystem[{b, -a}];


As we should we get diagonal matrices

vecs . a . Transpose[vecs] // Chop // MatrixForm
vecs . b . Transpose[vecs] // Chop // MatrixForm


Hope that helps.

• Thank you very much. This helped alot! Commented Oct 21, 2022 at 13:14