OK, so the subs matrix represent the rows and columns of the matrix A. (1, 1) represent the first row and first column. (3, 2) represents the 3rd row and 2nd column. If you don't get this, look at the document on ind2sub on this page. Now the vals are associated with each row in subs (e.g. the first (1, 1) associated with 101, up to (4, 1) with 106). And the values get accumulated in the groups determined by the rows of subs. Therefore, the two rows of (1, 1) have the values 101 and 104. This adds to 205, which is why A(1,1) = 205. (2, 2) are 102 and 105. Thus, A(2,2) = 207. (3,2) --> 103, A(3,2) = 103. (4,1) --> 106, A(4,1) = 106. It's pretty clear if you look at the last image on the accumarray page. https://www.mathworks.com/help/matlab/ref/accumarray.png

You could always use |accumarray| with two dimensions of _subs_ rather than _vals_: A=[1;1;2;3;3;3]; B=rand(6,10); [xx, yy] = ndgrid(A,1:size(B,2)); C=accumarray([xx(:) yy(:)],B(:)); Frankly, I'll bet the |for|-loop will be fastest. *More* And my suspicions are right depending on how big A and B actually are. I am seeing the |for|-loop be faster with small arrays and accumarray with two-dimensional subs being faster for larger ones. |arrayfun/cell2mat| etc are really slow since they require the conversion to cells which are slow and arrayfun is just slower than a |for|-loop anyway. function timeAccumarray A=[1;1;2;3;3;3]; B=rand(6,10); [t1,t2,t3] = deal(0); for ii = 1:100 tic [xx, yy] = ndgrid(A,1:size(B,2)); C=accumarray([xx(:) yy(:)],B(:)); t1 = t1+toc; tic for jj=size(B,2):-1:1 %dynamically preallocate, avoiding slow unique() C2(:,jj)=accumarray(A,B(:,jj)); end t2=t2+toc; tic C3=cell2mat(arrayfun(@(x) accumarray(A,B(:,x)),1:size(B,2),'un',0)); t3 = t3+toc; end isequal(C,C2,C3) [t1 t2 t3] _ans = 0.0072 0.0289 0.1381_

I think the unexpected values are because you are expecting the values to be passed to the anonymous function in the same order as the relevant indices appears in subs.

However, the documentation for accumarray says:

Note If the subscripts in subs are not sorted, fun should not depend on the order of the values in its input data.

