How to create N+1 dimensional array by taking exterior product of 1st dimension of two N dimensional arrays?

I would like to efficiently create an N+1 dimensional array C, whose first 2 dimensions are the exterior product of the 1st dimension of the N dimensional arrays A and B. The exterior product of m by 1 column vectors x and y is the m by m matrix x*y'.

As an example, given the 2 by 3 by 5 arrays A and B, I would like to create the 2 by 2 by 3 by 5 array C such that C(i,j,k,l) = A(i,k,l)*B(j,k,l).

For efficiency pruposes, I would like to do this without for loops, to the extent possible (I know how to do this with slow, ugly, brute force for loops). Given that I will apply this to YALMIP sdpvar arrays, implicit expansion can't be used. The following (non-exhaustive list) can be used in any combination:

reshape
repmat
vec'ing (i.e., A(:))
.*
kron
bsxfun (if needed)

Thanks.

Answer 1

[Bruno Luong]

Solution without expansion as requested

sizeA = size(A);
reshapeB = reshape(B,[1 size(B)]);
reshapeA = reshape(A,[sizeA(1) 1 sizeA(2:end)]);
rA = ones(size(size(A))); rA(2)=size(B,1);
rB = ones(size(size(B))); rB(1)=size(A,1);
C = repmat(reshapeA,rA) .* repmat(reshapeB,rB) ;

[Mark Stone]

Thanks. And thanks also to @James Tursa who made a good first stab.

[James Tursa]

@Mark: I fixed the typos in my first post. I think you missed the point I was trying to make about the 1's, so I wrote it out explicitly for you.

[Bruno Luong]

I believe James's answer deserves the acceptation. My contribution is next to nothing.

[Mark Stone]

@Bruno Luong Perhaps you built on the answer by @James Tursa, but your answer fulfills my criteria, and his, evren after editing does not. Thanks to you both,.

Answer 2

[James Tursa]

Maybe try this:

sizeB = size(B);
reshapeA = reshape(A,[1 size(A)]);
reshapeB = reshape(B,[sizeB(1) 1 sizeB(2:end)]);
C = reshapeA .* reshapeB; % would do this if implicit expansion works for your case
C = bsxfun(@times,reshapeA,reshapeB); % try this if implicit expansion doesn't work for your case

Or you can swap where the 1's go in A and B if you want the 2D transpose in those first two dimensions.

EDIT

Corrected last lines above to use reshapeA and reshapeB.

Also my comment about the 1's was intended to steer you towards this solution if it was more appropriate:

sizeA = size(A);
reshapeA = reshape(A,[sizeA(1) 1 sizeA(2:end)]);
reshapeB = reshape(B,[1 size(B)]);
etc.

[Mark Stone]

Thanks for answering, but perhaps another iteration is needed to get what I need.

I presume the first solution was intended to form C with reshapeA and reshapeB, rather than A and B. I.e.,

sizeB = size(B);
reshapeA = reshape(A,[1 size(A)]);
reshapeB = reshape(B,[sizeB(1) 1 sizeB(2:end)]);
C = reshapeA .* reshapeB; 

I think this is correct for A and B being double, except for the matrix in the first two dimensions being the transpose of what it should be. In any event, it does not work for my inteneded use due to lack of proper support for implicit expansion. Could judicious use of repmat avoid the need for implicit expansion?

Instead of

C = bsxfun(@times,A,B);

do you mean

C = bsxfun(@times,reshapeA,reshapeB);

I think this latter formula produces the same result as

reshapeA .* reshapeB

and has the same transposition of the matrix in the 1st two diemnsions as

reshapeA .* reshapeB

It seems to require the same explict expansion as well, and therefore does not work for my inteneded use.

[Bruno Luong]

sizeA = size(A);
reshapeB = reshape(B,[1 size(B)]);
reshapeA = reshape(A,[sizeA(1) 1 sizeA(2:end)]);
C = bsxfun(@times,reshapeA,reshapeB);

or

sizeA = size(A);
reshapeB = reshape(B,[1 size(B)]);
reshapeA = reshape(A,[sizeA(1) 1 sizeA(2:end)]);
rA = ones(size(size(A))); rA(2)=size(B,1);
rB = ones(size(size(B))); rB(1)=size(A,1);
C = repmat(reshapeA,rA) .* repmat(reshapeB,rB) ;

[Mark Stone]

@Bruno Luong Thanks. I believe that corrects the matrix transposition issue. So now I am down to the matter of not havimng implicit expansion available.

[Bruno Luong]

My second solution does not contain any expansion, implicit or explicit.

[Mark Stone]

@Bruno Luong. Somehow I guess you edited in that 2nd solution after my reply. Yes, that 2nd solution, which does not require explicit expansion, seems to be exactly what I want. Thanks.

If you make that into a separate answer I can accept it.

[Mark Stone]

Yes, the edit by @James Tursa does fix the transposition issue, but I believe does still rely on implicit expansion, and therefore does not handle my situation. The answer by @Bruno Luong does not use implicit expansion and meets my needs. Thansk to you botj.

[James Tursa]

Sorry I apparently misunderstood. But according to your post, using bsxfun was acceptable so that is what I gave you with this line, which doesn't use implicit expansion:

C = bsxfun(@times,reshapeA,reshapeB);

However, as long as you have a solution that works for you from Bruno, all is good.