Q: extracting coefficients and monomial degrees from multivariate polynomials

[Audrius-St]

Hello,

Please find below a Sympy MWE that describes my question in the title:

import sympy as sp

def main():

    x, y = sp.symbols('x, y')
    a22, a21, a12, a11, a10, a01, a00 = sp.symbols('a22, a21, a12, a11, a10, a01, a00')
    b60, b53, b32, b21, b10, b00 = sp.symbols('b60, a53, a32, b21, b10, b00')

    expr_1 = a22*x**2*y**2 + a21*x**2*y + a12*x*y**2 + a11*x*y + a10*x + a01*y + a00
    expr_2 = b60*x**6 + b53*x**5*y**3 + b32*x**3*y**2 + b21*x**2*y + b10*x + b00

    expr = expr_1 + expr_2

    poly_coeffs = sp.poly(expr, [x, y]).coeffs(order='grevlex')[::-1]

    poly_monoms = sp.poly(expr, [x, y]).monoms(order='grevlex')[::-1]

    print(f'expr = {expr} \n')
    print(f'expr coefficients: {len(poly_coeffs)} {poly_coeffs} \n')
    print(f'expr monomials: {len(poly_monoms)} {poly_monoms} \n')


if __name__ == "__main__":
    main()

gives the results

expr = a00 + a01*y + a10*x + a11*x*y + a12*x*y**2 + a21*x**2*y + a22*x**2*y**2 + a32*x**3*y**2 + a53*x**5*y**3 + b00 + b10*x + b21*x**2*y + b60*x**6 

expr coefficients: 10 [a00 + b00, a01, a10 + b10, a11, a12, a21 + b21, a22, a32, b60, a53] 

expr monomial degrees: 10 [(0, 0), (0, 1), (1, 0), (1, 1), (1, 2), (2, 1), (2, 2), (3, 2), (6, 0), (5, 3)] 

Are similar functions currently available in Symbolics.jl?
A search of Symbolics.jl issues did turn up the following discussion:

#677

However, a search for "coefficient" in the Symbolics.jl documentation only gave a link to Gröbner bases.

[bowenszhu]

using Symbolics
@variables x, y, a22, a21, a12, a11, a10, a01, a00, b60, b53, b32, b21, b10, b00
expr_1 =
    a22 * x^2 * y^2 +
    a21 * x^2 * y +
    a12 * x * y^2 +
    a11 * x * y +
    a10 * x +
    a01 * y +
    a00
expr_2 =
    b60 * x^6 +
    b53 * x^5 * y^3 +
    b32 * x^3 * y^2 +
    b21 * x^2 * y +
    b10 * x +
    b00
expr = expr_1 + expr_2
dict, residual = polynomial_coeffs(expr, (x, y))

result:

julia> dict
Dict{Any, Any} with 10 entries:
  y           => a01
  (x^3)*(y^2) => b32
  x*y         => a11
  (x^5)*(y^3) => b53
  1           => a00 + b00
  x*(y^2)     => a12
  x           => a10 + b10
  y*(x^2)     => a21 + b21
  x^6         => b60
  (x^2)*(y^2) => a22

julia> residual
0

[bowenszhu]

julia> using Symbolics

julia> @doc semipolynomial_form
  semipolynomial_form(expr, vars, degree::Real; consts) -> Tuple{Any, Any}
  

  Returns a tuple of two objects:

    1. A dictionary of coefficients keyed by monomials in vars upto the
       given degree,

    2. A residual expression which has all terms not represented as a
       product of monomial and a coefficient

  degree should be a nonnegative number.

  If consts is set to true, then the returned dictionary will contain a key 1
  and the corresponding value will be the constant term. If false, the
  constant term will be part of the residual.

  semipolynomial_form(exprs::AbstractArray, vars, degree::Real; consts) -> Tuple{Any, Any}
  

  For every expression in exprs computes the semi-polynomial form and returns
  a tuple of two objects – a vector of coefficient dictionaries, and a vector
  of residual terms.

  If consts is set to true, then the returned dictionary will contain a key 1
  and the corresponding value will be the constant term. If false, the
  constant term will be part of the residual.

julia> @doc polynomial_coeffs
  polynomial_coeffs(expr, vars)
  

  Find coefficients of a polynomial in vars.

  Returns a tuple of two elements:

    1. A dictionary of coefficients keyed by monomials in vars

    2. A residual expression which is the constant term

  (Same as semipolynomial_form(expr, vars, Inf))

julia> @doc semilinear_form
  semilinear_form(exprs::AbstractArray, vars) -> Tuple{SparseArrays.SparseMatrixCSC{Num, Int64}, Any}
  

  Returns a tuple of a sparse matrix A, and a residual vector c such that, A * vars + c is the same as exprs.

julia> @doc semiquadratic_form
  semiquadratic_form(exprs, vars) -> Tuple{SparseArrays.SparseMatrixCSC{Num, Int64}, SparseArrays.SparseMatrixCSC{Num, Int64}, Any, Any}
  

  Returns a tuple of 4 objects:

    1. a matrix A of dimensions (m x n)

    2. a matrix B of dimensions (m x (n+1)*n/2)

    3. a vector v2 of length (n+1)*n/2 containing monomials of vars upto degree 2 and zero where they are not required.

    4. a residual vector c of length m.

  where n == length(exprs) and m == length(vars).

  The result is arranged such that, A * vars + B * v2 + c is the same as exprs.

[bowenszhu]

using Symbolics
@variables x, y
expr = (x^(4//3) + y^(5//2))^3
d, r = semipolynomial_form(expr, [y], 5)

result:

julia> d
Dict{Any, Any} with 2 entries:
  y^(5//1) => 3(x^(4//3))
  1.0      => x^(4//1)

julia> r
y^(15//2) + 3(x^(8//3))*(y^(5//2))

[bowenszhu]

using Symbolics
@variables x, y, z
exprs = [3x + tan(z),
    y / z + 5z,
    x * y + y * z / x]
A, c = semilinear_form(exprs, [x, y, z])

result:

julia> A
3×3 SparseArrays.SparseMatrixCSC{Num, Int64} with 2 stored entries:
 3  ⋅  ⋅
 ⋅  ⋅  5
 ⋅  ⋅  ⋅

julia> c
3-element Vector{Num}:
          tan(z)
           y / z
 x*y + (y*z) / x

[bowenszhu]

For more showcases, see
https://github.com/JuliaSymbolics/Symbolics.jl/blob/master/test/semipoly.jl

[bowenszhu]

You need Symbolics v4.10.4 or newer versions

[bowenszhu]

This feature hasn't been added to the documentation of Symbolics.

[Audrius-St]

Thank you for the new Symbolics.jl functions detailed documentation and examples. Much appreciated.

Applied to the initial MWE

# test_poly_coeffs.jl

using Symbolics

begin
    
    @variables x, y
    @variables a22, a21, a12, a11, a10, a01, a00
    @variables b60, b53, b32, b21, b10, b00

    expr_1 = a22*x^2*y^2 + a21*x^2*y + a12*x*y^2 + a11*x*y + a10*x + a01*y + a00
    expr_2 = b60*x^6 + b53*x^5*y^3 + b32*x^3*y^2 + b21*x^2*y + b10*x + b00

    expr = expr_1 + expr_2

    # Polynomial coefficients
    expr_coeffs, _ = polynomial_coeffs(expr, [x, y])

    # Useful to have the polynomial coefficients sorted by the degree of the multivariate terms
    expr_coeffs_sorted = sort(collect(expr_coeffs), by = term->degree(term[1]))

    @show expr_coeffs_sorted

end

gives the result

10-element Vector{Pair{Any, Any}}:
           1 => a00 + b00
           y => a01
           x => a10 + b10
         x*y => a11
     x*(y^2) => a12
     y*(x^2) => a21 + b21
 (x^2)*(y^2) => a22
 (x^3)*(y^2) => b32
         x^6 => b60
 (x^5)*(y^3) => b53

For the second part of my question, how to access the powers of x and y?

{m, n | x^m*y^n where m, n = 0, 1, 2, . . .}

In the case of the MWE

expr monomial degrees: 10 [(0, 0), (0, 1), (1, 0), (1, 1), (1, 2), (2, 1), (2, 2), (3, 2), (6, 0), (5, 3)]

[bowenszhu]

julia> map(expr -> Symbolics.degree.(expr[1], (x, y)), expr_coeffs_sorted)
10-element Vector{Tuple{Int64, Int64}}:
 (0, 0)
 (0, 1)
 (1, 0)
 (1, 1)
 (1, 2)
 (2, 1)
 (2, 2)
 (3, 2)
 (6, 0)
 (5, 3)

Note that degree will probably not be exported in Symbolics after v4.10.4

[Audrius-St]

Great.
Good to know that Symbolics.degree has this functionality.
Thanks again for your explanations and examples.

[xiang-yu]

I am switching to Symbolics.jl now and new to Julia.
How do you @Audrius-St @bowenszhu get all the terms as a vector, i.e., [x2*y2, x2y, xy2, x*y, x, y, 1]?
Same for the expr_coeffs as [ a22, a21, a12, a11, a10, a01, a00].

[Audrius-St]

# test_symbols_vector.jl

using Symbolics

begin
    @variables x, y, x2, y2, x2y, xy2

    # [x2*y2, x2y, xy2, x*y, x, y, 1]
    variables_vec = Vector{Num}(undef, 7)

    # Below you would instead loop overthe 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