Maximizing the area of intersection of rectangles
D.B.A. Epstein∗

M.S. Paterson†

24 March 2017

arXiv:1703.08658v1 [] 25 Mar 2017

Abstract
This paper attacks the following problem. We are given a large number N of
rectangles in the plane, each with horizontal and vertical sides, and also a number
r < N . The given list of N rectangles may contain duplicates. The problem is
to find r of these rectangles, such that, if they are discarded, then the intersection
of the remaining (N − r) rectangles has an intersection with as large an area as
possible. We will find an upper bound, depending only on N and r, and not on the
particular data presented, for the number of steps needed to run the algorithm on
(a mathematical model of) a computer. In fact our algorithm is able to determine,
for each s ≤ r, s rectangles from the given list of N rectangles, such that the
remaining (N − s) rectangles have as large an area as possible, and this takes
hardly any more time than taking care only of the case s = r. Our algorithm
extends to d-rectangles—analogues of rectangles, but in dimension d instead of in
dimension 2. Our method is to exhaustively examine all possible intersections—this

is much faster than it sounds, because we do not need to examine all Ns subsets in
order to find all possible intersection rectangles. For an extreme example, suppose
the rectangles are nested, e.g., concentric squares of distinct sizes, then the only
intersections examined are the smallest s + 1 rectangles.

1

Background

This problem arose from the use of a novel robotically controlled microscope technology
(Toponome Imaging System, abbreviated to TIS ), invented by Walter Schubert (see [4]).
TIS records, pixel by pixel, the location and abundance of proteins in a tissue section,
thus co-locating many different proteins in the same pixel. Since proximity is necessary
(though not sufficient) for interaction, this gives a powerful new method for discovering
protein complexes and protein networks, particularly since the anatomy of the section
is not disrupted by the TIS process. A TIS run results in a large number (typically,
hundreds) of different images of the same tissue section.
Such sequences of images may issue from the TIS machine approximately aligned,
but they will seldom be completely correctly aligned without further processing, because
physical reality, for example changes in temperature or vibrations from a passing heavy
goods vehicle, ensures that camera measurements are never perfect. In order to obtain
good information about co-location of proteins and other biomolecules, alignments that
∗

Mathematics Institute, University of Warwick. Partial support from BBSRC Grant BB/K018868/1.
Computer Science Department, University of Warwick. Partial support from Centre for Discrete
Mathematics and its Applications (DIMAP).
†

1

s=2
s=0
s=1

Figure 1: N = 6, r = 2. Shading shows maximal intersections for s = 0, 1, 2.

are as accurate as possible must first be achieved. In [3], it is explained how this can be
done. The only adjustments necessary to achieve good alignments turn out to be rigid
translation of the different images against one fixed image.
After alignment, some pixels will appear in all N images, and others may appear in
only (N − 1) images, or in (N − 2) images, and so on, as shown in Figure 1. If we insist
that a pixel can be analyzed only if it records all these signals simultaneously, then we
have to restrict to the region where all N images overlap. However, if one is prepared to
discard completely one or more images, then the set of pixels which figure in all remaining
images will be larger, much larger if one is discarding very seriously misaligned images.
In practice one would normally decide to jettison poorest quality images first, so the
algorithm presented here, though useful, is not the final word.
For a similar scenario in more conventional photography, suppose one wants an image
of St. Mark’s Basilica in Venice, but without the tourists and pigeons. To achieve this,
one might take snapshots periodically from the same spot over an extended period. The
tourists and pigeons can then be eliminated by aligning the images and then averaging
the results over each pixel in the intersection of the images. Using our algorithm one can
quickly determine which images to discard, in order to obtain a final image of desired
area.

2

Statement of result

We start by fixing our notation and clarifying the meaning of the terms we use. In
general, we will work in dimension d. In the discussion above, it was assumed that d = 2,
but our algorithm will work well in general.
Definition 2.1. Given an ordered pair (λ, ρ) of reals numbers, we can form the space
{x : λ ≤ x ≤ ρ} ⊆ R. If λ < ρ, this is a closed interval; if λ = ρ, it is a point; and if λ > ρ
it is the empty set. For the purposes of this paper, it is convenient to have a concept that
includes the endpoints λ and ρ as well as the underlying space. Refraining from pedantic
rigour, we will use the word interval as though it is nothing more the underlying space,
but nevertheless, we will feel free to extract the “endpoints” whenever necessary, even
when λ > ρ. We will denote such an “interval” by [λ, ρ], even if λ > ρ.
Definition 2.2. By a d-rectangle, we mean a product
[λ1 , ρ1 ] × · · · × [λd , ρd ] ⊆ Rd .
2

With the same level of informality as in Definition 2.1, we are able to extract the 2d
endpoints, even when the underlying space is the empty set. Given a d-rectangle R as
above, suppose that λi < ρi for exactly k values of i, and that λi = ρi for the other values
of i. Then we call the R a k-dimensional rectangle and define dim(R) = k.
Note that a d-rectangle may or may not be a d-dimensional rectangle. In general the
underlying space of a d-rectangle may be the empty set or a point or a k-dimensional
rectangle in Rd . A d-rectangle can be presented as a 2d-dimensional row vector


R = λ1 ρ1 . . . λd ρd .
(1)
Definition 2.3. Given a d-rectangle R as above, we extend the definition of its volume
slightly. The d-volume of R is a pair (dim(R), vol(R)) where:
1. If, for each i, λi < ρi ), then
dim(R) = d and vol(R) = (ρ1 − λ1 ) × · · · × (ρd − λd ) > 0.
2. If d > k > 0 and R is a k-dimensional rectangle in Rd (as in Definition 2.2), then
dim(R) = k and vol(R) is equal to the volume of R, regarded as a k-rectangle in
Rk . For example, if k = 1, then R is an interval in Rd and vol(R) is its usual
length.
3. If R is a 0-rectangle (that is, a point), then (dim(R), vol(R)) = (0, 0).
4. If, for some i, λi > ρi , so that the underlying point set of R is ∅, then
(dim(R), vol(R)) = (−1, 0).
If d = 1, this gives the length of the interval, and if d = 2 the area of the rectangle. If d-volumes need to be ordered, then we use lexicographical ordering on the pair
(dim(R), vol(R)).
Convention 2.4. Throughout this paper (unless otherwise stated) d > 0 will denote the
dimension of the euclidean space Rd in which we are working, and R will denote a specific
ordered list of N d-rectangles. This list can equally well be represented by an N × 2d
real matrix. We fix an integer r with 0 < r < N .
Our result can now be stated. Its proof will occupy most of the rest of this paper.
Theorem 2.5. We are given data as in Convention 2.4. Then there is an algorithm
taking at most
 


