Sometimes in programming you need to do a pairwise comparison of some elements coming from two collections, for example, checking possible collisions between particles (which may be embedded inside a quadtree representation for efficiency). A handy operation is then the Cartesian product of the two sets of elements, to get the set of all pairs, which can then be traversed.

Whenever I need a Cartesian product of two lists in Haskell, I whip out the list monad to generate the Cartesian product:

cartesianProduct :: [a] -> [b] -> [(a, b)] cartesianProduct as bs = as >>= (\a -> bs >>= (\b -> [(a, b)]))

Thus for every element `a`

in `as`

we get every element `b`

in `bs`

and form the singleton list of the pair `(a, b)`

, all of which get concatenated to produce the result. For example:

*Main> cartesianProduct [1..4] [1..4] [(1,1),(1,2),(1,3),(1,4),(2,1),(2,2),(2,3),(2,4),(3,1),(3,2),(3,3),(3,4),(4,1),(4,2),(4,3),(4,4)]

This traverses the Cartesian space in row order. That is, if we imagine the square grid of possibilities here (where the element at the ith row and jth column corresponds to element `(i, j)`

), then `cartesianProduct`

generates the pairs in the following order:

In mathematics, the Cartesian product is defined on sets, so the order of these pairs is irrelevant. Indeed, if we want to examine all pairs, then it may not matter in what order. However, if we learn something as we look at the pairs (i.e., accumulating some state), then the order can be important.

In some recent research, I was building an automated algebra tool to find axioms for some algebraic structure, based on an interpretation. This involved generating all pairs of syntax trees `t`

and `t'`

of terms over the algebraic structure, up to a particular depth, and evaluating whether `t = t'`

. I also had some automated proof search machinery to make sure that this equality wasn’t already derivable from the previously generated axioms. I could have done this derivability check as a filtering afterwards, but I was exploring a very large space, and did not expect to even be able to generate all the possibilities. I just wanted to let the thing run for a few hours and see how far it got. Therefore, I needed Cartesian product to get all pairs, but the order in which I generated the pairs became important for the effectiveness of my approach. The above ordering (row major order) was not so useful as I was unlikely to find interesting axioms quickly by traversing the span of a (very long) row; I needed to unpack the space in a more balanced way.

My first attempt at a more balanced Cartesian product was the following:

cartesianProductBalanced :: [a] -> [b] -> ([(a, b)]) cartesianProductBalanced as bs = concatMap (zip as) (tails bs) ++ concatMap (flip zip bs) (tail (tails as))

`tails`

gives the list of successively applying `tail`

, e.g., `tails [1..4]`

= `[[1,2,3,4],[2,3,4],[3,4],[4],[]]`

This definition for `cartesianProductBalanced`

essentially traverses the diagonal of the space and then the lower triangle of the matrix, progressing away from the diagonal to the bottom-left, then traversing the upper triangle of the matrix (the last line of code), progressing away from the diagonal to the top-right corner. The ordering for a 4×4 space is then:

*Main> cartesianProductBalanced [1..4] [1..4] [(1,1),(2,2),(3,3),(4,4),(1,2),(2,3),(3,4),(1,3),(2,4),(1,4),(2,1),(3,2),(4,3),(3,1),(4,2),(4,1)]

This was more balanced, giving me a better view of the space, but does not scale well: traversing the elements of this Cartesian product linearly means first traversing all the way down the diagonal, which could be very long! So, we’ve ended up with a similar problem to the row-major traversal.

Instead, I finally settled on a “tiling” Cartesian product:

cartesianProductTiling :: [a] -> [b] -> [(a, b)] cartesianProductTiling [] _ = [] cartesianProductTiling _ [] = [] cartesianProductTiling [a] [b] = [(a, b)] cartesianProductTiling as bs = cartesianProductTiling as1 bs1 ++ cartesianProductTiling as2 bs1 ++ cartesianProductTiling as1 bs2 ++ cartesianProductTiling as2 bs2 where (as1, as2) = splitAt ((length as) `div` 2) as (bs1, bs2) = splitAt ((length bs) `div` 2) bs

This splits the space into four quadrants and recursively applies itself to the upper-left, then lower-left, then upper-right, then lower-right, e.g.,

*Main> cartesianProductTiling [1..4] [1..4] [(1,1),(2,1),(1,2),(2,2),(3,1),(4,1),(3,2),(4,2),(1,3),(2,3),(1,4),(2,4),(3,3),(4,3),(3,4),(4,4)]

