Usage

Note

Also check out the API documentation or the code.

Let’s look at the example multivariate polynomial:

\(p(x) = 5 + 1 x_1^3 x_2^1 + 2 x_1^2 x_3^1 + 3 x_1^1 x_2^1 x_3^1\)

Which can also be written as:

\(p(x) = 5 x_1^0 x_2^0 x_3^0 + 1 x_1^3 x_2^1 x_3^0 + 2 x_1^2 x_2^0 x_3^1 + 3 x_1^1 x_2^1 x_3^1\)

A polynomial is a sum of monomials. Our example polynomial has \(M = 4\) monomials and dimensionality \(N = 3\).

The coefficients of our example polynomial are: 5.0, 1.0, 2.0, 3.0

The exponent vectors of the corresponding monomials are:

  • [0, 0, 0]

  • [3, 1, 0]

  • [2, 0, 1]

  • [1, 1, 1]

To represent polynomials this package requires the coefficients and the exponent vectors as input.

import numpy as np

coefficients = np.array(
    [[5.0], [1.0], [2.0], [3.0]], dtype=np.float64
)  # numpy (M,1) ndarray
exponents = np.array(
    [[0, 0, 0], [3, 1, 0], [2, 0, 1], [1, 1, 1]], dtype=np.uint32
)  # numpy (M,N) ndarray

Note

by default the Numba jit compiled functions require these data types and shapes

Horner factorisation

to create a representation of the multivariate polynomial \(p\) in Horner factorisation:

from multivar_horner import HornerMultivarPolynomial

horner_polynomial = HornerMultivarPolynomial(coefficients, exponents)

the found factorisation is \(p(x) = x_1^1 (x_1^1 (x_1^1 (1.0 x_2^1) + 2.0 x_3^1) + 3.0 x_2^1 x_3^1) + 5.0\).

pass rectify_input=True to automatically try converting the input to the required numpy data structures and types

coefficients = [5.0, 1.0, 2.0, 3.0]
exponents = [[0, 0, 0], [3, 1, 0], [2, 0, 1], [1, 1, 1]]
horner_polynomial = HornerMultivarPolynomial(
    coefficients, exponents, rectify_input=True
)

pass keep_tree=True during construction of a Horner factorised polynomial, when its factorisation tree should be kept after the factorisation process (e.g. to be able to compute string representations of the polynomials later on)

horner_polynomial = HornerMultivarPolynomial(coefficients, exponents, keep_tree=True)

Note

for increased efficiency the default for both options is False

canonical form

it is possible to represent the polynomial without any factorisation (refered to as ‘canonical form’ or ‘normal form’):

from multivar_horner import MultivarPolynomial

polynomial = MultivarPolynomial(coefficients, exponents)

use this if …

  • the Horner factorisation takes too long

  • the polynomial is going to be evaluated only a few times

  • fast polynomial evaluation is not required or

  • the numerical stability of the evaluation is not important

Note

in the case of unfactorised polynomials many unnecessary operations are being done (internally uses naive numpy matrix operations)

string representation

in order to compile a string representation of a polynomial pass compute_representation=True during construction

Note

the number in square brackets indicates the number of multiplications required to evaluate the polynomial.

Note

exponentiations are counted as exponent - 1 operations, e.g. x^3 <-> 2 operations

polynomial = MultivarPolynomial(coefficients, exponents)
print(polynomial)  # [#ops=10] p(x)


polynomial = MultivarPolynomial(coefficients, exponents, compute_representation=True)
print(polynomial)
# [#ops=10] p(x) = 5.0 x_1^0 x_2^0 x_3^0 + 1.0 x_1^3 x_2^1 x_3^0 + 2.0 x_1^2 x_2^0 x_3^1 + 3.0 x_1^1 x_2^1 x_3^1

horner_polynomial = HornerMultivarPolynomial(
    coefficients, exponents, compute_representation=True
)
print(horner_polynomial.representation)
# [#ops=7] p(x) = x_1 (x_1 (x_1 (1.0 x_2) + 2.0 x_3) + 3.0 x_2 x_3) + 5.0

the formatting of the string representation can be changed with the parameters coeff_fmt_str and factor_fmt_str:

polynomial = MultivarPolynomial(
    coefficients,
    exponents,
    compute_representation=True,
    coeff_fmt_str="{:1.1e}",
    factor_fmt_str="(x{dim} ** {exp})",
)

the string representation can be computed after construction as well.

Note

for HornerMultivarPolynomial: keep_tree=True is required at construction time

