Selection on X₁ + X₂ + ⋯ + X_m via Cartesian product trees

Layer-ordered heaps

In a standard binary heap, the only known relationships between a parent and a child is A_i ≤ A_children(i). A layer-ordered heap (LOH) has stricter ordering than the standard binary heap, but is able to be created in $\Theta \left(nlog\left(\frac{n}{n\cdot \left(\alpha -1\right)+1}\right)+\frac{n\cdot \alpha \cdot log\left(\alpha \right)}{\alpha -1}\right)=\Theta \left(n\right)$ time for constant α > 1 (Pennington et al., 2020). α is the rank of the LOH and determines how fast the layers grow. A LOH partitions the array into several layers, L_i, which grow exponentially such that |L₁| = 1 and $\frac{\left|L_{i+1}\right|}{\left|L_i\right|}\approx \alpha$. Every value in a layer L_i is ≤ every value in proceeding layers L_i+1, L_i+2, … which we denote as L_i ≤ L_i+1, this is the “layer-ordering property.” If α = 1,  then all layers are size one and the LOH is sorted; therefore, to be constructed in Θ(n) time the LOH must have α > 1.

Pairwise selection

Serang’s method of selection on A + B utilizes LOHs to be both optimal in theory and fast in practice. In this section we provide a summary of Serang’s method, for a more detailed description with analysis see Serang (2020). The method has four phases. Phase 0 is simply to LOHify (make into a layer-ordered heap) the input arrays which can be done in Θ(n) time.

Phase 1 finds which layer products may be necessary for the selection. A layer product, A^(u) + B^(v), is the Cartesian product of layers A^(u) and B^(v): $A_1^{\left(u\right)}+B_1^{\left(v\right)},A_2^{\left(u\right)}+B_1^{\left(v\right)},\dots ,A_1^{\left(u\right)}+B_2^{\left(v\right)},\dots$. Finding which layer products are necessary for the selection can be done using a standard binary heap. A layer product is represented in the binary heap in two separate ways: a min tuple ⌊u, v⌋ = (min(A^(u) + B^(v)), (u, v), false) and a max tuple ⌈(u, v)⌉ = (max(A^(u) + B^(v)), (u, v), true). Creating the tuples does not require calculating the Cartesian product of A^(u) + B^(v) since min(A^(u) + B^(v)) = min(A^(u)) + min(B^(v)) which can be found in a linear pass of A and B separately. The same argument applies for ⌈(u, v)⌉. false and true note that the tuple contains the minimum or maximum value in the layer product, respectively. Also, let false = 0 and true = 1 so that a min tuple is popped before a max tuple even if they contain the same value.

Phase 1 uses a binary heap to retrieve the tuples in sorted order. When a min tuple, ⌊(u, v)⌋, is popped, then its neighbors (⌊(u + 1, v)⌋, ⌊(u, v + 1)⌋) and the corresponding max tuple are pushed, assuming ⌊(u + 1, v)⌋, ⌊(u, v + 1)⌋ are in bounds (a set is used to ensure a layer product is inserted only once). When a max tuple is popped, a variable s is increased by |A^(u) + B^(v)| = |A^(u)|⋅|B^(v)| and (u, v) is appended to a list q. This continues until s ≥ k. Figure 1 shows an example of phase 1 when k = 10.

In phases 2 and 3 all max tuples still in the heap have their index appended to q, then the Cartesian product of all layer products in q are generated. A linear time one-dimensional k-select is performed on the values in the Cartesian products to produce only the top k values in A + B. The time complexity of the algorithm is linear in the overall number of values produced which is O(k).

In this paper we efficiently perform selection on X₁ + X₂ + ⋯ + X_m by utilizing Serang’s pairwise selection method.

Tree construction

The tree has height ⌈log₂(m)⌉ with m leaves, each one is a wrapper around one of the input arrays which are unsorted and have no restrictions on the type of data they contain. Upon construction, the input arrays are LOHified in Θ(n⋅m) time, which is amortized out by the cost of loading the data. Each node in the tree above the leaves performs pairwise selection on two LOHs: one generated by its left child and one generated by its right child. All nodes in the tree generate their own LOH, but this is done differently for the leaves vs the pairwise selection nodes. When a leaf generates a new layer it simply allows its parent to have access to the values in the next layer of the LOHified input array. For a pairwise selection node, generating a new layer is more involved.

Pairwise selection nodes

Each node above the leaves is a pairwise selection node, each of which has two children that may be leaves or other pairwise selection nodes. In contrast to the leaves, the pairwise selection nodes will have to calculate all values in their LOHs by generating an entire layer at a time. Generating a new layer requires performing selection on A + B, where A is the LOH of its left child and B is the LOH of its right child. Due to the combinatorial nature of this problem, simply asking a child to generate their entire LOH can be exponential in the worst case so they must be generated one layer a time and only as necessary.

