Adaptive Geospatial Joins for Modern Hardware

02/26/2018 ∙ by Andreas Kipf, et al. ∙ Centrum Wiskunde & Informatica 0

Geospatial joins are a core building block of connected mobility applications. An especially challenging problem are joins between streaming points and static polygons. Since points are not known beforehand, they cannot be indexed. Nevertheless, points need to be mapped to polygons with low latencies to enable real-time feedback. We present an adaptive geospatial join that uses true hit filtering to avoid expensive geometric computations in most cases. Our technique uses a quadtree-based hierarchical grid to approximate polygons and stores these approximations in a specialized radix tree. We emphasize on an approximate version of our algorithm that guarantees a user-defined precision. The exact version of our algorithm can adapt to the expected point distribution by refining the index. We optimized our implementation for modern hardware architectures with wide SIMD vector processing units, including Intel's brand new Knights Landing. Overall, our approach can perform up to two orders of magnitude faster than existing techniques.



There are no comments yet.


page 4

page 5

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Connected mobility companies need to process vast amounts of location data in near real-time to run their businesses. For instance, Uber needs to join locations of passenger requests with a set of predefined polygons to display available products (e.g., Uber X) and to enable dynamic pricing111 Another example are traffic use cases where the positions of vehicles need to be joined with street segments to enable real-time traffic control. With a future of connected (and possibly self-driving) cars and drones, high-performance geospatial joins will become a key feature for increasingly demanding applications.

While geospatial joins have been studied for decades, many of them optimize for I/O operations, which is an important performance factor for disk-based systems. With the large main-memory capacities of modern hardware, however, it is for the first time possible to maintain fine-grained index structures purely in main memory.

We propose an adaptive geospatial join that addresses workloads with streaming points and static polygons. Our approach is optimized for modern hardware with large high-bandwidth memory and our implementation fully utilizes many-core processors and their Single Instruction Multiple Data (SIMD) units.

Fig. 1: We adaptively compute a fine-grained (true hit) filter that identifies most or even all join partners. We can treat candidate hits as true hits while guaranteeing a user-defined precision. The arrow width reflects the amount of data passing through each of the respective paths.

Figure 1 illustrates our approach. We make use of true hit filtering proposed by [1]. The authors of this work used maximum enclosed rectangles and circles as progressive approximations of single polygons. In contrast, we use an adaptive grid represented by a specialized radix tree, Adaptive Cell Trie (ACT), to not only index single polygons but entire sets of polygons.

