The Algorithm.

Given the graph representation of the experimental sample, G=(V,E), and a prespecified threshold k, the approximate blocking algorithm proceeds as follows:

• 1. Construct a (k−1)-nearest neighbor subgraph of G. Denote this graph Gnn=(V,Enn).

• 2. Find a maximal independent set of vertices, S, in the second power of the (k−1)-nearest neighbor subgraph, Gnn2. Vertices in S are referred to as the block seeds.

• 3. For each seed i∈S, create a block comprising its closed neighborhood in Gnn, Vi=NGnn[i].

• 4. For each yet unassigned vertex, assign it to any block that contains one of its adjacent vertices in Gnn.

When the algorithm terminates, the collection of blocks, balg={Vi}i∈S, is a valid threshold blocking of the experimental units that satisfies the optimality bound.

Informally, the algorithm constructs the blocking by selecting suitable vertices, the seeds, from which the blocks are grown. Seeds are spaced sufficiently far apart so as not to interfere with each other’s growth, but they are dense enough so that all nonseed vertices have a seed nearby. Specifically, the second step of the algorithm prevents any vertex from being adjacent to two distinct seeds in the (k−1)-nearest neighbor subgraph, but also never more than a walk of two edges away from a seed. This ensures that the seeds’ closed neighborhoods do not overlap, and that vertices assigned to the same block are at a close geodesic distance. Fig. 1 illustrates how the algorithm constructs blocks in an example sample.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

An illustration of the approximation algorithm for a sample with 2D covariate data when a minimum block size of two is desired (k=2). (A) The algorithm is provided with a set of data points and forms the graph by drawing an edge between all possible pairs of units. The edges are here omitted to ease presentation. (B) A (k−1)-nearest neighbor subgraph is constructed. (C) The second power of the nearest neighbor subgraph is derived, as shown by the edges, and a maximal independent set is found, as shown by the red vertices (the seeds). (D) All vertices adjacent to a seed in the nearest neighbor subgraph are included in the blocks formed by the seeds, as shown by the edges marked in red. (E) The two yet unassigned vertices are assigned to the blocks that contain one of their adjacent vertices in the nearest neighbor subgraph. (F) The final blocking.

Validity and Complexity.

We first prove that the algorithm is guaranteed to produce a valid threshold blocking, and then examine its time and space complexity.

Lemma 1. For any nonseed vertex, i∉S:

• 1. There exist no two seeds both adjacent to i in Gnn.

• 2. There exists a walk in Gnn of two or fewer edges from i to the seed of the block that i is assigned to.

Proof: The lemma follows from S being a maximal independent set in Gnn2. Refer to Proof of Lemma 1 for a complete proof.

Theorem 1. (Validity) The blocking algorithm produces a threshold blocking: balg∈Bk.

Proof: By Lemma 1, each vertex assigned in the third step is adjacent to exactly one seed; thus it will be in exactly one block. In the fourth step, vertices are assigned to exactly one block each. This ensures that the blocks are disjoint and span V; thus balg is a partition of V.

All seeds have at least k−1 adjacent vertices in Gnn. In the third step, these vertices and the seeds themselves will form the blocks, ensuring that each block contains at least k vertices. This satisfies Definition 1.

Theorem 2. (Complexity) The blocking algorithm terminates in polynomial time using O(kn) space.

Proof: Naively, the (k−1)-nearest neighbor subgraph can be constructed by sorting each vertex’s edge costs and finding its k−1 nearest neighbors. Thus, Gnn is constructed in, at most, O(n2⁡log⁡n) time (18). To enable constant time access to the neighbors of any vertex, store the nearest neighbor subgraph in n lists containing each vertex’s edges. There can be, at most, (k−1)n edges in Gnn. This implies an O(kn) space complexity for the edge lists.

Using the edge lists, a maximal independent set in the second power of Gnn can be found in O(kn) time without changing the space complexity. See Subroutine for the Second Step of the Algorithm for additional details. The third step is completed within O(n) time as Lemma 1 ensures that, at most, n units will be assigned to blocks in this step and the edge lists enable constant time access to the seeds’ neighbors. In the fourth step, it will never be necessary to search through all edge lists more than once, implying a complexity of O(kn).

Remark 1: After the initial construction of the (k−1)-nearest neighbor subgraph, the algorithm terminates in O(kn) time. As the nearest neighbor search problem is well studied, the naive subroutine in the proof can be improved on in most applications. In particular, most experiments will have reasonably low-dimensional metric spaces. Using specialized algorithms, the subgraph can, in that case, be constructed in O(kn⁡log⁡n) expected time (19) or worst-case time (20). If the covariates are not few to begin with, it is often advisable to use some dimensionality reduction technique before blocking so as to extract the most relevant information.

