Does anyone know how to get coefficients for ALL combination of the variables of a multivariable polynomial using sympy.jl or other Julia package for symbolic computation?

Here is an example from MATLAB,

syms a b y
[cxy, txy] = coeffs(ax^2 + by, [y x], ‘All’)
cxy =
[ 0, 0, b]
[ a, 0, 0]
txy =
[ x^2y, xy, y]
[ x^2, x, 1]

As you can see it returns all combination of the variables with 0 coefficients.

I know the Python package Sympy has a function “all_coeffs”, but it cannot handle multivariable polynomials and is not available at sympy.jl as far as I know.

The choice of appropriate CAS (computer algebra system) really depends on the type of math you will be working with. Nemo/Oscar in Julia is good for abstract math projects, Symbolics (or Sympy) is useful for generic symbolic manipulations over real/complex numbers.

In Nemo you can do what you want as follows:

julia> using Nemo # or the basic, unoptimized, pure-julia AbstractAlgebra
julia> R, (x,y) = PolynomialRing(RealField(), ["x", "y"])
(Multivariate Polynomial Ring in x, y over Real Field with 256 bits of precision and error bounds, AbstractAlgebra.Generic.MPoly{arb}[x, y])
julia> f = x*x*3+3*y*x;
julia> [coeff(f, [x,y], [i,j]) for i in 0:2 for j in 0:2]
9-element Vector{AbstractAlgebra.Generic.MPoly{arb}}:
0
0
0
0
3.00000000000000000000000000000000000000000000000000000000000000000000000000000
0
3.00000000000000000000000000000000000000000000000000000000000000000000000000000
0
0

You can use a different field, to specify Rational numbers, or only Integers, or finite precision Floats. Here I am using arbitrary precision BigFloats (the default in RealField).

Thank you for the prompt and informative reply.
I am dealing with derivatives, getting coefficients and terms and other manipulations of a 500*500 symbolic matrix. The MATLAB code is way too slow for “gradient”, “coeffs”, and for loop for such a big symbolic matrix. That is why I am shifting to sympy.jl or symengine.jl.

I see that Nemo can deal with coefficients for multivariable polynomials. How do you get ALL the terms x^3y, x^2, xy, x, y, 1 using Nemo?
Also, based on your experience, would Nemo be the candidate for the manipulation of 500*500 symbolic matrix with polynomial elements? The MATLAB couldn’t return a result after 24 hours wall-clock time with 64 CPUs.

If you are going to do code generation out of that matrix, I would probably suggest Symbolics.jl.

If you know that all you are going to do is polynomials, especially to high orders and of not too many variables, and only derivatives on them, then Nemo.jl probably would be the fastest. Especially if you know the polynomials would only be over integers or rationals.

To get all terms, you probably need to make your own function using the degrees function to find the highest term in the polynomial and then the coeff function to extract all the coefficients.

If you use Symbolics.jl instead, the API is probably different.

I am still not sure what “getting all the terms” mean. Given f = 2*x^3 + x^2 * y + y, what do you want it to return? Unfortunately I don’t speak MATLAB. Do you just want the list [2*x^3, x^2 * y, y]? Or do you want the coefficients of f for some fixed list of monomials like [1, y, x, x*y, x^2, x^2*y]? If so, for which monomials?

Nemo is happy with 500x500 matrices, depending on what one wants to do.

Only now I am noticing that you seem to want polynomials with symbolic coefficients. This is not really the type of task most “abstract math” algebra systems are best at: for them polynomials are over a specific ring of coefficients, not over wildcard symbols. That enables all the interesting fast mathematical operations (e.g. in cryptography) that people care about.

If you just want expressions over many symbols where only addition, multiplication, and exponentiation by an integer are used, then Symbolics is probably going to be easier to use. But that might not actually use a special-purpose “polynomial” datastructure, rather a generic “symbolic expression tree” datastructure.

The whole “symbolic coefficients in my polynomial” kinda throws a wrench in the common tools (at least if you want to be fast – if you do not care about speed, there is always a solution). Maybe you can share a bit more about the actual problem you want to solve and why you need symbolic coefficients in your polynomials?

I don’t know how well this scales, but if you want to try SymPy for this task, I think this is doing what you want:

