在 Mathematica 中计算这种递推关系的更有效方法 [英] More efficient way of calculating this recurrence relation in Mathematica
问题描述
Verbeia 就 Mathematica 中函数式编程风格的性能展开了一个相当有趣的讨论.可在此处找到:在 Mathematica 中构建大型块矩阵的最有效方法是什么?)
我正在解决一个问题，在对我的代码进行了一些计时之后，一个特别耗时的部分是我通过递归关系计算矩阵的条目:
I'm working on a problem, and after doing some timing of my code one particularly time consuming portion is where I calculate entries of a matrix through a recurrence relation:
c = Table[0, {2L+1}, {2L+1}];
c[[1, 1]] = 1;
Do[c[[i, i]] = e[[i - 1]] c[[i - 1, i - 1]], {i, 2, 2 L + 1}];
Do[c[[i, 1]] = (1 - e[[i - 1]]) c[[i - 1, 1]], {i, 2, 2 L + 1}];
Do[c[[i, j]] = (1 - e[[i - 1]]) c[[i - 1, j]] +
e[[i - 1]] c[[i - 1, j - 1]], {i, 2, 2 L + 1}, {j, 2, i - 1}];
其中 e 是一些外部定义的列表.有什么办法可以更有效地写这个吗?我似乎找不到任何明显的方法来使用内置函数以更惯用和有效的方式完成此任务.
Where e is some externally defined list. Is there any way I could write this in a more efficient manner? I can't seem to find any obvious way of using the built in functions to accomplish this in a more idiomatic and efficient way.
我意识到我只能做这么多，因为这段代码是O(n^2)，但我有一系列矩阵乘法(总共大约 6 个)，组合起来，需要更少时间比这个语句运行.我可以做的任何事情来加快速度，哪怕只是一点点，都会对我的运行时间产生明显的影响.
I realize I can only do so much, since this code is O(n^2), but I have a series of matrix multiplications (About 6 in all) that, combined, take less time to run than this statement. Anything I can do to speed this up even slightly would make an appreciable difference in run times for me.
更新:根据 acl 的建议，我尝试使用 Compile 来加速我的表达式.对于相对较小的 L = 600，我在 naive Do[...] 上得到 3.81 秒，对于普通的 Compile，我得到 1.54 秒，Compile[..., CompilationTarget->"C"] 为 0.033 秒.
Update:
In line with what acl recommended, I tried using Compile to speed up my expressions. For a relatively small L = 600, I get 3.81 seconds on the naive Do[...], 1.54 seconds for plain old Compile, and 0.033 seconds for Compile[..., CompilationTarget->"C"].
对于 L = 1200 的更实际大小，Do、Compile 和 的时间分别为 16.68、0.605 和 0.132>Compile[.., CompilationTarget->"C"] 分别.我能够实现与 acl 在他的帖子中提到的相同的 2 个数量级的加速.
For a more realistic size of L = 1200, the timings become 16.68, 0.605, and 0.132 for Do, Compile and Compile[.., CompilationTarget->"C"] respectively. I'm able to achieve the same 2 orders of magnitude speedup that acl mentioned in his post.
推荐答案
尝试 Compile.在这里我定义了 3 个函数:f 正如你定义的那样，fc 编译(到某种字节码)和 fcc 编译到 C(查找有关如何检查生成的代码的文档).
Try Compile. Here I define 3 functions: f as you defined it, fc compiled (to some sort of bytecode) and fcc compiled to C (look up the documentation as to how to examine the generated code).
首先，让 mma 告诉我们是否有无法编译的内容:
First, make mma tell us if something can't be compiled:
SetSystemOptions["CompileOptions"->"CompileReportExternal"->True]
然后是函数:
ClearAll[f];
f=Function[{ell,e},
Module[{c=Table[0,{2ell+1},{2ell+1}]},
c[[1,1]]=1;
Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
{i,2,2 ell+1},{j,2,i-1}];
c
]
];
ClearAll[fc];
fc=Compile[{{ell,_Integer},{e,_Integer,1}},
Module[{c},
c=Table[0,{2ell+1},{2ell+1}];
c[[1,1]]=1;
Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
{i,2,2 ell+1},{j,2,i-1}];
c
]
];
ClearAll[fcc];
fcc=Compile[{{ell,_Integer},{e,_Integer,1}},
Module[{c},
c=Table[0,{2ell+1},{2ell+1}];
c[[1,1]]=1;
Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
{i,2,2 ell+1},{j,2,i-1}];
c
],
CompilationTarget->"C",
RuntimeOptions->"Speed"
];
没有错误，所以没关系.现在测试(这些在 2.4GHz core 2 duo 的 macbook 上使用电池运行):
no errors, so it's OK. And now test (these on a macbook with a 2.4GHz core 2 duo running on battery):
ell=400;
e=RandomInteger[{0,1},2*ell];
f[ell,e];//Timing
fc[ell,e];//Timing
fcc[ell,e];//Timing
给予
{2.60925, Null}
{0.092022, Null}
{0.022709, Null}
所以编译成 C 的版本在这里要快两个数量级.
so the version compiled to C is two orders of magnitude faster here.
如果您更改类型并出现编译错误，请询问.
If you change the types and get compilation errors, ask.
如果 e 包含实数，请尝试
If e contains reals, try
ClearAll[fc];
fc=Compile[{{ell,_Integer},{e,_Real,1}},
Module[{c},
c=Table[0.,{2ell+1},{2ell+1}];
c[[1,1]]=1;
Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
{i,2,2 ell+1},{j,2,i-1}];
c
]
];
ClearAll[fcc];
fcc=Compile[{{ell,_Integer},{e,_Real,1}},
Module[{c},
c=Table[0.,{2ell+1},{2ell+1}];
c[[1,1]]=1;
Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
{i,2,2 ell+1},{j,2,i-1}];
c
],
CompilationTarget\[Rule]"C",
RuntimeOptions\[Rule]"Speed"
];
相反.
人们可以通过说来感受一下这是如何工作的
One can get a feel for how this works by saying
Needs["CompiledFunctionTools`"]
CompilePrint[fc]
并获得
2 arguments
18 Integer registers
6 Real registers
3 Tensor registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}
I0 = A1
T(R1)0 = A2
I12 = 0
I1 = 2
I3 = 1
I14 = -1
R0 = 0.
Result = T(R2)2
1 I11 = I1 * I0
2 I11 = I11 + I3
3 I7 = I1 * I0
4 I7 = I7 + I3
5 I17 = I12
6 T(R2)2 = Table[ I11, I7]
7 I15 = I12
8 goto 13
9 I16 = I12
10 goto 12
11 Element[ T(R2)2, I17] = R0
12 if[ ++ I16 < I7] goto 11
13 if[ ++ I15 < I11] goto 9
14 R1 = I3
15 Part[ T(R2)2, I3, I3] = R1
16 I4 = I1 * I0
17 I4 = I4 + I3
18 I5 = I3
19 goto 26
20 I10 = I5 + I14
21 I8 = I10
22 R4 = Part[ T(R1)0, I8]
23 R5 = Part[ T(R2)2, I8, I8]
24 R4 = R4 * R5
25 Part[ T(R2)2, I5, I5] = R4
26 if[ ++ I5 < I4] goto 20
27 I4 = I1 * I0
28 I4 = I4 + I3
29 I5 = I3
30 goto 40
31 I10 = I5 + I14
32 I13 = I10
33 R4 = Part[ T(R1)0, I13]
34 R5 = - R4
35 R4 = I3
36 R4 = R4 + R5
37 R5 = Part[ T(R2)2, I13, I3]
38 R4 = R4 * R5
39 Part[ T(R2)2, I5, I3] = R4
40 if[ ++ I5 < I4] goto 31
41 I4 = I1 * I0
42 I4 = I4 + I3
43 I5 = I3
44 goto 63
45 I6 = I5 + I14
46 I17 = I3
47 goto 62
48 I16 = I5 + I14
49 I9 = I16
50 R4 = Part[ T(R1)0, I9]
51 R2 = R4
52 R4 = - R2
53 R5 = I3
54 R5 = R5 + R4
55 R4 = Part[ T(R2)2, I9, I17]
56 R5 = R5 * R4
57 I16 = I17 + I14
58 R4 = Part[ T(R2)2, I9, I16]
59 R3 = R2 * R4
60 R5 = R5 + R3
61 Part[ T(R2)2, I5, I17] = R5
62 if[ ++ I17 < I6] goto 48
63 if[ ++ I5 < I4] goto 45
64 Return
寄存器的名称表明它们的类型等.如果您使用C"选项，您还可以查看生成的 C 代码(但这有点难以阅读).
with the names of the registers indicating their type etc. You can also look at the generated C code if you use the "C" option (but that is a bit harder to read).
这篇关于在 Mathematica 中计算这种递推关系的更有效方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
