Simplify function with polynomial via least-squares

Question

I want to "adjust" (simplify) $f(x)$, a function, by $g(x)$, a polynomial, via least-squares. I want to write code for that. Apperently my code is issuing wrong results, so I was wondering if my mistake lies on the math I had to do in order to simplify the problem into something "codeable" or on the code itself. My "simplification" was as follows:

$r^2 = \int_{x_0}^{x_n}(f(x)-g(x))^2dx$ must be as low as possible, i.e.

$$\nabla \int_{x_0}^{x_n}(f(x)-g(x))^2dx = 0$$

$g(x)$ being a polynomial of degree $n$ then there's a linear system of $n$ partials to be solved, such as

$$\begin{cases}\int_{x_0}^{x_n}\frac{\partial r^2}{\partial c_i}(f(x)-\sum\limits_{j=0}^{n}c_jx^j)^2dx=0\end{cases}$$

$c_j$ the $j$th coefficient of the polynomial we need to determine, with the partial derivative in respect of $c_i$, $i$ being the system's $i$th equation (the $i$th partial).

Solving the partial of $r^2$ (chain rule) gives $\int_{x_0}^{x_n}2r\frac{\partial r}{\partial c_i}rdx=0$. Taking the $2$ out of the integral, solving the partial of $r$ and replacing $r$ gives out

$$\int_{x_0}^{x_n}(f(x)-\sum\limits_{j=0}^{n}c_j x^j)(-x^i)dx=0$$

Since $x^i$ is the only term of $r$ with $c_i$ (and it's $-g(x)$). Taking the $-1$ out of the integral, expanding and adjusting the equation such that $c_j$ gets out of the integral

$$\int_{x_0}^{x_n}x^if(x)dx-\sum\limits_{j=0}^{n}c_j\int_{x_0}^{x_n}x^jx^idx=0$$

Alternatively (so it's easier to write code)

$$\begin{cases}\int_{x_0}^{x_n}x^if(x)dx=\sum\limits_{j=i}^{n+i}c_{j-i}\int_{x_0}^{x_n}x^jdx\end{cases}$$

I tested it with the following case:

$$f(x) = sin(2 \pi x) + cos(3 \pi x) \\ [x_0, x_n] = [-1, 1] \\ n (degree) = 5$$

Using the discrete least-squares approach, with $h=10^{-6}$, that is, 2000000 points, I got

$$g(x) = 14.52221837x^5 - 3.90881979x^4 - 17.86276751x^3 + 3.09711194x^2 + 4.01638587x - 0.25060696$$

Which looks good. Using the continuous approach described above I got

$$g(x) = -0.57270982x^5 - 3.90880894x^4 + 2.17959533x^3 + 3.0971047x^2 - 0.74400025x - 0.25060645$$

Which seems off. So... are there any issues with the math?

Answer 1

I don't see anything wrong with your math.

There is a nice simple description here which agrees with what you did (as far as I can see).

But there is a much simpler solution. If you express your approximation as a linear combination of Legendre polynomials, instead of as a linear combination of $1,x,x^2, \ldots, x^n$, then there is a simple closed-form formula for the coefficients. The Legendre polynomials are orthonormal (with respect to the inner product you're using), so solving the linear equations to get the coefficients becomes trivial. Also, numerical conditioning is much better. The technique is covered in the document cited above.

Answer 2

I think you need to construct a Vandermonde matrix. Given $N$ points on a function, to create a $M$-th order polynomial with least squares you need a [$N$×$M+1$] Vandermonde coefficient matrix.

Given $y_1,\ y_2,\ \ldots \ y_N$ values corresponding to $x_1,\ x_2,\ \ldots \ x_N$ points make $ y = A(x) c$

$$ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} = \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^M \\ 1 & x_2 & x_2^2 & \cdots & x_2^M \\ \vdots \\ 1 & x_N & x_N^2 & \cdots & x_N^M \end{bmatrix} \begin{pmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_M \end{pmatrix} $$

The coefficients are found with the least squares method

$$ y = A c \\ A^\top y = (A^\top A) c \\ c = (A^\top A)^{-1} A^\top y$$

This is also known as a pseudo-inverse.

Your polynomial inerpolation function is now $y(x) = c_0 + c_1 x + c_2 x^2 + \ldots c_M x^M$

Answer 3

I checked your results and they are correct to me. The beauty of the integral approach is that everything can be done analytically; as usual, you end with $n$ linear equations for $n$ unknowns (for a polynomial of degree ($n-1$)). So, your coefficients are in fact $$a_0=-\frac{35 \left(2 \pi ^2-3\right)}{24 \pi ^4}$$ $$a_1=-\frac{105 \left(1485-612 \pi ^2+16 \pi ^4\right)}{256 \pi ^5}$$ $$a_2=\frac{35}{\pi ^2}-\frac{175}{4 \pi ^4}$$ $$a_3=\frac{363825}{128 \pi ^5}-\frac{39375}{32 \pi ^3}+\frac{315}{8 \pi }$$ $$a_4=\frac{1225}{24 \pi ^4}-\frac{175}{4 \pi ^2}$$ $$a_5=-\frac{654885}{256 \pi ^5}+\frac{72765}{64 \pi ^3}-\frac{693}{16 \pi }$$ whose values are identical to those your reported for the discrete approach.

Concerning the second part, as already pointed out by bubba, it seems that there is a mistake somewhere since the curve looks much more to a quartic polynomial. I repeated the same calculations and got the same results. You probably have some errors in the integrals.