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)
# ]