from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

9.4.Orthogonal polynomials¶

Interpolation is not the only way to use polynomials for global approximation of functions. In Fitting functions to data we saw how to find least-squares polynomial fits to data by solving linear least-squares matrix problems. This idea can be extended to fitting functions.

Example 9.4.1 (Continuous least-squares fitting)

Let’s approximate $e^x$ over the interval $[−1,1]$ . We can sample it at, say, 15 points, and find the best-fitting straight line to that data.

from numpy.linalg import lstsq
t = linspace(-1, 1, 15)
y = exp(t)
plot(t, y, label="function")

V = [[ti**j for j in range(2)] for ti in t]
c = lstsq(V, y, rcond=None)[0]
print("fit coeffs:", c)

x = linspace(-1, 1, 600)
plot(x, c[1] + c[0] * x, label="fit")
xlabel("x"),  ylabel("value")
legend(),  title("Least squares fit of exp(x)");

fit coeffs: [1.20159125 1.11828384]

There’s nothing special about 15 points. Choosing more doesn’t change the result much.

t = linspace(-1, 1, 150)
y = exp(t)
plot(t, y, label="function")

V = [[ti**j for j in range(2)] for ti in t]
c = lstsq(V, y, rcond=None)[0]
print("fit coeffs:", c)

x = linspace(-1, 1, 600)
plot(x, c[1] + c[0] * x, label="fit")
xlabel("x"),  ylabel("value")
legend(),  title("Least squares fit of exp(x)");

fit coeffs: [1.17767125 1.10507318]

This situation is unlike interpolation, where the degree of the interpolant increases with the number of nodes. Here, the linear fit is apparently approaching a limit that we may think of as a continuous least-squares fit.

n = arange(40, 420, 60)
results = PrettyTable(["n", "intercept", "slope"])
slope = zeros(n.size)
intercept = zeros(n.size)

for k in range(n.size):
    t = linspace(-1, 1, n[k])
    y = exp(t)
    V = [[ti**j for j in range(2)] for ti in t]
    c = lstsq(V, y, rcond=None)[0]
    results.add_row([n[k], c[1], c[0]])

print(results)

+-----+--------------------+--------------------+
|  n  |     intercept      |       slope        |
+-----+--------------------+--------------------+
|  40 | 1.1090551764670034 | 1.1846492800335244 |
| 100 | 1.1057932881809656 | 1.1789195568501278 |
| 160 | 1.1049832955976733 | 1.1775158384490083 |
| 220 | 1.104615859586127  | 1.176881503540092  |
| 280 | 1.104406105963981  | 1.176520063299236  |
| 340 | 1.1042704646318244 | 1.1762865906944793 |
| 400 | 1.104175553929563  | 1.1761233467177028 |
+-----+--------------------+--------------------+

We can extend least-squares fitting from data to functions by extending several familiar finite-dimensional definitions. The continuous extension of a sum is an integral, which leads to the following.

9.4.1Quasimatrices¶

If we are extending our notion of vectors to include continuous functions, what should serve as an extension of a matrix? One of our most important interpretations of a matrix is the connection to linear combinations. For instance, the Vandermonde-like system $\mathbf{V} \mathbf{c} \approx \mathbf{y}$ from (3.1.2) is the statement

\mathbf{y} \approx \mathbf{V} \mathbf{c} = c_0 \mathbf{v}_0 + c_1 \mathbf{v}_1 + \cdots + c_n \mathbf{v}_n,

(9.4.4)

in which the $i$ th row of $\mathbf{v}_j$ is $t_i^j$ . This was derived as a discrete approximation for $j=0,\ldots,n$ .

\mathbf{y} \approx c_0 + c_1 x + \cdots + c_n x^n,

(9.4.5)

which we want to abbreviate in “matrix”-vector form.

Definition 9.4.2 (Quasimatrix)

Given functions $f_1,\ldots,f_n$ in inner product space $S$ , define the quasimatrix

\mathbf{F} = \begin{bmatrix} \underline{f_1(x)} & \underline{f_2(x)} & \cdots & \underline{f_n(x)} \end{bmatrix}.

(9.4.6)

For a vector $\mathbf{z} \in \real^n$ , define the quasimatrix-vector product

\mathbf{F}\mathbf{z} = z_1f_1(x) + z_2f_2(x) + \cdots + z_n f_n(x).

(9.4.7)

For another function $g\in S$ , define the adjoint product

\mathbf{F}^T g = \begin{bmatrix} \langle f_1,g \rangle \\ \vdots \\ \langle f_n,g \rangle \end{bmatrix}.

(9.4.8)

