1

I have a 5D matrix A, and I need to multiply the 3rd-5th dimensions with a vector. For example, see the following sample code:

A=rand(50,50,10,8,6);
B=rand(10,1);
C=rand(8,1);
D=rand(6,1);

for i=1:size(A,3)
    for j=1:size(A,4)
        for K=1:size(A,5)
            A(:,:,i,j,K)=A(:,:,i,j,K)*B(i)*C(j)*D(K);
        end
    end
end

I wonder if there's a better \ vectorized \ faster way to do this?

Max
  • 85
  • 6

1 Answers1

4

Firstly, as a note, these days in Matlab, with JIT compilation, vectorised code is not necessarily faster/better. For big problems the memory usage in particular can cause performance issues.

Nevertheless, here is a vectorised solution that seems to give the same results as your code:

A=rand(3,4,5,6,7);
B=rand(5,1);
C=rand(6,1);
D=rand(7,1);   

s=size(A);
[b,c,d]=ndgrid(B,C,D);
F=b.*c.*d;
G=zeros(1,1,prod(s(3:5)));
G(1,1,:)=F(:);
A=reshape(A,s(1),s(2),[]);
A=bsxfun(@times,A,G);
A=reshape(A,s);

EDIT: An alternate solution:

A=bsxfun(@times,A,permute(B,[2 3 1]));
A=bsxfun(@times,A,permute(C,[2 3 4 1]));
A=bsxfun(@times,A,permute(D,[2 3 4 5 1]));
David
  • 8,449
  • 1
  • 22
  • 32
  • Did you check both results with `isequal` method? It does not say that they are the same. They are different even by dimensions. – Cheery Oct 02 '14 at 00:41
  • @Cheery try now, `isequal` wont work because we are doing floating point operations (the multiplication), but the maximum difference between my code and OP's code is roughly `1e-16`. – David Oct 02 '14 at 00:44
  • Checked using `randi` function.. yes, now they are the same ) – Cheery Oct 02 '14 at 00:47
  • thanks Dan, for some reason when I use my real matrix (size (50,50,10,8,6)) I get that `Error using bsxfun Non-singleton dimensions of the two input arrays must match each other.` Any reason why this happens? I've edited the question with the proper sizes... – Max Oct 02 '14 at 00:52
  • @Max clear `G` from your workspace before running the code and see if that helps. – David Oct 02 '14 at 01:02
  • Or add `G=zeros(1,1,prod(s(3:5)));` on the line before `G(1,1,:)=F(:);`. – David Oct 02 '14 at 01:10
  • 2
    +1, FWIW, you can replace the G lines with just `A=bsxfun(@times,A,permute(F(:),[3 2 1]))` ... – bla Oct 02 '14 at 06:56
  • thanks, @David you've been very helpful! I've just asked a second question that you may be interested to see. http://stackoverflow.com/questions/26167321/vectorizing-a-nested-loop-where-one-loop-variable-depends-on-the-other – Max Oct 02 '14 at 18:39
  • @David, I like the new solution, it's short and fast. Do you think it can be applied somehow to a slightly more elaborate loop, like the one mentioned here http://stackoverflow.com/questions/26167321/vectorizing-a-nested-loop-where-one-loop-variable-depends-on-the-other ? – Max Oct 03 '14 at 00:55
  • 2
    @Max I've spent several hours trying, but no luck so far. – David Oct 03 '14 at 03:29
  • thanks david, I deeply appreciate it. @natan gave an answer that works based on your solution. I just hoped you can so something like that using only bsxfun and permute... – Max Oct 03 '14 at 16:31