Is fractal interpolation effective to recover missing data?

Posted on March 28, 2019April 10, 2019 by admin

Introduction

Unfortunately, it is a common occurrence that the data a data analyst has is incomplete.

This missing data has to be filled out a certain way and then the work is done with the “new” data. However, it may happen that the resulting data significantly distorts the real one. Thus, it will obviously lead to a loss in quality for the models based on this data. Therefore, finding the optimal solution to the missing values problem is a very important task. In this article, we will talk about different ways of filling in the missing data and how to compare the results obtained by different methods and on various open datasets.

A few words about interpolation

According to Wikipedia, interpolation is a method of constructing new data points within the range of a discrete set of known data points.

We will take a look at one-dimensional interpolation, which will be in finding the function of one real variable passing through a given set of points in the plane.

Next, we briefly touch upon the basic interpolation methods. Let’s note that the interp1d class is implemented in the scipy.interpolate package to perform one-dimensional interpolation using various methods.

Nearest-neighbor interpolation is the simplest interpolation method. Its essence is that it selects the nearest known function value as an intermediate value:

*Nearest neighborhood interpolation example*

Polynomial interpolation is another interpolation method. It is well known from the higher mathematics course that if there is a set of N points of the plane (with pairwise distinct abscissas), then there is a single polynomial of degree N-1 passing through the indicated points. For example, if there are two points, then there is a unique line (a first-degree polynomial) passing through these points. There are formulas for finding such polynomials (for instance, Lagrange polynomial).

In reality, a situation referred to as “overfitting” will occur for large enough N points. The following illustration demonstrates this (data marked with blue dots and generated from y=sin(x), the red curve is an polynomial interpolation with degree 1, 2, 3 and 7; we see that the best interpolation for polynomial degree 3 and observe overfitting for the polynomial degree 7):

Spline interpolation can be used in order to avoid overfitting. In this interpolation method, at each interval between the interpolation points, the interpolation function is a polynomial, which the degree is limited (for example, in cubic interpolation, this is a polynomial of degree less than 3). Since this kind of interpolation is ambiguous, additional conditions are imposed on the resulting function, for example, differentiability at interpolation points.

*Interpolation methods, determined by the value of the argument “kind” for the class scipy.interpolate.interp1*

Fractal interpolation

Fractal interpolation is another method that was proposed by the American mathematician Michael Barnsley in the mid-1980s. This method is described in a series of his articles and the wonderful book “Fractals Everywhere”. There the interpolation occurs by the so-called fractal self-affine functions. Let us briefly describe the nature of such mathematical objects.

About self-similar, self-affine sets, and fractals

What is a self-similar set? It’s a set which consists of several parts, each of which is a reduced copy of the set itself. Mathematically speaking, each of these parts is an image of the original set with a certain similarity transformation.

Cantor set is the most famous self-similar set (probably because historically it was the very first). It has the following construction. From a unit interval $\left [ 0,1 \right ]$ let’s remove the middle third, which is the interval $\left ( \frac{1}{3},\frac{2}{3} \right )$ . The remaining set is the union of the two intervals. Now, we remove the middle third from each segment. Thus, obtaining a new set, which is the union of four intervals. Then repeat this procedure. By removing the middle thirds of all four intervals, we get a new set. Those points that are never left out of a single segment form a Cantor set.

It is easy to understand that the part of the Cantor set that belongs to the segment $\left [ 0,\frac{1}{3} \right ]$ (like the one that belongs $\left [ \frac{2}{3},1 \right ]$ ) is similar to the original set (i.e. it is obtained from the original set by compression with the coefficient $\frac{1}{3}$ ).

The specified property is a characteristic of the Cantor set. More strictly $f_{1}(x)=\frac{1}{3}x, f_{2}(x)=\frac{1}{3}x+\frac{2}{3}$ – two contraction mappings of $\left [ 0, 1 \right ]$ . Then there’s the unique set $K$ such that $K=f_{1}(K)\cup f_{2}(K)$ . This set is a Cantor set.

However, we are more interested in sets on the plane. Sierpinski carpet is the simplest example of a self-similar set on a plane where a similar procedure, as in the case of a Cantor set, occurs with a unit square:

The first four steps to construct a Sierpinski carpet

The Sierpinski carpet is a self-similar set, for it can be represented as a union of eight reduced copies of itself.

A rather general idea of specifying self-similar sets was probably first described in the article J. Hutchinson “Fractals and self similarity” in 1981. Some time later, Hutchinson’s ideas were developed by M. Barnsley. To some extent, he replaced similarity transformations with affine transformations. As a result, self-affine sets became the natural generalization of self-similar sets. Now such sets are called fractals. An outstanding mathematician Benoit Mandelbrot introduced this term into scientific use in the 70s of the twentieth century. By the fractal, Mandelbrot meant objects that have to some extent self-similar properties. Often, fractal sets are graphs of some continuous functions. Such functions, which graphs are fractal sets, are usually called fractal functions.