Run time can also be improved by using an approximate nearest neighbor search algorithm. However, approximate optimality is not guaranteed in that case.

Remark 2: It is rarely motivated to increase the block size as the sample grows; thus k can be considered fixed in the typical experiment. When k is fixed and one can use a specialized procedure to derive the nearest neighbor subgraph, the algorithm has O(n⁡log⁡n) time and O(n) space complexity.

Approximate Optimality.

To prove the optimality bound, we will first show that the edge costs in the (k−1)-nearest neighbor subgraph are bounded. As the algorithm ensures that vertices in the same block are at a close geodesic distance in that subgraph, approximate optimality follows from the triangle inequality.

Lemma 2. No edge cost in Gnn can be greater than the maximum cost in the optimal blocking,∀ij∈Enn,cij≤λ.[6]Proof: Consider the graphGλ=(V,Eλ={ij∈E:cij≤λ}).[7]For all edges in an optimal blocking, ij∈E(b∗), we have cij≤λ from optimality. It follows that E(b∗)⊆Eλ.

Let c+=max{cij:ij∈Enn} and considerG+=(V,E+={ij∈E:cij<c+}).[8]The minimum degree of this graph, δ(G+), must be less than k−1. If not, a (k−1)-nearest neighbor graph exists as a subgraph of G+. As this new graph does not contain c+, it is contradictory that c+ is the maximum edge cost in Gnn.

Suppose that c+>λ. It then follows that Eλ⊆E+; thusδ[G(b∗)]≤δ(Gλ)≤δ(G+)<k−1.[9]That is, there exists a vertex in G(b∗) with fewer than k−1 edges. It follows that there must exist a block in G(b∗) with fewer than k vertices and, as Definition 1 then is violated, it cannot be a valid blocking. The contradiction proves that c+≤λ which bounds all edges in Enn.

Theorem 3. (Approximate optimality) The blocking algorithm is a 4-approximation algorithm,maxij∈E(balg)cij≤4λ.[10]Proof: Let balg denote the blocking produced by the algorithm. Consider any within-block edge ij∈E(balg). We must show that cij is bounded by 4λ.

If ij∈Enn, we have cij≤λ by Lemma 2. If ij∉Enn and i∉S, j∈S, then, by Lemma 1, there exists some ℓ so that iℓ,ℓj∈Enn. Lemma 2 applies to both these edges. By Eq. 1, the triangle inequality, it follows,cij≤ciℓ+cℓj≤λ+λ=2λ.[11]If ij∉Enn and i,j∉S, let ℓ∈S be the seed in the block that vertices i and j are assigned to. From above, we have ciℓ,cℓj≤2λ, and, by the triangle inequality,cij≤ciℓ+cℓj≤2λ+2λ=4λ.[12]As there is exactly one seed in each block, i,j∈S is not possible, and we have considered all edges in E(balg).

Remark 3: In some settings, a slight reduction in the sample size is acceptable or required, e.g., for financial constraints or when blocks are constructed before units are sampled using secondary data sources. In these cases, the algorithm can easily be altered into a 2-approximation algorithm. By terminating at the end of the third step and disregarding the unassigned vertices, one ensures that all remaining vertices are, at most, a distance of λ from the seed (where λ refers to the maximum distance in the optimal blocking of the selected subsample). Applying the triangle inequality proves that all edge costs in the blocking of the subsample are bounded by 2λ. It is also possible to apply a caliper to the blocking so as to restrict the maximum possible edge cost by excluding some hard-to-block vertices.

A concern when using a bottleneck objective is that densely populated regions of the sample space will be ignored, as the blocks in these regions will not affect the maximum edge cost. This is especially worrisome when there are a few hard-to-block vertices that result in a large λ. This can lead to poor performance, as covariate balance often can be improved by ensuring good block assignments for all vertices. However, as the presented algorithm does not directly use the bottleneck objective to form the blocks, it avoids this issue. Instead, its optimality follows from the use of the nearest neighbor subgraph as the basis of blocking, and from this graph’s connection with the optimal edge cost as shown in Lemma 2.

The following theorem shows that our algorithm leads to approximate optimality not only in the complete sample but also in all subsamples. Thus, if there is a densely populated region, the algorithm ensures that the blocking is near optimal also within that region.

Theorem 4. (Local approximate optimality) Let bsub⊆balg be any subset of blocks from a blocking constructed by the algorithm. Define Vsub=∪Vx∈bsubVx as the set of all vertices contained in the blocks of bsub. Let λsub denote the maximum edge cost in an optimal blocking of Vsub. The subset of blocks is an approximately optimal blocking of Vsub,maxij∈E(bsub)cij≤4λsub.[13]Proof: Theorem 4 is proven in Proof of Theorem 4.

