Eigensolver for Complex Companion Matrix

Specific Domains Numerics
#1

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

Thanks,

-Arrigo

#2

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

#3

Have you looked at

3 Likes
#4

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

#5

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

#6

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.

#7

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.

#8

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
#9

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

#10

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
#11

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