multivar_horner: a python package for computing Horner factorisations of multivariate polynomials

Many applications in the sciences require numerically stable and computationally efficient evaluation of multivariate polynomials. Finding beneficial representations of polynomials, such as Horner factorisations, is therefore crucial. multivar_horner, the python package presented here, is the first open source software for computing multivariate Horner factorisations. This work briefly outlines the functionality of the package and puts it into reference to previous work in the field. Benchmarks additionally prove the advantages of the implementation and Horner factorisations in general.


Introduction
Polynomials are a central concept in mathematics and find application in a wide range of fields.(Multivariate) polynomials have different possible mathematical representations and the beneficial properties of some representations are in great demand in many applications [2]- [4].
The so called Horner factorisation is such a representation with beneficial properties.Compared to the unfactorised representation of a multivariate polynomial, in the following called canonical form, this representation offers some important advantages.First of all the Horner factorisation is more compact in the sense that it requires less mathematical operations in order to evaluate the polynomial (cf.fig.4).Consequently, evaluating a multivariate polynomial in Horner factorisation is faster and numerically more stable [5]- [7] (cf.fig.2).These advantages come at the cost of an initial computational effort required to find the factorisation.
The multivar horner python package implements a multivariate Horner scheme ("Horner's method", "Horner's rule") [8] and thereby allows computing Horner factorisations of multivariate polynomials given in canonical form.Representing multivariate polynomials of arbitrary degree also in canonical form, computing derivatives of polynomials and evaluating polynomials at a given point are further features of the package.Accordingly the package presented here can be helpful always when (multivariate) polynomials have to be evaluated efficiently, the numerical error of the polynomial evaluation has to be small or a compact representation of the polynomial is required.This holds true for many applications applying numerical analysis.One example use case where this package is already being employed are novel response surface methods [9] based on multivariate Netwon interploation [4].

