如何在Matlab中计算两个矩阵的外积减去一个公共矩阵的平方和

Question

假设有三个n*n矩阵X、Y、S。如何快速计算以下标量b

for i = 1:n
b = b  + sum(sum((X(i,:)' * Y(i,:) - S).^2));
end

计算成本为O(n^3(。存在一种快速计算两个矩阵的外积的方法。具体来说，矩阵C

for i = 1:n
C = C + X(i,:)' * Y(i,:);
end

可以在没有仅为O(n^2(的循环CCD_ 6的情况下计算。有没有一种更快的方法来计算b？

Answer 1

您可以使用：

X2 = X.^2;
Y2 = Y.^2;
S2 = S.^2;
b = sum(sum(X2.' * Y2 - 2 * (X.' * Y ) .* S + n * S2));

以为例

b=0;
for i = 1:n
b = b  + sum(sum((X(i,:).' * Y(i,:) - S).^2));
end

我们可以首先将总和从循环中取出：

b=0;
for i = 1:n
b = b  + (X(i,:).' * Y(i,:) - S).^2;
end
b=sum(b(:))

知道我们可以把(a - b)^2写成a^2 - 2*a*b + b^2

b=0;
for i = 1:n
b = b  + (X(i,:).' * Y(i,:)).^2 - 2.* (X(i,:).' * Y(i,:)) .*S + S.^2;
end
b=sum(b(:))

我们知道(a * b) ^ 2和a^2 * b^2:是一样的

X2 = X.^2;
Y2 = Y.^2;
S2 = S.^2;
b=0;
for i = 1:n
b = b  + (X2(i,:).' * Y2(i,:)) - 2.* (X(i,:).' * Y(i,:)) .*S + S2;
end
b=sum(b(:))

现在我们可以分别计算每个项：

b = sum(sum(X2.' * Y2 - 2 * (X.' * Y ) .* S + n * S2));

以下是Octave中的一个测试结果，该测试比较了我的方法和@AndrasDeak提供的另外两种方法，以及针对大小为500*500:的输入的原始基于循环的解决方案

===rahnema1 (B)===
Elapsed time is 0.0984299 seconds.
===Andras Deak (B2)===
Elapsed time is 7.86407 seconds.
===Andras Deak (B3)===
Elapsed time is 2.99158 seconds.
===Loop solution===
Elapsed time is 2.20357 seconds

n=500;
X= rand(n);
Y= rand(n);
S= rand(n);
disp('===rahnema1 (B)===')
tic
X2 = X.^2;
Y2 = Y.^2;
S2 = S.^2;
b=sum(sum(X2.' * Y2 - 2 * (X.' * Y ) .* S + n * S2));
toc
disp('===Andras Deak (B2)===')
tic
b2 = sum(reshape((permute(reshape(X, [n, 1, n]).*Y, [3,2,1]) - S).^2, 1, []));
toc
disp('===Andras Deak (B3)===')
tic
b3 = sum(reshape((reshape(X, [n, 1, n]).*Y - reshape(S.', [1, n, n])).^2, 1, []));
toc
tic
b=0;
for i = 1:n
b = b  + sum(sum((X(i,:)' * Y(i,:) - S).^2));
end
toc

Answer 2

您可能无法节省时间复杂性，但您可以使用矢量化来摆脱循环，并尽可能多地使用低级代码和缓存。它是否真的更快取决于你的尺寸，所以你需要做一些时间测试，看看这是否值得：

% dummy data
n = 3;
X = rand(n);
Y = rand(n);
S = rand(n);
% vectorize
b2 = sum(reshape((permute(reshape(X, [n, 1, n]).*Y, [3,2,1]) - S).^2, 1, []));
% check
b - b2 % close to machine epsilon i.e. zero

发生的情况是，我们在其中一个数组中插入一个新的单例维度，最终得到一个大小为[n, 1, n]的数组，而另一个数组为[n, n]，后者隐含地与[n, n, 1]相同。重叠的第一个索引对应于循环中的i，其余两个索引对应着每个i的二元乘积的矩阵索引。然后，我们对索引进行排列，以便将"i"索引放在最后，这样我们就可以使用(隐式(大小为[n, n, 1]的S再次广播结果。然后我们有一个大小为[n, n, n]的矩阵，其中前两个索引是原始索引中的矩阵索引，最后一个对应于i。然后，我们只需要取平方并对每个项求和(而不是求两次和，我将数组重新整形为一行并求和一次(。

上面的一个微小变化转置了S，而不是可能更快的3d阵列(同样，你应该计时(：

b3 = sum(reshape((reshape(X, [n, 1, n]).*Y - reshape(S.', [1, n, n])).^2, 1, []));

就性能而言，reshape是免费的(它只重新解释数据，不复制(，但permute/transpose在复制数据时通常会导致性能下降。