Eigensolver for Complex Companion Matrix

I’m wondering if there is any Julia eigensolver specialized for complex companion matrices of moderate size, say up to 30x30.

Thanks,

-Arrigo

It would be helpful if you defined what you mean by ‘companion matrices’. Is this what you’re referring to?

Have you looked at

3 Likes

Yes, here is it for clarity:

C=\begin{bmatrix} 0 & 0 & \dots & 0 & -c_0 \\ 1 & 0 & \dots & 0 & -c_1 \\ 0 & 1 & \dots & 0 & -c_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & -c_{n-1} \end{bmatrix}

where c_i \in \mathbb{C}^n

Thanks for the link Sheehan. If this package supports trigonometric polynomials that’s exactly what I’m looking for.

Here’s a toy implementation:

using LinearAlgebra, Polynomials
struct CompanionMatrix{T} <: AbstractMatrix{T}
p::Poly{T}
function CompanionMatrix(p::Poly{T}) where {T}
@assert (coeffs(p))[end] ≈ one(T)
new{T}(p)
end
end
function Base.size(cm::CompanionMatrix)
L = length(coeffs(cm.p)) - 1
(L, L)
end
Base.size(cm::CompanionMatrix, i::Integer) = size(cm)[i]

function Base.length(cm::CompanionMatrix)
(length(coeffs(cm.p)) - 1)^2
end

function Base.getindex(cm::CompanionMatrix{T}, i::Integer, j::Integer) where {T}
if i == j + 1
one(T)
elseif j == size(cm, 1)
coeffs(cm.p)[i]
else
zero(T)
end
end

function LinearAlgebra.eigvals(cm::CompanionMatrix)
roots(cm.p)
end

now at the REPL

julia> cm = CompanionMatrix(Poly([1,2,3,4,5,1]))
5×5 CompanionMatrix{Int64}:
0  0  0  0  1
1  0  0  0  2
0  1  0  0  3
0  0  1  0  4
0  0  0  1  5

julia> eigvals(cm)
5-element Array{Complex{Float64},1}:
-4.192725723692124 + 0.0im
-0.5640990957045247 - 0.3909026378974719im
-0.5640990957045247 + 0.3909026378974719im
0.1604619575505869 - 0.6932715588679903im
0.1604619575505869 + 0.6932715588679903im

and as @dlfivefifty suggests, using FastPolynomialRoots.jl would be advisable as this doesn’t seem to be any faster than using eigvals on collect(cm) for the sizes you care about.

I don’t understand why you say trigonometric polynomials. The eigenvalues of a companion matrix are precisely the roots of a polynomial.

I believe FastPolynomialRoots.jl uses the AMVW algorithm, which is based on orthogonal conjugation of companion matrices.

Alternatively, you can just use a Hessenberg matrix and the relevant LAPACK routine, see

This may be faster at smaller dimensions (the complexity is worse than AMVW though).

1 Like

It may be appropriate to add any of the implementations suggested by @dlfivefifty and @Mason to Companion in SpecialMatrices.jl.

I’m not sure if anything I wrote above is worth contributing to SpecialMatrices.jl, but I hereby give permission to anyone to use any of the code I wrote above in this thread for whatever purpose they desire with no obligation to credit me. I offer no warranty, guarantee or recourse that I didn’t do something stupid or wrong.

2 Likes

I see your point, I just mentioned “trigonometric” polynomials since that’s the kind of polynomials I’m working with.