6

I have a relatively large matrix NxN (N~20,000) and a Nx1 vector identifying the indices that must be grouped together.

I want to sum together parts of the matrix, which in principle can have a different number of elements and non-adjacent elements. I quickly wrote a double for-loop that works correctly but of course it is inefficient. The profiler identified these loops as one of the bottlenecks in my code.

I tried to find a smart vectorization method to solve the problem. I explored the arrayfun, cellfun, and bsxfun functions, and looked for solutions to similar problems... but I haven't found a final solution yet.

This is the test code with the two for-loops:

M=rand(10); % test matrix
idxM=[1 2 2 3 4 4 4 1 4 2]; % each element indicates to which group each row/column of M belongs
nT=size(M,1);
sumM=zeros(max(idxM),max(idxM));
for t1=1:nT
    for t2=1:nT
        sumM(t1,t2) = sum(sum(M(idxM==t1,idxM==t2)));
    end
end
anothermh
  • 6,482
  • 3
  • 20
  • 39
user9998992
  • 125
  • 6
  • 2
    Perhaps worth noting that `arrayfun` and `cellfun` are basically loops in disguise, it's quite likely that optimised `for` loops would be just as quick (if not as concise). – Wolfie Jun 27 '18 at 10:45
  • Thank you for pointing this out. What about bsxfun that you used in your answer? – user9998992 Jun 27 '18 at 10:57
  • 2
    I assume you mean in Luis answer... `bsxfun` is a different beast and can often improve speed (depending how it is used). Think of `cellfun` and `arrayfun` as *"I want to loop through an array and apply a function to each element"*, whereas `bsxfun` is *"expand these vectors to equivalent size then perform a matrix operation"* so isn't element-by-element. They have different useful applications, `bsxfun` is useful for numerical operations between arrays and matrices, whereas the looping functions are good for concise code without having to initialise outputs etc. – Wolfie Jun 27 '18 at 11:03
  • 1
    @user9998992 `bsxfun` is [very fast](https://stackoverflow.com/questions/29719674/comparing-bsxfun-and-repmat), as is [implicit expansion](https://stackoverflow.com/q/42559922/2586922) – Luis Mendo Jun 27 '18 at 11:10
  • In your example you're calculating the values for 100 sums, but only 16 are actually relevant. If you have the same level of duplication in `idxM` in your actual data, I'd strongly suggest limiting the calculations to only the unique values of `idxM`. – beaker Jun 27 '18 at 17:37
  • 1
    Yes, agreed. The lines: `nT=size(M,1); sumM=zeros(max(idxM),max(idxM));` should be exchanged with the lines: `nT=max(idxM); sumM=zeros(nT,nT);` (in my original code, they were correct, but I put typos in my test code) – user9998992 Jun 27 '18 at 19:33
  • You could also use `idxU = unique(idxM); for t1=idxU, for t2=idxU...`, which would skip any unused indices in the range `1:max(idxM)`. – beaker Jun 27 '18 at 21:20

3 Answers3

3

You can use accumarray as follows:

nT = size(M,1); % or nT = max(idxM)
ind = bsxfun(@plus, idxM(:), (idxM(:).'-1)*nT); % create linear indices for grouping
sumM = accumarray(ind(:), M(:), [nT^2 1]); % compute sum of each group
sumM = reshape(sumM, [nT nT]); % reshape obtain the final result
Luis Mendo
  • 106,541
  • 12
  • 66
  • 138
  • 1
    And since R2016b you can rely on MATLAB's implicit expansion and don't need `bsxfun`... line 2 becomes `ind = idxM(:) + (idxM(:).'-1)*nT` – Wolfie Jun 27 '18 at 10:43
  • @Wolfie I find it hard to let go of `bsxfun` :-) – Luis Mendo Jun 27 '18 at 10:44
  • 1
    Before redrafting I did have some snark in my comment about how much you love it! In fairness it does add explicit clarity as to what the line does. – Wolfie Jun 27 '18 at 10:46
  • 2
    @Wolfie don't forget backward compatibility :) – Dev-iL Jun 27 '18 at 10:50
  • 2
    @Wolfie Yes, that's an important reason too. In fact, when I've used implicit expansion in my code I've felt obliged to include a comment, as it's easy to overlook what's going on. So `bsxfun` is clearer – Luis Mendo Jun 27 '18 at 10:53
  • 1
    @LuisMendo: This is the first time I see "`bsxfun` is clearer" -- this is the most obscure, obtuse function name they could come up with. It took me, literally, three years to memorize that name. I'm glad we don't need it any more. :) – Cris Luengo Jun 27 '18 at 12:52
  • 1
    @CrisLuengo I see your point. But when you see `bsxfun` you know that something a bit tricky is going on. It prompts you to look carefully. With implicit expansion the same tricky thing happens, only without any signal about it. If you see `x.*y`, is implicit expanstion taking place there? And if so, was it really the author's intent? – Luis Mendo Jun 27 '18 at 13:22
  • 1
    @LuisMendo: I know what you mean. I guess it's about what you're used to. I've had implicit singleton expansion for images in my toolbox for many, many years. So I'm used to it and it feels totally natural. It usually is not all that tricky what is going on. – Cris Luengo Jun 27 '18 at 13:31
  • The proposed solution is definitely elegant. However, for matrices M of dimensions smaller than about 10,000x10,000 on my laptop it is only as fast as the original ugly for-loops, but using ~50% more RAM. For larger matrices M, the computation time explodes on my laptop, probably due to RAM limitations. I can try with the office PC... – user9998992 Jun 27 '18 at 19:39
2

A solution using cumsum and diff.

[s,is] = sort(idxM);
sumM  = M(is,is);
idx = [diff(s)~=0 ,true];
CS = cumsum(sumM);
CS = cumsum(CS(idx,:),2);
n=sum(idx);
result = diff([zeros(n,1) diff([zeros(1,n); CS(:,idx)])],1,2);
sumM (:)=0;
sumM (s(idx),s(idx))=result;
rahnema1
  • 14,144
  • 2
  • 12
  • 26
  • The proposed solution is quite clever. However, as observed also for the other proposed solution, for matrices M of dimensions smaller than about 10,000x10,000 on my laptop it is only as fast as the original ugly for-loops, but using ~50% more RAM. For larger matrices M, the computation time explodes on my laptop, probably due to RAM limitations. – user9998992 Jun 27 '18 at 19:40
1

I'd like to point those who are interested to this answer provided on another forum

S=sparse(1:N,idxM,1); sumM=S.'*(M*S);

Credits (and useful discussion):

https://www.mathworks.com/matlabcentral/answers/407634-how-to-sum-parts-of-a-matrix-of-different-sizes-without-using-for-loops

user9998992
  • 125
  • 6