Accelerating the XGBoost algorithm using GPU computing

Decision trees

Decision tree learning is a method of predictive modelling that learns a model by repeatedly splitting subsets of the training examples (also called instances) according to some criteria. Decision tree inducers are supervised learners that accept labelled training examples as an input and generate a model that may be used to predict the labels of new examples.

In order to construct a decision tree, we start with the full set of training instances and evaluate all possible ways of creating a binary split among those instances based on the input features in $$\overrightarrow{x}$$. We choose the split that produces the most meaningful separation of the target label y. Different measures can be used to evaluate the quality of a split. After finding the ‘best’ split, we can create a node in the tree that partitions training instances down the left or right branch according to some feature value. The subsets of training instances can then be recursively split to continue growing the tree to some maximum depth or until the quality of the splits is below some threshold. The leaves of the tree will contain predictions for the target label y. For categorical labels, the prediction can be set as the majority class from the training instances that end up in that leaf. For regression tasks, the label prediction can be set as the mean of the training instances in that leaf.

To use the tree for prediction, we can input an unlabelled example at the root of the tree and follow the decision rules until the example reaches a leaf. The unlabelled example can be labelled according to the prediction of that leaf.

Figure 1 shows an example decision tree that can predict whether or not an individual owns a house. The decision is based on their age and whether or not they have a job. The tree correctly classifies all instances from Table 1.

Decision tree algorithms typically expand nodes from the root in a greedy manner in order to maximise some criterion measuring the value of the split. For example, decision tree algorithms from the C4.5 family (Quinlan, 2014), designed for classification, use information gain as the split criterion. Information gain describes a change in entropy H from some previous state to a new state. Entropy is defined as $$H(T)=-{\displaystyle \sum _{y\in Y}P}(y){\text{log}}_{b}P(y)$$ where T is a set of labelled training instances, y ∈ Y is an instance label and P(y) is the probability of drawing an instance with label y from T. Information gain is defined as $$\text{IG}(T,T_{\text{left}},T_{\text{right}})=H_T-(n_{\text{left}}/n_{\text{total}})*H(T_{\text{left}})-(n_{\text{right}}/n_{\text{total}})*H(T_{\text{right}})$$ Here T_left and T_right are the subsets of T created by a decision rule. n_total, n_left and n_right refer to the number of examples in the respective sets.

Many different criteria exist for evaluating the quality of a split. Any function can be used that produces some meaningful separation of the training instances with respect to the label being predicted.

In order to find the split that maximises our criterion, we can enumerate all possible splits on the input instances for each feature. In the case of numerical features and assuming the data has been sorted, this enumeration can be performed in O(nm) steps, where n is the number of instances and m is the number of features. A scan is performed from left to right on the sorted instances, maintaining a running sum of labels as the input to the gain calculation. We do not consider the case of categorical features in this paper because XGBoost encodes all categorical features using one-hot encoding and transforms them into numerical features.

Another consideration when building decision trees is how to perform regularisation to prevent overfitting. Overfitting on training data leads to poor model generalisation and poor performance on test data. Given a sufficiently large decision tree it is possible to generate unique decision rules for every instance in the training set such that each training instance is correctly labelled. This results in 100% accuracy on the training set but may perform poorly on new data. For this reason it is necessary to limit the growth of the tree during construction or apply pruning after construction.

Gradient boosting

Decision trees produce interpretable models that are useful for a variety of problems, but their accuracy can be considerably improved when many trees are combined into an ensemble model. For example, given an input instance to be classified, we can test it against many trees built on different subsets of the training set and return the mode of all predictions. This has the effect of reducing classifier error because it reduces variance in the estimate of the classifier.

Figure 2 shows an ensemble of two decision trees. We can predict the output label using all trees by taking the most common class prediction or some weighted average of all predictions.

Ensemble learning methods can also be used to reduce the bias component in the classification error of the base learner. Boosting is an ensemble method that creates ensemble members sequentially. The newest member is created to compensate for the instances incorrectly labelled by the previous learners.

