How can I calculate dyadics in matlab without usin

Question

I was wondering if someone could help me with my problem.

Let say that I have the coordinates of MxN vectors in a tensor r of dimensions [M,N,3]. I would like to save in a 3M-by-3N block matrix all dyadic products r_0'*r_0, where r_0 is the vector r_0 = r(m,n,:) for some m and n, and I would like to do this without using for loops.

If haven't explain myself correctly, here is an example code that shows what I would like to obtain (but using for loops, of course):

N=10;
M=5;
r=rand(M,N,3);
Dyadic=zeros(3*M,3*N);
for m=1:M
    a1=3*m-2;
    a2=3*m;
    for n=1:N
        b1=3*n-2;
        b2=3*n;
        aux(3)=r(m,n,3);
        aux(2)=r(m,n,2);
        aux(1)=r(m,n,1);
        Dyadic(a1:a2,b1:b2)=transpose(aux)*aux
    end
end

Thanks in advance!

Answer 1

You need to use bsxfun(@times and then re-arrange elements to have the desired output -

%// Get the multipliication result
mat_mult = bsxfun(@times,permute(r,[1 2 4 3]),r);
 %// OR if you would like to keep mat_mult as 3D that could be potentially faster -
 %//     mat_mult = bsxfun(@times,reshape(r,[],3),permute(reshape(r,[],3),[1 3 2]));

%// Re-arrange elements to have them the way you are indexing in nested loops
Dyadic = reshape(permute(reshape(mat_mult,M,N,3,[]),[3 1 4 2]),M*3,N*3);

The major play about this solution is really the re-arrangement of elements after we have the multiplication result.

Quick runtime tests with the input r as 1000 x 1000 x 3 sized array, show that this bsxfun based approach gives over 20x speedup over the nested loop code listed in the question!