Case Study 1 — Orthogonal Polynomials: Geometry That Powers Numerical Computing
Field: numerical analysis / signal processing. This case study cashes in §34.4: the function inner product, plus Gram–Schmidt from Chapter 20, manufactures the orthogonal polynomial families that sit inside numpy, scipy, and essentially every scientific-computing library on Earth.
The problem: approximating and integrating with polynomials
Polynomials are the workhorses of computation — fast to evaluate, easy to differentiate and integrate. So two everyday tasks of numerical computing are: approximate a complicated function by a polynomial, and integrate a function you can only evaluate at a few points. The naive approach to both uses the monomial basis $1, x, x^2, x^3, \dots$ — and the naive approach is a numerical disaster. The reason, and the cure, are pure inner-product-space geometry.
Here is the trouble with monomials. On the interval $[-1,1]$ with the function inner product $\langle f,g\rangle=\int_{-1}^1 fg\,dx$, the monomials are far from orthogonal. We computed in §34.4 that the angle between $1$ and $x^2$ is about $42^\circ$ — nowhere near the $90^\circ$ of perpendicularity. Worse, as you go to higher degree, $x^4$, $x^6$, $x^8$ all look almost identical on $[-1,1]$: they are nearly flat near the origin and shoot up near the endpoints, so they point in nearly the same direction in function space. A basis of nearly-parallel vectors is the numerical equivalent of a wobbly, near-degenerate coordinate system. Expressing a function in it requires solving a linear system whose matrix (the famous Hilbert matrix, whose $(i,j)$ entry is $\langle x^{i},x^{j}\rangle$) is catastrophically ill-conditioned — the condition number explodes exponentially with the degree, a phenomenon Chapter 38 will diagnose precisely. Tiny rounding errors in the data become enormous errors in the coefficients.
The cure is the one this chapter is built on: replace the skewed monomial basis with an orthogonal one. Orthogonal directions are independent — perpendicular axes do not contaminate one another — so coordinates become trivial to compute (one inner product each, §34.7) and numerically stable. And we already own the machine that builds orthogonal bases from arbitrary ones: Gram–Schmidt, run in the function inner product space.
Building the Legendre polynomials
Run Gram–Schmidt on $1, x, x^2, x^3, \dots$ under $\langle f,g\rangle=\int_{-1}^1 fg\,dx$ and out come the Legendre polynomials (up to scaling), exactly as we did for the first three in §34.4. The recipe is identical to Chapter 20's — subtract from each new monomial its shadow on everything already straightened — only the inner product changed:
$$ P_0 = 1,\quad P_1 = x,\quad P_2 \propto x^2-\tfrac13,\quad P_3 \propto x^3-\tfrac35 x,\quad\dots $$
These are mutually orthogonal: $\int_{-1}^1 P_j P_k\,dx = 0$ whenever $j\neq k$. Let us verify the orthogonality numerically, treating each polynomial as a vector of samples (the Chapter 5 picture, "a function is a vector with very many components"):
# Legendre polynomials are Gram-Schmidt on monomials under <f,g> = integral_{-1}^{1} f g dx.
import numpy as np
x = np.linspace(-1, 1, 4000); dx = x[1] - x[0]
def ip(f, g): return float(np.sum(f * g) * dx)
def gram_schmidt(vs, ip):
Q = []
for v in vs:
u = np.array(v, float)
for q in Q:
u = u - (ip(v, q) / ip(q, q)) * q
Q.append(u)
return Q
mono = [x**k for k in range(4)] # 1, x, x^2, x^3
P = gram_schmidt(mono, ip)
G = np.array([[ip(P[i], P[j]) for j in range(4)] for i in range(4)])
print(np.round(G, 4))
# [[ 2. 0. -0. 0. ] <- ||1||^2 = 2
# [ 0. 0.6667 -0. 0. ] <- ||x||^2 = 2/3
# [-0. -0. 0.1778 -0. ] <- ||x^2 - 1/3||^2 = 8/45
# [ 0. 0. -0. 0.0457]] <- ||x^3 - (3/5)x||^2 = 8/175
print(round(P[2][-1], 4), round(P[3][-1], 4)) # 0.6667 0.4 -> values at x=1: x^2-1/3, x^3-(3/5)x
The Gram matrix $G_{ij}=\langle P_i,P_j\rangle$ comes out diagonal — every off-diagonal entry is $0$ to four decimals — which is the numerical signature of an orthogonal set. The diagonal entries are the squared norms $\lVert P_k\rVert^2$. The last line confirms the monic forms: $P_2(1)=1-\tfrac13=\tfrac23$ and $P_3(1)=1-\tfrac35=\tfrac25$, matching the hand formulas. We rebuilt a cornerstone of numerical analysis with the Chapter 20 algorithm and a one-line inner product.
Why orthogonality buys numerical stability
The payoff is that, in the orthogonal basis, the coordinates of any function decouple. To approximate a function $f$ by a degree-$n$ polynomial, you project onto each $P_k$ independently: the coefficient is $c_k = \langle P_k,f\rangle/\langle P_k,P_k\rangle$ (the generalized Fourier coefficient of §34.7, scaled because the $P_k$ are orthogonal but not orthonormal). There is no ill-conditioned linear system to solve — each coefficient is a single, stable inner product. Adding a higher-degree term never disturbs the coefficients already computed, exactly the independence that made Fourier coefficients clean in Chapter 22. Orthogonality converts a near-degenerate problem into $n$ separate, well-conditioned ones.
This is also the principle behind Gaussian quadrature, the gold standard for numerical integration. The roots of the degree-$n$ Legendre polynomial, used as evaluation points with suitable weights, integrate every polynomial of degree up to $2n-1$ exactly — a near-miraculous efficiency that falls directly out of the orthogonality of the $P_k$. When scipy.integrate.quad integrates your function to machine precision in a handful of evaluations, it is exploiting the geometry of this chapter.
Chebyshev polynomials: orthogonality in a weighted inner product
Legendre is the right family for the plain integral, but a different application wants a different geometry — and this is where the weighted inner product of §34.3 meets the function inner product. For polynomial approximation that minimizes the worst-case error (rather than the average squared error), the right basis is the Chebyshev polynomials $T_0=1$, $T_1=x$, $T_2=2x^2-1$, $T_3=4x^3-3x$, and they are orthogonal not under the plain integral but under the weighted function inner product
$$ \langle f,g\rangle_w = \int_{-1}^1 \frac{f(x)\,g(x)}{\sqrt{1-x^2}}\,dx. $$
The weight $1/\sqrt{1-x^2}$ piles up near the endpoints $x=\pm1$, which is exactly why Chebyshev approximation controls the endpoint error that ruins naive polynomial interpolation (the Runge phenomenon). This is a single inner product space combining both generalizations from the chapter: functions (so the "vectors" are polynomials) and a weight (so the geometry is reshaped). The weight is positive on $(-1,1)$, so $\langle\cdot,\cdot\rangle_w$ is a genuine inner product and all of Part IV applies.
Verifying Chebyshev orthogonality numerically takes one clean trick: the substitution $x=\cos\theta$ turns the awkward endpoint-singular integral into an ordinary one, $\int_{-1}^1 \frac{f g}{\sqrt{1-x^2}}dx = \int_0^\pi f(\cos\theta)\,g(\cos\theta)\,d\theta$. (This substitution is itself the secret identity of the Chebyshev polynomials — $T_k(\cos\theta)=\cos(k\theta)$, so they are cosines in disguise, which is why their orthogonality mirrors the Fourier cosine orthogonality of Chapter 22.)
# Chebyshev T_k are orthogonal under the weighted IP <f,g> = integral f g / sqrt(1-x^2) dx.
# Substitution x = cos(theta) removes the endpoint singularity: integral_0^pi f(cos t) g(cos t) dt.
import numpy as np
t = np.linspace(0, np.pi, 200_000)
xt = np.cos(t)
def wip(f, g): return float(np.trapezoid(f * g, t)) # weighted IP via the cos-theta substitution
T = [np.ones_like(xt), xt, 2*xt**2 - 1, 4*xt**3 - 3*xt] # T0..T3
G = np.array([[wip(T[i], T[j]) for j in range(4)] for i in range(4)])
print(np.round(G, 4))
# [[3.1416 0. 0. 0. ] <- <T0,T0> = pi
# [0. 1.5708 0. 0. ] <- <T1,T1> = pi/2
# [0. 0. 1.5708 0. ] <- <T2,T2> = pi/2
# [0. 0. 0. 1.5708]] <- <T3,T3> = pi/2
Again the Gram matrix is diagonal — the Chebyshev polynomials are orthogonal under the weighted inner product — with $\langle T_0,T_0\rangle_w=\pi$ and $\langle T_k,T_k\rangle_w=\pi/2$ for $k\ge 1$, the known normalization constants. The off-diagonal zeros confirm mutual orthogonality. Switching the weight changed which polynomials count as perpendicular, and a different famous family emerged from the same Gram–Schmidt idea, run in a different geometry.
What this case study shows
Two of the most important tools in scientific computing — Legendre polynomials for stable approximation and exact quadrature, Chebyshev polynomials for minimax approximation — are not special inventions. They are Gram–Schmidt (Chapter 20) applied in a function inner product space (§34.4), with Legendre using the plain integral and Chebyshev using a weighted one (§34.3). The orthogonality that makes them numerically stable is the same orthogonality that decoupled Fourier coefficients in Chapter 22 and that will let us read off qubit amplitudes in Case Study 2. The lesson is the chapter's thesis, sharpened by an application: choose your inner product to fit the problem, and the geometry of Part IV — orthogonal bases, independent coordinates, stable projection — does the rest. When a library hands you numpy.polynomial.legendre or a scipy quadrature rule, it is handing you the geometry of inner product spaces, compiled.