using SymPy
@syms x[1:10]
Σ = zero(Sym)
for i ∈ 1:9
Σ += i * x[i] * x[i+1] # some arbitrary polynomial expression
end
p = Σ.as_poly()
xs = p.gens
ds = p.as_dict()
d = Dict( prod(filter(!iszero, k .* xs)) => v for (k, v) ∈ ds)

I know the Python package Sympy has a function “all_coeffs”, but it cannot handle multivariable polynomials

Sympy.py can handle multivariate polynomials in the manner you require.
Please refer to the link below for an example.

and is not available at sympy.jl as far as I know.

Have used Sympy.py for a number of years, recently switched to the all-Julia package Symbolics.jl for symbolic computation. Speed of execution was one deciding factor.

Indeed, I want polynomials with symbolic coefficients. For simplicity, let me share a TOY version of my problem. Considering solving an equation of dx/dt = Ax, where x is a symbolic vector with polynomial elements, x = [a1x1^2+a2x2^2+a3x3^2+a4x1x2+a5x2x3+1, b1x1^2+b2x2^2+b3x3^2+b4x1x2+b5x2x3+1,
c1x1^2+c2x2^2+c3x3^2+c4x1x2+c5x2*x3+1 ]

First, I need to get ALL the symbolic coefficients, a1,…a5, b1,…b5, and c1, …, c5 and the terms x1^2x2^2x3^2, …, 1.

Then I define another vector y=[y1, y2, y3, …,y9] such that y1=x1^2, y2=x2^2, y3=x3^2, y4=x1*x2, …, y9=1.

I construct a new linear equation of y using the chain rule , dy/dt=(dy/dx)*(dx/dt) = By, where B = (dy/dx).

My goal is to obtain a symbolic matrix B, this means I need to derive dy/dx symbolically. The main operations are calculating the symbolic derivatives and obtaining the coefficients.

My real problem deals with a much larger x and higher order.
The parallel MATLAB is way too slow, so I am switching to Julia symbolic computation. @Krastanov, @thofma, speed (parallelisation) is exactly what I need.

Thanks. I am also experimenting Symbolics.jl. Speed and parallelisation are exactly what I need. I’ve now posted my exact problem below (sorry I don’t know how to point you to my post, hope you can see it), any tips will be highly appreciated.

Thanks. This seems to leave out coefficients of first order, i.e., x₁, x₁₀, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉.
I need to try out the scaling. The speed is what I am after.

Because of the need for symbolic coefficients, I would suggest using Audrius-St suggestion above with Symbolics.jl. Its general purpose symbolic expression structure would not be as fast as the polynomial structures from Nemo, but it would probably be necessary in order to deal with the unconstrained coefficients.

My coefficients are all symbolic. Hope that works with PolynomialRing. I’ve used sympy.jl to derive my vector f. I guess I need to convert f to PolynomialRing and then perform collect(coefficients(f)) and collect(exponent_vectors(f)).

I am switching to Symbolics.jl now.
In your GitHub post,
How do you get all the terms as a vector, i.e., [x2*y 2, x2y, xy 2, x*y, x, y, 1]?
Same for the expr_coeffs as [ a22, a21, a12, a11, a10, a01, a00].
I’ve got my matrix using Julia and now trying to get the symbolic coefficients and the corresponding terms in the form of a vector, respectively.

# test_symbols_vector.jl
using Symbolics
begin
# [x2*y2, x2y, xy2, x*y, x, y, 1]
@variables x, y, x2, y2, x2y, xy2
variables_vec = Vector{Num}(undef, 7)
# Below you would instead loop over the indices of the matrix
variables_vec[1] = x2*y2
variables_vec[2] = x2y
variables_vec[3] = xy2
variables_vec[4] = x*y
variables_vec[5] = x
variables_vec[6] = y
variables_vec[7] = 1
@show variables_vec
end

Output:

variables_vec = Num[x2*y2, x2y, xy2, x*y, x, y, 1]
7-element Vector{Num}:
x2*y2
x2y
xy2
x*y
x
y
1

Thanks. I just figured it out. The terms and coefficients can be collected from the returned dictionary “coeffs” from
Julia> coeffs, _ = polynomial_coeffs()