Поэлементное умножение матриц для многомерного массива - PullRequest
Новая
       94

Поэлементное умножение матриц для многомерного массива

2 голосов
/ 30 апреля 2019

Я хочу реализовать компонентное умножение матриц в MATLAB, что можно сделать с помощью numpy.einsum в Python, как показано ниже: 

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

Пока я использую цикл fori, j и n, но я не хочу, по крайней мере, цикл for n.

1 Ответ

8 голосов
/ 30 апреля 2019

Вариант 1: вызов numy из MATLAB

Если ваша система настроена в соответствии с документацией , и у вас установлен пакет 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: Native MATLAB

Здесь наиболее важной частью является получение перестановок.правильно, поэтому мы должны отслеживать наши измерения.Мы будем использовать следующий порядок: 

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

Теперь я объясню, как я получил правильный порядок перестановок (давайте рассмотрим пример A): einsum ожидает, что порядок измерения равенбыть mki, что согласно нашей нумерации 5 3 1.Это говорит нам о том, что 1 st измерение A должно быть 5 th , 2 nd должно быть 3 rd и 3 rd должно быть 1 st (короче 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].Наконец: 

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 внизу поста.)

Другой способ сделать это в MATLAB, , как упомянуто @ AndrasDeak , является следующим: 

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, shiftdim.

Вот полный пример, который показывает, проверяет, являются ли эти методы эквивалентными: 

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

Обратите внимание, что эти два метода вызова 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 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

В моей системе это приводит к: 

>> q55913093_benchmark
ans =
    1.3964    0.1864    0.2428

Это означает, что метод 2 nd предпочтителен (по крайней мере для входных размеров по умолчанию).

...