contraction Tensor dans Matlab

Question

Possible en double:
MATLAB: How to vector-multiply two arrays of matrices?

Est-il possible de contracter tenseurs plus dimensions dans Matlab?

Par exemple, supposons que j'ai deux tableaux 3 dimensions, avec ces dimensions:

size(A) == [M,N,P] 
size(B) == [N,Q,P] 

Je veux contracter A et B sur les deuxième et premiers indices, respectivement. En d'autres termes, je veux considérer A comme un tableau de matrices de taille [M,N] et B être un tableau de longueur égale de matrices [N,Q]; Je veux multiplier ces tableaux élément-par-élément (matrix-by-matrix) pour obtenir quelque chose de la taille [M,Q,P].

Je peux le faire via une boucle for:

assert(size(A,2) == size(B,1)); 
assert(size(A,3) == size(B,3)); 

M = size(A,1); 
P = size(A,3); 
Q = size(B,2); 

C = zeros(M, Q, P); 
for ii = 1:size(A,3) 
    C(:,:,ii) = A(:,:,ii) * B(:,:,ii); 
end 

Est-il possible de le faire qui évite la boucle for? (Et peut-être travaille avec des tableaux d'un nombre arbitraire de dimensions?)

Answer 1

Voici une solution (similaire à ce qui a été fait here) qui calcule le résultat dans une seule opération de multiplication matricielle, bien qu'il implique une manipulation lourde des matrices pour les mettre dans la forme désirée. Je compare alors à la simple boucle de calcul (qui je l'avoue est beaucoup plus lisible)

%# 3D matrices 
A = rand(4,2,3); 
B = rand(2,5,3); 
[m n p] = size(A); 
[n q p] = size(B); 

%# single matrix-multiplication operation (computes more products than needed) 
AA = reshape(permute(A,[2 1 3]), [n m*p])';  %'# cat(1,A(:,:,1),...,A(:,:,p)) 
BB = reshape(B, [n q*p]);       %# cat(2,B(:,:,1),...,B(:,:,p)) 
CC = AA * BB; 
[mp qp] = size(CC); 

%# only keep "blocks" on the diagonal 
yy = repmat(1:qp, [m 1]); 
xx = bsxfun(@plus, repmat(1:m,[1 q])', 0:m:mp-1); %' 
idx = sub2ind(size(CC), xx(:), yy(:)); 
CC = reshape(CC(idx), [m q p]); 

%# compare against FOR-LOOP solution 
C = zeros(m,q,p); 
for i=1:p 
    C(:,:,i) = A(:,:,i) * B(:,:,i); 
end 
isequal(C,CC) 

Notez que ce qui précède effectue plus multiplications que nécessaire, mais parfois "Anyone who adds, detracts (from execution time)". Malheureusement, ce n'est pas le cas, comme boucle FOR est beaucoup plus rapide ici :)

Mon but était de montrer que la vectorisation n'est pas facile, et que les solutions à base de boucles ne sont pas toujours mauvais ...