3
$\begingroup$

I want to have common subexpression elimination of comlicated functions where the elimination is done by factoring out the expressions. The result is not to be digested by a compiler, it should remain symbolic. I want this to work in general, without choosing things by hand.

As an example of what I want to do, take to following expression:

-3 a - 2 a^3 + 4 Sqrt[1 + a^2] (5 - 9 Log[2]) + 
 4 a^2 Sqrt[1 + a^2] (5 - 9 Log[2]) + 
 12 (1 + a^2)^(3/2) Log[1 + Sqrt[1 + 1/a^2]] - 
 6 (4 (Sqrt[1 + a^2] - a (2 + a^2 - a Sqrt[1 + a^2])) Log[a] + 
    a Log[1 + a^2])

$-2 a^3+4 \sqrt{a^2+1} a^2 (5-9 \log (2))+12 \left(a^2+1\right)^{3/2} \log \left(\sqrt{\frac{1}{a^2}+1}+1\right)-6 \left(4 \left(\sqrt{a^2+1}-a \left(a^2-\sqrt{a^2+1} a+2\right)\right) \log (a)+a \log \left(a^2+1\right)\right)+4 \sqrt{a^2+1} (5-9 \log (2))-3 a$

This will not change by doing FullSimplify[ExpandAll[%], Assumptions -> _Symbol \[Element] Reals].

What I am after is to transform this expression to something where the complicated terms (e.g. measured in LeafCount) are nested. What I came up with until now is something which mostly works. However in terms of programming it is probably the ugliest function written in recent human history. I would like to ask if somebody has an idea how to do things nicer. Perhaps using OptimizeExpression is actually the wrong thing to do.

nestComplicatedTerms[x_] := Block[
  {expr, compiled, exprList, exprList2, countList}, {};
  compiled = 
   Experimental`OptimizeExpression[x, OptimizationLevel -> 2];
  compiled = Map[Hold, compiled, {2}][[1, 2]];
  compiled = ToString[InputForm[compiled]];
  compiled = 
   StringReplace[
    compiled, {"=" -> "->", ";" -> ",", "Hold" -> "List"}];
  exprList = ToExpression[compiled];
  exprList2 = exprList;
  countList = Table[(
     exprList = exprList /. exprList[[i]];
     {exprList2[[i]][[1]], 
      LeafCount[exprList[[i]][[2]]]*
       Boole[Count[exprList2[[-1]], exprList2[[i]][[1]], {-1}] > 0]}
     ), {i, 1, Length[exprList] - 1}];
  countList = SortBy[countList, -Last[#] &];
  countList = Table[countList[[i]][[1]], {i, 1, Length[countList]}];
  exprList = exprList2;
  expr = HornerForm[exprList[[-1]], countList];
  Do[(
    expr = expr /. exprList[[i]];
    exprList = exprList /. exprList[[i]];
    ), {i, 1, Length[exprList] - 1}];
  expr
  ]

What it does is to first produce the output from OptimizeExpression. Then comes the ugly step. The output is transformed to a string, replacing some characters and then transfomed back to an expression. I had to do this, because I was unable to extract the replacement rules without actually replacing the variables in the function.

This is how the output from OptimizeExpression looks like:

Experimental`OptimizedExpression[
 Block[{Compile`$938, Compile`$948, Compile`$954, Compile`$955, 
   Compile`$956, Compile`$957, Compile`$958, Compile`$961, 
   Compile`$962, Compile`$963, Compile`$964},
  Compile`$938 = a^2; 
      Compile`$948 = a Compile`$938;
  Compile`$954 = 1 + Compile`$938; 
      Compile`$955 = Sqrt[Compile`$954];
  Compile`$956 = Log[2]; 
      Compile`$957 = -9 Compile`$956;
  Compile`$958 = 5 + Compile`$957; 
      Compile`$961 = Compile`$955 Compile`$954; 
      Compile`$962 = 1/Compile`$938;
  Compile`$963 = 1 + Compile`$962; 
      Compile`$964 = Sqrt[Compile`$963];
  -3 a - 2 Compile`$948 + 
       4 Compile`$955 Compile`$958 + 
   4 Compile`$938 Compile`$955 Compile`$958 + 
   12 Compile`$961 Log[1 + Compile`$964] - 
       6 (4 (Compile`$955 - a (2 + Compile`$938 - a Compile`$955)) Log[
            a] + a Log[Compile`$954])]]

It is a mathematica coding expression, which actually is being processed when I want to extract the replacement rules, changing its own content when I do.. How could I convert this output to a list of replacement rules, like

