change monomial basis/coefficients in another basis

[tweisser]

Is there a reason why I cannot find a function like:

using MultivariatePolynomials
const MP = MultivariatePolynomials
using MultivariateBases
const MB = MultivariateBases

function MB.change_basis(p::MP.AbstractPolynomialLike, Basis::Type{<:MB.AbstractPolynomialBasis})
    basis = MB.maxdegree_basis(Basis, variables(p), maxdegree(p))
    coeffs = Float64[]
    rem = p
    for i = 1:length(basis)
        c, rem = divrem(rem, basis.polynomials[i])
        push!(coeffs, c)
    end
    return coeffs, basis
end

which takes a polynomial and returns its coefficients with respect to a certain basis. I can only find this for ScaledMonomialBasis.

using DynamicPolynomials
@polyvar x y
c, m = MB.change_basis(x^2+y, ChebychevBasis)

LinearAlgebra.dot(c, m.polynomials) == x^2+y

[blegat]

This is the role of:

MultivariateBases.jl/src/scaled.jl

Lines 48 to 50 in 17d642d

	function MP.coefficients(p, ::Type{<:ScaledMonomialBasis})
	return unscale_coef.(MP.terms(p))
	end

We should implement it for other basis too. For orthogonal basis, it is probably easier to use the scalar product. With your method, are you sure that rem is zero at the end ?

[tweisser]

I haven't checked, but p is an element of the vector space generated by basis, i.e., there exist unique coefficients p[a] such that p = sum(p[a]*a for a in basis). Now, when I take b in basis of highest degree, c, rem = divrem(p, a) should return rem = sum(p[a]*a for a in basis if !(a==b)) and c = p[b], because b does not divide any of the other a. To see the last assertion, assume b divides a. Since b has maximal degree, a has maximal degree, too, and a = bc for some constant c. This however, contradicts that p has a unique representation p = sum(p[a]*a for a in basis), because

p[b]*b+ p[a]*a = p[b]*b + p[a]*b*c = (p[b]+ p[a]*c)*b + 0*a.

Am I missing something?

Agreed for the scalar product for orthogonal bases though...

[blegat]

No, that's convincing, we should do that for non-orthogonal polynomial basis like FixedPolynomialBasis.

[tweisser]

Coming back to your idea for orthogonal polynomials - we would have to define the scalar product for each basis. How to do this efficiently is actually not 100%clear to me...

[blegat]

Closed by #5