Gradient boosting is a variation on boosting which represents the learning problem as gradient descent on some arbitrary differentiable loss function that measures the performance of the model on the training set. More specifically, the boosting algorithm executes M boosting iterations to learn a function F(x) that outputs predictions $\widehat{y}=F(x)$ minimising some loss function $L(y,\widehat{y})$. At each iteration we add a new estimator f(x) to try to correct the prediction of y for each training instance: $$F_{m+1}(x)=F_m(x)+f(x)=y$$ We can correct the model by setting f(x) to: $$f(x)=y-F_m(x)$$ This fits the model f(x) for the current boosting iteration to the residuals y − F_m(x) of the previous iteration. In practice, we approximate f(x), for example by using a depth-limited decision tree.

This iterative process can be shown to be a gradient descent algorithm when the loss function is the squared error: $$L(y,F(x))=\frac{1}{2}{(y-F(x))}^2$$ To see this, consider that the loss over all training instances can be written as $$J={\displaystyle \sum _{i}L}(y_i,F(x_i))$$ We seek to minimise J by adjusting F(x_i). The gradient for a particular instance x_i is given by $$\frac{dJ}{dF(x_i)}=\frac{d{\displaystyle \sum _{i}L}(y_i,F(x_i))}{dF(x_i)}=\frac{dL(y_i,F(x_i))}{dF(x_i)}=F_m(x_i)-y_i$$ We can see that the residuals are the negative gradient of the squared error loss function: $$f(x)=y-F_m(x)=-\frac{dL(y,F(x))}{dF(x)}$$ By adding a model that approximates this negative gradient to the ensemble we move closer to a local minimum of the loss function, thus implementing gradient descent.

Generalised gradient boosting and XGBoost

Herein, we derive the XGBoost algorithm following the explanation in Chen & Guestrin (2016). XGBoost is a generalised gradient boosting implementation that includes a regularisation term, used to combat overfitting, as well as support for arbitrary differentiable loss functions.

Instead of optimising plain squared error loss, an objective function with two parts is defined, a loss function over the training set as well as a regularisation term which penalises the complexity of the model: $$\text{Obj}={\displaystyle \sum _{i}L}(y_i,{\widehat{y}}_i)+{\displaystyle \sum _k\Omega }(f_k)$$ $$L(y_i,{\widehat{y}}_i)$$ can be any convex differentiable loss function that measures the difference between the prediction and true label for a given training instance. Ω(f_k) describes the complexity of tree f_k and is defined in the XGBoost algorithm (Chen & Guestrin, 2016) as (1) $$\Omega (f_k)=\gamma T+\frac{1}{2}\lambda w^2$$ where T is the number of leaves of tree f_k and w is the leaf weights (i.e. the predicted values stored at the leaf nodes). When Ω(f_k) is included in the objective function we are forced to optimise for a less complex tree that simultaneously minimises $L(y_i,{\widehat{y}}_i)$. This helps to reduce overfitting. γT provides a constant penalty for each additional tree leaf and λw² penalises extreme weights. γ and λ are user configurable parameters.

Given that boosting proceeds in an iterative manner we can state the objective function for the current iteration m in terms of the prediction of the previous iteration ${\widehat{y}}_i{}^{(m-1)}$ adjusted by the newest tree f_k: $${\text{Obj}}^m={\displaystyle \sum _{i}L}(y_i,{\widehat{y}}_i{}^{(m-1)}+f_k(x_i))+{\displaystyle \sum _k\Omega }(f_k)$$ We can then optimise to find the f_k which minimises our objective.

