Multivariate Polynomial Library

PolyMat is a Python library designed for the representation and manipulation of multivariate polynomial matrices.

Features

• Expression Building: Create polynomial expressions using various operators provided by the library.
• Efficient Internal Representation: Uses a sparse internal structure to optimize intermediate computations.
• Stateful Computation: The sparse internal structures are computed based on a state object. This eliminates any dependency on global variables and provides control over the sparse intermediate structures stored in memory for reuse.
• Performance Optimized: Designed for speed, it outperforms symbolic computation tools like sympy, while remaining lightweidgth and purely written in Python, unlike CasADi or Drake.

Installation

You can install PolyMat via pip:

pip install polymat

Basic Usage

Polynomial expressions are built from constants and variables. The variable $x$ can be defined as follows:

import polymat

x = polymat.define_variable('x')

This variable can then be used to define the polynomial matrix expression:

p = polymat.from_((
    (1 + x, 1 + 2*x + x**2),
    (1, 1 - x**2),
))

The package provides various operators to construct new polynomial expressions from existing ones. For instance, matrix multiplication can be performed as:

ones = polymat.from_(np.ones((2, 1)))

# matrix multiplication
p @ ones

Polynomial expressions in PolyMat stores the necessary information to compute the polynomial monomials and coefficients. Internally, they are represented as tree-like structures, where each node corresponds to an operation (e.g. addition or multiplication) or a terminal value (e.g. a constant or a variable). For example, the polynomial expression $p_{12}(x) = 1 + 2 x + x^2$ can be visualized as the following tree stucture:

graph TD
    A[add] --> B[add]
    A --> C[mult]
    B --> D[1]
    B --> E[mult]
    E --> F[2]
    E --> G[x]
    C --> H[x]
    C --> I[x]

Loading

Operations

PolyMat defines the following operations to create or manipulate polynomial expressions.

Creating Polynomial Expressions

• Polynomial Variable: Define a polynomial variable.

  x = polymat.define_variable('x')

• From Data: Create a polynomial expression from:
  • Tuple of numbers and polynomial variables

    j = polymat.from_(((0, -1), (1, 0)))
    # Matrix([[0, -1], [1, 0]])

  • numpy arrays (possibly containing polynomial variables)

    i = polymat.from_(np.eye(2))
    # Matrix([[1, 0], [0, 1]])

  • sympy expressions (symbols are automatically converted to polynomial variables).

Combining Polynomial Expressions

• Block Diagonal: Combine expression into block diagonal matrices.

  xblk = polymat.block_diag((x, x))
  # Matrix([[x, 0], [0, x]])

• Horizontal Stacking: Stack multiple polynomial expressions horizontally.

  xhstack = polymat.h_stack((x, x))
  # Matrix([[x, x]])

• Vertical Stacking: Stack multiple polynomial expressions vertically.

  xvstack = polymat.v_stack((x, x))
  # Matrix([[x], [x]])

Polynomial Expression Manipulation

• Arithmetic operations: Perform addition (+), subtraction (-), scalar multiplication and division (*, /), matrix multiplication (@), and exponentiation (**).

  f = j * x**2 - i
  # Matrix([[-1.0, -1.0*x**2], [x**2, -1.0]])

• Caching: Cache the polynomial expression to store intermediate results and speed up computation.

  x = x.cache()

• Combinations: Stack multiplied combinations of elements from polynomial vectors. For instance, the vector of monomials $\begin{bmatrix}1&x&x^2\end{bmatrix}$ can be created as follows:

  m = x.combinations((0, 1, 2))
  # Matrix([[1], [x], [x**2]])

• Diagonalization: Extract a diagonal or construct diagonal matrix.

  mdiag = m.diag()
  # Matrix([[1, 0, 0], [0, x, 0], [0, 0, x**2]])

• Differentiation: Compute the Jacobian matrix of a polynomial vector.

  mdiff = m.diff(x)
  # Matrix([[0], [1], [2.0*x]])

