The method I'll present below, is described (briefly) in the paper "A Revisit to Least Squares Orthogonal Distance Fitting of Parametric Curves and Surfaces" by Liu and Wang. Other methods for solving this problem are also presented in that paper and its references.

We wish to minimize the following sum-of-squares error of the cubic Bezier curve $C(u)$:

\begin{equation} min \sum_{k=0}^{n-1}(C(u_k) - S_k)^2 = min \sum_{k=0}^{n-1}(\sum_{i=0}^{3}P_i B_i(u_k) - S_k)^2 \end{equation}

Where $P_i$ are the Bezier control points to be computed (the output), $B_i(u)$ are the Bezier basis functions, $S_k$ are the input points and $u_k$ are the parameters of the curve corresponding to the input points.

If $u_k$ are given as input, then this can be solved as a linear least square problem. The equation can be interpreted as finding the minimal-norm solution (in the least square sense) of an over-determined set of linear equations of the form:

\begin{equation} \left[ \begin{array}{cccc} B_0(u_0) & B_1(u_0) & B_2(u_0) & B_3(u_0) \\ & \ldots & & \\ & \ldots & & \\ B_0(u_k) & B_1(u_k) & B_2(u_k) & B_3(u_k) \\ & \ldots & & \\ & \ldots & & \\ B_0(u_{n-1}) & B_1(u_{n-1}) & B_2(u_{n-1}) & B_3(u_{n-1}) \end{array} \right] % \left( \begin{array}{c} P_0 \\ P_1 \\ P_2 \\ P_3 \end{array} \right) % \mbox{=} % \left( \begin{array}{c} S_0 \\ \ldots \\ \ldots \\ S_k \\ \ldots \\ \ldots \\ S_{n-1} \end{array} \right) \end{equation}

Denoting $B$ as the left-hand matrix, $P$ as the column vector $(P_0,P_1,P_2,P_3)^T$, and $S$ as the right hand column vector (of size $n$), the standard least squares solution to this over-determined linear system is the solution of: \begin{equation} (B^T B) P = B^T S \end{equation}

However, in the context of the question, $u_k$ are not given and are part of the unknowns. They are the parameters where the minimal distance should be attained. Thus, the problem is non-linear (as @fedja noted in the comments).

A standard heuristic of fitting a curve is to assign (normalized) chord-length values to the $u_k$. While a practical solution, it only gives a coarse approximation to the real minimum attaining curve. However, it can be used as an initial approximation for an iterative process that converges to the minimum. In each iteration $j$, one alternates between solving the linear least square problem to attain the current Bezier curve $C_j(u)$, and then projecting the points $S_k$ onto $C_j$ to attain improved $u_k$ parameters. By executing these two steps iteratively, improved parameters are obtained.

A note on projecting a point onto a Bezier curve (sometimes called point-inversion, see for example Chapter 6.1 of The NURBS Book). Given a point $S_k$ and a curve $C(u)$, the parameter $u_k$ can be attained by finding the roots of the equation $f(u) = C'(u) \cdot (C(u) - S_k)$. In the case of a cubic Bezier curve $f(u)$ is a degree-5 polynomial and numerical methods such as Newton-Raphson (or any other polynomial solving method) can be used to find the root $u_k$ corresponding to the closest projected point.

The following figure demonstrates using this method. The blue crosses are the input points, the orange curve is the initial guess (achieved with chord-length parameter assignment) and the green curve is the result after ten iterations (the start point is fixed as requested in the original question, which means $P_0$ is set and therefore some small modifications to the formulations above were needed).

An interesting thing to notice is that for some configurations there is no minimaal solution. The following figure shows a "corner" configuration where except for the first point, all other points are on a line. The curves in the figure correspond to results after 100, 200 and 300 iterations. It can be seen that, since there is no constraint on the length of the curve, as the curve extends more and more to the left, it will achieve better and better approximation.

The method I'll present below, is described (briefly) in the paper "A Revisit to Least Squares Orthogonal Distance Fitting of Parametric Curves and Surfaces" by Liu and Wang. Other methods for solving this problem are also presented in that paper and its references.

We wish to minimize the following sum-of-squares error of the cubic Bézier curve $C(u)$:

\begin{equation} \min \sum_{k=0}^{n-1}(C(u_k) - S_k)^2 = \min \sum_{k=0}^{n-1}(\sum_{i=0}^{3}P_i B_i(u_k) - S_k)^2 \end{equation}

where $P_i$ are the Bézier control points to be computed (the output), $B_i(u)$ are the Bézier basis functions, $S_k$ are the input points and $u_k$ are the parameters of the curve corresponding to the input points.

If $u_k$ are given as input, then this can be solved as a linear least square problem. The equation can be interpreted as finding the minimal-norm solution (in the least square sense) of an over-determined set of linear equations of the form:

\begin{equation} \left[ \begin{array}{cccc} B_0(u_0) & B_1(u_0) & B_2(u_0) & B_3(u_0) \\ & \vdots & & \\ & \vdots & & \\ B_0(u_k) & B_1(u_k) & B_2(u_k) & B_3(u_k) \\ & \vdots & & \\ & \vdots & & \\ B_0(u_{n-1}) & B_1(u_{n-1}) & B_2(u_{n-1}) & B_3(u_{n-1}) \end{array} \right] % \left( \begin{array}{c} P_0 \\ P_1 \\ P_2 \\ P_3 \end{array} \right) = \left( \begin{array}{c} S_0 \\ \vdots \\ \vdots \\ S_k \\ \vdots \\ \vdots \\ S_{n-1} \end{array} \right). \end{equation}

Denoting $B$ as the left-hand matrix, $P$ as the column vector $(P_0,P_1,P_2,P_3)^T$, and $S$ as the right hand column vector (of size $n$), the standard least squares solution to this over-determined linear system is the solution of: \begin{equation} (B^T B) P = B^T S. \end{equation}

However, in the context of the question, $u_k$ are not given and are part of the unknowns. They are the parameters where the minimal distance should be attained. Thus, the problem is non-linear (as @fedja noted in the comments).

A standard heuristic of fitting a curve is to assign (normalized) chord-length values to the $u_k$. While a practical solution, it only gives a coarse approximation to the real minimum attaining curve. However, it can be used as an initial approximation for an iterative process that converges to the minimum. In each iteration $j$, one alternates between solving the linear least square problem to attain the current Bézier curve $C_j(u)$, and then projecting the points $S_k$ onto $C_j$ to attain improved $u_k$ parameters. By executing these two steps iteratively, improved parameters are obtained.

A note on projecting a point onto a Bézier curve (sometimes called point-inversion, see for example Chapter 6.1 of The NURBS Book). Given a point $S_k$ and a curve $C(u)$, the parameter $u_k$ can be attained by finding the roots of the equation $f(u) = C'(u) \cdot (C(u) - S_k)$. In the case of a cubic Bézier curve $f(u)$ is a degree-5 polynomial and numerical methods such as Newton–Raphson (or any other polynomial solving method) can be used to find the root $u_k$ corresponding to the closest projected point.

The following figure demonstrates using this method. The blue crosses are the input points, the orange curve is the initial guess (achieved with chord-length parameter assignment) and the green curve is the result after ten iterations (the start point is fixed as requested in the original question, which means $P_0$ is set and therefore some small modifications to the formulations above were needed).

An interesting thing to notice is that for some configurations there is no minimal solution. The following figure shows a "corner" configuration where except for the first point, all other points are on a line. The curves in the figure correspond to results after 100, 200 and 300 iterations. It can be seen that, since there is no constraint on the length of the curve, as the curve extends more and more to the left, it will achieve better and better approximation.