Let’s take a look at one of the simplest examples of a fractal function. Construct a sequence of functions on the segment $[0,1]$ . Suppose $f_{0}(x)$ is the diagonal of a unit square (see figure below). Function $f_{1}(x)$ is obtained from $f_{0}(x)$ by replacing a graph segment with a broken line of three links connecting the points $\left ( 0, 0 \right ),\left ( \frac{1}{3}, \frac{2}{3} \right ),\left ( \frac{2}{3}, \frac{1}{3} \right ),\left ( 1, 1 \right )$ . With each of the resulting graph of the function $f_{1}(x)$ repeat the same procedure and get the function graph $f_{2}(x)$ , etc.

The first three steps to construct a Katsuura curve

From the course of mathematical analysis, it is known that the resulting functional sequence has a limit, and it is also a continuous function. Furthermore, this is one of the simplest examples of a continuous nowhere differentiable function, first described in the monograph N. Bourbaki “Functions of a Real Variable”. Now it’s also known as Katsuura function.

The described function is a particular example of self-affine functions. In fact, if we designate $\Gamma$ as a graph of the limit function, then

$\Gamma = f_{1}(\Gamma)\cup f_{2}(\Gamma) \cup f_{3}(\Gamma)$ ,

where

$f_1\left(\begin{bmatrix} x \\ y\end{bmatrix}\right) = \begin{bmatrix} \frac13x \\[3pt] \frac23y\end{bmatrix},$ $f_2\left(\begin{bmatrix} x \\ y\end{bmatrix}\right) = \begin{bmatrix} - \frac13x+\frac23 \\[3pt] \frac13y + \frac13\end{bmatrix},$ $f_3\left(\begin{bmatrix} x \\ y\end{bmatrix}\right) = \begin{bmatrix} \frac13x+\frac23 \\[3pt] \frac23y + \frac13\end{bmatrix}.$

As in the case of self-similar sets, in order to define a self-affine function, it suffices to specify a certain set of contracting affine transformations of the plane. M. Barnsley used this idea, defining fractal interpolation functions. More details on how to define such functions can be found in the article Curve Fitting by Fractal Interpolation (P. Manousopoulos, V. Drakopoulos, T. Theoharis).

The interpolation method is as follows: the data set $D$ matches the set of plane transformations $w_1, w_2, \ldots, w_n \colon \mathbb{R}^2\to\mathbb{R}^2$ of the type $w_i\left(\begin{bmatrix}x\\y\end{bmatrix} \right )=\begin{bmatrix} a_i & c_i \\ 0 & d_i \end{bmatrix}\begin{bmatrix} x \\ y \end{bmatrix}+\begin{bmatrix} b_i \\ e_i \end{bmatrix}.$ Here $d_i$ are free parameters, and the remaining coefficients are uniquely determined by the $d_i$ from the conditions

$w_i\left(\begin{bmatrix}x_0\\y_0\end{bmatrix} \right )=\begin{bmatrix} x_{i-1} \\ y_{i-1} \end{bmatrix}$ and $w_i\left(\begin{bmatrix}x_n\\y_n\end{bmatrix} \right )=\begin{bmatrix} x_{i} \\ y_{i} \end{bmatrix}$ .

The Hutchinson operator $W$ is constructed on compact subsets of $\mathbb{R}^2$ :

$W(A)=\bigcup_{i=1}^nw_i(A)$ .

It is known, that if $w_i$ are contraction mappings, then the operator W is contractive and has a unique fixed point in the space of compact sets of with the Hausdorff metric. Barnsley showed that if $\max_{i=1,...,n}|d_i|<1$ , then mappings are contractive in some metric on $\mathbb{R}^2$ . In this case, the fixed point of the constructed operator, i.e. compact set $G$ satisfying the condition $W(G)=G$ , will be a graph of continuous function f, which interpolates the initial data.

So, if we have a set of $N + 1$ points on the plane, then there is a whole class of fractal interpolation functions interpolating them. Various sets of free parameters will spawn various functions.

Finding methodology of the optimal fractal interpolation function

