Умножьте 3D-матрицу на 2D-матрицу - PullRequest
Новая
       55

Умножьте 3D-матрицу на 2D-матрицу

25 голосов
/ 17 ноября 2009

Предположим, у меня есть AxBxC матрица X и BxD матрица Y.

Существует ли нецикличный метод, с помощью которого я могу умножить каждую из C AxB матриц на Y?

Ответы [ 9 ]

16 голосов
/ 17 ноября 2009

В качестве личного предпочтения мне нравится, чтобы мой код был максимально лаконичным и читабельным.

Вот что я бы сделал, хотя это не соответствует вашему требованию «без петель»: 

for m = 1:C

    Z(:,:,m) = X(:,:,m)*Y;

end

В результате получается A x D x C матрица Z .

И, конечно, вы всегда можете заранее выделить Z, чтобы ускорить процесс, используя Z = zeros(A,D,C);.

15 голосов
/ 17 ноября 2009

Вы можете сделать это в одной строке, используя функции NUM2CELL , чтобы разбить матрицу X на массив ячеек и CELLFUN для работы с ячейками: 

Z = cellfun(@(x) x*Y,num2cell(X,[1 2]),'UniformOutput',false);

Результат Z представляет собой массив ячеек 1-на-C , где каждая ячейка содержит матрицу A-by-D . Если вы хотите, чтобы Z была матрицей A-by-D-by-C , вы можете использовать функцию CAT : 

Z = cat(3,Z{:});


ПРИМЕЧАНИЕ: Мое старое решение использовало MAT2CELL вместо NUM2CELL , что было не так кратко: 

[A,B,C] = size(X);
Z = cellfun(@(x) x*Y,mat2cell(X,A,B,ones(1,C)),'UniformOutput',false);
8 голосов
/ 17 ноября 2009

Вот решение в одну строку (два, если вы хотите разделить на 3-е измерение): 

A = 2;
B = 3;
C = 4;
D = 5;

X = rand(A,B,C);
Y = rand(B,D);

%# calculate result in one big matrix
Z = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])' * Y;