Functionality
In its core multivar horner implements a multivariate Horner scheme with the greedy heuristic presented in [7].In the following the key functionality of this package is being outlined.For a more details on polynomials and Horner factorisations please refer to the literature, e.g.[10].
A polynomial in canonical form is a sum of monomials.For a univariate polynomial, which can be written as f (x) = a 0 + a 1 x + a 2 x 2 + ... + a d x d (canonical form), the Horner factorisation is unique: f (x) = a 0 + x(a 1 + x(...x(a d )...) In the multivariate case however the factorisation is ambiguous, as there are multiple possible factors to factorise with.The key functionality of multivar horner is finding a good instance among the many possible Horner factorisations of a multivariate polynomial.
Let's consider the example multivariate polynomial in canonical form p(x) = 5 + 1x 3  1 The polynomial p is the sum of 5 monomials, has dimensionality 3 and can also be written as p The coefficients of the monomials are 5, 1, 2 and 3 respectively.It is trivial but computationally expensive to represent this kind of formulation with matrices and vectors and to evaluate it in this way.In this particular case for example a polynomial evaluation would require 27 operations.Due to its simplicity this kind of formulation is being used for defining multivariate polynomials as input.The following code snipped shows how to use multivar horner for computing a Horner factorisation of p and evaluating p at a point x: The factorisation computed by multivar horner is p(x) = x 1 (x 1 (x 1 (1x 2 )+2x 3 )+3x 2 x 3 )+ 5 and requires 10 operations for every polynomial evaluation.It should be noted that the implemented factorisation procedure is coefficient agnostic and hence does not for example optimise multiplications with 1.This design choice has been made in order to have the ability to change the coefficients of a computed polynomial representation a posteriori.
With the default settings a Horner factorisation is being computed by recursively factorising with respect to the factor most commonly used in all monomials.When no leaves of the resulting binary "Horner Factorisation Tree" can be factorised any more, a "recipe" for evaluating the polynomial is being compiled.This recipe encodes all operations required to evaluate the polynomial in numpy arrays [11].This enables the use of functions just in time compiled by Numba [12], which cause the polynomial evaluation to be computationally efficient.The just in time compiled functions are always being used, since polynomial evaluation in pure python would to some extent outweigh the benefits of Horner factorisation representations.

Degrees of multivariate polynomials
It is important to note that in contrast to the one dimensional case, several concepts of degree exist for polynomials in multiple dimensions.Following the notation of [13] the usual notion of degree of a polynomial, the maximal degree, is the maximal sum of exponents of all monomials.This is equal to the maximal l 1 -norm of all exponent vectors of the monomials.Accordingly the euclidean degree is the maximal l 2 -norm and the maximal degree is the maximal l ∞ -norm of all exponent vectors.Refer to [13] for precise definitions.
A polynomial is called fully occupied with respect to a certain degree if all possible monomials having a smaller or equal degree are present.The occupancy of a polynomial can then be defined as the amount of existing monomials relative to the fully occupied polynomial of this degree.A fully occupied polynomial hence has an occupancy of 1.The amount of coefficients (equal to the amount of possible monomials) in multiple dimensions highly depends on the type of degree a polynomial has (cf.fig.1).This effect intensifies as the dimensionality grows.

Benchmarks
For benchmarking our method the following procedure is used: In order to draw polynomials with uniformly random occupancy, the probability of monomials being present is picked randomly.For a fixed maximal degree n in m dimensions there are (n + 1) m possible exponent vectors corresponding to monomials.Each of these monomials is being activated with the chosen probability.
For each maximal degree up to 7 and until dimensionality 7, 5 polynomials were drawn randomly.In order to compute the numerical error, each polynomial has been evaluated at the point of all ones.The true result in this case should always be the sum of all coefficients.Any deviation of the evaluation value from the sum of coefficients hence is numerical error.In order to receive more representative results, the obtained numerical error is being averaged over 100 tries with uniformly random coefficients in the range [−1; 1].Note that even though the original monomials are not actually present in a Horner factorisation, the amount of coefficients however is identical to the amount of coefficients of its canonical form.With increasing size in terms of the amount of included coefficients the numerical error of both the canonical form and the Horner factorisation found by multivar horner grow exponentially (cf.fig.2).However, in comparison to the canonical form, the Horner factorisation is more numerically stable as it has also been visualised in fig. 3.
Even though the amount of operations required for evaluating the polynomials grow exponentially with their size irrespective of the representation, the rate of growth is lower for the Horner factorisation (cf.fig.4).As a result, the Horner factorisations are computationally easier to evaluate.
Figure 4: amount of operations required to evaluate randomly generated polynomials.

Related work
The package has been created due to the recent advances in multivariate polynomial interpolation [4], [14].High dimensional interpolants of large degrees create the demand for evaluating multivariate polynomials computationally efficient and numerically stable.Among others, these advances enable modeling the behaviour of (physical) systems with polynomials.Obtaining an analytical, multidimensional and nonlinear representation of a system opens up many possibilities.With so called interpolation response surface methods [9] for example a system can be analysed and optimised.NumPy [11] offers functionality to represent and manipulate polynomials of dimensionality up to 3. SymPy offers the dedicated module sympy.polysfor symbolically operating with polynomials.As stated in the documentation, SymPy does not use Horner factorisations for polynomial evaluation in the multivariate case.Sage covers the algebraic side of polynomials.
multivar horner has no functions to directly interoperate with other software packages.
The generality of the required input parameters (coefficients and exponents) however still ensures the compatibility with other approaches.It is for example easy to manipulate a polynomial with other libraries and then compute the Horner factorisation representation of the resulting output polynomial with multivar horner afterwards, by simply transferring coefficients and exponents.Some intermediary operations to convert the parameters into the required format might be necessary.

Further reading
The documentation of the package is hosted on readthedocs.io.Any bugs or feature requests can be issued on GitHub.The contribution guidelines can be found there as well.
The underlying basic mathematical concepts are being explained in numerical analysis text books like [10].The Horner scheme at the core of multivar horner has been theoretically outlined in [7].
Instead of using a heuristic to choose the next factor, one can allow a search over all possible Horner factorisations in order to arrive at a minimal factorisation.The amount of possible factorisations, however, is increasing exponentially with the degree and dimensionality of a polynomial (the amount of monomials).One possibility to avoid computing each factorisation is to employ a version of A * search [15] adapted for factorisation trees.multivar horner also implements this approach, which is similar to the branch-and-bound method suggested in [16, ch. 3.1].[17] shows how factorisation trees can be used to evaluate multivariate polynomials and their derivatives.In [18] Monte Carlo tree search has been used to find more performant factorisations than with greedy heuristics.Other beneficial representations of polynomials are for example being specified in [2] and [3].

Figure 1 :
Figure 1: the amount of coefficients of fully occupied polynomials of different degrees in 3 dimensions.

Figure 2 :
Figure 2: numerical error of evaluating randomly generated polynomials of varying sizes.

Figure 3 :
Figure 3: numerical error of evaluating randomly generated polynomials in canonical form relative to the Horner factorisation.
The commercial software Maple offers the ability to compute multivariate Horner factorisations.multivar horner however is the first open source implementation of a multivariate Horner scheme.The Wolfram Mathematica framework supports univariate Horner factorisations.The Julia package StaticPolynomials has a functionality similar to multivar horner, but does not support computing Horner factorisations.