Gathering monomial coefficients of orthogonal polynomials using SpecialPolynomials.jl

General Usage Performance
1

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 ?

2

Maybe

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

uses the generic methods a bit better.

3

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
4

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.

5

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
6

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 : The Society for Industrial and Applied Mathematics. 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
7

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
8

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

9

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

10

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

11

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.

12

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))
    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]
    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)

A,B = generalized_comrade(coeffs,Legendre)
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)
A,B = generalized_comrade(C.coeffs,Chebyshev)
rr = eigvals(Matrix(A),Matrix(B))
rel(a,b) = abs(a-b)/abs(b)
plot(rel.(-r,Cnodes(n)))
plot!(rel.(-rr,Cnodes(n)))
13

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

14

Sure. That would be great.

15

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

16

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
17

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

18

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…

19

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?)

20

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)