如何在 Matlab 中更有效地进行以下矩阵乘法？

Question

我想知道是否有任何方法可以更有效地执行以下矩阵乘法，例如对其进行矢量化：

m=1000;n=500;
a=zeros(n,1);
b=rand(n,1);
A=rand(m,n);
B=rand(m,m);
for i=1:n
    a(i)=b'*(A'*B(i,:)'*B(i,:)*A)*b;
end

提前致谢

Answer 1

我建议重新排序矩阵乘法，使 B 和 B' 位于两侧，并使用 diag 获得对角线元素：

m=1000;n=500;
a=zeros(n,1);
b=rand(n,1);
A=rand(m,n);
B=rand(m,m);

%original for comparison
for i=1:n
    a(i)=b'*(A'*B(i,:)'*B(i,:)*A)*b;
end

%new version
%a2=diag(B*A*b*b'*A'*B');
%a2=a2(1:n); %emulate original cut-off

%new new version
a2=diag(B(1:n,:)*A*b*b'*A'*B(1:n,:)');

%compare difference
max(abs(a-a2)) %absolute error
max(abs(a-a2))./max(abs(a)) %relative error

与rand()相比，绝对误差似乎很大：

>> max(abs(a-a2)) %absolute error

ans =

   8.1062e-06

但相对误差表明我们对机器精度是正确的：

>> max(abs(a-a2))./max(abs(a)) %relative error

ans =

   1.9627e-15

重新排列背后的数学原理是原始总和中的一个项，

b'*(A'*B(i,:)'*B(i,:)*A)*b

是一系列矩阵乘积，如果把向量想成row/column矩阵。由于第一个和最后一个维度（通过 b' 和 b）都是 1，所以每个 i 都会得到一个标量结果。如果您将此标量视为 1 x 1 矩阵，则可以用它自己的迹代替它。并且，矩阵可以在迹线下循环置换：

Tr (b' * A * Bi' * Bi * A * b) = Tr (Bi * A * b * b' * A * Bi')

为简单起见，Bi 代表 B 的第 i 行。但是我们现在在跟踪中拥有的是又是一个标量，因此我们可以删除跟踪。所以对于每个 i 你需要

a(i)=B(i,:) * A * b * b' * A * B(i,:)';

这显然是矩阵的 (i,i)（对角线）分量

B * A * b * b' * A * B'

Answer 2

使用associativity and the transpose-product property可以大大减少操作次数:

t = B*A*b;
a = abs(t).^2;

这是有效的，因为原始表达式 b'*(A'*B(i,:)'*B(i,:)*A)*b 等于 b'*A'*B(i,:)'*B(i,:)*A*b（根据结合律），等于 (B(i,:)*a*b)'*B(i,:)*A*b（根据转置积属性）；后者可以为所有 i 计算为 (B*a*b)'*B*A*b.

如果您的矩阵包含实数，删除 abs.

或许可以加快速度

如何在 Matlab 中更有效地进行以下矩阵乘法？

How to do the following matrix multiplication more efficient in Matlab?

arrays

matlab

matrix

vectorization

matrix-multiplication