我想知道是否有任何方法可以更有效地执行以下矩阵乘法，例如对其进行矢量化：
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
提前致谢
【问题讨论】：
标签： arrays matlab matrix vectorization matrix-multiplication
我建议将矩阵乘法重新排序，让 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
是一系列矩阵乘积，如果您将向量视为行/列矩阵。由于第一个和最后一个维度（通过 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'
【讨论】：
您可以使用associativity 和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（通过 transpose-product 属性）；而后者可以为所有i 计算为(B*a*b)'*B*A*b。
如果您的矩阵包含实数，您也许可以加快删除abs 的速度。
【讨论】：