Understanding polyharmonic splines

Question

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

Answer 1

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:

$$\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.

Answer 2

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