32
$\begingroup$

I am working with some unpleasantly tedious polynomials, which need to be manipulated in various ways (integrate with respect to some variable, differentiate with respect to another). Since these will eventually be part of numeric routines, I'd like to pull out common subexpressions that I can reuse, e.g.

(G u^2 (6 p (2 h + p) - 8 (h + p) u + 3 u^2))/(12 h^2)
(G (3 h + 3 p - 2 u) u^2)/(3 h^2)

should at least note that both have a common factor of

(G u^2)/(3 h^2)

Is there a convenient way to instruct Mathematica to look for this sort of computational-expense-reduction in pairs of expressions? Ideally it would notice even in the case where it's not just a factor multiplied by both, e.g. if I added + 1 to the second equation, I'd still like it to find that common subexpression.

Just for clarification, this is what I do by hand:

e1 = (G u^2)/(12 h^2) * (6 p (2 h + p) - 8 (h + p) u + 3 u^2))
e2 = (G u^2)/( 3 h^2) * (3 h + 3 p - 2 u)

A = (G u^2)/(12 h^2)
e1 = A * (6 p (2 h + p) - 8 (h + p) u + 3 u^2))
e2 = A * (12 (h + p) - 6 u)

B = u^2
C = h + p
A = (G B)/(12 h^2)
e1 = A * (6 p (2 h + p) - 8 C u + 3 B)
e2 = A * (12 C - 6 u)

to go from 6 additions, 20 multiplications, 2 divisions, and 2 assignments, to 5 additions, 14 multiplications, 1 division, and 5 assignments (the extra three of which are temporary so can be in registers and cost essentially nothing).

Improve this question
$\endgroup$
3
  • 6
    $\begingroup$ Leonid put together a common sub-expression eliminator for another question. I have not tested it on polynomials, hence this is a comment, not an answer, but it may be useful, here. $\endgroup$
    – rcollyer
    May 13, 2013 at 16:33
  • 1
    $\begingroup$ Could you clarify how this result will be used, i.e. whether the numerical routine will be implemented in Mathematica or another language? In the former case I wouldn't worry so much about strict optimality because if your goal is to get peak hardware performance then you'll probably end up disappointed. In the latter, my instinct would be to use Simplify et al. and then let your compiler do the CSE and other transformations: it's better placed to know what optimizations are possible at that stage and by doing it prematurely you may disrupt this process. $\endgroup$
    – Oleksandr R.
    May 14, 2013 at 11:40
  • 1
    $\begingroup$ @OleksandrR. - I will probably port it to another language since I doubt that Mathematica will be able to do it fast enough--and it will probably be Java (due to other concerns) where the optimization is deferred to the JIT compiler, and it doesn't have time to do CSE over a very broad scope (the example above is one of the smallest cases), plus I will have to do a fair bit of work by hand, so CSE would also save me time by having less work to check/convert. $\endgroup$
    – Rex Kerr
    May 14, 2013 at 15:27

3 Answers 3

Reset to default
31
$\begingroup$

The engine behind this inside Compile is a well-hidden function called OptimizeExpression. it has two levels, 1 and 2. Setting to 2 makes it work harder to find CSEs.

e1 = (G u^2 (6 p (2 h + p) - 8 (h + p) u + 3 u^2))/(12 h^2);
e2 = (G (3 h + 3 p - 2 u) u^2)/(3 h^2);

Experimental`OptimizeExpression[{e1, e2}, 
 OptimizationLevel -> 2]

(* Out[40]= Experimental`OptimizedExpression[
 Block[{Compile`$7, Compile`$9, Compile`$10}, Compile`$7 = h^2; 
  Compile`$9 = 1/Compile`$7; 
      Compile`$10 = 
   u^2; {1/12 G Compile`$9 Compile`$10 (6 p (2 h + p) - 8 (h + p) u + 
          3 Compile`$10), 
   1/3 G Compile`$9 (3 h + 3 p - 2 u) Compile`$10}]] *)
Improve this answer
$\endgroup$
6
  • $\begingroup$ This is close, but why doesn't it understand that G should be part of Compile$9`? Surely I don't want to do an extra multiplication for no reason. Also, can one take this apart to get access to the block and individual expressions? $\endgroup$
    – Rex Kerr
    May 14, 2013 at 0:15
  • 2
    $\begingroup$ Probably worth noting that, with "OptimizationLevel" -> 2, it does CSE completely without regard to context. Thus, it will more than likely produce a non-working result if there are subexpressions whose side effects are relied upon. $\endgroup$
    – Oleksandr R.
    May 14, 2013 at 4:00
  • $\begingroup$ Not sure why it misses multiplying G*Compile`$9*Compile`$10 as a CSE. Overall it would save a multiplication if I am counting correctly. $\endgroup$
    – Daniel Lichtblau
    May 14, 2013 at 13:42
  • 1
    $\begingroup$ @DanielLichtblau Do you happen to know why OptimizeExpression fails to extract the common expression in the second example in this post? $\endgroup$
    – Silvia
    Apr 3, 2014 at 3:07
  • $\begingroup$ @Silvia No, I do not know what is going on there. Will send that example to relevant parties. $\endgroup$
    – Daniel Lichtblau
    Apr 3, 2014 at 15:05
11
$\begingroup$

Level may provide a means of getting started. Level breaks down your polynomials into their constituent parts, with various "levels" of complexity. As @Rex Kerr noted in a comment, Level 1 happens to be interesting in the present example. 

a=(G u^2 (6 p (2 h+p)-8 (h+p) u+3 u^2))/(12 h^2);
b=(G (3 h+3 p-2 u) u^2)/(3 h^2);
Intersection[Level[a,1],Level[b,1]]

level

Improve this answer
$\endgroup$
2
  • $\begingroup$ This is interesting at e.g. level 1! Thanks! $\endgroup$
    – Rex Kerr
    May 14, 2013 at 0:10
  • $\begingroup$ Yes. I made us of your observation in the latest version. $\endgroup$
    – DavidC
    May 14, 2013 at 2:35
6
$\begingroup$

Since you are going to work with numeric functions, Compile will optimize your functions along those lines. If you define :

f = Compile[{{g, _Real}, {u, _Real}, {h, _Real}, {p, _Real}},
 {(g u^2 (6 p (2 h + p) - 8 (h + p) u + 3 u^2))/(12 h^2), 
  (g (3 h + 3 p - 2 u) u^2)/(3 h^2)}
    ]

you can already see some of it in the output. To get an even more in depth look try : 

Needs["CompiledFunctionTools`"]
CompilePrint[f]
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ +1 since this is better than nothing, but it has two defects: it is no longer symbolic, so you lose the opportunity to do any further manipulation; and it is not very good at identifying those commonalities useful for numeric computation (for example, it divides by h^2 twice instead of once, even though division is much more expensive than multiplication, and it multiplies by g twice even though there's no reason to). $\endgroup$
    – Rex Kerr
    May 14, 2013 at 0:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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