Heuristic Improvements.

The algorithm allows for several improvements of heuristic character. Although the guaranteed optimality bound or complexity level remains unchanged, we expect these changes to improve general performance. In particular, the algorithm has a tendency to construct blocks that are too large. Although flexibility in the block size is beneficial—the main idea behind threshold blocking—the current version tends to overuse that liberty.

The first improvement exploits an asymmetry in the nearest neighbor subgraph that currently is disregarded. The cardinality condition is met as each seed and its k−1 nearest neighbors are assigned to the same block. However, in addition to those necessary neighbors, the current version assigns vertices that have the seed as their nearest neighbor to the block. With some minor alterations, detailed in Algorithm Using Nearest Neighbor Digraphs, the algorithm can use a (k−1)-nearest neighbor digraph to form the blocks. This digraph is such that an arc (i.e., directed edge) is drawn from i to j if j is among the (k−1) closest neighbors of i. Using the digraph, one can differentiate whether an edge indicates that a vertex is a neighbor of the seed or vice versa. Fig. S1 illustrates the difference between the undirected and directed versions of the algorithm.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S1.

An illustration of the directed version of the blocking algorithm with a comparison with the undirected version. (A) The directed algorithm creates a nearest neighbor digraph by drawing an arc from each vertex to its closest neighbor. (A′) The undirected algorithm draws the same graph but disregards the directionality of the edges. (B) The directed version finds seeds (red vertices) so that no two seeds point toward the same vertex. (B′) As the undirected version disregards the direction of the edges, it requires that the shortest path between any two seeds in the nearest neighbor graph is at least three edges regardless of whether that is needed to satisfy the block-size threshold. It, thereby, misses one possible seed in this example. (C) Blocks are formed with the seeds’ closest neighbors. (C′) Blocks are formed both with the seeds’ closest neighbors and with vertices that have the seeds as their closest neighbors. (D and D′) Remaining vertices are assigned to the block containing their closest neighbor.

There is rarely a unique maximal independent set in the second power of the nearest neighbor graph (i.e., the seeds). The current version selects one arbitrarily. The second improvement is to choose the seeds more deliberately. As each seed assigns at least k−1 vertices to its block, a straightforward way to reduce the block sizes is to maximize the number of seeds—a larger set is expected to produce better blockings. The ideal may be the maximum independent set, but deriving such a set is an NP-hard problem. Most heuristic algorithms are, however, expected to perform well.

Third, despite the above improvements, the algorithm will occasionally produce blocks that are much larger than k. Whenever a block contains 2k or more vertices, it can safely be split into two or more blocks, each containing at least k vertices. As the algorithm ensures that all edge costs satisfy the optimality bound and no edges are added by splitting, this can only lower the maximum within-block cost. In Greedy Threshold Algorithm, we describe a greedy threshold blocking algorithm for splitting blocks that runs fast and is expected to perform well. This algorithm can also be used to block the complete sample, but will not perform on par with the approximation algorithm.

A fourth improvement changes how vertices are assigned in the fourth step. With larger desired block sizes, some blocks may contain peripheral vertices that are far from their seeds. In these cases, it is often beneficial to assign the remaining vertices in the fourth step to the block containing their closest seed instead of the block containing a closest neighbor. This avoids the situation where a vertex is assigned to distant block due to having a peripheral vertex close by. The optimality bound is maintained as Theorem 3 ensures that at least one seed exists at a distance of at most 2λ. If a vertex is not assigned to that seed, it must have been assigned to a seed that is at a closer distance.

Finally, once a blocking is derived, searching for moves or swaps of vertices between blocks can lead to improvements as in other partitioning problems (21). It is, however, seldom feasible to let such searches continue until no additional improvements are possible (i.e., until a local optimal is found), as the flexible block structure allows for a vast number of moves and swaps.

Run Time and Memory.

To investigate the resource requirements, we let each algorithm block samples containing between 100 and 100 million data points generated from the model above. Each setting was replicated 250 times. As time and memory use are quite stable over the runs, these replications suffice to investigate even small differences. In these simulations, the directed version performs almost identically to the original version, and it is omitted from the graphs to ease presentation.

