I would like to make a multivariate polynomial with given degree.
For example, the bivariate polynomial basis of a vector x = [x1, y2] with degree of 2 would be x -> [1, x1, x1^2, x2, x2^2, x1*x2].

It seems that there are some nice packages for handling multivariate polynomials such as MultivariatePolynomials.jl but I don’t know how to create a polynomial basis functions easily.

Well, I don’t wanna define polyvar objects in this case.

I actually find a way as follows:

using Combinatorics, Transducers
"""
n ∈ N: length of array, i.e., x ∈ Rⁿ
d ∈ N: degree
"""
function polynomial_basis(n, d; with_bias=true)
_n = with_bias ? n+1 : n
exponents = multiexponents(_n, d)
return function (x)
_x = with_bias ? [x..., 1] : x
exponents |> Map(exponent -> prod(_x.^exponent)) |> collect
end
end