The problem of curve fitting is fundamental to font technology, as we want to make Béziers which most closely resemble the "true" shape of the glyph. Font tools need to apply curve fitting to simplify outlines, apply transformations such as offset curve, delete a smooth on-curve point, and other applications. However, the problem is surprisingly tricky and there is no definitive solution in the literature. I recently worked out a much better solution and wrote about it in my blog: https://raphlinus.github.io/curves/2021/03/11/bezier-fitting.html
The main reason it's so tricky is that C-shaped Bézier curves come in triplets of similar shape, causing local minima so that an optimizer may fall into one of these minima, even if one of the other two has lower error. My blog post gives a relatively simple formula for finding all three minima.
One subtle detail is that the goal of optimizing Bézier curves for shape accuracy is in tension with making them interpolation compatible. To me, this is a strike against cubic Béziers as a fundamental shape representation for fonts, and I continue to explore techniques that have smoothness guarantees by construction. But hopefully that's a topic for a future blog post.