How can I separate a separable function

Question

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 — Niki Estner, 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. — Niki Estner, Jul 10 '12 at 14:59
Another version: https://mathematica.stackexchange.com/a/224888/69875 — SnzFor16Min, Jun 29 '20 at 05:29

Leonid Shifrin · Accepted Answer · 2012-07-10T15:05:53.253

24

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.

edited Jul 10 '12 at 15:05

answered Jul 10 '12 at 14:59

Leonid Shifrin

114,335
15
329
420

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. – Niki Estner 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. – Niki Estner Jul 10 '12 at 18:19
1

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? – Nasser Dec 23 '17 at 07: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? – Nasser Dec 23 '17 at 07:14
@Nasser "when the expression can't be separated, the call does not even go through. Why is that?" - because of the way I have written the function, particularly the condition FreeQ[dlogg, yvar]. "How do you suggest the correct way one should check the result for the cases when the expression can't be separated?" - I wouldn't know off hand, and alas have no time right now to look closer into it. – Leonid Shifrin Dec 23 '17 at 12:08
fyi, I have a case where it does not seem to produce correct splitting. Here it is: expr = (x^2*y + y)*y^3*Exp[-x - y + 1]*3^(-x - y)*Sqrt[x^2*y - y]; now getHY[expr,x,y] gives expression which contains both x and y in it: (3^-y E^(1 - y) y^4 Sqrt[(-1 + x^2) y])/Sqrt[1 - x^2]. The splitting of this expression should be: xPart = (x^2 + 1) Sqrt[x^2 - 1] Exp[-x] 3^(-x) and yPart = y^(9/2) Exp[-y + 1] 3^(-y) May be I should open a new question on this. Not sure now. – Nasser Dec 24 '17 at 09:10
@Nasser The dependence on x in your above result is merely a complex phase. If both multipliers under the square root are non-negative, you can do FullSimplify[expr/getGX[expr, x, y], Assumptions -> y > 0 && x >= 1 ], which results in an expression that depends on y only. Such issues are inevitable given that variables are generally considered as complex numbers in Mathematica, and for some expressions one would need extra conditions to ensure they are real and can be factored. Perhaps I should've mentioned that in the post. – Leonid Shifrin Dec 24 '17 at 13:22
You are right, Here is a small example of where assumptions are needed. f=Sqrt[x y], splitting gives Sqrt[x] and Sqrt[x y]/Sqrt[x] but the second one can be simplified to Sqrt[y] only under assumption that x>0. I found paper that takes about method, very similar to what you show but was not too clear for me to follow all the details. You might be interested in it. Here is the link https://arxiv.org/ftp/math/papers/0210/0210167.pdf – Nasser Jun 12 '18 at 09:33
@Nasser Thanks, I will have a look when I get a spare moment. – Leonid Shifrin Jun 12 '18 at 10:15

J. M.'s missing motivation · Answer 2 · 2012-07-10T15:50:54.323

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.

edited Jul 10 '12 at 15:50

answered Jul 10 '12 at 14:53

J. M.'s missing motivation

124,525
11
401
574

1

Great! I had to add Together[f[x,y]] to make it work for D[gauss, x, x], but this looks very promising. – Niki Estner Jul 10 '12 at 15:00
Ah, yes. I've added it to my snippet. – J. M.'s missing motivation 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.'s missing motivation 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
@Leonid: I don't know, it's just the basic idea of take a logarithm, pick out terms, and exponentiate... I didn't even have to differentiate/integrate here, and the *Expand[] and Together[] functions just expedite the separation for Select[]. – J. M.'s missing motivation Jul 10 '12 at 16:25
1

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). – Leonid Shifrin Jul 10 '12 at 16:41
3

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. – Niki Estner Jul 10 '12 at 17:52

Michael E2 · Answer 3 · 2021-02-17T19:12:13.333

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, Internal`DependsOnQ[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|>
*)

How can I separate a separable function

3 Answers3

Linked

Related