Taking the Taylor expansion of the above function to the second order allows us to easily accommodate different loss functions: $${\text{Obj}}^m\simeq {\displaystyle \sum }_i[L(y_i,{\widehat{y}}_i{}^{(m-1)})+g_i{f}_k(x)+\frac{1}{2}{h}_i{f}_k{(x)}^2]+{\displaystyle \sum }_k\Omega (f_k)+\text{constant}$$ Here, g_i and h_i are the first and second order derivatives respectively of the loss function for instance i: $$g_i=\frac{dL(y_i,{\widehat{y}}_i{}^{(m-1)})}{d{\widehat{y}}_i{}^{(m-1)}}{h}_i=\frac{d^{2}L(y_i,{\widehat{y}}_i{}^{(m-1)})}{d{({\widehat{y}}_i{}^{(m-1)})}^2}$$ Note that the model ${\widehat{y}}_i{}^{(m-1)}$ is left unchanged during this optimisation process. The simplified objective function with constants removed is $${\text{Obj}}^m={\displaystyle \sum }_i[g_i{f}_k(x)+\frac{1}{2}{h}_i{f}_k{(x)}^2]+{\displaystyle \sum }_k\Omega (f_k)$$ We can also make the observation that a decision tree predicts constant values within a leaf. f_k(x) can then be represented as w_q(x) where w is the vector containing scores for each leaf and q(x) maps instance x to a leaf.

The objective function can then be modified to sum over the tree leaves and the regularisation term from Eq. (1): $${\text{Obj}}^m={\displaystyle \sum _{j=1}^T\left[\left({\displaystyle \sum _{i\in I_j}{g}_i}\right)w_{q(x)}+\frac{1}{2}\left({\displaystyle \sum _{i\in I_j}{h}_i}\right)w_{q(x)}^2\right]}+\gamma T+\frac{1}{2}\lambda {\displaystyle \sum _{j=1}^T{w}^2}$$ Here, I_j refers to the set of training instances in leaf j. The sums of the derivatives in each leaf can be defined as follows: $$G_j={\displaystyle \sum _{i\in I_j}{g}_i}{H}_j={\displaystyle \sum _{i\in I_j}{h}_i}$$ Also note that w_q(x) is a constant within each leaf and can be represented as w_j. Simplifying we get (2) $${\text{Obj}}^m={\displaystyle \sum _{j=1}^T\left[G_j{w}_j+\frac{1}{2}(H_j+\lambda )w_j^2\right]}+\gamma T$$ The weight w_j for each leaf minimises the objective function at $$\frac{\partial {\text{Obj}}^m}{\partial w_j}=G_j+(H_j+\lambda )w_j=0$$ The best leaf weight w_j given the current tree structure is then $$w_j=-\frac{G_j}{H_j+\text{λ}}$$ Using the best w_j in Eq. (2), the objective function for finding the best tree structure then becomes (3) $${\text{Obj}}^m=-\frac{1}{2}{\displaystyle \sum _{j=1}^T}\frac{G_j^2}{H_j+\text{λ}}+\text{γ}T$$ Eq. (3) is used in XGBoost as a measure of the quality of a given tree.

Growing a tree