r + 2d
O d
+ dN
2d
steps, that determines, for each s with 0 < s ≤ r, exactly s entries to delete from the list
R, so that the list of the remaining (N − s) d-rectangles, represented by the remaining
(N − s)
 elements of R, has an intersection whose d-volume is maximal, as we vary over
N
the s possible choices of s elements of R. The constant in our bound for the number of
steps is independent of N , r and d.
Corollary 2.6. For any fixed dimension d, the running time of the algorithm is O(r2d +
N ).
3

In applications such as that of microscopy, which was our primary motivation, many
rectangles may coincide in some coordinates. Our algorithm is designed to deal efficiently
with this situation. Several of the complications result from this.
We have been thinking mainly of the case when r is small compared with N . For r
near N , it is possible to achieve better results with a brute force search. For example, if
r = N − 1, then one deletes all except one d-rectangle. So one need only compute the
volume of each of the given d-rectangles, and one then takes the maximum of these. This
takes O(dN ) steps. We will not bother to make explicit other similar brute force results
when r is nearly as big as N .
To be explicit, we need a model of computation: we will assume the usual model of
a RAM machine. This means that the number of memory cells is as large as we need for
the computation, and that each of these is accessed in unit time. Moreover, memory cells
are assumed to be large enough so that any particular integer that needs to be stored
can be stored in a single cell. Formally, we do not work with floating point numbers. If
necessary, these can be converted to integers by a suitable change of scale. However, it is
convenient and less confusing to pretend that we can work with arbitrary real numbers.
Arithmetical operations including simple comparisons are assumed to take unit time.
This gives a reasonably realistic mathematical model of a computer for most practical
purposes.

3

Some notations and definitions

We can think of a d-rectangle as an element in R2d , a row vector as in Equation 1. R is
then a map R : {1, · · · , N } → R2d . In particular, repetition is allowed.
Convention 3.1. Formally speaking, a list is a map, like R. However, it is convenient
to use the same notation also for the image of the map. The advantage of this abuse of
notation is that it allows us to write R ∈ R for any element R in the list. This paper
will contain several lists, and we allow the same abuse of notation in each case. We will
disambiguate if we think there is a possibility of confusion.
For 1 ≤ j ≤ N , we write Rj instead of R(j), and, following Equation 1, we will
consistently use the notation:.

 

Rj = λj,1 ρj,1 . . . λj,d ρj,d = rj,1 · · · rj,2d
(2)
Definition 3.2. We often refer to elements of {1, · · · , N } as rectangle indices, to remind
ourselves that they are indices for our list R of d-rectangles.
Definition 3.3. Let L ⊆ {1, · · · , N } be non-empty. The intersection of the family
(Rj )j∈L is defined as follows. For each i with 1 ≤ i ≤ d, let Mi = maxj∈L λji and let
mi = minj∈L ρji . Then the intersection is the d-rectangle


R(L) = M1 m1 . . . Md md .
(3)
Equation 3 is the result of a trivial calculation when R(L) is not empty, and is a definition
when the point set of R(L) is empty. R(L) is called an intersection rectangle. We denote
by Ir the set of non-empty intersection rectangles (R(L) 6= ∅), where L ⊆ {1, · · · , N }
and |L| ≥ N − r.
4

We will find an efficient way to produce, one by one, the elements of Ir . We will then
be able to calculate, for each R ∈ Ir , the d-volume of R and compare these results to
find the best one. The number of P
subsets like
 L above, giving rise to elements of Ir in
Definition 3.3, could be as large as rs=0 Ns , potentially a large number. To reduce this
number, note that a given element R ∈ Ir can in general satisfy R = R(L) for many
different L ⊂ {1, · · · , N }. For any such L, we define L(R), by
L(R) = {j|R ⊆ Rj } ,

(4)

which is the largest subset of the given set of rectangles with intersection R. Then
L ⊆ L(R) = L(R(L)), and |L(R(L))| ≥ |L| ≥ N − r.
We are thinking of r as small compared with N , and so it is easier to work with the
complement of L(R) in {1, · · · , N }. We define
U(R) = {j|R * Rj }

(5)

so that |U(R)| ≤ r.
Definition 3.4. We define
Ur = {U(R) : R ∈ Ir } .
The map U : Ir → Ur is clearly a bijection.
An essential feature of our solution will be the use of different orderings on the set
{1, · · · , N }, one for each i with 1 ≤ i ≤ 2d. To understand better why and how orderings
are important, we discuss the special case d = 1. In this case, we are given a list of N
intervals R = (I1 , . . . , IN ), with Ij = [λj , ρj ]. To avoid irrelevant detail, we assume for the
moment that the 2N endpoints of these intervals are all distinct, and that the intersection
of the N intervals is not empty. We define u, v ∈ {1, · · · , N } so that λu = maxj λj = m
and ρv = minj ρj = M , as in Definition 3.3. Unless Iu is removed from R, the lefthand
endpoint of the intersection will continue to be m. Unless Iv is removed from R, the
righthand endpoint of the intersection will continue to be M . So Iu or Iv must be
removed (and remember that u = v is a possibility) if the length of the intersection is to
strictly increase. If r ≥ 2, and we are omitting r intervals from R, then similar reasoning
will apply to subsequent removal of intervals. We have to look repeatedly for the largest
lefthand endpoint and, separately, for the smallest righthand endpoint, in the remaining
intervals.
Given a list R of N d-rectangles, we now define, for each i with 1 ≤ i ≤ 2d, orderings
≺i and 4i on {1, · · · , N }, modelling the definitions on the discussion in the previous
paragraph.
Definition 3.5. Let j, k ∈ {1, · · · , N } and let 1 ≤ t ≤ d. If i = 2t − 1, we define j ≺i k
if λj,t > λk,t , and j 4i k if λj,t > λk,t or if (λj,t = λk,t and j ≤ k). If i = 2t, we define
j ≺i k if ρj,t < ρk,t , and j 4i k if ρj,t < ρk,t or if (ρj,t = ρk,t and j ≤ k). Then ≺i is a
partial order and 4i is a total order on {1, · · · , N }. For odd i, the orderings are (weakly)
decreasing in the λ’s, and, for even i, they are (weakly) increasing in the ρ’s.

