7
$\begingroup$

I have a large symbolic expression with six variables a,b,c,d,e,f per term, roughly like this:

expr=a^2 b c d^2 + 1/2 a b (b + I d) e^3 - 1/8 I a (c + 1/2 I (c + I d))^2 (I b + 1/2 (c + I d)) e f

I want to test whether expr contains terms with c and d and e and f. In the example above, there is one term that fulfills the requirement: (-(1/8)-I/16) a b c d e f.

I want a function g[x], which returns a new expression, only containing terms with c*d*e*f.

I can solve it using Expand and symbolic replacements, but this is very slow, and I was hoping for a much faster way. Does anybody have a suggestion for a faster implementation?

Here is the example code - which first creates the expression (such that we can compare our solutions), then defines g[x], runs it and prints the result and time:

(* Same seed for fair comparison *)
SeedRandom[1];

(* Creating one large expression that has to be analysed *)
func = 0;
For[ii = 1, ii <= 100, ii++,
  rndFull = Product[RandomChoice[{a, b, c, d, e, f}], {i, 6}];
  For[jj = 1, jj <= 60, jj++,
   rndVar0 = RandomChoice[{a, b, c, d, e, f}];
   rndVar1 = RandomChoice[{a, b, c, d, e, f}];
   rndVar2 = RandomChoice[{a, b, c, d, e, f}];
   rndFull = rndFull /. {rndVar0 -> 1/2*(I*rndVar1 + rndVar2)};
   ];
  func += Simplify[Exp[(I*\[Pi])/2*RandomInteger[4]]]*rndFull;
  ];

(* g[x] works correctly but is very slow *)
g[expr_] := (Return[Expand[ZERO*expr] /. {ZERO*c*d*e*f -> c*d*e*f} /. {ZERO -> 0}])


CurrTime = AbsoluteTime[];
Print[g[func]]
(* (-(3066695705460323281748465708319/340282366920938463463374607431768211456)-(592396087826201433092643851215 I)/170141183460469231731687303715884105728) a^2 c d e f-(28061855884788386235225/1267650600228229401496703205376-(1188923985339435970275 I)/158456325028528675187087900672) a b c d e f-(23306256125181506635725/2535301200456458802993406410752-(47084620130455206162825 I)/5070602400912917605986812821504) b^2 c d e f *)
Print[AbsoluteTime[] - CurrTime] (* 4.5930638 sec *)

Does anybody have a suggestion for a faster implementation?

$\endgroup$
  • $\begingroup$ So in your example expr, a b c d e f counts but a d c d e f doesn't? This will be very dependent on how the terms in the expression are arranged. Would you want it to generalise to, for example, picking out terms involving b d e f? $\endgroup$ – aardvark2012 Oct 31 '17 at 2:37
  • $\begingroup$ If you're just interested in the literal occurrence of c d e f as a string in the expanded form of expr, and you're not worried about the order depending on Mathematica's canonical ordering, you could use Pick[#, StringContainsQ["c d e f"] /@ ToString /@ #] &@(List @@ Expand@expr). $\endgroup$ – aardvark2012 Oct 31 '17 at 2:48
  • 1
    $\begingroup$ @aardvark2012 No, the order does not matter when doing /. Mathematica will match a b or b a for the pattern a b. Compare a b c/.a b->x and a b c/.b a->x they give same answer which is c x. Your string solution will not work, since only a b will be matched and not b a as in the case of general pattern $\endgroup$ – Nasser Oct 31 '17 at 2:53
  • 2
    $\begingroup$ As you note, all the time seems to be being spent on Expand (Simplify is even worse). Coefficient[func, c d e f] c d e f returns the same as g and is almost order of magnitude faster when they're both passed already Expanded (or Simplifyed) expressions (and marginally faster when the argument isn't pre-Expanded, but there's not a lot to choose between them). $\endgroup$ – aardvark2012 Oct 31 '17 at 4:05
  • $\begingroup$ @aardvark2012 thanks for your solution; unfortunatly it also uses Expand, thus its as fast/slow as my original solution. Also it gives a slightly different result, as it also returns expressions like b c d e f^2 which I do not want. $\endgroup$ – NicoDean Oct 31 '17 at 12:06
3
+50
$\begingroup$

Here is another faster solution:

prune[X_][expr_] := Select[expr, MemberQ[#, X, \[Infinity]] &]
h[expr_] := Plus@@Cases[Expand@prune[c]@prune[d]@prune[e]@prune[f]@expr, c d e f X_] 
r1 = g@func; //Timing
r2 = h@func; //Timing
r1 == r2

{12.0938, Null}

{1.39063, Null}

True

We see a speed up of almost a factor 10.

$\endgroup$
  • $\begingroup$ An operator version of prune would be prune[x_] = Select[Not@*FreeQ[x]]. $\endgroup$ – Carl Woll Nov 8 '17 at 23:09
  • $\begingroup$ Yes, and I'm sure the pruning could be made more efficient. However, the step after the prune is the most time consuming. $\endgroup$ – mmeent Nov 9 '17 at 10:09
  • $\begingroup$ Very nice solution, and indeed in my test more than 5 times faster. Thanks +50 $\endgroup$ – NicoDean Nov 10 '17 at 12:26
4
$\begingroup$

The following is slightly faster:

r1 = With[{dd = D[func, c, d, e, f]},
    Expand[c d e f Block[{c = 0, d = 0, e = 0, f = 0}, dd]]
]; //AbsoluteTiming

r2 = g[func]; //AbsoluteTiming

r1 === r2

{2.8915, Null}

{3.49264, Null}

True

$\endgroup$
  • $\begingroup$ Wow that is very clever, using a derivative :-) Didn't think about that. Thanks for the solution. $\endgroup$ – NicoDean Oct 31 '17 at 15:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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