Polynomial Bases

The foundation of this library is a set of polynomial basis classes. Each class defines a set of functions \(\lbrace \phi_j \rbrace_j\) that form a basis for the desired polynomial space; i.e., a mononomial basis for polynomials of degree 3 of one variable has basis functions

\[\phi_0(x) = 1, \quad \phi_1(x) = x, \quad \phi_2(x) = x^2, \quad \phi_3(x) = x^3.\]

Here we consider two types of polynomial bases: total degree polynomial bases and maximum degree polynomial bases:

\[\begin{split}\mathcal P_{m}^{\text{tot}} & := \text{Span} \left\lbrace f: f(\mathbf x) = \prod_{i=1}^d x_i^{\alpha_i} \right\rbrace_{|\boldsymbol \alpha| \le m}, \quad & \text{where} & \quad |\boldsymbol \alpha| = \sum_{i=1}^d \alpha_i; \\ \label{eq:max} \mathcal P_{\mathbf m}^{\text{max}} & := \text{Span} \left \lbrace f: f(\mathbf x) = \prod_{i=1}^d x_i^{\alpha_i} \right\rbrace_{ \boldsymbol \alpha \le \mathbf m}, & \text{where} & \quad \boldsymbol \alpha \le \mathbf m \Leftrightarrow \alpha_i \le m_i.\end{split}\]

Note

Although these polynomial bases ideally represent the same quantities, depending on the input coordinates their numerical conditioning can be very different. The most robust approach, although expensive, is to use ArnoldiPolynomialBasis. If performance is desired, choose a tensor product polynomial basis for which the input coordinates (after scaling) approximately sample the measure under which the polynomials are orthogonal. As a rule of thumb:

• if \(x_i\) is uniformly distributed on \([a,b]\) use LegendrePolynomialBasis;
• if \(x_i\) is normally distributed with mean \(\mu\) and standard deviation \(\sigma\) use HermitePolynomialBasis;
• if \(\mathbf x_i \in \mathbb{C}^d\) for \(d\ge 2\) and degree greater than 10, use ArnoldiPolynomialBasis.

Most bases are initialized using a set of input coordinates \(\mathbf x_i\) and a degree. We require the input coordinates to determine a scaling to make the basis well conditioned; i.e., for the Legendre polynomial basis we perform an affine transformation (which does not alter polynomial basis) such that these input coordinates are in the \([-1,1]^m\) hypercube.

class polyrat.PolynomialBasis(X, degree)

An abstract base class for polynomial bases.

__init__(X, degree)

Parameters:

• X (array_like) – Array of dimensions (M, dim) specifying discrete points on which the polynomial basis will be evaluated.
• degree (int or list of int) –

  The degree of the polynomial:

  • if an int, this specifies a total degree polynomial;
  • if a list of int, this specifies a maximum degree polynomial (must match the number of columns in X)

X

The points on which the basis was defined.

degree

The degree of the polynomial.

dim

The dimension of the input space.

mode

What type of polynomial basis is used

roots(coef)

Compute the roots of a univariate, scalar valued polynomial.

Given coefficients \(c_j\) that define a polynomial

\[p(x) = \sum_{j=1}^N c_j \phi_j(x)\]

compute the roots \(r_k\) such that \(p(r_k) = 0\).

Parameters:coef (array-like, (N,)) – Coefficients \(c_j\)

vandermonde(X)

Construct the generalized Vandermonde matrix associated with the polynomial basis.

For a given input X whose rows are the vectors \(\mathbf x_i\), construct the generalized Vandermonde matrix whose entries are

\[[\mathbf V]_{i,j} = \phi_j(\mathbf x_i).\]

Parameters:X (array_like) – Coordinates at which to evaluate the basis; of dimensions (M, dim). Returns:V – generalized Vandermonde matrix; of dimensions (M, N) Return type:ndarray

vandermonde_X

Alias for vandermonde(X) where X is the points where the basis was initialized

