Minimum error linear interpolation of arbitrary curve using fixed number of points

Question

Problem

A finite curve, $r$, defined by $y = f(x)$, is defined by the linear interpolation between the 2 closest points (immediately higher and immediately lower) of a set points $\{(x_0,y_0), (x_1,y_1),...,(x_{N-1},y_{N-1})\}$, where the X coordinates are always increasing (ie: if $i>j \implies x_i>x_j$).

That is

1. the value of $r$ at $x=x_j$ is $y_j$ if $(x_j,y_j)$ belongs to the set of points,
2. the value of $r$ at $x=x_j$ is the linear interpolation of the two closest points (immediately higher, $(x_k,y_k)$, and immediately lower, $(x_i,y_i)$). $$y_j=f(x_j)=y_i+\frac{y_k-y_i}{x_k-x_i} (x-x_j)$$ if $(x_j,y_j)$ does not belong to the set of points, but $x_i<x_j<x_k$ and $k-i=1$
3. the curve is not defined outside $x_0 \leq x \leq x_{N-1}$.

Intuitively, there exists at least 1 reduced curve, $y_r=g(x)$, formed by only $M$ points ($M < N$) that minimizes the error, $E$, when defined as

$$E = \int_{x_0}^{x_{N-1}} (y(x)-y_r(x))^2 dx$$

All values are real values.

Things I have tried

For reference my original curve was around N = 10,000 points, and the reduced curve had to have M = 12 points.

I know there are a few orders of magnitude of difference between both, but by visualizing the curve with Python it looked like it could be done.

I had to solve this problem a few times, so i needed it to be more efficient, but I ended up settling for method 2.

1. The brute-force method.

I started with an algorithm that would create "the initial curve" using ${(x_0,y_0)}$ and ${(x_{N-1},y_{N-1})}$ as the initial and final points, and then divide then pick the $M-2$ remaining points to be equally spaced in X.

That would look like this $$\{(x_0,y_0), (x_i, f(x_i)), (x_j, f(x_j)), ..., (x_{N-1},y_{N-1})\}$$

Then I would iterate over every point changing it a small amount on all 4 directions and compare it with the unchanged version, until the error did not get smaller. Then I might reduce the small amount even further or be done.

This worked OK, but it took forever.

2. The slightly better (third derivative) yet still brute-forced method

I noticed that if I filtered the data and calculated the x-coordinates where the third derivative was zero ($\frac{d^3}{dx^3} = 0$), those points $\{(x_{d_0},f(x_{d_0})), ..., (x_{d_{M-2}},f(x_{d_{M-2}}))\}$, fell closer to the optimized result than the previous method did.

This reduced the number of step to achieve minimum $E$, but it still relied on manually calculating the $E$ over and over.

Question

Is anyone familiar with more "optimal" (aka: less iterative based) approach to solve this (or similar) problem?