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

Verbeia 对 Mathematica 中函数式编程风格的表现展开了一场相当有趣的讨论。在这里能找到它：在 Mathematica 中构建大型分块矩阵最有效的方法是什么？ https://stackoverflow.com/q/6867079/312124)

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

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

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

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

Update:符合什么acl推荐，我尝试使用Compile来加快我的表达速度。对于一个相对较小的L = 600，我天真地得到了 3.81 秒Do[...], 普通旧版 1.54 秒Compile，以及 0.033 秒Compile[..., CompilationTarget->"C"].

为了获得更真实的尺寸L = 1200，时间变为 16.68、0.605 和 0.132Do, Compile and Compile[.., CompilationTarget->"C"]分别。我能够实现 acl 在他的帖子中提到的相同 2 个数量级的加速。

Try Compile。这里我定义了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

giving

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

instead.

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

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