Is the Ordering of the Output from combinations() Documented, Repeatable, and Sensible?

Suppose I want to evaluate a function of three variables over a grid of points in 3D space.

For example, using three nested loops

x = 1:2; y = 10:11; z = 100:102;
f = @(x,y,z) x*y*z;
for ii = 1:numel(x)
    for jj = 1:numel(y)
        for kk = 1:numel(z)
            a(ii,jj,kk) = f(x(ii),y(jj),z(kk));
        end
    end
end

However, I'd rather use a single @doc:parfor loop.

One approach is to @doc:ndgrid the inputs, then reshape to rows, evaluate, and reshape the result. Or, preallocate because we know the size and ordering of the output and use linear indexing

[X,Y,Z] = ndgrid(x,y,z);
c = [X(:),Y(:),Z(:)];
aa = nan(size(X));
for ii = 1:height(c)  % would be parfor in reality
    aa(ii) = f(c(ii,1),c(ii,2),c(ii,3));
end
isequal(aa,a)

ans = logical
   1

The "modern" approach avoids forming the temporary variables X,Y, and Z and instead uses @doc:combinations and then iterate over the rows. I could try the same approach as above, but I can't find a deifnitive statement on ordering of the output of combinations(), so don't know the dimensions needed for preallocation. Instead, preallocate and fill in a 1D vector for starters

c = combinations(x,y,z);
aaa = nan(height(c),1);
for ii = 1:height(c) % would be parfor in reality
    aaa(ii) = f(c{ii,1},c{ii,2},c{ii,3});
end

Now i want to @doc:reshape aaa back to the same shape as the underlying 3D grid, which would be most naturally represented in @doc:ndgrid form, so I can index into the reshaped aaa to extract the values of the function.

The first questions are: what is the ordering of the rows of c, and is that ordering guaranteed to be the same for all calls to combinations()? As far as I can tell, doc@combinations is silent on these issues.

For this particular problem, the rows of c are ordered in reverse ndgrid

[Z,Y,X] = ndgrid(z,y,x);
[X(:),Y(:),Z(:)],table2array(combinations(x,y,z))

ans = 12×3

1 10 100 1 10 101 1 10 102 1 11 100 1 11 101 1 11 102 2 10 100 2 10 101 2 10 102 2 11 100

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = 12×3

