Moltiplicare una matrice 3D con una matrice 2D

Question

Supponiamo di avere un AxBxC matrice X e un BXD matrice Y.

Esiste un metodo non loop che posso moltiplicare ciascuna delle C AxB matrici con Y?

Answer 1

È possibile farlo in una sola riga utilizzando le funzioni NUM2CELL per rompere il X matrice in una matrice di celle e CELLFUN di operare attraverso le celle:

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

Il Z risultato è un 1-by-C matrice di celle in cui ogni cella contiene un A-by-D matrice. Se si desidera Z di essere un A-by-D-by-C a matrice, è possibile utilizzare il funzione CAT :

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

---

Nota: La mia vecchia soluzione utilizzata MAT2CELL invece di NUM2CELL , che non era così succinta:

[A,B,C] = size(X);
Z = cellfun(@(x) x*Y,mat2cell(X,A,B,ones(1,C)),'UniformOutput',false);

Answer 2

Come una preferenza personale, mi piace il mio codice di essere il più sintetico e leggibile possibile.

Ecco quello che avrei fatto, anche se non soddisfa le vostre esigenze 'no-loop':

for m = 1:C

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

end

Ciò si traduce in un A x D x C matrice Z .

E, naturalmente, si può sempre pre-allocare Z per velocizzare le cose utilizzando Z = zeros(A,D,C);.

Answer 3

Ecco una soluzione one-line (due se si desidera dividere in 3 dimensioni):

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

Quindi ora: Z(:,:,i) contiene il risultato di X(:,:,i) * Y

---

Spiegazione:

È possibile che questo può sembrare confuso, ma l'idea è semplice. Per prima cosa ho Iniziare prendo la terza dimensione della X e fare una concatenazione verticale lungo la prima fioca:

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

... la difficoltà è stata che C è una variabile, quindi non si può generalizzare che espressione utilizzando cat o vertcat . Successivo moltiplichiamo questo Y:

ZZ = XX * Y;

Finalmente ho diviso di nuovo nella terza dimensione:

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

Così si può vedere che richiede solo una moltiplicazione di matrici, ma è necessario reshape la matrice prima e dopo.

Answer 4

sto avvicinando lo stesso problema esatto, con un occhio per il metodo più efficiente. Ci sono circa tre approcci che vedo in giro, a corto di utilizzare librerie esterne (ad esempio, mtimesx ):

1. Loop attraverso fette di matrice 3D
2. repmat-e-permute magia
3. cellfun moltiplicazione

Recentemente ho confrontato tutti e tre i metodi per vedere che era più veloce. La mia intuizione è che (2) sarebbe stato il vincitore. Ecco il codice:

% 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

Tutti e tre gli approcci hanno prodotto lo stesso risultato (uff!), Ma, a sorpresa, il ciclo è stato il più veloce:

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

Si noti che i tempi variano parecchio da un processo ad un altro, e talvolta (2) viene fuori il più lento. Queste differenze diventano più drammatico con i dati più grandi. Ma con molto di dati più grandi, (3) batte (2). Il metodo ciclo è ancora meglio.

% 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.

Ma il metodo ciclo possono essere più lento di (2), se la dimensione loop è molto più grande rispetto agli altri.

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) vince da un grande fattore, in questo caso (forse estremo). Non ci può essere un approccio che è ottimale in tutti i casi, ma il loop è ancora piuttosto buona, e meglio in molti casi. E 'anche meglio in termini di leggibilità. Loop via!

Answer 5

No. Ci sono diversi modi, ma si tratta sempre fuori in un ciclo, diretto o indiretto.

Proprio per compiacere la mia curiosità, perché si vuole che in ogni caso?

Answer 6

Per rispondere alla domanda, e per migliorare la leggibilità, si veda:

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

ingresso

• 2 array
• dim

Esempio

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

sorgente

Fonte originale. Ho aggiunto commenti in linea.

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

Answer 7

Mi raccomando di utilizzare il MMX toolbox di mATLAB. Si può moltiplicare le matrici n-dimensionali il più velocemente possibile.

I vantaggi di MMX sono:

1. È semplice da usare.
2. Moltiplicare matrici n-dimensionale (realmente può moltiplicare le matrici di matrici 2-D)
3. Si esegue altre operazioni di della matrice (trasposizione, quadratica Multiply, Chol decomposizione e more)
4. Si utilizza C compilatore e multi-thread calcolo per accelerare.

Per questo problema, è sufficiente scrivere questo comando:

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

qui è un punto di riferimento per tutti i metodi possibili. Per maggiori dettagli fare riferimento a questo domanda .

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

Answer 8

Vorrei pensare ricorsione, ma questo è l'unico altro metodo non-ciclo si può fare

Answer 9

Si potrebbe "unroll" del ciclo, cioè scrivere tutte le moltiplicazioni in sequenza che si verificherebbe nel loop

Answer 10

Vorrei condividere la mia risposta ai problemi di:

1) rendendo il prodotto tensoriale di due tensori (di qualsiasi valenza);

2) rendendo la contrazione di due tensori lungo qualsiasi dimensione.

Ecco le mie subroutine per il primo e secondo compiti:

1) prodotto tensoriale:

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

2) contrazione: Qui A e B sono i tensori da contrarre lungo rispettivamente le dimesions i e j. Le lunghezze di queste dimensioni devono essere uguali, naturalmente. Non c'è alcun controllo per questo (questo sarebbe oscurare il codice), ma a parte questo funziona bene.

   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's 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

Tutti i critici?