Find subexpression to minimize leafcount after replacment with temporary variable

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.