3

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$

Fraccalo
  • 6,057
  • 13
  • 28
  • 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 08:16

1 Answers1

5

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->Reduce`FreeVariables[f],
        Except[_List]->{vars}
        }
    ];
    res = Exp @ MapThread[
        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
]

For your examples:

decompose[Exp[-x^2-y^2]]
decompose[Exp[-x^2-y^2](x-y)]

{{1, E^-x^2, E^-y^2}}

{{1, E^-x^2 x, E^-y^2}, {-1, E^-x^2, E^-y^2 y}}

Note that the constant terms are different from what you expected.

Carl Woll
  • 130,679
  • 6
  • 243
  • 355