Definition

A spline function \(s(x)\) of degree \(k > 0\) (order \(k+1\)) is defined on the finite interval \([a,b]\) containing monotonically increasing knots \(\lambda_j\) where \(j=0,1,\cdots,g+1\), \(\lambda_0 = a\), and \(\lambda_{g+1}=b\). For each knot interval \([\lambda_j, \lambda_{j+1}]\) the spline function is given by a polynomial of degree \(\leq k\), and must be continues along with all its derivatives up to \(k-1\) on \([a,b]\). For further information on the theory and computational implementation for univariate and bivariate spline fitting, please see Paul Dierckx, Curve and Surface Fitting with Splines, Oxford University Press, 1993.

Implementation - Spline representation

Spline coefficients and knots are derived iteratively starting with only two knots per dimension by default at \([min(x), max(x)]\). At each iteration, a new knot will be added until:

\begin{eqnarray} \left( \sum_i^n{(d_i - s(x_i))^2} \right) - \alpha \leq \alpha * tolerance \end{eqnarray}

where \(d_i\) is the sample value at point \(i\) and \(\alpha\) is the smoothing factor for the reduction. When \(\alpha = 0\), the reduction will terminate once the resulting function represents an interpolating spline. The \(tolerance\) parameter is a user defined value, typically of the order \(10^{-3}\).

During each new reduction, the current knots are used to divide the entire sample space into panels. Each panel is defined as an n-dimensional volume with neighbouring knots at each vertex. Fits are performed to all samples in a panel, and the panel with the greatest \(w = \sum{(d - s(x))^2}\) is taken as the candidate in which to place a new knot. If no valid knot can be placed at that location, the algorithm moves onto the next greatest panel \(w\) and so on.

Note that in rare circumstances, iterations will halt when the location of the new knot coincides with an existing one or the maximum number of knots has been reached. Iterations may also terminate once the maximum number of iterations has been reached, the solution no longer converges, or the required storage space exceeds a certain limit.

Implementation - Spline evaluation

Once a spline representation has been derived, it may be evaluated at any given coordinate \(x\) in the range \(a \leq x \leq b\). B-splines are generated for \(x\) in each dimension using the Cox-de Boor recursion formula (2007):

\begin{eqnarray} B_{j,0}(x) & = & \left\{ \begin{matrix} 1 & \text{if} \quad \lambda_j \leq x < \lambda_{j+1} \\ 0 & \text{otherwise} \end{matrix} \right\} \\ B_{j,k}(x) & = & \frac{x - \lambda_j} {\lambda_{j+k} - \lambda_j} B_{j,k-1}(x) + \frac{\lambda_{j+k+1} - x} {\lambda_{j+k+1} - \lambda_{j+1}} B_{j+1,k-1}(x) \end{eqnarray}

If \(j\) defines the knot interval \(\lambda_j \leq x < \lambda_{j+1}\), a solution at \(x\) is given by:

\begin{eqnarray} s(x) = \sum_{i=j-k}^{j} c_i B_{i,k}(x) \end{eqnarray}

where \(c\) are the spline coefficients which are determined during the spline representation phase.