Returns:V – generalized Vandermonde matrix of dimensions (M, N) evaluated using the points the basis was initialized with Return type:ndarray

vandermonde_derivative(X)

Construct the derivative of the generalized Vandermonde matrix.

For a given input X whose rows are the vectors \(\mathbf x_i\), construct the 3-tensor whose entires correspond to the derivative of the basis polynomial \(\phi_j\) with respect to each coordiante evaluated at \(\mathbf x_i\).

\[[\mathbf V]_{i,j,k} = \left. \frac{\partial}{\partial x_k} \phi_j(\mathbf x) \right|_{\mathbf x = \mathbf x_i}.\]

Parameters:X (array-like (M, dim)) – Coordinates at which to evaluate the basis; of dimensions (M, dim). Returns:DV – slice-wise derivative of generalized Vandermonde matrix; of dimensions (M, N, dim) Return type:ndarray

Tensor Product Bases

Several polynomial bases are constructed from univariate polynomial bases defined in Numpy numpy.polynomial. These bases are each constructed by subclassing polyrat.TensorProdutPolynomialBasis which performs two tasks:

• multiplying univariate polynomials in each coordinate to construct the desired multivariate basis and
• using an affine transformation in each coordinate to improve conditioning.

Here we provide an overview of the internals of this class to show how new bases could be implemented in the same way.

class polyrat.TensorProductPolynomialBasis(X, degree)

Abstract base class for polynomial bases constructed from univariate bases

_vander(X, degree)

Construct a univariate Vandermonde matrix.

When subclassing, this function should implement the same API as functions like polyvander().

Parameters:

• X (array_like) – A one-dimensional array of coordinates
• degree (int) – Nonnegative integer for the degree

Returns:

V – Array of shape (len(X), degree+1).

Return type:

ndarray

_der(coef)

Given univariate polynomial coefficients, yield the coefficients of the deriative polynomial.

When subclassing, this function should implement the same API as functions like polyder().

Parameters:coef (array_like) – Polynomial coefficients Returns:der – Coefficients of deriative polynomial Return type:ndarray

roots(coef)

Given univariate polynomial coefficients, compute the roots in this basis.

When subclassing, this function should implement the same API as functions like polyroots().

Parameters:coef (array_like) – Polynomial coefficients Returns:roots – Roots of the polynomial Return type:ndarray

Notes

Although computing the roots is more a property of a polynomial, and not a basis, how we compute these roots depends strongly on the choice of basis.

_set_scale()

Perform any precomputation needed to scale the basis.

To better condition the polynomial, we perform an affine transformation on the input data; i.e.,

\[\hat{\mathbf{x}}_i = \text{diag}(\mathbf a) \mathbf x_i + \mathbf b.\]

By default we scale the basis such that \(\hat{\mathbf{x}}_i \in [-1,1]^d\). However other bases may implement a different scaling.

_scale(X)

Scale coordinates using an affine transform

Parameters:X (ndarray) – Coordinates to transform Returns:Xhat – Transformed coordinates Return type:ndarray

_inv_scale(X)

Perform the inverse of the affine scaling

Parameters:Xhat (ndarray) – Transformed coordinates Returns:X – Coordinates in original space Return type:ndarray

Monomial Basis

class polyrat.MonomialPolynomialBasis(X, degree)

A tensor product polynomial bases constructed from monomials

Univariate monomials take the form of

\[\begin{split}\phi_0(x) &= 1 \\ \phi_1(x) &= x \\ \phi_2(x) &= x^2 \\ \phi_3(x) &= x^3 \\ \vdots & \phantom{=}\vdots \\ \phi_k(x) &= x^k\end{split}\]

To improve conditioning we scale the input to the \([-1,1]^d\) hypercube.

Legendre Basis

class polyrat.LegendrePolynomialBasis(X, degree)

A tensor product polynomial basis constructed from Legendre polynomials