4

Properties of X ∈ Ur .

We want to produce the elements X ∈ Ur efficiently one by one. With this in mind,
let X ⊂ {1, · · · , N } with |X| ≤ r < N , in which case we find necessary and sufficient
conditions for X to be an element of Ur .
5

Definition 4.1. An endpoint list is a non-empty subset of {1, · · · , N }, arranged in order,
smallest first, with respect to one of the total orderings 4i described in Definition 3.5.
We set Ei equal to the set {1, · · · , N }, arranged in order according to 4i , so that Ei is
an endpoint list.
Definition 4.2. Let X ⊂ {1, · · · , N }, with |X| ≤ r < N . For each i, let j(i, X) be
the smallest element of Ei that is not in X (which must exist since X has fewer than N
elements). We define bv(i, X) = rj(i,X),i (where bv stands for barrier value). We set
bv(X) = (bv(1, X), · · · , bv(2d, X)).

(6)

Lemma 4.3. Suppose X ⊂ {1, · · · , N }, with |X| ≤ r < N . X ∈ Ur if and only if the
following two conditions are satisfied:
4.3.1: For each i with 1 ≤ i ≤ d, bv(2i − 1, X) ≤ bv(2i, X).
4.3.2: For each j ∈ X, there is an i such that j ≺i j(i, X).
Under these conditions, X = U(bv(X)).
Proof. Suppose first that X ∈ Ur , or equivalently that X = U(R) for some R ∈ Ir , where


∅=
6 R = r1 r2 . . . r2d−1 r2d .
For any j ∈
/ U(R), we have R ⊆ Rj , and so, for 1 ≤ i ≤ d,
rj,2i−1 ≤ r2i−1 ≤ r2i ≤ rj,2i .
This holds in particular with j = j(2i − 1, X) and with j = j(2i, X), giving rise to
bv(2i − 1, X) ≤ r2i−1 ≤ r2i ≤ bv(2i, X),

(7)

proving the first condition in the statement of Lemma 4.3.
Continuing with the hypotheses R ∈ Ir and X = U(R), we next show T
that bv(X) = R,
when we regard both as elements of R2d . Suppose 1 ≤ i ≤ 2d. Since R = {Rj : R ⊆ Rj },
we know that ri = rj,i , for some j such that R ⊆ Rj , that is, for some j ∈
/ X. By Equation 7, if i is odd, we have
ri = rj,i ≥ bv(i, X) = rj(i,X) .
But j(i, X) is, by definition, the smallest element not in X, and the ordering is reversed for
odd i, according to the ordering 4i (see Definition 3.5). So a strict inequality bv(i, X) <
rj,i is impossible, and we must have the equality bv(i, X) = ri . Similarly if i is even. This
proves that R = bv(X). Therefore X = U(R) = U(bv(X)).
We continue to assume that X = U(R), and let j ∈ X. Then R * Rj . It follows that,
for some i, [r2i−1 , r2i ] * [rj,2i−1 , rj,2i ]. This implies that r2i−1 < rj,2i−1 or r2i > rj,2i . We
have shown that R = bv(X), so this means that j ≺2i−1 j(2i − 1, X) or j ≺2i j(2i, X),
proving the second condition in the statement of Lemma 4.3.
Now let us assume the two conditions in the statement of the lemma, and prove that
X ∈ Ur . Let R = bv(X) be a d-rectangle. By Condition 4.3.1, R 6= ∅. If j ∈
/ X, then,
for each i, j(i, X) 4i j. Therefore bv(2i − 1, X) ≥ rj,2i−1 and bv(2i, X) ≤ rj,2i . It follows
that R ⊆ R
/ X and each i with 1 ≤ i ≤T2d, R ⊆ Rj(i,X) , so
Tj . In particular, for each j ∈
that R ⊆ i Rj(i,X)
T . On the other hand, Definition 3.3 shows that i Rj(i,X) ⊆ R. We
deduce that R = j ∈X
/ Rj , which shows that R ∈ Ir .
If, on the other hand, j ∈ X, then, by hypothesis, for some i, j ≺i j(i, X). Then
rj,i > bv(i, X) if i is odd, and rj,i < bv(i, X) if i is even. It follows that R * Rj .
Therefore X = {j : R * Rj } and X ∈ Ur .
6

5

Program paths

