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,
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]]

2 Answers 2


Clearly, OP did not even try to read the answer I linked to in my previous answer.

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:


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,

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.

  • $\begingroup$ 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. $\endgroup$ Commented Nov 29, 2015 at 15:17

I had trouble following the step going from

enter image description here


enter image description here

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

enter image description here

Step 2. Apply matrix inversion by parts

J.M. shows how to compute matrix inversion by parts here.

enter image description 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

enter image description here

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.

enter image description here


  • $\begingroup$ 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}$. $\endgroup$ Commented Dec 1, 2015 at 2:24
  • $\begingroup$ @J.M. Thank you. I think I got it now. $\endgroup$ Commented Dec 6, 2015 at 15:32
  • $\begingroup$ Looks right to me. :) $\endgroup$ Commented Dec 9, 2015 at 2:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.