1

I'm trying to calculate the outer product of all rows in a matrix. So far I'm doing the following:

A = rand(10,8);
[N J] = size(A);
for i = 1:N,
    B(((i-1)*J+1):J+((i-1)*J),:) = A(i,:)'*A(i,:)
end

I then take the mean of the same elements of the same row in each submatrix created above.

for j=1:J,
    C(1:J,j) = accumarray(repmat((1:J)',N,1),B(:,j)).*(1/N); 
end

Now this works, but my input matrix will eventually be a number of magnitudes larger, so a vectorized version would be nice. I assume there is some way to do this with permute and bsxfun but the solutions for outer product I have seen so far don't seem to apply here.

TiredHornet
  • 49
  • 4

1 Answers1

1

Sure. You can compute C directly as

C = mean(bsxfun(@times, permute(A, [2 3 1]), permute(A, [3 2 1])), 3);

Or, if you really need the B variable, you can do it this way:

B_3D = bsxfun(@times, A.', permute(A,[3 1 2])); %'// 3D array
B = reshape(B_3D, numel(A), []);
C = mean(permute(B_3D, [1 3 2]), 3);
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • This answer is great and results look to be equal. But if I try to verify with iseqaul() then it returns 0. Why could this be the case? – TiredHornet Jul 24 '14 at 12:54
  • @TiredHornet Most probably floating point precision issues. How about manually inspect some of the values? – Divakar Jul 24 '14 at 13:37
  • You should always [__beware of comparing reals__](http://stackoverflow.com/questions/13699596/is-this-a-matlab-bug-do-you-have-the-same-issue), for the reason Divakar has said. Instead, try `yourC ./ myC` and check it gives all entries equal to `1` except for finite precission issues – Luis Mendo Jul 24 '14 at 15:39
  • @Luis Mendo: Very nice idea. I did that and it I get all ones. – TiredHornet Jul 25 '14 at 11:18
  • @TiredHornet Note they are not actually ones, but very close (try for example `yourC./myC-1`). That's normal; it happens because of finite numerical precision – Luis Mendo Jul 25 '14 at 11:53
  • @Luis Mendo Okay, I'll remember that. Now, it is time to learn how to use permute on my own. Hard to get your head this flipping of dimensions... – TiredHornet Jul 25 '14 at 12:24
  • @TiredHornet Yes, it's hard initially, but you kind-of get used :-) The basic idea is to make dimensions non-coincident if you want `bsxfun` to compute all "combinations" involving those dimensions – Luis Mendo Jul 25 '14 at 12:26