# Probabilistic interpretation of Thin Plate Smoothing Splines

TLDR: Do thin plate regression splines have a probabilistic/Bayesian interpretation?

Given input-output pairs $(x_i,y_i)$, $i=1,...,n$; I want to estimate a function $f(\cdot)$ as follows $$f(x)\approx u(x)=\phi(x_i)^T\beta +\sum_{i=1}^n \alpha_i k(x,x_i),$$ where $k(\cdot,\cdot)$ is a kernel function and $\phi(x_i)$ is a feature vector of size $m<n$. The coefficients $\alpha_i$ and $\beta_i$ can be found by solving $${\displaystyle \min _{\alpha\in R^{n},\beta \in R^{m}}{\frac {1}{n}}\|Y-\Phi\beta -K\alpha\|_{R^{n}}^{2}+\lambda \alpha^{T}K\alpha},$$ where the rows of $\Phi$ are given by $\phi(x_i)^T$ and, with some abuse of notation, the $i,j$'th entry of the kernel matrix $K$ is ${\displaystyle k(x_{i},x_{j})}$. This gives $$\alpha^*=\lambda^{-1}(I+\lambda^{-1}K)^{-1}(Y-\Phi\beta^*)$$ $$\beta^*=\{\Phi^T(I+\lambda^{-1}K)^{-1}\Phi\}^{-1}\Phi^T(I+\lambda^{-1}K)^{-1}Y.$$ Assuming that $k(\cdot,\cdot)$ is a positive definite kernel function, this solution can be seen as the Best Linear Unbiased Predictor for the following Bayesian model: $$y~\vert~(\beta,h(\cdot))~\sim~N(\phi(x)\beta+h(x),\sigma^2),$$ $$h(\cdot)~\sim~GP(0,\tau k(\cdot,\cdot)),$$ $$\beta\propto1,$$ where $\sigma^2/\tau=\lambda$ and $GP$ denotes a Gaussian process. See for example https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2665800/

