Gathering monomial coefficients of orthogonal polynomials using SpecialPolynomials.jl

General Usage Performance
29

I did some quick tests, and get the following:

julia> DataFrame(case1.(T, ns))  # T=Float64; ns = [10,25,30]
3×6 DataFrame
 Row │ n      η              roots        proots       amrvw        gcomp       
     │ Int64  Float64        Float64      Float64      Float64      Float64     
─────┼──────────────────────────────────────────────────────────────────────────
   1 │    10    107.518      5.71565e-15  7.75204e-16  4.82674e-15  5.0533e-15
   2 │    25  11169.5        4.17834e-13  3.2743e-14   1.03356e-13  3.64146e-13
   3 │    50      1.05848e8  1.58139e-9   6.67526e-11  1.23539e-10  2.50662e-10

All roughly the same, though PolynomialRoots a tad better (smaller maximum polynomial value over the identified roots).

julia> DataFrame(case2.(T, ns))
3×4 DataFrame
 Row │ n      η             conv         cmp         
     │ Int64  Float64       Float64      Float64     
─────┼───────────────────────────────────────────────
   1 │    10    57.6049     1.94289e-14  7.56062e-14
   2 │    25  2548.86       2.95504e-9   2.6897e-6
   3 │    50     4.78515e7  4.40678      0.226592

As anticipated in this thread, conversion from the the orthogonal basis to the standard basis loses accuracy compared to the companion matrix approach

julia> DataFrame(case3.(T, ns))
3×6 DataFrame
 Row │ n      η              avw          avw_✓        conv         cmp         
     │ Int64  Float64        Float64      Float64      Float64      Float64     
─────┼──────────────────────────────────────────────────────────────────────────
   1 │    10    160.907      2.93324e-11  4.00999e-17  9.55756e-10  1.57792e-10
   2 │    25  77960.7        2.45473e-13  4.4044e-15   2.74531e-12  9.89575e-12
   3 │    50      1.7372e10  3.4879e-14   2.98266e-14  1.06142e-9   2.04649e-8

The AVW algorithm in the SpecialPolynomials PR performs better here.

These runs are typical, but a fair amount of the time the random nature produces poorly behaved polynomials for the larger n values

The above is based on the following slight modifications to your example:

import PolynomialRoots
using DataFrames
# should use a different norm...
function condition_number(p, z)
    n = degree(p)
    sqrt(n) * (1 + norm(z)^2)^((n-2)/2) * norm(coeffs(p), 1)/norm(derivative(p)(z))
end

## Case 1
function roots_generalized_companion(ps)
    A,B = generalized_companion(ps)
    eigvals!(Matrix(A),Matrix(B))
end

function case1(T,n)
    P = Polynomial
    ps = rand(T,n)
    p = P(ps);

    (n=n,
     η =maximum(z -> condition_number(p, z), roots(p)),
     roots=maximum(norm∘p, roots(p)),
     proots=maximum(norm∘p, PolynomialRoots.roots(ps)),
     amrvw = maximum(norm∘p, AMRVW.roots(ps)),
     gcomp=maximum(norm∘p, roots_generalized_companion(ps))
     )
end

function case2(T,n)
    OP = Legendre
    ps = rand(T, n);
    p = OP(ps); # roots(p) will fall back to roots(convert(Polynomial, p))
    bp = basis_pond(n,OP);
    co = bp*p.coeffs;

    (n=n,
     η = maximum(z -> condition_number(p, z), roots(p)),
     conv = maximum(norm∘p, roots(p)),
     cmp= maximum(norm∘p, roots_generalized_companion(co))
     )
end

function case3(T,n)
    OP = MonicLegendre
    ps = rand(T, n);
    p = OP(ps); # roots(p) using PR
    q = convert(Polynomial, p);

    bp = basis_pond(n,OP);
    co = bp*p.coeffs;

    (n=n,
     η = maximum(z -> condition_number(p, z), roots(p)),
     avw=maximum(norm∘p, roots(p)),
     avw_✓=maximum(SpecialPolynomials.a_posteriori_check(roots(p), p)),
     conv=maximum(norm∘p, roots(q)),
     cmp=maximum(norm∘p, roots_generalized_companion(co))
     )
end