Julia implementation of the `quadeig` algorithm for the solution of quadratic matrix pencils
Author pablosanjose
Popularity
1 Star
Updated Last
1 Year Ago
Started In
March 2021

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).

### Used By Packages

No packages found.