如何在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

提前感谢您。

- Ignacio

---

如果这些答案解决了您的问题，请考虑标记最适合您需求的答案为已接受:)（我相信Luis的解决方案不仅更简单，而且更快速和更节省内存。） - Andras Deak -- Слава Україні

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

如果你将向量视为行/列矩阵，那么这是一系列的矩阵乘积。由于第一个和最后一个维度（通过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'

- Andras Deak -- Слава Україні

---

太棒了！运行得很好！ 谢谢你 - Ignacio

Answer 2

你可以通过结合律和转置乘积性质大大减少操作数量：

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来提高一些速度。

- Luis Mendo

---

@AndrasDeak 谢谢！我相信这种方法会比你的方法更快。我总是尽量避免计算矩阵，只为保留对角线 :-) - Luis Mendo

我非常确定它会的 :) 更不用说节省内存了。 - Andras Deak -- Слава Україні

是的，它确实更加高效。 - Ignacio