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

Question

受Mike Bantegui 关于构造定义为递归关系的矩阵的问题的启发,我想知道是否可以在最少的计算时间内设置大块矩阵.根据我的经验,构建块然后将它们放在一起可能效率很低(因此我的答案实际上比Mike的原始代码慢).Join并且可能ArrayFlatten效率低于它们.

显然,如果矩阵是稀疏的,可以使用SparseMatrix构造,但有时你构造的块矩阵不是稀疏的.

这类问题的最佳做法是什么？我假设矩阵的元素是数字.

Answer 1

下面的代码可以在这里找到: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循环替换它,但我保留它因为它看起来更干净.

天真Do循环的当前计数为9.06秒,使用Outer和的下一个版本为1.389秒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:似乎对于这些情况,功能方式优于程序方式.但是在编译时,过程代码优于功能代码.