As I've already mentioned in my previous answer, the method of using the QZ algorithm for matrix pencils on the Frobenius companion linearization of the polynomial eigenproblem is not always the most efficient approach. To illustrate this, I'll outline a general method for solving a hyperbolic quadratic eigenvalue problem, which is known to have all its eigenvalues real. (Encapsulating the strategy as a working Mathematica routine is left as an exercise.)
The method also makes use of a linearization; this linearization, as described in this paper by Higham, Mackey, Mackey, and Tisseur, makes use of block symmetric matrices whose blocks are block Hankel and block antitriangular. A Mathematica routine for constructing the matrix $X_m(P(\lambda))$ (in the notation of section 3.3 of that paper) follows:
BlockSymmetricBasis[k_Integer?NonNegative, matCof : {__?MatrixQ}] :=
Module[{p = Length[matCof] - 1, lm, um},
Switch[k,
0, -ArrayFlatten[Partition[PadRight[Take[matCof, -p], 2 p - 1], p, 1]],
p, ArrayFlatten[Partition[PadLeft[Take[matCof, p], 2 p - 1], p, 1]],
_,
lm = ArrayFlatten[Partition[PadLeft[Take[matCof, k], 2 k - 1], k, 1]];
um = ArrayFlatten[Partition[PadRight[Take[matCof, k - p], 2 (p - k) - 1], p - k, 1]];
ArrayFlatten[{{lm, 0}, {0, -um}}]]] /;
SameQ @@ (Dimensions /@ matCof) && k < Length[matCof]
From here, one now has a family of pencils $X_{m-1}(P(\lambda))-\lambda X_m(P(\lambda))$ at disposal; the key is to choose among these pencils such that $X_m$ is positive definite (i.e., we want to pick out which of the $X_{m-1}(P(\lambda))-\lambda X_m(P(\lambda))$ is a symmetric-definite pencil).
To continue further with the discussion, here is a concrete example of a hyperbolic quadratic eigenvalue problem:
polycof = N@{{{3, 2, 1}, {2, 3, 2}, {1, 2, 3}},
{{-2, -1, -1}, {-1, -3, 2}, {-1, 2, -1}},
{{-5, 1, -2}, {1, -4, -3}, {-2, -3, -5}}};
We check which of the $X_m$ are positive definite:
Table[PositiveDefiniteMatrixQ[BlockSymmetricBasis[k, polycof]], {k, 0, 2}]
{False, True, False}
We thus continue further with the pencil $X_0-\lambda X_1$. Now, we check the condition number $\kappa$ of $X_1$:
X0 = BlockSymmetricBasis[0, polycof]; X1 = BlockSymmetricBasis[1, polycof];
LinearAlgebra`MatrixConditionNumber[X1]
22.2727
The value of the condition number obtained is rather modest in size, so we can continue with one of the usual approaches for symmetric-definite pencils, which uses Cholesky decomposition. First, build the intermediate matrix:
ℳ = CholeskyDecomposition[X1]; lft = LinearSolve[Transpose[ℳ]];
ℋ = lft[Transpose[lft[X0]]];
Here are the eigenvalues:
λ = Eigenvalues[ℋ]
{6.61035, -1.8856, 1.37725, 1.21165, -1.06445, -0.124207}
Check the eigenvalues:
Table[Det[Fold[#1 λ + #2 &, 0, polycof]] // Chop, {λ, %}]
{3.47379*10^-10, 0, 0, 0, 0, 0}
Here are the eigenvectors:
lf = LinearSolve[ℳ];
\[ScriptV] = Take[lf[#], 3] & /@ Eigenvectors[ℋ]
{{-0.46788, 0.903438, -0.600705}, {0.420022, -0.335861, -0.176273},
{-0.425712, -0.113381, 0.331892}, {-0.0699161, -0.208497, -0.204744},
{0.196805, -0.0360833, 0.26503}, {0.069413, 0.115626, -0.103772}}
Check eigenvalues and eigenvectors:
MapThread[Function[{λ, \[ScriptV]}, Chop[Fold[#1 λ + #2 &, 0, polycof].\[ScriptV]]],
{λ, \[ScriptV]}]
{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}
Had the result of LinearAlgebra`MatrixConditionNumber[X1] been rather large (say, $\approx 10^7$), an alternative route would have been the eigendecomposition of X1, which would have gone like this:
ℳ = #2.#1 & @@ MapThread[#1@#2 &, {{Composition[DiagonalMatrix, Sqrt], Transpose},
Eigensystem[X1]}];
lf = LinearSolve[ℳ];
ℋ = lf[Transpose[lf[X0]]]
λ = Eigenvalues[ℋ] (* eigenvalues *)
lft = LinearSolve[Transpose[ℳ]];
\[ScriptV] = Take[lft[#], 3] & /@ Eigenvectors[ℋ] (* eigenvectors *)
See this reference for more information on definite matrix polynomials, which is the general class of matrix polynomials with real eigenvalues.
To summarize: the Frobenius linearization + QZ route works generally, but one should be on the lookout for methods that exploit structure if you will be solving a lot of structured problems, since they will usually be more efficient with space, computational effort, or both.