Given that it is intractable to enumerate through all possible tree structures, we greedily expand the tree from the root node. In order to evaluate the usefulness of a given split, we can look at the contribution of a single leaf node j to the objective function from Eq. (3): $${\text{Obj}}_{\text{leaf}}=-\frac{1}{2}\frac{G_j^2}{H_j+\text{λ}}+\text{γ}$$ We can then consider the contribution to the objective function from splitting this leaf into two leaves: $${\text{Obj}}_{\text{split}}=-\frac{1}{2}\left(\frac{G_{j\text{L}}^2}{H_{j\text{L}}+\text{λ}}+\frac{G_{j\text{R}}^2}{H_{j\text{R}}+\text{λ}}\right)+2\text{γ}$$ The improvement to the objective function from creating the split is then defined as $$\text{Gain}={\text{Obj}}_{\text{leaf}}-{\text{Obj}}_{\text{split}}$$ which yields (4) $$\text{Gain}=\frac{1}{2}\left[\frac{G_{\text{L}}^2}{H_{\text{L}}+\text{λ}}+\frac{G_{\text{R}}^2}{H_{\text{R}}+\text{λ}}-\frac{{(G_{\text{L}}+G_{\text{R}})}^2}{H_{\text{L}}+H_{\text{R}}+\text{λ}}\right]-\text{γ}$$ The quality of any given split separating a set of training instances is evaluated using the gain function in Eq. (4). The gain function represents the reduction in the objective function from Eq. (3) obtained by taking a single leaf node j and partitioning it into two leaf nodes. This can be thought of as the increase in quality of the tree obtained by creating the left and right branch as compared to simply retaining the original node. This formula is applied at every possible split point and we expand the split with maximum gain. We can continue to grow the tree while this gain value is positive. The γ regularisation cost at each leaf will prevent the tree arbitrarily expanding. The split point selection is performed in O(nm) time (given n training instances and m features) by scanning left to right through all feature values in a leaf in sorted order. A running sum of G_L and H_L is kept as we move from left to right, as shown in Table 3. G_R and H_R are inferred from this running sum and the node total.

Table 2 shows an example set of instances in a leaf. We can assume we know the sums G and H within this node as these are simply the G_L or G_R from the parent split. Therefore, we have everything we need to evaluate Gain for every possible split within these instances and select the best.

XGBoost: data format

Tabular data input to a machine learning library such as XGBoost or Weka (Hall et al., 2009) can be typically described as a matrix with each row representing an instance and each column representing a feature as shown in Table 3. If f2 is the feature to be predicted then an input training pair $({\overrightarrow{x}}_i\mathrm{,}{y}_i)$ takes the form ((f0_i, f1_i), f2_i) where i is the instance id. A data matrix within XGBoost may also contain missing values. One of the key features of XGBoost is the ability to store data in a sparse format by implicitly keeping track of missing values instead of physically storing them. While XGBoost does not directly support categorical variables, the ability to efficiently store and process sparse input matrices allows us to process categorical variables through one-hot encoding. Table 4 shows an example where a categorical feature with three values is instead encoded as three binary features. The zeros in a one-hot encoded data matrix can be stored as missing values. XGBoost users may specify values to be considered as missing in the input matrix or directly input sparse formats such as libsvm files to the algorithm.

XGBoost: handling missing values

Representing input data using sparsity in this way has implications on how splits are calculated. XGBoost’s default method of handling missing data when learning decision tree splits is to find the best ‘missing direction’ in addition to the normal threshold decision rule for numerical values. So a decision rule in a tree now contains a numeric decision rule such as f0 ≤ 5.53, but also a missing direction such as missing = right that sends all missing values down the right branch. For a one-hot encoded categorical variable where the zeros are encoded as missing values, this is equivalent to testing ‘one vs all’ splits for each category of the categorical variable.

The missing direction is selected as the direction which maximises the gain from Eq. (4). When enumerating through all possible split values, we can also test the effect on our gain function of sending all missing examples down the left or right branch and select the best option. This makes split selection slightly more complex as we do not know the gradient statistics of the missing values for any given feature we are working on, although we do know the sum of all the gradient statistics for the current node. The XGBoost algorithm handles this by performing two scans over the input data, the second being in the reverse direction. In the first left to right scan the gradient statistics for the left direction are the scan values maintained by the scan, the gradient statistics for the right direction are the sum gradient statistics for this node minus the scan values. Hence, the right direction implicitly includes all of the missing values. When scanning from right to left, the reverse is true and the left direction includes all of the missing values. The algorithm then selects the best split from either the forwards or backwards scan.

Languages and libraries

