fitting - 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.

Answer

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.