The easy explanation goes as follows: for each missing value take the average over the 4 surrounding values, i.e. of the entries above, below, left and right. This holds for the inner points of the matrix, special rules apply to the borders, see below. Of course, maybe some or even all of the surrounding values are missing, too. In general this leads to a system of linear equations one has to solve.

The more involved explanation: search for a harmonic function interpolating the non-missing data. A harmonic function satisfies

from which the mean value property follows. This says that each function value is the average of the function values on or in a surrounding sphere

where . That's deeper math of course. To find the harmonic function we have to solve its defining partial differential equation, which we discretize on the grid given by the data matrix itself assuming equidistant nodes (with distance say ). For an inner point we have

which is the same as

matching the easy explanation from the beginning. For each inner point of the matrix we either get this averaging equation, or if the data point is known to be simply

This way we collect equations, if the original data matrix has rows and columns. The borders have modified equations, which I do not list here explicitly, but which are based on the conditions for the upper and lower border, for the left and right border and for the four corners.

Here is a tentative implementation in QuantLib using the BiCGStab solver at its heart.

Let’s do some examples (the code can be found here, together with a gnuplot script generating the pictures.

The original data is generated by filling a matrix sample like this

    for (Size i = 0; i < N; ++i) {
        for (Size j = 0; j < N; ++j) {
            sample(i, j) =
                sin((double)(i + 1) / N * 4.0) * cos((double)(j + 1) / N * 4.0);
        }
    }

By the way, the function generating the test data is not harmonic. Plotted as a heatmap it looks like this:

Now we delete of the points, randomly.

We reconstruct the missing values by applying Laplace interpolation, which is done as easy as this

    laplaceInterpolation(sample2);

where sample2 is the destructed matrix with missing values marked as Null. The function laplaceInterpolation replaces these values by interpolated ones and writes them back to the input matrix sample2.

The reconstructed matrix looks is this:

As beautiful as the original, isn’t it. Let’s delete more points and look and the reconstruction.

For we get

for it is

for

and for

Quite amazing. Let’s look at how well the method works if we do not delete points randomly, but erase a whole stripe at the center of the matrix.

If the stripe is of the height

for

and for

What about extrapolation? So let’s delete a stripe at the top of the matrix.

For we get

for

and for

Of course the method can not see what is no longer there. Still it is working surprisingly well if you ask me.