多维数组的逐元素矩阵乘法

多维数组的逐元素矩阵乘法

本文介绍了多维数组的逐元素矩阵乘法的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我想在MATLAB中实现基于组件的矩阵乘法,可以使用 numpy.einsum 在Python中如下所示:

I want to realize component-wise matrix multiplication in MATLAB, which can be done using numpy.einsum in Python as below:

import numpy as np
M = 2
N = 4
I = 2000
J = 300

A = np.random.randn(M, M, I)
B = np.random.randn(M, M, N, J, I)
C = np.random.randn(M, J, I)

# using einsum
D = np.einsum('mki, klnji, lji -> mnji', A, B, C)

# naive for-loop
E = np.zeros(M, N, J, I)
for i in range(I):
    for j in range(J):
        for n in range(N):
            E[:,n,j,i] = B[:,:,i] @ A[:,:,n,j,i] @ C[:,j,i]

print(np.sum(np.abs(D-E))) # expected small enough

到目前为止,我使用ijn的for循环,但是我不想,至少是n的for循环.

So far I use for-loop of i, j, and n, but I don't want to, at least for-loop of n.

推荐答案

选项1:从MATLAB调用numpy

假设您的系统已设置文档,并且您已经安装了numpy软件包,则可以执行以下操作(在MATLAB中):

Option 1: Calling numpy from MATLAB

Assuming your system is set up according to the documentation, and you have the numpy package installed, you could do (in MATLAB):

np = py.importlib.import_module('numpy');

M = 2;
N = 4;
I = 2000;
J = 300;

A = matpy.mat2nparray( randn(M, M, I) );
B = matpy.mat2nparray( randn(M, M, N, J, I) );
C = matpy.mat2nparray( randn(M, J, I) );

D = matpy.nparray2mat( np.einsum('mki, klnji, lji -> mnji', A, B, C) );

matpy的位置此处.

最重要的部分是正确排列,因此我们需要跟踪尺寸.我们将使用以下顺序:

Here the most important part is to get the permutations right, so we need to keep track of our dimensions. We'll be using the following order:

I(1) J(2) K(3) L(4) M(5) N(6)

现在,我将说明如何获得正确的置换顺序(以A为例):einsum期望尺寸顺序为mki,根据我们的编号,其为5 3 1.这告诉我们A的第1 维需要为第5 ,第2 必须为3 和3 必须是1 (简称1->5, 2->3, 3->1).这也意味着无源尺寸"(意味着没有原始尺寸的尺寸;在这种情况下为2 4 6)应为单件.使用ipermute这真的很容易编写:

Now, I'll explain how I got the correct permute order (let's take the example of A): einsum expects the dimension order to be mki, which according to our numbering is 5 3 1. This tells us that the 1 dimension of A needs to be the 5, the 2 needs to be 3 and the 3 needs to be 1 (in short 1->5, 2->3, 3->1). This also means that the "sourceless dimensions" (meaning those that have no original dimensions becoming them; in this case 2 4 6) should be singleton. Using ipermute this is really simple to write:

pA = ipermute(A, [5,3,1,2,4,6]);

在上面的示例中,1->5表示我们首先编写5,其他两个维度也是如此(产生[5,3,1]).然后,我们只需在末尾添加单例(2,4,6)即可得到[5,3,1,2,4,6].最后:

In the above example, 1->5 means we write 5 first, and the same goes for the other two dimensions (yielding [5,3,1]). Then we just add the singletons (2,4,6) at the end to get [5,3,1,2,4,6]. Finally:

A = randn(M, M, I);
B = randn(M, M, N, J, I);
C = randn(M, J, I);

% Reference dim order: I(1) J(2) K(3) L(4) M(5) N(6)
pA = ipermute(A, [5,3,1,2,4,6]); % 1->5, 2->3, 3->1; 2nd, 4th & 6th are singletons
pB = ipermute(B, [3,4,6,2,1,5]); % 1->3, 2->4, 3->6, 4->2, 5->1; 5th is singleton
pC = ipermute(C, [4,2,1,3,5,6]); % 1->4, 2->2, 3->1; 3rd, 5th & 6th are singletons

pD = sum( ...
  permute(pA .* pB .* pC, [5,6,2,1,3,4]), ... 1->5, 2->6, 3->2, 4->1; 3rd & 4th are singletons
  [5,6]);

(请参阅帖子底部有关sum的注释.)

(see note regarding sum at the bottom of the post.)

@AndrasDeak 提到的在MATLAB中执行此操作的另一种方法是: /p>

Another way to do it in MATLAB, as mentioned by @AndrasDeak, is the following:

rD = squeeze(sum(reshape(A, [M, M, 1, 1, 1, I]) .* ...
                 reshape(B, [1, M, M, N, J, I]) .* ...
... % same as:   reshape(B, [1, size(B)]) .* ...
... % same as:   shiftdim(B,-1) .* ...
                 reshape(C, [1, 1, M, 1, J, I]), [2, 3]));

另请参阅: squeeze reshape permute ipermute .

下面是一个完整的示例,该示例显示了测试这些方法是否等效的方法:

Here's a full example that shows that tests whether these methods are equivalent:

function q55913093
M = 2;
N = 4;
I = 2000;
J = 300;

mA = randn(M, M, I);
mB = randn(M, M, N, J, I);
mC = randn(M, J, I);

%% Option 1 - using numpy:
np = py.importlib.import_module('numpy');

A = matpy.mat2nparray( mA );
B = matpy.mat2nparray( mB );
C = matpy.mat2nparray( mC );

D = matpy.nparray2mat( np.einsum('mki, klnji, lji -> mnji', A, B, C) );

%% Option 2 - native MATLAB:
%%% Reference dim order: I(1) J(2) K(3) L(4) M(5) N(6)

pA = ipermute(mA, [5,3,1,2,4,6]); % 1->5, 2->3, 3->1; 2nd, 4th & 6th are singletons
pB = ipermute(mB, [3,4,6,2,1,5]); % 1->3, 2->4, 3->6, 4->2, 5->1; 5th is singleton
pC = ipermute(mC, [4,2,1,3,5,6]); % 1->4, 2->2, 3->1; 3rd, 5th & 6th are singletons

pD = sum( permute( ...
  pA .* pB .* pC, [5,6,2,1,3,4]), ... % 1->5, 2->6, 3->2, 4->1; 3rd & 4th are singletons
  [5,6]);

rD = squeeze(sum(reshape(mA, [M, M, 1, 1, 1, I]) .* ...
                 reshape(mB, [1, M, M, N, J, I]) .* ...
                 reshape(mC, [1, 1, M, 1, J, I]), [2, 3]));

%% Comparisons:
sum(abs(pD-D), 'all')
isequal(pD,rD)

运行上述操作,我们得出的结果确实是等效的:

Running the above we get that the results are indeed equivalent:

>> q55913093
ans =
   2.1816e-10
ans =
  logical
   1

请注意,这两种调用sum的方法是在最近的发行版中引入的,因此,如果您的MATLAB相对较旧,则可能需要替换它们:

Note that these two methods of calling sum were introduced in recent releases, so you might need to replace them if your MATLAB is relatively old:

S = sum(A,'all')   % can be replaced by ` sum(A(:)) `
S = sum(A,vecdim)  % can be replaced by ` sum( sum(A, dim1), dim2) `


根据评论中的要求,这是比较方法的基准:


As requested in the comments, here's a benchmark comparing the methods:

function t = q55913093_benchmark(M,N,I,J)

if nargin == 0
  M = 2;
  N = 4;
  I = 2000;
  J = 300;
end

% Define the arrays in MATLAB
mA = randn(M, M, I);
mB = randn(M, M, N, J, I);
mC = randn(M, J, I);

% Define the arrays in numpy
np = py.importlib.import_module('numpy');
pA = matpy.mat2nparray( mA );
pB = matpy.mat2nparray( mB );
pC = matpy.mat2nparray( mC );

% Test for equivalence
D = cat(5, M1(), M2(), M3());
assert( sum(abs(D(:,:,:,:,1) - D(:,:,:,:,2)), 'all') < 1E-8 );
assert( isequal (D(:,:,:,:,2), D(:,:,:,:,3)));

% Time
t = [ timeit(@M1,1), timeit(@M2,1), timeit(@M3,1)];

function out = M1()
  out = matpy.nparray2mat( np.einsum('mki, klnji, lji -> mnji', pA, pB, pC) );
end

function out = M2()
  out = permute( ...
          sum( ...
            ipermute(mA, [5,3,1,2,4,6]) .* ...
            ipermute(mB, [3,4,6,2,1,5]) .* ...
            ipermute(mC, [4,2,1,3,5,6]), [3,4]...
          ), [5,6,2,1,3,4]...
        );
end

function out = M3()
out = squeeze(sum(reshape(mA, [M, M, 1, 1, 1, I]) .* ...
                  reshape(mB, [1, M, M, N, J, I]) .* ...
                  reshape(mC, [1, 1, M, 1, J, I]), [2, 3]));
end

end

在我的系统上,结果为:

On my system this results in:

>> q55913093_benchmark
ans =
    1.3964    0.1864    0.2428

这意味着最好使用2 方法(至少对于默认输入大小而言).

Which means that the 2 method is preferable (at least for the default input sizes).

这篇关于多维数组的逐元素矩阵乘法的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!

08-19 23:41