# Gathering monomial coefficients of orthogonal polynomials using SpecialPolynomials.jl

**URL:** <https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554>\
**Category:** Performance\
**Tags:** polynomials\
**Created:** [June 25, 2021, 2:00pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554 "2021-06-25T14:00:08Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 25, 2021, 2:00pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/1 "2021-06-25T14:00:09Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 25, 2021, 2:35pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/2 "2021-06-25T14:35:33Z")

</div>

Maybe

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

```

uses the generic methods a bit better.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 25, 2021, 3:04pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/3 "2021-06-25T15:04:44Z")

</div>

> [@BambOoxX](#):
>
> While using `SpecialPolynomials.jl` I want to gather the monomial coefficients of multiple orthogonal polynomials of increasing size in a matrix.

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?

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 25, 2021, 3:09pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/4 "2021-06-25T15:09:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 25, 2021, 3:24pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/5 "2021-06-25T15:24:40Z")

</div>

> [@BambOoxX](#):
>
> 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](http://eprints.ma.man.ac.uk/2281/1/covered/MIMS_ep2015_24.pdf)

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 25, 2021, 4:46pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/6 "2021-06-25T16:46:20Z")

</div>

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](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 !

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2021, 7:52am UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/7 "2021-06-28T07:52:06Z")

</div>

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](https://academic.oup.com/imajna/article/29/1/141/670606?login=true).  
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](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.590.1869&rep=rep1&type=pdf) 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](https://www.sciencedirect.com/science/article/pii/0024379575900415), 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.

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 28, 2021, 11:50am UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/8 "2021-06-28T11:50:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2021, 12:19pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/9 "2021-06-28T12:19:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 28, 2021, 12:28pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/10 "2021-06-28T12:28:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2021, 12:31pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/11 "2021-06-28T12:31:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 28, 2021, 4:09pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/12 "2021-06-28T16:09:14Z")

</div>

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

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

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

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

```

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 29, 2021, 11:08am UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/13 "2021-06-29T11:08:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 29, 2021, 11:16am UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/14 "2021-06-29T11:16:37Z")

</div>

Sure. That would be great.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 29, 2021, 12:37pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/15 "2021-06-29T12:37:56Z")

</div>

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

In principle, yes, but it’s not implemented yet: [eigensolvers for Hessenberg matrices · Issue #41406 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/41406)

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 29, 2021, 12:41pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/16 "2021-06-29T12:41:56Z")

</div>

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 !

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [July 1, 2021, 8:55am UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/17 "2021-07-01T08:55:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [July 1, 2021, 8:56am UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/18 "2021-07-01T08:56:56Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 1, 2021, 12:42pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/19 "2021-07-01T12:42:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 1, 2021, 12:46pm UTC](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554/20 "2021-07-01T12:46:55Z")

</div>

> [@BambOoxX](#):
>
> I just opened the issue. Tell me if you need anything else.

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](https://github.com/jverzani/SpecialPolynomials.jl/issues/23))

[Next page](https://discourse.julialang.org/t/gathering-monomial-coefficients-of-orthogonal-polynomials-using-specialpolynomials-jl/63554.md?page=2)