# Understanding polyharmonic splines

Rahul Narain and J.M. gives implementation of the thin plate polyharmonic splines in this post. Rahul Narain's method is clear. J.M.'s is convenient to reuse. But J.M. used too many matrix(ls,ws,lx,bb,etc.) that I cannot follow his idea. Who can tell me what's the meaning of the last three lines of coding?

ls = LinearSolve[Φ[N[Function[point, Apply[Abs[#1 + I #2] &, point - tx]] /@ xa,
Precision[data]]]];
ws = ls[wa]; lx = ls[xap];
xap = Transpose[xap]; bb = LinearSolve[xap.lx, xap.ws];
(ws - lx.bb).Φ[EuclideanDistance[vars, #] & /@ xa] + bb.Append[vars, 1]]


In any event: I merely exploited the block structure of the underlying linear system for the polyharmonic spline.

We start with data looking like this:

$$\begin{pmatrix}x_1&y_1&z_1\\x_2&y_2&z_2\\&\vdots&\\x_n&y_n&z_n\end{pmatrix}$$

wa = data[[All, -1]];
xa = Drop[data, None, -1];


This results in

$$\mathtt{xa = }\begin{pmatrix}x_1&y_1\\x_2&y_2\\\vdots&\\x_n&y_n\end{pmatrix}\quad \mathtt{wa = }\begin{pmatrix}z_1\\z_2\\\vdots\\z_n\end{pmatrix}$$

xap = PadRight[xa, {n, p}, 1];


$$\mathtt{xap}=\mathbf W=\begin{pmatrix}x_1&y_1&1\\x_2&y_2&1\\&\vdots&\\x_n&y_n&1\end{pmatrix}$$

ls = LinearSolve[Φ[N[Function[point, Sqrt[Total[(point - tx)^2]]] /@ xa,
Precision[data]]]];


I actually did two things at once here. First, I formed the matrix of RBFs:

$$\mathbf A=\small \begin{pmatrix}0&\phi(\|(x_1,y_1)-(x_2,y_2)\|)&\cdots&\phi(\|(x_1,y_1)-(x_n,y_n)\|)\\\phi(\|(x_1,y_1)-(x_2,y_2)\|)&0&&\vdots\\\vdots&&\ddots&\\\phi(\|(x_1,y_1)-(x_n,y_n)\|)&\cdots&&0\end{pmatrix}$$

and then precomputed the LU decomposition as ls for later use. (I used this construction instead of Outer[], per Leonid.)

Now, we have to solve the block linear system

$$\begin{pmatrix}\mathbf A&\mathbf W\\\mathbf W^\top&\mathbf 0\end{pmatrix}\begin{pmatrix}\mathbf v\\\mathtt{bb}\end{pmatrix}=\begin{pmatrix}\mathtt{wa}\\\mathbf 0\end{pmatrix}$$

ws = ls[wa]; lx = ls[xap];


We then have $\mathtt{ws}=\mathbf A^{-1}\mathtt{wa}$ and $\mathtt{lx}=\mathbf A^{-1}\mathbf W$.

xap = Transpose[xap]; bb = LinearSolve[xap.lx, xap.ws];


We now have $\mathtt{xap}=\mathbf W^\top$, and then we solve the linear system $\mathbf W^\top\mathbf A^{-1}\mathbf W\cdot\mathtt{bb}=\mathbf W^\top\mathbf A^{-1}\mathtt{wa}$. (Again, refer to the math.SE answer on how these expressions were derived.)

Finally, $\mathbf v=$ ws - lx.bb is just $\mathbf A^{-1}\mathtt{wa}-\mathbf A^{-1}\mathbf W\cdot\mathtt{bb}$, and we now form the polyharmonic spline function from the derived coefficients.

• I have followed your link and am unable to get from one particular step to another. I am posting an answer below to show where I am stuck. I will delete the answer after I have an understanding. – Jack LaVigne Nov 29 '15 at 15:17

I had trouble following the step going from to Thanks to J.M. for helping me through the steps.

I thought other users might find the following breakdown helpful.

Step 1. Multiply both sides by the inverse Step 2. Apply matrix inversion by parts

J.M. shows how to compute matrix inversion by parts here. Step 3. Simplify the bottom row

Since we now have a vector on each side of the equality we can equate and simplify the bottom row Step 4. Equate and simplify the top row

In the simplification here the right hand side of the first equation in step three is replaced with bb. QED

• You're supposed to multiply the expression for the partitioned inverse with $\begin{pmatrix}\bf \tt{wa}\\\mathbf 0\end{pmatrix}$, thus resulting in expressions for $\begin{pmatrix}\mathbf v\\\bf \tt{bb}\end{pmatrix}$. – J. M.'s torpor Dec 1 '15 at 2:24
• @J.M. Thank you. I think I got it now. – Jack LaVigne Dec 6 '15 at 15:32
• Looks right to me. :) – J. M.'s torpor Dec 9 '15 at 2:49