The two main languages for general purpose GPU programming are CUDA and OpenCL. CUDA was chosen for the implementation discussed in this paper due to the availability of optimised and production ready libraries. The GPU tree construction algorithm would not be possible without a strong parallel primitives library. We make extensive use of scan, reduce and radix sort primitives from the CUB (Merrill & NVIDIA-Labs, 2016) and Thrust (Hoberock & Bell, 2017) libraries. These parallel primitives are described in detail in ‘Parallel primitives.’ The closest equivalent to these libraries in OpenCL is the boost compute library. Several problems were encountered when attempting to use Boost Compute and the performance of its sorting primitives lagged considerably behind those of CUB/Thrust. At the time of writing this paper OpenCL was not a practical option for this type of algorithm.

Execution model

CUDA code is written as a kernel to be executed by many thousands of threads. All threads execute the same kernel function but their behaviour may be distinguished through a unique thread ID. Listing 1 shows an example of kernel adding values from two arrays into an output array. Indexing is determined by the global thread ID and any unused threads are masked off with a branch statement.

Threads are grouped according to thread blocks that typically each contain some multiple of 32 threads. A group of 32 threads is known as a warp. Thread blocks are queued for execution on hardware streaming multiprocessors. Streaming multiprocessors switch between different warps within a block during program execution in order to hide latency. Global memory latency may be hundreds of cycles and hence it is important to launch sufficiently many warps within a thread block to facilitate latency hiding.

A thread block provides no guarantees about the order of thread execution unless explicit memory synchronisation barriers are used. Synchronisation across thread blocks is not generally possible within a single kernel launch. Device-wide synchronisation is achieved by multiple kernel launches. For example, if a global synchronisation barrier is required within a kernel, the kernel must be separated into two distinct kernels where synchronisation occurs between the kernel launches.

Memory architecture

CUDA exposes three primary tiers of memory for reading and writing. Device-wide global memory, thread block accessible shared memory and thread local registers.

• Global memory: Global memory is accessible by all threads and has the highest latency. Input data, output data and large amounts of working memory are typically stored in global memory. Global memory can be copied from the device (i.e. the GPU) to the host computer and vice versa. Bandwidth of host/device transfers is much slower than that of device/device transfers and should be avoided if possible. Global memory is accessed in 128 byte cache lines on current GPUs. Memory accesses should be coalesced in order to achieve maximum bandwidth. Coalescing refers to the grouping of aligned memory load/store operations into a single transaction. For example, a fully coalesced memory read occurs when a warp of 32 threads loads 32 contiguous 4 byte words (128 bytes). Fully uncoalesced reads (typical of gather operations) can limit device bandwidth to less than 10% of peak bandwidth (Harris, 2013).

• Shared memory: 48 KB of shared memory is available to each thread block. Shared memory is accessible by all threads in the block and has a significantly lower latency than global memory. It is typically used as working storage within a thread block and sometimes described as a ‘programmer-managed cache.’

• Registers: A finite number of local registers is available to each thread. Operations on registers are generally the fastest. Threads within the same warp may read/write registers from other threads in the warp through intrinsic instructions such as shuffle or broadcast (Nvidia, 2017).

Reduction

A parallel reduction reduces an array of values into a single value using a binary-associative operator. Given a binary-associative operator ⊕ and an array of elements the reduction returns $(a_0\oplus a_1\oplus \cdots \oplus a_{n-1})$. Note that floating point addition is not strictly associative. This means a sequential reduction operation will likely result in a different answer to a parallel reduction (the same applies to the scan operation described below). This is discussed in greater detail in ‘Floating point precision.’ The reduction operation is easy to implement in parallel by passing partial reductions up a tree, taking O(logn) iterations given n input items and n processors. This is illustrated in Fig. 3.

