Matlab tensor reduction

Question

Possible duplicate:
MATLAB: How is the vector multiplication of two arrays of matrices?

Is there a way to compress multidimensional tensors in Matlab?

For example, suppose I have two 3-dimensional arrays with the following sizes:

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

I want to enclose A and B in the second and first indices respectively. In other words, I want to consider A as an array of matrices of size [M,N] and B as an equal array of matrix lengths [N,Q] ; I want to multiply these arrays by elements (matrix by matrix) to get something of size [M,Q,P] .

I can do this through a for loop:

 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 

Is there a way to do this to avoid the for loop? (And maybe it works with arrays of an arbitrary number of dimensions?)

Answer 1

Here is a solution (similar to what was done here ) that calculates the result in a single matrix multiplication operation, although it involves heavy manipulation of the matrices with place them in the desired shape. Then I compare it with a simple calculation for a loop (which, I admit, is more readable)

 %# 3D matrices A = rand(4,2,3); B = rand(2,5,3); [mnp] = size(A); [nqp] = size(B); %# single matrix-multiplication operation (computes more products than needed) AA = reshape(permute(A,[2 1 3]), [nm*p])'; %'# cat(1,A(:,:,1),...,A(:,:,p)) BB = reshape(B, [nq*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), [mqp]); %# compare against FOR-LOOP solution C = zeros(m,q,p); for i=1:p C(:,:,i) = A(:,:,i) * B(:,:,i); end isequal(C,CC) 

Note that the above performs more multiplications than necessary, but sometimes "Anyone who adds distracts (from runtime) . " Unfortunately, this is not the case, since the FOR loop is much faster here :)

My point was to show that vectorization is not easy, and that loop-based solutions are not always bad ...