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

如果您将向量视为行/列矩阵,则是一系列矩阵产品.并且由于第一个和最后一个维度(通过b'和b)都是1,因此每个维度都得到一个标量结果i.如果您将此标量视为1 x 1矩阵,则可以使用自己的跟踪替换它.并且,矩阵可以在跟踪下循环置换:

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

简单的Bi地方代表第一i排B.但是我们现在在跟踪内部的内容再次是标量,因此我们可以删除跟踪.所以i你需要的每一个

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

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

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