Multiplier une matrice 3D avec une matrice 2D

Question

Supposons que j'ai une matrice AxBxC X et une matrice BxD Y.

Y at-il une méthode non-boucle par laquelle je peux multiplier chacun des C AxB matrices avec Y?

Answer 1

Vous pouvez le faire en une seule ligne en utilisant les fonctions NUM2CELL pour briser la matrice X dans une matrice de cellules et CELLFUN faire fonctionner à travers les cellules:

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

Le résultat est un Z 1 par C réseau de cellules où chaque cellule contient un A-par-D matrice. Si vous voulez Z être A par D par C matrice, vous pouvez utiliser le la fonction de CAT:

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

---

Remarque: Mon ancienne solution utilisée MAT2CELL au lieu de NUM2CELL , qui n'a pas été aussi succincte:

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

Answer 2

Comme une préférence personnelle, j'aime mon code soit aussi succinct et lisible que possible.

Voici ce que je l'aurais fait, même si elle ne répond pas à vos besoins « sans boucles »:

for m = 1:C

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

end

Il en résulte une A x D x C de la matrice Z .

Et bien sûr, vous pouvez toujours pré-allouer Z pour accélérer les choses en utilisant Z = zeros(A,D,C);.

Answer 3

Voici une solution d'une ligne (deux si vous voulez diviser en 3ème dimension):

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

Par conséquent maintenant: Z(:,:,i) contient le résultat de X(:,:,i) * Y

---

Explication:

Le dessus peut paraître confus, mais l'idée est simple. Prenons d'abord je commence par la troisième dimension de X et faire une concaténation verticale le long de la première dim:

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

... La difficulté était que C est une variable, donc vous ne pouvez pas généraliser cette expression en utilisant chat ou vertcat . Ensuite, nous multiplier par Y:

ZZ = XX * Y;

Enfin je partage le dans la troisième dimension:

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

Vous pouvez le voir ne nécessite qu'une seule multiplication de matrices, mais vous devez Reshape la matrice avant et après.

Answer 4

J'aborde exactement le même problème, avec un oeil pour la méthode la plus efficace. Il y a à peu près trois approches que je vois autour, à court d'utiliser des bibliothèques externes (c.-à- mtimesx ):

1. boucle par tranches de la matrice 3D
2. repmat-et permute wizardry
3. multiplication cellfun

J'ai récemment comparé les trois méthodes pour voir ce qui était le plus rapide. Mon intuition est que (2) serait le vainqueur. Voici le code:

% 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

Les trois approches ont produit la même sortie (ouf!), Mais, de façon surprenante, la boucle a été le plus rapide:

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

Notez que les temps peuvent varier beaucoup d'un essai à l'autre, et parfois (2) sort le plus lent. Ces différences sont plus spectaculaires avec des données plus importantes. Mais avec beaucoup plus grandes données, (3) temps (2). La méthode de la boucle est encore mieux.

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

Mais la méthode de la boucle peut être plus lent que (2), si la dimension est beaucoup plus grande boucle que les autres.

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) gagne par un facteur important, dans ce (peut-être extrême) cas. Il ne peut pas être une approche qui est optimale dans tous les cas, mais la boucle est encore assez bon, et le meilleur dans de nombreux cas. Il est également préférable en termes de lisibilité. Boucle magnétique!

Answer 5

Non. Il y a plusieurs façons, mais il revient toujours dans une boucle, directe ou indirecte.

Juste pour faire plaisir à ma curiosité, pourquoi voudriez-vous que de toute façon?

Answer 6

Pour répondre à la question, et pour une meilleure lisibilité, s'il vous plaît voir:

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

Entrée

• 2 tableaux
• dim

Exemple

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

Source

source originale. J'ai ajouté des commentaires inline.

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

Je recommande fortement d'utiliser le boîte à outils MMX de Matlab. On peut multiplier les matrices n dimensions aussi vite que possible.

Les avantages de MMX sont:

1. facile à utiliser.
2. Multiplier matrices n dimensions (en fait, on peut multiplier les ensembles de matrices 2-D)
3. Il effectue d'autres opérations matricielles (transposition, Quadratic Multiply, décomposition Chol et plus)
4. Il utilise compilateur C et multi-thread calcul jusqu'à la vitesse.

Pour ce problème, il vous suffit d'écrire cette commande:

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

ici est une référence pour toutes les méthodes possibles. Pour plus de détails se référer à cette 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  <===================

Answer 8

Je pense récursion, mais c'est la seule autre méthode non-boucle que vous pouvez faire

Answer 9

Vous pouvez « déroulez » la boucle, à savoir écrire toutes les multiplications séquentielle qui se produirait dans la boucle

Answer 10

Je voudrais partager ma réponse aux problèmes de:

1) la fabrication du produit tensoriel de deux tenseurs (quelle valence);

2) rendant la contraction de deux tenseurs le long de toute dimension.

Voici mes sous-routines pour les première et deuxième tâches:

1) produit tensoriel:

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

2) la contraction: Ici, A et B sont des tenseurs d'être traitées le long des dimesions i et j respectivement. Les longueurs de ces dimensions doivent être égales, bien sûr. Il n'y a pas vérifier pour cela (ce masquerait le code), mais à part cela, il fonctionne bien.

   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

Les critiques?