# Disentangling the data

An experiment yields a functional dependence on two variables $F(x,y)$, can be imagined as a 2D map. It is known that the function can be factored as follows $F(x,y)=fx(x) gy(y) +gx(x) fy(y)$. Given the tabulated $F(x,y)$ I would like to numerically find the four unknown functions $fx$, $fy$, $gx$, $gy$. It is clear that in full generality the problem is ill posed and has no unique solution. There are however several constraints that are known a priori:

i) $fx(x)$, $fy(y)$ are non-zero only on a finite interval

ii) $gy(y)$ has a gaussian profile, $gx(x)$ decay exponentially for $x>0$.

Just to give an example let us use the following definitions:

fx[x_] := 1/(Exp[-x - 7] + 1) + 1/(Exp[x - 7] + 1) - 1
gy[y_] := Exp[-y^2]
fy[y_] := (1/(Exp[-y - 10] + 1) + 1/(Exp[y - 10] + 1) - 1) (5 + 0.5 Sin[y])
gx[x_] := Exp[-0.4 x]/(Exp[-5 x] + 5)
F[x_, y_] := fx[x] gy[y] + gx[x] fy[y]


looking as The resulting function has a cross shape

 Plot3D[F[x, y], {x, -15, 15}, {y, -15, 15}, PlotRange -> All, PlotPoints -> {50, 50}] I do not have many ideas how to find $fx,fy,gx,gy$ numerically. I was trying to use a fitting procedure, however, it works like you get what you put in. I am seeking for a more general solution, perhaps with the help of Fourier or wavelet analysis. I would be grateful for any ideas.

• there is a typo here gx[x_] := gx[x_] – chris May 1 '14 at 18:10
• you could use BSplines to parametrize the rest? as in mathematica.stackexchange.com/a/10997/1089 and then use maximum likelihood on the data (which I assume is made of pairs of {x_i,y_i}s ? – chris May 1 '14 at 18:12
• @b.gatessucks and @ chris: thank you for correcting a typo – yarchik May 2 '14 at 7:43
• Just to make the problem statement more precise: experimental data looks like the 3d plot above, i.e. a set of points $x_,y_i, F(x_i,y_i)$; based on this input the 4 single-argument functions (like on the Fig. 1) need to be numerically reconstructed. – yarchik May 2 '14 at 7:46
• A general solution might be quite difficult. Therefore I fell a need to formulate 2 simpler problems: i) numerically find coordinates of the intersection point; ii) numerically verify if the angle between 2 waves is strictly 90 degrees. – yarchik May 2 '14 at 7:53

So, you have a function $F(x,y) = f_x(x)g_y(y) + g_x(x)f_y(y)$, and you want to recover $f_x,g_y,g_x,f_y$.

If you've tabulated the values of $F(x,y)$ in a matrix $\mathbf F$ with entries $f_{ij} = F(x_i,y_j)$, then this amounts to decomposing the matrix as $$\mathbf F \approx \mathbf f_x\mathbf g_y^T + \mathbf g_x\mathbf f_y^T,$$ where $\mathbf f_x,\mathbf g_y,\mathbf g_x,\mathbf f_y$ are column vectors with non-negative entries. (I'm using $\approx$ instead of $=$ because your data presumably has some noise in it.) Stick the vectors into two matrices, $\mathbf W = \begin{bmatrix}\mathbf f_x & \mathbf g_x\end{bmatrix}$ and $\mathbf H = \begin{bmatrix}\mathbf g_y^T \\ \mathbf f_y^T\end{bmatrix}$, and you have $$\underbrace{\mathbf F}_{n\times n} \approx \underbrace{\mathbf W}_{n\times2}\underbrace{\mathbf H}_{2\times n}$$ where all three matrices have non-negative entries. This is precisely the problem of non-negative matrix factorization. And look, there's a Mathematica implementation in the open source Mathematica for Prediction project.

Let's try it!

f = Table[F[x, y], {x, -15, 15, 0.1}, {y, -15, 15, 0.1}];
Needs["NonNegativeMatrixFactorization"];
{w, h} = GDCLS[f, 2];
fx = w[[All, 1]];
gx = w[[All, 2]];
gy = h[[1, All]] // Normal;
fy = h[[2, All]] // Normal;
ListLinePlot[{fx, gx}, PlotRange -> All]
ListLinePlot[{gy, fy}, PlotRange -> All]  There's a bit of cross-talk between the components, and the results appear to be a little different every time you run it (maybe because of random initialization), but overall it looks pretty good.

I am not sure where your problem lies exactly. If you have a set of points, x_i,y_i which obey the PDF F[x,y], you could do maximum likelihood analysis. A parametric model could be

