Element-wise matrix multiplication for multi-dimensional array

Question

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

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

Answer 1

选项1：从MATLAB调用numpy

假设您的系统是根据文档设置的，并且已经安装了numpy软件包，则可以执行以下操作（在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可以找到这里。

选项2：本机MATLAB

这里最重要的部分是正确排列，因此我们需要跟踪尺寸。我们将使用以下顺序：

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

现在，我将说明如何获得正确的置换顺序（以的示例为例A）：einsum预期尺寸顺序为mki，根据我们的编号为5 3 1。这告诉我们，1 ^日的尺寸A需要是5 ^次，2 ^个需要是3 ^次和3 ^个需要是1 ^日（简称1->5, 2->3, 3->1）。这也意味着“无源尺寸”（意味着没有原始尺寸的尺寸；在这种情况下为2 4 6）应为单件。使用ipermute它真的很简单：

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

在上面的示例中，1->5意味着我们5首先进行编写，而其他两个维度也是如此（产生[5,3,1]）。然后，我们只需在末尾添加单例（2,4,6）即可[5,3,1,2,4,6]。最后：

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

（请参阅有关sum帖子底部的注释。）

如@AndrasDeak所述，在MATLAB中执行此操作的另一种方法是：

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

另请参阅：squeeze，reshape，permute，ipermute，shiftdim。

这是一个完整的示例，显示了测试这些方法是否等效的示例：

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]);

运行以上操作，我们得出的结果确实是等效的：

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]));

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

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

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

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)

在我的系统上，这导致：

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

这意味着最好使用^第二种方法（至少对于默认输入大小而言）。