Efficient method to sort very large set of integer vectors by all coordinates simultaneously

I have a set $E$ which is the set of all possible $d$-tuples ($d$-dimensional vectors) of integers between $1$ and $n$.

Typically $d=3$ and $n\approx1000$, but for the sake of making a small example, suppose $d = 2$ and $n = 4$, so

$$E = \{(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)\}\,.$$

Given an arbitrary nonempty $J\subseteq E$, I would like to find the product of the lengths of its projections on the axes, in other words, the area of the minimal axis-aligned bounding box of the 2d points in $J$.

For example, given $J = \{(2,3),(3,2),(4,1),(4,2)\}$, my function would return

$$f(J) = (\max(2,3,4,4)-\min(2,3,4,4))(\max(3,2,1,2)-\min(3,2,1,2)) = (4 - 2)(3 - 1) = 4\,.$$

I want to perform this calculation for every possible non-empty $J\subseteq E$, of which there are $2^{16} - 1 = 65,535$. At present I call built-in "Max" and "Min" commands every time, i.e. four such calls for each evaluation of $f$, hence $65,535\times 4 = 262,140$ calls. Presumably each time Max or Min is called, it sorts the positive integers fed to it.

This seems horribly inefficient, since for $J_1$ and $J_2$ having many elements in common, the lists to be sorted for $J_1$ are very similar to the lists to be sorted for $J_2$, so there may be great repetition of comparisons, even if the built in Max or Min function is efficient within itself.

What is a good way to solve this problem which achieves a balance between efficiency of time and memory, given the typical parameters $d = 3$ and $n = 1000$?

Should I store the nonempty subsets $J$ of $E$ in some kind of graph structure, and have my algorithm traverse this graph somehow?

EDIT:

So what I’m actually trying to do is to evaluate

$R(s_1, \ldots, s_d, n_1, \ldots, n_d; q) = 1 + \sum_{J \in \mathcal{P}(E) \setminus \emptyset} (-1)^{|J|} \prod_{J' \in \mathcal{P}(J) \setminus \emptyset} \exp \left [ (-1)^{|J'|} \ln \left ( \frac{1}{q} \right ) \prod_{r=1}^d \max \left ( 0, s_r - \left ( \max_{e \in J'} e_r - \min_{e \in J'} e_r \right ) \right ) \right ],$

where $d \in \mathbb{N}$, $n_1, \ldots, n_d \in \mathbb{N}$, and $s_1, \ldots s_d \in \mathbb{N}$ with $1 \leq s_r \leq n_r$ for all $r \in \{1, \ldots, d\}$, and $E = \prod_{r=1}^d \{1, \ldots, n_r - s_r + 1 \}$, and $q$ is a symbol. This expression gives a polynomial in $q$. The sequence of coefficients of the polynomial is what I would like to obtain.

This formula itself might well not be a very good formula for computing this polynomial, but my present task is to make an algorithm based upon the formula.

The quantity

$\prod_{r=1}^d \max \left ( 0, s_r - \left ( \max_{e \in J'} e_r - \min_{e \in J'} e_r \right ) \right )$

represents the volume of the intersection of several contiguous subarrays of dimensions $s_1 \times \ldots \times s_d$, within a large array of dimensions $n_1 \times \ldots \times n_d$. The reason for taking the maximum of 0 and $s_r - \left ( \max_{e \in J'} e_r - \min_{e \in J'} e_r \right )$ is that some such intersections will be empty.

Re. Tom’s comment about enumerating bounding boxes instead of subsets, I did have the following idea, which I think might be in that spirit:

Since the volume of the intersection of several contiguous subarrays is invariant with respect to translation of that set of subarrays within the large array, I could choose only to compute the volumes of a certain number $M$ of representative cases, and then copy the result to all the translates. I think

$M= 2^ { \left ( \prod_{r=1}^d \min(s_r, n_r – s_r + 1) \right ) -1 }$,

which I computed by supposing that the “top left” contiguous subarray was included in the set of subarrays, and asking which other contiguous subarrays could intersect that one.

EDIT:

I do apologise. In my original question, $f(J)$ should not be the product of the lengths of the projections of the subset on the axes, it should simply be the sequence of lengths of the projections on the axes. As can be seen in my first edit, above, I do not compute the product of those lengths $l_r$, rather, the product $\prod_{r=1}^d \max(0,(s_r - l_r))$. When the $n_r$ are large compared with the $s_r$, many of those products $\prod_{r=1}^d \max(0,(s_r - l_r))$ will be zero, which should allow me to check many fewer than the $2^{(1000^3)}$ cases mentioned in Tom's comment.

• For $d=3$ and $n=1000$ there are $1000^3$ elements, and you want to do something with all possible subsets (of which there are, well, $2^{(1000^3)}$). It seems that the inefficiency of something is the last thing you'd want to worry about. Have you considered enumerating all possible bounding boxes, and for each bounding box counting the number of subsets that would give rise to that bounding box? You could presumably do that more efficiently. – Tom van der Zanden Sep 29 '15 at 15:22
• A sane max or min implementation definitely doesn't sort its inputs: sorting takes time $n\log n$ but you can determine the max (min) of a list in linear time by checking each element against the biggest (smallest) you've found so far. – David Richerby Sep 29 '15 at 16:06
• I think there's a massive X-Y problem here. You've decided that the solution to some problem is to compute bounding boxes for $2^{n^d}$ sets of points and you're asking how to do that efficiently. That can't be done, since $2^{n^d}$ is too big. Why do you think you want to calculate all those bounding boxes? – David Richerby Sep 29 '15 at 16:09
• Dear @TomvanderZanden and David Richerby thank you both very much for your great answers ! Thank you to David for fixing my post, too, and for explaining about max/sort. Tom you make an excellent point - I suppose as stated my calculation would take astronomically long. I have thought of something similar to what you suggest, but to explain it I must first confess to the massive X-Y sin of which David rightly accuses me ! I had thought it better to provide a minimal example, but having read the FAQ (thank you for pointing me that way, David) I see that its better to tell you the whole story. – Simon Sep 29 '15 at 18:50
• I have now edited my original question to explain more about what I am trying to do, i.e. to try to resolve the X-Y problem. – Simon Sep 29 '15 at 18:57