Generate a matrix containing all combinations of elements taken from n vectors

The ndgrid function almost gives the answer, but has one caveat: n output variables must be explicitly defined to call it. Since n is arbitrary, the best way is to use a comma-separated list (generated from a cell array with ncells) to serve as output. The resulting n matrices are then concatenated into the desired n-column matrix:

vectors = { [1 2], [3 6 9], [10 20] }; %// input data: cell array of vectorsn = numel(vectors); %// number of vectorscombs = cell(1,n); %// pre-define to generate comma-separated list[combs{end:-1:1}] = ndgrid(vectors{end:-1:1}); %// the reverse order in these two%// comma-separated lists is needed to produce the rows of the result matrix in%// lexicographical order combs = cat(n+1, combs{:}); %// concat the n n-dim arrays along dimension n+1combs = reshape(combs,[],n); %// reshape to obtain desired matrix

A little bit simpler ... if you have the Neural Network toolbox you can simply use combvec:

vectors = {[1 2], [3 6 9], [10 20]};combs = combvec(vectors{:}).' % Use cells as arguments

which returns a matrix in a slightly different order:

combs =     1     3    10     2     3    10     1     6    10     2     6    10     1     9    10     2     9    10     1     3    20     2     3    20     1     6    20     2     6    20     1     9    20     2     9    20

If you want the matrix that is in the question, you can use sortrows:

combs = sortrows(combvec(vectors{:}).')% Or equivalently as per @LuisMendo in the comments: % combs = fliplr(combvec(vectors{end:-1:1}).') 

which gives

combs =     1     3    10     1     3    20     1     6    10     1     6    20     1     9    10     1     9    20     2     3    10     2     3    20     2     6    10     2     6    20     2     9    10     2     9    20

If you look at the internals of combvec (type edit combvec in the command window), you'll see that it uses different code than @LuisMendo's answer. I can't say which is more efficient overall.

If you happen to have a matrix whose rows are akin to the earlier cell array you can use:

vectors = [1 2;3 6;10 20];vectors = num2cell(vectors,2);combs = sortrows(combvec(vectors{:}).')

I've done some benchmarking on the two proposed solutions. The benchmarking code is based on the timeit function, and is included at the end of this post.

I consider two cases: three vectors of size n, and three vectors of sizes n/10, n and n*10 respectively (both cases give the same number of combinations). n is varied up to a maximum of 240 (I choose this value to avoid the use of virtual memory in my laptop computer).

The results are given in the following figure. The ndgrid-based solution is seen to consistently take less time than combvec. It's also interesting to note that the time taken by combvec varies a little less regularly in the different-size case.

Benchmarking code

Function for ndgrid-based solution:

function combs = f1(vectors)n = numel(vectors); %// number of vectorscombs = cell(1,n); %// pre-define to generate comma-separated list[combs{end:-1:1}] = ndgrid(vectors{end:-1:1}); %// the reverse order in these two%// comma-separated lists is needed to produce the rows of the result matrix in%// lexicographical ordercombs = cat(n+1, combs{:}); %// concat the n n-dim arrays along dimension n+1combs = reshape(combs,[],n);

Function for combvec solution:

function combs = f2(vectors)combs = combvec(vectors{:}).';

Script to measure time by calling timeit on these functions:

nn = 20:20:240;t1 = [];t2 = [];for n = nn;    %//vectors = {1:n, 1:n, 1:n};    vectors = {1:n/10, 1:n, 1:n*10};    t = timeit(@() f1(vectors));    t1 = [t1; t];    t = timeit(@() f2(vectors));    t2 = [t2; t];end