Legendre polynomials are an orthogonal basis on the interval \([-1,1]\):

\[\begin{split}\phi_0(x) &= 1 \\ \phi_1(x) &= x \\ \phi_2(x) &= \frac12 (3x^2 -1) \\ \phi_3(x) &= \frac12 (5x^3 -3x) \\ \vdots & \phantom{=}\vdots \\ \phi_k(x) &= \frac{1}{2^k k!} \frac{\partial^k}{\partial x^k} (x^2 - 1)^k\end{split}\]

Chebyshev Basis

class polyrat.ChebyshevPolynomialBasis(X, degree)

A tensor product polynomial basis constructed from Chebyshev polynomials

Chebyshev polynomials are an orthogonal basis on the interval \([-1,1]\) with the weight \((1 - x^2)^{-1/2}\):

\[\begin{split}\phi_0(x) &= 1 \\ \phi_1(x) &= x \\ \phi_2(x) &= 2x^2 - 1 \\ \phi_3(x) &= 4x^3 - 3x \\ \vdots & \phantom{=}\vdots \\ \phi_k(x) &= 2x \phi_{k-1}(x) - \phi_{k-2}(x)\end{split}\]

Hermite Basis

class polyrat.HermitePolynomialBasis(X, degree)

A tensor product polynomial basis constructed from Hermite polynomials

Hermite polynomials are an orthogonal basis on the interval \((-\infty, \infty)\) with weight \(e^{-x^2/2}\). Following hermite(), we use physicists Hermite polynomials

\[\begin{split}\phi_0(x) &= 1 \\ \phi_1(x) &= 2x \\ \phi_2(x) &= 4x^2 - 1 \\ \phi_3(x) &= 8x^3 - 12x \\ \vdots & \phantom{=}\vdots \\ \phi_k(x) &= (-1)^k e^{x^2} \frac{\partial^k}{\partial x^k} e^{-x^2}.\end{split}\]

Here we scale the coordinates such that they have zero mean and unit variance to improve the conditioning of the generalized Vandermonde matrix.

Laguerre Basis

class polyrat.LaguerrePolynomialBasis(X, degree)

A tensor product polynomial basis constructed from Laguerre polynomials

Laguerre polynomials are an orthogonal basis on the interval \([0, \infty)\) with weight \(e^{-x}\).

\[\begin{split}\phi_0(x) &= 1 \\ \phi_1(x) &= -x+1 \\ \phi_2(x) &= \frac12(x^2 - 4x + 1) \\ \phi_3(x) &= \frac16(-x^3 + 9x^2 - 18x + 6) \\ \vdots & \phantom{=}\vdots \\ \phi_k(x) &= \sum_{\ell=0}^k {k \choose \ell} \frac{ (-1)^\ell}{\ell!} x^\ell.\end{split}\]

Here we scale the coordinates to best approximate a unit exponential distribution.

Arnoldi Polynomial Basis

class polyrat.ArnoldiPolynomialBasis(X, degree, weight=None)

A polynomial basis constructed using Vandermonde with Arnoldi.

A polynomial basis constructed with Vandermonde with Arnoldi; see [BNT19x], [AKL+20] for discussion.

Notes

• In order to compute the roots in the univariate case we convert to Legendre polynomial.

For low-level access, the following functions are available.

polyrat.vandermonde_arnoldi_CGS(X, degree, weight=None, mode=None)

Multivariate Vandermode with Arnoldi using classical Gram-Schmidt with reorthogonalization

This function uses the Arnoldi proceedure to construct an orthogonal basis for polynomials. In the univariate case, this corresponds to constructing a Krylov subspace for the starting vector \(\mathbf{b}\) and matrix \(\text{diag}(\mathbf{x})\); see [BNT19x]. The multivariate case uses a similar approach, but care is taken to use the right combination of input coordinates.; see [AKL+20].

Notes

