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

Question

This question pops up quite often in one form or another (see for example here or here). So I thought I'd present it in a general form, and provide an answer which might serve for future reference.

Given an arbitrary number n of vectors of possibly different sizes, generate an n-column matrix whose rows describe all combinations of elements taken from those vectors (Cartesian product) .

For example,

vectors = { [1 2], [3 6 9], [10 20] }

should give

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 ]

Answer 1

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 vectors

n = numel(vectors); %// number of vectors
combs = 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+1
combs = reshape(combs,[],n); %// reshape to obtain desired matrix

Answer 2

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{:}).')

Answer 3

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 vectors
combs = 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+1
combs = 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