In practice, GPU implementations of reductions do not launch one thread per input item but instead perform parallel reductions over ‘tiles’ of input items then sum the tiles together sequentially. The size of a tile varies according to the optimal granularity for a given hardware architecture. Reductions are also typically tiered into three layers: warp, block and kernel. Individual warps can very efficiently perform partial reductions over 32 items using shuffle instructions introduced from Nvidia’s Kepler GPU architecture onwards. As smaller reductions can be combined into larger reductions by simply applying the binary associative operator on the outputs, these smaller warp reductions can be combined together to get the reduction for the entire tile. The thread block can iterate over many input tiles sequentially, summing the reduction from each. When all thread blocks are finished the results from each are summed together at the kernel level to produce the final output. Listing 2 shows code for a fast warp reduction using shuffle intrinsics to communicate between threads in the same warp. The ‘shuffle down’ instruction referred to in Listing 2 simply allows the current thread to read a register value from the thread d places to the left, so long as that thread is in the same warp. The complete warp reduction algorithm requires five iterations to sum over 32 items.

Reductions are highly efficient operations on GPUs. An implementation is given in Harris (2007) that approaches the maximum bandwidth of the device tested.

Parallel prefix sum (scan)

The prefix sum takes a binary associative operator (most commonly addition) and applies it to an array of elements. Given a binary associative operator ⊕ and an array of elements the prefix sum returns $[a_0,(a_0\oplus a_1),\mathrm{...},(a_0\oplus a_1\oplus \mathrm{...}\oplus a_{n-1})]$. A prefix sum is an example of a calculation which seems inherently serial but has an efficient parallel algorithm: the Blelloch scan algorithm.

Let us consider a simple implementation of a parallel scan first, as described in Hillis & Steele (1986). It is given in Algorithm 1. Figure 4 shows it in operation: we apply a simple scan with the addition operator to an array of 1’s. Given one thread for each input element the scan takes ${{\displaystyle \text{log}}}_{2}n=3$ iterations to complete. The algorithm performs O(nlog₂n) addition operations.

Given that a sequential scan performs only n addition operations, the simple parallel scan is not work efficient. A work efficient parallel algorithm will perform the same number of operations as the sequential algorithm and may provide significantly better performance in practice. A work efficient algorithm is described in Blelloch (1990). The algorithm is separated into two phases, an ‘upsweep’ phase similar to a reduction and a ‘downsweep’ phase. We give pseudocode for the upsweep (Algorithm 2) and downsweep (Algorithm 3) phases by following the implementation in Harris, Sengupta & Owens (2007).

Figures 5 and 6 show examples of the work efficient Blelloch scan, as an exclusive scan (the sum for a given item excludes the item itself). Solid lines show summation with the previous item in the array, dotted lines show replacement of the previous item with the new value. O(n) additions are performed in both the upsweep and downsweep phase resulting in the same work efficiency as the serial algorithm.

A segmented variation of scan that processes contiguous blocks of input items with different head flags can be easily formulated. This is achieved by creating a binary associative operator on key value pairs. The operator tests the equality of the keys and sums the values if they belong to the same sequence. This is discussed further in ‘Scan and reduce on multiple sequences.’

A scan may also be implemented using warp intrinsics to create fast 32 item prefix sums based on the simple scan in Fig. 4. Code for this is shown in Listing 3. Although the simple scan algorithm is not work efficient, we use this approach for small arrays of size 32.

Radix sort

Radix sorting on GPUs follows from the ability to perform parallel scans. A scan operation may be used to calculate the scatter offsets for items within a single radix digit as described in Algorithm 4 and Fig. 7. Flagging all ‘0’ digits with a one and performing an exclusive scan over these flags gives the new position of all zero digits. All ‘1’ digits must be placed after all ‘0’ digits, therefore the final positions of the ‘1’s can be calculated as the exclusive scan of the ‘1’s plus the total number of ‘0’s. The exclusive scan of ‘1’ digits does not need to be calculated as it can be inferred from the array index and the exclusive scan of ‘0’s. For example, at index 5 (using 0-based indexing), if our exclusive scan shows a sum of 3 ‘0’s, then there must be two ‘1’s because a digit can only be 0 or 1.

