User guide

Theory

  • Multivariate splines

Spline interpolation is well studied in the 1-dimensional case – [HTF09] is a good reference. The natural cubic splines, in particular, arise as the optimal estimator for the least squares problem

(1)\[\sum_{i=1}^N \left( y_i - f(x_i) \right)^2 - \lambda \int_{-\infty}^\infty \left( f''(x) \right)^2 dx,\]

where \((y_i, x_i)_{i=1}^N\) are the sample points and \(f\) is a candidate estimator. The penalty \(\lambda > 0\) is on the convexity of the estimator, which prevents oscillation. The knots of such cubic spline are exactly \((x_i)_{i=1}^N\).

The criterion (1) can be generalised to multiple dimensions, although the maths become quite involved – see [AdST04]. The generalisation is not only for multiple dimensions, but also on regularity, including different orders of differentiation. It is shown that the solution to the generalised version of (1) is of the form

(2)\[\begin{split}\begin{align*} f(x) & = \sum_{i=1}^N \beta_i K_\nu \left(\lVert x - x_i \rVert\right) + p(x), \\ K_\nu(r) & = \begin{cases} r^\nu, & \nu \text{is odd,} \\ r^\nu \log r, & \nu \text{is even,} \end{cases} \end{align*}\end{split}\]

where \(\nu\) is linked to the regularisation parameters, \((\beta_i)_{i=1}^N\) are coefficients to be found and \(p\) is a polynomial with degree \(\nu - 1\). This function \(f\) is called a polyharmonic spline. Although it is not a trivial task, the 1D natural cubic splines can be shown to a particular case of (2). The coefficients \((\beta_i)_{i=1}^N\) and polynomial \(p\) can be found via a linear system of equations.

The kernel form (2) can be extended to other kernels – e.g. multiquadric, Gaussian – with \(p \equiv 0\). Of course, this is not anymore the solution of the generalised version of the criterion (1), but can be useful in practice.

  • Hermite splines

It is useful to highlight two types of interpolation: Lagrange-type and Hermite-type. The Lagrange-type interpolation is the one treated in the previous section, where we only have information of levels of the true function \((y_i, x_i)_{i=1}^N\). Interpolation of the Hermite type uses, in addition, information about the derivative of the true function \((y_i, (\partial_j y_i)_{j=1}^M, x_i)_{i=1}^N\). The corresponding generalisation of criterion (1) for the Hermite case yields the solution

(3)\[\begin{split}\begin{align*} f(x) & = \sum_{i=1}^N \beta_{i, 0} K_\nu \left(\lVert x - x_i \rVert\right) + \sum_{i=1}^N \sum_{j=1}^M \beta_{i, j} \partial_j K_\nu \left(\lVert x - x_i \rVert\right) + p(x), \\ K_\nu(r) & = \begin{cases} r^\nu, & \nu \text{is odd,} \\ r^\nu \log r, & \nu \text{is even,} \end{cases} \end{align*}\end{split}\]

where, compared to (2), the partial derivatives of the kernels are added. As before, the coefficients \(((\beta_{i, j})_{j=1}^M)_{i=1}^N\) and polynomial \(p\) can be found via a linear system of equations.

The kernel form (3) can be extended to other kernels with \(p \equiv 0\). Our implementation is based in the MATLAB code available in [Fas07].

  • Example

In this example, we compare the natural cubic spline interpolation – which corresponds to the 1-D polyharmonic spline with \(\nu = 3\) – and the Hermite spline interpolation using the multiquadric kernel.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from parma import polyharmonic_interpolator, polyharmonic_hermite_interpolator


def f(x):
    return np.sin(x)

def df(x):
    return np.cos(x)

data_locs = -2/3 + 2*np.arange(3)/3
data_vals = f(data_locs)
data_hermite_vals = df(data_locs)

x = np.linspace(-1, 1, 20)
hermite_splines, _ = polyharmonic_hermite_interpolator(
    (data_locs,), data_vals, (0,), (data_hermite_vals,), 5)
lagrange_splines = polyharmonic_interpolator((data_locs,), data_vals, 2)

fig, ax = plt.subplots()
pd.Series(f(x), x, name='sin').plot(ax=ax)
pd.DataFrame([lagrange_splines((x,)), hermite_splines((x,))],
             ['lagrange', 'hermite'], x).T.plot(ax=ax, linestyle='--')
pd.Series(data_vals, data_locs, name='samples'
          ).plot(ax=ax, marker='o', color='k', linewidth=0)
ax.legend()

(Source code)

Sampling

When the true function is costly to compute, it is important to find a scheme that minimises the number of sampled points. Mesh grids can be applied when the dimensionality is low, otherwise the number of points in the grid grows exponentially with the number of dimensions. For the high-dimensional cases, simplexes could be applied.

A simplex is a geometrical object that generalises the notion of a triangle and tetrahedron to arbitrary dimensions. The important property to be exploited is that simplexes in \(\mathbb R^n\) have \(n + 1\) vertices. Therefore, by sampling the vertices of a simplex, the number of samples grow linearly with the number of dimensions.

In the example below, we start with correlated data points in 2 dimensions. We orthonormalise the data and construct a standard simplex in this orthonormal space. In order to obtain the sample points, we invert the orthonormalisation. The origin is also added.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from parma.regularisation import orthonormalise_data
from parma.utils import standard_simplex, right_inverse


z = np.random.randn(1000, 2)
x = z @ np.array([[1, 0], [0.75, 0.66]]).T

mean_vec = np.mean(x, axis=0)
cov_mat = np.cov(x.T)
total_var = np.trace(cov_mat)
orth_mat = orthonormalise_data(cov_mat)
simplex = (np.sqrt(total_var)*standard_simplex(2)
           @ right_inverse(orth_mat).T + mean_vec)

fig, ax = plt.subplots()
pd.DataFrame(*x.T[::-1]).plot(marker='o', alpha=0.25, linewidth=0,
                              legend=False, ax=ax)
pd.Series(*np.append(simplex, simplex[:1], axis=0).T[::-1]).plot(ax=ax, c='k')
pd.Series([0], [0]).plot(ax=ax, c='k', marker='o')

(Source code)

References

AdST04

Rémi Arcangeli, María Cruz López de Silanes, and Juan José Torrens. Multidimensional minimizing splines: Theory and applications. Springer Science & Business Media, 2004.

Fas07

Gregory E Fasshauer. Meshfree approximation methods with MATLAB. Volume 6. World Scientific, 2007.

HTF09

Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction. Springer Science & Business Media, 2009.