Dot products of functions¶

We can apply the Gram–Schmidt process to any vector space as long as we define a dot product (also called an inner product). (Technically, a continuous ("complete") vector space equipped with an inner product is called a Hilbert space.)

For column vectors, the usual dot product is to multiply the components and add them up.

But (real-valued) functions $f(x)$ also define a vector space (you can add, subtract, and multiply by constants). In particular, consider functions defined on the interval $x \in [-1,1]$. The "components" of $f$ can be viewed as its values $f(x)$ at each point in the domain, and the obvious analogue of "summing the components" is the integral. Hence, the most obvious "dot product" of two functions in this space is:

$$ f \cdot g = \int_{-1}^1 f(x) g(x) \, dx $$

Such a generalized inner product is commonly denoted $\langle f, g \rangle$ (or $\langle f | g \rangle$ in physics).

This satisfies the key properties of dot products that make linear algebra "work":

• symmetry: $f \cdot g = g \cdot f$
• linearity: $f \cdot (\alpha g + \beta h) = \alpha (f\cdot g) + \beta (f \cdot h)$
• positivity: $f \cdot f = \Vert f \Vert^2 \ge 0$, and $=0$ only if $f = 0$ (almost everywhere).

As long as the dot product has these properties, much of what we do in 18.06 will "just work" for functions too. Gram–Schmidt will just work.

Orthogonal polynomials¶

In particular, let us consider a subspace of functions defined on $[-1,1]$: polynomials $p(x)$ (of any degree).

One possible basis of polynomials is simply:

$$ {1, x, x^2, x^3, \ldots} $$

(There are infinitely many polynomials in this basis because this vector space is infinite-dimensional.)

Instead, let us apply Gram–Schmidt to this basis in order to get an orthogonal basis of polynomials known as the Legendre polynomials.

Julia code¶

I'll use the Polynomials package to do polynomial arithmetic for me.

However, I'll need to define a few extra methods to perform my dot products from above, and I also want to display ("pretty print") the polynomials a bit more nicely than the default.

Gram–Schmidt on polynomials¶

Now, let's apply Gram–Schmidt on the polynomials $a_i = x^i$ for $i = 0, 1, \ldots$.

Ordinarily, in Gram–Schmidt, I would normalize each result $p(x)$ by dividing by $\Vert p \Vert = \sqrt{p \cdot p}$, but that will result in a lot of annoying square roots. Instead, I will divide by $p(1)$ to result in the more conventional Legendre polynomials. This is not an "orthonormal" basis because we have chosen a different normalization, but it is still an orthogonal basis.

That means that to get $p_i(x)$, I will do:

$$ v_i(x) = a_i(x) - \sum_{j=0}^{i-1} p_j(x) \frac{p_j \cdot a_i}{p_j \cdot p_j} $$

We would get an ordinary orthnormal basis of polynomials with

$$ q_i(x) = v_i(x) / \Vert v_i \Vert $$

but instead we will use the Legendre normalization

$$ p_i(x) = v_i(x) / v_i(1) $$

where I explicitly divide by $p_i \cdot p_i$ in the projections to compensate for the lack of normalization.

In Julia, I will use the special syntax 2 // 3 to construct the exact rational $\frac{2}{3}$, etc. This will allow me to see the exact Legendre polynomials without any roundoff errors or annoying decimals.

Orthogonalization didn't change $x$, because $x$ and $1$ are already orthogonal under this dot product. In fact, any even power of $x$ is orthogonal to any odd power (because the dot product is the integral of an even function times an odd function).

On the other hand, $x^2$ and $1$ are not orthogonal, so orthogonalizing them leads to a different polynomial of degree 2:

It quickly gets tiresome to type in these expressions one by one, so let's just write a function to compute the Legendre polynomials $p_0, \ldots, p_n$:

Let's display them more nicely with LaTeX:

Key things to notice:

• The polynomials contain only even or only odd powers of $x$, but not both. The reason is that the even and odd powers of $x$ are already orthogonal under this dot product, as noted above.

• A key property of Gram–Schmidt is that the first k vectors span the same space as the original first k vectors, for any k. In this case, it means that $p_0, \ldots, p_k$ span the same space as $1, x, \ldots, x^k$. That is, the $p_0, \ldots, p_k$ polynomials are an orthogonal basis for all polynomials of degree k or less.

