问题 在Mathematica中构造大块矩阵的最有效方法是什么？

下面显示的代码可在此处获得： http://pastebin.com/4PWWxGhB。只需将其复制并粘贴到笔记本中即可进行测试。

我实际上是在尝试做几种计算矩阵的函数方法，因为我 认为功能方式（在Mathematica中通常是惯用的）更有效。

作为一个例子，我有这个由两个列表组成的矩阵：

In: L = 1200;
e = Table[..., {2L}];
f = Table[..., {2L}];

h = Table[0, {2L}, {2L}];
Do[h[[i, i]] = e[[i]], {i, 1, L}];
Do[h[[i, i]] = e[[i-L]], {i, L+1, 2L}];
Do[h[[i, j]] = f[[i]]f[[j-L]], {i, 1, L}, {j, L+1, 2L}];
Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];

我的第一步是计时。

In: h = Table[0, {2 L}, {2 L}];
AbsoluteTiming[Do[h[[i, i]] = e[[i]], {i, 1, L}];]
AbsoluteTiming[Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}];]
AbsoluteTiming[
 Do[h[[i, j]] = f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}];]
AbsoluteTiming[Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];]

Out: {0.0020001, Null}
{0.0030002, Null}
{5.0012861, Null}
{4.0622324, Null}

DiagonalMatrix[...] 比do循环慢，所以我决定只使用 Do 最后一步循环。如你所见，使用 Outer[Times, f, f] 在这种情况下要快得多。

然后我写了等效的使用 Outer 对于矩阵右上角和左下角的块，和 DiagonalMatrix 对角线：

AbsoluteTiming[h1 = ArrayPad[Outer[Times, f, f], {{0, L}, {L, 0}}];]
AbsoluteTiming[h1 += Transpose[h1];]
AbsoluteTiming[h1 += DiagonalMatrix[Join[e, e]];]


Out: {0.9960570, Null}
{0.3770216, Null}
{0.0160009, Null}

该 DiagonalMatrix 实际上比较慢。我可以用这个代替 Do 循环，但我保留它，因为它看起来更干净。

天真的当前计数是9.06秒 Do 循环，我的下一个版本使用1.389秒 Outer 和 DiagonalMatrix。大约6.5倍的加速，还不错。

听起来快很多，现在不是吗？我们试试吧 Compile 现在。

In: cf = Compile[{{L, _Integer}, {e, _Real, 1}, {f, _Real, 1}},
   Module[{h},
    h = Table[0.0, {2 L}, {2 L}];
    Do[h[[i, i]] = e[[i]], {i, 1, L}];
    Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}];
    Do[h[[i, j]] = f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}];
    Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];
    h]];

AbsoluteTiming[cf[L, e, f];]

Out: {0.3940225, Null}

现在它的运行速度比我上一版本快3.56倍，比第一版快23.23倍。下一个版本：

In: cf = Compile[{{L, _Integer}, {e, _Real, 1}, {f, _Real, 1}},
   Module[{h},
    h = Table[0.0, {2 L}, {2 L}];
    Do[h[[i, i]] = e[[i]], {i, 1, L}];
    Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}];
    Do[h[[i, j]] = f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}];
    Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];
    h], CompilationTarget->"C", RuntimeOptions->"Speed"];

AbsoluteTiming[cf[L, e, f];]

Out: {0.1370079, Null}

大多数速度来自 CompilationTarget->"C"。在这里，我获得了超过最快版本的2.84加速，比第一版高66.13倍。但我所做的只是编译它！

现在，这是一个非常简单的例子。但这是我用来解决凝聚态物理问题的真实代码。因此，不要将其视为“玩具示例”。

关于我们可以使用的技术的另一个例子怎么样？我有一个相对简单的矩阵我必须建立起来。我有一个矩阵，它只包含从开始到某个任意点的矩阵。天真的方式可能看起来像这样：

In: k = L;
AbsoluteTiming[p = Table[If[i == j && j <= k, 1, 0], {i, 2L}, {j, 2L}];]
Out: {5.5393168, Null}

相反，让我们使用它来构建它 ArrayPad 和 IdentityMatrix：

In: AbsoluteTiming[ArrayPad[IdentityMatrix[k], {{0, 2L-k}, {0, 2L-k}}
Out: {0.0140008, Null}

这实际上不适用于k = 0，但你可以特殊情况，如果你需要它。此外，根据k的大小，这可以更快或更慢。它总是比Table [...]版本更快。

你甚至可以用这个来写 SparseArray：

In: AbsoluteTiming[SparseArray[{i_, i_} /; i <= k -> 1, {2 L, 2 L}];]
Out: {0.0040002, Null}

我可以继续谈论其他一些事情，但我担心如果我这样做，我会让这个答案不合理地大。在我尝试优化某些代码时，我已经积累了许多用于形成这些不同矩阵和列表的技术。我使用的基本代码需要6天才能运行一个计算，现在只需要6个小时就可以完成相同的操作。

我会看看我是否可以选择我想出的一般技术，然后将它们粘在笔记本上即可使用。

TL; DR：似乎对于这些情况，功能方式优于程序方式。但是在编译时，过程代码优于功能代码。