MATLAB中非常大矩阵的高效乘法

Question

我的内存不足以简单地创建一个大小为D的对角线D矩阵，因为D很大。我一直收到'out of memory'错误。

我尝试通过在第一次乘法中执行M x D操作而不是M x D x D操作来减少计算量，但是我的代码仍然运行缓慢。

有没有人能找到更有效的方法来执行乘法A'*B*A？以下是我迄今为止尝试过的内容：

D=20000
M=25

A = floor(rand(D,M)*10);
B = floor(rand(1,D)*10);

for i=1:D
    for j=1:M
        result(i,j) = A(i,j) * B(1,j);
    end
end    

manual = result * A';
auto = A*diag(B)*A';
isequal(manual,auto)

- matcheek

---

我有点困惑。矩阵B应该是D乘D还是M乘M？你的图片显示前者，但你的代码却表明后者。 - gnovice

另外，您是在尝试计算A'BA吗？这将给您一个M乘以M的结果。 - gnovice

是的，那就是我想要实现的。 - matcheek

Answer 1

解决您问题的一个选项是使用稀疏矩阵。以下是一个示例：

D = 20000;
M = 25;
A = floor(rand(D,M).*10);    %# A D-by-M matrix
diagB = rand(1,D).*10;       %# Main diagonal of B
B = sparse(1:D,1:D,diagB);   %# A sparse D-by-D diagonal matrix
result = (A.'*B)*A;         %'# An M-by-M result

另一种选择是使用REPMAT函数将D元素沿着B的主对角线复制，创建一个M x D矩阵，然后使用逐元素乘法与A.'相乘。

B = repmat(diagB,M,1);   %# Replicate diagB to create an M-by-D matrix
result = (A.'.*B)*A;    %'# An M-by-M result

还有另一种选择是使用函数BSXFUN：

result = bsxfun(@times,A.',diagB)*A;  %'# An M-by-M result

- gnovice

---

如果内存是一个问题，bsxfun优于repmat，因为它不会生成复制的矩阵。然而，bsxfun直到Matlab 2006才可用... - shabbychef

Answer 2

也许我有点脑抽了，但是你不能将你的DxD矩阵转换为一个DxM矩阵（其中包含给定向量的M个副本），然后.*最后两个矩阵而不是将它们相乘（然后，当然，正常地将第一个乘以找到的乘积数量）吗？

- Yonatan N

---

我不会创建 D x D 矩阵，因为会出现“内存错误”，如果那样的话就太容易了。你的解决方案几乎立即就给我报了“内存错误”，而我的需要更长一些时间，但确实完成了相同的任务。它们都是正确的，只是矩阵太大了。 - matcheek

我的解决方案和gnovice的解决方案都只需要O(DxM)的存储空间。我不明白为什么我的解决方案是错误的，而他的是正确的。事实上，他的repmat解决方案与我的完全相同（乘法顺序除外，但这并不重要）。 - Yonatan N

谢谢你的回答。为了明确，我尝试了你的解决方案，但是无法解决“内存不足”的错误。 - matcheek

Answer 3

1. 您之所以出现“内存不足”的情况，是因为MATLAB找不到足够大的内存块来容纳整个矩阵。有不同的技术可以避免这个错误，在MATLAB文档中有描述。

2. 在MATLAB中，通常情况下不需要编写显式循环，因为可以使用运算符*。如果使用显式循环进行矩阵乘法，存在一种加速技术，这里有C#示例。它有一个好主意，即如何将（可能很大的）矩阵分成较小的矩阵。要在MATLAB中包含这些较小的矩阵，可以使用单元矩阵。系统更有可能找到足够的RAM来容纳两个较小的子矩阵，而不是最终的大矩阵。

- Mikhail Poda

---

我是一个大无知，甚至从未想过缓存未命中、循环展开或分支预测，但感谢提醒。 - matcheek

@matcheek：将大矩阵分割成较小的矩阵的技术可以像应用于缓存一样应用于 RAM。对于缓存来说这无关紧要，但对于 RAM 来说很重要。 - Mikhail Poda