Imagine \(n\) objects placed in a multidimensional space, stored across a distributed network of \(p\) compute nodes. The problem is simple:

  • how to divide the universe in such a way that each compute node gets an equal number of objects \(n/p\)? and
  • how to ensure that the objects of each compute node are as close as possible to each other? and
  • how to ensure that the regions of overlap of objects across compute nodes are minimal or nonexistent?

This is a very common problem setup, and one of high importance in large-scale scientific applications. The main rationale is that a good spatial decomposition should allow close objects to be stored in the same memory region (compute node), so that inter-node dependency (communication) is reduced. This post compares the simplest spatial decompositions, and details a distributed algorithm for the most efficient of them.

The simplest spatial decomposition algorithms perform an orthogonal spatial decomposition (slicing). In a 2D space, there are several slicing alternatives, numbered 1 to 4 from left to right:

Single Axis
 Non-Uniform Grid
 Orthogonal Recursive Bisection
 Sort Tile Recursive

A small note on efficiency. Slicing layout 1 incurs a lot of object duplication when objects are not points in space but have volume, and when there is a high number of slices. Layout 2 works better by reducing cloning, yet it does not guarantee an accurate spatial decomposition, as the cuts of each dimension are shared by all slices, which may then hold different numbers of objects.

Layout 3 is the orthogonal recursive bisection (ORB), introduced by Berger and Bokhari in 1987 to balance the work of nonuniform problems on multiprocessors. Layout 4 is the sort tile recursive (STR) layout, introduced by Leutenegger et al. in 1997 to pack the R-trees of spatial databases. Layouts 3 and 4 are similar, but layout 4 outperforms layout 3 in parallelism. This is because layout 3 performs a tree-based recursive bisectional decomposition (i.e. at every level of the tree it computes 2 slices), while layout 4 performs several slices at once for each of the 2 dimensions.

The following table details the complexity of each algorithm (assuming sorting with an average complexity of \(n \log n\)):

Scheme Computation steps Properties / Limitations

Single Axis		 \(= n \cdot \log_2 n\) One sorting of all elements on a single direction, independently of the number of dimensions.

Non-Uniform Grid		 \(= d \cdot n \cdot \log_2 n\) One sorting of all elements applied to each dimension, iteratively; slicing configuration for the \(d\) dimensions must decompose final number of slices \(s\) i.e. \(\prod_{k=1}^d S_k = s\).

Sort Tile Recursive		 \(\begin{array}{l} = n \cdot \log_2 n \\ +n \cdot \log_2 \frac{n}{S_1} \\ +n \cdot \log_2 \frac{n}{S_1 S_2} \\ +... \\ +n \cdot \log_2 \frac{n}{S_1 S_2 ... S_{d-1}} \\ = \sum_{i=1}^d n \cdot \log_2 \frac{n}{\prod_{k=0}^{i-1} S_k} \end{array}\) One sorting for each subslice’s data on each dimension, applied iteratively; sub-slices computations are independent; slicing configuration for the \(d\) dimensions must decompose final number of slices \(s\) i.e. \(\prod_{k=1}^d S_k = s\); Faster than NUG as each new dimension requires sorting across smaller datasets; the sum assumes \(S_0 = 1\) i.e. the initial universe is a single slice including all elements.

Orthogonal Recursive Bisection		 \(\begin{array}{l} = n \cdot \log_2 n \\ +n \cdot \log_2 \frac{n}{2} \\ +n \cdot \log_2 \frac{n}{4} \\ +n \cdot \log_2 \frac{n}{8} \\ +... \\ +n \cdot \log_2 \frac{n}{s/2} \\ = \sum_{i=1}^{\log_2 s} n \cdot \log_2 \frac{n}{2^{i-1}} \end{array}\) Total number of slices \(s\) must be a power of 2; requires, at every iteration, one sorting of each slice, which holds half of the elements of its parent slice, for a total of \(\log_2 s\) steps, therefore a slow method on large networks; sub-slices computations are independent.