• Evaluation: Replace variable symbols within tuple of floats.

  meval = m.eval({x.symbol: (2,)})
  # Matrix([[1], [2.0], [4.0]])

• Kronecker Product: Compute the Kronecker products.

  fkron = f.kron(i)
  # Matrix([[-1.0, 0, -1.0*x**2, 0], [0, -1.0, 0, -1.0*x**2], [x**2, 0, -1.0, 0], [0, x**2, 0, -1.0]])

• Repmat: Repeat polynomial expressions.

  xrepmat = x.rep_mat(3, 1)
  Matrix([[x], [x], [x]])

• Reshape: Modify the shape of polynomial matrices.

  freshape = f.reshape(-1, 1)
  # Matrix([[-1.0], [x**2], [-1.0*x**2], [-1.0]])

• Summation: Sum the rows of the polynomial expression.

  fsum = f.sum()
  # Matrix([[-1.0*x**2 - 1.0], [x**2 - 1.0]])

Specialized methods:

• Monomial Vector: Construct a monomial vector $Z(x)$ appearing in a polynomial expression.

  p = x**3 - 2*x + 3
  
  p_monom = p.monomial_vector(x)
  # Matrix([[1], [x], [x**3]])

• Coefficient Vector: Compute a coefficient matrix $Q$ associated with a vector of monomials $Z(x)$ and a polynomial vector $p(x) = Q Z(x)$.

  p_coeff = p.coefficient_vector(x, monomials=p_monom)
  # Matrix([[3, -2.0, 1]])

Internal Sparse Representation

PolyMat uses an internal sparse polynomial structure to optimize intermediate computations. This internal structure can be accessed by calling the apply method of the polynomial expression using a state object:

state = polymat.init_state()

state, internal_repr = p.apply(state)

# prints the internal sparse structure
# FromPolynomialMatrixImpl(
#   data={
#     (0, 0): {(): 1, ((0, 1),): 1.0}, 
#     (0, 1): {(): 1}, 
#     (1, 0): {(): 1, ((0, 1),): 2.0, ((0, 2),): 1.0},
#     (1, 1): {(): 1, ((0, 2),): -1.0}}, 
#  shape=(2, 2))
print(internal_repr)

The internal sparse structure is built using nested Python dictionaries and tuples, allowing for efficient construction and modifications of sparse representations. While polynomial expressions only maintain a tree-structured representation of how to construct the polynomial matrix, the internal sparse structure explicitly computes and stores the coefficient for each monomial in the matrix. As a result, generating an internal sparse structure is computationally expensive and has been carefully designed to optimize performance.

Output Functions

PolyMat provides output functions that convert polynomial expressions into structured objects for use outside of PolyMat. The following list summarizes these output functions:

• Sympy Representation: Convert a polynomial expression to a sympy expression.

  state, sympy_repr = polymat.to_sympy(f).apply(state)
  # Matrix([[-1.0, -1.0*x**2], [x**2, -1.0]])

• Array Representation: Convert polynomial expressions to an array representation (implemented through numpy and scipy array).

  state, farray = polymat.to_array(f, x).apply(state)
  # {0: array([[-1.], [ 0.], [ 0.], [-1.]]), 2: array([[ 0.], [ 1.], [-1.], [ 0.]])}

• Tuple Representation: Output the constant parts of each component of the polynomial matrix as nested tuples.

  # Setting assert_constant=False will prevent an exception form being raised, even if f is not a constant polynomial expression
  state, ftuple = polymat.to_tuple(f, assert_constant=False).apply(state)
  # ((-1.0,), (-1.0,))

• Polynomial Degrees: Obtain degrees of each component of the polynomial matrix.

  state, fdegree = polymat.to_degree(f).apply(state)
  # ((0, 2), (2, 0))

