Matlab中的块矩阵内积答案

【问题标题】：Block matrix inner products in MatlabMatlab中的块矩阵内积
【发布时间】：2019-02-20 09:45:12
【问题描述】：

我一直在使用以下自定义函数来执行向量与矩阵的乘法运算，其中向量的每个元素乘以 (3xN)x(3) 矩阵内的 3x3 块：

function [B] = BlockScalar(v,A)

 N=size(v,2);
 B=zeros(3*N,3);

 for i=1:N
     B(3*i-2:3*i,:) = v(i).*A(3*i-2:3*i,:);
 end

end

类似地，当我想将 3x3 矩阵的集合乘以 3x3 向量的集合时，我使用以下代码

function [B] = BlockMatrix(A,u)

 N=size(u,2);
 B=zeros(N,3);

 for i=1:N
     B(i,:) = A(3*i-2:3*i,:)*u(:,i);
 end

end

由于我经常调用它们，不幸的是，它们显着减慢了我的代码运行速度。我想知道上述操作是否有更高效（也许是矢量化）的版本。

【问题讨论】：

标签： matlab matrix linear-algebra

【解决方案1】：

在这两种情况下，您都可以取消 for 循环（尽管未经测试，我无法确认这是否一定会加快您的计算速度）。

对于第一个功能，你可以这样做：

function [B] = BlockScalar(v,A)
% We create a vector N = [1,1,1,2,2,2,3,3,3,...,N,N,N]
N=ceil((1:size(A,1))/3); 

% Use N to index v, and let matlab do the expansion
B = v(N).*A;

end

对于第二个函数，我们可以制作一个块对角矩阵。

function [B] = BlockMatrix(A,u)

 N=size(u,2)*3;
 % We use a little meshgrid+sparse magic to convert A to a block matrix
 [X,Y] = meshgrid(1:N,1:3);

% Use sparse matrices to speed up multiplication and save space
 B = reshape(sparse(Y+(ceil((1:N)/3)-1)*3,X,A) * (u(:)),3,size(u,2))';

end

请注意，如果您能够访问单个 3x3 矩阵，则可以通过使用本机 blkdiag 来更快/更简单：

function [B] = BlockMatrix(a,b,c,d,...,u)
% Where A = [a;b;c;d;...];

% We make one of the input matrices sparse to make the whole block matrix sparse
% This saves memory and potentially speeds up multiplication by a lot
% For small enough values of N, however, using sparse may slow things down.
reshape(blkdiag(sparse(a),b,c,d,...) * (u(:)),3,size(u,2))';

end

【讨论】：

谢谢！对于第一个函数，我不清楚 i 是什么，因为没有 for 循环。
另外，第二个函数不适用于 3x3 矩阵和 3 分量列向量。
对不起，我在这两种情况下都犯了错误的复制和粘贴。我相信已经解决了这些问题。你能再试一次吗？
错误使用重塑要重塑元素的数量不能改变。
好的，我想我也应该解决这个问题。

【解决方案2】：

这里是矢量化解决方案：

function [B] = BlockScalar(v,A)
    N = size(v,2);
    B = reshape(reshape(A,3,N,3) .* v, 3*N, 3);
end


function [B] = BlockMatrix(A,u)
    N = size(u,2);
    A_r = reshape(A,3,N,3);
    B = (A_r(:,:,1) .* u(1,:) + A_r(:,:,2) .* u(2,:) + A_r(:,:,3) .* u(3,:)).';
end

function [B] = BlockMatrix(A,u)
    N = size(u,2);
    B = sum(reshape(A,3,N,3) .* permute(u, [3 2 1]) ,3).';
end

【讨论】：