summing over a multi-dimensional cartesian product using nested arrayfun commands

I have an anonymous function with multiple arguments and I want to sum over all of them. For example, consider the simple function

U  = @(i,j) i+j;
I  = 1:3;  J = 1:3;

The command

sumJ = @(I,J) sum(arrayfun(@(j)sum(U(I,j)),J)

will sum over all pairs in the cartesian product IxJ

I'd like to do the analogous thing for

V  = @(i,j,k) i + j + k

I realize that there are many other ways of doing this simple summation, it's just an easy example. For my more complex problem, I believe I would need to use nested arrayfun commands, but can't figure out how to do this. Thanks very much for any help.

Answer 1

[Guillaume]

If I understood correctly, it may be simpler to ndgrid your arrays and have a single arrayfun

V = @(i,j,k,l) i*j/k+l;                %demo function
I = 1:10; J = 1:20; K = 2:5; L = 3:7;  %demo arrays

[II, JJ, KK, LL] = ndgrid(I, J, K, L);
sumall = sum(arrayfun(V, II(:), JJ(:), KK(:), LL(:)))

The downside of this method is that if your arrays are big enough, you're going to need a lot of memory for the ndgrid output.

A completely generic version of the above which works for n input arrays stored in a cell array:

fun = @(i,j,k,l) i*j/k+l; %some function with |n| input arguments.
arrays = {1:10, 1:20, 2:5, 3:7};  %|n| input arrays stored in cell array

griddedarrays = cell(size(arrays)); %variable to receive output of ndgrid
[griddedarrays{:}] = ndgrid(arrays{:});
griddedarrays = cellfun(@(m) reshape(m, [], 1), griddedarrays, 'UniformOutput', false);
sumall = sum(arrayfun(fun, griddedarrays{:}));

PS: It's a shame bsxfun is limited to two inputs as it would avoid the ndgrid entirely.

Answer 2

[Stephen23]

Maybe this can help:

U = @(m,n)m+n;
V = @(M,n)sum(arrayfun(@(m)U(m,n),M));
W = @(M,N)sum(arrayfun(@(n)V(M,n),N));

A = 1:3;
B = 1:3;
W(A,B)

But note that loops are often faster than arrayfun or cellfun. And even faster would be to consider using bsxfun instead of arrayfun or a loop:

X = bsxfun(@plus,A,B(:));
sum(X(:))

[Leo Simon]

Thanks, Steve, for the suggestion. The trouble is that this doesn't seem to generalized to 3 dimensions. For example, the obvious extension of what you did is:

U = @(m,n,o)m+n+o;
V = @(M,n,o)sum(arrayfun(@(m)U(m,n,o),M));
W = @(M,N,o)sum(arrayfun(@(n)V(M,N,o),N));
X = @(M,N,O)sum(arrayfun(@(n)V(M,N,P),P));
A = 1:3;
B = 1:3;
C = 1:3;
X(A,B,C)

Throws an error about nonscalar in UniformOutput, but I'm unable to fix it using UniformOutput. Can you show me how to make the generalization work?

[Stephen23]

You could combine bsxfun and cellfun:

>> U = @(m,n,o)m+n+o;
>> V = @(M,N,o)sum(sum(bsxfun(@(m,n)U(m,n,o),M,N(:))));
>> W = @(M,N,O)sum(cellfun(@(o)V(M,N,o),num2cell(O)));
>> W(1:3,1:3,1:4)
ans =  234

[Leo Simon]

Thanks Steven, could you please explain the role of N(:) in the definition of V, I tried replacing it with N, i.e.,

V = @(M,N,o)sum(sum(bsxfun(@(m,n)U(m,n,o),M,N)));

keeping the definitions of U and W the same, and then

W(1:3,1:3,1:4)

returned 78, not 234. So obviously, both definitions are legal, but do different things. Thanks for your patience! Leo

[Guillaume]

The purpose of N(:) is simply to transpose N. It is the same as N.'.

In Stephen's answer bsxfun is used to apply your function to the cartesian product of the vectors. Hence one vector needs to be along one dimension, the other vector along another. bsxfun is equivalent to repmat'ing each vector along the dimension of the other (without needing the memory for the repmat'ed matrix). My answer's basically does this repmat'ing explicitly, using ndgrid.

That's why I said it was a shame that bsxfun is limited to two inputs. An hypothetical bsxfun that allowed a arbitrary number of inputs would have allowed you to write:

W = bsxfun(V, I, J', permute(K, [3 2 1]), permute(L, [4 3 2 1]);
sumall = sum(W:));