The basic radix sort implementation only sorts unsigned integers but this can be extended to correctly sort signed integers and floating point numbers through simple bitwise transformations. Fast implementations of GPU radix sort perform a scan over many radix bits in a single pass. Merrill & Grimshaw (2011) show a highly efficient and practical implementation of GPU radix sorting. They show speedups of 2× over a 32 core CPU and claim to have the fastest sorting implementation for any fully programmable microarchitecture.

Segmented scan

A scan can be performed on the sequences from Table 6 using the conventional scan algorithm described in ‘Parallel prefix sum (scan)’ by modifying the binary associative operator to accept key value pairs. Listing 4 shows an example of a binary associative operator that performs a segmented summation. It resets the sum when the key changes.

Segmented reduce

A segmented reduction can be implemented efficiently by applying the segmented scan described above and collecting the final value of each sequence. This is because the last element in a scan is equivalent to a reduction.

Interleaved sequences: multireduce

A reduction operation on interleaved sequences is commonly described as a multireduce operation. To perform a multireduce using the conventional tree algorithm described in ‘Reduction’ a vector of sums can be passed up the tree instead of a single value, with one sum for each unique sequence. As the number of unique sequences or ‘buckets’ increases, this algorithm becomes impractical due to limits on temporary storage (registers and shared memory).

A multireduce can alternatively be formulated as a histogram operation using atomic operations in shared memory. Atomic operations allow multiple threads to safely read/write a single piece of memory. A single vector of sums is kept in shared memory for the entire thread block. Each thread can then read an input value and increment the appropriate sum using atomic operations. When multiple threads contend for atomic read/write access on a single piece of memory they are serialised. Therefore, a histogram with only one bucket will result in the entire thread block being serialised (i.e. only one thread can operate at a time). As the number of buckets increases this contention is reduced. For this reason the histogram method will only be appropriate when the input sequences are distributed over a large number of buckets.

Interleaved sequences: multiscan

A scan operation performed on interleaved sequences is commonly described as a multiscan operation. A multiscan may be implemented, like multireduce, by passing a vector of sums as input to the binary associative operator. This increases the local storage requirements proportionally to the number of buckets.

General purpose multiscan for GPUs is discussed in Eilers (2014) with the conclusion that ‘multiscan cannot be recommended as a general building block for GPU algorithms.’ However, highly practical implementations exist that are efficient up to a limited number of interleaved buckets, where the vector of sums approach does not exceed the capacity of the device. The capacity of the device in this case refers to the amount of registers and shared memory available for each thread to store and process a vector.

Merill and Grimshaw’s optimised radix sort implementation (Merrill & NVIDIA-Labs, 2016; Merrill & Grimshaw, 2011), mentioned in ‘Radix sort,’ relies on an eight-way multiscan in order to calculate scatter addresses for up to 4 bits at a time in a single pass.

Data layout

To facilitate enumeration through all split points, the feature values should be kept in sorted order. Hence, we use the device memory layout shown in Tables 9 and 10. Each feature value is paired with the ID of the instance it belongs to as well as the leaf node it currently resides in. Data are stored in sparse column major format and instance IDs are used to map back to gradient pairs for each instance. All data are stored in arrays in device memory. The tree itself can be stored in a fixed length device array as it is strictly binary and has a maximum depth known ahead of time.

Block level parallelism

Given the above data layout notice that each feature resides in a contiguous block and may be processed independently. In order to calculate the best split for the root node of the tree, we greedily select the best split within each feature, delegating a single thread block per feature. The best splits for each feature are then written out to global memory and are reduced by a second kernel. A downside of this approach is that when the number of features is not enough to saturate the number of streaming multiprocessors—the hardware units responsible for executing a thread block—the device will not be fully utilised.

Calculating splits

In order to calculate the best split for a given feature we evaluate Eq. (4) at each possible split location. This depends on (G_L, H_L) and (G_R, H_R). We obtain (G_L, H_L) from a parallel scan of gradient pairs associated with each feature value. (G_R, H_R) can be obtained by subtracting (G_L, H_L) from the node total which we know from the parent node.

