Minimize distance between curves

Question

I have a dataset with values of multiple curves. An example plot is shown below. I want to shift the curves (up/down) so that all curves overlap. This would mean the data points in each curve is scaled up/down by a factor. I am not sure how to frame this as a mathematical problem (to minimize the vertical distance between each pair of curves) and determine the scaling factor for each curve. I tried to start by computing the pair-wise distance matrix but I am not sure what to do next.

Suggestions will be really appreciated.

EDIT:

The following is a sample dataset which includes that datapoints corresponding to 5 curves and coordinate inputs are below

scale =  1.5;
x1 = [0,4,6,10,15,20]*scale;
y1 =  [18,17.5,13,12,8,10];
x2 = [0,10.5,28]*scale;
y2= [18.2,10.6,10.3];
x3 = [0,4,6,10,15,20]*scale;
y3 = [18,13,15,12,11,9.6];
x4 = [9,17,28]*scale;
y4 = [5,5.5,7];
x5 = [1,10,20]*scale;
y5 = [3,0.8,2];

plot(x1,y1, '*-', x2, y2,  '*-', x3, y3,  '*-', x4, y4,  '*-', x5, y5,  '*-')

And the plot looks like below. I tried to do a piecewise interpolation using a linear approximation to find the function values. Since the curves don't have a common x, I am not sure how the target function has to be selected.