gxa[x_, a_] = Exp[-a x]/(Exp[-5 x] + 5);
G[x_, y_, a_] = fx[x] gy[y] + gxa[x, a] fy[y];


with the corresponding normalization (so that its a PDF)

norma[a_] =
Table[{a,
NIntegrate[G[x, y, a], {x, -20, 0, 20}, {y, -20, 0, 20},
PrecisionGoal -> 2]}, {a, 0.1, 0.9, 0.1}] //
Interpolation[#, a] &;


So that the various PDFs look like this

Table[ContourPlot[G[x, y, a]/norma[a], {x, -8, 8}, {y, -8, 8},
PlotRange -> All, PlotPoints -> {50, 50},
MeshFunctions -> Function[{x, y, z}, z],
ContourStyle -> ColorData[a*10], ContourShading -> False],
{a, 0.1, 0.9, 0.2}] // Show A draw would be drawn from and set of points

sample = RandomVariate[UniformDistribution[{-15, 15}], {10000, 3}];
sample = sample /. {x_, y_, z_} :> {x, y, Abs[z]/10/norma[0.4]};


using the 'keep below PDF' prescription (for a=0.4)

ok = Select[sample, #[] < G[#[], #[], 0.4]/norma[0.4] &];
ok = ok /. {x_, y_?NumberQ, _} -> {x, y}; ok // Length


so that we can check that the maximum likelihood value of our draw corresponds to a=0.4

Table[{a, Plus @@ Log@Map[(G[#[], #[], a]/norma[a]) &, ok]}, {a,
0.1, 0.9, 0.025}] // ListLinePlot which indeed peaks near 0.4.

Note that the above parametric model could be made more general, using e.g. BSplines at the expense of a more complex optimization problem.

Note finally that mathematica has its own MaximumLikelihood function.

• thank you for your answer. It provides a partial solution and demonstrates how difficult a general solution might be. Your solution shows how to find a single parameter when 3 other functions are known as well as the functional dependence on the parameter $a$. I still hope that all 4 functions can be reconstructed by imposing some constraints. – yarchik May 2 '14 at 7:51

Keep in mind this is only a partial solution and it is supposed to give you ideas and some insight.

First to give you an idea of where I am going consider the following example

img = Import["http://i.stack.imgur.com/yV8FW.png"]
WaveletImagePlot[DiscreteWaveletTransform[img]] Now on to the more interesting problem at hand

fx[x_] := 1/(Exp[-x - 7] + 1) + 1/(Exp[x - 7] + 1) - 1
gy[y_] := Exp[-y^2]
fy[y_] := (1/(Exp[-y - 10] + 1) + 1/(Exp[y - 10] + 1) - 1) (5 + 0.5 Sin[y])
gx[x_] := Exp[-0.4 x]/(Exp[-5 x] + 5)
F[x_, y_] := fx[x] gy[y] + gx[x] fy[y]

dwd = DiscreteWaveletTransform[Table[F[x, y], {x, -15, 15}, {y, -15, 15}]]

WaveletMatrixPlot[dwd] Now do you see why I started with an example ?

Now we shall recover the waves

ListPlot3D[InverseWaveletTransform[
WaveletMapIndexed[#1 0.0 &, dwd, {___, 1 | 4}]], PlotRange -> All] ListPlot3D[InverseWaveletTransform[
WaveletMapIndexed[#1 0.0 &, dwd, {___, 2 | 4}]], PlotRange -> All] Or directly from the wavelet transform

ListPlot3D[Abs@Reverse@dwd[{2}, "Values"], PlotRange -> All,
Boxed -> False, ImageSize -> 500] ListPlot3D[Abs@Reverse@dwd[{1}, "Values"], PlotRange -> All,
Boxed -> False, ImageSize -> 500]
` You can and shall always play with the different parameter settings, use different wavelet families, etc.

• It seems there are many solutions to my problem! – yarchik May 12 '14 at 15:00
• Can you explain a little about how your method works? – user484 May 12 '14 at 15:29
• Yes, of course. Well, as simple as I can, you can see the wavelet transform resolves some features differently. Going back to the intersecting lines example we can clearly see a line in every part of the image pyramid and another one that is present only in the bottom half. If analyse two perpendicular lines we would also get the intersection point in the bottom right corner. The crucial point here is that we are actually using 1D wavelets to analyse 2D data. What this means is we can't "cover" the whole data set (2D) with the wavelet, instead we are "covering" in 1D. – Sektor May 12 '14 at 21:26
• @RahulNarain I know it is vague, but decided to keep it simple :D If you have any questions do not hesitate :) – Sektor May 12 '14 at 21:31