# Eigensolver for Complex Companion Matrix

**URL:** <https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587>\
**Category:** Numerics\
**Created:** [April 14, 2020, 7:56pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587 "2020-04-14T19:56:54Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Arrigo\_Benedetti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arrigo_benedetti/32/25545_2.png) [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)\
**Post date:** [April 14, 2020, 7:56pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/1 "2020-04-14T19:56:54Z")

</div>

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

Thanks,

-Arrigo

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [April 14, 2020, 7:58pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/2 "2020-04-14T19:58:44Z")

</div>

It would be helpful if you defined what you mean by ‘companion matrices’. Is [this](https://en.wikipedia.org/wiki/Companion_matrix) what you’re referring to?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 14, 2020, 8:03pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/3 "2020-04-14T20:03:32Z")

</div>

Have you looked at

> **[GitHub - andreasnoack/FastPolynomialRoots.jl: Fast and backward stable...](https://github.com/andreasnoack/FastPolynomialRoots.jl)**
>
> Fast and backward stable computation of roots of polynomials in Julia - GitHub - andreasnoack/FastPolynomialRoots.jl: Fast and backward stable computation of roots of polynomials in Julia

---

<div class="post-metadata">

**Author:** ![Arrigo\_Benedetti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arrigo_benedetti/32/25545_2.png) [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)\
**Post date:** [April 14, 2020, 8:03pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/4 "2020-04-14T20:03:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Arrigo\_Benedetti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arrigo_benedetti/32/25545_2.png) [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)\
**Post date:** [April 14, 2020, 8:15pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/5 "2020-04-14T20:15:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [April 14, 2020, 8:20pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/6 "2020-04-14T20:20:02Z")

</div>

Here’s a toy implementation:

```julia
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
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.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 14, 2020, 8:26pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/7 "2020-04-14T20:26:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 14, 2020, 8:29pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/8 "2020-04-14T20:29:36Z")

</div>

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

[https://github.com/JuliaApproximation/ApproxFunBase.jl/blob/master/src/LinearAlgebra/hesseneigs.jl](https://github.com/JuliaApproximation/ApproxFunBase.jl/blob/master/src/LinearAlgebra/hesseneigs.jl)

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

---

<div class="post-metadata">

**Author:** ![eliassno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eliassno/32/18917_2.png) [@eliassno](https://discourse.julialang.org/u/eliassno)\
**Post date:** [April 14, 2020, 8:48pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/9 "2020-04-14T20:48:28Z")

</div>

It may be appropriate to add any of the implementations suggested by @dlfivefifty and @Mason to `Companion` in [SpecialMatrices.jl](https://github.com/JuliaMatrices/SpecialMatrices.jl/blob/master/src/companion.jl).

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [April 14, 2020, 8:53pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/10 "2020-04-14T20:53:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Arrigo\_Benedetti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arrigo_benedetti/32/25545_2.png) [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)\
**Post date:** [April 14, 2020, 9:22pm UTC](https://discourse.julialang.org/t/eigensolver-for-complex-companion-matrix/37587/11 "2020-04-14T21:22:01Z")

</div>

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