矢量化和嵌套矩阵乘法

Question

这是原始代码:

K = zeros(N*N)
for a=1:N
    for i=1:I
        for j=1:J
            M = kron(X(:,:,a).',Y(:,:,a,i,j));

            %A function that essentially adds M to K. 
        end
    end
end

目标是向量化kroniker乘法调用.我的直觉是将X和Y视为矩阵的容器(作为参考,X和Y的切片被馈送到kron是7x7阶的矩阵).在此容器方案下,X显示为1-D容器,Y显示为3-D容器.我的下一个猜测是将Y重塑为2-D容器或更好的1-D容器,然后进行X和Y的元素乘法.问题是:如何以保留M和M的轨迹的方式重塑matlab甚至可以在这个容器构思中处理这个想法,还是需要进一步重塑容器以进一步暴露内部矩阵元素？

Answer 1

方法#1:使用6D置换的矩阵乘法

% Get sizes
[m1,m2,~] =  size(X);
[n1,n2,N,n4,n5] =  size(Y);

% Lose the third dim from X and Y with matrix-multiplication
parte1 = reshape(permute(Y,[1,2,4,5,3]),[],N)*reshape(X,[],N).';

% Rearrange the leftover dims to bring kron format
parte2 = reshape(parte1,[n1,n2,I,J,m1,m2]);

% Lose dims correspinding to last two dims coming in from Y corresponding
% to the iterative summation as suggested in the question
out = reshape(permute(sum(sum(parte2,3),4),[1,6,2,5,3,4]),m1*n1,m2*n2)

方法#2:简单的7D置换

% Get sizes
[m1,m2,~] =  size(X);
[n1,n2,N,n4,n5] =  size(Y);

% Perform kron format elementwise multiplication betwen the first two dims
% of X and Y, keeping the third dim aligned and "pushing out" leftover dims
% from Y to the back
mults = bsxfun(@times,permute(X,[4,2,5,1,3]),permute(Y,[1,6,2,7,3,4,5]));

% Lose the two dims with summation reduction for final output
out = sum(reshape(mults,m1*n1,m2*n2,[]),3);

验证

这是运行原始方法和建议方法的设置 -

% Setup inputs
X = rand(10,10,10);
Y = rand(10,10,10,10,10);

% Original approach
[n1,n2,N,I,J] =  size(Y);
K = zeros(100);
for a=1:N
    for i=1:I
        for j=1:J
            M = kron(X(:,:,a).',Y(:,:,a,i,j));
            K = K + M;
        end
    end
end

% Approach #1
[m1,m2,~] =  size(X);
[n1,n2,N,n4,n5] =  size(Y);
mults = bsxfun(@times,permute(X,[4,2,5,1,3]),permute(Y,[1,6,2,7,3,4,5]));
out1 = sum(reshape(mults,m1*n1,m2*n2,[]),3);

% Approach #2
[m1,m2,~] =  size(X);
[n1,n2,N,n4,n5] =  size(Y);
parte1 = reshape(permute(Y,[1,2,4,5,3]),[],N)*reshape(X,[],N).';
parte2 = reshape(parte1,[n1,n2,I,J,m1,m2]);
out2 = reshape(permute(sum(sum(parte2,3),4),[1,6,2,5,3,4]),m1*n1,m2*n2);

跑完后,我们看到了最大值.与原始方法的拟议方法的绝对偏差 -

>> error_app1 = max(abs(K(:)-out1(:)))
error_app1 =
   1.1369e-12
>> error_app2 = max(abs(K(:)-out2(:)))
error_app2 =
   1.1937e-12

价值观对我来说很好!

标杆

使用与验证相同的大数据集对这三种方法进行定时,我们得到类似的结果 -

----------------------------- With Loop
Elapsed time is 1.541443 seconds.
----------------------------- With BSXFUN
Elapsed time is 1.283935 seconds.
----------------------------- With MATRIX-MULTIPLICATION
Elapsed time is 0.164312 seconds.

看起来像矩阵乘法对这些大小的数据集做得相当好!