Constructing Orthogonal Polynomials with Minterpy

This is a parent issue for extending Minterpy such that it can be used to compute orthogonal polynomials more conveniently.

About orthogonal polynomials

A sequence of polynomials \{ \pi_n \}_{n=0}^{\infty} each of which having a degree \mathrm{deg}[\pi_n] = n is called orthogonal with respect to a weight function w on the interval (a, b);\; a < b if and only if

\int_a^b \pi_k(x) \, \pi_l(x) \, w(x) \, dx = \lVert \pi_k \rVert^2 \delta_{k, l},

where

\delta_{k, l} = \begin{cases}
0, k \neq l \\
1, k = l
\end{cases}

and

\lVert \pi_k \rVert^2 = \langle \pi_k, \pi_k \rangle_{w} = \int_a^b \pi_k^2(x) \, w(x) \, dx.

The quantity \lVert \pi_k \rVert is called the norm of the function \pi_k.

The polynomials may be normalized:

\tilde{\pi}_k = \frac{\pi_k}{\lVert \pi_k \rVert}

such that:

\lVert \tilde{\pi}_k \rVert^2 = \langle \tilde{\pi}_k, \tilde{\pi}_k \rangle_{w} = \int_a^b \tilde{\pi}_k^2(x) \, w(x) \, dx = 1.

The sequence of orthonormal polynomials \{ \tilde{\pi}_n \}_{n=0}^{\infty} also satisfies:

\int_a^b \tilde{\pi}_k(x) \, \tilde{\pi}_l(x) \, w(x) \, dx = \lVert \tilde{\pi}_k \rVert^2 \delta_{k, l}.

Several known and classic orthogonal polynomials are Legendre, Hermite, Jacobi, and Laguerre each of which has a particular form of weight function. These polynomials are known analytically. However, at least in the context of uncertainty quantification, several common weight functions (correspond to probability density function) are not classics and therefore, the polynomials must be constructed from scratch.

Construction: three-term recurrence relation and Stieltjes procedure

The three-term recurrence relation is a fundamental theorem useful for its construction. The theorem reads as follows for monic orthogonal polynomials:

\pi_{k + 1}(x) = (x - \alpha_k) \pi_k(x) - \beta_k \pi_{k - 1} (x), \; k = 0, \ldots

with \pi_{-1} (x) = 0 and \pi_{0} (x) = 1.

The sequence of pairs \{ (\alpha_k, \beta_k) \}_{k = 0}^\infty is called the recurrence coefficients and it may be obtained using the Stieltjes procedure:

\alpha_k = \alpha_k (w) = \frac{ \langle \pi_k, x \pi_k \rangle_{w}}{\langle \pi_k, \pi_k \rangle_{w}} = \frac{\int_{\mathcal{D}_X} x \, \pi_k^2(x) \, w(x) \, dx }{\int_{\mathcal{D}_X} \pi_k^x(x) \, w(x) \, dx}

and

\beta_k = \beta_k (w) = \frac{ \langle \pi_k, \pi_k \rangle_{w}}{\langle \pi_{k - 1}, \pi_{k - 1} \rangle_{w}} = \frac{\int_{\mathcal{D}_X} \pi_k^2(x) \, w(x) \, dx }{\int_{\mathcal{D}_X} \pi_{k-1}^x(x) \, w(x) \, dx}.

The coefficient \beta_0 is arbitrary but by convention is usually set to \beta_0 = \int_{\mathcal{D}_X} w(x) \, dx.

For orthonormal polynomials, the three-term recurrence relation reads:

\sqrt{\beta_{k + 1}} \tilde{\pi}_{k + 1}(x) = (x - \alpha_k) \tilde{\pi}_k(x) - \sqrt{\beta_k} \tilde{\pi}_{k - 1} (x), \; k = 0, \ldots

with \pi_{-1} (x) = 0 and \pi_{0} (x) = \frac{1}{\sqrt{\beta_0}}.

Except for taking the square root of \beta_k's, the recurrence coefficients are the same.

In other words, the computation of the recurrence coefficients requires the computation of inner products of polynomials with respect to a given weight function w.

Open issues

I would like to carry out the Stieltjes procedure with Minterpy and it requires a couple of extensions to the current codebase so that we can manipulate Minterpy polynomials more conveniently. The core computation, computing inner products, is an integral operation; we have implemented polynomial integration in Minterpy for all supported bases.

Here are several proposed extensions; note that the construction of orthogonal polynomials as described above is only in one dimension:

Scalar multiplication of polynomials (currently, only polynomials can be multiplied)
Scalar division of polynomials (i.e., dividing a polynomial with a scalar); in essence this is a scalar multiplication with the scalar reciprocal.
Scalar addition and subtraction

Some of the operations above are possible for the Lagrange basis if a constant polynomial is first constructed from the scalar; but this is not convenient. Moreover, the operations are not supported yet for all polynomial bases so there are still inconsistencies in the higher-level interface.

Furthermore, while polynomial integration via interpolatory quadrature is possible in Minterpy, some common weight functions are not polynomials (e.g., Hermite, Jacobi with non-integer parameters, Laguerre weight functions). A naive approach would be to create an interpolant polynomial for the weight function first and use the polynomial in the Stieltjes procedure. However, this may not be optimal because we don't know in advance the degree requires to interpolate a given weight function. Moreover, weight functions may even have integrable singularities that cannot be interpolated nicely.

A more sensible approach is perhaps to use another interpolatory quadrature rule; Gautschi recommends the Fejer-1 quadrature rule in [-1, 1] (the natural domain of Minterpy) and to deal with other bounds using fractional transformation.

Investigate the extension of the method integrate_over() to allow passing a weight function, either a polynomial or a generic callable. When it's generic callable, use Fejer-1 rule to compute the integral (what if it's multiple dimension?).

The goal of these proposed extensions is not a direct computation of orthogonal polynomials and their recurrence relations by Minterpy per se, but for Minterpy to better support such computations via the Stieltjes procedure. In other words, the extensions are made such that Minterpy can handle computations involving polynomials more naturally and consistently (scalar addition, subtraction, multiplication, weighted integration, etc.).