So, there’s a set of points $(x_0,y_0), (x_1,y_1), \ldots, (x_N,y_N)$ , which completely describes some data. Let’s call this set the initial data. Assume that we only know the values at some n of these points: $(x_{i_0},y_{i_0})=(x_0,y_0),$ $(x_{i_1},y_{i_1}),$ …, $(x_{i_n},y_{i_n})=(x_{N},y_{N})$ . This set of points is called known data. Our goal is to restore initial data according to known data as accurately as possible. Next, assume that $n>3″ alt=”n>3″ align=”absmiddle”>.</div> <p>We proceed from the hypothesis (which we want to check) that fractal functions are more accurate to recover missing data than the functions used in traditional interpolation methods. However, as we saw earlier, there are infinitely many fractal self-affine functions interpolating this data. How do we choose the one that best suits us?</p> <p>There is no answer to the question because the lack of information about the process that generated the data cannot lead us to an unequivocal solution in the task of finding such a function. However, we have constructed an algorithm that constructs a certain function, and we will test this hypothesis on whether it interpolates better than other interpolation functions using the real data. It is obtained as the averaging of several fractal interpolation functions built on different parts of the known data.</p> <p>So, we had the next idea of constructing a fractal interpolation function.</p> <ol> <li> <div>Select several consecutive points of the original data. Let’s call this set of such points a “window”. Let’s call the number of such points the width of the window (we denote it as WINDOW). This value is not less than 4, and the default value is 5 (the value is 5 in the examples described below).</div> </li> <li> <div>We will successively exclude every point of the window, except the first and last, and each time, look for an interpolation function which value in the missing point is closest to the true one.</div> </li> <li> <div>We will calculate the values of each constructed interpolation function in the missing points of the window.</div> </li> <li> <div>At the missing points, the totals are calculated as the arithmetic average of the values of all interpolation functions at that point. The obtained points are smoothed with the <a href=$ Savitzky-Golay filter. For each $x_i$ , $i=0,1,...,N$ we get the predicted values $\hat{y}_i$ .

Now let’s compare the interpolation results on various open datasets. In this case, we will use root-mean-square error to estimate the accuracy of the prediction:

$RMSE = \sqrt{\frac1N\sum_{i=1}^N\left(y_i-\hat{y}_i \right )^2}$ .

Example 1. Initial data: Daily total female births in California in 1959. The initial data set contains 365 points:

Daily total female births in California, 1959

We include randomly selected 20 points in the known data:

The result of the fractal interpolation is presented in the figure:

After applying the Savitzky-Golay filter, we obtain such a smooth curve:

Next table presents the root-mean-square error, mean and standard deviation, calculated for different interpolation methods:

Example 2. The same data set, values at 100 points are known:

The result of fractal interpolation:

and after smooting:

As you can see, the result of smoothing does not meet expectations. However, let’s take a look at the statistical indicators and the error of various interpolation methods:

Example 3. Initial data is a fragment of electrocardiogram data of one of the patients (data is taken from CEBS database).

The data consists of 3333 points. We leave only the uniformly distributed 18 points and interpolate the data at these points:

The results of fractal interpolation (before and after smoothing) are presented in the following figures:

For comparison, let’s take a look at the standard interpolation methods:

Let’s compare the results of this interpolation with the traditional ones:

Example 4. Initial data: a fragment of electrocardiogram data from one of the patients (data is taken from CEBS database). The initial data consisted of 20,000 consecutive values, but we left only every sixth point:

As before, we leave only the evenly distributed 18 points and interpolate the data at these points:

The results of fractal interpolation (before and after smoothing) are presented in the following figures:

Again, we compare the results of this interpolation with the traditional ones:

Conclusion

The examples show slightly better total error values with fractal interpolation than with conventional methods. At the same time, a “miracle” didn’t happen, for the error value, although better, it’s not “way better”.

We also tried to interpolate points of a randomly chosen fractal interpolation function (with random free parameter initialization). The result is worse and is comparable with traditional interpolation methods. The result worsens and becomes comparable with the traditional.

If the data values at each next step in general do not depend on the value at the previous point (as in examples 1 and 2), then it is better to leave out the smoothing. Such a process is just more logical to think of a non-differentiable, and therefore fractal function. If the nature of the data is such that it is convenient to describe them with a continuous and differentiable function, then after the fractal interpolation it is better to perform smoothing.

At the same time, we note the drawbacks of the described algorithm to construct a fractal interpolation function:

Time consuming. The search for the optimal interpolation function is resource-intensive because at each step a search is made for a set of free parameters in which the function best suits the artificially missing value.
The fractal interpolation function that best describes the value at the artificially missing point may result in “overfitting”.

The described cons attract our attention to building the best algorithms to find fractal interpolation functions. It is possible that this will require new restrictions on them. Particularly, such restriction could be the fractal dimension of the resulting graph, i.e. the search for the fractal interpolation function is performed in the class of functions with a predetermined fractal dimension (for example, box-counting dimension). Preventing the fractal dimension of the graph from exceeding the predetermined number can serve as an element of regularization and prevent overfitting.

One thought on “Is fractal interpolation effective to recover missing data?”

Sebastian says:

January 31, 2020 at 12:53 pm

Hy!
I’m very interested in the topic and want to use some of the ideas mentioned here in my PhD research.
Can you maybe send me the programs you used for doing these interpolations?
And do you maybe have some more corresponding references?
All the best Sebastian