• The use of Classical Gram-Schmidt with reorthogonalization was suggested by Yuji Nakatsukasa to improve performance. This uses matrix operations rather than the vector operations of modified Gram-Schmidt allowing the use of more efficient BLAS3 operations.

Parameters:

• X (np.array) – Input coordinates of size (M, dim)
• degree (int or list of ints) – Polynomial degree. If an int, a total degree polynomial is constructed; if a list, the list must be length m and a maximum degree polynomial is constructed.
• weight (None or np.array (M,)) – Initial vector in the Arnoldi iteration
• mode (None or ['total', 'max']) – What type of polynomial basis to construct; only matters if dim>1. If None, the type of basis will be automatically detected.

Returns:

• Q (ndarray) – Orthonormal basis for the desired polynomials on the input coordinates.
• R (ndarray) – Matrix containing orthogonalization information; needed to evaluate polynomial at new coordinates.
• indices (ndarray) – List of indices showing the order in which the basis was constructed.

polyrat.vandermonde_arnoldi_eval(X, R, indices, weight=None)

Evaluate a Vandermonde with Arnoldi polynomial

Parameters:

• X (array_like) – Coordiantes at which to evaluate the polynomial
• R (ndarray) – Upper triangular matrix containing orthogonalization information from vandermonde_arnoldi_CGS.
• indices (ndarray) – The ordering of polynomial powers returned from vandermonde_arnoldi_CGS.
• weight (None or array_like) – Starting vector for Arnoldi; by default is the ones vector. In general, this should not change without good reason.

Returns:

W – Polynomial basis evaluated at the specified points.

Return type:

ndarray

polyrat.vandermonde_arnoldi_eval_der(X, R, indices, weight=None, V=None)

Evaluate the derivative of Vandermonde with Arnoldi polynomial basis

Parameters:

• X (array_like) – Coordiantes at which to evaluate the polynomial
• R (ndarray) – Upper triangular matrix containing orthogonalization information from vandermonde_arnoldi_CGS.
• indices (ndarray) – The ordering of polynomial powers returned from vandermonde_arnoldi_CGS.
• weight (None or array_like) – Starting vector for Arnoldi; by default is the ones vector. In general, this should not change without good reason.
• V (ndarray) – Evaluation of the polynomial without the derivative at X.

Returns:

W – Polynomial basis deriative evaluated at the specified points.

Return type:

ndarray

Lagrange Polynomial Basis

class polyrat.LagrangePolynomialBasis(nodes)

Constructs a Lagrange polynomial basis in barycentric form.

Here we construct a univariate polynomial in barycentric form as described in [BT04]:

\[p(x) = \sum_{j=1}^d \frac{ f_j w_j (x - \xi_j)^{-1}}{w_j (x - \xi_j)^{-1}}, \quad w_j := \prod_{k\ne j} (\xi_j - \xi_k)^{-1}\]

such that \(p(\xi_j) = f_j\).

Parameters:nodes (array_like) – List of the nodes \(\xi_j\) specifying the basis.

polyrat.lagrange_roots(nodes, weights, coef, deflation=True)

Compute the roots of a Lagrange polynomial

This implements the deflation algorithm from [LC14] to compute the roots of a Lagrange polynomial without changing basis and to high accuraccy.

Parameters:

• nodes (numpy.array (n,)) – Nodes \(x_j\) where the polynomial value is defined.
• weights (numpy.array (n,)) – The barycentric weights \(\prod_{k\ne j} (\xi_j - \xi_k)^{-1}\)
• coef (numpy.array (n,)) – The coeffients \(f_j\) defining the value of the polynomial at each point; i.e., \(p(\xi_j) = f_j\).
• deflation (bool) – In the standard formulation to compute these roots we solve a generalized eigenvalue problem (GEP) which has two infinite poles. If True, we explicitly remove these infinite poles by shrinking the matrices; if False we do not.

Returns:

roots – Roots of the polynomial.

Return type:

ndarray