# Gathering monomial coefficients of orthogonal polynomials using SpecialPolynomials.jl

While using `SpecialPolynomials.jl` I want to gather the monomial coefficients of multiple orthogonal polynomials of increasing size in a matrix.

My current implementation is

``````using LinearAlgebra,SpecialPolynomials,Polynomials
function basis_pond(p::Int64)
L = zeros(Float64,p+1,p+1)
for r = 1:p+1
v = zeros(r)
v[r] = 1
coeffs = convert(Polynomial,MonicLaguerre{0}(v)).coeffs
L[1:r,r] .= coeffs
end
return UpperTriangular(L)
end
``````

Is there a better way to access these coefficients ?

Maybe

``````p = basis(MonicLaguerre{0}, r)
q = convert(Polynomial, p)
L[1:r,r] .= coeffs(q)
``````

uses the generic methods a bit better.

Note that expressing polynomials in a monomial basis is often a disaster in finite precision. For example, you may get very ill-conditioned matrices. Is this just for educational purposes, or are you trying to actually solve some practical problem, and if so what is your underlying problem?

1 Like

Yes, this is sort of the end issue with this strategyâ€¦

The practical problem is that I have to compute zeros of a matrix valued polynomial (the coefficients of the polynomial are matrices, the variable is a scalar).

I know that I can compute the zeros using a companion matrix under the assumption that

• the coefficients are associated with a canonical polynomial
• the polynomial is monic

Basically, what I canâ€™t figure otherwise is : how to get the zeros of my polynomial expressed in an orthogonal basis.

There are â€ścompanion matrixâ€ť (or â€ścolleague matrixâ€ť or â€śconfederate matrixâ€ť) methods that operate directly on the coefficients of the orthogonal basis. See, for example: http://eprints.ma.man.ac.uk/2281/1/covered/MIMS_ep2015_24.pdf

1 Like

Well, I guess it is a bit late on a friday to anlyze this completely, but it feels like a great starting point.
I am currently trying out multiple basis for my proble so I found this paper that may be of help for any kind of orthogonal basis : https://epubs.siam.org/doi/pdf/10.1137/040609847. Would you care to have a quick look at it @stevengj ? I am not sure I can discriminate good or bad strategies on such topic. Thanks anyway for the paper !

1 Like

I guess I found exactly what I need in ref 2 of the paper you linked : 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.
Another paper 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 seems to cover this as well.
Side note, here is an accessible paper (on science direct though) from Stephen Barnett about this 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, which is probably a seminal paper on the topic.

Thanks again for the help ! I have to see if and where I can contribute once I finish the implementation.

1 Like

For the Aurentz et. al. paper, the AMRVW.jl package has a lot of the necessary plumbing implemented.

My end goal is to work on block-comrade matrices, is it supported ?

Not immediately. Iâ€™ll have a look this week to see if it can be incorporated.

All right, thanks. I started a small implementation based on SpecialPolynomials.
The idea would be to pass two arguments to the function

• A vector of coefficients (scalar or matrices)
• An orthogonal polynomial of SpecialPolynomials

The output would be either a comrade matrix if the leading term is identity, or two matrices forming the generalized comrade problem otherwise.

Here is a quick implementation based on the paper by Amiraslami, Corless and Lancaster

``````using LinearAlgebra, SpecialPolynomials
function generalized_comrade(coeffs::Vector{T},OP) where T <: Number
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))

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

A[1,1] = bn_p[1]
A[2,1] = an_p[1]
@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]
A[end-1,end] += kn_p[2] * cn_p[end] * coeffs[end]
A[end  ,end] += kn_p[2] * bn_p[end] * coeffs[end]
return UpperHessenberg(A),B
end
``````

There is however a non-negligible difference between this implementation and what is found by back substitution and companion matrix computation (which is what is currently done by SpecialPolynomials AFAICT).
I canâ€™t figure if this comes from a bug on my part (I checked again,but I may have missed something) or from the back substitution.

See :

``````using Polynomials
coeffs = rand(6)
L = Legendre(coeffs)
P = convert(Polynomial,L)

ref_roots = roots(P)

comp_roots = eigvals(Matrix(A),Matrix(B)) # Don't know if diagonal and hessenberg structure can be leveraged there)
``````

EDIT, I checked for Chebyshev polynomials, and the results check out w.r.t the analystical expression of the Chebyshev nodes

``````using Plots
n = 30
Cnodes(n) = cos.((2 .* (1:n) .- 1)./(2 .* n) .* pi)
C = basis(Chebyshev,n)
r = roots(C)
rr = eigvals(Matrix(A),Matrix(B))
rel(a,b) = abs(a-b)/abs(b)
plot(rel.(-r,Cnodes(n)))
plot!(rel.(-rr,Cnodes(n)))
``````

Do you want me to open an issue with the developments Iâ€™ve done so far ?

Sure. That would be great.

In principle, yes, but itâ€™s not implemented yet: eigensolvers for Hessenberg matrices Â· Issue #41406 Â· JuliaLang/julia Â· GitHub

I see, thanks for the interest anyway ! I donâ€™t know if you checked but I managed the computations thanks to the information you provided. I still do not understand why the results differ when using a Legendre basis but otherwise it works fine !

1 Like

I just opened the issue. Tell me if you need anything else.

When I switch to complex valued polynomial coefficients, the eigensolver seems to have much more troubleâ€¦ Would you know a way to circumvent this? Because of that, for my end application back substitution works better than using comrade matrices at the momentâ€¦

You have to be more specific. What do you mean by â€śmore troubleâ€ť? (What do you mean by â€śback substitution?â€ť Conversion back to a monomial basis and using a companion matrix?)

If you have code, ideally you can open a pull request. (There already seems to be an issue comrade matrix Â· Issue #23 Â· jverzani/SpecialPolynomials.jl Â· GitHub)