The pairwise selection performed is Serang’s method with a few modifications (though Fig. 1 is still representative of the method). The size of the selection is always the size of the next layer, k = |L_i|, to be generated by the parent. The selection begins in the same way as Serang’s: a heap is used to pop min and max layer product tuples. When a min tuple, ⌊(u, v)⌋, is popped the values in the Cartesian product are generated and appended to a list of values to be considered in the k-selection. The neighboring layer products inserted into the heap are determined using the scheme from Kaplan et al. which differs from Serang’s method: ⌈(u, v)⌉ and ⌊(u, v + 1)⌋ are always inserted and, if v = 1, ⌊(u + 1, v)⌋ is inserted as well. This insertion scheme will not repeat any indices and therefore does not require the use of a set to keep track of the indices in the heap. When any min tuple is proposed, the parent asks both children to generate the layer if it is not already available. If one or both children are not able to generate the layer (i.e., the index is larger than the full Cartesian product of the child’s children) then the parent does not insert the tuple into its heap. The newly generated layer is simply appended to the parent’s LOH and may now be accessed by the parent’s parent. An example of a pairwise selection node generating a new layer which requires a child to generate a new layer is shown in Fig. 2.

The dynamically generated layers should be kept in individual arrays, then a list of pointers to the arrays may be stored. This avoids resizing a single array every time a new layer is generated.

Theorem 1 in Serang (2020) proves that Serang’s method performs selection in O(k) time. lemmas 6 and 7 show that the number of all values generated is O(n + k); however, lemma 7 may be amended to show that any layer product of the form (u, 1) or (1, v) will generate ≤α⋅|(u − 1, 1)| ∈ O(k) or ≤α⋅|(1, v − 1)| ∈ O(k) values, respectively, to show that the total values generated is O(k). Thus, the total number of values generated when a parent adds a new layer L_i is O(|L_i|).

Selection from the root

In order to select the top k values from X₁ + X₂ + ⋯ + X_m, the root is continuously asked to generate new layers until the cumulative size of the layers in their LOH exceeds k. Then a k-selection is performed on the layers to retrieve only the top k. Due to the layer-ordering property, a selection only needs to be performed on the last layer, all values in previous layers will be in the top k.

The Cartesian product tree is constructed in the same way as the FastSoftTree (Kreitzberg, Lucke & Serang, 2020a) as both methods dynamically generate new layers in a similar manner with the same theoretical runtime. In both methods, the pairwise selection creates at most O(α²⋅k) values. Thus the theoretical runtime of both methods is O(n⋅m + k⋅m^log₂(α2)) with space usage O(n⋅m + klog(m)).

Wobbly version

In Serang’s pairwise selection, after enough layer product tuples are popped from the heap to ensure they contain the top k values, there is normally a selection performed to remove any excess values. Strictly speaking, this selection is not necessary anywhere in the tree except for the root when the final k values are returned. When the last max tuple ⌈(u, v)⌉ is popped from the heap, max(A^(u) + B^(v)) is an upper bound on the k^th value in the k-selection. Instead of doing a k-selection and returning the new layer, which requires a linear time selection followed by a linear time partition, we can simply do a value partition on max(A^(u) + B^(v)).

A new layer generated from only a value partition and not a selection is not guaranteed to be size k, it is at least size k but contains all values ≤max(A^(u) + B^(v)). In the worst case, this may cause layer sizes to grow irregularly with a constant larger than α. For example, if k = 2 and |L₁| = |L₂| = 1 then in the worst case every parent will ask their children to each generate two layers and the value partition will not remove any values. Each leaf will generate two values, their parents will then have a new layer of size 2² = 4, their parent will have a new layer of size (2²)² = 16, etc. Thus the root will have to perform a 2-selection on 2^m values which will be quite costly.

In an application like calculating the most probable isotopologues of a compound, this version can be quite beneficial. For example, to generate a significant amount of the isotopologues of the titin protein may require k to be hundreds of millions. Titin is made of only carbon, hydrogen, nitrogen, oxygen, and sulfur so it will only have five leaves and a tree height of three. The super-exponential growth of the layers for a tree with height three is now preferential because it will still not create so many more than k values but it will do so in many fewer layers with only value partitions and not the more costly linear selections. We call this the “wobbly” Cartesian product tree.

Cartesian product trees vs FastSoftTree

In a Cartesian product tree, replacing the pairwise A + B selection steps from Kaplan et al.’s soft heap-based algorithm with Serang’s optimal LOH-based method provides the same o(n⋅m + k⋅m) theoretical performance for the Cartesian product tree but is practically much faster (Table 1). This is particularly true when k⋅m^log₂(α2) ≫ n⋅m, where popping values dominates the cost of loading the data. When k ≥ 2¹⁰, k⋅m^0.2750 < n⋅m which is reflected in our results where for k = 2²⁰ we get a 630.4 × speedup, significantly larger than for k = 2¹⁰ which only has a 10.27 × speedup. For any reasonable ϵ and α, the difference in number of values generated (and thus in memory usage) will not differ significantly. Regardless of which algorithm generates less values, the performance gain due to generating less overall values will be negligible compared to the difference in cache performance between soft heaps and LOHs.

Standard vs Wobbly

As we see in Table 2, for small m the Cartesian product tree can gain significant increases in performance when there are no linear selections performed in the tree and the layers are allowed to grow super-exponentially. As k grows, the speedup of the wobbly version continues to grow, resulting in a 1.786 × speedup for k = 2³⁰. When m ≫ 5 the growth of the layers near the root start to significantly hurt the performance. For example, if n = 32, m = 256 and k = 256 the wobbly version takes 0.5805 seconds and produces 149,272 values at the root compared to the non-wobbly version which takes 1.8810 × 10⁻⁰³ seconds and produces just 272 values at the root.