Creating and manipulating three dimensional matrices in Matlab

Question

I'm desperately trying to avoid a for loop in Matlab, but I cannot figure out how to do it. Here's the situation:

I have two m x n matrices A and B and two vectors v and w of length d. I want to outer multiply A and v so that I get an m x n x d matrix where the (i,j,k) entry is A_(i,j) * v_k, and similarly for B and w.

Afterward, I want to add the resulting m x n x d matrices, and then take the mean along the last dimension to get back an m x n matrix.

I'm pretty sure I could handle the latter part, but the first part has me completely stuck. I tried using bsxfun to no avail. Anyone know an efficient way to do this? Thanks very much!

EDIT: This revision comes after the three great answers below. gnovice has the best answer to the question I asked without a doubt. However,the question that I meant to ask involves squaring each entry before taking the mean. I forgot to mention this part originally. Given this annoyance, both of the other answers work well, but the clever trick of doing algebra before coding doesn't help this time. Thanks for the help, everyone!

Answer 1

EDIT:

Even though the problem in the question has been updated, an algebraic approach can still be used to simplify matters. You still don't have to bother with 3-D matrices. Your result is just going to be this:

output = mean(v.^2).*A.^2 + 2.*mean(v.*w).*A.*B + mean(w.^2).*B.^2;

If your matrices and vectors are large, this solution will give you much better performance due to the reduced amount of memory required as compared to solutions using BSXFUN or REPMAT.

Explanation:

Assuming M is the m-by-n-by-d matrix that you get as a result before taking the mean along the third dimension, this is what a span along the third dimension will contain:

M(i,j,:) = A(i,j).*v + B(i,j).*w;

In other words, the vector v scaled by A(i,j) plus the vector w scaled by B(i,j). And this is what you get when you apply an element-wise squaring:

M(i,j,:).^2 = (A(i,j).*v + B(i,j).*w).^2;
            = (A(i,j).*v).^2 + ...
              2.*A(i,j).*B(i,j).*v.*w + ...
              (B(i,j).*w).^2;

Now, when you take the mean across the third dimension, the result for each element output(i,j) will be the following:

output(i,j) = mean(M(i,j,:).^2);
            = mean((A(i,j).*v).^2 + ...
                   2.*A(i,j).*B(i,j).*v.*w + ...
                   (B(i,j).*w).^2);
            = sum((A(i,j).*v).^2 + ...
                  2.*A(i,j).*B(i,j).*v.*w + ...
                  (B(i,j).*w).^2)/d;
            = sum((A(i,j).*v).^2)/d + ...
              sum(2.*A(i,j).*B(i,j).*v.*w)/d + ...
              sum((B(i,j).*w).^2)/d;
            = A(i,j).^2.*mean(v.^2) + ...
              2.*A(i,j).*B(i,j).*mean(v.*w) + ...
              B(i,j).^2.*mean(w.^2);

Answer 2

Try reshaping the vectors v and w to be 1 x 1 x d:

  mean (bsxfun(@times, A, reshape(v, 1, 1, [])) ...
        + bsxfun(@times, B, reshape(w, 1, 1, [])), 3)

Here I am using [] in the argument to reshape to tell it to fill that dimension in based on the product of all the other dimensions and the total number of elements in the vector.

Answer 3

Use repmat to tile the matrix in the third dimension.

A =

     1     2     3
     4     5     6

>> repmat(A, [1 1  10])

ans(:,:,1) =

     1     2     3
     4     5     6


ans(:,:,2) =

     1     2     3
     4     5     6

etc.

Answer 4

You still don't have to resort to any explicit loops or indirect looping using bsxfun et al. for your updated requirements. You can achieve what you want by a simple vectorized solution as follows

output = reshape(mean((v(:)*A(:)'+w(:)*B(:)').^2),size(A));

Since OP only says that v and w are vectors of length d, the above solution should work for both row and column vectors. If they are known to be column vectors, v(:) can be replaced by v and likewise for w.

---

You can check if this matches Lambdageek's answer (modified to square the terms) as follows

outputLG = mean ((bsxfun(@times, A, reshape(v, 1, 1, [])) ...
        + bsxfun(@times, B, reshape(w, 1, 1, []))).^2, 3);

isequal(output,outputLG)

ans =

     1