Efficiently compute a 3D matrix of outer products - MATLAB

Question

Suppose I have a matrix of elements like so:

A = reshape(1:25, 5, 5)

A =

 1     6    11    16    21
 2     7    12    17    22
 3     8    13    18    23
 4     9    14    19    24
 5    10    15    20    25

I would like to efficiently compute a 3D matrix of outer products, such that the i^th slice of this output matrix is the outer product of the i^th column of A with itself. The outer product between two vectors u and v is simply u*v.' if u and v are both column vectors.

Therefore, each slice of this output matrix B should be structured such that:

B(:,:,1) = A(:,1) * A(:,1).';
B(:,:,2) = A(:,2) * A(:,2).';
        ...
        ...
B(:,:,5) = A(:,5) * A(:,5).';

My current method is the following. I tried doing it this way using arrayfun and cell2mat:

cellmatr = arrayfun(@(x) A(:,x) * A(:,x).', 1:size(A,2), 'uni', 0);
out = reshape(cell2mat(cellmatr), size(A,1), size(A,1), size(A,2));

I simply loop over a linear index array between 1 and as many columns we have in A, and for each element in this array, I access the corresponding column and compute the outer product. The output will thus give a 1D grid of cells, which I then convert back into a 2D matrix, then reshape into a 3D matrix to find the 3D matrix of outer products.

However, for large matrices, this is quite slow. I've also tried replacing the matrix product with kron (i.e. kron(A(:,x), A(:,x))) inside my arrayfun call, but this is still quite slow for my purposes.

---

Does anyone know of an efficient way to compute this 3D matrix of outer products in this fashion?

Answer 1

This is just a minor improvement over Divakar's answer. It is a little faster because it replaces a 3D-array permute with a 2D-array permute:

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

Answer 2

To state the obvious, have you tried a simple for-loop:

[m,n] = size(A);
B = zeros(m,m,n);
for i=1:n
    B(:,:,i) = A(:,i) * A(:,i).';
end

You'll be surprised how competitively fast it is.

Answer 3

This -

B = permute(bsxfun(@times,A,permute(A,[3 2 1])),[1 3 2])