Introduction

In numpy the concept of array is generalized to imply arrays of arbitrary dimension, overlapping with the concept of scalars, matrices and tensors. To allow arrays of various dimensions to operate together it defines unambiguous broadcasting rules for what to expect. The results is a library that is used as the reference for almost all of the numerical Python community.

In mathematical literature the term polynomial expansions is used to denote a collection of polynomials. Though they strictly do not need to, they are often indexed, giving each polynomial both a label and a position for where to locate a polynomial relative to the others. Assuming that there always is an index, one could say that polynomial expansions could just as well be termed polynomial array. And using the rules defined in numpy, there no reason not to also start talking about multi-dimensional polynomial arrays.

The main idea here is that in the same way as numpy.ndarray are composed of scalars, numpoly.ndpoly – the baseclass for the polynomial arrays – are composed of simpler polynomials. This gives us a mental model of a polynomial that looks like this:

\[\Phi(q_1, \dots, q_D) = [\Phi_1(q_1, \dots, q_D), \cdots, \Phi_N(q_1, \dots, q_D)]\]

where \(\Phi\) is polynomial vector, \(N\) is the number of terms in the polynomial sum, and \(q_d\) is the \(d\)-th indeterminant name. This mental model is shown in practice in how numpoly displays its polynomials in the REPL:

>>> q0, q1 = numpoly.variable(2)
>>> expansion = numpoly.polynomial([1, q0, q1**2])
>>> expansion
polynomial([1, q0, q1**2])

Here numpoly.variable() creates to simple indeterminants, and the numpoly.polynomial() constructor joins an array of polynomials into a polynomial array, much like numpy.array() does for numeric in numpy.

Another way to look at the polynomials is to keep the polynomial array as a single polynomial sum: A multivariate polynomial can in the case of numpoly be defined as:

\[\Phi(q_1, \dots, q_D) = \sum_{n=1}^N c_n q_1^{k_{1n}} \cdots q_D^{k_{Dn}}\]

where \(c_n\) is a multi-dimensional polynomial coefficients, and \(k_{nd}\) is the exponent for the \(n\)-th polynomial term and the \(d\)-th indeterminant name.

Neither of the two ways of representing a polynomial array is incorrect, and serves different purposes. The former works well for visualisation, while the latter form gives a better mental model of how numpoly handles its polynomial internally.

Modelling polynomials by storing the coefficients as multi-dimensional arrays is deliberate. Assuming few \(k_{nd}\) and large dimensional \(c_n\), all numerical operations that are limited to the coefficients, can be done fast, as numpy can do the heavy lifting.

This way of representing a polynomial also means that to uniquely defined a polynomial, we only need the three components:

  • coefficients – the polynomial coefficients \(c_n\) as multi-dimensional arrays.

  • exponents – the exponents \(k_{nd}\) as a 2-dimensional matrix.

  • indeterminants – the names of the variables, typically q0, q1, etc.

We can access these three defining properties directly from any numpoly.ndpoly polynomial. For example, for a simple polynomial with scalar coefficients:

>>> q0, q1 = numpoly.variable(2)
>>> poly = numpoly.polynomial(4*q0+3*q1-1)
>>> poly
polynomial(3*q1+4*q0-1)
>>> indet = poly.indeterminants
>>> indet
polynomial([q0, q1])
>>> coeff = poly.coefficients
>>> coeff
[-1, 3, 4]
>>> expon = poly.exponents
>>> expon
array([[0, 0],
       [0, 1],
       [1, 0]], dtype=uint32)

Because these three properties uniquely define a polynomial array, they can also be used to reconstruct the original polynomial:

>>> terms = coeff*numpoly.prod(indet**expon, axis=-1)
>>> terms
polynomial([-1, 3*q1, 4*q0])
>>> poly = numpoly.sum(terms, axis=0)
>>> poly
polynomial(3*q1+4*q0-1)

Here numpoly.prod() and numpoly.sum() is used analogous to their numpy counterparts numpy.prod() and numpy.sum() to multiply and add terms together over an axis. See Numpy functions for more details on how this works.

Note

As mentioned the chosen representation works best with relatively few \(k_{nd}\) and large \(c_n\). for large number \(k_{nd}\) and relatively small \(c_n\) however, the advantage disappears. And even worse, in the case where polynomial terms \(q_1^{k_{1n}} \cdots q_D^{k_{Dn}}\) are sparsely represented, the numpoly representation is quite memory inefficient. So it is worth keeping in mind that the advantage of this implementation depends a little upon what kind of problems you are working on. It is not the tool for all problems.