The thread block moves from left to right across a given feature, consuming ‘tiles’ of input. A tile herein refers to the set of input items able to be processed by a thread block in one iteration. Table 11 gives an example of a thread block with four threads evaluating a tile with four items. For a given tile, gradient pairs are scanned and all splits are evaluated.

Each 32 thread warp performs a reduction to find the best local split and keeps track of the current best feature value and accompanying gradient statistics in shared memory. At the end of processing the feature, another reduction is performed over all the warps’ best items to find the best split for the feature.

Missing values

The original XGBoost algorithm accounts for missing values by scanning through the input values twice as described in ‘XGBoost: handling missing values’—once in the forwards direction and once in the reverse direction. An alternative method used by our GPU algorithm is to perform a sum reduction over the entire feature before scanning. The gradient statistics for the missing values can then be calculated as the node sum statistics minus the reduction. If the sum of the gradient pairs from the missing values is known, only a single scan is then required. This method was chosen as the cost of a reduction can be significantly less than performing the second scan.

Node buckets

So far the algorithm description only explains how to find a split at the first level where all instances are bucketed into a single node. A decision tree algorithm must, by definition, separate instances into different nodes and then evaluate splits over these subsets of instances. This leaves us with two possible options for processing nodes. The first is to leave all data instances in their current location, keeping track of which node they currently reside in using an auxiliary array as shown in Table 12. When we perform a scan across all data values we keep temporary statistics for each node. We therefore scan across the array processing all instances as they are interleaved with respect to their node buckets. This is the method used by the CPU XGBoost algorithm. We also perform this method on the GPU, but only to tree depths of around 5. This interleaved algorithm is fully described in ‘Interleaved algorithm: finding a split.’

The second option is to radix sort instances by their node buckets at each level in the tree. This second option is described fully in ‘Sorting algorithm: finding a split.’ Briefly, data values are first ordered by their current node and then by feature values within their node buckets as shown in Table 13. This transforms the interleaved scan (‘multiscan’) problem described above into a segmented scan, which has constant temporary storage requirements and thus scales to arbitrary depths in a GPU implementation.

In our implementation, we use the interleaved method for trees of up to depth 5 and then switch to the sorting version of the algorithm. Avoiding the expensive radix sorting step for as long as possible can provide speed advantages, particularly when building small trees. The maximum number of leaves at depth 5 is 32. At greater depths there are insufficient shared memory resources and the exponentially increasing run-time begins to be uncompetitive.

Interleaved algorithm: finding a split

In order to correctly account for missing values a multireduce operation must be performed to obtain the sums within interleaved sequences of nodes. A multiscan is then performed over gradient pairs. Following that, unique feature values are identified and gain values calculated to identify the best split for each node. We first discuss the multireduce and multiscan operations before considering how to evaluate splits.

Sorting algorithm: finding a split

The sorting implementation of the split finding algorithm operates on feature value data grouped into node buckets. Given data sorted by node ID first and then feature values second we can perform segmented scan/reduce operations over an entire feature only needing a constant amount of temporary storage.

The segmented reduction to find gradient pair sums for each node is implemented as a segmented sum scan, storing the final element from each segment as the node sum. Another segmented scan is then performed over the input feature to get the exclusive scan of gradient pairs. After scanning each tile, the split gain is calculated using the scan and reduction as input and the best splits are stored in shared memory.

The segmented scan is formulated by performing an ordinary scan over key value pairs with a binary associative operator that resets the sum when the key changes. In this case the key is the current node bucket and the value is the gradient pair. The operator is shown in Eq. (6).

An overview of the split finding algorithm for a single thread block processing a feature is provided in Algorithm 8. The output of this algorithm, like that of the interleaved algorithm, consists of the best splits for each feature, and each node. This is reduced by a global kernel to find the best splits for each node, of any feature.