Support block companion / comrade matrices ?

[BambOoxX]

This issue is based on the discussion initiated in this Discourse thread

Here are some toy implementations to compare different strategies to compute roots of polynomials expressed in terms of orthogonal or canonical polynomials. Theoretical background is provided in papers cited on the discourse thread (recalled here for convenience)

• A. Amiraslami, R. M. Corless, P. Lancaster. Linearization of matrix polynomials expressed in polynomial bases .IMA. J. Numer. Anal.29(1), pp. 141–157, 2009
• Jared L. Aurentz·Raf Vandebril·David S. Watkins. Fast computation of eigenvalues of companion,comrade, and related matrices. BIT Numer Math (2014) 54:7–30
• Stephen Barnett, A companion matrix analogue for orthogonal polynomials, Linear Algebra and its Applications, Volume 12, Issue 3, 1975, Pages 197-202, ISSN 0024-3795, https://doi.org/10.1016/0024-3795(75)90041-5

using Polynomials,SpecialPolynomials,LinearAlgebra,AMRVW,Plots
plotlyjs()

# Computation of companion matrix
function compmatrix(X::AbstractVector{T}) where T
    @assert X[end] == 1 "Leading coefficient must be unitary, use a comrade matrix otherwise"
    # Get dimensions of X, following the adopted convention
    p = length(X)
    # Initialize the companion matrix
    comp = zeros(T,p-1,p-1)
    # Fill comp with X coefficients
    comp[:,end] .= .-view(X,1:p-1)
    # Fill comp with I block of appropriate dimensions
    iblock!(comp, 2:p, 1:(p-1))
    return comp
end

function generalized_companion(coeffs::Vector{T}) where T <: Number
    p = length(coeffs)
    B = Diagonal([ones(T,p-2);coeffs[end]])
    A = zeros(T,p-1,p-1)
    A[1,end] = -coeffs[1]
    @inbounds for r = 2:p-1
        A[r,r-1] = 1
        A[r,end] = -coeffs[r]
    end
    return UpperHessenberg(A),B
end

function iblock!(matrix, rowiter, coliter)
    @inbounds for k in eachindex(rowiter) 
        matrix[rowiter[k], coliter[k]] = one(eltype(matrix))
    end
    return matrix
end

function basis_pond(p::Int64, OP)
    L = zeros(Float64, p, p)
    for r = 1:p
        P = basis(OP, r-1)
        Q = convert(Polynomial, P)
        L[1:r,r] .= coeffs(Q)
    end
    return UpperTriangular(L)
end

function generalized_comrade(coeffs::Vector{T},OP) where T <: Union{Float64,ComplexF64}
    p = length(coeffs)
    an_p = SpecialPolynomials.an.(OP,collect(0:p-2))
    bn_p = SpecialPolynomials.bn.(OP,collect(0:p-1))
    cn_p = SpecialPolynomials.cn.(OP,collect(1:p-1))
    kn_p = SpecialPolynomials.leading_term.(OP,[p-1,p])

    B = Diagonal([ones(T,p-2);kn_p[2]*coeffs[end]])
    A = zeros(T,p-1,p-1)

    A[1,1] = bn_p[1]
    if p > 2
    A[2,1] = an_p[1]
    end
    @inbounds for r in 2:p-2
        A[r-1,r] = cn_p[r-1]
        A[r  ,r] = bn_p[r]
        A[r+1,r] = an_p[r]
    end
    @views A[:,end] .= .-kn_p[1] * coeffs[1:end-1] 
    if p > 2
    A[end-1,end] = A[end-1,end] + kn_p[2] * cn_p[end] * coeffs[end]
    end
    A[end  ,end] = A[end  ,end] + kn_p[2] * bn_p[end] * coeffs[end]
    return UpperHessenberg(A),B
end

# Define polynomial order
n = 10
# Generate vector of random polynomial coefficients
X = rand(n)

# Case 1 : polynomial is canonical
P = Polynomial(X)
# Default root computation
roots_P = roots(P)
# Computation with companion matrix
co = P.coeffs
co = co./co[end]
roots_comp = eigvals!(compmatrix(co))
# Computation with generalized companion matrix
A,B = generalized_companion(P.coeffs)
roots_gencomp = eigvals!(Matrix(A),Matrix(B))
# Computation with AMRVW
roots_amrvw = AMRVW.roots(co)

plot(roots_P,lt=scatter,label="Polynomials")
plot!(roots_comp,lt=scatter,label="Companion")
plot!(roots_amrvw,lt=scatter,label="AMRVW")
plot!(roots_gencomp,lt=scatter,label="GenCompanion")

# Case 1 : polynomial is Legendre based
P = Legendre(X)
# Default root computation
roots_P = roots(P)
# Computation with companion matrix
bp = basis_pond(n,Legendre)
co = bp*P.coeffs
co = co./co[end]
roots_comp = eigvals!(compmatrix(co))
# Computation with generalized companion matrix
A,B = generalized_companion(co)
roots_gencomp = eigvals!(Matrix(A),Matrix(B))
# Computation with AMRVW
roots_amrvw = AMRVW.roots(co)
# Computation with generalized comrade matrix
A,B = generalized_comrade(P.coeffs,Legendre)
roots_gencomr = eigvals!(Matrix(A),Matrix(B))

plot(roots_P,lt=scatter,label="Polynomials")
plot!(roots_comp,lt=scatter,label="Companion")
plot!(roots_amrvw,lt=scatter,label="AMRVW")
plot!(roots_gencomp,lt=scatter,label="GenCompanion")
plot!(roots_gencomr,lt=scatter,label="GenComrade")

This only tests a comrade implementation versus existing roots evaluation. At the moment, I do not know which way would be best to

• Test a comrade matrix implementation (I do not know if there are analytical results to test against for polynomials with matrix coefficients)
• Ensure the implementation is scalable or at least little-sensitive to numerical error propagation
• Ensure the implementation works well even for complex coefficients (this seems to be somewhat of a difficult case for the default eigensolver)

[jverzani]

No, this is okay to leave open. It turned out that I couldn't leverage so much of this code to implement the Auretz paper, so it made more sense to just add it to SpecialPolynomials.

[BambOoxX]

@jverzani All right. I will add the matrix implementations I got so far tomorrow. As I said earlier, numerical issues seem to occur faster for complex coefficients even for scalar polynomials it seems. I wonder if it is due to the default solver, or else...

[jverzani]

Yes, I've noticed that too when implementing AMRVW. I assumed it was due to more computations creating more opportunity for floating point drift. But it might be something else.

The code in SpecialPolynomials in the PR is such that this might be reverse. The authors note that the single shift case (which is the one complex values use) has better properties than the double shift case.

I'm not sure that handles the block comrade matrices you ulitmately desire though.