两个数组的对应二维切片的乘法和数组切片的反转

Question

我有两个数组A 和B 具有相同的维度1000 x 3 x 20 x 20。我想生成维度为3 x 3 x 20 x 20 的第三个数组C，这将是A 和B 的对应切片矩阵相乘的结果，即C(:,:,i,j) = A(:,:,i,j)'*B(:,:,i,j)。然后我需要通过反转相应的3 x 3 矩阵，即D(:,:,i,j) = inv(C(:,:,i,j))，将数组C 转换为新数组D。同样，很清楚如何使用循环来做到这一点。有没有办法避免循环遍历 400 项目？

编辑：比较不同解决方案性能的基准代码是 -

%// Inputs
n1 = 50;
n2 = 200;
A = rand(n1,3,n2,n2);
B = rand(n1,3,n2,n2);

%// A. CPU loopy code
tic
C = zeros(3,3,n2,n2);
for ii = 1:n2
    for jj = 1:n2
        C(:,:,ii,jj) = A(:,:,ii,jj)'*B(:,:,ii,jj); %//'
    end
end
toc

%// B. Vectorized code (using squeeze)
tic
C1 = squeeze(sum(bsxfun(@times,permute(A,[2 1 5 3 4]),permute(B,[5 1 2 3 4])),2));
toc

%// C. Vectorized code (avoiding squeeze)
tic
C2 = sum(bsxfun(@times,permute(A,[2 5 3 4 1]),permute(B,[5 2 3 4 1])),5);
toc

%// D. GPU vectorized code
tic
A = gpuArray(A);
B = gpuArray(B);
C3 = sum(bsxfun(@times,permute(A,[2 5 3 4 1]),permute(B,[5 2 3 4 1])),5);
C3 = gather(C3);
toc

运行时结果 -

Elapsed time is 0.287511 seconds.
Elapsed time is 0.250663 seconds.
Elapsed time is 0.337628 seconds.
Elapsed time is 1.259207 seconds.

Answer 1

代码

%// Part - 1
C = sum(bsxfun(@times,permute(A,[2 5 3 4 1]),permute(B,[5 2 3 4 1])),5);

%// Part - 2: Use MATLAB file-exchange tool multinv
D = multinv(C);

multinv 的功能码可在here 使用，它声称非常高效。

第一部分，你也可以试试这个-

C = squeeze(sum(bsxfun(@times,permute(A,[2 1 5 3 4]),permute(B,[5 1 2 3 4])),2));

这似乎是重新排列元素不像上面代码中提到的那样“破坏性”，但缺点是需要squeeze，这可能会减慢它的速度。我会把它留给您，并鼓励您进行基准测试并选择更好的。

---

为什么是bsxfun + GPU？

我增加了循环限制，因为这可能是循环代码和矢量化代码之间的真正测试。所以，这是第 1 部分的修改代码 -

%// Inputs
n1 = 50;
n2 = 200;
A = rand(n1,3,n2,n2);
B = rand(n1,3,n2,n2);

%// A. CPU loopy code
tic
C = zeros(3,3,n2,n2);
for ii = 1:n2
    for jj = 1:n2
        C(:,:,ii,jj) = A(:,:,ii,jj)'*B(:,:,ii,jj); %//'
    end
end
toc

%// B. GPU vectorized code
tic
A = gpuArray(A);
B = gpuArray(B);
C1 = sum(bsxfun(@times,permute(A,[2 5 3 4 1]),permute(B,[5 2 3 4 1])),5);
C1 = gather(C1);
toc

我系统上的运行时结果是 -

Elapsed time is 0.310056 seconds.
Elapsed time is 0.172499 seconds.

所以，你看！