How do I apply the Fast Multipole Method to Thin Plate Spline interpolation?

Question

I have n scattered measurements of elevation, z, as a function of x and y coordinate: $ \{(x_i,y_i,z_i)\}_{i=1..n}$ that I want to interpolate so that I find z(x,y) for all x and y.

Using Thin Plate Splines:

$$ z(x,y) = \sum_i l_i\phi(r)\tag{1}$$

where $\ \phi = r^2\ln(r) $ and $ r = \sqrt{(x_i-x)^2+(y_i-y)^2} $.

The interpolation has two steps:

I Find the coefficients

Demand that the interpolated curve goes exactly trough the measurements: $$ z_i = \Phi_{ij}l_j \tag{2}$$ where $ \Phi_{ij}=r_{ij}ln(r_{ij})$ and $r_{ij}=\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}$.

Solve this for l.

II Add up contributions from all "sources"

For any point (x,y) use equation (1) to find z.

This works fine but does not scale well with increasing n since the matrix in (2) is dense. I quickly run out of RAM and solving the matrix equation is slow.

To the rescue The Fast Multipole Method!?

Consider 10 random xy points. Partition the xy space into boxes using a Quadtree:

Approximate the points within a box by the center: $$ l_c= \Sigma_{i} l_i $$ where i is over the points within the box only:

Here the contributions from the 3 "sources" in the box 2<=x<=4, 0<=y<=2 are considered to be far away, R >> r, from the point in question: x=6.8,y=5 and can therefore be approximated by a combined "source" in the center of the box.

Consider one box of size $\Delta x * \Delta y $.

What is the max error from using the center point as a function of distance from point to center, R?

Introduce $$ \beta = r/R \tag{3}$$ Taylor approximation around $ \beta = 0 $: $$ (R+\beta R)^2ln(R+\beta R) = R^2ln(R) + \beta R^2(2ln(R)+1) + O(\beta^2)$$
Thus error of using a "source" in R instead of R + r can be approximated:

$$ \Delta z_i \approx l_i*\beta R^2(2ln(R)+1) $$

Relative error: $$ \Delta z_i/l_i \approx \beta R^2(2ln(R)+1) \tag{4}$$

Demand that relative error is less than some threshold $ \epsilon $ e.g.: 1e-3.

Worst case: all points in one corner: $$ r = \sqrt{(\Delta x/2)^2+(\Delta y/2)^2} \tag{5}$$

With a formulae for the max error I can speed-up step II by the following algorithm:

function visitBox(box):
  r = sqrt(0.25*dx*dx+0.25*dy*dy) 
  beta = r/R
  relError = beta*R*R*(2*ln(R)+1)
  if relError > epsilon:
     foreach child:
        visitBox(child) 
  else
      Z = Z + lc*R*ln(R)
epsilon = 1e-3
Z = 0
visitBox(root)

But how do I use this approach to overcome the problems with step I?

I understand that I can simplify a given row in (1) so that the number of elements are far less than n. However the simplification will vary from row to row and as a result the matrix $ \Phi_{ij} $ will still be nxn. In fact it will be even larger since it must include centroid contributions. If there are m centroids; the size of this matrix will be (m+n)^2. But it will be sparse. And that is the key I guess?

Answer 1

Ok I'll answer why is not good idea to do Barnes-hut style approximation with thin plates.

Barnes-hut is used for faster galaxy simulation so we(in 2d) want to calculate this sum fast:

$$ \phi(x) = - \sum_{t=1}^N m_i \log(\|x-x_i\|)$$

in FMM(in 2d) we use complex numbers so

$$ \phi(x) = -\Re\left( \sum_i m_i \log(x-x_i)\right)$$

where $x,x_i\in \mathbb{C}$.

Now assume that $\|x\|>\|x_i\|$, this corresponds to the situation when you want to evaluate the field caused by points($x_i$) in certain box(which is centered at origin) and you($x$) are far from it, so:

$$\sum_i m_i\log(x-x_i)= $$ $$\sum_i m_i\log(x) +m_i \log(1-\frac{x_i}{x})= $$ $$ \log(x)\sum_i m_i - \sum_i m_i \sum_{n=1}^\infty \frac{x_i^n}{x^nn} =$$ $$\log(x)\sum_i m_i - \sum_{n=1}^\infty \frac1{x^nn}\sum_i m_i x_i^n$$

The idea of FMM is to truncate the infinite sum to $K$ terms and calculate sums $\sum_i m_ix_i^n$ beforehead. So when evaluating $\phi$ you do work proportional to $K$ not to $N$.

so when you want to calculate $\phi$ calculate only

$$\log(x)\sum_i m_i - \sum_{n=1}^K \frac1{x^nn}\sum_i m_i x_i^n$$

Now if you do very rough approximation you pick $K=0$ so you get:

$$\phi(x) = \log(x)\sum_i m_i$$

So this is your $\sum_i l_i$

But now look what happens for thin plates:

$$ f(x) = \sum_i w_i \|x-x_i\|^2 \log \|x-x_i\|$$

using complex numbers you can rewrite it like ($x,x_i \in \mathbb{C}$)

$$ f(x) = \Re\left( \sum_i w_i (x-x_i)(x-x_i)^*\log(x-x_i)\right)$$

Assume that $\|x\|>\|x_i\|$:

$$ \sum_i w_i (x-x_i)(x-x_i)^* \log(x-x_i) = $$ $$ \sum_i w_i (x-x_i)(x-x_i)^* (\log(x)+ \log(1-\frac{x_i}{x})=$$ $$ \sum_i w_i (xx^* -x^*x_i - xx_i^* + x_ix_i^*)\left(\log(x) - \sum_{n=1}^\infty \frac{x_i^n}{x^nn}\right)= $$

Now get rid of the infinite sum(picking $K=0$) as before and you get

$$f(x) \approx \Re\left( \log(x)\sum_i w_i (xx^* -x^*x_i - xx_i^* + x_ix_i^*) \right)$$

but this is different from

$$f(x) \approx \Re\left( xx^*\log(x)\sum_i w_i \right)$$