Data approximation
==================

In order to build :math:`r`, one can either use a pre-defined sequence of control points :math:`\{ c[k] \}_{k=0,...,M-1}` or of knots :math:`\{ n[k] \}_{k=0,...,M-1}`. Alternatively, one can also attempt to retreive the control points that best approximate a set of data points, as in the classical spline approximation setting.

The problem is framed as follows. We consider a set :math:`\{ p[i] \}_{i=0,...,N-1}` of :math:`N` points to be approximated with the spline model :ref:`(1) <theory:eq:1>` of :math:`M` control points. Hereafter, we will assume a periodic spline model, but a similar derivation can easily be done for the non-periodic case.

We obtain an approximation by ensuring that the samples of the spline :math:`r` match the data points :math:`p`, which translates to

.. math::
   :name: approx:eq:1

   p[i] = \sum_{k=0}^{M-1}c[k]\phi\left(\frac{Mi}{N}-k\right).

Since :math:`\phi` is of finite support, we can re-write :ref:`(1) <approx:eq:1>` as

.. math::
   :name: approx:eq:2

   \mathbf{\Phi}\mathbf{C} = \mathbf{P},

with the basis matrix :math:`\mathbf{\Phi}` (size :math:`N \times M`), the control point matrix :math:`\mathbf{C}` (size :math:`M \times 1`), and the data points matrix :math:`\mathbf{P}` (size :math:`N \times 1`) given by

.. math::
   :name: approx:eq:3

   \mathbf{\Phi} = \begin{bmatrix}
    \phi(0) &  \phi(-1) & \dots & \ \phi(-(M-1)) \\
    \phi\left(\frac{M}{N}\right) &  \phi\left(\frac{M}{N}-1\right) & \dots & \ \phi\left(\frac{M}{N}-(M-1)\right) \\
    \vdots & \vdots & \ddots & \vdots \\
    \phi\left(\frac{(N-1)M}{N}\right) &  \phi\left(\frac{(N-1)M}{N}-1\right) & \dots & \ \phi\left(\frac{(N-1)M}{N}-(M-1)\right)
   \end{bmatrix}

.. math::
   :name: approx:eq:4

   \mathbf{C}  =  \begin{bmatrix}
    c[0] \\
    \vdots  \\
    c[M-1]
   \end{bmatrix}

.. math::
   :name: approx:eq:5

   \mathbf{P}  =  \begin{bmatrix}
    p[0] \\
    \vdots  \\
    p[N-1]
   \end{bmatrix}.

The control points :math:`\mathbf{C}` can then be retrieved by finding the least-square best solution that minimizes

.. math::
   :name: approx:eq:6

   \| \mathbf{P} - \mathbf{\Phi} \mathbf{C} \|^2_2.


Boundary conditions
-------------------

Open splines are padded with additional control points at the ends (:ref:`getting_started/padding:Padding`).
They are located at parameters values :math:`t = -1, -2, \ldots` and :math:`t = M, M+1, \ldots`.
Since the data points are only positions on the parameters interval :math:`[0, M-1]`, the additional control points can cause erratic behavior of the spline when fitting noisy data.
To control this behaviour, one of the following boundary conditions can be enforced:

* *clamped*: :math:`r'(0)=0` and :math:`r'(M-1)=0`
* *natural*: :math:`r''(0)=0` and :math:`r''(M-1)=0`

With out loss of generality, we will describe how to fit a spline with the *clamped* boundary condition.

The spline can be written as:

.. math::
   :name: approx:eq:7

   r(t) = \sum_{k=-p}^{M-1+p}c[k]\Phi(t-k),

where :math:`p` is the amount of padding.

The boundary condition allows us to express the first control point as a combination of the other control points:

.. math::

   r'(0) &= 0 \\
       0 &= \sum_{k=-p}^{M-1+p}c[k]\Phi'(-k) \\
   c[-p] &= \sum_{k=-p+1}^{M-1+p} -c[k] \frac{\Phi'(-k)}{\Phi'(p)}

*Note*: We assume that :math:`\Phi'(p) \neq 0`.

Plugging this into equation :ref:`(7) <approx:eq:7>` and grouping the summand by control point yields:

.. math::
   :name: approx:eq:8

   r(t) = \sum_{k=-p+1}^{M-1+p} c[k] (\Phi(t-k) - \frac{\Phi'(-k)}{\Phi'(p)}\Phi(t+p))

*Note*: :math:`\Phi(t+p)` will be zero for most :math:`t` since the support of the basis function :math:`\Phi` end in the interval :math:`[p, p+1)`.

We can do the same for the last control point, starting from equation :ref:`(8) <approx:eq:8>`:

.. math::

   r'(M-1) &= 0 \\
         0 &= \sum_{k=-p+1}^{M-1+p} c[k] (\Phi'(M-1-k) - \frac{\Phi'(-k)}{\Phi'(p)}\Phi'(M-1+p)) \\

Because of the way we choose :math:`p` and the fact that :math:`M > 1`, we know that is outside the support of :math:`\Phi`
and :math:`\Phi'(M-1+p)=0`.
Therefor we get:

.. math::

          0 &= \sum_{k=-p+1}^{M-1+p} c[k] \Phi'(M-1-k) \\
   c[M-1+p] &= \sum_{k=-p+1}^{M-2+p} -c[k] \frac{\Phi'(M-1-k)}{\Phi'(-p)}

Plugging this into equation :ref:`(8) <approx:eq:8>` yields:

.. math::

   r(t) = \sum_{k=-p+1}^{M-2+p} c[k] (\Phi(t-k) - \frac{\Phi'(-k)}{\Phi'(p)}\Phi(t+p) - \frac{\Phi'(M-1-k)}{\Phi'(-p)}\Phi(t-M+1-p))

Like equation :ref:`(1) <approx:eq:1>`, this equation can be written as a matrix-vector multiplication and can be solved using least-squares.