Suggestions will be really appreciated. Also, I only want to scale the datapoints $$\tilde f_i(a) = a \cdot f_i$$ and do not want to shift it along the x-dimension (I'm not sure if b $\tilde f_i(a,b) = a \cdot f_i$ + b shifts the function value along the x-dimension; may be I am wrong and b shifts the function value in the y-direction). And I also prefer to define a tolerance window/interval $\epsilon$ (user-defined) above and below the target function which will allow scaling the other curves such that they don't completely overlap on the target curve.

EDIT2: I would like to ask for a suggestion. If I want to scale $f_i$ (move it a bit up) and also scale $g_i$ (move it a bit down) would it make sense if the objective function is expressed as $(a⋅fi − gi)^2 + (b.gi - fi)^2$. Basically, I want to try and move both the curves towards each other instead of moving one curve towards the other.

EDIT3: @Vladislav Gladkikh

I like this idea, this preserves the shape of all $x_k$'s i.e the trend of the $x_k$'s before and after shifting doesn't change and I find this great. Could you please add a pictorial explanation of the third point?

I think overall we do, $f(x_k) - \bar{f(x_k)} + \bar{g(x_k)}$ Here, f is the non-main curve and g is the main curve (i.e y values of the curves at the corresponding x values), and $\bar{f(x_k)}$ is the average.

Then shift all curves toward the main curve so that their averages coincide within the corresponding domain segments.

I think I am a bit confused here since we are subtracting/adding the mean the $\bar{f(x_k)}$'s/ $\bar{g(x_k)}$'s computed at the $x_k$'s in the entire domain. I am not sure if my interpretation of segment is wrong, please correct me if I am wrong.

Could you please clarify if we refer to a segment as the domain between any 2 points that lie on the curve xk, yk here. Could you please illustrate this in a figure, if possible?

Answer 1

Let's assume you have a set of abscissas $x_i$ and two sets of function values on these grid points $f_i, g_i$ representing the functions $f$ and $g$.

As mentioned in the comments, you'll need a model how you adjust the curves, so let's combine the shift-and-scale approach and assume $$\tilde f_i(a,b) = a \cdot f_i + b$$

which corresponds to an affine transformation of function $f$. You then have to minimizr some kind of distance, for which usually the sum-of-squares deviation is used, which yields the problem

$$ \min_{a,b} \left(\sum_i (\tilde f_i(a,b) -g_i)^2\right)\\ = \min_{a,b} \left(\sum_i ( a \cdot f_i + b -g_i)^2\right) $$

As usual you can now set the partial derivatives w.r.t. $a$ and $b$ to zero and thus obtain two equations for the parameter which depend on your data $f_i$ and $g_i$.

This will yield you the optimal shift-and-scale parameters for two selected curves. As you seem to want to overlap any of the five curves in the picture, just repeat the above algorithm for the other three curves.

---

So how does one solve the optimization problem given above? Evaluate the partial derivatives with respect to the parameters and set them to zero: $$ 0 = \partial_a \left(\sum_i ( a \cdot f_i + b -g_i)^2\right) = 2 \sum_i f_i (a \cdot f_i + b -g_i) \\ 0 = \partial_b \left(\sum_i ( a \cdot f_i + b -g_i)^2\right) = 2 \sum_i (a \cdot f_i + b -g_i)\\ $$ and thus obtain the following linear system $$ \begin{pmatrix}\sum_i f_i^2 & \sum_i f_i \\ \sum_i f_i & \sum_i 1\end{pmatrix} \begin{pmatrix}a \\ b\end{pmatrix} = \begin{pmatrix}\sum_i f_ig_i\\ \sum_i g_i\end{pmatrix}\\ $$

You can either solve this system numerically or analytically. The analytical solution can be found, e.g., via the explicit inverse of a 2x2 matrix,

$$ \begin{pmatrix}a \\ b\end{pmatrix} = \frac1{(\sum_i f_i^2)(\sum_i 1) - (\sum_i f_i)^2} \begin{pmatrix} \sum_i 1& -\sum_i f_i \\ -\sum_i f_i &\sum_i f_i^2\end{pmatrix} \begin{pmatrix}\sum_i f_ig_i\\ \sum_i g_i\end{pmatrix} $$

which exists only if the discriminant (the denominator of the prefactor) is not zero.

---

EDIT2: If you only want to scale the curves, things get even easier.

You should then minimize $$ \min_{a} \left(\sum_i ( a \cdot f_i -g_i)^2\right) $$ which leads to a single equation $$ 0 = \partial_a \left(\sum_i ( a \cdot f_i -g_i)^2\right) = 2 \sum_i f_i (a \cdot f_i -g_i) $$ and thus $$ a = \frac{\sum f_i g_i}{\sum f_i^2} $$

Answer 2

Here is a simple solution.

1. Find the curve $(x_m, y_m)$ with the largest domain $x$. In your case it is (x2, y2).
2. Assign it to be the main curve and shift all other (non-main) curves toward it.
3. Each non-main curve $(x_k, y_k)$ should be shifted in such a way that its average value equals the average value of the segment of the main curve within the same domain segment as the domain $x_k$ of the non-main curve.

Because the man curve, $(x_m, y_m)$, has the largest domain, make an interpolation for it (e.g. linear), find the values for this curve at all points for other, non-main curves, $(x_k, y_k)$. Then subtract from each non-main curve, $(x_k, y_k)$, its average, and add the average of the main curve values interpolated at $x$-points of the non-main curve, $(x_k, y_k)$.

It is easier to show in code.

My code is in Julia but you can do analogously in other languages.

using PyPlot, Statistics, Interpolations

scale =  1.5;
x1 = [0,4,6,10,15,20]*scale;
y1 =  [18,17.5,13,12,8,10];
x2 = [0,10.5,28]*scale;
y2= [18.2,10.6,10.3];
x3 = [0,4,6,10,15,20]*scale;
y3 = [18,13,15,12,11,9.6];
x4 = [9,17,28]*scale;
y4 = [5,5.5,7];
x5 = [1,10,20]*scale;
y5 = [3,0.8,2];

plot(x1, y1, "*-", label="y1")
plot(x2, y2, "*-", label="y2")
plot(x3, y3, "*-", label="y3")
plot(x4, y4, "*-", label="y4")
plot(x5, y5, "*-", label="y5")
legend();

This gives your original plot

Then interpolate the curve with the largest domain

interp_linear = LinearInterpolation(x2, y2);

Then shift all curves toward the main curve so that their averages coincide within the corresponding domain segments

plot(x1, y1 .+ mean(interp_linear.(x1)).-mean(y1) , "*-", label="y1")
plot(x2, y2, "*-", label="y2")
plot(x3, y3 .+ mean(interp_linear.(x3)).-mean(y3) , "*-", label="y3")
plot(x4, y4 .+ mean(interp_linear.(x4)).-mean(y4) , "*-", label="y4")
plot(x5, y5 .+ mean(interp_linear.(x5)).-mean(y5) , "*-", label="y5")
legend();

Answer 3

Taking a function as reference $f_r$ the scale-translation transformations for each remaining functions can be handled by minimizing

$$ E(a,b,r) = \sum_{k\ne r}^{m}\sum_{j=1}^n \left(f_k(j)a_k+b_k - f_r(j)\right)^2 $$

or

$$ \min_{a,b,r}E(a,b,r) $$

Here, for each $r$, the minimization can be accomplished by solving a linear system. As the stationary condition is $\nabla E(a,b,r) = 0$ we have

$$ \cases{ \frac{\partial}{\partial a_k}E(a,b,r) = \sum_j^n f_k(j)^2 a_k + \sum_j^n f_k(j)b_k - \sum_j^n f_k(j)f_r(j) = 0\\ \frac{\partial}{\partial b_k}E(a,b,r) = \sum_j^n f_k(j) a_k + \sum_j^n n b_k - \sum_j^n f_r(j) = 0 } $$

so at each $r$ we should solve $m-1$, $(2\times 2)$linear systems or instead one of $(2(m-1),2(m-1)).$

$$ \left(\matrix{\sum_j^n f_k(j)^2 & \sum_j^n f_k(j)\\ \sum_j^n f_k(j) & m}\right) \left(\matrix{a_k \\ b_k} \right)=\left(\matrix{\sum_j^n f_k(j)f_r(j)\\ \sum_j^n f_r(j)}\right) $$

In the proposed case study, the best arrangement is shown in the attached figure. Here the base function is the third from top to bottom.

NOTE

Code attached: MATHEMATICA.

data = {{{157, 60}, {140, 57}, {121, 52}, {103, 47}, {86, 39}, {67, 28}, {50,21}, {32, 15}}, {{159, 70}, {140, 67}, {122, 63}, {105, 58}, {85, 45}, {68, 42}, {50, 34}, {33, 29}}, {{157, 74}, {140, 71}, {121, 70}, {105, 65}, {86, 55}, {69, 52}, {52, 43}, {32, 37}}, {{158, 100}, {141, 100}, {123, 95}, {105, 90}, {88, 78}, {69, 72}, {51, 57}, {33, 46}}, {{160, 143}, {140, 146}, {123, 146}, {105, 145}, {86, 136}, {70, 121}, {51, 97}, {34, 60}}}
gr1 = ListPlot[data, AspectRatio -> 1.3];
gr2 = ListLinePlot[data, AspectRatio -> 1.3];
Show[gr1, gr2]

min = 1000000;
For[r = 1, r <= 5, r++,
 obj = Total[Table[If[k != r, Sum[(data[[k]][[j, 2]] a[k] + b[k] - data[[r]][[j, 2]])^2, {j, 1, 8}], 0], {k, 1, 5}]];
 vars = Flatten[Table[If[k != r, {a[k], b[k]}], {k, 1, 5}]];
 sol = Solve[Grad[obj, vars] == 0, vars][[1]];
 min0 = obj /. sol;
 If[min0 < min, rmin = r; min = min0; coefs = sol]
]

sdata = Table[Table[{data[[k]][[j, 1]], data[[k]][[j, 2]] a[k] + b[k]} /. coefs, {j, 1, 8}], {k, 1, 5}] /. {a[rmin] -> 1, b[rmin] -> 0}

gr1b = ListPlot[sdata, AspectRatio -> 1.3];
gr2b = ListLinePlot[sdata, AspectRatio -> 1.3];
Show[gr1b, gr2b]

Answer 4

At the beginning of the finite element method books, there is a section called vector-function approximations. There is what you are looking for.

Since there is no target function to approach, you can take a weighted average of all curves to find a target $T$ curve. Then find $A$ values to minimize the total difference. Minimize this $(f-A-T)^2$.

Minimization is at the beginning of FEM books, in the approximation section. However, it is the subject of optimization methods and can be found even in Wikipedia.