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

【问题标题】：Mathematica: using simplify to do common sub-expression elimination and reduction in strengthMathematica：使用简化来做常见的子表达式消除和强度降低
【发布时间】：2010-11-17 08:56:45
【问题描述】：

所以最近我一直在研究如何在编译器优化中很好地利用 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

【问题讨论】：

嗨！你介意分享你的提供可能替换的函数吗？
Tnx！我发现“Experimental`OptimizeExpression”......这似乎消除了重复的计算。不过，仍然不适合我。
顺便说一句，Simplify 更像是启发式搜索，它尝试不同的转换规则以使表达式更短。它的强大之处在于知道许多特殊的函数关系，但我经常遇到只有 exp,+,-,/,* 的表达式，Simplify 找不到最简单的形式，而一些自定义重写规则 a->b 可以完成这项工作

标签： wolfram-mathematica compiler-optimization

【解决方案1】：

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

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

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

【讨论】：

工作链接stoney.sb.org/wordpress/…

【解决方案2】：

要识别重复的子表达式，您可以使用类似这样的方法

(*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

【讨论】：

不错！我想知道构建一个基本的优化器还剩下多少工作。

【解决方案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] &]

【讨论】：