Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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]
share|improve this question
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
up vote 12 down vote accepted

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:

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

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

A test function:

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



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


(* 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.

share|improve this answer
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.

share|improve this answer
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. – J. M. 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...) – J. M. 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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.