1 10 100 1 10 101 1 10 102 1 11 100 1 11 101 1 11 102 2 10 100 2 10 101 2 10 102 2 11 100

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now I need to @doc:permute after the reshape to recover the function in the desired form, which is a bit counterintuitive (and I shouldn't have had to experimentally determine the ordering of c to figure it out)

aaa = permute(reshape(aaa,size(X)),ndims(X):-1:1);
isequal(aaa,a)

ans = logical
   1

Suppose aaa is a very large array. Does that call to permute incur a lot of additional overhead?

I suppose that if I knew a priori how combinations() ordered it's output I could preallocate to the correct size, fill in with linear indexing, and then there would be no need to reshape. But I don't see an easy way to get around the need to permute.

1. Is the ordering of the output rows from combinations documented?

2. Will the ordering be the same on every call to combinations?

3. Is the ordering shown for this problem sensible? Frankly, I was quite surprised; I was expecting the ordering to be consistent with ndgrid(x,y,z). I wonder what led to this design choice.

These questions apply regardless of the data types of the inputs to combinations() and the output of f(); numeric just shown here for simplicity.

Answer 1

[Stephen23]

1. "Is the ordering of the output rows from combinations documented?"

Implicitly.

The COMBINATIONS description states "These combinations are the same as the Cartesian product...":

https://en.wikipedia.org/wiki/Cartesian_product

Although sets have no order, in practice the common, naive implentation of the cartesian product defines mathematical tuples which iterate fastest over the last operand's values and slowest over the first operand's values, which is exactly what COMBINATIONS delivers. This order is clearly illustrated in the examples given on Wikipedia, for example a deck of playing cards is the cartesian product of

• Ranks × Suits = {(A, ♠), (A, ♥), (A, ♦), (A, ♣), (K, ♠), ..., (3, ♣), (2, ♠), (2, ♥), (2, ♦), (2, ♣)}
• Suits × Ranks = {(♠, A), (♠, K), (♠, Q), (♠, J), (♠, 10), ..., (♣, 6), (♣, 5), (♣, 4), (♣, 3), (♣, 2)}

or the cartesian product of some numbers:

• A × B = {1,2} × {3,4} = {(1,3), (1,4), (2,3), (2,4)}
• B × A = {3,4} × {1,2} = {(3,1), (3,2), (4,1), (4,2)}

2. "Will the ordering be the same on every call to combinations?"

Yes, see the code.

3. "Is the ordering shown for this problem sensible?"

It matches the common naive implementations of the cartesian product that I can find. For example:

• Python's ITERTOOLS.PRODUCT
• Mathematica's TUPLES
• WolframAlpha

So even though set theory does not require this, in common practice the cartesian product is returned in https://en.wikipedia.org/wiki/Lexicographic_order of the operand element indices (which is the same order returned by SORTROWS). This is also the order used by various educational websites that my AI deepsearch tool found when searching online, advising me that "its effectively the de-facto standard for listing A×B in textbooks and tools" and also "across programming languages, math software, and educational platforms, Cartesian products are listed in consistent lexicographic order by default". It would therefore be rather odd/awkward for MATLAB to use another order for the cartesian product than all(?) other tools and resources.

Perhaps the function should have been named CartesianProduct.

[Stephen23]

Just for fun:

R = ["A","K","Q","J","10",num2cell('9':-1:'2')];
S = ["♠","♥","♦","♣"];
combinations(R,S)

ans = 52x2 table
     R      S 
    ___    ___

    "A"    "♠"
    "A"    "♥"
    "A"    "♦"
    "A"    "♣"
    "K"    "♠"
    "K"    "♥"
    "K"    "♦"
    "K"    "♣"
    "Q"    "♠"
    "Q"    "♥"
    "Q"    "♦"
    "Q"    "♣"
    "J"    "♠"
    "J"    "♥"
    "J"    "♦"
    "J"    "♣"

combinations(S,R)

ans = 52x2 table
     S      R  
    ___    ____

    "♠"    "A" 
    "♠"    "K" 
    "♠"    "Q" 
    "♠"    "J" 
    "♠"    "10"
    "♠"    "9" 
    "♠"    "8" 
    "♠"    "7" 
    "♠"    "6" 
    "♠"    "5" 
    "♠"    "4" 
    "♠"    "3" 
    "♠"    "2" 
    "♥"    "A" 
    "♥"    "K" 
    "♥"    "Q" 

[Stephen23]

Summary

NDGRID is for computation on grids (F-order / MATLAB-style):

• Purpose: designed for numerical computing — evaluating functions over a structured N-D grid (e.g., for interpolation, PDEs, image processing).
• Order rationale: uses Fortran-style (F-order, "column-major") layout — the first dimension varies fastest when iterating. This aligns with MATLAB's memory layout, where columns (1st index) are stored contiguously.
• Efficient memory access when looping through arrays in MATLAB.
• Grid coordinates are generated in a way that matches how MATLAB stores and traverses arrays in memory.

COMBINATIONS is for enumeration / logic (math-based):

• Purpose: built for enumerating all possible tuples from multiple sets — used in combinatorics, logic, and table construction.
• Order rationale: follows mathematical lexicographic convention — the last dimension varies fastest.
• Matches other languages, tools, learning resources.
• Matches the order of nested loops.

[Paul]

When I read that Wikipedia page before posting I was concerned that, as you say, "sets have no order." I was unaware that there is a de facto standard.

The order should be explicitly documented. Perhaps as you explained it in a single sentence.

I thought that one purpose of COMBINATIONS is to refactor multiple nested loops into a single loop to operate on all possible tuples. I suppose it can be used for that purpose, but maybe other approaches are better.

[Steven Lord]

@Paul Please submit this feedback about the information you feel is missing from the documentation to Technical Support or submit it using the star rating section at the end of the combinations function documentation page (after you give a star rating you will be able to provide text feedback.)

@Stephen23 "Perhaps the function should have been named CartesianProduct."

As someone who has been involved in a LOT of discussions / negotiations / arguments / brainstorming sessions about function naming here at MathWorks, choosing the right name (descriptive, discoverable, "not scary", accurate, etc.) can be difficult and time consuming. [And we don't always succeed at satisfying all of those considerations. The first example I think someone will mention is bsxfun. We spent a long time arguing about that name!]

I don't remember if we specifically discussed cartesianProduct as a name when designing this function, but would most users who were looking to generate combinations of elements from vectors necessarily think to look for "Cartesian"? I'm not so sure. [Even if they did, because that phrase is used in the Description section of the documentation page, you can find it with a search.] Would they look for something like "how to create combinations of elements"? That seems more likely IMO.

Answer 2

[Walter Roberson]

The source code for combinations is readable.

The output order is not documented.

The calculation is determinstic, so the order will be the same on every call. (Well, unless the undocumented algorithm changes...)

[Paul]

Should it be documented?

Is it sensible?