# Find subexpression to minimize leafcount after replacment with temporary variable

Given an expression expr, is there an automated way to find a subexpression subexpr such that replacing subexpr with a temporary variable v minimizes LeafCount[expr//.subexpr->v]+LeafCount[subexpr] over all possible subexpr?

To provide a concrete example to work with, here is a cleaned up small fragment of an expression that someone asked an unrelated question about recently:

(1/Sqrt[2](Sqrt[((Sqrt[((-x3 - (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) + Tan[a]) - Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] + Tan[b])^2)] + Sqrt[((-x3 - Cos[a] (m-h Cos[a])/(2 (Tan[a]-Tan[b])))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] -
(Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)] + Sqrt[((-(Cos[a] (m-h Cos[a]))/(2 (Tan[a]-Tan[b])) + (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) + Tan[a]) - Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-Tan[b] -
(Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)])(-Sqrt[((-x3 - (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) + Tan[a]) - Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] +
Tan[b])^2)] + 1/2 (Sqrt[((-x3 - (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (- m Sec[a]/(x2-x3) + Tan[a])-Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-1/2 Cos[a] (m-h Cos[a]) -x3 Tan[a] + Tan[b])^2)] + Sqrt[((-x3 - (Cos[a] (m-h Cos[a]))/
(2 (Tan[a]-Tan[b])))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] -(Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)] +Sqrt[((-(Cos[a] (m-h Cos[a]))/(2 (Tan[a]-Tan[b])) +(1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) +
Tan[a])-Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-Tan[b] -(Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)])))]))


The original problem had a LeafCount of almost 6000. Fumbling with manual substitutions I was finally able to make the structure clear and that happened to enable Simplify to reduce the total LeafCount to about 250 in a few seconds.

Is there a simple method of automating this? It would be even better if it could recognize and handle -subexpr and 1/subexpr as containing subexpr, but that may be too much to ask.

Note: I'm trying to actually see the result, not just have Mathematica do subexpression sharing behind the scenes where I cannot see what happened.

-
If your subexpression == expression the count comes down to 1. So perhaps other restrictions are in place ... and ought to be specified – Dr. belisarius Jul 27 '12 at 20:41
@belisarius Well, you could say that the leafcount to minimize is the leafcount of the expression along with any expression needed to define the used subexpressions, then the result makes slightly more sense. – jVincent Jul 27 '12 at 21:11

The compiler does a pretty good job of finding and eliminating common subexpressions. We can leverage that functionality.

Starting from the provided expression...