Thus, we explore the upper-left portion first, without having to traverse down a full row, column, or diagonal. This turned out to be much better for my purposes of searching the space of possible axioms for an algebraic structure.

Note that this visits the elements in Z-order (also known as the Lebesgue space-filling curve) [Actually, it is a reflected Z-order curve, but that doesn’t matter here].

The only downside to this tiling approach is that it does not work for infinite spaces (as it calculates the length of `as`

and `bs`

), so we cannot exploit Haskell’s laziness here to help us in the infinite case. I’ll leave it as an exercise for the reader to define a tiling version that works for infinite lists.

Of course, there are many other possibilities depending on your application domain!

Later edit: In the case of the example I mentioned, I actually do not need the full Cartesian product since the order of pairs was not relevant to me (equality is symmetric) and neither do I need the pairs of identical elements (equality is reflexive). So for my program, computing just the lower triangle of Cartesian product square is much more efficient, i.e,:

cartesianProductLowerTriangle :: [a] -> [b] -> [(a, b)] cartesianProductLowerTriangle as bs = concatMap (zip as) (tail (tails bs))

However, this does not have the nice property of tiling product where we visit the elements in a more balanced fashion. I’ll leave that one for another day (I got what I needed from my program in the end anyway).

Expanding on my comment on twitter:

1 2 4 7

3 5 8 11

6 9 12 14

10 13 15 16

(or its transpose) is a much fairer interleaving in that it handles the infinite cases rather nicely.

You can think of it as a form of diagonalization.

For reference, http://hackage.haskell.org/package/control-monad-omega-0.3.1/docs/src/Control-Monad-Omega.html#diagonal implements this kind of search. Calling it a monad there gives me pause, but it does do the kind of thing you want.

Note how even your balanced version starves when given an infinite grid, because it has to reach all the way to the end before it starts to wander. Here we don’t.

Great! Thanks for that; very interesting. Yes the balanced version has the problem of chasing the (infinite) diagonal first in the infinite setting- bad! I will try that out in the place I was using it and see it if performs better (in terms of getting useful results faster).

Another option is to use LogicT from the logic package, which supplies explicit “fair” interleaving and bind operators.

The downside to a LogicT-style search is that because, at least effectively, the first item gets half of its attention and the second half of that, etc, you wind up with lots of searches at different stages of the process. This is great if you can stack-rank your choices at each turn, but don’t have exact probabilities.

In guanxi I go further and offer a monad for weighted non-deterministic choice, and then draw terms from that without replacement.

This is a great idea! I recently tried to integrate `testing-feat` with `hedgehog`: https://github.com/etorreborre/registry-hedgehog/blob/testing-feat/test/Test/Data/Registry/FeatSpec.hs. My idea was to start enumerating data structures before resorting to random sampling of values. I also stumbled on the fact that a naive cartesian product was not so nice in terms of “diversity” of the first generated values, so I implemented a “breadth-first” version of a cartesian product of n lists of elements instead of a “depth-first” one: https://github.com/etorreborre/registry-hedgehog/blob/testing-feat/src/Data/Registry/Internal/Feat.hs#L55. But now that makes me want to implement a “diagonal” version as well! I don’t know how hard it is to implement though, I already found my version hard to implement, my first attempt was both complicated and wrong :-).

Cool! Interesting to hear you stumbling on the same kind of idea elsewhere. Yes these things are a bit subtle. I still need to sit down and do the “Omega” version of the lower-triangle. I tried the other day and it was harder than expected (I kept getting it wrong).

Good to know I am not the only one struggling with this kind of code :-). And I suppose that writing the version which goes across n lists will be even more fun!

This is a great idea! I recently tried to integrate `testing-feat` with `hedgehog`: https://github.com/etorreborre/registry-hedgehog/blob/testing-feat/test/Test/Data/Registry/FeatSpec.hs. My idea was to start enumerating data structures before resorting to random sampling of values. I also stumbled on the fact that a naive cartesian product was not so nice in terms of “diversity” of the first generated values, so I implemented a “breadth-first” version of a cartesian product of n lists of elements instead of a “depth-first” one: https://github.com/etorreborre/registry-hedgehog/blob/testing-feat/src/Data/Registry/Internal/Feat.hs#L55. But now that makes me want to implement a “diagonal” version as well! I don’t know how hard it is to implement though, I already found my version hard to implement, my first attempt was both complicated and wrong :-).