polynomial.compute_string_representation(
    coeff_fmt_str="{:1.1e}", factor_fmt_str="(x{dim} ** {exp})"
)
print(polynomial)
# [#ops=10] p(x) = 5.0e+00 (x1 ** 0) (x2 ** 0) (x3 ** 0) + 1.0e+00 (x1 ** 3) (x2 ** 1) (x3 ** 0)
#                   + 2.0e+00 (x1 ** 2) (x2 ** 0) (x3 ** 1) + 3.0e+00 (x1 ** 1) (x2 ** 1) (x3 ** 1)

change the coefficients of a polynomial

in order to access the polynomial string representation with the updated coefficients pass compute_representation=True with in_place=False a new polygon object is being generated

Note

the string representation of a polynomial in Horner factorisation depends on the factorisation tree. the polynomial object must hence have keep_tree=True

new_coefficients = [
    7.0,
    2.0,
    0.5,
    0.75,
]  # must not be a ndarray, but the length must still fit
new_polynomial = horner_polynomial.change_coefficients(
    new_coefficients,
    rectify_input=True,
    compute_representation=True,
    in_place=False,
)

optimal Horner factorisations

use the class HornerMultivarPolynomialOpt for the construction of the polynomial to trigger an adapted A* search to find the optimal factorisation.

See this chapter for further information.

Note

time and memory consumption is MUCH higher!

from multivar_horner import HornerMultivarPolynomialOpt

horner_polynomial_optimal = HornerMultivarPolynomialOpt(
    coefficients,
    exponents,
    compute_representation=True,
    rectify_input=True,
)

Caching

by default the instructions required for evaluating a Horner factorised polynomial will be cached either as .c file or .pickle file in the case of numpy+numba evaluation.

One can explicitly force the compilation of the instructions in the required format:

horner_polynomial = HornerMultivarPolynomial(
    coefficients, exponents, store_c_instr=True, store_numpy_recipe=True
)

If you construct a Horner polynomial with the same properties (= exponents) these cached instructions will be used for evaluation and a factorisation won’t be computed again. Note that as a consequence you won’t be able to access the factorisation tree and string representation in these cases.

the cached files are being stored in <path/to/env/>multivar_horner/multivar_horner/__pychache__/

horner_polynomial.c_file
horner_polynomial.c_file_compiled
horner_polynomial.recipe_file

you can read the content of the cached C instructions:

instr = horner_polynomial.get_c_instructions()
print(instr)

you can also export the whole polynomial class (including the string representation etc.):

path = "file_name.pickle"
polynomial.export_pickle(path=path)

to load again:

from multivar_horner import load_pickle

polynomial = load_pickle(path)

evaluating a polynomial

in order to evaluate a polynomial at a point x:

# define a query point and evaluate the polynomial
x = np.array([-2.0, 3.0, 1.0], dtype=np.float64)  # numpy ndarray with shape [N]
p_x = polynomial(x)  # -29.0

or

p_x = polynomial.eval(x)  # -29.0

or

x = [-2.0, 3.0, 1.0]
p_x = polynomial.eval(x, rectify_input=True)  # -29.0

As during construction of a polynomial instance, pass rectify_input=True to automatically try converting the input to the required numpy data structure.

Note

the default for both options is False for increased speed

Note

the dtypes are fixed due to the just in time compiled Numba functions

evaluating a polynomial at complex query points

x = [np.complex(-2.0, 1.0), np.complex(3.0, -1.0), np.complex(1.0, 0.5)]
p_x = polynomial.eval_complex(x, rectify_input=True)

computing the partial derivative of a polynomial

Note

BETA: untested feature

Note

partial derivatives will be instances of the same parent class

Note

all given additional arguments will be passed to the constructor of the derivative polynomial

Note

dimension counting starts with 1 -> the first dimension is #1!

deriv_2 = polynomial.get_partial_derivative(2, compute_representation=True)
# p(x) = x_1 (x_1^2 (1.0) + 3.0 x_3)

computing the gradient of a polynomial

Note

BETA: untested feature

Note

all given additional arguments will be passed to the constructor of the derivative polynomials

grad = polynomial.get_gradient(compute_representation=True)
# grad = [
#     p(x) = x_1 (x_1 (3.0 x_2) + 4.0 x_3) + 3.0 x_2 x_3,
#     p(x) = x_1 (x_1^2 (1.0) + 3.0 x_3),
#     p(x) = x_1 (x_1 (2.0) + 3.0 x_2)
# ]