$expr = (1/Sqrt[2](Sqrt[((Sqrt[((-x3 - (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) + Tan[a]) - Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] + Tan[b])^2)] + Sqrt[((-x3 - Cos[a] (m-h Cos[a])/(2 (Tan[a]-Tan[b])))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] - (Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)] + Sqrt[((-(Cos[a] (m-h Cos[a]))/(2 (Tan[a]-Tan[b])) + (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) + Tan[a]) - Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-Tan[b] - (Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)])(-Sqrt[((-x3 - (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) + Tan[a]) - Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] + Tan[b])^2)] + 1/2 (Sqrt[((-x3 - (1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (- m Sec[a]/(x2-x3) + Tan[a])-Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-1/2 Cos[a] (m-h Cos[a]) -x3 Tan[a] + Tan[b])^2)] + Sqrt[((-x3 - (Cos[a] (m-h Cos[a]))/ (2 (Tan[a]-Tan[b])))^2 + (-1/2 Cos[a] (m-h Cos[a]) - x3 Tan[a] -(Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)] +Sqrt[((-(Cos[a] (m-h Cos[a]))/(2 (Tan[a]-Tan[b])) +(1/2 Cos[a] (m-h Cos[a]) + x3 Tan[a] - x3 (-m Sec[a]/(x2-x3) + Tan[a])-Tan[b])/(-m Sec[a]/(x2-x3) + Tan[a]))^2 + (-Tan[b] -(Cos[a] (m-h Cos[a]) Tan[b])/(2 (Tan[a]-Tan[b])))^2)])))]));  We extract the variable names: $variables = Cases[$expr, _Symbol, Infinity] // Union (* {a, b, h, m, x2, x3} *)  Now we compile the expression. The result is an object which can be destructured. We are interested in the component which represents the compiled Function expression: $function =
Cases[
Compile[##]&[$variables,$expr]
, x_Function :> x
] // First

(*
Function[{a,b,h,m,x2,x3},
Block[{Compile$398,Compile$406,[...snip...],Compile$441}, Compile$398=-3;
Compile$406=Cos[a]; Compile$403=Tan[a];
Compile$399=x2+Compile$398;
[...snip...]
Compile$441=Compile$423+Compile$433+Compile$440;
Sqrt[Compile$441 (-Compile$423+Compile$441/2)]/Sqrt[2]]] *)  In this result, we can see that the compiler has found all of the common subexpressions and assigned them to variables with generated names. If the output is intended for machine consumption, we are done. All of the symbols have been safely localized. But as humans, we might find all of the generated symbol names distracting. We can apply a few more transformations in the interest of readability. We will change all of the variables of the form Compile$nnn to p[n], where p is carefully chosen as a symbol that does not appear in the original expression:

FreeQ[$expr, p] (* True *)  The following transformations do the trick: $cse =
Cases[$function, s_Symbol /; "Compile" === Context[s], Infinity] // DeleteDuplicates // MapIndexed[#1 -> p@@#2 &, #]& //$function /.
# /.
HoldPattern[Function[_, Block[_, b_]]] :> Hold@\[FormalM][{p}, b] /.
\[FormalM] -> Module &

(*
Hold[Module[{p},
p[1] = -x3;p[2] = Cos[a]; p[3] = Tan[a];
p[4] = x2 + p[1]; p[5] = 1/p[4]; p[6] = Sec[a]; p[7] = -m p[5] p[6];
[...snip...]
p[44] = p[30] + p[36] + p[43];
Sqrt[p[44] (-p[30] + p[44]/2)]/Sqrt[2]]]
*)


The final expression expands to the original expression when evaluated:

$expr === ReleaseHold[$cse]

(* True *)


If desired, another transformation will give us a base expression and replacement rules:

$rep =$cse /. Hold[Module[_, CompoundExpression[s___, f_]]] :> Hold[f /. {s}] /. Set -> Rule

(*
Hold[Sqrt[p[44](-p[30]+p[44]/2)]/Sqrt[2] /. {p[1]->-x3, p[2]->Cos[a], [...snip...]}]
*)


Inspection of this result reveals that the compiler sometimes introduces temporary variables even for expressions that are only used once. Such variables appear exactly twice in the expression -- once to define it and once for its use:

$unnecessary = Cases[Cases[$rep, p[_], Infinity] // Tally, {_, 2}][[All, 1]]

(*
{p[4],p[5],p[6],p[7],p[9],p[15],p[19],p[20],p[21],p[24],p[25],p[26],p[27],p[28],
p[29],p[31],p[32],p[33],p[34],p[35],p[36],p[37],p[38],p[39],p[40],p[41],p[42],p[43]}
*)


We can further simplify the compiler's expression by removing these unnecessary variables:

$cse2 = Verbatim[Rule][Alternatives @@$unnecessary, _] //
DeleteCases[$rep, #, Infinity] //. Cases[$rep, #, Infinity] &


Once again, this expression is equivalent to the original when evaluated:

ReleaseHold[$cse2] ===$expr

(* True *)


Just For Fun

Sometimes a layered graph plot can come in handy when trying to visualize the structure of complex expressions:

edge[x:_[___]] := (Scan[If[!AtomQ[#], Sow[x -> #]]& , x]; x)
edge[x_] := x
Reap[edge //@ $expr][[2, 1]] // DeleteDuplicates // LayeredGraphPlot[#, Top, DirectedEdges -> False]&  The expression at each node is revealed as a tooltip when the mouse hovers over a vertex. - you could also use ExperimentalOptimizeExpression[$expr] – acl Jul 28 '12 at 13:59
+2 Beachcombing the compiled code is a great idea! – Dr. belisarius Jul 28 '12 at 18:12
Great idea and answer! – Thies Heidecke Jul 28 '12 at 19:10
Is there a limit to how large an expression that Mathematica will do common subexpression elimination on? I have an expression with LeafCount of 98701 with 6 variables and {List, Log, Plus, Power, Times}  as heads (from Cases[$expr, x_[__] -> x, {0, Infinity}, Heads -> True] // Union ) – Åsmund Hj Feb 26 at 9:21 @ÅsmundHj Sorry, but I have no knowledge concerning the upper limits of this machinery. – WReach Feb 26 at 15:05 Assuming that we assign your test expression to a variable expr, and using my common subexpression eliminator gives immediately without any tweaking: res = csub[Hold[Evaluate[expr]], {}, 10] (* Hold[ Let[{$3:=-((m Sec[a])/(x2-x3))+Tan[a],
$4:=1/2 Cos[a] (m-h Cos[a])+x3 Tan[a]-x3$3-Tan[b],
$5:=-(1/2) Cos[a] (m-h Cos[a]),$6:=(Tan[b] $5)/(Tan[a]-Tan[b]),$7:=$5/(Tan[a]-Tan[b]),$8:=(-x3-$4/$3)^2+($5-x3 Tan[a]+Tan[b])^2,$9:=Sqrt[$8]+Sqrt[(-x3+$7)^2+($5-x3 Tan[a]+$6)^2]+
Sqrt[($7+$4/$3)^2+(-Tan[b]+$6)^2]
},
Hold[Sqrt[$9 (-Sqrt[$8]+$9/2)]/Sqrt[2]] ] ] *)  Please see the mentioned post for the definitions of csub and Let. The LeafCount of this is 195, and that includes 2 extra Hold-s and one Let, all of which are strictly not parts of the expression. Note that calling ReleaseHold@ReleaseHold@res  will bring you back to your expression. - Hi Leonid! Just out of laziness ... How do I invoke it to treat something like expr = (Sin[(x + y + z)]/(x + y + z))^(x + y + z)? – Dr. belisarius Jul 27 '12 at 22:20 Hi @belisarius! You will have to modify code for csub a bit, and tweak the limitCount parameter. In the code for csub, I put a restriction Depth[Unevaluated[x]] > 2. Change this to Depth[Unevaluated[x]] > 1, and set limitCount to 2 (say): csub[Hold[Evaluate[expr]], {}, 2]. I get then: Hold[Let[{$25 := x + y + z}, Hold[(Sin[$25]/$25)^\$25]]]. I should probably make a minimal depth as a parameter. I could also add a LeafCount minimizer which would automatically try different values for parameters and optimize minimal depth and limitCount. – Leonid Shifrin Jul 27 '12 at 22:30
Ok. Got it working. Thanks! – Dr. belisarius Jul 27 '12 at 22:43

I don't think you can have such things fully automated, but as for tips to help you identify subexpressions, I think this is helpfull:

 PatternTable[
expression_] := {(#[[1]] - 1) (#[[2]] - 1), #[[1]], #[[2]], #[[3]]} & /@
({Count[exp, #, \[Infinity]], LeafCount[#], #} & /@
Cases[exp, _, \[Infinity]] // Union) // Sort // Reverse


I know it's messy, but it's just suggestive. It returns a list of lists containing

 {score number(to be explained), number of appearances, leafcount, and sub-expression}


Then you can look through it and see which expressions seem worth it for you to define separately. The score number is just an estimate of how much you can reduce the leaf count by making the replacement, it goes by the logic, that if you have n appearances of a subexpression with l leafs, then you can at most replace n-1, since you still need it once in your definition, and you still leave 1 leaf everywhere you replace it, so the score is (n-1)(l-1). It serves nicely to sort the expressions into a suggestive order.

You can then skim through the entire list, or sort it depeding on what you deam suitable criteria for a replacement being worth it, for example for your given expression.

((replacements =
Cases[PatternTable[expression],
{score_ /; score > 20, count_ /; count > 3, leaf_ /; leaf > 5, exp_}
])) // Column

 {
{156, 5, 40, 1/2 Cos[a] (m - h Cos[a]) + x3 Tan[a] -  x3 (-((m Sec[a])/(x2 - x3)) +Tan[a]) - Tan[b]},
{126, 10, 15, -((m Sec[a])/(x2 - x3)) + Tan[a]},
{102, 18, 7, m - h Cos[a]},
{99, 10, 12, -((m Sec[a])/(x2 - x3))},
{69, 4, 24, -((Cos[a] (m - h Cos[a]) Tan[b])/(2 (Tan[a] - Tan[b])))}},
{68, 5, 18, -x3 (-((m Sec[a])/(x2 - x3)) + Tan[a])},
.... + 7 more
}


If you take these and substitute into the expression it ends up reducing the leafcount, but not alot if you take into acount the leafcount needed to define the subesxpressions.

 subexpressions = Rule @@@ ({replacements[[1 ;;, -1]],
shortexp /@ Range[Length[replacements]]}\[Transpose])

subexpressions // LeafCount
exp /. subexpressions // LeafCount

244
200


So down from 693 to 444, but it may help as an aide in further simplifications.

-
Nice, but of course you can't overcome with this Mma's particular idiosyncrasy for representing expressions. Try e = (Sin[(x + y + z)]/(x + y + z))^(x + y + z) – Dr. belisarius Jul 27 '12 at 21:33
@belisarius I'm not sure what you mean, with that particular expression the only replacement with a positive score is x+y+z, which seems a reasonable substitution. – jVincent Jul 28 '12 at 17:51
Yes, but I am not getting it from your code above. Perhaps my fault. Rechecking – Dr. belisarius Jul 28 '12 at 18:07
You won't get it from using that definition of replacements, no. But the sorting is as I stated just a way to reduce the number of replacements, in your case the obvious thing would be to examine the entire list since it's just around 10 entries and notice that all but one are useless. – jVincent Jul 28 '12 at 22:17