# Decompose many-variables function as sum of single-variable functions - Schmidt decomposition

Is there any built-in function for decomposing a many-variables function as sum of products of single-variable functions? (I usually call this procedure Schmidt decomposition, but I'm not 100% that it is the correct name, and it is in some way analogous to the singular value decomposition in linear algebra)

$f(x,y)=\sum_k g_k u_k(x) v_k(y)$

where the functions u and v should be orthonormal.

Here a few examples of what I mean:

If the function is: $f(x,y)=\exp(-x^2-y^2)$ the output should be:

$g_1=1$

$u_1=\exp(-x^2)$

$v_1=\exp(-y^2)$

If the function is: $f(x,y)=\exp(-x^2-y^2)(x-y)$ the output should be:

$g_1=\frac{1}{\sqrt{2}}$

$u_1=\exp(-x^2)\ x$

$v_1=\exp(-y^2)$

$g_2=\frac{1}{\sqrt{2}}$

$u_2=\exp(-x^2)$

$v_2=-\exp(-y^2)\ y$

• Why aren't $g_1$ and $g_2$ just 1 in the second example? – Carl Woll Jul 11 '17 at 16:36
• that's not really crucial: the important thing is that the scaling is the same. $1/\sqrt{2}$ derives from the normalization of the coefficients: $\sum_k |g_k|^2 =1$ – Fraccalo Jul 11 '17 at 19:01
• Well, your expected outputs don't satisfy your equation $f(x,y)=\sum_k g_k u_k(x) v_k(y)$, hence my question. – Carl Woll Jul 11 '17 at 19:08
• yes you are right, I put those weight because I'm used to use those objects after having normalised them, but that's a step further :) – Fraccalo Jul 12 '17 at 8:16

Here is a variant of @LeonidShifrin's answer to How can I separate a separable function:

separate[f_, vars_:Automatic] := Module[{v, res, other, free},
v = Replace[vars,
{
Automatic->ReduceFreeVariables[f],
Except[_List]->{vars}
}
];
Integrate,
{Simplify @ D[Log@f,{v}], v}
];
other = Activate @ Thread[Inactive[DeleteCases][Alternatives @@ v, v]];
res = Replace[res Boole@MapThread[FreeQ,{res, other}], 0->1, {1}];
Prepend[res, Simplify[f / Times@@res]]
]


Examples:

separate[Exp[-x^2-y^2-z^2](x+x^2)Sin[y]]
separate[Expand[(x (y+1))^5]]


{1, E^-x^2 x (1 + x), E^-y^2 Sin[y], E^-z^2}

{1, x^5, (1 + y)^5}

For your application, I think you can use the following function:

decompose[expr_] := Module[{terms},
terms = Replace[Expand[expr],
{
a_Plus :> List @@ a,
a_ :> {a}
}
];
separate /@ terms
]


decompose[Exp[-x^2-y^2]]
`