These polynomials are very special in many ways. To get a hint of that, let's plot them:

Note that $p_n(x)$ has exactly $n$ roots in the interval $[-1,1]$!

This is essentially required by the fact that they are orthogonal: $p_n$ has to oscillate in sign faster and faster in $[-1,1]$ as $n$ increases in order to integrate to zero against the previous polynomials.

Expanding a polynomial in the Legendre basis.¶

Now that we have an orthogonal (but not orthonormal) basis, it is easy to take an arbitrary polynomial $p(x)$ and write it in this basis:

$$ p(x) = \alpha_0 p_0(x) + \alpha_1 p_1(x) + \cdots = \sum_{i=0}^\infty \alpha_i p_i(x) $$

because we can get the coefficients $\alpha_i$ merely by projecting:

$$ \alpha_i = \frac{p_i \cdot p}{p_i \cdot p_i} $$

Note, however, that this isn't actually an infinite series: if the polynomial $p(x)$ has degree $d$, then $\alpha_i = 0$ for $i > d$. The polynomials $p_0, \ldots, p_d$ are a basis for the subspace of polynomials of degree $d$ (= span of $1, x, \ldots, x^d$)!

Let's see how this works for a "randomly" chosen $p(x)$ of degree 5:

Here are the coefficients α:

Let's check that the sum of $\alpha_i p_i(x)$ gives $p(x)$:

Polynomial fits¶

Review: Projections and Least-Square¶

Given a matrix $Q$ with $n$ orthonormal columns $q_i$, we know that the orthogonal projection

$$ p = QQ^T b = \sum_{i=1}^n q_i q_i^T b $$

is the closest vector in $C(Q)$ to $b$. That is, it minimizes the distance:

$$ \min_{p \in C(Q)} \Vert p - b \Vert \; . $$

Closest polynomials¶

Now, suppose that we have some function $f(x)$ on $x\in[-1,1]$ that is not a polynomial, and we want to find the closest polynomial of degree $n$ to $f(x)$ in the least-square sense. That is, we want to find the polynomials $p(x)$ of degree $n$ that minimizes

$$ \min_{p\in \mathcal{P}_n} \int_{-1}^1 |f(x)-p(x)|^2 dx = \min_{p\in \mathcal{P}_n} \Vert f(x) - p(x) \Vert^2 $$

where

$$ \mathcal{P}_n = \operatorname{span} \{ 1, x, x^2, \ldots, x^n \} = \operatorname{span} \{ p_0(x), p_1(x), \ldots, p_n(x) \} $$

is the space of polynomials of degree $\le n$, spanned by our Legendre polynomials up to degree $n$.

Presented in this context, we can see that this is the same problem as our least-square problem above, and the solution should be the same: $p(x)$ is the orthogonal projection of $f(x)$ onto $\mathcal{P}_n$, given by:

$$ p(x) = p_0(x) \frac{p_0 \cdot f}{p_0 \cdot p_0} + \cdots p_n(x) \frac{p_n \cdot f}{p_n \cdot p_n} \; . $$

Let's try this out for $f(x) = e^x$. Because we're lazy, we'll have Julia compute the integrals numerically using its quadgk function, and fit it to polynomials of degree 5 using our Legendre polynomials from above.

Now, let's use dot products to compute the coefficients in the $p_i(x)$ expansion above for $f(x) = e^x$ (the exp function in Julia):

One thing to notice is an important fact: expanding functions, especially smooth functions, in orthogonal bases like Legendre polynomials or Fourier series tends to converge very rapidly.

Let's write out the resulting polynomial $p(x)$:

Let's plot it:

They are so close that you can hardly tell the difference!

Let's plot the fits for degree 0, 1, …, 3 so that we can watch it converge:

By degree 3, it is hard to tell the difference from $e^x$.

Fitting a non-smooth function¶

It may not be so surprising that we can fit $e^x$ to polynomials; after all, $e^x$ has a convergent Taylor series, which is also a polynomial. But what we try to fit a non-smooth function, such as a triangle-shape function $t(x)$? That is,

$$ t(x) = 1 - |x| $$

(a triangle peaked at $x=0$ and zero at $x=\pm 1$). Let's see:

It is still converging, just much more slowly!