recursion - 在 Mathematica 中计算此递推关系的更有效方法

Question

Verbeia 就 Mathematica 中函数式编程风格的性能展开了一场相当有趣的讨论。可以在这里找到:What is the most efficient way to construct large block matrices in Mathematica? )

我正在解决一个问题，在对代码进行一些计时后，一个特别耗时的部分是我通过递归关系计算矩阵的条目:

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是一些外部定义的列表。有什么办法可以更有效地写这个吗？我似乎找不到任何明显的方法来使用内置函数以更惯用和更有效的方式完成此任务。

我意识到我只能做这么多，因为这段代码是O(n^2)，但是我有一系列矩阵乘法(总共大约6次)，组合起来需要更少的时间比这个语句运行的时间。我能做的任何事情来加快速度，哪怕只是一点点，都会对我的运行时间产生明显的影响。

更新:根据acl的建议，我尝试使用Compile来加速我的表达式。对于相对较小的 L = 600，我在简单的 Do[...] 上得到了 3.81 秒，在普通的旧 Compile 上得到了 1.54 秒， Compile[..., CompilationTarget->"C"] 为 0.033 秒。

对于更实际的 L = 1200 大小，Do、Compile 和 的计时变为 16.68、0.605 和 0.132 >分别编译[.., CompilationTarget->"C"]。我能够实现 acl 在他的帖子中提到的相同 2 个数量级的加速。

Answer 1

尝试编译。这里我定义了 3 个函数:f 正如你所定义的那样，fc 编译为(某种字节码)，fcc 编译为 C(查找有关如何检查生成的代码的文档)。

首先，让 mma 告诉我们是否有东西无法编译:

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 上进行测试):

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 的版本在这里快了两个数量级。

如果您更改类型并出现编译错误，请询问。

编辑:如果e包含实数，请尝试

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"
];

相反。

人们可以通过说来感受一下它是如何工作的

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 代码(但这有点难以阅读)。