The QuadEig package implements the `quadeig`

algorithm to deflate zero and infinite
eigenvalues of quadratic pencils. The algorithm is published in:

"An Algorithm for the Complete Solution of Quadratic Eigenvalue Problems" S. Hammarling, C. J. Munro, and F. Tisseur. ACM Trans. Math. Softw. 39, (2013)

Given a quadratic pencil `Q(λ) = A₀ + λ A₁ + λ² A₂`

, where `Aᵢ`

are square matrices of size
`N`

, we want to solve the quadratic right-eigenvalue problem `Q(λ)φ = 0`

. A powerful approach is
to linearize `Q(λ)`

into an equivalent linear pencil `L(λ) = A - λB`

of size `2N`

, and solve
the generalized eigenvalue problem `L(λ)φ´ = (A - λ B)φ´ = 0`

instead. There are many
possible linearizations. The eigenvectors `φ´`

of `L`

will be related to the original `φ`

in
a way that depends on the chosen linearization.

The `quadeig`

algorithm helps to more efficiently solve the `L(λ)φ´ = 0`

problem by
transforming `L(λ)`

into a smaller `L₋(λ) = Q L(λ) V`

with orthogonal `Q`

and `V`

operators.
`L₋(λ)`

shares the same finite eigenvalues as `L`

, but has less (or no) `λ = 0`

and `λ = ∞`

eigenvalues which one wants to discard. This process is called "deflation".

The algorithm relies on the specific structure of the so-called second companion
linearization, defined by matrices `A = [A₁ -I; A₀ 0], B = [-A₂ 0; 0 -I]`

of size `2N`

. The
right-eigenvectors of the original problem `Q`

are obtained from those of `L`

(deflated or
not) by `φ = V * φ´[1:N]`

, where `V`

is the deflation transformation on the right. For
undeflated linearizations, `V`

is not the identity, because a non-deflating transformation
of the second companion linearization is performed for performance reasons.

The QuadEig package exports a `linearize`

function to build `L`

, and a `deflate`

function to
transform an `L`

into a deflated `L₋`

. The `A`

, `B`

and `V`

matrices of a linearization `l`

can be accessed by `l.A`

, `l.B`

and `l.V`

, or through destructuring `A, B, V = l`

.

Example

```
julia> using QuadEig, LinearAlgebra
julia> A₀ = rand(6,6); A₁ = rand(6, 6); A₂ = rand(6, 6);
julia> A₀[:, 1:3] .= A₀[:, 4:6]; # This creates 3 λ = 0 eigenvalues
julia> A₂[:, 2:3] .= A₂[:, 4:5]; # This creates 2 λ = ∞ eigenvalues
julia> l = linearize(A₀, A₁, A₂)
Linearization{T}: second companion linearization of quadratic pencil
Matrix size : 20 × 20
Matrix type : Matrix{ComplexF64}
Scalings γ, δ : (1.0, 1.0)
Deflated : false
julia> eigvals(l.A, l.B) # Note the 3 zero (within machine precision) and 2 infinite (NaN) eigenvalues
12-element Vector{ComplexF64}:
-10.86932670379268 - 7.26632452718938e-15im
-0.9585605368704543 - 1.3660371264720934im
-0.9585605368704528 + 1.3660371264720919im
-0.4087545396926745 + 0.8719854788559378im
-0.4087545396926742 - 0.8719854788559386im
-8.591313021709173e-16 + 0.0im
-7.180754871717689e-17 + 0.0im
2.0338957932144944e-17 - 0.0im
0.32412827713495224 - 4.9149302039555125e-17im
1.281810969835058 + 7.879489919456524e-16im
NaN + NaN*im
NaN + NaN*im
julia> d = deflate(l) # or deflate(A₀, A₁, A₂)
Linearization{T}: second companion linearization of quadratic pencil
Matrix size : 7 × 7
Matrix type : Matrix{ComplexF64}
Scalings γ, δ : (1.0, 1.0)
Deflated : true (12 -> 7)
julia> eigvals(d.A, d.B) # The finite eigenvalues are the same, within machine precision
7-element Vector{ComplexF64}:
-10.869326703792948 - 4.237043003816544e-20im
-0.9585605368704537 + 1.366037126472092im
-0.9585605368704532 - 1.366037126472091im
-0.40875453969267517 + 0.8719854788559379im
-0.40875453969267417 - 0.8719854788559379im
0.32412827713495357 - 0.0im
1.281810969835057 + 1.0657465373859974e-16im
```

The `deflate`

function admits an `atol`

keyword argument to specify a threshold for
eigenvalues to deflate (`|λ| < atol`

for zeros and `|λ| > atol⁻¹`

for infinities).