@j_verzani @stevengj Iβve done some benchmarks against real and complex scalar polynomials using
roots from Polynomialsroots from AMRVW- generalized companion matrix computation (i.e for non-monic polynomials)
- generalized comrade matrix computation (for polynomial expressed in terms of orthogonal polynomials)
I tested it for polynomials of order 10 and 50 for the moment. Also, as advised by @stevengj , I evaluated the polynomial at the predicted roots to see the residuals.
Depending on the use case, there seem to be large variations in results quality, but I am at the very edge of my understanding of this. Some feedback would be greatly appreciated.
@j_verzani the articles are posted on a public archive, but are under review, I will send you the link ASAP.
The code used for the bechmark is the following
using Polynomials,SpecialPolynomials,LinearAlgebra,AMRVW,Plots
plotlyjs()
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 iblock!(matrix, rowiter, coliter)
@inbounds for k in eachindex(rowiter)
matrix[rowiter[k], coliter[k]] = one(eltype(matrix))
end
return matrix
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 generalized_companion(coeffs::Vector{Matrix{T}}) where T <: Number
n,m = size(coeffs[1])
@assert n == m "Matrix polynomial coefficients should be square matrices"
p = length(coeffs)
B = zeros(T,n*(p-1),n*(p-1))
B[diagind(B)[1:end-n]] .= 1
B[end-n+1:end,end-n+1:end] .= coeffs[end]
A = similar(B)
A[diagind(A,-n)] .= 1
@inbounds for r = 0:p-2
A[r*n.+1:n,end-n+1:end] = -coeffs[r]
end
return A,B
end
function generalized_companion(coeffs::Matrix{T}) where T <: Number
pni,ni = size(coeffs)
p = pni/ni
@assert isinteger(p) "Number of rows of the coefficient matrix should be a multiple of the number of columns"
p = floor(Int64,p)
B = zeros(T,ni*(p-1),ni*(p-1))
B[diagind(B)[1:end-ni]] .= 1
B[end-ni+1:end,end-ni+1:end] .= coeffs[end-ni+1:end,:]
A = similar(B)
A[diagind(A,-ni)] .= 1
A[:,end-ni+1:end] = -coeffs[1:end-ni,:]
return A,B
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
function generalized_comrade(coeffs::Vector{Matrix{T}},OP) where T <: Number
n,m = size(coeffs[1])
@assert n == m "Matrix polynomial coefficients should be square matrices"
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)
diagAa = diagind(A,-n)
diagAb = diagind(A)
diagAc = diagind(A,n)
@views A[diagAa[1:end-n]] .= repeat(an_p,n)
@views A[diagAb[1:end-n]] .= repeat(bn_p[1:end-n],n)
@views A[diagAc[1:end-n]] .= repeat(cn_p[1:end-n],n)
for r = 0:p-2
A[r.+1:n,end-n+1:end] .= .-kn_p[1] .* coeffs[r+1]
end
A[end-2*n+1:end-n,end-n+1:end] .+= kn_p[2] .* cn_p[end] .* coeffs[end]
A[end- n+1:end ,end-n+1:end] .+= kn_p[2] .* bn_p[end] .* coeffs[end]
return A,B
end
function generalized_comrade(coeffs::Matrix{T},OP) where T <: Number
ni = size(coeffs,2)
p = size(coeffs,1)/ni
# Ensure X is of an acceptable format
if !isinteger(p)
error("Coefficient matrix must be of dimensions ni*p Γ ni")
else
p = floor(Int64,p)
end
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])
Cn = @view coeffs[end-ni+1:end,:]
B = zeros(T,ni*(p-1),ni*(p-1))
iblock!(B,1:(p-2)*ni,1:(p-2)*ni)
B[end-ni+1:end,end-ni+1:end] .= kn_p[2]*Cn
A = similar(B)
imat = Matrix(I,ni,ni)
A[1:ni,1:ni] .= bn_p[1]*imat
if p > 2
A[ni+1:2*ni,1:ni] .= an_p[1]*imat
end
shift = 1:ni
for r in 2:p-2
A[ni*(r-2).+shift,ni*(r-1).+shift] .= cn_p[r-1].*imat
A[ni*(r-1).+shift,ni*(r-1).+shift] .= bn_p[r].*imat
A[ni*(r ).+shift,ni*(r-1).+shift] .= an_p[r].*imat
end
@views A[:,end-ni+1:end] .= .-kn_p[1] .* coeffs[1:end-ni,:]
if p > 2
A[end-2*ni+1:end-ni,end-ni+1:end] .+= kn_p[2] .* cn_p[end] .* Cn
end
A[end-ni+1:end ,end-ni+1:end] .+= kn_p[2] .* bn_p[end] .* Cn
return A,B
end
# Define polynomial order
for n in [10,50], t in [Float64,ComplexF64], OP in [Chebyshev,Legendre]
# Generate vector of random polynomial coefficients
X = rand(t,n)
# Case 1 : polynomial is canonical
P = Polynomial(X)
# Default root computation
roots_P = roots(P)
# 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",ms=5,title="Roots_order=$(n)_Canonical_coeffstype=$(eltype(X))")
plot!(roots_amrvw,lt=scatter,label="AMRVW",ms=4)
plot!(roots_gencomp,lt=scatter,label="GenCompanion",ms=3)
savefig("Roots_order=$(n)_Canonical_coeffstype=$(eltype(X)).png")
plot(P.(roots_P),lt=scatter,label="Polynomials",ms=5,title="PVals_order=$(n)_Canonical_coeffstype=$(eltype(X))")
plot!(P.(roots_amrvw),lt=scatter,label="AMRVW",ms=4)
plot!(P.(roots_gencomp),lt=scatter,label="GenCompanion",ms=3)
savefig("PVals_order=$(n)_Canonical_coeffstype=$(eltype(X)).png")
# Generate vector of random polynomial coefficients
X = rand(t,n)
# Case 2 : polynomial is orthogonal polynomial based
P = OP(X)
# Default root computation
roots_P = roots(P)
# Computation with companion matrix
bp = basis_pond(n,OP)
co = bp*P.coeffs
# 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,OP)
roots_gencomr = eigvals!(Matrix(A),Matrix(B))
plot(roots_P,lt=scatter,label="Polynomials",ms=6,title="Roots_order=$(n)_op=$(OP)_coeffstype=$(eltype(X))")
plot!(roots_amrvw,lt=scatter,label="AMRVW",ms=5)
plot!(roots_gencomp,lt=scatter,label="GenCompanion",ms=4)
plot!(roots_gencomr,lt=scatter,label="GenComrade",ms=3)
savefig("Roots_order=$(n)_op=$(OP)_coeffstype=$(eltype(X)).png")
plot(P.(roots_P),lt=scatter,label="Polynomials",ms=6,title="PVals_order=$(n)_op=$(OP)_coeffstype=$(eltype(X))")
plot!(P.(roots_amrvw),lt=scatter,label="AMRVW",ms=5)
plot!(P.(roots_gencomp),lt=scatter,label="GenCompanion",ms=4)
plot!(P.(roots_gencomr),lt=scatter,label="GenComrade",ms=3)
savefig("PVals_order=$(n)_op=$(OP)_coeffstype=$(eltype(X)).png")
end
And the results :
. I also checked using the 
