Polynomial expansions#

Consider a function \(f(x)\) expressed as a finite polynomial expansion

\[f(x) = \sum_{k=0}^d c_k p_k(x)\]

Smooth functions can often be accurately approximated by truncated expansions over a complete polynomial basis \(\{P_k(x)\}\). SeeMPS provides tools for both direct encoding of functions and function composition using polynomial expansions.

SeeMPS provides two distinct but related sets of polynomial approximation tools:

Exact construction of MPS for 1D polynomials, when the coefficients are known explicitly.
Polynomial expansions of operators or states, evaluating \(f(A)\) where \(A\) is either an MPS or an MPO.

Exact polynomial MPS constructions#

The function mps_from_polynomial() constructs the exact MPS corresponding to a one-dimensional polynomial expressed in the monomial basis \(p_k(x)=x^k\) over some equispaced discretization \([a, b]\), represented by the class RegularInterval.

Any degree-\(d\) polynomial can be constructed using an MPS with bond dimension \(\chi \leq d+1\) (in practice, often smaller if compression routines are used). This function can either take the coefficients \(\{p_i\}\) or convert a NumPy representation of a polynomial to MPS form.

Polynomial expansions of MPS and MPOs#

This implements the remaining functionality using the to_mps() and to_mpo() methods. These apply the Clenshaw evaluation method, a numerically stable technique to evaluate polynomials in situations of finite precision.

Function composition#

Given an MPS \(\mathbf{v}\) encoding a function \(g(x)\), we can compute the encoding \(\mathbf{w}\) of the composed function through the expansion:

\[(f \circ g)(x) = f(g(x)) \to \mathbf{w} = \sum_k c_k P_k(\mathbf{v})\]

The same technique applies straightforwardly to operator-valued functions, allowing the evaluation of \(f(O)\) for an operator \(O\) represented as an MPO.

Available expansion classes#

SeeMPS currently provides the following expansion classes:

PowerExpansion: A polynomial expansion in the monomial basis \(p_k(x)=x^k\). In this scenario, Clenshaw’s evaluation formula reduces to Horner’s method, an efficient and robust technique for evaluating polynomial expansions. The user must explicitly provide the coefficients \(\{c_k\}\).
ChebyshevExpansion: An expansion in the orthogonal basis of Chebyshev polynomials. The use of Chebyshev polynomials yields an approximation framework analogous to Matlab’s “ChebFun” package, but formulated entirely within the MPS/MPO formalism.
LegendreExpansion: An expansion in the orthogonal basis of Legendre polynomials.

Three-term recurrence#

Orthogonal polynomial expansions rely on three-term recurrence relations:

\[P_{k+1}(x) = (\alpha_k x + \beta_k) P_k(x) - \gamma_k P_{k-1}(x)\]

which are evaluated using numerically stable Clenshaw formulas. To completely determine the basis, two additional coefficients determining the linear term and fixing affine translations are required:

\[P_{1}(x) = \sigma x + \mu\]

The user provides the target function \(f\) together with an initial MPS or MPO encoding the argument.

This expansion framework is easily extensible to any classical orthogonal polynomial family by subclassing PolynomialExpansion, requiring only the three-term recurrence relation, the affine fixing coefficients, and the orthogonality domain of the basis.

Coefficient computation#

All expansion objects can be constructed by providing the coefficients \([c_0,c_1,...]\) explicitly. Alternatively, the coefficients can be computed by projecting the target function onto the orthogonal basis through projection methods, which estimate a finite-order expansion of a scalar function using numerical quadratures.

Limitations#

The applicability of this technique is limited by the regularity of the target function and the properties of the basis. Generally, highly differentiable functions present favorable convergence rates, while functions with discontinuities or sharp features—such as Heaviside functions—are poorly approximated by polynomial expansions, requiring prohibitively large expansion orders.

Example#

An example demonstrating the use of these functions for the case of Chebyshev polynomials is shown in Chebyshev.ipynb.

The following example encodes a multivariate function \(f(x, y) = e^{x + y}\) using a Chebyshev expansion:

>>> import numpy as np
>>> from seemps.state import mps_tensor_sum
>>> from seemps.analysis.mesh import RegularInterval
>>> from seemps.analysis.factories import mps_interval
>>> from seemps.analysis.expansion import ChebyshevExpansion
>>>
>>> interval = RegularInterval(-1, 1, 10)
>>> mps_x = mps_interval(interval)
>>> mps_xy = mps_tensor_sum([mps_x] * 2)
>>> f = lambda x: np.exp(x)
>>> expansion = ChebyshevExpansion.project(f, (-1, 1))
>>> mps_f = expansion.to_mps(argument=mps_xy)

`PolynomialExpansion`	Abstract base class for polynomial expansions.
`PowerExpansion`	Polynomial expansion in the monomial basis {1, x, x^2, ...}, coresponding to a standard power or Taylor series.
`ChebyshevExpansion`	Expansion in the Chebyshev basis.
`LegendreExpansion`	Expansion in the Legendre basis.
`project`(func[, approximation_domain, order])	Project a scalar function onto the Chebyshev basis on the given approximation domain.
`mps_from_polynomial`(p, domain[, first, strategy])	Construct a tensor representation of a polynomial.