# 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. Commented Jul 10, 2012 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. Commented Jul 10, 2012 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 Commented Jul 10, 2012 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. Commented Jul 10, 2012 at 14:59
• OK, I see. Makes sense. Commented Jul 10, 2012 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. Commented Jul 10, 2012 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. Commented Jul 10, 2012 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. Commented Jul 10, 2012 at 18:19
• when the expression can't be separated, the call does not even go through. Why is that? Should not these function return {} or such in this case? For example getGX[1/(x +y),x,y] Mathematica just echo the input back to the screen. This makes it little hard to check if the expression was separated or not? Commented Dec 23, 2017 at 7:04
• And here is another strange result for something that can't be separated getHY[3 y + Sin[x], x, y] now Mathematica echos this: (3 y+Sin[x])/getGX[3 y+Sin[x],x,y]. How do you suggest the correct way one should check the result for the cases when the expression can't be separated? Commented Dec 23, 2017 at 7:14

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. Commented Jul 10, 2012 at 15:00
• Ah, yes. I've added it to my snippet. Commented Jul 10, 2012 at 15:21
• (I wish I knew why my answer seems to be less optimal than Leonid's for readers, given the voting patterns...) Commented Jul 10, 2012 at 16:02
• What I mean is that my algorithm can, for example, be also applied numerically, for function which is defined numerically, where the notions of Expand, Together etc make no sense. Actually, Expand and Together, as well as other operations such as PowerExpand, FreeQ, don't have a well-defined mathematical meaning, and are only defined in Mathematica. Again, this does not answer which solution is more optimal, but it seems that mine is more universal in the sense that its semantics largely does not depend on a particular tool used to implement it (Mathematica here). Commented Jul 10, 2012 at 16:41
• 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. Commented Jul 10, 2012 at 17:52

While the two existing answers are good, I've always approached this sort of problem with FactorList and there should be an answer showing how to use this old Mathematica function in this way. Much time having passed since this question was asked, I'll use current tools to process the list of factors.

I decided to separate the task into two parts, which makes the solution a little more complicated than necessary. OTOH, I've sometimes wanted to know the factors, including ones that cannot put in the form $$g(x)h(y)$$, separated by type.

ClearAll[sepFac, sepParts, dependsOn];
sepParts[expr_, vars_List] := With[{
factors = GroupBy[FactorList[expr],
Function[fac,
Pick[dependsOn @@ vars, InternalDependsOnQ[fac, #] & /@ vars]]
],
apply = Function[{h, lev}, Apply[h, #, lev] &]},
apply[Times, {1}]@ apply[Power, {2}]@ factors
];
sepFac::nonsep = "The factor  could not be separated.";
sepFac[expr_, vars_List] := With[{parts = sepParts[expr, vars]},
Message[sepFac::nonsep, #] & /@
Values@KeySelect[parts, Length[#] >= 2 &];
(MapAt[dependsOn[] # &, (* multiply first by constant factor *)
dependsOn /@ vars,
1] /. parts
) /; Max[Length /@ Keys[parts]] <= 1
];


Examples:

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

(*
{(E^(-(x^2/(2 σ^2))) x)/(2 π σ^6),
E^(-(y^2/(2 σ^2))) y}
*)

sepFac[(x + y) x y, {x, y}]

sepFac::nonsep :  The factor x+y could not be separated.

(*  sepFac[x y (x + y), {x, y}]  *)

sepParts[(x + y) x y, {x, y}]

(*
<|dependsOn[] -> 1, dependsOn[x] -> x, dependsOn[y] -> y,
dependsOn[x, y] -> x + y|>
*)
`