{Compile`$938 -> a^2, Compile`$948 -> a Compile`$938, ..}) 

other than doing the forth and back string conversion?

The rest of the function is straightforward and does what I want. It is these first couple of lines which I do not like.

Oh and what it does to the expression I actually wanted to process? It decreases the LeafCount from 114 down to 89:

-3 a - 2 a^3 + 12 (1 + a^2)^(3/2) Log[1 + Sqrt[1 + 1/a^2]] + 
 48 a Log[a] + 24 a^3 Log[a] + 
 Sqrt[1 + a^2] ((4 + 4 a^2) (5 - 9 Log[2]) - 24 Log[a] - 
    24 a^2 Log[a]) - 6 a Log[1 + a^2]

$-2 a^3+24 a^3 \log (a)-6 a \log \left(a^2+1\right)+12 \left(a^2+1\right)^{3/2} \log \left(\sqrt{\frac{1}{a^2}+1}+1\right)+\sqrt{a^2+1} \left(-24 a^2 \log (a)+\left(4 a^2+4\right) (5-9 \log (2))-24 \log (a)\right)-3 a+48 a \log (a)$

Doing a FullSimplify on it further simplifies the expression to a LeafCount of 72:

-3 a - 2 a^3 + 12 (1 + a^2)^(3/2) Log[1 + Sqrt[1 + 1/a^2]] + 
 48 a Log[a] + 24 a^3 Log[a] - 
 4 (1 + a^2)^(3/2) (-5 + Log[512] + 6 Log[a]) - 6 a Log[1 + a^2]

$-2 a^3+24 a^3 \log (a)-6 a \log \left(a^2+1\right)+12 \left(a^2+1\right)^{3/2} \log \left(\sqrt{\frac{1}{a^2}+1}+1\right)-4 \left(a^2+1\right)^{3/2} (6 \log (a)-5+\log (512))-3 a+48 a \log (a)$

And if one repeats the procedure on this last output by FullSimplify[nestComplicatedTerms[%]] the expression is again reduced to a LeafCount of 61.

-3 a + 48 a Log[a] - 
 4 (1 + a^2)^(
  3/2) (-5 + Log[512] - 3 Log[1 + Sqrt[1 + 1/a^2]] + 6 Log[a]) + 
 a^3 (-2 + 24 Log[a]) - 6 a Log[1 + a^2]

$a^3 (24 \log (a)-2)-6 a \log \left(a^2+1\right)-4 \left(a^2+1\right)^{3/2} \left(-3 \log \left(\sqrt{\frac{1}{a^2}+1}+1\right)+6 \log (a)-5+\log (512)\right)-3 a+48 a \log (a)$

I do not know why this is not directly done the first time, though. Anyhow comparing to the initial expression it is much shorter and has way less repeating terms.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey May 29 '15 at 21:38
3
$\begingroup$

Ok so I decided to change my initial code with something less hacky. I found the relevant idea on a page. I am not sure I am allowed to post the link here, so I will put the reference in a comment. The first step is now the following function.

getCommonSubexpressions[expr_, OptionsPattern[]] := Block[
  {vnum, modExpr, subExprs, subExpr, mapping, inactivemap, 
   my`$Inactive$}, {
   vnum = 0,
   modExpr = expr,
   subExprs = 
    Join[Level[expr, {1, -2}], 
     Cases[expr, 
      Except[__Symbol?(Context@# === "System`" &), _Symbol], -1]]
   };
  If[subExprs == {}, Return[{modExpr, {}}]];

  (*This part is for being sure that functions do not get evaluated \
while containing the temporary variables. 
  my`$Inactive$ must be not defined as a function*)
  inactivemap = 
   DeleteDuplicates[
    Map[Head[#] -> my`$Inactive$[Head[#]] &, 
     Cases[subExprs, 
      Except[_Symbol | _Integer | _Real | _Rational | _Complex | \
_Plus | _Times | _Power | _List], {1}]]];
  inactivemap = inactivemap /. my`$Inactive$ -> Inactive;
  modExpr = modExpr /. inactivemap;
  subExprs = subExprs /. inactivemap;

  subExprs = Sort[subExprs, LeafCount[#1] < LeafCount[#2] &];
  Reap[While[subExprs != {},
    subExpr = First[subExprs];
    subExprs = Rest[subExprs];
    If[MemberQ[subExprs, subExpr] || 
      MemberQ[inactivemap, Head[subExpr], {-2}],
     mapping = subExpr -> Symbol["my`$" <> IntegerString[vnum++]];
     Sow[mapping[[2]] -> mapping[[1]]];
     modExpr = modExpr /. mapping;
     subExprs = DeleteCases[subExprs, subExpr] /. mapping]];
   modExpr]
  ]

The second step ist mostly like it was before, with the change that it seems advantageous to rather use Collect instead of my initial HornerForm, because HornerForm does fail when the variables contain roots while collect still works.

nestComplicatedTerms[x_] := Block[
  {expr, exprList, exprList2, countList}, {
   exprList = getCommonSubexpressions[x];
   };
  If[exprList[[2]] == {}, Return[x]];
  exprList = Flatten[Extract[exprList, {{2, 1}, {1}}]];
  exprList2 = exprList;
  countList = Table[(
     exprList = exprList /. exprList[[i]];
     {exprList2[[i]][[1]], 
      LeafCount[exprList[[i]][[2]]]*
       Boole[Count[exprList2[[-1]], exprList2[[i]][[1]], {-1}] > 0]}
     ), {i, 1, Length[exprList] - 1}];
  countList = Transpose[SortBy[countList, -Last[#] &]][[1]];
  exprList = exprList2;
  expr = Collect[exprList[[-1]], countList];
  Do[(
    expr = expr /. exprList[[i]];
    exprList = exprList /. exprList[[i]];
    ), {i, 1, Length[exprList] - 1}];
  Activate[expr]
  ]

This code will work as-is only in v10 and above, because I make use of the Inactivate mechanism. It could however be easily removed, which will make it usefull also for lowever versions, with the consrain that it will not work if the expression contains some rather unusual functions like MeijerG.

$\endgroup$
0
$\begingroup$

What constitutes simplification is in the eye of the beholder. With the initial expression in the Question designated z,

Collect[z, {(_)^(3/2), Sqrt[_]}, FullSimplify]
(* -4*(1 + a^2)^(3/2)*(-5 + Log[512] - 3*Log[1 + Sqrt[1 + a^(-2)]] + 6*Log[a]) 
   + a*(-3 - 2*a^2 + 24*(2 + a^2)*Log[a] - 6*Log[1 + a^2]) )*

produces an expression with a LeafCount of 60. Of course, some human insight went into selecting which functions to Collect. A more automated approach might identify all functions, Collect all permutations of them, and selection the one with smallest LeafCount.

$\endgroup$

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.