# Surface fitting

I have a sparse matrix of points in a volume. I know a surface passing by these points exists and this surface is mostly flat and relatively smooth with some small harmonic contents. See the two images below (this surface of the first one is a bit too complicated and de second one is a bit too simple).

I am looking for a method (Laplacian surfaces / Splines / etc.) to fit my points to the best surface of minimal order. Actually I am expecting do describe my surface with less than 30 coefficients.

Any advice of the method I can use (with Matlab or Mathematica) ? If I can get a working example it would be awesome.  EDIT:

Ok, It seems my expectations were too high. In fact, it is very difficult to fit a 2D surface without a model. I recently tried with a pretty simple model:

$f(x,y) = a_0 + a_1\cdot x^3 + a_2\cdot x^2y + a_3\cdot xy + a_4\cdot y^2x a_5\cdot y^3 + b_1 \sin(b_2\cdot x) + b_3 \cos(b_4\cdot x) + c_1 \sin(c_2\cdot x) + c_3 \cos(c_4 \cdot x)$

With this one I can fit my surface with a least squares approach and normal equations. The result is not bad but it would be better to find a "smart" algorithm that can build the output function with an iterative method. I still haven't find anything like this.

• Is the first one "too complicated" because your data is in fact simpler than this, or because you are willing to accept a worse fit (of smaller order)? If the latter, is 30 your upper-bound for the acceptable order? Aug 12, 2014 at 20:03

Along with Kriging using sparse matrix techniques, you could also look up literature on sparse approximation (aka sparse decomposition) methods or sparse Gaussian graphical models.

• It might be helpful to explain why you think sparse approximation / gausian graphical models would be appropriate.
– Paul
Jul 14, 2014 at 0:32

I would first look into the thin plate splines (TPS), which are widely available in MATLAB. They are simple to use, and were found to perform fairly good in many cases. Here is the documentation: http://www.mathworks.de/de/help/curvefit/tpaps.html?refresh=true

More detailed information is here: http://user.engineering.uiowa.edu/~aip/papers/bookstein-89.pdf

• That MathWorks page is weird... it doesn't ever seem to state over which space the function $f$ is minimized.
– cfh
Apr 22, 2015 at 7:33

I guess you already tried the matlab "fit" function.

f = fit( [x, y], z, 'poly23' )


Looks like it would be a good candidate for Gaussian Process Regression (also known as Kriging).

Edit: I'll explain a little more. Consider the space $X$ of all functions on a grid with $N$ nodes, $X := \mathbb{R}^N$. For example, given a vector $x \in X$, the first component $x_1$ is the value of the function at the first gridpoint, $x_2$ is the value at the second gridpoint, and so forth.

The basic idea is to put a gaussian prior distribution on $X$, $$\pi_\text{prior} \sim e^{-x^T C^{-1} x},$$ and consider making measurements $y_k = x_k$ of the function at some nodes of interest $k$.

We can combine the prior information and the measurement information to construct a posterior distribution on $X$ using Bayes' theorem. Some care must be taken since the likelihood is an improper distribution, but it can be done.

Then one takes the mean of the posterior to get a good function that fits the data. This mean function is equivalent to what you get from Kriging.

Entries of the prior covariance, $C_{ij}$, measure how similar you expect the function to be at nodes $i$ and $j$. A couple typical examples are,

• Exponential decay covariance: The similarity between the function at two points decreases as the points get further apart with decay rate $q$ (eg., $q=2$ is the most common) and characteristic length $\rho$: $$C_{ij} \propto e^{-\frac{||x_i - x_j||^q}{\rho}}.$$ This matrix is a dense $n^2$-by-$n^2$ matrix, which may present some computational difficulties.
• Laplacian covariance: The covariance falls off in the same manner as the greens function of an elliptic PDE, $\Delta^q C(x-x_0) = \delta_x$: $$C \propto \Delta^{-q},$$ where $\Delta$ is the (graph) Laplacian. This covariance is particularly nice since applying $\Delta^{-q}$ can be done efficiently with multigrid. So, you can compute the action of a dense decaying covariance matrix whenever needed (or it's square root), except without actually forming the dense matrix.
• Nick: It might be helpful if you explain why you think kriging would be a good idea.
– Paul
Jul 14, 2014 at 0:30
• @Paul Ok, I added an explanation. Nov 5, 2014 at 2:56