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

Question

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

我的第一个想法是，我将使用 Mathematica 的 simplify 函数来解决开发问题。可以指定一个复杂性函数来比较两个表达式的相对简单性。我打算使用相关算术运算的权重创建一个，并将表达式的 LeafCount 添加到其中，以考虑所需的赋值操作。这解决了强度方面的降低问题，但正是消除了常见的子表达式让我陷入了困境。

我正在考虑将公共子表达式消除添加到可能的转换函数中以简化使用。但对于大型表达式，可能有许多可能的子表达式可以被替换，并且在看到表达式之前不可能知道它们是什么。我编写了一个提供可能替换的函数，但您指定的转换函数似乎只需要返回一个可能的转换，至少从文档中的示例来看是这样。关于如何绕过这一限制有什么想法吗？有谁更好地了解简化如何使用可能暗示前进方向的转换函数？

我想象 Simplify 在幕后正在进行一些动态编程，尝试对表达式的不同部分进行不同的简化，并返回复杂度分数最低的简化。我尝试使用常见的代数简化（例如因子和集合）自己进行动态规划会更好吗？

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

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

通过添加选项来过滤掉列表中作为其他元素的子表达式的元素，对 @Yaroslav 的答案进行了轻微的补充：

subexpressionQ[a_, b_] := MemberQ[Cases[a, _, All], b]

filterOutSubexpressions[x_] := 
 Select[x, 
  s |-> Not@
    AnyTrue[x[[All, 1]] // 
      DeleteCases[s[[1]]], subexpressionQ[#, s[[1]]] &]]

Yaroslav 封装到函数中的代码：

subExpressionCount[s_] :=
 Module[{index,value,indices,values,items,indexQ,counts},
  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];
  Map[(counts[#]=If[indexQ[counts,#],counts[#]+1,1];#)&,s,Infinity];
  items[counts]
  ]

注意： 我不清楚排序后 subExpressionCount 的输出是否可以与 Tally@Cases[#, _, All] & 不同。

示例：

subExpressionCount[
    x /. Solve[a*x^3 + b*x^2 + c*x + d == 0, x][[1]]] // 
   Select[#[[2]] > 1 && 
      LeafCount[#] > 10 &] // filterOutSubexpressions // 
 ReverseSortBy[Last]

{{-2 b^3 + 9 abc - 27 a^2 d + Sqrt[ 4 (-b^2 + 3 ac)^3 + (-2 b^3 + 9 abc - 27 a^ 2 d)^2], 2}}

Answer 2

我试图模仿此博客上出现的字典压缩功能： 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] &]

Answer 3

要识别重复的子表达式，您可以使用类似的东西

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

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

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

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