My question is as follows. Suppose that I let $k(x,x'):=|x-x'|^2 \ln(|x-x'|)$ and $\phi(x)^T=(1,x)$, i.e. thin plate spline regression. Now, $k(\cdot,\cdot)$ is not a positive semidefinite function and the above interpretation doesn't work. Does the above model and its solution still have a probabilistic interpretation as for the case the $k(\cdot,\cdot)$ is positive semidefinite?

• You seem to assume the $x$ is in a $d$-dimensional space with $d=2$ or at least that the integer $d$ is even.
– Yves
Mar 10, 2018 at 10:03
• Ok, so what are the implications?
– MthQ
Mar 11, 2018 at 7:28
• This was just a side remark because in the question one may think that $x_i$ are scalars. But in this case Duchon's kernel has the form $|x - x'|^{2m-1}$ with $m$ integer, and $m=2$ for the usual smoothing spline. I think that the probabilistic interpretation remains nearly unchanged but the GP is non-stationary: it is an Intrinsic Random Function. For the usual smoothing spline, this turns out to be an integrated Wiener process.
– Yves
Mar 11, 2018 at 8:56
• @Yves that sounds interesting. You may want to expand your comment to an answer, explaining a bit more what an intrinsic random function is, and adding the classic example of the smoothing spline. If you worry about proving that the TPS kernel gives rise to a non-stationary GP, maybe a simulation could be a useful compromise, especially if you add a non parametric estimate of the variance of the posterior predictive distribution. Mar 12, 2018 at 7:37
• @DeltaIV. Thank you. I will try to do it, not an easy task yet. I am quite sure this holds when the functions $\phi_j$ are suitable polynomials related to the kernel, but this might no longer be true with arbitrary $\phi_j$ as in the more classical GP context.
– Yves
Mar 12, 2018 at 10:34

Let the model of the question be written as $$$$\tag{1} Y_i = \boldsymbol{\phi}(\mathbf{x}_i)^\top\boldsymbol{\beta} + h(\mathbf{x}_i) + \varepsilon_i$$$$ where $$h(\mathbf{x})$$ is an unobserved GP with index $$\mathbf{x} \in \mathbb{R}^d$$ and $$\varepsilon_i$$ is a normal noise term with variance $$\sigma^2$$. The GP is usually assumed to be centered, stationary and non-deterministic. Note that the term $$\boldsymbol{\phi}(\mathbf{x})^\top \boldsymbol{\beta}$$ can be regarded as a (deterministic) GP with kernel $$\boldsymbol{\phi}(\mathbf{x})^\top \mathbf{B}\, \boldsymbol{\phi}(\mathbf{x})$$ where $$\mathbf{B}$$ is a infinite-valued'' covariance matrix. Indeed, by taking $$\mathbf{B} := \rho \, \mathbf{I}$$ with $$\rho \to \infty$$ we get the kriging equations of the question. This is often named the diffuse prior for $$\boldsymbol{\beta}$$. A proper posterior for $$\boldsymbol{\beta}$$ results only when the matrix $$\boldsymbol{\Phi}$$ has full rank. So the model writes as well as $$$$\tag{2} Y_i = \zeta(\mathbf{x}_i) + \varepsilon_i$$$$ where $$\zeta(\mathbf{x})$$ is a GP. The same Bayes interpretation can be used with restrictions when $$\zeta(\mathbf{x})$$ is no longer a GP but rather is an Intrinsic Random Function (IRF). The derivation can be found in the book of G. Wahba. Readable presentations of the concept of IRF are e.g. in the book by N. Cressie and the article by Mardia et al cited below. IRFs are similar to the well-known integrated processes in the discrete-time context (such as ARIMA): an IRF is transformed into a classical GP by a kind of differencing operation.

Here are two examples of IRF for $$d=1$$. Firstly, consider a Wiener process $$\zeta(x)$$ with its initial condition $$\zeta(0) = 0$$ replaced by a diffuse initial condition: $$\zeta(0)$$ is normal with an infinite variance. Once a value $$\zeta(x)$$ is known, the IRF can be predicted as is the Wiener GP. Secondly, consider an integrated Wiener process given by the equation $$\text{d}^2 \zeta(x) / \text{d}x^2 = \text{d} W(x)/\text{d}x$$ where $$W(x)$$ is a Wiener process. To get a GP we now need two scalar parameters: two values $$\zeta(x)$$ and $$\zeta(x')$$ for $$x \neq x'$$, or the values $$\zeta(x)$$ and $$\text{d}\zeta(x) / \text{d}x$$ at some chosen $$x$$. We may consider that the two extra parameters are jointly Gaussian with an infinite $$2 \times 2$$ covariance matrix. In both examples, as soon as a suitable finite set of observations is available, the IRF is nearly coped with as a GP. Moreover we used a differential operator: $$L := \text{d}/ \text{d}x$$ and $$L := \text{d}^2/ \text{d}x^2$$ respectively. The nullspace is a linear space $$\mathcal{F}$$ of functions $$\phi(x)$$ such that $$L \phi = 0$$. It contains the constant function $$\phi_1(x)=1$$ in the first case and the functions $$\phi_1(x)=1$$ and $$\phi_2(x) = x$$ in the second case. Note that in the first example $$\zeta(x) - \zeta(x + \delta)$$ is GP for any fixed $$\delta$$ in the first example and similarly $$\zeta(x-\delta) - 2 \zeta(x) + \zeta(x + \delta)$$ is a GP in the second case.

For a general dimension $$d$$, consider a linear space $$\mathcal{F}$$ of functions defined on $$\mathbb{R}^d$$. We call an increment relative to $$\mathcal{F}$$ a finite collection of $$s$$ locations $$\mathbf{x}_i \in \mathbb{R}^d$$ and $$s$$ real weights $$\nu_i$$ such that $$\sum_{i=1}^s \, \nu_i \,\phi(\mathbf{x}_i) = 0 \text{ for all } \phi \in \mathcal{F}.$$ Consider $$\mathcal{F}$$ as being the nullspace of our examples. For the first example we can take e.g. $$s=2$$ with $$x_1$$ and $$x_2$$ arbitrary and $$[1, \, -1]$$. For the second example we can take $$s = 3$$ equally spaced $$x_i$$s and $$\boldsymbol{\nu} = [1,\,-2,\,1]$$. The definition of an IRF involves a space of functions $$\mathcal{F}$$ and a function $$g(\mathbf{x}, \, \mathbf{x}')$$ which is conditionally positive w.r.t. $$\mathcal{F}$$, which means that $$\sum_{i=1}^s \sum_{j=1}^s \nu_i \nu_j \, g(\mathbf{x}_i, \, \mathbf{x}'_j) \geq 0$$ holds as soon as $$[\nu_i,\,\mathbf{x}_i]_{i=1}^s$$ is an increment w.r.t. $$\mathcal{F}$$. From $$\mathcal{F}$$ and $$g(\mathbf{x},\,\mathbf{x}')$$ we can make a covariance kernel hence a GP as in Mardia et al. We can start from a linear differential operator $$L$$ and use the nullspace as $$\mathcal{F}$$; the IRF will then have connection with the equation $$L \zeta =$$ a Gaussian noise.

The computation of the prediction of the IRF is nearly the same as in the question, with $$k(\mathbf{x},\,\mathbf{x}')$$ replaced by $$g(\mathbf{x},\,\mathbf{x}')$$, but with the $$\phi_i(\mathbf{x})$$ now forming a basis of $$\mathcal{F}$$. The extra constraint $$\boldsymbol{\Phi}^\top \boldsymbol{\alpha} = \mathbf{0}$$ must be added in the optimisation problem, which will grant that $$\boldsymbol{\alpha}^\top \mathbf{K} \boldsymbol{\alpha} \geq 0$$. We still can add more basis functions which are not in $$\mathcal{F}$$ if needed; this will have the effect of adding a deterministic GP, say $$\boldsymbol{\psi}(\mathbf{x})^\top\boldsymbol{\gamma}$$ to the IRF $$\zeta(\mathbf{x})$$ in (2).

The thin-plate spline depends on an integer $$m$$ such that $$2m> d$$, the space $$\mathcal{F}$$ contains polynomials with low degree, with dimension $$p(m)$$ depending on $$m$$ and $$d$$. It can be shown that if $$E(r)$$ is the following function for $$r \geq 0$$ $$E(r) := \begin{cases} (-1)^{m + 1 + d /2} \, r^{2m-d} \log r & d \text{ even},\\ r^{2m-d} & d \text{ odd,} \end{cases}$$ then $$g(\mathbf{x},\,\mathbf{x}') := E(\|\mathbf{x} - \mathbf{x}'\|)$$ defines a conditionally positive w.r.t. $$\mathcal{F}$$. The construction relates to a differential operator $$L$$. It turns out that for $$d=1$$ and $$m=2$$ the thin plate spline is nothing than the usual natural cubic spline, which relates to the integrated Wiener example above, with $$g(x,\,x') = |x - x'|^3$$. So (2) is nothing then than the usual smoothing spline model. When $$d=2$$ and $$m=2$$ the nullspace has dimension $$p(m)=3$$ and is generated by the functions $$1$$, $$x_1$$ and $$x_2$$.

Cressie N Statistics for Spatial Data. Wiley 1993.

Mardia KV, Kent JT, Goodall CR and Little JA. Kriging and splines with derivative information. Biometrika (1996), 83,1, pp. 207-221.

Wahba G Spline Models for Observational Data. SIAM 1990.

Wang, Y Smoothing Splines, Methods and Applications. Chapman and Hall, 2011.

• Thanks a lot for the effort that you put in. Extremely useful. I have a one additional question. So adding additional basis functions to $\boldsymbol{\phi}(\cdot)$ (on top of the basis functions of $\mathcal{F}$) does not change the interpretation of $\zeta(\cdot)$. What I did notice, however, is that the solution $\alpha^*$ given in my question above, always satisfies $\boldsymbol{\Phi}^\top \boldsymbol{\alpha} = \mathbf{0}$, not only if $\boldsymbol{\phi}(\cdot)\in \mathcal{F}$. How can this be interpreted?
– MthQ
Mar 19, 2018 at 10:26
• Yes. In both cases there are $n+p$ basis functions in the approximation of $f(x)$, while only $n$ observations are used. So we have something like a rank-deficient regression with coefficients $\beta_i$ and $\alpha_j$. Since the $\beta$ part is not penalised, it tends to 'absorb' more of the variation of $y$ than does the $\alpha$ part which brings $p$ linear constraints. Note that nothing forbids the use of some of the $n$ "kernel shifts" functions $x \mapsto k(x, x_i)$ as $\phi_j(x)$. If we use all of them then all $\alpha_j^\star$ are zero, which seems sensible.
– Yves
Mar 19, 2018 at 16:34
• Probably mean $m \ge 2d$ given the cubic spline example? Sep 18, 2021 at 17:25
• @mathtick Thank you for pointing out this problem: the condition is $2m > d$.
– Yves
Sep 18, 2021 at 17:40