This means that the start of subs ( 1 2; 1 2; ) can result in 102 101 or 101 102 being passed to the anonymous function (because the documentation says that fun should not depend on the order of values.

sum(diff([102 101])) is -1 which is the value written to element 1,2 in the example

val = 101:106;
subs=[1 2; 1 2; 3 1; 4 1; 4 4; 4 1];
B = accumarray(subs,val,[],@(x)sum(diff(x)))

B =

     0    -1     0     0
     0     0     0     0
     0     0     0     0
     2     0     0     0

I'm not sure why Mathworks would choose to use the diff function in an example when the result of diff is obviously dependent on the order of values passed to it.

It looks like this input value ordering is the reason for the other output values being different from your expectations.

It is something of a challenge to figure out when your first encounter it, but when you think about it and work with it, it’s not all that difficult to understand. In this example, rows *|1|* and *|4|* in *|‘subs’|* are both *|(1,1)|*. They correspond to rows *|1|* and *|4|* in *|‘val’|*, those being *|101|* and *|104|*. The function takes those values, adds them, (or ‘accumulates’ them), and puts that value, *|205|* in the *|(1,1)|* position in the output matrix. The *|‘subs’|* argument does not define the *|(1,2)|* position anywhere, so that element remains at *|0|*. The others work the same way. This isn’t a full explanation, but it should get you started in understanding an extremely useful function.

What exactly are you trying to do? Typically, the first input arg is a column vector and you get a group count, as you already know. If you pass a matrix in the first input arg, you are asking it to shape the output in a different ways. In this case, by passing [1,2,3] , you are asking it to return the result in 1x2x3 3D matrix. Why on earth anyone wants it I have no idea, but you can do it if you want to.

In your accumarray calls you're specifying the third input (sz) as an empty array, which tells accumarray basically "you figure out how large the result should be." If you know how large you want it to be, provide that information to accumarray. locations = [1 2; 2 1]; values = [2; 3]; A = accumarray(locations, values, []) B = accumarray(locations, values, [3 3]) In the example above, A is 2-by-2. It doesn't need any more rows and columns to contain the values at the coordinates given in the locations matrix. B is 3-by-3, as while the values do technically fit in a 2-by-2 matrix you explicitly asked for the result to be 3-by-3.

Yes, |accumarray| will do it easily. You just need to generate the correct subscripts for it, which is the cartesian product of all the indices in all the dimension, replacing the appropriate dimensions with your grouping variables: A = randi(100, [3 4 2 7]); %demo data %generate indices in all the dimension, and replace appropriate one with grouping variables: indices = arrayfun(@(s) 1:s, size(A), 'UniformOutput', false); %vector of indices in each dimension indices{2} = [2 1 2 3]; %grouping in 2nd dimension indices{4} = [3 3 2 1 3 2 4]; %grouping in 4th dimension %calculate cartesian product of indices [indices{:}] = ndgrid(indices{:}); indices = cell2mat(cellfun(@(v) v(:), indices, 'UniformOutput', false)); %apply accumarray with default sum function B = accumarray(indices, A(:)); _Edit_: note that the above *will work* with arrays of any dimension (within reason, ndgrid will run out of memory for too many dimensions).

For vals you can set it as 1:size(vals2,1) and use it to extract rows of vals2. Also it is required for the function to return cell.

result2 = accumarray(subs, 1:size(vals2,1), [], @(x) {prod(1+vals2(x,:),1)-1})

You can concatenate cell elements:

result3 = vertcat(result2{:})

Or all in one line:

result3 = cell2mat( accumarray(subs, 1:size(vals2,1), [], @(x) {prod(1+vals2(x,:),1)-1}))

result3 =

   0.38600   0.38600
   0.76635   0.76635
   0.18720   0.18720
   0.27008   0.27008

Result of a test in Octave comparing three proposed methods using a [10000 x 200] matrix as input:

subs = randi(1000,10000,1);
vals2 = rand(10000,200);

=========CELL2MAT========
Elapsed time is 0.130961 seconds.
=========NDGRID========
Elapsed time is 3.96383 seconds.
=========FOR LOOP========
Elapsed time is 6.16265 seconds.

Online Demo

One solution is to replicate the row indices in labels and add another column of column indices. Then you can reshape X into a column vector and apply accumarray once:

labels = [repmat(labels(:),nCols,1) ...            % Replicate the row indices
          kron(1:nCols,ones(1,numel(labels))).'];  % Create column indices
totals = accumarray(labels,X(:));  % I used "totals" instead of "means"

How it works...

A = accumarray(subs,val) for a column vector subs and vector val works by adding the number in val(i) to the total in row subs(i) in the output column vector A. However, subs can contain more than just row indices. It can contain subscript indices for multiple dimensions to assign values to in the output. This feature is what allows you to handle an input val that is a matrix instead of a vector.

First, the input for val can be reshaped into a column vector using the colon operator X(:). Next, in order to keep track of which column in the output the values in X(:) should be placed, we can modify the input subs to include an additional column index. To illustrate how this works, I'll use these sample inputs:

labels = [3; 1; 1];
X = [1 2 3; ...
     4 5 6; ...
     7 8 9];
nCols = 3

And here are what the variables in the above code end up looking like:

labels = 3 1    X(:) = 1    totals = 11 13 15
         1 1           4              0  0  0
         1 1           7              1  2  3
         3 2           2
         1 2           5
         1 2           8
         3 3           3
         1 3           6
         1 3           9

Notice, for example, that the values 1 4 7 that were originally in the first column of X will only be accumulated in the first column of the output, as denoted by the ones in the first three rows of the second column of labels. The resulting output should be the same as what you would have gotten by using the code in the question where you loop over each column to perform the accumulation.

This wasn't an option in the release the original poster was using (the two functions required were released in release R2016b, which came out a few months after this question was asked) but today I suggest turning your table into a timetable using table2timetable then call retime on the resulting timetable.