%'# split into third dimension
Z = permute(reshape(Z',[D A C]),[2 1 3]);

Следовательно, теперь: Z(:,:,i) содержит результат X(:,:,i) * Y

Пояснение:

Выше может показаться запутанным, но идея проста. Сначала я беру третье измерение X и выполняю вертикальную конкатенацию вдоль первого затемнения: 

XX = cat(1, X(:,:,1), X(:,:,2), ..., X(:,:,C))

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

ZZ = XX * Y;

Наконец я разделил его на третье измерение: 

Z(:,:,1) = ZZ(1:2, :);
Z(:,:,2) = ZZ(3:4, :);
Z(:,:,3) = ZZ(5:6, :);
Z(:,:,4) = ZZ(7:8, :);

Таким образом, вы можете видеть, что требуется только одно умножение матриц, но вам нужно изменить матрицу до и после.

5 голосов
/ 22 октября 2015

Я подхожу к точно такой же проблеме, ориентируясь на наиболее эффективный метод. Есть примерно три подхода, которые я вижу, за исключением использования внешних библиотек (то есть mtimesx ):

  1. Цикл срезов трехмерной матрицы
  2. Колдовско-перестановочное волшебство
  3. умножение на ячейки

Недавно я сравнил все три метода, чтобы увидеть, какой из них был самым быстрым. Моя интуиция заключалась в том, что (2) будет победителем. Вот код: 

% generate data
A = 20;
B = 30;
C = 40;
D = 50;

X = rand(A,B,C);
Y = rand(B,D);

% ------ Approach 1: Loop (via @Zaid)
tic
Z1 = zeros(A,D,C);
for m = 1:C
    Z1(:,:,m) = X(:,:,m)*Y;
end
toc

% ------ Approach 2: Reshape+Permute (via @Amro)
tic
Z2 = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])' * Y;
Z2 = permute(reshape(Z2',[D A C]),[2 1 3]);
toc


% ------ Approach 3: cellfun (via @gnovice)
tic
Z3 = cellfun(@(x) x*Y,num2cell(X,[1 2]),'UniformOutput',false);
Z3 = cat(3,Z3{:});
toc

Все три подхода дали одинаковый результат (фу!), Но, что удивительно, цикл оказался самым быстрым: 

Elapsed time is 0.000418 seconds.
Elapsed time is 0.000887 seconds.
Elapsed time is 0.001841 seconds.

Обратите внимание, что время может сильно варьироваться от одного испытания к другому, и иногда (2) получается самое медленное. Эти различия становятся более существенными с большими данными. Но с намного большими данными, (3) удары (2). Метод цикла все еще лучший. 

% pretty big data...
A = 200;
B = 300;
C = 400;
D = 500;
Elapsed time is 0.373831 seconds.
Elapsed time is 0.638041 seconds.
Elapsed time is 0.724581 seconds.

% even bigger....
A = 200;
B = 200;
C = 400;
D = 5000;
Elapsed time is 4.314076 seconds.
Elapsed time is 11.553289 seconds.
Elapsed time is 5.233725 seconds.

Но метод цикла может быть медленнее, чем (2), если зацикленное измерение намного больше, чем другие. 

A = 2;
B = 3;
C = 400000;
D = 5;
Elapsed time is 0.780933 seconds.
Elapsed time is 0.073189 seconds.
Elapsed time is 2.590697 seconds.

Итак (2) выигрывает с большим фактором, в этом (может быть, крайнем) случае. Может не быть подхода, который был бы оптимальным во всех случаях, но цикл все еще довольно хорош, а во многих случаях лучше. Это также лучше с точки зрения читабельности. Отойди!

1 голос
/ 12 августа 2017

Я настоятельно рекомендую вам использовать MMX toolbox из Matlab. Он может максимально быстро умножить n-мерные матрицы.

Преимущества MMX :

  1. легко использовать.
  2. Умножить n-мерные матрицы (на самом деле он может умножать массивы двумерных матриц)
  3. Выполняет другие матричные операции (транспонирование, квадратичное умножение, декомпозиция Chol и другие)
  4. Он использует C компилятор и многопоточность вычисления для ускорения.

Для этой проблемы вам просто нужно написать следующую команду: 

C=mmx('mul',X,Y);

вот эталон для всех возможных методов. Более подробно см. Этот вопрос . 

    1.6571 # FOR-loop
    4.3110 # ARRAYFUN
    3.3731 # NUM2CELL/FOR-loop/CELL2MAT
    2.9820 # NUM2CELL/CELLFUN/CELL2MAT
    0.0244 # Loop Unrolling
    0.0221 # MMX toolbox  <===================
1 голос
/ 09 января 2014

Чтобы ответить на вопрос, и для удобочитаемости, см .:

  • ndmult , ajuanpi (Хуан Пабло Карбахал), 2013, GNU GPL

Input

  • 2 массива
  • тусклый

Пример 

 nT = 100;
 t = 2*pi*linspace (0,1,nT)’;

 # 2 experiments measuring 3 signals at nT timestamps
 signals = zeros(nT,3,2);
 signals(:,:,1) = [sin(2*t) cos(2*t) sin(4*t).^2];
 signals(:,:,2) = [sin(2*t+pi/4) cos(2*t+pi/4) sin(4*t+pi/6).^2];

 sT(:,:,1) = signals(:,:,1)’;
 sT(:,:,2) = signals(:,:,2)’;
   G = ndmult (signals,sT,[1 2]);

Источник

Оригинальный источник. Я добавил встроенные комментарии. 

function M = ndmult (A,B,dim)
  dA = dim(1);
  dB = dim(2);

  # reshape A into 2d
  sA = size (A);
  nA = length (sA);
  perA = [1:(dA-1) (dA+1):(nA-1) nA dA](1:nA);
  Ap = permute (A, perA);
  Ap = reshape (Ap, prod (sA(perA(1:end-1))), sA(perA(end)));

  # reshape B into 2d
  sB = size (B);
  nB = length (sB);
  perB = [dB 1:(dB-1) (dB+1):(nB-1) nB](1:nB);
  Bp = permute (B, perB);
  Bp = reshape (Bp, sB(perB(1)), prod (sB(perB(2:end))));

  # multiply
  M = Ap * Bp;

  # reshape back to original format
  s = [sA(perA(1:end-1)) sB(perB(2:end))];
  M = squeeze (reshape (M, s));
endfunction
1 голос
/ 17 ноября 2009

Неа. Есть несколько способов, но это всегда происходит в цикле, прямом или косвенном.

Просто чтобы порадовать мое любопытство, зачем ты вообще этого хочешь?

0 голосов
/ 17 ноября 2009

Вы можете «развернуть» цикл, т.е. выписать все умножения последовательно, которые будут происходить в цикле

0 голосов
/ 17 ноября 2009

Я бы подумал, что рекурсия, но это единственный другой нецикличный метод, который вы можете сделать

...