Mathematica：使用简化来进行常见的子表达式消除和强度降低

Question

因此，最近我一直在研究如何在编译器优化中很好地利用Mathematica的模式匹配和术语重写...试图高度优化作为循环内部部分的短代码块。减少评估表达式所需的工作量的两种常用方法是，识别多次出现的子表达式并存储结果，然后在随后的点上使用存储的结果来保存工作。另一种方法是尽可能使用便宜的操作。例如，我的理解是，与加法和乘法相比，求平方根需要更多的时钟周期。明确地说，我对评估表达式所需的浮点运算的成本感兴趣，而不是对Mathematica评估它需要多长时间。

[我的最初想法是，我将解决使用Mathematica的simplify函数开发的问题。可以指定比较两个表达式的相对简单性的复杂度函数。我将使用权重为相关的算术运算创建一个，并为表达式添加LeafCount，以说明所需的赋值运算。这解决了强度方面的降低，但是消除常见的子表达式使我绊倒了。

我当时正在考虑将常见的子表达式消除添加到可能简化使用的转换函数中。但是对于较大的表达式，可能会有许多可能的子表达式可以被替换，并且直到您看到该表达式，才可能知道它们是什么。我已经编写了一个函数，该函数给出了可能的替换，但是看来您指定的转换函数仅需要返回单个可能的转换，至少从文档中的示例中可以。关于如何克服这一限制的任何想法？有谁对如何简化使用可能暗示前进方向的转换函数有更好的想法？

我以为Simplify在幕后进行了一些动态编程，试图对表达式的不同部分进行不同的简化，并返回具有最低复杂度得分的表达式。我会更好地尝试使用常见的代数简化（例如factor和collect）独自进行这种动态编程吗？

编辑：我添加了生成可能要删除的子表达式的代码

(*traverses entire expression tree storing each node*)
AllSubExpressions[x_, accum_] := Module[{result, i, len},
  len = Length[x];
  result = Append[accum, x];
  If[LeafCount[x] > 1,
   For[i = 1, i <= len, i++,
     result = ToSubExpressions2[x[[i]], result];
     ];
   ];
  Return[Sort[result, LeafCount[#1] > LeafCount[#2] &]]
  ]

CommonSubExpressions[statements_] := Module[{common, subexpressions},
subexpressions = AllSubExpressions[statements, {}];
  (*get the unique set of sub expressions*)
  common = DeleteDuplicates[subexpressions];
  (*remove constants from the list*)
  common = Select[common, LeafCount[#] > 1 &];
  (*only keep subexpressions that occur more than once*)
  common = Select[common, Count[subexpressions, #] > 1 &];
  (*output the list of possible subexpressions to replace with the \
number of occurrences*)
  Return[common];
  ]

一旦从CommonSubExpressions返回的列表中选择了一个公共子表达式，下面将执行替换操作。

eliminateCSE[statements_, expr_] := Module[{temp},
temp = Unique["r"];
Prepend[ReplaceAll[statements, expr -> temp], temp[expr]]
]

冒这个问题的风险，我将举一些示例代码。我认为尝试优化的一个合适的表达式是用于求解微分方程的经典Runge-Kutta方法。

Input:
nextY=statements[y + 1/6 h (f[t, n] + 2 f[0.5 h + t, y + 0.5 h f[t, n]] + 
    2 f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]] + 
    f[h + t, 
     y + h f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]]])];
possibleTransformations=CommonSubExpressions[nextY]
transformed=eliminateCSE[nextY, First[possibleTransformations]]

Output:
{f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]], 
y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]], 
0.5 h f[0.5 h + t, y + 0.5 h f[t, n]], 
f[0.5 h + t, y + 0.5 h f[t, n]], y + 0.5 h f[t, n], 0.5 h f[t, n], 
0.5 h + t, f[t, n], 0.5 h}

statements[r1[f[0.5 h + t, y + 0.5 h f[0.5 h + t, y + 0.5 h f[t, n]]]], 
y + 1/6 h (2 r1 + f[t, n] + 2 f[0.5 h + t, y + 0.5 h f[t, n]] + 
 f[h + t, h r1 + y])]

最后，下面是判断不同表达式的相对代价的代码。权重在这一点上是概念性的，因为这仍然是我正在研究的领域。

Input:
cost[e_] := 
 Total[MapThread[
   Count[e, #1, Infinity, Heads -> True]*#2 &, {{Plus, Times, Sqrt, 
     f}, {1, 2, 5, 10}}]]
cost[transformed]

Output:
100

Answer 1

要确定重复的子表达式，可以使用类似这样的东西

(*helper functions to add Dictionary-like functionality*)

index[downvalue_, 
   dict_] := (downvalue[[1]] /. HoldPattern[dict[x_]] -> x) // 
   ReleaseHold;
value[downvalue_] := downvalue[[-1]];
indices[dict_] := 
  Map[#[[1]] /. {HoldPattern[dict[x_]] -> x} &, DownValues[dict]] // 
   ReleaseHold;
values[dict_] := Map[#[[-1]] &, DownValues[dict]];
items[dict_] := Map[{index[#, dict], value[#]} &, DownValues[dict]];
indexQ[dict_, index_] := 
  If[MatchQ[dict[index], HoldPattern[dict[index]]], False, True];


(*count number of times each sub-expressions occurs *)
expr = Cos[x + Cos[Cos[x] + Sin[x]]] + Cos[Cos[x] + Sin[x]];
Map[(counts[#] = If[indexQ[counts, #], counts[#] + 1, 1]; #) &, expr, 
  Infinity];
items[counts] // Column

Answer 2

作者在这里还实现了一些例程：http://stoney.sb.org/wordpress/2009/06/converting-symbolic-mathematica-expressions-to-c-code/

我将其打包到* .M文件中，并修复了一个错误（如果该表达式没有重复的子表达式，它将终止。），我正在尝试查找作者的联系信息，以查看是否可以将其修改后的代码上传到pastebin或任何地方

编辑：我已经获得作者的上传许可，并将其粘贴到此处：http://pastebin.com/fjYiR0B3

Answer 3

我试图模仿字典压缩功能出现在此博客上：https://writings.stephenwolfram.com/2018/11/logic-explainability-and-the-future-of-understanding/

这是我的作品：

DictionaryCompress[expr_, count_, size_, func_] := Module[
  {t, s, rule, rule1, rule2},
  t = Tally@Level[expr, Depth[expr]];
  s = Sort[
    Select[{First@#, Last@#, Depth[First@#]} & /@ 
      t, (#[[2]] > count && #[[3]] > size) &], #1[[2]]*#1[[3]] < #2[[
        2]]*#2[[2]] &];
  rule = MapIndexed[First[#1] -> func @@ #2 &, s];
  rule = (# //. Cases[rule, Except[#]]) & /@ rule;
  rule1 = Select[rule, ! FreeQ[#, Plus] &];
  rule2 = Complement[rule, rule1];
  rule = rule1 //. (Reverse /@ rule2);
  rule = rule /. MapIndexed[ Last[#1] -> func @@ #2 &, rule];
  {
   expr //. rule,
   Reverse /@ rule
   }
  ];

poly = Sum[Subscript[c, k] x^k, {k, 0, 4}];
sol = Solve[poly == 0, x];
expr = x /. sol;
Column[{Column[
     MapIndexed[
      Style[TraditionalForm[Subscript[x, First[#2]] == #], 20] &, #[[
       1]]], Spacings -> 1],
    Column[Style[#, 20] & /@ #[[2]], Spacings -> 1, Frame -> All]
    }] &@DictionaryCompress[expr, 1, 1, 
  Framed[#, Background -> LightYellow] &]