# Data Interpolation by Near-Optimal Splines with Free Knots Using Linear Programming
## 3 Problem Formulation
### Linear Optimization Formulation for Near-Optimal Solution to (1) and (2)  
This section presents a numerical approach for finding near-optimal splines in maximum norm space (typically referred to as \(L_{\infty}\) space) since we minimize \(\left\|f^{(r)}\right\|_{\infty}=\max\left|f^{(r)}\right|\). The main idea of the approach is to formulate a linear optimization problem, the solution of which approximates \(\left\|f^{(r)}\right\|_{\infty}\), with certain tolerance, of the optimal spline interpolating given \(n\) data points \(\{x_{i},y_{i}\}_{i=1}^{n}\). The idea is presented by using \(n\) equidistant data points (note that in the following development of the formulation, there are no requirements of the point equidistance and the data points are not restricted to be equidistant) such that \(A=x_{1}\), \(B=x_{n}\) and \(x_{i}-x_{i-1}=H\), \(H>0\), \(i=2,\ldots,n\), and \([A,B]=\bigcup_{i=1}^{n-1}[x_{i},x_{i+1}]\). In order to discretize the problem, we chose to divide each interval \([x_{i},x_{i+1}]\), \(i=1,\ldots,n-1\) in \(m\) subintervals of equal width. Thus, we partition interval \([A,\,B]\) into \(l\) evenly spaced subintervals, where \(l=(n-1)m\). Note that the number of points defining the \(l\) subintervals is \(l+1\) and that the points are denoted by \(\left\{\hat{x}_{j}\right\}_{j=1}^{(n-1)m+1=l+1}\), such that the relation \(\hat{x}_{(i-1)m+1}=x_{i}\), \(i=1,\ldots,n-1\) holds between original data points \(x_{i}\) and data points \(\hat{x}_{j}\) created by subdivisions. For example, \(\hat{x}_{1}=x_{1}=A\), \(\hat{x}_{11}=x_{2}\), \(\hat{x}_{21}=x_{3}\),..., \(\hat{x}_{(n-1)m+1=l+1}=x_{n}=B\). It should also be pointed out that the general formula for \(\hat{x}_{j}\) is \(\hat{x}_{j}=B\frac{(j-1)}{l}+A\frac{(1-j+l)}{l}\), so that \(\hat{x}_{1}=A\), \(\hat{x}_{2}=\frac{B}{l}+A\frac{(1+l)}{l}\), \(\hat{x}_{3}=\frac{2B}{l}+A\frac{(2+l)}{l},\ldots,\hat{x}_{l+1}=B\). Additionally, while \(H\) denotes the width of each interval \([x_{i},x_{i+1}]\) specified in the data, \(h\) denotes the width of each smaller interval \([\hat{x}_{j},\hat{x}_{j+1}]\) created after \(m\) subdivisions of each interval. In other words, \(H=mh\). After the interval \([A,\,B]\) is discretized, the derivatives are replaced by their discrete counterparts. Derivatives of a function are often approximated by finite differences. A discrete analog of the \(r\)th derivative of a continuous function is the \(r\)th order difference divided by \(h^{r}\). A general formula of the \(r\)th order _forward_ difference is based on \(r\) consecutive (sub)intervals and is represented as follows [55; 56; 57]  
\[\Delta_{h}^{r}[f]\left(x\right)=\sum_{i=0}^{r}(-1)^{i}\binom{r}{i}f(x+(r-i)h).\]  
_Central_ differences \(\delta_{h}^{r}[f](x)=\sum_{i=0}^{r}(-1)^{i}\binom{r}{i}f\left(x+\left(\frac{ r}{2}-i\right)h\right)\) provide more accurate approximations as compared to forward (or backward) differences: \(\frac{\Delta_{h}^{r}[f]\left(x\right)}{h^{r}}=f^{(r)}(x)+O(h)\) and \(\frac{\delta_{h}^{r}[f]\left(x\right)}{h^{r}}=f^{(r)}(x)+O\left(h^{2}\right)\); thus, errors in forward differences are of the order \(h\), while in central differences, they are of the order \(h^{2}\). Essentially, this is due to considering the values of the function at \(h/2\) distance away both to the left and to the right of the function in central differences rather than the entire \(h\) distance only on one side of the function in forward differences, resulting in a better approximation of the function derivative by central differences.