Matlab 中的块矩阵内积
Block matrix inner products in Matlab
我一直在使用以下自定义函数执行向量与矩阵的乘法,其中向量的每个元素与 (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
由于我经常调用它们,不幸的是,它们显着降低了我代码的 运行 速度。我想知道是否有上述操作的更有效(可能是矢量化)版本。
在这两种情况下,您都可以取消 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
:
来实现 faster/simpler
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
以下是矢量化解决方案:
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
我一直在使用以下自定义函数执行向量与矩阵的乘法,其中向量的每个元素与 (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
由于我经常调用它们,不幸的是,它们显着降低了我代码的 运行 速度。我想知道是否有上述操作的更有效(可能是矢量化)版本。
在这两种情况下,您都可以取消 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
以下是矢量化解决方案:
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