• Shape of the Matrix: Retrieve the shape of the polynomial matrix.

  state, fshape = polymat.to_shape(f).apply(state)
  # (2, 2)

Array Representation

The first output function polymat.to_array returns an object of type ArrayRepr, a class included in the \textit{PolyMat} package. This array representation serves two purposes. First, it enables efficient evaluation of a polynomial at multiple points. Second, it provides a structured format suitable for input to an SDP solver. For example, when applied to a polynomial vector expression $r(x) \in (R[x])^n$, the resulting ArrayRepr object contains a matrix $R_d$ for each degree $d$ of the polynomial expression, such that

$$r(x) = R_0 + R_1 x + R_2 (x \otimes x) + ...$$

where $R_0 \in \mathbb R^n$, $R_1 \in \mathbb R^{n \times n}$, and $R_2 \in \mathbb R^{n \times n^2}$. To optimize performance, low-degree coefficient matrices ($R_0$ and $R_1$) are stored as dense matrices using numpy, whereas higher degree matrices ($R_2$, ...) are stored as sparse matrices using scipy. The following code snippet demonstrates how to evaluate the polynomial $r(x)$ at different points:

import numpy as np

context, array_repr = polymat.to_array(r).apply(context)

for i in range(5):
    x_eval = np.array(i, 0, 0).reshape(-1, 1)

    # evaluation is implemented by calling the array object
    print(array_repr(x_eval))

Example

The following example computes a vector field defined by the gradient of a scalar polynomial

$$f(x_1, x_2) := (0.1 x_1^2 + 0.2 x_1 x_2 + 0.1 x_2^2 - 1) (x_1 x_2 + x_1 + x_2- 1).$$

The vector field is visualized by first converting the polnyomial expression into its array representation and then evaluating it at each point within the plot domain.

import polymat

# Define polynomial variables and stack them into a vector
names = ('x1', 'x2')
x1, x2 = (polymat.define_variable(n) for n in names)
x = polymat.v_stack((x1, x2))

# Create a polynomial expressions using arithmetic operations
f1 = 0.1*x1**2 + 0.2*x1*x2 + 0.1*x2**2 - 1
f2 = x1*x2 + x1 + x2 - 1
f = f1 * f2

# compute the gradient of the scalar polynomial expression
df = f.diff(x)

# Print the human-readable string representation
# diff(mul(add(add(add(mul(0.1, mul(x1, x1)), mul(mul(0.2, x1), x2)), mul(0.1, mul(x2, x2))), -1), add(add(add(mul(x1, x2), x1), x2), -1)), v_stack(x1,x2))
print(f'{df}')

# Print the internal Python representation of the expression
print(f'{df=}')

# Initialize state
state = polymat.init_state()

# sympy representation
state, sympy_repr = polymat.to_sympy(f).apply(state)
print(f'{sympy_repr}')

# array representation
context, df1_array = polymat.to_array(df[0, 0], x).apply(context)
context, df2_array = polymat.to_array(df[0, 1], x).apply(context)

def plot_vector_field(ax):
    x1 = np.linspace(-2, 2, 20)
    x2 = np.linspace(-2, 2, 20)
    X1, X2 = np.meshgrid(x1, x2)

    def to_array(x1, x2):
        return np.array((x1, x2)).reshape(-1, 1)

    U1 = np.vectorize(lambda x1, x2: df1_array(to_array(x1, x2)))(X1, X2)
    U2 = np.vectorize(lambda x1, x2: df2_array(to_array(x1, x2)))(X1, X2)

    ax.quiver(X1, X2, U1, U2)

References

Here are some references related to the probject:

• State-Monad is a Python library that encapsulates stateful computations into a monadic structure.
• CasADi is a tool for nonlinear optimization and algorithmic differentiation.
• Drake is a tool for model-based design and verification for robotics.
• DynamicPolynomials.jl is a Julia library that provides a sparse dynamic representation of multivariate polynomials.