I have written a graphics library that has a module that plots function lines. Internally, software goes from a starting x-value to an ending x-value in steps that equal a pixel. This produces substantially more points than necessary to faithfully represent the line. An algorithm in the spirit of Ramer–Douglas–Peucker (but much faster) is used to reduce, for example, 248 points to 43 without a loss of accuracy. The latter are then plotted as a polyline (see the first two images below).
The results are quite satisfactory given the objectives of the library. Reducing the points needed by 87% decreases the size of SVGs and significantly improves image generation performance. Wanting to push some optimizations further, I attempted curve fitting. The third image below shows the application of the algorithm (from "Graphics Gems" text) that fits one or more cubic Bezier curves to a polyline. fit-curve npm package was used. The fit line in the image was generated with an error setting of 1 and consists of 16 coordinates.
Having 16 coordinates is certainly an improvement over 43, but there are two main issues. First, for some functions the fitted line has about as many coordinates as the input coordinates it is meant to simplify. And second, fitting is time consuming. On my machine performing 1000 fittings for the line presented takes 2 seconds. The library currently can generate 10,000 graphs shown in Figure 2 in about 1.2 seconds. Adding 20 seconds of fitting processing obviously makes no sense.
Is there anything that I am missing? Is it possible to mathematically derive the Bezier curves? For example, would it make sense to break the line below into two segments (i.e., from origin and down and from origin and up) and try to fit one quadratic Bezier into each segment? When I play with Bezier curves manually, I get an intuitive sense that this may be possible.
Here is another way to frame the problem. The following are the first six coordinates for the polyline in figure 2: [20.9, 284.1], [27, 263], [34, 241.3], [41, 222.4], [48, 205.9], [54, 193.7], [60, 183.1], [66, 173.9]. I can use the following Bezier command as part of the path tag to represent that section almost the same (while using half the characters): M 20.9 284.1 C 32.3 244.6, 43.4 208.6, 66 173.9. I would like, if possible, to non-iteratively (i.e., deterministically) compute that Bezier segment.
Figure 1: Dots that Compose a Function Line
Figure 2: Function Line Plotted from Dots
Figure 3: Bezier Curve Fitted to Dots



In mathematical terms, you’re trying to approximate your original function by a string of cubic polynomials (which you choose to represent in Bézier form).
Approximation by polynomials is an old and rich area of mathematics, so you have many techniques to choose from. I will mention two, specifically:
You can Google these two terms to get started.
For both techniques, there are formulas for the errors involved in the approximations. These formulas typically involve fourth derivatives of your original function, which you will need to estimate (by sampling, typically).