Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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).

share|improve this question
6  
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. –  rcollyer May 13 '13 at 16:33
1  
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. –  Oleksandr R. May 14 '13 at 11:40
1  
@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. –  Rex Kerr May 14 '13 at 15:27

3 Answers 3

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}]] *)
share|improve this answer
    
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? –  Rex Kerr May 14 '13 at 0:15
2  
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. –  Oleksandr R. May 14 '13 at 4:00
    
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. –  Daniel Lichtblau May 14 '13 at 13:42
    
@DanielLichtblau Do you happen to know why OptimizeExpression fails to extract the common expression in the second example in this post? –  Silvia Apr 3 at 3:07
    
@Silvia No, I do not know what is going on there. Will send that example to relevant parties. –  Daniel Lichtblau Apr 3 at 15:05

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

share|improve this answer
    
This is interesting at e.g. level 1! Thanks! –  Rex Kerr May 14 '13 at 0:10
    
Yes. I made us of your observation in the latest version. –  David Carraher May 14 '13 at 2:35

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]  
share|improve this answer
1  
+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). –  Rex Kerr May 14 '13 at 0:09

Your Answer

 
discard

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.