Before describing the algorithm that produces all the elements of Ur one by one, we
describe how to generate one particular U ∈ Ur . From Lemma 4.3, U = U(bv(U )).S
For 1 ≤ i ≤ 2d, we write Ui = {j|j ≺i j(i, U )}, so that, by Lemma 4.3, U = i Ui .
For 1 ≤ i ≤ 2d let
[
U i = Ui \
Uk .
(8)
1≤k<i

U is equal to the disjoint union of the U i . We generate first the elements of U 1 , then the
elements of U 2 , and so on.
We prefer to think of these rectangle indices in U as being discarded or deleted,
rather than generated, since we are thinking of generating the complement of U , L =
{1, · · · , N } \ U , and then taking the intersection of the d-rectangles in L. We start with
{1, · · · , N }, discard the rectangle indices in U 1 , then the rectangle indices in U 2 , and so
on. To completely specify the order in which the rectangle indices in U are discarded, we
insist that the elements of U i are discarded in increasing 4i order, that is, least first.
Definition 5.1. By a ∆-move we mean one of the deletions just mentioned. While
working with U i , ∆ deletes the 4i -smallest element of Ei remaining. Until all elements
of U i are deleted, this smallest element is an element of U i . While carrying out such
deletions, we say that the algorithm is focussed at i. An S-move moves from the focus
index i to (i + 1). In our process, there are only these two types of move, and the process
starts with focus at i = 1. ∆ and S are mnemonic for “delete” and “shift index by one”.
At the end of the process, a final S-move causes the process to stop.
During the algorithm the Ei keep on changing. In order to keep track of what is going
on, we introduce a time variable t which increases by one at each move. Let δ(t, U ) be
the number of ∆-moves and s(t, U ) the number of S-moves during the first t units of
time, so that t = δ(t, U ) + s(t, U ). At time t, the focus is at i = s(t, U ) + 1. Recall
from Definition 4.1 that Ei consists of the elements of {1, · · · , N } listed in ascending
order according to 4i . We write Ei (t, U ) for Ei at time t, when all the lists Ei (t, U )
have exactly the same set of N − δ(t, U ) entries, though the orders in which these entries
appear are, in general, different. When a ∆-move deletes an element of Ei (t, U ) during
focus at i (always the smallest element according to 4i ), then that element is also deleted
from all the Ek (t, U ), for 1 ≤ k ≤ 2d.
The process starts at t = 0 with focus at i = 1. For each k with 1 ≤ k ≤ 2d,
Ek (0, U ) = {1, · · · , N }. For 1 ≤ i ≤ 2d, let ti be the first time at which the focus moves
to i (so t1 = 0), and let t2d+1 be the time when the process ends. At time ti , the set
of
indices so far deleted from {1, · · · , N } and from the various E` , is equal to
S rectangle
k
1≤k<i U . From the definition of Ui immediately preceding Equation 8, we know that
U i = {j|j ∈ Ei (ti , U ) and j ≺i j(i, U )}

(9)

(see Convention 3.1 for the meaning of j ∈ Ei ), so that the 4i -smallest element of U i ,
which is the next rectangle number to be discarded (assuming U i is not empty), is also
the 4i -smallest element of Ei (t, U ). Therefore, |U i | successive ∆-moves will delete exactly
the elements of U i , and, following this, an S-move transfers focus to (i + 1), or, if i = 2d,
S causes the process to stop.

7

Rules 5.2. It is convenient to think of the process π(U ) just described, when U ∈ Ur ,
as a word W(U ) in the terms S and P
∆, read from left to right. Let u(i) = |U i | be the
number of elements in U i . Note that i u(i) ≤ r. We define
W(U ) = ∆u(1) S∆u(2) S · · · S∆u(2d) S.

(10)

We will now present several rules on words w in S and ∆ that every w = W(U ) has
to obey. Once the rules have all been made explicit, we will define Wr to be the set of
legal words, that is, the set of all words satisfying all the rules. We will then show that
the map FWU : Ur → Wr given by FWU (U ) = W(U ) is a bijection, so that the rules
will be not only necessary but also sufficient for a word w to be equal to W(U ) for some
U ∈ Ur . This will reduce our task to producing an algorithm that produces one by one
the set of all legal words.
We already have three rules, satisfied by any word w = W(U ) for some U ∈ Ur :
Rule A: w contains exactly 2d terms equal to S.
Rule B: w ends with an S.
Rule C: w contains at most r terms equal to ∆.
Rules A and B mean that each coordinate direction is considered and Rule C ensures
that at most r deletions are used.
Definition 5.3. Now suppose a word w in ∆ and S satisfies Rule A, Rule B and Rule C.
We obtain a process π(w), analogous to the process π(U ) (see first sentence of 5.2), that
reads w from left to right, taking action as appropriate for each term read. More precisely,
we proceed as follows. If t ≤ length(w), we factorize w = wt .vt , where the length of wt is
t. We denote by δ(t, w) the number of terms equal to ∆ in wt , and by s(t, w) the number
of terms equal to S in wt . After reading wt , the focus is at i = i(t, w) = s(t, w) + 1. In
particular, the focus is at i = 1 when t = 0. We denote by Ek (t, w) the list Ek after
the δ(t, w) deletions that occur while wt is being read. With the focus at i, let j be the
4i -smallest element of Ei (t, w). If the first term of vt is ∆, then j is deleted from each
Ek (t, w), giving the endpoint list Ek (t + 1, w). If the first term of vt is S, then focus is
transferred from i to (i + 1), and, for each k with 1 ≤ k ≤ 2d, Ek (t + 1, w) = Ek (t, w). Let
U (w) be the set of all the rectangle indices j deleted while w is read. We are interested
only when U (w) ∈ Ur , which is not in general the case for an arbitrary word w satisfying
only the three rules above.
We will now describe further rules that take account of coincident coordinates, and
will prove in Lemma 5.4 that these rules are satisfied by w = W(U ), whenever U ∈ Ur .
It will turn out that, conversely, if w satisfies all the rules we present, then there is a
unique U ∈ Ur such that w = W(U ).
Rule D: As rectangles are discarded, barrier values may change. A lefthand endpoint
barrier value will not increase and may decrease. A righthand endpoint barrier value will
not decrease and may increase. This rule will ensure that any rectangle whose removal
would change the barrier value at i must be removed when the focus is no later than
i. Let j(t, k, w) be the 4k -smallest element of Ek (t, w), and let bv(t, k, w) = rj(t,k,w),k .
Suppose that 2 ≤ i ≤ 2d and that ti is the first time at which the focus is at i. The
condition on w is that, if 1 ≤ k < i, then bv(t, k, w) is independent of t provided that
ti ≤ t ≤ length(w). We denote this constant value by bv(k, w).
8

Rule E: Fix i with 1 ≤ i ≤ 2d. Let A be an endpoint list with ordering 4i . By an
A-block we mean a non-empty set of entries in A of the form
B = {j|j ∈ A and rji = x}
for some x ∈ R. Any j ∈ A is contained in a unique A-block. Such a block can be
of length one, but blocks can also be of any length up to N , as would happen, for
example, if all of the d-rectangles in the list R were equal to each other. We denote
by Bk (t, w) the initial Ek (t, w)-block, and suppose t is such that the focus is at k. Let
b(t, w) = |Bk (t, w)|. Then b(t, w) > 0. The condition on w = wt .vt is that, if vt starts
with ∆, then δ(t, w) + b(t, w) ≤ r and vt starts with ∆b(t,w) . In other words, if vt starts
with a ∆, then all the elements of the initial block of Ek (t, w) must be deleted before the
focus increases.
Rule F: For each k with 1 ≤ k ≤ d, bv(2k − 1, w) ≤ bv(2k, w). This rule ensures that
only non-empty intersections are reached, and this can be checked as the focus leaves any
even integer.
Lemma 5.4. Let U ∈ Ur and w = W(U ). Then w satisfies all the rules above.
Proof. From Equation 10 we immediately see that Rule A, Rule B and Rule C are true.
To prove Rule D, recall that j(k, U ) is the 4k smallest element of Ek (0, w) \ U (see
Definition 5.3 for the meaning of Ek (t, w)). Since only elements of U are discarded while
w is being read, j(k, U ) is never discarded. Therefore, for all t, j(k, U ) is the 4k smallest
element of Ek (t, w) \ U . Therefore, for all t,
j(t, k, w) 4k j(k, U ).

(11)

bv(t, k, w) = rj(t,k,w),k = rj(k,U ),k = bv(k, U ),

(12)

We claim that, for t ≥ tk+1 ,

where the first and third equalities are definitions. For, if Equation 12 were false, then
Equation 11 would give j(t, k, w) ≺k j(k, U ). Recall from the definition just before Equation 8, that Uk = {j : j ≺k j(k, U )}, so that j(t, k, w) ∈ Uk . However, for t ≥ tk+1 , all
elements of Uk have been discarded, in particular j(t, k, w), This contradiction proves our
claim, which proves Rule D.
Now we prove Rule E. Choose k with 1 ≤ k ≤ 2d, and suppose the focus is at k
at time t. Let B be the initial block of Ek (t, w), and let x ∈ R be the real number
such that B = {j ∈ Ek (t, w)|rjk = x}. The initial ∆ in vt deletes some j ∈ B ∩ U k , so
that j ≺k j(k, U ). If k is even, this means that x = rj,k < rj(k,U ),k , and, if k is odd,
x = rj,k > rj(k,U ),k . It follows that, for all j ∈ B, j ≺k j(k, U ). But this means that
B ⊆ Uk . Since the focus is on k, B ∩ U j = ∅ for j < k. We deduce that B ⊆ U k , and so
all elements of B are discarded immediately. This completes the proof of Rule E.
Rule F follows from Equation 12 and Lemma 4.3, and this completes the proof of
Lemma 5.4.
Proposition 5.5. The map FWU : Ur → Wr , defined by FWU (U ) = W(U ), is a bijection.
Proof. Lemma 5.4 shows that we do indeed have W(U ) ∈ Wr , for each U ∈ Ur .
Conversely, given w ∈ Wr , we obtain a subset U ⊂ {1, · · · , N } of discards by means
of the process π(w) (Definition 5.3). We need to prove the two conditions of Lemma 4.3.
9

Let Vi = {j : j ≺i j(i, U )}. By the definition of j(i, U ) (see Definition 4.2), Vi ⊆ U , so,
in due course, all elements of Vi are discarded during the process π(w). It follows that,
for i fixed and t sufficiently large (of course t ≤ r + 2d throughout), j(t, i, w) = j(i, U ),
and therefore bv(i, w) = bv(i, U ). Rule F implies 4.3.1.
Let j ∈ U . Then, for some i, j is discarded while the focus is at i. While the focus
is at i, only elements k ∈ U with k ≺i j(i, U ) are discarded. Therefore j ≺i j(i, U ). We
claim that rj,i 6= rj(i,U ),i . For otherwise j and j(i, U ) would be in the same Ei -block. But
then Rule E would imply that j(i, U ) is also discarded, implying that j(i, U ) ∈ U . This
contradiction implies that j ≺i j(i, U ), which is 4.3.2. We have proved the two conditions
of Lemma 4.3, and so U ∈ Ur .

6

A brief sketch of the algorithm

We have reduced our problem to that of finding an algorithm that produces the words
of Wr one-by-one. Consider the finite set W(r) of words w satisfying Rule A, Rule B
and Rule C. In the standard fashion, we form these words into a rooted tree T(r), with
directed edges, each edge labelled with either ∆ or S, so that there is exactly one path
in T from the root to a leaf for each word in W(r). There is a subtree Tr corresponding
to the words of Wr .
Our algorithm traverses a certain subtree of the tree T(r), in depth-first order, choosing to follow a ∆-edge rather than an S-edge, where there is a choice. For efficiency
reasons, we would prefer to arrange for this subtree to be exactly the subtree Tr of
interest. However, staying inside Tr would require a substantial amount of additional
computation—this is because Rule F cannot be easily checked term by term. Instead we
explore a subtree of T(r) that is somewhat larger than Tr . As soon as we discover that
we are outside Tr , we backtrack.
Our algorithm uses a recursive function call to perform a depth-first search of Tr . We
maintain a global data structure enabling us to check that the rules are being followed
during the search. We denote the recursive call by F (p, s). Here p is the maximum
number of times we permit further rectangle indices to be discarded from {1, · · · , N },
that is, the maximum number of terms ∆ that the algorithm may still generate, and p
is mnemomic for permissible discard. Also s is the number of endpoint lists still to be
examined, that is, the number of terms S still to be generated, and s is mnemonic for shift.
The behaviour of F (p, s) depends on the state of the global data structure, for example
deciding whether to call F (p − 1, s) or F (p, s − 1) or both. The depth-first aspect of the
search comes about since F (p, s) calls F (p − 1, s) before F (p, s − 1) if the data structure
tells us to call both. A principle followed by our algorithm is that F (p, s) alters the most
important part of the data structure before calling F (p − 1, s) and/or F (p, s − 1). On its
return, F (p, s) changes this part of the data structure back again before exiting. Thus,
(the most important part of) the global data structure, when F (p, s) starts, is the same
as the structure when F (p, s) finishes.
The depth-first search starts by calling F (r, 2d).

7

Selection versus sorting

We will later describe in detail the data structures used by our algorithm to check, as
rapidly as possible, that the rules are being followed. In this section we discuss a minor
10

improvement that can be made on the most obvious approach. The first three rules are
enforced by the notation F (p, s) and so have no implications for the data structure. The
procedure described above requires us to order {1, · · · , N } with respect to the 2d-different
orderings 4i described in Definition 3.5. This takes O(dN log N ) steps.
However, we can do a little better than this. For each k with 1 ≤ k ≤ 2d, let
Gk ⊆ {1, · · · , N } be the ordered set of the first r + 1 elements of {1, · · · , N } with respect
to the ordering 4k . It turns out (see Lemma 7.1) that our algorithm can be confined to
consideration only of the ordered sets Gk , and that it will then run in the same number
of steps (or slightly fewer). To compute the Gk we proceed as follows. For each k, find
the (r + 1)-st element pk ∈ Ek . An algorithm doing this is called a selection algorithm,
and it requires O(N ) steps, as described in the next paragraph. As an unordered set
Gk = {j : j 4k pk } is found in O(N ) steps. Gk is then sorted in O(r log r) steps. To
obtain all the Gk as ordered sets takes
O(d(N + r log r))

(13)

steps.
There is a randomized selection algorithm [2], due to Hoare, with an expected number
of steps equal to O(N ), though the worst-case bound is O(N 2 ). [1] contains an algorithm
where the number of steps is bounded by O(N ) steps even in the worst case, though
the constant hidden by the O-notation is not small. The latter algorithm only rarely
outperforms the Hoare algorithm, and the Hoare algorithm also has the advantage that
it can be carried out in place, that is, without needing more space than is needed by
the input data. There is a good article about this topic, giving Hoare’s code, at http:
//en.wikipedia.org/wiki/Selection_algorithm.
Given a word w in ∆ and S satisfying Rule A, Rule B and Rule C, we have the process
π(w) (see Definition 5.3). We change this slightly to the process πG (w), which is the same
as π(w), except that, when deletion of j from Ek is required in π(w), then no action is
taken unless j ∈ Gk .
Lemma 7.1. Rule D, Rule E and Rule F can be checked within the process πG (w).
Proof. We fix k, with 1 ≤ k ≤ 2d. By the time the process π(w) has completed, α
elements have been deleted from Gk , and β elements from {1, · · · , N } \ Gk , where α ≤
α + β ≤ r, so at most r elements have been deleted. Since |Gk | = r + 1, one or more
elements of Gk is never deleted during π(w). This means that we can compute j(t, k, w)
(defined in Rule D) within πG (w). So we can also compute bv(t, k, w). We are therefore
able to check Rule D and Rule F.
We assume now that w fails Rule E but passes all the other rules. We show that the
failure can be seen within the process πG (w). Using the notation in which Rule E was
stated, we suppose the rule fails at time t, when the focus is at k. We are assuming that
w = wt .vt , where wt has length t, and that vt deletes the first element of Bk (t, w), but
not all of it. This means that vt starts with ∆α .S, where 1 ≤ α < b(t, w). The number
of discards from Gk caused by wt is bounded above by δ(t, w), and δ(t, w) + α ≤ r, by
Rule C. Therefore w = wt .vt deletes no more than r elements from Gk , so that at least
one element remains in Gk . Since b(t, w) > α, the smallest such element must be in the
block Bk (t, w), so that we observe the failure of Rule E within Gk .

11

8

Data structure

We start with the given input: the N × 2d-matrix R as in 2.4. Given r with 0 < r < N ,
we want to generate all words in Wr . In particular, this means that we need the ability
to efficiently check the rules, avoiding unnecessary computation. We now describe the
different parts of the data structure. As stated above, significant parts of the data
structure, used and changed by F (p, s), are restored as F (p, s) exits. We will label these
with the label restore.
8.1 Tree: The tree T(r) of all words satisfying our first three rules was introduced in
Section 6. The subtree Tr ⊆ T(r) is formed from Wr , the set of words satisfying
all six of our rules. Tree is a subtree of T(r) that normally grows inside Tr as our
algorithm progresses, and various words are explored. From time to time we may
have to explore words outside Wr , because it may not be immediately apparent
that Rule F must inevitably fail, no matter how the word is extended. The error is
eventually discovered, and the illegal branch of Tree is snipped off.
Each node ν of the tree corresponds to a unique word w(ν) in the terms ∆ and S,
read from left to right as we descend the tree. Each node of the tree is provided
with data under the following headings:
8.2 Nodes and pointers ∆, S, parent: We establish two symbolic values UNDEF and
ILLEGAL. The node ν has pointers ν.∆, ν.S and ν.parent, each pointing to another
node of the tree or set equal to UNDEF or to ILLEGAL. As soon as we find out
that w(ν) is not a proper prefix of a word in Wr , we set ν.∆ = ILLEGAL, and
ν.S = ILLEGAL.
Initially, Tree has only one node, namely root, representing the null word, and
initially all its three pointers point to UNDEF. For every node ν the parent pointer
is never changed, so that the parent pointer of root always points to UNDEF. If
ν 6= root, then its parent pointer always points to a node, not to UNDEF, so this
can be used to recognize root.
8.3 CurrentNode (restore): is a global variable, pointing to the node ν of Tree
currently reached by the algorithm. It fixes the word w = w(ν) corresponding
to this node. Initially, CurrentNode points to root, the only initial node of Tree.
8.4 BSZ (restore) and Q (constant): Let Bi be the set of Gi -blocks, arranged in order
using ≺i . We set bi = |Bi |, and number the elements of Bi in order, using block
numbers {1, . . . , bi }, with 1 numbering the ≺i -smallest block and bi the ≺i -largest
block. Q is an N × 2d matrix, whose entry in position (j, i) gives the block number
of the Gi -block containing the rectangle index j, provided j ∈ Gi , and is equal to
UNDEF if j ∈
/ Gi . Constructing Q takes O(d N ) steps and Q does not change during
the algorithm. The changing block sizes are recorded in the (r + 1) × 2d-matrix
BSZ, such that BSZ(k, i) is the block size of the k-th block of Gi if k ≤ bi . If k > bi ,
then BSZ(k, i) = UNDEF. Initializing BSZ takes O(d r) steps. Updating BSZ when the
rectangle Rj is discarded or restored takes O(d) steps.
8.5 DLL (restore Doubly Linked List): In Section 7 we defined, for each i with 1 ≤
i ≤ 2d, Gi as the 4i -ordered set of the r + 1 rectangle indices that are smallest
with respect to 4i . DLL is an N × 2d-matrix, whose k-th column is denoted by
12

DLLk . Each entry of DLL is a pair of rectangle indices, that initially records Gi and
its order as follows. If j ∈
/ Gi , the entry at row j and column i is UNDEF and this
entry is unchanged throughout the algorithm. If j ∈ Gi , the (j, i) entry is initially a
pair (prev,next) of rectangle indices, specifying the ordering 4i on Gi in the usual
manner. We adjoin a row 0 to DLL to accommodate, for each i, a start entry at (0, i)
with prev = UNDEF and next equal to the 4i -smallest surviving element of Gi . If
j is the 4i -largest rectangle index surviving, then next = UNDEF at position (j, i).
The entries in DLL are changed in the usual way (taking constant time) as elements
j ∈ Gi are deleted. When j is deleted from Gi , we do not change the temporarily
meaningless values (prev,next) at position (j, i); instead these values are retained
for use when j is restored. DLL takes O(d (N + r log r)) steps to initialize. Deleting
or restoring a rectangle index j takes O(d) steps.
8.6 ν.vol: Let m be the number of terms in the word w = w(ν) that are equal to S. For
k ≤ m, bv(k, w) = rj,k , where j is the 4k -smallest element remaining in the doubly
linked list DLLk —this is because Rule D is enforced at each step. For 2i ≤ m, let
α(i) = bv(2i, w) − bv(2i − 1, w).
This is already known, unless m = 2i and the final term of the word w(ν) is S. If
α(i) < 0, we set ν.parent.S = ILLEGAL and snip off the part of the tree hanging
from ν. Otherwise set
Y
ν.vol =
{α(i)| 2i ≤ m and α(i) > 0} ,
which can be evaluated inductively, each step of the induction taking constant time.
8.7 ν. dim: We set ν.dim equal to the number of i such that α(i) > 0.
8.8 Results: For 0 ≤ p ≤ r, Results(p) has three subcomponents, namely dim, vol
and address. The first two subcomponents give the lexicographically largest result
of the form (dim,vol) so far achieved by discarding p d-rectangles from the list R,
then taking the intersection of the remaining d-rectangles. The third subcomponent
is the address of the Tree node νp that was current at the time that this best result
was found. Initially each dim entry is −1, each vol entry is UNDEF, and each address
entry is equal to UNDEF. Initialization takes O(r) steps.

9

Pseudocode

Here is pseudocode for the algorithm sketched in Section 6. We will use & and * for the
address and indirection operators, as in C. The address can be, as in C, an address in
memory, or, in languages that do not attempt to mimic machine architecture, the number
of a row in a table or matrix. The code starts with initialization, which consists of the
sorting processes described above in Section 8 and setting up the various data structures
in the obvious way. This takes O(dN + dr log r) steps.
The execution path of the program F (p, s), whose description follows, depends on the
state in which it starts. By the state we mean the values contained in the global data
structures described above. So F (p, s) stands for many different possible execution paths.
13

However, We will be able to give a reasonable upper bound for the number of steps that
F (p, s) requires, and this upper bound is independent of the particular execution path.
Rule A is enforced by the inequality 0 ≤ s ≤ 2d. We need to check that F (p, s)
enforces the other five rules.

9.1

F (p, s): deciding how the tree can be extended
Input: The global variables described above; p, an upper bound for the
number of characters ∆ still to be generated, or, equivalently, the number
of rectangles still to discard; s, the number of characters S still to be
generated, or, equivalently, 2d minus the number of characters S already
generated; CurrentNode, the address of the current node ν in the tree.
Output: Assignment of the pointers ν.∆ and ν.S to ILLEGAL, where appropriate. Storage of best results so far.
(1)

f ← 2d − s + 1; /* focus now at f */

(2)
(3)

j ← smallest element of DLLf ;
ν ← ∗CurrentNode; /* ν is the current node */

(4)
(5)

if s == 0
/* Exploration of this word is complete and all rules are
* satisfied. */
set ν.∆ = ν.S = ILLEGAL; /* Enforcing Rule B */

(6)
(7)
(8)
(9)
(10)

if Results(r − p).(dim, vol) < ν.(dim, vol) (lexicographically)
Results(r − p).((dim, vol) ← (ν.(dim, vol), &ν);
endif
else
/* s > 0. Check all rules, but without altering the main
* data structures. */

(11)

if p == 0
ν.∆ ← ILLEGAL; /* Enforcing Rule C */

(12)

else /* p > 0 */

(13)
(14)

if BSZ(Q(j, f ), f ) > p
/* Note that BSZ(Q(j, f ), f ) is the block containing
* j. When p = 0, the inequality is automatically sat* isfied, we have already discarded r rectangle indices,
* and further discards are not allowed. If p > 0 and the
* inequality is satisfied then Rule E would be violated
* by a further discard and this rules out continuations
* of w(ν) with next letter ∆. */

(15)
(16)

ν.∆ ← ILLEGAL
else

14

/* Rule E is satisfied and we check Rule D */
(17)
(18)
(19)
(20)

foreach i such that 1 ≤ i < f
k ← smallest element of DLLi ;
if Q(j, i) == Q(k, i) and BSZ(Q(j, i), i) == 1
/* j, k in same Gi -block and Rule D would be
* violated by a further deletion. */

(21)
(22)
(23)
(24)
(25)
(26)
(27)
(28)
(29)
(30)

ν.∆ ← ILLEGAL;
break from foreach loop
endif
endfor
endif
endif
if f is even
k ← smallest element of DLLf −1
if rk,f −1 > rj,f
/* Rule F would be violated by following with S. */

(31)
(32)
(33)
(34)

9.2

ν.S ← ILLEGAL;
endif
endif
endif

F (p, s): extending the tree, and recursion

At this stage we have determined all failures in the rules that would be
immediately apparent on following w(ν) with either ∆ or S. As a result,
we know which new nodes to construct, and we proceed with this task. We
also alter the global data structures as we explore the tree.
Input: The global variables; p, an upper bound for the number of characters ∆ still to be generated, or, equivalently, the number of rectangles still to
discard; s, the number of characters S still to be generated, or, equivalently,
2d minus the number of characters S already generated; CurrentNode, the
address of the current node ν in the tree.
Output: Construction of new nodes of the tree. Determination of the best
results possible so far, using recursion.
(1)

f ← 2d − s + 1; /* same as Line 1 of 9.1 */

(2)

j ← smallest element of DLLf ;

15

/* same as Line 2 of 9.1 */
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)

ν ← ∗CurrentNode
if ν.∆ 6= ILLEGAL
Create node α;
Treelocation ← &α;
ν.∆ ← &α;
α.parent ← &ν;
foreach i such that j ∈ DLLi
Delete j from DLLi ;
Decrease BSZ(Q(j, i), i) by 1;
endfor
F (p − 1, s); /* p > 0 by Line 11 of 9.1 */

(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
(24)
(25)
(26)
(27)
(28)
(29)
(30)
(31)
(32)
(33)
(34)

if α.∆ == ILLEGAL and α.S == ILLEGAL
set ν.∆ = ILLEGAL;
endif
Reinsert j in the doubly linked lists DLLi ;
Increase by 1 the various BSZ(Q(j, i), i) (see Line 11 above);
CurrentNode ← &ν;
endif
if ν.S 6= ILLEGAL
Create node β;
Treelocation ← &β;
ν.S ← &β;
β.parent ← &ν;
foreach i such that j ∈ DLLi
Delete j from DLLi ;
Decrease BSZ(Q(j, i), i) by 1;
endfor
if f is even
k ← smallest element of DLLf −1 ;
if rk,f −1 == rj,f
β.(dim, vol) = ν.(dim, vol);
else /* By Line 31 of Subsection 9.1 rk,f −1 < rj,f . */

(35)
(36)
(37)
(38)
(39)
(40)
(41)
(42)
(43)
(44)
(45)
(46)

β.vol ← ν.vol × (rj,f − rk,f −1 );
β.dim ← ν.dim + 1;
endif
endif
F (p, s − 1); /* s > 0 by Line 5 of Subsection 9.1 */
if β.∆ == ILLEGAL and β.S == ILLEGAL
set ν.S = ILLEGAL;
endif
Reinsert j in the doubly linked lists DLLi ;
Increase by 1 the various BSZ(Q(j, i), i) (see Line 11 above);
CurrentNode ← &ν;
endif

16

10

Estimate for number of steps

By induction on p + s, there is a unique function f (p, s) with the following properties:
f (p, 0) = 1 for p ≥ 0 ;
f (0, s) = s + 1 for s ≥ 0 ;
f (p, s) = f (p − 1, s) + f (p, s − 1) + d for p > 0 and s > 0 .

(14)

From the pseudocode in Section 9, one sees that there is a constant k > 0, independent
of N , d and r, such that the number of steps necessary to execute F (p, s) is bounded by
k.f (p, s).
Using the addition formula for binomial coefficients
  
 

n
n−1
n−1
=
+
,
s
s
s−1
one checks that f (p, s) is given by

 

p+s
p+s+1
f (p, s) = d
+
− d,
s
s
since it satisfies Equation 14. In particular, we are interested in

 

r + 2d
r + 2d + 1
f (r, 2d) = d
+
− d.
2d
2d

(15)

(16)

To complete the computation of the upper bound for the number of steps involved
in our algorithm, we need to extract the information in the data structure Results (see
8.8). From ρ = Results(p), we find the set of discarded d-rectangles by following parent
pointers, starting at ρ.address. This takes at most O(p + 2d) steps. Summing over
0 ≤ p ≤ r, we obtain a bound
O((r + 1)2d + r(r + 1)/2).

(17)

steps, where the constant implicit in the O notation is independent of N , d and r.
Combining this with the bounds of Equations 13 and 16, we get an overall bound
 


r + 2d
O d
+ dN
(18)
2d
since the omitted terms are dominated by those that remain in Equation 18.

11

The dual situation

Let S be the set of closed non-empty d-rectangles in Rd , together with the null-set. We
have a lattice: the greatest lower bound of a finite subset of S is given by intersection.
The least upper bound of any finite subset of S exists, but is not in general equal to the
union.
Theorem 2.5 applies equally well to the dual situation. Here we are given N drectangles and a number r with 0 ≤ r ≤ N . The task is to find which r rectangles to
discard, so that the smallest rectangle containing all except these r rectangles is as small
as possible.
We can then prove the dual result to Theorem 2.5 by reversing all inequalities and
replacing intersection by least upper bound. It is possible to sketch scenarios in which
such a result might be used, but we spare the reader rather than labour the point.
17

References
[1] M. Blum, R.W. Floyd, V. Pratt, R.L. Rivest, and R.E. Tarjan. Time bounds for
selection. Journal of Computer and System Sciences, 7(4):448–461, 1973.
[2] C.A.R. Hoare. Algorithm 63: partition algorithm 65: Find. Communications of the
ACM, 4(7):321–322, 1961.
[3] Shan-e-Ahmed E. Raza, Ahmad Humayun, Sylvie Abouna, Tim W. Nattkemper,
David B.A. Epstein, Michael Khan, and Nasir M. Rajpoot. Ramtab: Robust alignment of multi-tag bioimages. Plos One, 7(2):e30894–e30894, 2012.
[4] Walter Schubert, Bernd Bonnekoh, Ansgar J. Pommer, Lars Philipsen, Raik Bckelmann, Yanina Malykh, Harald Gollnick, Manuela Friedenberger, Marcus Bode, and
Andreas W. M. Dress. Analyzing proteome topology and function by automated multidimensional fluorescence microscopy. Nature biotechnology, 24(10):1270–8, 2006.

18