Across all possible slicing layouts, (4) is the most efficient: it guarantees low data duplication, allows for accurate slicing, and requires only \(d\) (dimensions) recursive iterations, regardless of the number of compute nodes (unlike (3)). However, its implementation in a distributed memory environment is not trivial. Possible options are:

  • A histogram/binning algorithm that clusters elements in spatial bins and uses the number of elements per bin as the metric for slicing;
  • A method that performs serial slicing on one compute node that holds only a sample of the data, and uses the slicing of the sampled dataset as the slicing of the final dataset.

Distributed Slicing via Sort-Balance-Split

A technique published by the author in this paper, called the Sort-Balance-Split (SBS), implements a parallel multi-dimensional slicing based on layout (4). The algorithm follows three steps:

  • The sorting step performs a distributed sort of all elements on a given dimension (the implementation uses a sample sort);
  • The load balancing step equalizes the number of elements across compute nodes, while respecting the initial sorting, as illustrated in the following picture:

Load balancing of 19 sorted elements on 3 ranks, colored by rank. The first 19 mod 3 = 1 rank holds one more element. In an MPI implementation, step 1 is one MPI_Allgather of the element counts, and step 2 is one MPI_Alltoall of the send counts followed by one MPI_Alltoallv of the elements. Adapted from: An Efficient Parallel Load-balancing Framework for Orthogonal Decomposition of Geometrical Data

  • The splitting step creates smaller networks, called subnetworks, from subsets of the ranks of the main network, as illustrated in the following picture:

The split of 8 ranks into \(S_k = 2\), \(S_k = 3\) and \(S_k = 4\) subnetworks of consecutive ranks, after the sort and the balance, colored by rank. The data stays in place, and each rank gets a new rank in its subnetwork. In an MPI implementation, the split is one collective operation, MPI_Comm_split. For \(S_k = 3\), which does not divide the number of ranks, the subnetworks have 3, 3 and 2 ranks. Adapted from: An Efficient Parallel Load-balancing Framework for Orthogonal Decomposition of Geometrical Data

The three steps explain the name Sort-Balance-Split. The SBS algorithm is performed recursively at every level, where each level operates on the subnetworks created by the split step of the previous level.

This method has two advantages over existing methods: (1) a fixed complexity, as the number of communication steps is fixed for any network and input size; and (2) an accurate spatial decomposition, i.e. not an approximation. However, SBS moves the objects at every level, while the histogram and sampling methods compute an approximate slicing first and then move each object only once: these are preferable when an approximate balance is enough.

A sample application of the algorithm on a 2D universe of 16 shapes distributed across 4 compute nodes (ranks) is presented below:

Adapted from: An Efficient Parallel Load-balancing Framework for Orthogonal Decomposition of Geometrical Data

From left to right:

  1. Initial data layout;
  2. The initial sorting operation sorts all elements on the first (X) dimension, without caring about the other coordinates. Data is load balanced on the X axis, so that each rank holds 4 shapes;
  3. A network split divides the initial problem of 16 shapes across 4 nodes on a single network into 2 problems of 8 shapes, each on a network of 2 nodes;
  4. Each network performs a distributed sort of its objects on the next axis (Y), and balances the workload;
  5. Final data layout.

Here is another illustration of the same problem in a 3D universe, now emphasizing the memory layout across nodes:

Sort-Balance-Split of 18 objects in 3D on \(p = 8\) ranks, with \(S_1 = S_2 = S_3 = 2\), colored by rank. Each column is one object, with its coordinates \(x\), \(y\) and \(z\). Red frames mark the coordinate that steps 1.1, 2.1 and 3.1 sort in each slice, and blue lines the cuts between the subnetworks. In step 2.1, two splitters of the second subnetwork are equal, so rank 5 receives no object until the balance.

And here is an application of distributed slicing to the N-body problem (left), which allows one to simulate cosmic structure formation after the Big Bang and whose simulations have won the ACM Gordon Bell Prize several times. The picture on the right shows the slicing of a small network of neurons:

Limitations, benchmarks and a more thorough analysis of the problem are available in the original publication. The C++ implementation, with MPI, is available in slicing.tar.gz.

To find more about the topic, the Multi-Jagged algorithm of Deveci et al. computes a layout similar to layout 4, with several cuts per dimension, and runs faster and scales better than the recursive bisection of layout 3.