# How can I separate a separable function

I have a separable function $f[x,y]$, and I would like to find two functions $g[x]$ and $h[y]$ with

$f[x,y]=g[x] h[y]$

where $g[x]$ doesn't depend on $y$ and $h[y]$ doesn't depend on $x$. Ideally, $g$ and $h$ should have the same magnitude, to prevent overflows/underflows. I have a hackish approach that works, but involves a lot of manual labor.

Background: $f[x,y]$ is a filter kernel I want to apply to an image, and using two separate 1d-filters is much more efficient.

My first approach was to start with $g[x]=f[x,0]$. But that doesn't work for e.g. $f[x,y]=\frac{e^{-\frac{x^2+y^2}{2 \sigma ^2}} x y}{2 \pi \sigma ^6}$

Currently, I have a function that "removes" $x$ or $y$ from $f[x,y]$ using pattern matching:

removeSymbol[f_, s_] := f //. {s^_ + a_ -> a, s^_.*a_ -> a}


but that means I have to manually adjust this pattern for different f's.

Is there a more elegant way to do this? $f[x,y]$ is usually a derivative of a gaussian, e.g.

gaussian[x_,y_] := 1/(2 π σ^2) Exp[-((x^2 + y^2)/(2 σ^2))]
f[x_,y_] := D[gaussian[x,y], x, y]

-
Since you "have to manually adjust this pattern for different f's" it would be helpful to see some different f functions. –  Mr.Wizard Jul 10 '12 at 14:51
Did you check whether the built-in GaussianFilter[image,r,{nx,ny}] would work for you? It might be sufficiently fast that you don't need to separate your kernel. –  Sjoerd C. de Vries Jul 10 '12 at 14:55
@Mr.Wizard: Ideally it would work for any function, but all the functions I've used so far were various derivatives of gaussian, D[gaussian[x,y], x], D[gaussian[x,y], x,x,y] and so on –  nikie Jul 10 '12 at 14:58
@SjoerdC.deVries: I'm using GaussianFilter when I'm prototyping an algorithm in Mathematica. But I have to build the filter kernels manually when I write the final version in C. –  nikie Jul 10 '12 at 14:59
OK, I see. Makes sense. –  Sjoerd C. de Vries Jul 10 '12 at 15:03

I would take a logarithmic derivative with respect to one variable - it should be then independent of the other one, then integrate it back over the first variable and exponentiate. The second function is found by plain division. Here is the code:

ClearAll[getGX];
getGX[expr_, xvar_, yvar_] :=
With[{dlogg = D[Log[expr], xvar] // FullSimplify},
Exp[Integrate[dlogg, xvar]] /; FreeQ[dlogg, yvar]];

Clear[getHY];
getHY[expr_, xvar_, yvar_] := FullSimplify[(#/getGX[#, xvar, yvar]) &[expr]]


A test function:

ftest[x_, y_] := (x^2 + 1)*y^3 *Exp[-x - y]


Now,

getGX[ftest[x,y],x,y]

(* E^-x (1+x^2)  *)

getHY[ftest[x,y],x,y]

(* E^-y y^3 *)


The integration constant ambiguity translates into an ambiguity of how you split the function, since this operation is only defined up to a multiplicative constant factor by which you can multiply one function, and divide the other one.

-
Dumb question: How can this work if expr can be negative and Log is only defined for positive values? Obviously it works, but I can't see why. –  nikie Jul 10 '12 at 15:15
@nikie Log is defined in the complex plane, not just real line. So, you can think of it as if I have analytically continued our function. But, we don't have to be so fancy - just replace the logarithmic derivative by D[expr,xvar]/expr, and you don't have to think of this issue then. –  Leonid Shifrin Jul 10 '12 at 15:18
Of course, I didn't think of complex numbers. But it even works if expr is zero. And neither Log[expr] nor D[expr,xvar]/expr is defined for expr=0. And still, it gives the right result. –  nikie Jul 10 '12 at 18:19

What I'd do:

gaussian[x_, y_] := 1/(2 π σ^2) Exp[-((x^2 + y^2)/(2 σ^2))];
f[x_, y_] = D[gaussian[x, y], x, y]

Exp[Select[Expand[PowerExpand[Log[Together[f[x, y]]]]], #]] & /@
{FreeQ[#, x | y] &, ! FreeQ[#, x] &, ! FreeQ[#, y] &}
{1/(2 Pi σ^6), E^(-(x^2/(2 σ^2))) x, E^(-(y^2/(2 σ^2))) y}


The snippet separates out the constant factor, the factors with x, and the factors with y.

More examples:

f = -Pi Cos[x]^2 Sin[y]^3/E;
Exp[Select[Expand[PowerExpand[Log[Together[f]]]], #]] & /@
{FreeQ[#, x | y] &, ! FreeQ[#, x] &, ! FreeQ[#, y] &}
{-Pi/E, Cos[x]^2, Sin[y]^3}


We see that negative constant factors are reproduced.

f = w[x] z[y];
Exp[Select[Expand[PowerExpand[Log[Together[f]]]], #]] & /@
{FreeQ[#, x | y] &, ! FreeQ[#, x] &, ! FreeQ[#, y] &}
{1, w[x], z[y]}


The implicit constant factor of 1 is detected.

-
Great! I had to add Together[f[x,y]] to make it work for D[gauss, x, x], but this looks very promising. –  nikie Jul 10 '12 at 15:00
Ah, yes. I've added it to my snippet. –  Ｊ. Ｍ. Jul 10 '12 at 15:21
(I wish I knew why my answer seems to be less optimal than Leonid's for readers, given the voting patterns...) –  Ｊ. Ｍ. Jul 10 '12 at 16:02
I think that my answer is based more on Mathematics than Mathematica. Which is more optimal is of course subjective. –  Leonid Shifrin Jul 10 '12 at 16:22
I upvoted both. And I really like this solution, because selecting factors out of a list seems far more intuitive than all those derivatives and integrals. But I do see advantages in @Leonid's answer: For example, it works for sums, even if they can't be factored by Together (Example: f[x_, y_] := Expand[(x*(y + 1))^5]). Taking the log to get a list of the factors only works if the original expression is a product. Taking the logarithmic derivative works for sums, too. –  nikie Jul 10 '12 at 17:52