There are systems (e.g., Oracle Spatial [2] and Spark Magellan222 that have used the idea of true hit filtering. We take this idea one step further and leverage the large main-memory capacities available in modern hardware to maximize the likelihood of true hits. We increase this likelihood by training the index using historical data points.

In summary, we transform the traditionally compute-intensive problem of geospatial joins into a memory-intensive one. Our approach results in speedups of up to two orders of magnitude compared to state-of-the-art standalone libraries and geospatial database systems.

Our contributions include:

  • An adaptive geospatial join that uses hierarchical grids combined with radix indexing techniques to identify most join partners in the filter phase.

  • An approximate version of this approach that guarantees a user-defined precision.

  • An optimization of this technique for Intel’s latest Xeon Phi, named Knights Landing (KNL).

Ii Background

Unless stated otherwise, we follow the semantics of the —ST_Covers— join predicate supported by PostGIS and Oracle Spatial. —ST_Covers— evaluates whether one geospatial object (e.g., a polygon) covers another (e.g., a point). Points on the edges or vertices of a polygon are considered to be within the polygon. Other predicates, such as —ST_Intersects—, are not explicitly addressed in this paper but can be implemented with similar techniques.

Our geospatial join uses several building blocks:

PIP: A point-in-polygon (PIP) test determines whether a point lies within a polygon. Typically such a test is performed using complex geometric operations, such as the ray tracing algorithm

which involves drawing a line from the query point to a point that is known to be outside of the polygon and counting the number of edges that the line crosses. If the line crosses an odd number of edges, the query point lies within the polygon. The runtime complexity of this algorithm is


being the number of edges. While there are many conceptual optimizations to the PIP test, this operation remains computationally expensive since it processes real numbers (e.g., latitude/longitude coordinates) and thus involves floating point arithmetics. In geospatial joins, this test should thus be avoided whenever possible.

Space-filling curves: Geographical coordinates can either be indexed using multi-dimensional data structures, such as R-trees, or they can be linearized and indexed using one-dimensional data structures, such as B-trees or tries. Linearization methods include the Hilbert space-filling curve and the Z curve. Our approach relies on discretization and linearization but does not depend on a concrete space-filling curve. For our approach to work, the cell enumeration must only fulfill the property that children cells share a common prefix with their parent cell.

Google S2: We use the Google S2 library333 that provides primitives that our approach builds upon. These primitives include a PIP test algorithm and adaptive grid approximations of polygons. In general, S2 offers many features for processing geospatial data in both discrete and non-discrete space. To transform coordinates from the non-discrete space (the latitude/longitude coordinate system) into the discrete space, S2 maps points on earth onto a surrounding unit cube. Then it recursively subdivides each of the six faces of the cube in a quadtree fashion and enumerates the cells using the Hilbert space-filling curve. The enumerated cells are 64 bit integers, called cell ids, that uniquely identify a cell. The three most significant bits represent one of the six faces of the cube. The following zero to 30 bit pairs identify the quadtree cell. The next bit is always set and is used to identify the level of a cell. The remaining bits (if any) are set to zero. The smallest cells use all 64 bits.

Figure 2 illustrates the hierarchical decomposition of a cube face and shows how many bits the corresponding cell ids require depending on their level. The colored bits encode the levels. The cell ids in the example share a common prefix with their ancestors up to the colored bit. For example, is an ancestor of and (i.e., fully contains and ). A containment relationship between cells can thus be efficiently computed using bitwise operations.

Fig. 2: Cells at three different levels enumerated by the Hilbert curve and their bit representations

As mentioned above, S2 offers methods for approximating polygons. Figure 2(a) shows a covering of a polygon. A covering is a collection of cells (possibly at different levels) fully covering a polygon. Figure 2(b) shows an interior covering. As the name suggests, an interior covering approximates the interior of a polygon. All cells of an interior covering (interior cells) are completely contained within the polygon. When a point lies in an interior cell, we know that it also lies within the polygon. Both types of coverings can be utilized to speed up PIP tests.

(a) Covering
(b) Interior covering
Fig. 3: A covering and an interior covering of a polygon

To allow for an efficient search, S2 stores the cell ids of a covering in a sorted vector. Besides sorting the cell id vector, it allows for normalizing the covering. A normalized covering neither contains conflicting nor duplicate cells. A conflict between two cells exists, if one cell contains the other. Only when the covering is normalized, cell containment checks can be efficiently () implemented using a binary search on the sorted vector. The normalization of an covering does not lead to a precision loss.

Iii Geospatial Join

With our approach, we address geospatial joins between streaming points and static polygons. An example of such a workload is mapping Uber or DriveNow cars and passenger requests to predefined zones for allocation purposes.

Iii-a Overview

Our geospatial join takes a set of simple polygons, builds an index, and probes points against the index. The high level idea is to index adaptive grid approximations of the polygons in a highly-efficient and specialized radix tree. In contrast to techniques that first reduce the number of polygons using an index, e.g., an R-tree on the polygons’ minimum bounding rectangles (MBRs), and then refine candidates using geometric operations, our approach implements true hit filtering [1] and identifies most or even all join partners in the filter phase. Given a memory budget and a precision bound, we adapt our index to satisfy these constraints while maximizing probe performance.

Our algorithm consists of five phases:

Build logical index: We compute coverings and interior coverings of all polygons and merge them to form a logical index.

Build physical index: Then we build a physical index on the cells of the logical index.

Training: When our approximate approach cannot guarantee a user-defined precision without exceeding a specified memory budget, we use an exact approach and train the index with historical data points to increase probe performance. Popular areas that expect more hits are approximated using a finer-grained grid than less popular areas.

Probe index: In this phase, we probe the points against the physical index. An index probe can either return a false hit or a list of hits, where each hit can either be a true hit or a candidate hit.

Refine candidates: Finally, we perform exact geometric computations for all candidate hits444Often, this phase can be skipped. In fact, it is only required if exact results are needed and if there are candidate hits at all..

We use the S2 library to compute coverings and interior coverings. Note that our technique does not rely on S2. In fact, any other adaptive grid index would work as well. We then merge the individual coverings to create a logical index, before we build a physical index. Note that the first two phases can be performed in one step, i.e., the individual coverings can be computed and merged while inserting their cells into the index. In our approach, the probe (filter) phase is the performance critical part. We therefore highly parallelize this phase using thread- and instruction-level parallelism to accelerate lookups in the physical index. In the refinement phase, we use S2’s PIP test, which implements the ray tracing algorithm (cf., [3] for performance numbers).

The overall strategy of our technique is to minimize the number of (expensive) refinements. We can decrease the number of refinements by using a more fine-grained index. Since we make use of true hit filtering, a precise index allows us to identify most join partners during the filter phase. We can always omit the refinement phase when approximate results are sufficient. The error of such an approximate algorithm is bound by the length of the diagonal of the largest covering cell, which we choose to be sufficiently small.

Our geospatial join takes two parameters: (i) a memory budget that must not be exceeded and (ii) an application-dependent precision555The maximum distance (in meters) between the two partners of a false positive join pair. that must be guaranteed. We first try using the approximate approach by refining the index until we can guarantee the user-defined precision bound. This involves replacing covering cells with their children cells until the largest covering cell is sufficiently small (i.e., its diagonal is less than the precision bound). It might occur that we exhaust the memory budget during this procedure. In that case, we use the exact approach and train the index until we have used the entire memory budget. Training the index reduces the likelihood for expensive geometric operations. We evaluate these two approaches in Section V.

In summary, we trade memory consumption with precision (approximate approach) and probe performance (exact approach). The approximate approach almost involves no computation (except result aggregation) while the exact approach reduces expensive computations by training the index. Thus, both strategies favor modern hardware with large main memory capacities and high memory bandwidths.

Iii-B Build Logical Index

We compute coverings and interior coverings of all polygons, and combine them into a logical index, which we call super covering. The precision of the super covering defines the selectivity of the index. As mentioned earlier, a normalized covering excludes conflicting and duplicate cells. In contrast to the lossless normalization of a single covering, we might experience a precision loss when merging two overlapping coverings and normalizing the resulting super covering666There is no precision loss if the overlap only consists of duplicate cells.. In other words, a normalized super covering may be less selective than two individual coverings. Figures 3(a) and 3(b) show the coverings of two individual polygons. The red cells have conflicts with cells of the other covering. Figure 3(c) shows the normalized super covering. The conflicting cells were expanded to larger cells, which lead to a precision loss.

(a) Covering
(b) Covering
(c) Super cov.
Fig. 4: A normalized super covering may be less selective than two individual coverings. The arrows indicate that the cells will be expanded.

When combining the covering and the interior covering of a single polygon, the situation is even more severe since there will be many overlapping cells. To retain the precision of the individual coverings, we can omit removing conflicting cells and only remove duplicates. However, there is a trade-off between the size of a super covering and its precision. We have experimented with three flavors of a super covering: (i) only covering cells, (ii) covering and interior cells with precision loss, and (iii) covering and interior cells without precision loss. The first approach only indexes covering cells. This approach cannot identify true hits during the filter phase. The second approach additionally indexes interior cells and can thus skip the refinement phase when hitting an interior cell. The first two approaches normalize the super covering when combining the cells leading to a precision loss. The third approach avoids the normalization. Instead of storing a cell and its descendant cell , we compute their difference and store and . This has the advantage that there will not be any overlap between the indexed cells and thus an index lookup will at most return a single cell. The side effect is that the total number of cells will increase since consists of at least three cells. This approach retains the precision and the type (covering or interior) of the individual cells as well as the mappings of cells to polygons777We expect the individual coverings (that serve as inputs to the merge and index phases) to be normalized. Instead, they could also be normalized “on the fly”.. Figure 5 illustrates this precision preserving conflict resolution. Assume that and are cells of two different coverings. First, we compute , which consists of six cells. We then copy all references of to and and omit . Note that the cell count is increased by five.

(a) and
(b) Difference
(c) and
Fig. 5: Precision preserving conflict resolution. is marked in blue, in green, and the cells in in purple. Note that contains .

Each cell in the super covering maps to a list of polygon references. A polygon reference has two attributes:

polygon id: The id of the polygon that this cell references.

interior flag: Whether the cell is an interior or a covering cell of the polygon.

[outerlinewidth=0.1pt,outerlinecolor=black, innerleftmargin=2pt,innerrightmargin=2pt,innertopmargin=-15pt,innerbottommargin=10pt]

  a list of coverings coverings // one per polygon
  a list of interior coverings interiors // one per polygon
  // a list of (cell, polygon references)
  the super covering superCovering
  for (covering in coverings) {
    for (cell in covering) {
      if (superCovering already contains cell) {
        add references of cell to existing cell
      if ( conflicts with existing cell in superCovering) {
        // cell is covered by existing cell or vice versa
        // resolve conflict
         = ascendant cell // may be cell or existing cell
         = descendant cell // may be cell or existing cell
         = difference of  and 
        add references of  to  and 
        remove  from superCovering // only required if the existing cell is the ascendant cell
        add  and  to superCovering
      add {cell, {covering.polygonId, interior flag=false}} to superCovering
  //  same code for interior coverings (with interior flag=true)
Listing 1: Build super covering

Listing 1 outlines the algorithm that builds a super covering of coverings and interior coverings and retains the precision of the individual coverings. We iterate over all input cells and try to insert them into the super covering. When a cell already exists, this means that it is also part of another covering that has already been processed. When a cell conflicts with another cell, this means that either the current cell covers the other cell or vice versa. These two cases may happen when polygons overlap or are close to each other. Conflicts also occur when we first insert the cells of a covering of a given polygon and then the cells of its interior covering. The interior cells always overlap some (if not all) covering cells. By appropriately resolving these conflicts, we retain the precision of the index. As mentioned earlier, this strategy increases the total number of cells. However, a more precise index reduces the number of refinements and thus increases the overall performance of our algorithm. In the remainder of this paper, we use the approach that retains the precision of the individual coverings.

(a) Coarse-grained
(b) Fine-grained
Fig. 6: Super coverings of the five NYC boroughs. Interior cells are marked in green, covering cells of single polygons in blue, and covering cells of multiple polygons in red.

Figure 6 illustrates a coarse-grained and a fine-grained super covering of the five NYC boroughs. For computing the illustrated super coverings, we use up to 50 and up to 200 cells per individual (interior) covering for the coarse-grained and the fine-grained super covering, respectively. The key message of this illustration is that the size of the green area increases with more precise individual coverings, while the size of the red and the blue area decreases. Increasing the size of the green area increases the likelihood that we can skip the refinement phase. Interior cells are not only more likely to be hit than covering cells but they also tend to be larger. Since larger cells use less bits, this leads to an additional performance improvement when probing ACT. The reason for this is that large cells are indexed high up in the tree and are thus found sooner.

Iii-C Build Physical Index

To store the logical index, we use a physical index consisting of two data structures, ACT and a lookup table. Both data structures are designed as in-memory data structures and are optimized for modern hardware. In particular, we utilize large main-memory capacities, avoid unnecessary cache misses, and have parallelized and vectorized the lookups in ACT (cf., Section IV). ACT stores the cell ids of the super covering and the lookup table stores the corresponding lists of polygon references. When a given cell maps to less than three polygon references, we inline the references into the tree nodes avoiding an additional indirection and a possible cache miss.

Adaptive Cell Trie: In S2, the cells of an individual covering are stored in a sorted vector that can be efficiently searched using a binary search (cf., Section II). While this approach has a runtime complexity of , it is not very space efficient, since the cell ids of the cells of a single covering often share a common prefix that is stored redundantly. For a single covering, the space consumption can be reduced by storing the common prefix separately and only storing the varying parts of the cells (e.g., using delta encoding). Since the cells of a single covering are spatially dense and thus likely to share a common prefix, this works well. However, when combining the coverings and interior coverings of multiple polygons, the cells are not necessarily dense anymore. Also, we only require prefix lookups, i.e., we need to search for the cell ids that share a common prefix with our search key (the cell id of a query point) to check whether the query point is contained in one of the indexed cells. A radix tree thus is the ideal data structure for our needs. It is more space-efficient than a sorted vector, since it avoids redundantly storing common prefixes (in a trie, the path to a leaf node implicitly defines the key). Another advantage of a radix tree is that lookups are in with being the key length. Besides the algorithmic complexity, lookup performance depends on the number of cache misses. Since we are generally interested in high selectivity and index many cells, the tree usually exceeds cache size (cf., Section V). Thus, traversing the tree results in many cache misses (at least one per non-cached node).

Let be the average key length and be the fanout of the tree. Then the average costs

of a lookup can be estimated as follows:

= * costs per node access

Increasing the number of indexed cells does not necessarily increase . Thus, the costs of a lookup in a non-cached tree are relatively independent of the number of indexed cells. The number of node accesses is bounded by the maximum key length , which is 60 when using S2’s adaptive grid index888The three face bits are indexed in a dedicated face node and the level encoding bit is not required.. Our tree uses a default fanout () of 256 (= 8 bits). Thus every level in the tree represents four S2 levels (recall that every S2 level is encoded with two bits). This has the side effect that we can only index cells at certain levels.

Let be the S2 level granularity of the tree. Then the following holds for indexed cells:

mod = 0

Thus, we need to denormalize999Denormalizing a cell to a given level means replacing the cell with all of its descendant cells at that level. cells upon insertion and replicate their payloads.

While a fanout of 256 results in sparse nodes and thus in a large space consumption, it allows for efficient lookups. An adaptive radix tree (ART) [4] is usually more space-efficient, however, experiments show that introducing a second node type with four children (Node4 in ART) would (i) only save a negligible amount of space in the case of our workload and (ii) switching between the different node types would have a significant impact on lookup performance. Also, lookups in compressed node types would be more expensive due to the additional indirection for accessing the payloads. Since ACT is designed as a transient data structure that is built at runtime, we favor performance over space consumption. With , the maximum number of node accesses is . In practice, a lower is often sufficient. For example, allows for indexing S2 cells up to level 24 (a level 24 cell represents at most 2  on earth) and limits the number of node accesses to 6.

(a) Adaptive Cell Trie
offset true candidate
0 {5} {3, 1}
1 {7, 2} {8}
(b) Lookup table
Fig. 7: Adaptive Cell Trie and the lookup table. Key parts (bit sequences) are marked in black and values are marked in blue (single or double payloads) and red (offsets).

Figure 6(a) illustrates the structure of ACT. Values (payloads or offsets) can be found in any node at any level of the tree. Every node consists of a fixed-sized array of 256 entries of 8 byte pointers. These pointers are tagged. By default, all entries point to a sentinel node indicating a false hit. Such a tagged entry can be:

  • An 8 byte pointer to a child or the sentinel node

  • An inlined payload (a 31 bit value)

  • Two inlined payloads (two 31 bit values)

  • An offset (a 31 bit value) into a lookup table

We use the two least significant bits of the 8 byte pointer to differentiate between these four possibilities. For an inlined payload, we differentiate between a true hit and a candidate hit using the least significant bit of the 31 bit payload. Thus, we can effectively only store bit payloads (i.e., index up to polygons). We decided to inline two payloads since it is very likely that a cell has less than three payloads in the case of disjoint polygons. We do not inline more than two payloads since this would mean further reducing the supported number of polygons.

Since our implementation uses S2, we need to maintain up to six radix trees, one for each face of the cube (cf., Section II). We call the node that stores the root nodes of these trees face node.

Lookup table: When a cell references more than two polygons, the tree contains an offset into a lookup table. Cells often reference the same set of polygons. To avoid redundancy, we therefore only store unique polygon reference lists. Figure 6(b)

shows an example of a lookup table. The reference lists are split into two parts, a list with true hits and a list with candidate hits. Both lists contain polygon ids. When a cell references at most two polygons, we inline its payloads into the tree. Assuming disjoint polygons and a uniform distribution of points, a more fine-grained super covering may reduce the chance that we need to access the lookup table. Let

be the area that is covered by cells that contain an offset into the lookup table (i.e., cells that reference more than two polygons). Increasing the granularity of the index may reduce the size of () and thus reduces the likelihood of having to access the lookup table. The lookup table is encoded as a single 32 bit unsigned integer array. The offsets stored in the tree are simply offsets into that array. Each encoded entry contains the number of true hits followed by the true hits, the number of candidate hits, and the candidate hits.

Iii-D Training

When we cannot build an index that satisfies a user-defined precision without exceeding the memory budget, we use an approach that may enter the expensive refinement phase. To minimize the likelihood of refinements, we train the index to adapt to the expected distribution of query points.

The index can either be trained online or offline. In this work, we only discuss and evaluate how to train the index in a dedicated training phase. The online case introduces additional concurrency and buffer management issues that we leave for future work.

Training the index minimizes the area that is covered by expensive cells. We define expensive cells as cells that map to polygon reference lists with at least one candidate hit. When we hit such a cell during the join, we need to enter the refinement phase. Expensive cells are marked blue and red in Figure 6.

We start off with an instance of ACT that does not exceed the memory budget. When a training point hits an expensive cell, we determine the logical representation of this cell101010Recall that ACT uses a fanout of 256 thus the physical cell in the tree might not be the same as the logical cell that was originally inserted into the tree.. Then we verify for each of its four children cells whether they intersect, are fully contained in, or do not intersect the referenced polygons at all, and update ACT accordingly. We repeat this procedure until we have exhausted the memory budget. We show the effect of training the index in Section V.

Iii-E Probe Index

A lookup in ACT returns at most one cell that maps to a list of polygon references. The probe algorithm is shown in Listing 2. In contrast to a search in a binary tree, a search in a radix tree is comparison-free. This means that we do not compare the value of the search key to the value(s) stored in the current node. We only need to extract the relevant bits of the search key and jump to the corresponding offset. However, we do a comparison for checking whether the entry contains a payload.

[outerlinewidth=0.1pt,outerlinecolor=black, innerleftmargin=2pt,innerrightmargin=2pt,innertopmargin=-15pt,innerbottommargin=10pt]

  face node faceNode
  the cell id of the point cellId
  tagged entry taggedEntry
  rootNode = extract face bits of cellId and look up root node in faceNode
  if (rootNode does not exist)
    // no indexed cells on that face
    return invalid entry
  if (common prefix of rootNode does not match)
    return invalid entry
  level = 0
  currNode = rootNode
  bits = getBits(cellId, level++) // extract relevant bits
  // traverse the tree until we either hit the sentinel node or found a payload
  while (taggedEntry = currNode.getEntry(bits) is a pointer) {
    if (taggedEntry points to the sentinel node)
      return false hit
    currNode = taggedEntry
    bits = getBits(cellId, level++)
Listing 2: Probe Adaptive Cell Trie

For the tagged entry returned by the tree probe, we need to differentiate between (i) one payload, (ii) two payloads, and (iii) an offset. In the first case, we need to check whether the payload is invalid, which indicates a false hit. Otherwise, we extract the interior flag (the least significant bit of the 31 bit payload) and the polygon id and return a polygon reference. In the second case, we extract and return both references. Only in the third case, we need to access the lookup table to retrieve the polygon references.

Listing 3 shows the complete probe algorithm. For a given point, we retrieve the cell that it is contained in (if such a cell exists) and go over all references of this cell. We can skip refinement checks for true hits and immediately output the join partners.

[outerlinewidth=0.1pt,outerlinecolor=black, innerleftmargin=2pt,innerrightmargin=2pt,innertopmargin=-15pt,innerbottommargin=10pt] Listing 3: The probe algorithm input:   points points // lat/lng coordinates and cell ids   polygons polygons // lat/lng coordinates of vertices   face node faceNode   lookup table lookupTable output:   list of join pairs pairs // point/polygon pairs begin   for (point in points) {     taggedEntry = probeAdaptiveCellTrie(faceNode, point.cellId) //  cf., Listing 2     if (taggedEntry is invalid)       continue     references = getPolygonReferences(lookupTable, taggedEntry) //  returns a list of polygon references     for (reference in references) {       polygonId = reference.polygonId       polygon = polygons[polygonId]       if (reference is true hit) {         add {point, polygon} to pairs       } else { // candidate hit         if (polygonCoversPoint(polygon, point)) // exact PIP test           add {point, polygon} to pairs   } } } end

Iv Modern Hardware

In this section, we present the optimizations of our proposed join algorithm for modern hardware architectures. The general hardware trend is having an increasing number of cores with ever larger vector processing units (VPUs). That applies to established server-class machines such as Intel’s Xeon and particularly to high-performance computing platforms such as Intel’s Xeon Phi. In this work, we focus on Intel’s second-generation Xeon Phi processor, Knights Landing.

Iv-a Knights Landing Processor

Just like its predecessor, the KNL processor is a many integrated core (MIC) architecture which draws its computational power from wide VPUs. The corresponding AVX-512 offer SIMD capabilities on 512-bit registers. In contrast to the first generation Phi, KNL is available as a co-processor expansion card and as a self-boot processor (socket on a mainboard) which is fully backward compatible as it also implements all legacy instructions. In this work, we focus on the self-boot processor and not the co-processor expansion card.

Compared to our Xeon processor (cf., Section V), equipped with 14 fast clocked brawny cores per socket, the KNL processor consists of 64 wimpy cores moderately clocked at 1.3 GHz. The register width as well as the number of vector registers has been doubled over AVX2 architectures. Thus, for a software to run efficiently on KNL, it is essential to exploit the offered SIMD capabilities.

Especially for data-intensive workloads like ours, KNL offers three promising features: (i) it is equipped with an on-chip High Bandwidth Memory (HBM) which is comparable to the memory on high-end GPUs in speed and size, (ii) the processor performs a much more aggressive memory prefetching, and (iii) the 4-way hyper threading can help to hide memory latencies.

Iv-B Parallelism

KNL is a massively parallel processor in terms of thread-level parallelism and data parallelism on the instruction level.

Iv-B1 Multithreading

Thread-level parallelism can be achieved using well-known techniques. In existing parallel joins, multiple threads probe an index structure (e.g., a hash table) in parallel. During the probe phase, the index structure is in an immutable state. Therefore, an arbitrary number of threads can perform read operations without the need for synchronization. As no contention points exist, the join performance typically scales linearly in the number of threads.

Iv-B2 Data Parallelism

Achieving data parallelism on instruction level is more involved and requires several modifications over the scalar implementation. The goal is again to perform multiple index lookups in parallel. We divide a 512 bit SIMD register into eight 64 bit SIMD lanes and can thus process up to eight elements simultaneously per vector unit. The major difference to thread parallelism is, that the same instruction is executed on all in-flight lookups, which is due to the lockstep nature of SIMD. Thus, for the SIMD implementation to be efficient, the in-flight lookups need to make progress with every issued instruction. This implies that every lookup must be in the same stage of the algorithm, i.e., it is no longer possible that one lookup performs a prefix comparison while simultaneously another issues a memory load. To achieve this, our probe algorithm from Listing 2 is decomposed into the following stages:

Determine tree root: As mentioned earlier, our implementation uses S2 for space partitioning the earth. Thus, a join index actually consists of up to six radix trees. In this stage, we extract the three face bits of the query point’s cell id to determine the radix tree that needs to be traversed. Further, we check if the query point and the polygons indexed in the radix tree share the predetermined common prefix. If the prefixes do not match, then it is guaranteed that there is no join partner. We want to point out that on the one hand this stage of the algorithm is S2 specific but on the other hand it is similar to a scenario where the data is radix partitioned by three bits.

Tree traversal: In this stage, we traverse the tree. In each iteration, we compute the offset of the entry in the current tree node according to the search key. If the entry is a pointer to a child tree node we follow the pointer and continue.

Produce output: In the third and last stage, we interpret the content of the entries where the search in the previous stage terminated. We have to distinguish the cases (i) false hit, (ii) single or double hit, and (iii) multiple hits which are either true or candidate hits (cf., Section III).

To illustrate the vectorized implementation, we use pseudo code instead of C++ for better readability. In the pseudo code, vectors are annotated with an arrow, bitmasks are named —m—, and for masked operations we use the following notation:

=  —operation——(——)—

The —operation— is applied to the components in specified by mask and the results are stored in . The remaining components are copied from to . Further, we make use of the —vpermd/q— instruction (for better readability we refer to it as —permute—) that shuffles vector components using the corresponding index vector, e.g.,

= —permute(——)—.

During index construction we initialize the following vectors: [outerlinewidth=0.1pt,outerlinecolor=black, innerleftmargin=2pt,innerrightmargin=2pt,innertopmargin=3pt,innerbottommargin=10pt]

   // pointers to the six root nodes
   // common prefixes
   // common prefix lengths

These vectors are used to determine the tree roots in the first stage (cf., Listing 4). – For performance reasons these vectors remain in CPU registers to avoid memory loads. – Given a vector of cell ids, we first extract the face bits of each cell id. The resulting vector is then used as a permutation index vector to move the pointer of the corresponding root nodes to the output vector . Similarly, we obtain the prefixes and the lengths of the prefixes. The actual prefix check, an equality comparison, results in a bitmask which is used as an execution mask in the later stage. I.e., if the cell id in SIMD lane does not match the common prefix, then the -th bit in the bitmask is set to zero.

[outerlinewidth=0.1pt,outerlinecolor=black, innerleftmargin=2pt,innerrightmargin=2pt,innertopmargin=-15pt,innerbottommargin=10pt]

  cell id vector 
  tree node pointers 
  active lane mask m
   =  // extract the 3 face bits
   = permute(, ) // get the node pointers
  // get common prefix and prefix length for each tree
   = permute(, )
   = permute(, )
   =   // extract the common prefixes
  m =  ==  // prefix check (results in a bitmask)
Listing 4: Stage 1 of the vectorized probe algorithm

The tree traversal, shown in Listing 5, is more involved. As it performs multiple traversals in lockstep, it is required to do some bookkeeping about the current traversal state. The mask —m— keeps track of the active SIMD lanes and in —m— the algorithm memorizes the lanes that produced an output. The output mask therefore represents the execution mask for the next stage. The output values itself are stored in .

In each iteration of the while loop, the algorithm first computes the entry offsets from which it subsequently gathers the contents. The gather is a masked operation for two reasons: (i) to avoid unnecessary memory accesses and (ii) to keep already existing output values in .

Once the bucket contents are loaded into a vector register, we have to check whether the values are payloads or pointers. In case of pointers, the addresses are compared to the address of the sentinel node. If the comparison evaluates to true, the search terminates without producing an output, otherwise the corresponding pointers are moved to . Based on the two masks —m— and —m—, the traverse mask —m— is updated and the search continues until —m— is zero.

[outerlinewidth=0.1pt,outerlinecolor=black, innerleftmargin=2pt,innerrightmargin=2pt,innertopmargin=-15pt,innerbottommargin=10pt]

  cell id vector 
  node pointer vector 
  mask m
  output mask 
  output payloads 
  level = 1 // start tree traversal at level 1
  m =  m // active lanes
  m = 0 // lanes that produced an output
   = [0,  , 0] // a place for the produced output
  while (m != 0) {
      compute_bucket_offset( level)
     // load bucket contents
     // check for pointer (the two LSBs are zero)
     m = (( & 11 0)
     m = ( sentinel) // check for sentinel
      // update current node pointers
     // identify lanes that produced an output
     m & 
     // identify trapped lanes and update traversal mask
      ^=  & 
Listing 5: Stage 2 of the vectorized probe algorithm

In the last stage, we have to distinguish between true and candidate hits. We have not vectorized this stage since its performance impact is neglectable. However, for the —select polygon_id, count(*) … group by polygon_id— query we use in the evaluation (cf., Section V), we partially optimized the aggregation using SIMD. Otherwise, the aggregation would become a bottleneck due to the wimpy CPU cores of KNL.

Our SIMDfied aggregation is optimized for the most likely cases, where either one polygon id or two polygon ids are returned by the ACT lookup. We have implemented a direct aggregation based on gather and scatter instructions and the efficient popcount implementation presented in [5]. In the less likely cases, where the query point matches more than two polygons (i.e., ACT returns an offset into the lookup table) or where the match is inconclusive, we materialize the results into a small L1 resident buffer, which is consumed sequentially later on.

V Evaluation

In this section, we show the performance of our join and the selectivity of our index for a join between points and polygons on two different machines. We compare the performance of our technique to an R-tree based approach that uses the same code in the refinement phase as ours. To give a point of reference, we also compare against PostGIS and Magellan. In addition to the results of the exact algorithm, we show that an approximate algorithm can significantly speed up the join while providing sufficiently precise results for many use cases.

Configuration: We evaluate our approach on a server-class machine (Xeon) and on KNL (cf., Section IV-A). Xeon runs Ubuntu 16.04 and is equipped with an Intel Xeon E5-2680 v4 CPU (2.40 GHz, 3.30 GHz turbo) and 256 GB DDR4 RAM. The machine has two NUMA sockets with 14 physical cores each, resulting in a total of 28 physical cores (56 hyperthreads). KNL runs CentOS Linux release 7.2.1511 and is equipped with a single Intel Xeon Phi CPU 7210 (1.30 GHz, 1.50 GHz turbo), 96 GB DDR4 RAM, and 16 GB on-chip HBM. We use GCC version 5.4.0 with O3 enabled in all experiments.

We join 1.23 B points from the NYC taxi dataset111111 against NYC’s boroughs (5 polygons), neighborhoods (289 polygons), and census blocks (39,184 polygons) and count the number of points per polygon. While there are only five boroughs, their polygons are significantly more complex.

We compare our technique to an approach that builds an R-tree on the polygons, probes the points against the index, and refines candidate hits. In the refinement phase, we use S2’s PIP test, the same algorithm that our technique uses. In the filter phase, we use the boost R-tree implementation (1.6.0) and evaluate different splitting strategies (—linear—, —quadratic—, and —rstar—) and maximum element counts per node (4, 8, 16, and 32). —rstar— with a maximum of 8 elements per node performs best in all workloads. We therefore omit the results for the other configurations. We have also evaluated the —STRtree— implementation in GEOS (3.5.0), however, omit it in the further evaluation as it cannot compete with —rstar— for our workloads. For example in the neighborhoods join on Xeon, —STRtree— achieves a throughput of 27.4 M points/s, while —rstar— processes 70.9 M points/s. The performance difference is caused by the selectivity of the respective index structure. —rstar— returns 1.60 candidates on average while —STRtree— returns 1.94. In addition, we benchmark PostgreSQL 9.6.1 (PostGIS 2.3.1) with a GiST index on the polygons. We configured PostgreSQL to use as many parallel workers121212The intra-query parallelism feature was introduced in PostgreSQL 9.6. as available hyperthreads. As another competitor, we evaluate a single-machine deployment of Magellan. Similar to our approach, Magellan uses true hit filtering. Our default configuration for computing the individual coverings is as follows: max covering cells = 128, max covering level = 30, max interior cells = 256, and max interior level = 20.

(a) Dual-socket Xeon
(b) Single-socket KNL
Fig. 8: Throughput of our approach (ACT) compared to the best competitor (boost R-tree) on the Xeon and on the KNL machine (approx100m and approx10m = approximate with 100 and 10 m precision, exact = refine candidate hits)

Results: Figure 8 shows the throughput of our approach compared to the boost R-tree using all threads (56 and 256 threads on Xeon and KNL, respectively).

approx always skips the refinement phase and treats all candidate hits as true hits. The precision of approx is bounded by the diagonal length of the largest covering cell. We refined the index to only contain sufficiently small covering cells to guarantee the denoted precision. approx10m has a false positive rate of almost 0% for boroughs, 3% for neighborhoods, and 17% for the census blocks dataset. The average distance of these false positives join partners is 1.31 m for boroughs, 1.62 m for neighborhoods, and 1.50 m for the census blocks dataset. This shows that our technique can achieve a significantly better precision than the user-defined precision bound (10 m in that case) for real workloads. approx10m consumes 1.2 GB of memory for the largest of the three datasets.

exact denotes the full join that refines candidate hits using S2. We train the index with 1 M historical data points for all datasets.

Fig. 9: Throughput after training ACT on boroughs with an increasing number of historical data points

Figure 9 shows the throughput increase when training ACT with an increasing number of points. An untrained ACT achieves a throughput of 2360 M points/s, whereas ACT trained with 1 M points achieves 3732 M points/s (+58%).

The Xeon binaries are compiled with AVX2 optimizations (march=core-avx2) while the KNL implementation is hand-optimized (cf., Section IV). Compared to an auto-vectorized version, the performance of our hand-optimized implementation is between 29% and 70% higher. Auto-vectorization with AVX-512 (march=knl) instead of AVX2 does not have a performance impact suggesting that hand-tuning is essential to fully utilize the wide vector processing units of KNL.

On Xeon (Figure 7(a)), ACT-exact achieves a throughput of 3735 M points/s for the boroughs dataset. With a higher number of polygons in the neighborhoods and census datasets, the throughput of ACT-exact decreases. The reason for the slowdown is the lower selectivity and the higher space consumption of our index for the larger datasets (cf., Table I). While the index of the boroughs dataset fits into the cache, the index of the neighborhoods and census blocks datasets exceed the cache size. The R-tree achieves a peak performance of 61.2 M points/s for the neighborhoods dataset. The reason for its slow performance for the boroughs dataset is as follows: S2’s PIP test implements the ray tracing algorithm. As mentioned earlier, this algorithm has a runtime complexity of with being the number of edges. Since the boroughs are complex polygons with many edges, the PIP tests in the refinement phase are very expensive. Here, our algorithm shines since it can identify most join partners in the filter phase and only needs to enter the refinement phase for 0.1% of the points.

As a point of reference, PostgreSQL (PostGIS) with a GiST index on the polygons only achieves a throughput of 0.39 M points/s for the boroughs, 1.09 M points/s for the neighborhoods, and 0.69 M points/s for the census blocks dataset. Similar to the R-tree based join, PostgreSQL’s performance suffers from the complex polygons in the boroughs dataset. Even though PostgreSQL parallelizes this query, its performance is still orders of magnitude lowers than ours. Magellan performs better than PostGIS and achieves a throughput of 0.88 M points/s for boroughs, 4.57 M points/s for neighborhoods, and 2.24 M points/s for census blocks.

With 10 meter precision, ACT achieves a throughput of 4532 M points/s for the boroughs dataset and still 874 M points/s for the almost 8000x larger census blocks dataset. The reason for the similar performance of ACT-approx100m and ACT-approx10m is that neither of the two performs any exact checks. Thus, their performance is dominated by the costs for the ACT node accesses and the final aggregation. For the census blocks dataset, approx10m is even faster than approx100m, since the finer-grained index of approx10m returns less candidates on average leading to less aggregations.

On KNL (Figure 7(b)), ACT-exact achieves a peak throughput of 4436 M points/s for the boroughs dataset outperforming the Xeon counterpart. For the census blocks dataset, we can see the same effect as on the Xeon machine with the more precise approx approach achieving the highest throughput. These results show that the single-socket KNL can compete with or even outperform a dual-socket Xeon machine when the code uses vector instructions.

metric boroughs neighborhoods census
tree nodes 413 13025 590147
false hits 0.07% 0.09% 0.07%
solely true hits 99.9% 87.1% 72.1%
candidates 1.06 1.66 1.72
TABLE I: Metrics of our index
metric boroughs neighborhoods census
tree nodes 456 22484 634005
false hits 0.07% 0.11% 0.07%
solely true hits 99.9% 97.7% 88.6%
candidates 1.15 1.57 1.70
TABLE II: Effect of training the index
(a) Dual-socket Xeon
(b) Single-socket KNL
Fig. 10: Scalability for the census blocks dataset on the Xeon and on the KNL machine

Table I shows different metrics of our index for the three polygon datasets. The quality metric solely true hits indicates the percentage of points that skipped the expensive refinement phase, which is clearly above 70% throughout all datasets. The points for which the index probe neither returns a false hit nor solely true hits and thus need to enter the refinement phase are compared against less than two candidate polygons on average. For the boroughs dataset, our index consumes less than one MB of memory and thus the top levels of the tree fit into the L2 cache of the Xeon machine. For the larger datasets, we use a larger number of cells and thus also a larger index structure. The exact index of the census blocks dataset requires more than one GB of memory. The high memory consumption is caused by the high fanout of 256 of our radix tree. An adaptive radix tree would be more space-efficient, however, switching between the different node types has a significant impact on probe performance, especially for cache-resident indices.

Table II shows the effect of training the index with 1 M training points (cf., Section III-B). This optimization improved the percentage of solely true hits for all three datasets (especially for census blocks) while increasing the number of tree nodes and thus the space consumption of the indices.

Finally, we evaluate the scalability of our technique on both machines. We choose the census blocks dataset for this experiment since its index significantly exceeds the cache sizes of both machines but still fits into the MCDRAM of KNL. Figure 10 shows the results.

On Xeon, all three approaches scale almost linearly with the number of physical cores and benefit from hyperthreading. The fact that an oversubscription of cores has a positive performance impact shows that our technique is bound by memory access latencies and having more threads than physical cores can hide these latencies.

On KNL, ACT-approx10 has its maximum performance with 2-way hyperthreading (128 threads) and decreases with 4-way hyperthreading (256 threads). ACT-exact scales from 13.0 M points/s with one thread to 417 M points/s with 56 threads on Xeon. On KNL, ACT-exact scales from 25.9 M points/s with eight threads to 260 M points/s with 256 threads.

Vi Related Work

Geospatial joins have been studied for decades and there is a lot of related work on algorithmic techniques. [6] gives a summary of spatial join techniques. [1] proposes true hit filtering in the form of maximum enclosed rectangles and circles allowing to skip the refinement phase in many cases. [7] follows up on this approach by using raster approximations in the form of uniform grids thereby improving the selectivity. Both approaches optimize not only for selectivity but also for I/O operations, which is an important performance factor for disk-based database systems. [2] recursively divides the MBR of a polygon into four cells until a certain granularity is reached, identifies interior cells, and indexes them in an R-tree to skip refinement checks during join processing. [8] presents an approach that leverages multiple grids with different resolutions to index geospatial objects using standard database indices.

In contrast to existing techniques, our approach proposes a new way of processing geospatial joins by identifying the majority or even all of the join partners in the filter phase using an adaptive grid while leveraging the large main-memory capacities, high memory bandwidths, multi-cores, and wide vector processing units of modern hardware. None of these techniques supports training the index with historical data points to improve probe performance.

Several database systems and geographical information systems (GIS) support geospatial joins. [9]

present a comprehensive survey of state-of-the-art systems and classifies the approaches that exist in literature to design and implement such systems. There are two main geospatial processing libraries used by database systems: GEOS

131313 and S2 (cf., Section II). MongoDB, a document-oriented database, uses S2 for geospatial processing and indexes the coverings of geospatial objects in a B-tree141414 HyPerSpace [3], a geospatial extension to the main-memory database system HyPer [10], also uses S2 for geospatial processing. PostGIS151515, a geospatial extension to PostgreSQL161616, uses GEOS for geospatial processing and an R-tree implemented on top of GiST [11] for indexing geospatial objects. MonetDB, a column-oriented database, also uses GEOS for geospatial processing171717 Magellan, a library for geospatial analytics that is based on Apache Spark, indexes geospatial objects using the Z curve. Similar to our technique, it is capable of identifying join partners in the filter phase. However, it uses a uniform grid and thus cannot adapt to the geometrical features of the polygons and the expected distribution of query points. Two popular Hadoop-based geospatial information systems are Hadoop-GIS [12] and SpatialHadoop [13]. Both systems partition the data and store it in HDFS blocks such that spatial objects that belong to a particular partition are stored in the same HDFS block and insert the MBRs of the partitions into a global index. Since these systems rely on offline partitioning of the data, they cannot efficiently handle the online case where points are streamed. GeoSpark [14] does not support the online case either since it does not allow to perform the join without an index on points, and Simba [15] does not yet support polygons. Most work on hardware optimizations for geospatial joins focuses on GPU offloading[16, 17, 18] and [19] propose a vision for a GPU-accelerated end-to-end system for performing spatial computations which decides whether to use CPUs or GPUs based on a cost model.

Vii Conclusions

We have presented an adaptive geospatial join that leverages true hit filtering using a hierarchical grid represented by a specialized radix tree. We have transformed a traditionally compute-intensive problem into a memory-intensive one. Our approach is enabled by adaptively using the large high-bandwidth main memory of modern hardware. We have shown that it is possible to refine the index up to a user-defined precision and identify all join partners in the filter phase. We have demonstrated that the exact version of our algorithm can adapt to the expected point distribution. We have optimized our implementation for the brand new KNL processor and have shown that our approach can outperform existing techniques by up to two orders of magnitude.

In future work, we plan to evaluate our approach on the next-generation Intel Xeon processor, named Skylake, which also supports AVX-512, and on GPUs. We also plan to address offline workloads by leveraging small materialized aggregates (SMAs) available in Data Blocks [20].


We want to thank Ram Sriharsha for his help in evaluating Magellan. This work has been sponsored by the German Federal Ministry of Education and Research (BMBF) grant 01IS12057 (FASTDATA and MIRIN). This work is further part of the TUM Living Lab Connected Mobility (TUM LLCM) project and has been funded by the Bavarian Ministry of Economic Affairs and Media, Energy and Technology (StMWi) through the Center Digitisation.Bavaria, an initiative of the Bavarian State Government.


  • [1] T. Brinkhoff, H. Kriegel, R. Schneider, and B. Seeger, “Multi-step processing of spatial joins,” in Proc. of SIGMOD, 1994, pp. 197–208. [Online]. Available:
  • [2] K. V. R. Kanth and S. Ravada, “Efficient processing of large spatial queries using interior approximations,” in Proc. of SSTD, 2001, pp. 404–424. [Online]. Available:
  • [3] V. Pandey, A. Kipf, D. Vorona, T. Mühlbauer, T. Neumann, and A. Kemper, “High-performance geospatial analytics in HyPerSpace,” in Proc. of SIGMOD.   ACM, 2016, pp. 2145–2148.
  • [4] V. Leis, A. Kemper, and T. Neumann, “The adaptive radix tree: ARTful indexing for main-memory databases,” in Proc. of ICDE, 2013, pp. 38–49. [Online]. Available:
  • [5] T. Gubner and P. Boncz, “Exploring Query Compilation Strategies for JIT, Vectorization and SIMD,” in ADMS, 2017.
  • [6] E. H. Jacox and H. Samet, “Spatial join techniques,” ACM Trans. Database Syst., vol. 32, no. 1, p. 7, 2007. [Online]. Available:
  • [7] G. Zimbrao and J. M. de Souza, “A raster approximation for processing of spatial joins,” in Proc. of VLDB, 1998, pp. 558–569. [Online]. Available:
  • [8] J. Roth, “The extended split index to efficiently store and retrieve spatial data with standard databases,” in Proc. of IADIS, 2009, pp. 85–92.
  • [9] A. Eldawy and M. F. Mokbel, “The era of big spatial data,” in ICDE Workshops, 2015, pp. 42–49. [Online]. Available:
  • [10] A. Kemper and T. Neumann, “HyPer: A hybrid OLTP&OLAP main memory database system based on virtual memory snapshots,” in Proc. of ICDE, 2011, pp. 195–206. [Online]. Available:
  • [11] J. M. Hellerstein, J. F. Naughton, and A. Pfeffer, “Generalized search trees for database systems,” in Proc. of VLDB, 1995, pp. 562–573. [Online]. Available:
  • [12] A. Aji, F. Wang, H. Vo, R. Lee, Q. Liu, X. Zhang, and J. H. Saltz, “Hadoop-GIS: A high performance spatial data warehousing system over MapReduce,” PVLDB, vol. 6, no. 11, pp. 1009–1020, 2013. [Online]. Available:
  • [13] A. Eldawy, “SpatialHadoop: Towards flexible and scalable spatial processing using MapReduce,” in SIGMOD, PhD Symposium, 2014, pp. 46–50. [Online]. Available:
  • [14] J. Yu, J. Wu, and M. Sarwat, “Geospark: a cluster computing framework for processing large-scale spatial data,” in Proc. of SIGSPATIAL, 2015, pp. 70:1–70:4. [Online]. Available:
  • [15] D. Xie, F. Li, B. Yao, G. Li, L. Zhou, and M. Guo, “Simba: Efficient in-memory spatial analytics,” in Proc. of SIGMOD, 2016, pp. 1071–1085. [Online]. Available:
  • [16] N. Bandi, C. Sun, D. Agrawal, and A. El Abbadi, “Fast computation of spatial selections and joins using graphics hardware,” Inf. Syst., vol. 32, no. 8, pp. 1073–1100, 2007. [Online]. Available:
  • [17] S. You, J. Zhang, and L. Gruenwald, “Parallel spatial query processing on gpus using r-trees,” in SIGSPATIAL Workshops, 2013, pp. 23–31. [Online]. Available:
  • [18] J. Zhang, S. You, and L. Gruenwald, “Parallel selectivity estimation for optimizing multidimensional spatial join processing on gpus,” in Proc. of ICDE, 2017, pp. 1591–1598. [Online]. Available:
  • [19] H. Chavan, R. Alghamdi, and M. F. Mokbel, “Towards a GPU accelerated spatial computing framework,” in ICDE Workshops, 2016, pp. 135–142. [Online]. Available:
  • [20] H. Lang, T. Mühlbauer, F. Funke, P. A. Boncz, T. Neumann, and A. Kemper, “Data Blocks: Hybrid OLTP and OLAP on compressed storage using both vectorization and compilation,” in Proc. of SIGMOD, 2016, pp. 311–326. [Online]. Available: