Dynamic Polynomials

Build Status	References to cite

Sparse dynamic representation of multivariate polynomials that can be used with MultivariatePolynomials (see the documentation there for more information). Both commutative and non-commutative variables are supported. The following types are defined:

• Variable{V,M}: A variable which is commutative with * when V<:Commutative. Commutative variables are created using the @polyvar macro, e.g. @polyvar x y, @polyvar x[1:8] and non-commutative variables are created likewise using the @ncpolyvar macro. The type parameter M is the monomial ordering.
• Monomial{V,M}: A product of variables: e.g. x*y^2.
• MultivariatePolynomials.Term{T,Monomial{V,M}}: A product between an element of type T and a Monomial{V,M}, e.g 2x, 3.0x*y^2.
• Polynomial{V,M,T}: A sum of Term{T,Monomial{V,M}}, e.g. 2x + 3.0x*y^2 + y.

All common algebraic operations between those types are designed to be as efficient as possible without doing any assumption on T. Typically, one imagine T to be a subtype of Number but it can be anything. This is useful for example in the package PolyJuMP where T is often an affine expression of JuMP decision variables. The commutativity of T with * is not assumed, even if it is the coefficient of a monomial of commutative variables. However, commutativity of T and of the variables + is always assumed. This allows to keep the terms sorted (Graded Lexicographic order is used) in polynomial and measure which enables more efficient operations.

Below is a simple usage example

julia> using DynamicPolynomials

julia> @polyvar x y # assigns x (resp. y) to a variable of name x (resp. y)
(x, y)

julia> p = 2x + 3.0x*y^2 + y # define a polynomial in variables x and y
y + 2.0x + 3.0xy²

julia> differentiate(p, x) # compute the derivative of p with respect to x
2.0 + 3.0y²

julia> differentiate.(p, (x, y)) # compute the gradient of p
(2.0 + 3.0y², 1.0 + 6.0xy)

julia> p((x, y)=>(y, x)) # replace any x by y and y by x
2.0y + x + 3.0x²y

julia> subs(p, y=>x^2) # replace any occurence of y by x^2
2.0x + x² + 3.0x⁵

julia> p(x=>1, y=>2) # evaluate p at [1, 2]
16.0

Below is an example with @polyvar x[1:n]

julia> n = 3;

julia> @polyvar x[1:n] # assign x to a tuple of variables x1, x2, x3
(Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}[x₁, x₂, x₃],)

julia> p = sum(x .* x) # compute the sum of squares
x₃² + x₂² + x₁²

julia> subs(p, x[1]=>2, x[3]=>3) # make a partial substitution
13 + x₂²

julia> A = reshape(1:9, 3, 3);

julia> p(x => A * vec(x))  # corresponds to dot(A*x, A*x), need vec to convert the tuple to a vector
194x₃² + 244x₂x₃ + 77x₂² + 100x₁x₃ + 64x₁x₂ + 14x₁²

To declare the created bindings as constants, use constant = true at global scope:

@polyvar x[1:10] constant = true

This option also works with @ncpolyvar and @complex_polyvar, and applies to every binding created by the macro. The macro still returns a tuple of the created bindings. constant must be a literal true or false and defaults to false.

Variable identity and ordering

By default, @polyvar uses variable_order = DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}. Each declaration creates a fresh variable, even when its name is reused:

@polyvar x
old_x = x
@polyvar x
old_x == x # false

Use NameOrder to identify variables by name, independently of creation order. Earlier names in lexicographic order have higher priority (x > y):

using DynamicPolynomials
import DynamicPolynomials as DP

@polyvar y x variable_order = DP.Commutative{DP.NameOrder}
old_x = x
@polyvar x variable_order = DP.Commutative{DP.NameOrder}
@assert old_x == x
@assert old_x + x == 2x
@assert x > y
@assert Dict(old_x => 1)[x] == 1

For noncommutative variables, use @ncpolyvar x variable_order = DP.NonCommutative{DP.NameOrder}. Complex variables, their conjugates, and their real and imaginary parts remain distinct. Use the same variable-order policy for variables combined in a polynomial.

Custom policies can subtype DynamicPolynomials.AbstractVariableOrdering and define DynamicPolynomials.instantiate, Base.cmp, and a matching Base.hash method. Variables compare equal when cmp returns zero, so equal variables must have equal hashes.

variable_order controls the variables; monomial_order independently controls how products of those variables are ordered.

Monomial ordering

The terms of a polynomial are ordered in increasing monomial order. The default ordering is the graded lex order but it can be modified using the monomial_order keyword argument of the @polyvar macro. We illustrate this below by borrowing the example p. 59 of "Ideals, Varieties and Algorithms" of Cox, Little and O'Shea:

julia> p(x, y, z) = 4x*y^2*z + 4z^2 - 5x^3 + 7x^2*z^2
p (generic function with 1 method)

julia> @polyvar x y z monomial_order = LexOrder
(x, y, z)

julia> p(x, y, z) 
4z² + 4xy²z + 7x²z² - 5x³

julia> @polyvar x y z
(x, y, z)

julia> p(x, y, z)
4z² - 5x³ + 4xy²z + 7x²z²

julia> @polyvar x y z monomial_order = Graded{Reverse{InverseLexOrder}}
(x, y, z)

julia> p(x, y, z)
4z² - 5x³ + 7x²z² + 4xy²z

Note that, when doing substitution, it is required to give the Variable ordering that is meant. The default ordering between variables follows their creation order, which can be nondeterministic with parallel computing. Therefore, this order cannot be used for substitution, even as a default (see here for a discussion about this).