We consider any other expressions involving a quasimatrix to be undefined. It might help to think of $\mathbf{F}$ as an $\infty\times n$ matrix, which is consistent with the definitions that $\mathbf{F}\mathbf{z}$ is a function ( $\infty\times 1$ ), $\mathbf{F}^T g$ is a vector ( $n\times 1$ ), and $\mathbf{F}^T\mathbf{F}$ is a matrix ( $n \times n$ ). When infinite dimensions combine in a product, we use integrals rather than sums.

Example 9.4.2

Let $\mathbf{F} = \bigl[ \,\underline{\cos(\pi x)} \quad \underline{\sin(\pi x)}\, \bigr]$ . Then

\begin{split} \mathbf{F} \begin{bmatrix} -2 \\ 1 \end{bmatrix} &= -2\cos(\pi x) + \sin(\pi x), \\ \mathbf{F}^T x & = \begin{bmatrix} \int_{-1}^1 x\cos(\pi x)\, dx \\ \int_{-1}^1 x\sin(\pi x)\, dx\end{bmatrix} = \begin{bmatrix} 0 \\ 2/\pi \end{bmatrix}, \\ \mathbf{F}^T\mathbf{F} &= \begin{bmatrix} \int_{-1}^1 \cos^2(\pi x)\, dx & \int_{-1}^1 \cos(\pi x)\sin(\pi x)\, dx \\ \int_{-1}^1 \cos(\pi x)\sin(\pi x)\, dx & \int_{-1}^1 \sin^2(\pi x)\, dx \end{bmatrix} = \mathbf{I}. \end{split}

(9.4.10)

9.4.2Normal equations¶

The discrete linear least-squares problem of minimizing $\| \mathbf{y} - \mathbf{V} \mathbf{c} \|_2$ over all possible $\mathbf{c}$ , given matrix $\mathbf{V}$ and data vector $\mathbf{y}$ , has a solution via the normal equations (3.2.3),

\mathbf{c} = \left(\mathbf{V}^T\mathbf{V}\right)^{-1} \mathbf{V}^T \mathbf{y}.

(9.4.11)

We can now reinterpret (9.4.11) in terms of quasimatrices.

There is no need to supply a proof of Theorem 9.4.1 because it will read exactly the same as for the discrete normal equations. All the effort has gone into making definitions that set up a perfect analogy. In retrospect, all we needed in the original discrete case were linear combinations and inner products.

Example 9.4.3

We revisit approximation of $e^x$ as suggested in Example 9.4.1. With the Vandermonde quasimatrix $\mathbf{V}= \begin{bmatrix} \underline{1} & \underline{x} \end{bmatrix}$ , we get

\mathbf{V}^Te^x = \begin{bmatrix} \langle 1,e^x \rangle \\[1mm] \langle x,e^x \rangle \end{bmatrix} = \begin{bmatrix} \int_{-1}^1 e^x\, dx \\[1mm] \int_{-1}^1 x e^x\, dx \end{bmatrix} = \begin{bmatrix} e-e^{-1} \\ 2 e^{-1} \end{bmatrix}

(9.4.14)

and

\mathbf{V}^T \mathbf{V} = \begin{bmatrix} \langle 1,1 \rangle & \langle 1,x \rangle \\[1mm] \langle x,1 \rangle & \langle x,x \rangle \end{bmatrix} = \begin{bmatrix} 2 & 0 \\ 0 & 2/3 \end{bmatrix}.

(9.4.15)

The normal equations (9.4.13) therefore have solution

\mathbf{c} = \begin{bmatrix} 2 & 0 \\ 0 & 2/3 \end{bmatrix}^{-1} \begin{bmatrix} e-e^{-1} \\ 2 e^{-1} \end{bmatrix} = \begin{bmatrix} \sinh(1) \\ 3e^{-1} \end{bmatrix} \approx \begin{bmatrix} 1.175201\\ 1.103638 \end{bmatrix},

(9.4.16)

which is well in line with the values found in Example 9.4.1.

If we extend $\mathbf{V}$ by an additional column for $x^2$ , then we need to calculate $\int_{-1}^1 x^2 e^x \, dx = e - 5e^{-1}$ and $\int_{-1}^1 x^4\, dx = 2/5$ to get

\mathbf{c} = \begin{bmatrix} 2 & 0 & 2/3\\ 0 & 2/3 & 0 \\ 2/3 & 0 & 2/5 \end{bmatrix}^{-1} \begin{bmatrix} e-e^{-1} \\ 2 e^{-1} \\ e - 5e^{-1} \end{bmatrix} \approx \begin{bmatrix} 0.9962940 \\ 1.103638 \\ 0.5367215 \end{bmatrix}.

(9.4.17)

9.4.3Legendre polynomials¶

