title: efficient implementation of a tensor dot product in matlab
body: <p>In MATLAB I have a series of 2x2 matrices stacked into a 3D tensor and I'd like to perform a matrix multiplication for each instance of matrices.</p>
  
  <p>so my C = A * B is defined as</p>
  
  <pre><code>C_ijk = sum(a_ilk * b_ljk, over all l)
  </code></pre>
  
  <p>my current implementation looks like this</p>
  
  <pre><code>    function mats = mul3D(A, B)
          % given a list of 2D matrices (e.g. rotation matrices) applies
          % the matrix product for each instance along the third dimension
          % mats(:,:,i) = A(:,:,i) * B(:,:,i) for all i
          % for this to succeed matrix dimensions must agree.
          mats = zeros(size(A,1), size(B,2), size(B,3));
          for i=1:size(B, 3)
              mats(:,:,i) = A(:,:,i) * B(:,:,i);
          end
      end
  </code></pre>
  
  <p>which is very easy to read but I remember someone saying that MATLAB doesn't like for-loops.</p>
  
  <p>So can you think of a better implementation that doesn't consume more memory than this one while being faster? my code spends about 50% of the run time in this for loop.</p>
  
  <p><strong>edit</strong></p>
  
  <p>thanks for your suggestions. Unfortunately I'm unable to introduce new dependencies to 3rd party code.</p>
  
  <p>Based on your questions I had the idea of exploiting the 2 x 2 x n structure of the tensors. My latest implementation looks like this:</p>
  
  <pre><code>    function mats = mul3D(A, B)
          % given a list of 2D matrices (e.g. rotation matrices) applies
          % the matrix product for each instance along the third dimension
          % mats(:,:,i) = A(:,:,i) * B(:,:,i) for all i
          % for this to succeed matrix dimensions must agree.
          mats = zeros(size(A,1), size(B,2), size(B,3));
  
          mats(1,1,:) = A(1,1,:) .* B(1,1,:) + A(1,2,:) .* B(2,1,:);
          mats(2,1,:) = A(2,1,:) .* B(1,1,:) + A(2,2,:) .* B(2,1,:);
          if(size(mats,2) &gt; 1)
              mats(1,2,:) = A(1,1,:) .* B(1,2,:) + A(1,2,:) .* B(2,2,:);
              mats(2,2,:) = A(2,1,:) .* B(1,2,:) + A(2,2,:) .* B(2,2,:);
          end
      end
  </code></pre>
  
  <p>any further suggestions are appreciated!</p>