multivar_horner: a python package for computing Horner factorisations of multivariate polynomials
mmultivar horner: a python package for computing Horner factorisations ofmultivariate polynomials
Jannik MichelfeitTechnische Universitt Dresden [email protected] keywords: multivariate polynomials, Horner factorisations, Horner factorizations,polynomial evaluation, python package, open source software
Abstract
Many applications in the sciences require numerically stable and computationallyefficient evaluation of multivariate polynomials. Finding beneficial representations ofpolynomials, such as Horner factorisations, is therefore crucial. multivar horner [1],the python package presented here, is the first open source software for computingmultivariate Horner factorisations. This work briefly outlines the functionality ofthe package and puts it into reference to previous work in the field. Benchmarksadditionally prove the advantages of the implementation and Horner factorisationsin general.
Polynomials are a central concept in mathematics and find application in a wide range offields. (Multivariate) polynomials have different possible mathematical representationsand the beneficial properties of some representations are in great demand in manyapplications[2]–[4]. 1 a r X i v : . [ c s . M S ] J u l he so called Horner factorisation is such a representation with beneficial properties.Compared to the unfactorised representation of a multivariate polynomial, in the followingcalled canonical form , this representation offers some important advantages. First of allthe Horner factorisation is more compact in the sense that it requires less mathematicaloperations in order to evaluate the polynomial (cf. fig. 4). Consequently, evaluating amultivariate 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 effortrequired 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 factori-sations of multivariate polynomials given in canonical form. Representing multivariatepolynomials of arbitrary degree also in canonical form, computing derivatives of polyno-mials and evaluating polynomials at a given point are further features of the package.Accordingly the package presented here can be helpful always when (multivariate) poly-nomials have to be evaluated efficiently, the numerical error of the polynomial evaluationhas to be small or a compact representation of the polynomial is required. This holdstrue for many applications applying numerical analysis. One example use case wherethis package is already being employed are novel response surface methods[9] based onmultivariate Netwon interploation[4]. In its core multivar horner implements a multivariate Horner scheme with the greedyheuristic presented in [7]. In the following the key functionality of this package is beingoutlined. For a more details on polynomials and Horner factorisations please refer to theliterature, 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 + a x + a x + ... + a d x d (canonical form), the Hornerfactorisation is unique: f ( x ) = a + x ( a + x ( ...x ( a d ) ... ) In the multivariate case howeverthe 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 manypossible Horner factorisations of a multivariate polynomial.Let’s consider the example multivariate polynomial in canonical form p ( x ) = 5 + 1 x x +2 x x + 3 x x x . The polynomial p is the sum of 5 monomials, has dimensionality 3 andcan also be written as p ( x ) = 5 x x x + 1 x x x + 2 x x x + 3 x x x . The coefficients ofthe monomials are 5, 1, 2 and 3 respectively. It is trivial but computationally expensive to2epresent 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 multivariatepolynomials 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 : from m u l t i v a r _ h o r n e r import H o r n e r M u l t i v a r P o l y n o m i a l c o e f f i c i e n t s = [5.0 , 1.0 , 2.0 , 3.0] e x p o n e n t s = [[0 , 0 , 0] , [3 , 1 , 0] , [2 , 0 , 1] , [1 , 1 , 1]] p = H o r n e r M u l t i v a r P o l y n o m i a l ( coefficients , exponents ,r e c t i f y _ i n p u t = True ) x = [ -2.0 , 3.0 , 1.0] p_x = p . eval (x , r e c t i f y _ i n p u t = True ) The factorisation computed by multivar horner is p ( x ) = x ( x ( x (1 x )+2 x )+3 x x )+5 and requires 10 operations for every polynomial evaluation. It should be noted thatthe implemented factorisation procedure is coefficient agnostic and hence does not forexample optimise multiplications with 1. This design choice has been made in order tohave the ability to change the coefficients of a computed polynomial representation aposteriori.With the default settings a Horner factorisation is being computed by recursively fac-torising with respect to the factor most commonly used in all monomials. When noleaves of the resulting binary “Horner Factorisation Tree” can be factorised any more,a “recipe” for evaluating the polynomial is being compiled. This recipe encodes alloperations required to evaluate the polynomial in numpy arrays[11]. This enables the useof functions just in time compiled by Numba [12], which cause the polynomial evaluationto be computationally efficient. The just in time compiled functions are always beingused, since polynomial evaluation in pure python would to some extent outweigh thebenefits of Horner factorisation representations.
It is important to note that in contrast to the one dimensional case, several conceptsof 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 sum3f exponents of all monomials. This is equal to the maximal l -norm of all exponentvectors of the monomials. Accordingly the euclidean degree is the maximal l -norm andthe maximal degree is the maximal l ∞ -norm of all exponent vectors. Refer to [13] forprecise definitions.A polynomial is called fully occupied with respect to a certain degree if all possiblemonomials having a smaller or equal degree are present. The occupancy of a polynomialcan then be defined as the amount of existing monomials relative to the fully occupiedpolynomial of this degree. A fully occupied polynomial hence has an occupancy of 1.Figure 1: the amount of coefficients of fully occupied polynomials of different degrees in3 dimensions.The amount of coefficients (equal to the amount of possible monomials) in multipledimensions highly depends on the type of degree a polynomial has (cf. fig. 1). This effectintensifies as the dimensionality grows. 4igure 2: numerical error of evaluating randomly generated polynomials of varying sizes. For benchmarking our method the following procedure is used: In order to draw poly-nomials with uniformly random occupancy, the probability of monomials being presentis 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 beingactivated with the chosen probability.For each maximal degree up to 7 and until dimensionality 7, 5 polynomials were drawnrandomly. In order to compute the numerical error, each polynomial has been evaluatedat the point of all ones. The true result in this case should always be the sum of allcoefficients. Any deviation of the evaluation value from the sum of coefficients hence isnumerical error. In order to receive more representative results, the obtained numericalerror is being averaged over 100 tries with uniformly random coefficients in the range[ −
1; 1]. 5igure 3: numerical error of evaluating randomly generated polynomials in canonicalform relative to the Horner factorisation.Note that even though the original monomials are not actually present in a Hornerfactorisation, the amount of coefficients however is identical to the amount of coefficientsof its canonical form. With increasing size in terms of the amount of included coefficientsthe 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 canonicalform, the Horner factorisation is more numerically stable as it has also been visualised infig. 3.Even though the amount of operations required for evaluating the polynomials growexponentially with their size irrespective of the representation, the rate of growth islower for the Horner factorisation (cf. fig. 4). As a result, the Horner factorisations arecomputationally easier to evaluate. 6igure 4: amount of operations required to evaluate randomly generated polynomials.
The package has been created due to the recent advances in multivariate polynomialinterpolation[4], [14]. High dimensional interpolants of large degrees create the demandfor evaluating multivariate polynomials computationally efficient and numerically stable.Among others, these advances enable modeling the behaviour of (physical) systems withpolynomials. Obtaining an analytical, multidimensional and nonlinear representationof a system opens up many possibilities. With so called interpolation response surfacemethods [9] for example a system can be analysed and optimised.The commercial software
Maple offers the ability to compute multivariate Horner fac-torisations. multivar horner however is the first open source implementation of amultivariate Horner scheme. The
Wolfram Mathematica framework supports univariateHorner factorisations. The
Julia package
StaticPolynomials has a functionality similarto multivar horner , but does not support computing Horner factorisations.
NumPy [11] offers functionality to represent and manipulate polynomials of dimensionalityup to 3.
SymPy offers the dedicated module sympy.polys for symbolically operating withpolynomials. As stated in the documentation,
SymPy does not use Horner factorisations7or polynomial evaluation in the multivariate case.
Sage covers the algebraic side ofpolynomials. multivar horner has no functions to directly interoperate with other software packages.The generality of the required input parameters (coefficients and exponents) however stillensures the compatibility with other approaches. It is for example easy to manipulatea polynomial with other libraries and then compute the Horner factorisation represen-tation of the resulting output polynomial with multivar horner afterwards, by simplytransferring coefficients and exponents. Some intermediary operations to convert theparameters into the required format might be necessary.
The documentation of the package is hosted on readthedocs.io . Any bugs or featurerequests can be issued on
GitHub . The contribution guidelines can be found there aswell.The underlying basic mathematical concepts are being explained in numerical analysistext books like [10]. The Horner scheme at the core of multivar horner has beentheoretically outlined in [7].Instead of using a heuristic to choose the next factor, one can allow a search overall possible Horner factorisations in order to arrive at a minimal factorisation. Theamount of possible factorisations, however, is increasing exponentially with the degreeand dimensionality of a polynomial (the amount of monomials). One possibility toavoid computing each factorisation is to employ a version of A ∗ search[15] adapted forfactorisation trees. multivar horner also implements this approach, which is similar tothe branch-and-bound method suggested in [16, ch. 3.1].[17] shows how factorisation trees can be used to evaluate multivariate polynomials andtheir derivatives. In [18] Monte Carlo tree search has been used to find more performantfactorisations than with greedy heuristics. Other beneficial representations of polynomialsare for example being specified in [2] and [3].8 Acknowledgements
Thanks to Michael Hecht and Steve Schmerler for valuable input enabling this publica-tion. References [1] J. Michelfeit. (2018). multivar horner , [Online]. Available: https://github.com/MrMinimal64/multivar_horner (visited on Jul. 24, 2020).[2] M. M.-D. Lee, “Factorization of multivariate polynomials”, PhD thesis, TechnischeUniversitt Kaiserslautern, Jun. 13, 2013.[3] C. E. Leiserson, L. Li, M. M. Maza, and Y. Xie, “Efficient evaluation of largepolynomials”, in
International Congress on Mathematical Software , Springer, 2010,pp. 342–353.[4] M. Hecht, K. B. Hoffmann, B. L. Cheeseman, and I. F. Sbalzarini, “MultivariateNewton interpolation”,
ArXiv preprint arXiv:1812.04256 , 2018.[5] J. M. Pea and T. Sauer, “On the multivariate Horner scheme”,
SIAM journal onnumerical analysis , vol. 37, no. 4, pp. 1186–1197, 2000.[6] ——, “On the multivariate Horner scheme II: Running error analysis”,
Computing ,vol. 65, no. 4, pp. 313–322, 2000.[7] M. Ceberio and V. Kreinovich, “Greedy algorithms for optimizing multivariateHorner schemes”, 2004.[8] W. G. Horner, “Xxi. a new method of solving numerical equations of all orders,by continuous approximation”,
Philosophical Transactions of the Royal Society ofLondon , no. 109, pp. 308–335, 1819.[9] J. Michelfeit, “A response surface method based on multivariate Newton interpo-lation”, Master’s thesis, Technische Universitt Dresden, 2019. [Online]. Available: https://mosaic.mpi- cbg.de/docs/Michelfeit2019.pdf (visited on Jul. 24,2020).[10] A. Neumaier,
Introduction to numerical analysis . Cambridge University Press, 2001.[11] S. v. d. Walt, S. C. Colbert, and G. Varoquaux, “The numpy array: A structurefor efficient numerical computation”,
Computing in Science and Engg. , vol. 13, no.2, pp. 22–30, Mar. 2011, issn : 1521-9615. doi : . [Online].Available: http://dx.doi.org/10.1109/MCSE.2011.37 . Max Planck Institute of Molecular Cell Biology and Genetics Helmholtz-Zentrum Dresden-Rossendorf
Proceedings of the Second Workshop on the LLVM Compiler Infrastructurein HPC , ser. LLVM ’15, Austin, Texas: Association for Computing Machinery,2015, isbn : 9781450340052. doi : . [Online]. Available: https://doi.org/10.1145/2833157.2833162 .[13] L. Trefethen, “Multivariate polynomial approximation in the hypercube”, Pro-ceedings of the American Mathematical Society , vol. 145, no. 11, pp. 4837–4844,2017.[14] M. Hecht and I. F. Sbalzarini, “Fast interpolation and Fourier transform in high-dimensional spaces”, in
Intelligent Computing. Proc. 2018 IEEE Computing Conf.,Vol. 2, , K. Arai, S. Kapoor, and R. Bhatia, Eds., ser. Advances in IntelligentSystems and Computing, vol. 857, London, UK: Springer Nature, 2018, pp. 53–75.[15] P. E. Hart, N. J. Nilsson, and B. Raphael, “A formal basis for the heuristicdetermination of minimum cost paths”,
IEEE transactions on Systems Science andCybernetics , vol. 4, no. 2, pp. 100–107, 1968.[16] M. Kojima, “Efficient evaluation of polynomials and their partial derivatives inhomotopy continuation methods”,
Journal of the Operations Research Society ofJapan , vol. 51, no. 1, pp. 29–54, 2008.[17] J. Carnicer and M. Gasca, “Evaluation of multivariate polynomials and theirderivatives”,
Mathematics of Computation , vol. 54, no. 189, pp. 231–243, 1990.[18] J. Kuipers, A. Plaat, J. Vermaseren, and H. J. van den Herik, “Improving multi-variate Horner schemes with Monte Carlo tree search”,