Fig. 2 presents the time and memory needed for the algorithms to successfully terminate. All versions of the approximation algorithm run fast and terminate within a minute for samples sizes up to a few millions. With 100 million data points, the original version terminates within 11 min, whereas the version will all three refinements does so in less than 16 min—all very manageable in real applications. Memory use is increasing linearly at a slow rate for both versions. A modern desktop computer would have enough memory to block samples with tens of millions of units.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Run time (A and B) and memory use (C and D) of five blocking algorithms with 2D input data over a range of sample sizes. Marker symbols are actual simulation results, and the connecting lines are interpolations. Results are presented with different scales due to the large differences in performance. Results are presented for all algorithms for sample sizes up to 40,000 data points (A and C), whereas results for sample sizes up to 100 million data points are only shown for the two approximation algorithms (B and D). No simulations were successful for the greedy algorithms and nonbipartite matching for sample sizes larger than 20,000, due to excessive run time or memory use. The undirected version of the approximation algorithm (shown in red) has almost identical run time and memory use as the directed version, as described in Heuristic Improvements, and its results are not shown in the figure. See Table S4 for detailed results.

The three comparison algorithms paint another picture altogether. For the rather modest sample size of 20,000 data points, these algorithms take more than 20 min—up to 2 h—to terminate. Even more problematic is their extensive memory use. For samples larger than 50,000 data points, all three algorithms try to allocate more than 48 gigabytes of memory; under these settings, these algorithms do not terminate successfully, and no results can be shown.

Detailed results for these simulations and simulations with input data with higher dimensionality are presented in Tables S4–S6.

Minimizing Distances.

To investigate how well the algorithms minimize distances, we increase the number of replications to 5,000; this statistic is less stable, and the differences between algorithms are smaller. However, with this number of replications, no difference can be attributed to simulation error.

The first three data columns of Table 1 show the maximum within-block distance, averaged over the simulation rounds, when the desired minimum block size is two. Refer to Table S1 for results when the desired minimum block size is four. We report values normalized by the performance of the approximation algorithm to ease interpretation. The two improved versions of the approximation algorithm lead to quite substantial decreases in the maximum distance. All versions of the approximation algorithm outperform the two greedy algorithms—for the fixed-sized version, drastically so. However, only the version with all three improvements performs better than nonbipartite matching.

In the fourth through sixth data columns of Table 1, the average within-block distance is presented. This measure is the objective of the greedy fixed-sized algorithm and nonbipartite matching, so it is not surprising that these algorithms show better performance. Nonbipartite matching is, here, the best performing algorithm, and only the approximation algorithm with all three improvements outperforms the fixed-sized greedy algorithm.

The seventh through ninth data columns in Table 1 present the average size of the blocks produced by the different algorithms. The improvements discussed in Heuristic Improvements are shown to be effective in taming the original algorithm’s tendency to construct blocks that are too large. However, smaller blocks do not automatically lead to better performance. This is evident from the greedy threshold algorithm, which produces smaller blocks than the approximation algorithms but has worse performance. The two fixed-sized blocking algorithms produce blocks of constant size by construction.

Reducing Uncertainty.

To investigate how the blockings affect an estimator’s performance, one must specify a data-generating process for the outcome. The results are highly sensitive to the details of this process. Even blockings that are optimal with respect to within-block distances need not lead to the lowest variance. An extensive investigation is beyond the scope of this paper, and we provide only indicative results.

We consider a simple setting with two treatment conditions when the desired minimum block size is two. Refer to Table S2 for results when the minimum block size is four. The outcome (y) is the product of the covariates with additive, normally distributed noise,y=x1x2+ε,ε≈N(0,1).[15]Note that the treatment does not enter into the model and thus has no effect—the potential outcomes are equal. The covariates are highly predictive in this model—more so than one can expect in a real application. This enables the methods to make the most out of the covariate information and helps us differentiate between the algorithms’ performances. For all blocking methods, we will use a difference-in-means estimator that is weighted by the size of each block to estimate treatment effects. Refer to Estimation Methods or ref. 23 for additional details on estimation with threshold blocking. We ran 5,000 replications in this setting, and results are presented relative to the approximation algorithm. Tables S3 and S7 present the results without normalization.

Table 2 presents results for the root of the average squared difference between estimates and the true treatment effect (RMSE) for each method’s estimator. All blocking methods with proven optimality level perform well. For the smaller sample sizes, the approximation algorithm with all three improvements performs best, and nonbipartite matching is slightly better in the larger samples. All blocking methods seem, however, to converge as the sample size grows.

In addition to the six blocking methods, Table 2 includes the results of two methods that do not use blocking. In both cases, the RMSE is markedly higher than in any of the methods using blocking. When controlling for imbalances using ordinary least square regression (OLS), the RMSE is at least twice as large as that of any blocking method with proven optimality. When the estimate is completely unadjusted for imbalances, the RMSE is up to 20 times higher than for the approximation algorithm. However, this difference certainly overstates the benefits of blocking that one can expect in real applications, as the covariates are unusually predictive in this simulation.