6.2 Continuous least squares

In discrete least squares, our starting point was a set of data points. Here we will start with a continuous function \(f\) on \([a,b]\) and answer the following question: how can we find the “best” polynomial \(P_{n}(x)=\sum_{j=0}^{n}a_{j}x^{j}\) of degree at most \(n\) that approximates \(f\) on \([a,b]?\) As before, “best” polynomial will mean the polynomial that minimizes the least squares error:

(52)\[\begin{equation} E=\int_{a}^{b}\left(f(x)-\sum_{j=0}^{n}a_{j}x^{j}\right)^{2}dx. \end{equation}\]

Compare this expression with that of the discrete least squares:

\[\begin{equation*} E=\sum_{i=1}^{m}\left(y_{i}-\sum_{j=0}^{n}a_{j}x_{i}^{j}\right)^{2}. \end{equation*}\]

To minimize \(E\) in (52) we set \(\frac{\partial E}{\partial a_{k}}=0\), for \(k=0,1,...,n,\) and observe

\[\begin{align*} \frac{\partial E}{\partial a_{k}} & =\frac{\partial}{\partial a_{k}}\left(\int_{a}^{b}f^{2}(x)dx-2\int_{a}^{b}f(x)\left(\sum_{j=0}^{n}a_{j}x^{j}\right)dx+\int_{a}^{b}\left(\sum_{j=0}^{n}a_{j}x^{j}\right)^{2}dx\right)\\ & =-2\int_{a}^{b}f(x)x^{k}dx+2\sum_{j=0}^na_{j}\int_{a}^{b}x^{j+k}dx=0, \end{align*}\]

which gives the \((n+1)\) normal equations for the continuous least squares problem:

(53)\[\begin{equation} \sum_{j=0}^{n}a_{j}\int_{a}^{b}x^{j+k}dx=\int_{a}^{b}f(x)x^{k}dx \end{equation}\]

for \(k=0,1,...,n.\) Note that the only unknowns in these equations are the \(a_{j}\)’s; hence this is a linear system of equations. It is instructive to compare these normal equations with those of the discrete least squares problem:

\[\begin{equation*} \sum_{j=0}^{n}a_{j}\left(\sum_{i=1}^{m}x_{i}^{k+j}\right)=\sum_{i=1}^{m}y_{i}x_{i}^{k}. \end{equation*}\]

Example 54

Find the least squares polynomial approximation of degree 2 to \(f(x)=e^{x}\) on \((0,2)\).

Solution.

The normal equations are:

\[\begin{equation*} \sum_{j=0}^{2}a_{j}\int_{0}^{2}x^{j+k}dx=\int_{0}^{2}e^{x}x^{k}dx \end{equation*}\]

\(k=0,1,2.\) Here are the three equations:

\[\begin{align*} a_{0}\int_{0}^{2}dx+a_{1}\int_{0}^{2}xdx+a_{2}\int_{0}^{2}x^{2}dx & =\int_{0}^{2}e^{x}dx\\ a_{0}\int_{0}^{2}xdx+a_{1}\int_{0}^{2}x^{2}dx+a_{2}\int_{0}^{2}x^{3}dx & =\int_{0}^{2}e^{x}xdx\\ a_{0}\int_{0}^{2}x^{2}dx+a_{1}\int_{0}^{2}x^{3}dx+a_{2}\int_{0}^{2}x^{4}dx & =\int_{0}^{2}e^{x}x^{2}dx \end{align*}\]

Computing the integrals we get

\[\begin{align*} 2a_{0}+2a_{1}+\frac{8}{3}a_{2} & =e^{2}-1\\ 2a_{0}+\frac{8}{3}a_{1}+4a_{2} & =e^{2}+1\\ \frac{8}{3}a_{0}+4a_{1}+\frac{32}{5}a_{2} & =2e^{2}-2 \end{align*}\]

whose solution is \(a_{0}=3(-7+e^{2})\), \(a_{1}=-\frac{3}{2}(-37+5e^{2}),\) \(a_{2}=\frac{15}{4}(-7+e^{2}).\) Then

\[\begin{equation*} P_{2}(x)=1.17+0.08x+1.46x^{2}. \end{equation*}\]

The solution method we have discussed for the least squares problem by solving the normal equations as a matrix equation has certain drawbacks:

• The integrals \(\int_{a}^{b}x^{i+j}dx=\left(b^{i+j+1}-a^{i+j+1}\right)/(i+j+1)\) in the coefficient matrix give rise to matrix equation that is prone to roundoff error.

• There is no easy way to go from \(P_{n}(x)\) to \(P_{n+1}(x)\) (which we might want to do if we were not satisfied by the approximation provided by \(P_{n}\)).

There is a better approach to solve the discrete and continuous least squares problem using the orthogonal polynomials we encountered in Gaussian quadrature. Both discrete and continuous least squares problem tries to find a polynomial \(P_{n}(x)=\sum_{j=0}^{n}a_{j}x^{j}\) that satisfies some properties. Notice how the polynomial is written in terms of the monomial basis functions \(x^{j}\) and recall how these basis functions caused numerical difficulties in interpolation. That was the reason we discussed different basis functions like Lagrange and Newton for the interpolation problem. So the idea is to write \(P_{n}(x)\) in terms of some other basis functions

\[\begin{equation*} P_{n}(x)=\sum_{j=0}^{n}a_{j}\phi_{j}(x) \end{equation*}\]

which would then update the normal equations for continuous least squares (53) as

\[\begin{equation*} \sum_{j=0}^{n}a_{j}\int_{a}^{b}\phi_{j}(x)\phi_{k}(x)dx=\int_{a}^{b}f(x)\phi_{k}(x)dx \end{equation*}\]

for \(k=0,1,...,n.\) The normal equations for the discrete least squares (43) gets a similar update:

\[\begin{equation*} \sum_{j=0}^{n}a_{j}\left(\sum_{i=1}^{m}\phi_{j}(x_{i})\phi_{k}(x_{i})\right)=\sum_{i=1}^{m}y_{i}\phi_{k}(x_{i}). \end{equation*}\]

Going forward, the crucial observation is that the integral of the product of two functions \(\int \phi_j(x) \phi_k(x) dx\), or the summation of the product of two functions evaluated at some discrete points, \(\sum \phi_j(x_i)\phi_k(x_i)\), can be viewed as an inner product \(\left\langle \phi_{j},\phi_{k}\right\rangle \) of two vectors in a suitably defined vector space. When the functions (vectors) \(\phi_{j}\) are orthogonal, the inner product \(\left\langle \phi_{j},\phi_{k}\right\rangle \) is 0 if \(j\neq k\), which makes the normal equations trivial to solve. We will discuss details in the next section.