Умножить трехмерную матрицу с 2D-матрицей

Ответ 1

Вы можете сделать это в одной строке, используя функции 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);

Ответ 2

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

Вот что я бы сделал, хотя он не соответствует вашим требованиям "no-loops":

for m = 1:C

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

end

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

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

Ответ 3

Здесь однострочное решение (два, если вы хотите разбить на 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 и выполняю вертикальную конкатенацию вдоль первого dim:

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, :);

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

Ответ 4

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

• Прокрутка по срезам трехмерной матрицы
• repmat-and-permut wizardry
• умножение сотовой сети

Недавно я сравнил все три метода, чтобы увидеть, что было быстрее. Моя интуиция заключалась в том, что (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). Метод loop еще лучше.

% 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) выигрывает большой фактор, в этом (возможно, экстремальном) случае. Во всех случаях не может быть оптимального подхода, но цикл по-прежнему очень хорош и во многих случаях лучше. Это также лучше с точки зрения удобочитаемости. Пройдите!

Ответ 5

Неа. Существует несколько способов, но он всегда появляется в цикле, прямо или косвенно.

Просто, чтобы угодить моему любопытству, зачем вам это все равно?

Ответ 6

Чтобы ответить на вопрос и прочитать его, читайте:

• ndmult, by ajuanpi (Juan Pablo Carbajal), 2013, GNU GPL

Ввод

• 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

Ответ 7

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

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

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

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

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

здесь является эталоном для всех возможных методов. Для более подробной информации обратитесь к question.

    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  <===================

Ответ 8

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

Ответ 9

Вы можете "развернуть" цикл, т.е. выписать все последовательные последовательности, которые произойдут в цикле

Ответ 10

У меня есть аналогичный вопрос, но трехмерная матрица X AxBxC и двухмерная матрица Y CxD, и вы хотите получить AxBxD. Размеры:

A = 30 B = 70 C = 300 D = 100

Трехмерная матрица - это фиктивная переменная, которая принимает value =

1 в каждом измерении C в экземплярах AxB, если (...) (и сумма всех Cs = 300), различная для каждого C.

0 в противном случае

Двухмерная матрица - это временные данные.

Моя самая большая проблема связана с фиктивной переменной.

Ответ 11

Я хотел бы поделиться своим ответом на проблемы:

1) создание тензорного произведения двух тензоров (любой валентности);

2) сокращение двух тензоров вдоль любого измерения.

Вот мои подпрограммы для первого и второго задания:

1) тензорное произведение:

function [C] = tensor(A,B)
   C = squeeze( reshape( repmat(A(:), 1, numel(B)).*B(:).' , [size(A),size(B)] ) );
end

2) сокращение: Здесь A и B - тензоры, которые должны быть сжаты по размерам i и j соответственно. Длина этих размеров, конечно, должна быть одинаковой. Для этого нет проверки (это может затенить код), но кроме этого он работает хорошо.

   function [C] = tensorcontraction(A,B, i,j)
      sa = size(A);
      La = length(sa);
      ia = 1:La;
      ia(i) = [];
      ia = [ia i];

      sb = size(B);
      Lb = length(sb);
      ib = 1:Lb;
      ib(j) = [];
      ib = [j ib];

      % making the i-th dimension the last in A
      A1 = permute(A, ia);
      % making the j-th dimension the first in B
      B1 = permute(B, ib);

      % making both A and B 2D-matrices to make use of the
      % matrix multiplication along the second dimension of A
      % and the first dimension of B
      A2 = reshape(A1, [],sa(i));
      B2 = reshape(B1, sb(j),[]);

      % here the implicit implication that sa(i) == sb(j),
      % otherwise - crash
      C2 = A2*B2;

      % back to the original shape with the exception
      % of dimensions along which we've just contracted
      sa(i) = [];
      sb(j) = [];
      C = squeeze( reshape( C2, [sa,sb] ) );
   end

Есть ли критики?