Equation (9.4.13) becomes much simpler if $\mathbf{V}^T\mathbf{V}$ is diagonal. By our definitions, this would imply that the columns of $\mathbf{V}$ are mutually orthogonal in the sense of the function inner product. This is not the case for the monomial functions $x^j$ . But there are orthogonal polynomials which do satisfy this property.

For what follows, let $\mathcal{P}_n \subset S$ be the set of polynomials of degree $n$ or less.

Here are some key facts that are straightforward to prove.

Now let us define the quasimatrix

\mathbf{L}_n(x) = \begin{bmatrix} \alpha_0^{-1} \underline{P_0} & \alpha_1^{-1} \underline{P_1} & \cdots & \alpha_{n}^{-1} \underline{P_{n}} \end{bmatrix}.

(9.4.20)

Then $\mathbf{L}_n^T\mathbf{L}_n=\mathbf{I}$ . The normal equations (9.4.13) thus simplify accordingly. Unraveling the definitions, we find the least-squares solution

\mathbf{L}_n \bigl( \mathbf{L}_n^T f \bigr) = \sum_{k=0}^n c_k P_k(x), \quad \text{where } c_k = \frac{1}{\alpha_k^2} \langle P_k,f \rangle.

(9.4.21)

9.4.4Chebyshev polynomials¶

Equation (9.4.1) is not the only useful way to define an inner product on a function space. It can be generalized to

\langle f,g \rangle = \int_{-1}^1 f(x)g(x)w(x)\,dx

(9.4.22)

for a positive function $w(x)$ called the weight function of the inner product. An important special case is

\langle f,g \rangle = \int_{-1}^1 \frac{f(x)g(x)}{\sqrt{1-x^2}}\,dx.

(9.4.23)

Chebyshev polynomials also have a startling alternative form,

T_k(x) = \cos\left( k \theta \right), \quad \theta = \arccos(x).

(9.4.25)

The results from Theorem 9.4.2 apply to Chebyshev polynomials as well, with orthogonality being in the sense of (9.4.23). Their Gram matrix is given by

\langle T_i,T_j \rangle = \begin{cases} 0, & i\neq j, \\ \gamma_0^2 = \pi, & i=j=0, \\ \gamma_i^2=\pi/2, & i=j>0. \end{cases}

(9.4.26)

The least-squares solution is not the same in the Legendre and Chebyshev cases: both find the nearest approximation to a given $f(x)$ , but the norm used to measure distances is not the same.

9.4.5Roots of orthogonal polynomials¶

Interesting properties can be deduced entirely from the orthogonality conditions. The following result will be relevant in Spectrally accurate integration. The same result holds for orthogonal polynomial families with different weight functions, such as the Chebyshev polynomials.

Proof

Let $x_1,\ldots,x_m$ be all the distinct roots of $P_n(x)$ between -1 and 1 at which $P_n(x)$ changes sign (in other words, all roots of odd multiplicity). Define

r(x) = \prod_{i=1}^m (x-x_i).

(9.4.27)

By definition, $r(x)P_n(x)$ does not change sign over $(-1,1)$ . Therefore

\int_{-1}^1 r(x)P_n(x) \, dx \neq 0.

(9.4.28)

Because $r$ is a degree- $m$ polynomial, we can express it as a combination of $P_0,\ldots,P_m$ . If $m<n$ , the integral (9.4.28) would be zero, by the orthogonality property of Legendre polynomials. So $m\ge n$ . Since $P_n(x)$ has at most $n$ real roots, $m=n$ . All the roots must therefore be simple, and this completes the proof.

The result of Theorem 9.4.3 holds for orthogonal families of polynomials for other weight functions. The Chebyshev case is unusual in that thanks to (9.4.25), the roots of $T_n$ are known explicitly:

t_k = \cos\left(\frac{2k-1}{2n}\pi\right), \qquad k=1,\ldots,n.

(9.4.29)

These are known as the Chebyshev points of the first kind. The chief difference between first-kind and second-kind points is that the latter type include the endpoints $\pm 1$ . Both work well for polynomial interpolation and give spectral convergence.

9.4.6Least squares versus interpolation¶

Both interpolation and the solution of a linear least-squares problem produce a projection of a function into the space of polynomials $\mathcal{P}_n$ . In the least-squares case, the close connection with inner products and orthogonality makes the 2-norm, perhaps with a weight function, a natural setting for analysis. Because a constant weight function is the simplest choice, Legendre polynomials are commonly used for least squares.

Interpolation has no easy connection to inner products or the 2-norm. With interpolants a different kind of approximation analysis is more fruitful, often involving the complex plane, in which the max-norm is the natural choice. For reasons beyond the scope of this text, Chebyshev polynomials are typically the most convenient to work with in this context.