Merge Hulls Without Copying Vertices: GPU Isotonic Regression

Merge Hulls Without Copying Vertices: GPU Isotonic Regression Elena Raikova Independent researcher Buenos Aires, Argentina lenabarretta@gmail.com

arXiv:2610.24038v1 [cs.DC] 16 Oct 2026

Abstract Isotonic regression finds the monotone sequence closest to the data in weighted least squares. It is used to calibrate probabilities and inside differentiable ranking. Its standard algorithm, pool adjacent violators, is sequential within a sequence, and torchsort, the exact GPU solver of differentiable sorting libraries, runs it with one thread per sequence, which is slow on long ones. The fit is the slope sequence of the lower convex hull of cumulative sums, and hulls of sorted points can be merged in parallel. Merging partial hulls level by level, however, copies the surviving vertices at every level, which moves Θ(𝑛 log 𝑛) data when most points survive. We store each chunk’s hull once and represent a merged hull as a list of segments of these chunk hulls, at most one per chunk; a merge then updates the list. We prove the representation correct, bound the work and depth of the merge phase, and measure a CUDA implementation on four NVIDIA GPUs against a control that copies the surviving vertices with the whole GPU at every level. On dense hulls the merge phase is 17 to 276 times faster with segment lists and the whole operator 3.9 to 19 times; on noisy data, where few vertices survive, the difference is small. In a differentiable soft rank, replacing torchsort’s isotonic solver makes forward and backward 2.3 to 230 times faster on a T4, and the training step of a compact model with a Spearman loss 1.6 to 9 times.

CCS Concepts • Theory of computation → Parallel algorithms; • Computing methodologies → Parallel computing methodologies; • Mathematics of computing → Mathematical software.

Keywords isotonic regression, convex hull, GPU, CUDA, parallel algorithms, PyTorch

1

Introduction

Click models in online advertising are often trained on data in which the non-clicks were downsampled, and their predicted probabilities are then too high [12]. Before the predictions are used to price impressions, they are calibrated: a monotone curve, fitted to a day of held-out traffic, maps each score to the click rate observed at that score. Isotonic regression fits this curve [18, 29]. A day of traffic has millions of scores, and calibration per campaign or per country needs hundreds of curves. The same problem, the non-decreasing sequence closest to the data in weighted least squares, appears in differentiable sorting and ranking [5], sparse top-𝑘 operators [23] and distributional regression [13], usually as a batch of many rows. The standard algorithm, pool adjacent violators, scans a sequence once and merges neighbouring blocks while they violate the order. torchsort [17], the exact GPU solver used by differentiable sorting

libraries, runs it with one thread per sequence. This is fast for many short rows. For one row of 107 elements such a solver takes 5.5 s on a T4 and 6.4 s on an A100, while SciPy takes 0.17 s on one CPU core. Iterative methods for accelerators [23] are parallel but approximate. The fit is also the slope sequence of the lower convex hull of the cumulative sum diagram [1]. Lower hulls of points sorted by abscissa can be built in parallel: split the points, build the hulls of the parts, and merge neighbouring hulls along the bridge between them [19, 21], with logarithmic depth on a PRAM [3, 9]. On a GPU the direct version merges partial hulls level by level and copies the surviving vertices at every merge. On noisy data few vertices survive and the copies are cheap. On data that are nearly monotone almost every point is a vertex, and the merge tree moves Θ(𝑛 log 𝑛) vertices. In our representation each chunk of the input writes its hull’s vertices once, and a merged hull is a list of segments, ranges of vertices of the chunk hulls with at most one per chunk. A merge finds the bridge, shortens the two segments it touches, and appends the right list to the left one (Section 3.3). Structures that merge hulls without copying exist for pointer machines and PRAMs (Section 8); ours is a flat array with one descriptor slot per chunk. The contributions are: • Segment-referenced hull merging (Sections 3 and 4): the data structure, the merge, and a proof that segment lists represent every hull of the merge tree (Theorem 4.4). • Bounds on the work and depth of the merge phase (Section 5), for the algorithm and for the code as executed. • A GPU implementation, TorchIsotonic,1 a weighted, differentiable PyTorch operator for batches of rows of any common length, measured from rows of 16 elements to a single row of 108 (Section 6). On our tests its float64 fits agree with a float64 reference to 1.3 · 10−10 and its float32 fits to 1.2 · 10−7 , in the largest error relative to 1 + |𝑧𝑖★ | (Section 7.6). • An evaluation on T4, L4, A100 and H100 GPUs (Section 7): segment lists against two controls that copy the surviving vertices, ablations of the other design choices, comparisons with a one-thread-per-row solver and with scikit-learn, a differentiable soft rank and a training step that uses it, and the calibration of a click model on the Criteo data. 1 https://github.com/LenaBarretta/torchisotonic; the code of the experiments is in its

paper folder.

Elena Raikova

2 Background 2.1 Isotonic regression as a lower convex hull

Algorithm 1 Merging the segment lists 𝐿 and 𝑅 of two neighbouring nodes.

Given 𝑦1, . . . , 𝑦𝑛 ∈ R and weights 𝑤 1, . . . , 𝑤𝑛 > 0, the (increasing) isotonic regression problem is 𝑛 ∑︁ min 𝑤𝑖 (𝑧𝑖 − 𝑦𝑖 ) 2 . (1) 𝑧 1 ≤𝑧 2 ≤ ··· ≤𝑧𝑛 𝑖=1 ∗ Its solution 𝑧 is unique and piecewise constant: it consists of blocks

of consecutive indices on which 𝑧 ∗ equals the weighted mean of 𝑦 over the block. The decreasing problem is the increasing one for −𝑦. The classical geometric characterisation (Barlow et al. [1]; Robertson et al. [22], Theorem 1.2.1) is the one this paper computes with. Let the cumulative sum diagram be the points ∑︁ ∑︁ 𝑃𝑘 = (𝑊𝑘 , 𝑆𝑘 ), 𝑊𝑘 = 𝑤 𝑖 , 𝑆𝑘 = 𝑤𝑖 𝑦𝑖 , 𝑘 = 0, . . . , 𝑛, 𝑖 ≤𝑘

𝑧𝑖∗ = the slope of the hull edge above [𝑊𝑖 −1,𝑊𝑖 ].

(3)

Indeed, the slope of the edge from vertex 𝑃𝑎 to vertex 𝑃𝑏 equals (𝑆𝑏 − 𝑆𝑎 )/(𝑊𝑏 −𝑊𝑎 ), the weighted mean of 𝑦 over (𝑎, 𝑏], and the vertices are exactly the block boundaries. The sequential pool-adjacentviolators algorithm (PAVA) builds this hull as the monotone chain does: pooling two blocks removes a non-convex vertex. Throughout, we take a point lying on the segment between its neighbours not to be a vertex (collinear points are removed). This changes which points are called vertices but not the slopes in (3), hence not 𝑧 ∗ . The lane and chunk hulls and the walks follow it; the bridge search of the merge levels may keep a collinear point (Section 4), which changes no slope either.

Parallel hulls of presorted points

The lower hull of 𝑛 points sorted by abscissa can be built in parallel with 𝑂 (𝑛) work and polylogarithmic depth by merging the hulls of neighbouring ranges of points (Section 8). Two 𝑥-separated lower hulls merge along their bridge, the common supporting line below both. Preparata and Hong [21] find it by walking from the facing ends, Overmars and van Leeuwen [19] by a simultaneous binary search over both hulls in 𝑂 (log 𝑛) steps, and Gæde et al. [8] reduce that search to three cases on edges. We apply these merges to the cumulative sum diagram on a GPU.

3

2: if the chain is convex at 𝑠 then 3:

append 𝑅, without 𝑠, to 𝐿

⊲ no bridge

4: else

(𝑎, 𝑏) ← walk from 𝑠’s neighbours, ≤ 𝑇 moves ⊲ thread levels only 6: if there was no walk, or it gave up then ⊲ bisection 7: 𝑎 ← first vertex of 𝐿 above its tangent to 𝑅 8: 𝑏 ← where that tangent touches 𝑅 9: end if 10: end 𝐿 at 𝑎: cut its segment, drop the later ones 11: start 𝑅 at 𝑏: cut its segment, drop the earlier ones 12: move 𝑅’s descriptors to follow 𝐿’s 13: end if 5:

𝑖 ≤𝑘

(2) with 𝑃0 = (0, 0); their abscissae strictly increase since the weights are positive. The greatest convex minorant of the diagram is its lower convex hull, and

2.2

1: 𝑠 ← the last vertex of 𝐿 (the first of 𝑅)

The Algorithm

The input is a batch of rows, each an independent problem with its cumulative sum diagram (Section 2.1). The algorithm works at three levels. A lane of a warp builds the hull of 32 consecutive points. The lanes of a chunk (one to four warps, 𝐶 = 1024 to 4096 points) merge their hulls in shared memory. The chunks of a row merge in a tree, in global memory, as lists of segments. A row that fits in one chunk is finished at the second level, and rows of up to 512 elements are packed several to a warp.

3.1

Lane hulls in a bit mask

Each lane takes 32 consecutive elements, forms their cumulative sums in shared memory, and runs the monotone chain over its points. The stack of the chain is a 32-bit mask in a register, with bit 𝑘 set while point 𝑘 is on the stack. The top and the entry under it are kept in registers, and the entry under those is found with a count-leading-zeros instruction on the mask. A pop clears a bit and a push sets one. Every iteration either pops or pushes and issues both loads it might need, so the loop has no data-dependent branch and the 32 lanes of a warp stay in step. When the chain ends, the set bits are the lane’s hull.

3.2

Merging lane hulls by walking

The lane hulls of a chunk merge pairwise in rounds, five within a warp and then across warps. Neighbouring hulls share their boundary point. The first lane of a merge tests whether the chain is convex at the shared point. If it is not, the lane walks from the point’s two neighbours, moving the left end left and the right end right until the line between them supports both hulls; this line is the bridge. Every lane of the merge then clears the mask bits strictly between the bridge’s ends. Each step of a walk except the last of each pass removes a vertex, and a vertex is removed only once, so all the walks of a chunk take 𝑂 (𝐶) steps (Lemma 5.1). The chunk’s hull is then the set bits of its lanes’ masks.

3.3

Chunk hulls as segment lists

A chunk that is not a whole row writes its hull to a slot of 𝐶 + 1 entries of a global vertex array. The vertices are ranked by a prefix sum of the lanes’ mask counts and stored as coordinates relative to the chunk’s first point, with the point’s index. These entries are written once; the merges above change only descriptors. The chunk’s segment list is one descriptor ( 𝑗, ℓ, ℎ, 𝑣) of 24 bytes: the chunk 𝑗, its first and last surviving vertex ℓ and ℎ, and the slot 𝑣 that holds its vertices. A node of the merge tree spanning chunks 𝑗0, . . . , 𝑗1 is the list of the segments of its hull, at most one per chunk (Theorem 4.4). The list is stored in the descriptor slots of the node’s chunks, starting at 𝑗0 ’s (Figure 1).

Merge Hulls Without Copying Vertices: GPU Isotonic Regression

4

(a) merging by copying chunk 0

chunk 1

chunk 2

level 1 level 2 (b) merging segment lists

chunk 0 (fixed)

chunk 1 (fixed)

chunk 2 (fixed)

chunk 3 (fixed)

(0, 0, 5)

(1, 0, 4)

(2, 0, 5)

(3, 0, 3)

(0, 0, 3) (1, 2, 4) (0, 0, 3) (1, 2, 3) (3, 2, 3)

(2, 0, 4) (3, 1, 3)

level 0 level 1 level 2

Figure 1: Merging four chunk hulls. (a) Every merge copies the surviving vertices into a new array, so a vertex that survives 𝑘 levels is written 𝑘 + 1 times. (b) Each chunk writes its vertices once; a node is a list of segments (chunk, first vertex, last vertex), and a merge cuts the two segments its bridge lands in and moves the right list’s descriptors. At level 2 the bridge passes over chunk 2, which keeps no segment. Dark: vertices of the final hull. Level 𝑟 of the tree merges pairs of nodes of 2𝑟 chunks (Algorithm 1). A vertex is addressed by its segment and its position in the chunk’s slot. Its coordinates in the row are its chunk-relative ones plus the chunk’s start, which comes from a prefix sum over the chunks’ totals. The bridge is found by walking, as within the chunks, now across segment boundaries. A walk that would take more than 𝑇 = 64 steps is replaced by a search: a bisection over the candidate vertices of 𝐿, each tested against its tangent to 𝑅, both bisecting first over segments and then within one. The merge then cuts two segments and moves the descriptors of 𝑅’s tail down by at most 2𝑟 places. A level with many merges runs one thread per merge, since its lists are short and a walk is serial. A level with few merges runs one warp per merge, which searches at once with 32 candidates per step and copies descriptors in coalesced rounds.

3.4

Reconstruction

At the root, the list holds the row’s hull. One kernel finds, for every chunk, the range of its surviving vertices and the slopes of the edges into and out of it, by a bisection over the root list. A chunk without a segment lies under a single edge. A second kernel assigns each of the chunk’s elements the slope of the edge above it, by a binary search in the span’s vertices, and writes the fit.

3.5

Correctness

chunk 3

Short rows

Rows of up to 512 elements are packed several to a warp. A row takes the smallest power of two 𝐿 of lanes with 32𝐿 ≥ 𝑛, and these lanes run the lane and chunk levels on their own, with the prefix sums and merges restricted to the row. The whole row stays in shared memory and is finished in one kernel.

All statements are in exact arithmetic; rounding is the subject of Section 7.6. Let 𝑃 0, . . . , 𝑃𝑛 be the cumulative sum diagram (Section 2.1), split into chunks: chunk 𝑗 covers the points 𝑃𝑐 𝑗 , . . . , 𝑃𝑐 𝑗 +1 , so neighbouring chunks share a point. Within a chunk the same holds for its lanes. We write 𝐻 (𝑄) for the lower hull of a point set 𝑄, as its vertex sequence by increasing abscissa. Lemma 4.1 (filtering). Let 𝑄 ⊆ {𝑃0, . . . , 𝑃𝑛 }. If 𝑃𝑘 ∈ 𝑄 is a vertex of 𝐻 ({𝑃0, . . . , 𝑃𝑛 }) then 𝑃𝑘 is a vertex of 𝐻 (𝑄). Proof. 𝑃𝑘 is a vertex of the lower hull of a set exactly when some non-vertical line through 𝑃𝑘 has every point of the set on or above it, with no two points of the set on the line on both sides of 𝑃𝑘 . The line that witnesses this for the whole diagram witnesses it for the subset 𝑄. □ Hence a point that is not a vertex of its lane’s hull, or of its chunk’s, is not a vertex of the diagram’s hull, and every stage may discard the non-vertices of its own level. The end points of a lane, chunk or node are always vertices of its hull, since they are extreme in 𝑥. Lemma 4.2 (merging, [19, 21]). Let 𝐴 and 𝐵 be point sets with 𝑥 (𝑎) ≤ 𝑥𝑠 ≤ 𝑥 (𝑏) for all 𝑎 ∈ 𝐴, 𝑏 ∈ 𝐵, sharing the point 𝑠 with 𝑥 (𝑠) = 𝑥𝑠 , which is the last vertex of 𝐻 (𝐴) and the first of 𝐻 (𝐵). Then 𝐻 (𝐴 ∪ 𝐵) = 𝐻 (𝐴) [0..𝛼] ∥ 𝐻 (𝐵) [𝛽..] for a unique bridge (𝐻 (𝐴) [𝛼], 𝐻 (𝐵) [𝛽]): if the chain is convex at 𝑠 (the slope into 𝑠 along 𝐻 (𝐴) is below the slope out of 𝑠 along 𝐻 (𝐵)) then 𝛼 and 𝛽 both denote 𝑠, which is kept once; otherwise 𝑠 is not a vertex of 𝐻 (𝐴 ∪ 𝐵), 𝛼 precedes 𝑠 and 𝛽 follows it, and the vertices of 𝐻 (𝐴) after 𝛼 and of 𝐻 (𝐵) before 𝛽 lie on or above the bridge. The implementation finds the bridge by walking from 𝑠: 𝛼 moves left while the turn (𝐻 (𝐴) [𝛼 − 1], 𝐻 (𝐴) [𝛼], 𝐻 (𝐵) [𝛽]) is not strictly left, 𝛽 moves right symmetrically, alternating until neither moves. The fixed point supports both chains, so it is the bridge. Every step except the last of each pass removes a vertex, and a vertex is removed once over all merges of a node, which bounds the walks of a node by 𝑂(its vertices) in total (Section 5). Definition 4.3 (segment lists). A segment ( 𝑗, ℓ, ℎ), ℓ ≤ ℎ, denotes the vertices ℓ..ℎ of chunk 𝑗’s hull 𝐻 𝑗 . A segment list is a sequence of segments with strictly increasing chunk indices; it denotes the concatenation of its segments’ vertices. Theorem 4.4 (segment representation). For every node of the merge tree, the vertex sequence of the hull of the node’s points, without collinear middle points (Section 2.1), equals the concatenation of its segment list, and each segment is a contiguous range of one chunk hull’s vertices. The list of a node of 𝑐 chunks has at most 𝑐 segments. Proof. By induction up the tree. A leaf, chunk 𝑗, has the list [( 𝑗, 0, |𝐻 𝑗 | − 1)]. Let nodes 𝐴 and 𝐵 be neighbours with lists 𝑆𝐴 , 𝑆 𝐵 satisfying the claim; their hulls share the point 𝑠 at the end of 𝑆𝐴 and the start of 𝑆 𝐵 , and 𝐴 lies to the left of 𝑠, 𝐵 to the right. By Lemma 4.2, 𝐻 (𝐴 ∪ 𝐵) = 𝐻 (𝐴) [0..𝛼] ∥ 𝐻 (𝐵) [𝛽..]. The prefix 𝐻 (𝐴) [0..𝛼] of the concatenation of 𝑆𝐴 is the sequence of segments of 𝑆𝐴 up to the one containing 𝛼, that one cut at 𝛼’s position, ( 𝑗, ℓ, pos𝛼 ) with ℓ ≤ pos𝛼 ; the suffix 𝐻 (𝐵) [𝛽..] is the segment of 𝑆 𝐵 containing 𝛽,

Elena Raikova

cut to ( 𝑗 ′, pos𝛽 , ℎ) with pos𝛽 ≤ ℎ, followed by the rest of 𝑆 𝐵 . Cut segments remain contiguous ranges of one chunk hull, and the chunk indices of 𝑆𝐴 are below those of 𝑆 𝐵 , so the concatenation is a segment list of 𝐴 ∪ 𝐵. Since the kept part of 𝑆𝐴 and of 𝑆 𝐵 use disjoint chunks, the list has at most as many segments as the node has chunks. □ By the last clause, a node’s list fits in one slot per chunk of the node. A merge can therefore write its result where the left child’s list was, moving the right child’s tail down by at most the number of its segments. Corollary 4.5 (reconstruction). At the root, the segment list denotes 𝐻 ({𝑃0, . . . , 𝑃𝑛 }). For chunk 𝑗 with a segment ( 𝑗, ℓ, ℎ) in it, the edges above the chunk’s elements are the edges between consecutive vertices ℓ..ℎ, the edge into vertex ℓ from the previous segment’s last vertex, and the edge out of ℎ into the next segment’s first vertex; a chunk with no segment lies under a single edge, from the last vertex of the last segment of an earlier chunk to the first vertex of the next segment. Assigning each element the slope of the edge over it gives 𝑧 ∗ by (3). Degenerate cases follow the same convention. A node always has at least two vertices, its ends, so the predecessor of 𝑠 in 𝐻 (𝐴) and its successor in 𝐻 (𝐵) exist. The shared point appears in both lists; it is kept once in the convex case and dropped from both otherwise. A cut cannot produce an empty segment. A turn of zero is treated as non-convex in the lane and chunk hulls and in the walks, so collinear middle points are dropped there. Theorem 4.4 is stated for this canonical hull. The implementation’s bridge search may stop at the first of several collinear points of contact and so keep a collinear middle point in a merged list. Such a list then also contains points that lie on edges of the hull. A point on an edge changes no slope, so the fit is the same.

5

Analysis

We analyse one row of 𝑛 elements; rows of a batch are independent. Let 𝐶 be the chunk size (𝐶 = 1024 for long rows), 𝑚 = ⌈𝑛/𝐶⌉ the number of chunks, ℎ 𝑗 the number of vertices of chunk 𝑗’s hull Í and 𝐻 = 𝑗 ℎ 𝑗 ≤ 𝑛 + 𝑚, and 𝑃 = 32 the points per lane. The merge tree has 𝑅 = ⌈log2 𝑚⌉ levels; a node of level 𝑟 spans 2𝑟 chunks and, by Theorem 4.4, its list has at most 2𝑟 segments. Work counts operations over all threads; depth counts the longest chain of dependent operations.

5.1

Chunk hulls

A lane processes its 𝑃 points with one push and at most one pop each, in a loop whose only data-dependent branch is the rare fallback of the turn predicate, so the lane hulls of a chunk cost Θ(𝐶) operations at depth Θ(𝑃). The lane and warp merges inside a chunk walk to their bridges: Lemma 5.1 (walks). In a merge by walking (Section 4), every step of a pass except its last removes a vertex of one of the two children’s hulls, one that is not a vertex of the merged hull, and the outer loop ends after the first pass without a step. Hence a merge that removes 𝑑 vertices takes at most 3𝑑 +2 steps, and over all merges of a node, whose vertices are each removed at most once, the walks take 𝑂(vertices of the node’s children) steps.

Thus a chunk’s merges take 𝑂 (𝐶) steps in total. A step may also skip lanes whose masks are empty to reach the next vertex, so it costs up to 𝑂 (𝐶/𝑃) operations, and the chunk’s merges at most 𝑂 (𝐶 2 /𝑃). Their depth is that of the longest single walk, Θ(𝐶) in the worst case, when one point at the end of the chunk pulls the bridge across all of it. 𝐶 is a constant of the implementation. For fixed 𝐶 the chunk phase therefore costs Θ(𝑛) work at depth 𝑂 (𝐶), and the exclusive scan of the chunk totals Θ(𝑚) work at depth 𝑂 (log 𝑚).

5.2

Merging the chunk hulls

A merge at level 𝑟 tests the shared point, finds the bridge, and moves the right child’s tail of at most 2𝑟 descriptors down. Moving costs 𝑂 (𝑚) descriptors per level, hence 𝑂 (𝑚 log 𝑚) in total. This is the only data movement of the phase that grows with 𝑛, and it depends on the number of chunks only. For a fixed 𝐶, 𝑚 log 𝑚 = (𝑛/𝐶) log(𝑛/𝐶) is not 𝑂 (𝑛). It is below 𝑛 when 𝐶 ≥ log2 (𝑛/𝐶), which holds for 𝐶 = 1024 on every input that fits in memory; the statements below keep the term. By Lemma 5.1, finding the bridges by walking alone is 𝑂 (𝐻 ) work over the whole tree. Its depth, however, is not bounded by any polylogarithm. In a row whose last point lies far below the rest, the merges along the right spine remove all other vertices, and the walk at the top level takes Θ(𝑛) steps. The implementation therefore caps each walk and then searches: a candidate vertex of the left list is tested against its tangent to the right list, and both are found by bisection, first over segments and then within one. Theorem 5.2 (merge phase, nested search). In the work-depth model, with the descriptor moves of a merge done as one parallel copy: with every walk capped at a constant 𝑇 steps (𝑇 = 64 in the implementation) and the nested search beyond it, the merge phase over 𝑚 = ⌈𝑛/𝐶⌉ chunk hulls does 𝑂 (𝑛 + 𝑚 log 𝑚) work at depth 𝑂 (𝑇 log 𝑛 + log3 𝑛). Proof. Capped walks take at most 3𝐻 +2𝑚 = 𝑂 (𝑛) steps in total (Lemma 5.1). A merge searches at most once. At level 𝑟 there are at most ⌈𝑚/2𝑟 +1 ⌉ merges, and the outer bisection of a search takes 𝑂 (𝑟 + log 𝐶) candidates, each tested with a tangent of 𝑂 (𝑟 + log 𝐶) Í steps; the searches therefore cost 𝑟 (𝑚/2𝑟 +1 + 1) 𝑂 ((𝑟 + log 𝐶) 2 ) = 𝑂 (𝑚 log2 𝐶 + log3 𝑚), which is 𝑂 (𝑛), as 𝑚 log2 𝐶 = 𝑛 log2 𝐶/𝐶. The descriptor moves cost 𝑂 (𝑚 log 𝑚). A level’s depth is at most 𝑇 + 𝑂 ((𝑟 + log 𝐶) 2 ), and the sum over the ⌈log2 𝑚⌉ levels is 𝑂 (𝑇 log 𝑚 + log3 𝑛). □ Theorem 5.2 assumes that a merge moves its descriptors in one parallel step. In the code, a merge run by one thread moves them one by one, and a merge run by a warp moves them in rounds of 32. This changes the depth. Proposition 5.3 (merge phase, as executed). The CUDA code’s merge phase does 𝑂 (𝑛+𝑚 log 𝑚) work at depth 𝑂 (𝑇 log 𝑛+log3 𝑛+𝑚): the descriptor moves add at most 𝑚/16 + log2 𝑚 sequential rounds to the bound of Theorem 5.2. Proof. A level runs one thread per merge only if it has at least 256 merges, so its nodes span 2𝑟 ≤ 𝑚/512 chunks, and a thread moves at most 2𝑟 descriptors; over these levels that is at most 𝑚/256 steps. A level with fewer merges runs a warp per merge, which moves at most 2𝑟 descriptors in ⌈2𝑟 /32⌉ rounds; over all

Merge Hulls Without Copying Vertices: GPU Isotonic Regression

Í levels at most 𝑟 <⌈log2 𝑚⌉ (2𝑟 /32 + 1) < 𝑚/16 + log2 𝑚 rounds. The warp levels do not walk, and their search is the nested one with 32 candidates per step, so their depth is at most that of Theorem 5.2. □

size is the largest whose shared memory (the chunk’s cumulative sums, lane starts and vertex stack) fits in the 48 KB of a launch. The lane hulls live in registers and use no local memory. A chunk’s vertex slot takes 12 bytes per vertex, 20 with weights.

The added term is linear in 𝑛 with a small constant, 𝑚/16 = 𝑛/16 384. For 𝑛 = 108 , 𝑚 = 97 657, and the moves take at most about 6 100 rounds. The term is largest on dense hulls, where every chunk keeps a segment. The measurements of Section 7 include it. The cube in the depth comes from the nested search, a bisection inside a bisection. The simultaneous bisection of both chains of Overmars and van Leeuwen [19], in the three-case form of Gæde et al. [8], finds a bridge in a logarithmic number of steps and lowers the depth of the model to 𝑂 (𝑇 log 𝑛+log2 𝑛) (Appendix A, Theorem A.2). We implemented it as a variant. It is slower in practice (Section 7.4), so the default keeps the nested search. Since the levels run one after another, neither search reaches the 𝑂 (log 𝑛) depth of Goodrich’s CREW hull [9].

Merge levels. The ⌈log2 𝑚⌉ levels are separate launches. A level runs one thread per merge if it has at least 256 merges, and one warp per merge otherwise. A thread walks for at most 64 steps before it switches to the nested search; a warp searches directly, 32 candidates per step, and moves the tail of the right list with coalesced copies. The three-case search of Theorem A.2 is implemented as a variant for the thread-per-merge levels.

5.3

Reconstruction

The ends kernel locates each chunk’s segments by bisection in the root list, which is 𝑂 (𝑚 log 𝑚) work. The fit kernel loads at most ℎ 𝑗 + 2 slopes per 1024 elements and finds each element’s edge by a binary search of 𝑂 (log ℎ 𝑗 ) steps, ℎ 𝑗 ≤ 𝐶 + 1. The phase costs 𝑂 (𝑛 log 𝐶) work in the worst case and Θ(𝑛) for small hulls, at depth 𝑂 (log 𝑛).

5.4

Totals, data movement and memory

Over all phases the work is 𝑂 (𝑛 log 𝐶 + 𝑚 log 𝑚), the searches included. The fit contributes a factor log 𝐶 per element and the descriptor moves log 𝑚/𝐶. The depth is 𝑂 (𝐶 + log3 𝑛) in the model and 𝑂 (𝐶 + log3 𝑛 + 𝑚) as executed (Proposition 5.3). The walks and searches depend on the hull’s geometry in either representation. The difference lies in the data moved. With materialized hulls a merge at level 𝑟 copies the right child’s surviving vertices, up to all of them, so the phase moves Θ(𝐻 log 𝑚) vertices of 12 to 20 bytes, which is Θ(𝑛 log 𝑛) when all points are vertices. With segment lists it moves 𝑂 (𝑚 log 𝑚) descriptors of 24 bytes, 40 MB for 𝑛 = 108 . Peak memory is dominated by the vertex slots, 𝐶 + 1 per chunk of 12 bytes (20 with weights). This is three times the size of a float32 input, or 2.5 times that of the values and weights together. The lists, counts and chunk starts take 𝑂 (𝑚).

6

Implementation

The operator has one CUDA code path for float32 and float64 and for weighted and unweighted problems. A call launches at most five kernels plus one per merge level. It also fills the rows’ NaN flags and, for weighted calls, checks that the weights are positive, which waits for the GPU. The constants below were tuned on the T4 (Section 7). Blocks and memory. Rows of up to 512 elements are packed several per warp and solved in registers and shared memory. Rows of up to 4 096 elements are one chunk, handled by a block of one, two or four warps; for weighted float64 the limit is 2 048, as a block of four warps would not fit in shared memory. Longer rows are split into chunks of 1 024 elements, one warp per chunk. The block

Coordinates. The coordinates form a hierarchy, so that differences of nearby points keep their precision. Within a lane, cumulative sums are taken in the working type from the lane’s first point. In float32 they are pairs (hi, lo) from a renormalized two-sum. The lo part is stored in shared memory as a 16-bit code in units of 2−38 of the power of two of hi, which suffices since a renormalized pair bounds it by 214 such units; a pair then takes 6 bytes instead of 8 and resolves 2−38 of a sum’s magnitude. Lane starts within a chunk are pairs from a two-sum scan across the warp. A chunk’s vertices are stored as doubles relative to the chunk, and the chunk starts along a row are double-double pairs from a scan over the chunk totals. A slope between vertices of different chunks is formed from the differences of the hi parts (exact for nearby chunks, by Sterbenz’s lemma), of the lo parts and of the chunk-relative coordinates. Absolute coordinates reach 1010 at 108 elements, and a difference of two of them would lose six digits. Predicates. The turns of the lane and chunk levels are first evaluated on the hi parts. If the result exceeds an a-priori error bound it is used; otherwise the turn is recomputed from the pairs with two-sum and two-product. The bound is in terms of the magnitudes of the summed terms, the lane’s start and the sum within the lane, rather than of the coordinate alone, since the two may cancel. The merge levels evaluate the cross product from absolute coordinates, one addition per point, and recompute it from the differences above only if its magnitude is within an error bound. That recomputation is double arithmetic at the scale of a chunk’s sums. On the T4 the merge levels are limited by double throughput, a thirty-second of the float rate. The cross product takes seven double operations, and evaluating its bound in double, or in float with conversions that run at the double rate, doubles the time of the merge levels. The code therefore replaces each term of the bound by a power of two read from the exponent field of a double’s high word, including the chunk starts, and compares exponents in integer arithmetic. The test is coarser than the bound by a constant factor, and zero or subnormal operands always take the second evaluation. Filters of this kind are classical [2, 24]; this one uses no floating-point operation. Appendix B gives the bounds and their derivation, and Section 7.6 the accuracy they lead to. Output. Every slope of the output is rounded to the working type once: in float32 from a quotient of the pairs corrected by its remainder, in float64 from double coordinates. The output is not guaranteed to be monotone (Section 7.6). The option exact=False drops the pairs and the filter of the lane and chunk levels for float32.

Elena Raikova

Gradient. The fit is constant on blocks, each the weighted mean of 𝑦 over the block, so its Jacobian is block averaging. With 𝑔 the gradient with respect to the fit and 𝐺, 𝑊 the sums of 𝑔 and 𝑤 over an element’s block, 𝜕/𝜕𝑦𝑖 = 𝑤𝑖 𝐺/𝑊 and 𝜕/𝜕𝑤𝑖 = (𝑦𝑖 − 𝑧𝑖∗ ) 𝐺/𝑊 , wherever the blocks do not change. The blocks are those of the hull. When a gradient is needed, the forward pass also writes one byte per element, set where the element’s point is a vertex of the row’s hull, that is where a block starts. The fit kernel knows this for every element it writes. An element at the start of a chunk is a vertex when it closes the previous chunk’s segment, which the ends kernel records. The backward pass sums over these blocks without atomics. A warp takes a tile of the row, sized to the row and at most 1 024 elements. A row that is one tile is finished in one kernel; a longer row takes three steps, and the sums of the blocks that cross tiles come from prefix sums over the tiles’ boundary sums. The sums are pairs of the working type. Neighbouring blocks whose exact means differ by less than the solver resolves may be pooled, so near such ties the gradient follows the solver’s partition; Appendix C describes how the gradients are validated.

Figure 2: One row of 𝑛 noisy float32 elements on a T4, median of ten calls: the default, the option without the float32 pairs, and torchsort (to 223 ).

7 Evaluation 7.1 Setup

from 106 rows of 16 elements through 104 × 103 , 103 × 104 and 102 × 105 to one row of 107 and one of 108 .

The measurements are on four NVIDIA GPUs of four generations: a Tesla T4 (Turing) hosted by Kaggle, on a host with an Intel Xeon at 2.00 GHz (2 cores, 4 threads) and 31 GB of memory, driver 580.178 and Python 3.13; and an L4 (Ada), an A100-SXM4-40GB (Ampere) and an H100 80GB HBM3 (Hopper) hosted by Modal, driver 610.57 and Python 3.12. All run CUDA 12.8 and PyTorch 2.11. The baselines are torchsort 0.1.10, built from PyPI against the installed PyTorch, and, on the T4’s host, SciPy 1.16.3 and scikit-learn 1.6.1 with NumPy 2.1.3. We measure TorchIsotonic 0.2.1. The constants of the implementation (chunk sizes, the walk cap, the scheduling threshold) were tuned on the T4 alone; the other three GPUs show whether the design depends on it. Inputs are resident on the GPU. A fit is timed with the host’s clock around the call and a torch.cuda.synchronize after it, over ten calls after a warm-up and a busy wait that raises the GPU’s clock, stopping early once the calls have taken ten seconds, so that the slowest torchsort cells are of two calls. We report the median and, where it matters, the interquartile range. The gradient benchmark times each call with CUDA events, and the per-stage times are GPU kernel times from the PyTorch profiler, means over five calls. On noisy data the per-stage times repeat to within a few percent across sessions, while on dense data the chunk stage of 108 increasing elements varies by up to 25% between calls of one session on the T4. Comparisons between variants are therefore taken within one session, with the default profiled beside each variant; the differences between representations in Section 7.3 are far larger than this variation, while some of the ablations are not. The data families determine how many points survive as hull vertices. In noise, a ramp from 0 to 3 plus standard normal noise, few do, and in strictly increasing data all points are vertices. In a sawtooth the bridges reach across many chunks, and in a cliff, increasing values followed by one far below, they reach across the whole row. Each number below states its family. The shapes run

one thread per row 103 102

ms

exact (default) exact=False

101 100

214

217

220

223

226

elements in one row

7.2

Comparison with a one-thread-per-row solver

Table 1 and Figure 2 compare the operator with torchsort [17], the exact GPU solver used by differentiable sorting libraries, on unweighted float32 rows. On many short rows, where a thread per row keeps every thread busy, torchsort is 7 to 43 times slower. The ratio grows with the row’s length and reaches three to four orders of magnitude for a row of 107 elements. torchsort’s time for such a row is about the same on the H100 as on the T4 and longer on the A100, as a single thread does not benefit from a larger GPU. With the gradient included the ratio is smaller: 4.2 against 6 917 ms for a row of 107 elements on the T4 (Section 7.8). On the T4’s host CPU (medians), SciPy’s pool adjacent violators on one core takes 171 ms for a row of 107 elements, and the library’s CPU path, parallel over rows on the host’s four threads, 380 ms for 106 rows of 16; the GPU operator is 80 and 169 times faster.

7.3

Segment-referenced versus materialized merging

Against torchsort the operator is faster on every shape (Section 7.2). This section isolates the representation, whose gain depends on how many vertices survive the merges. Table 2 and Figure 3 compare it with two controls. Both keep the pipeline and, after every merge, write the merged hull’s vertices out as one contiguous array. In permerge copying, the thread or warp that found a bridge also copies the right child’s surviving vertices down, in place, as the textbook merge does. At the top levels of a long row a single warp then copies tens of millions of vertices, so this control measures copying at the speed of one warp. Grid copying finds every merge’s bridge in one kernel and writes every merged hull into another buffer in a second one, with all of the GPU’s blocks and coalesced accesses, as stream compaction would [4]. Each surviving vertex then moves

Merge Hulls Without Copying Vertices: GPU Isotonic Regression

Table 1: Decreasing unweighted float32 fits, data on the GPU, ms, median of ten calls (two for the slowest torchsort cells). The interquartile range of ours is within 6% of the median except at 100 × 103 , where a call takes 0.1 ms (up to 38%, on the T4). torchsort runs pool adjacent violators with one thread per row; it was not run on 108 elements.

batch × 𝑛

T4 (Turing) ours torchsort

L4 (Ada) ours torchsort

A100 (Ampere) ours torchsort

H100 (Hopper) ours torchsort

100 × 103 106 × 16 104 × 103 103 × 104 1 × 106 1 × 107 1 × 108

0.12 2.25 2.30 1.77 0.54 2.13 20.0

0.11 0.73 0.65 0.88 0.40 1.23 8.38

0.10 0.49 0.53 0.58 0.29 0.83 4.33

0.07 0.28 0.35 0.37 0.21 0.56 2.74

3.04 72.2 99.8 197 550 5 493 -

2.42 6.95 14.0 139 415 4 947 -

3.09 3.40 18.7 190 588 6 377 -

2.42 1.91 12.5 123 432 4 683 -

Table 2: The merge stage alone, ms, one row of 𝑛 float32 elements: segment lists (seg.) against two controls that copy the surviving vertices of every merge into a flat array, with the whole GPU at every level (grid) or by the thread or warp of each merge (merge). Same chunk kernel, bridge search and schedule.

data, 𝑛

seg.

T4 grid merge

seg.

L4 grid merge

increasing, 107 increasing, 108 cliff, 108 noise, 107 noise, 108

0.19 1.85 2.21 0.50 2.21

25.1 351 343 0.49 3.76

0.15 1.17 1.32 0.40 0.80

21.6 323 320 0.38 1.55

increasing

171 2 350 2 312 0.39 3.92

noise

merge stage, ms

per-merge copy 103

grid copy

102 101

segments

100 10−1 10−2 214

220

226

elements

214

220

226

elements

Figure 3: The merge stage of one row of 𝑛 float32 elements on a T4: segment lists against copying the surviving vertices with the whole GPU (grid copy) and merge by merge, on increasing data (every point a vertex) and on noise. once per level at the memory’s bandwidth. Both controls evaluate their turns on absolute coordinates without the filter of Section 6; with plain predicates the default’s merge stage is up to 24% faster (Section 7.4). Where most points survive (increasing data, or a cliff whose bridges reach across the row), both controls move Θ(𝑛 log 𝑚) vertices. Grid copying moves them at half to two thirds of each GPU’s peak bandwidth. The merge stage with segment lists is 133 to 190 times faster on the T4, 141 to 276 times on the L4, 23 to 41 times on

seg.

162 0.19 1 644 1.25 1 611 1.32 0.32 0.33 0.99 0.60

A100 grid merge 4.34 50.7 50.2 0.32 0.96

seg.

223 0.14 2 175 0.95 2 128 1.01 0.27 0.23 0.56 0.40

H100 grid merge 2.40 27.5 27.3 0.21 0.60

173 1 717 1 676 0.17 0.32

the A100 and 17 to 29 times on the H100, at 107 and 108 elements. Per-merge copying is another 5 to 72 times slower. Figure 3 shows the growth on the T4: the copies grow as 𝑛 log 𝑛, the segment lists with the number of chunks. On noisy data only a few vertices of each chunk’s hull survive its first merges. There the three are close at 107 elements, with per-merge copying up to 27% faster, as a flat array saves one indirection per vertex. At 108 the segment lists are 1.5 to 2.0 times faster than grid copying, and between 1.8 times faster (T4) and 23% slower (H100) than per-merge copying. Table 3 gives the whole operator. Only the chunk stage is the same kernel in both. The control reconstructs the fit from its flat array with a kernel of its own, and before its first merge it tags every vertex slot with the vertex’s chunk, since a copied vertex leaves its chunk’s slot. On dense hulls the merge stage accounts for most of the difference: the operator with segment lists is 9.1 to 9.3 times faster on the T4, 12 to 19 times on the L4, 4.3 to 5.9 times on the A100 and 3.9 to 5.3 times on the H100. On noise at 107 elements the merge stages take the same time, and the operator’s difference of 10 to 13% comes from the control’s other kernels: on the T4 its reconstruction takes 0.23 against 0.19 ms and its tagging and buffers 0.28 against 0.10 ms. Part of this is the cost of a flat layout and part that of the control’s implementation, so we do not count it for the segment lists. At 108 noisy elements the merge stage contributes 1.6 of the T4’s 5.1 ms difference.

7.4

Ablations

Each variant changes one choice of the default. Table 4 gives its time relative to the default’s on rows of 108 elements, for the chunk stage in the first variant and the merge stage in the others. The

Elena Raikova

Table 3: The whole operator against grid copying, ms of GPU kernel time (all kernels), one row of 𝑛 float32 elements; same profiles as Table 2.

data, 𝑛

seg.

T4 grid

seg.

L4 grid

A100 seg. grid

H100 seg. grid

increasing, 107 increasing, 108 cliff, 108 noise, 107 noise, 108

3.25 42.8 41.6 2.01 18.2

29.4 396 386 2.22 23.3

2.01 19.4 17.9 1.10 7.99

23.8 345 340 1.21 10.4

1.33 11.9 10.4 0.75 4.25

0.83 7.28 6.39 0.49 2.67

Table 4: Ablations: time of a variant over the default’s, one row of 108 float32 elements; > 1 means the default is faster. variant (stage)

data

T4

L4 A100 H100

stack in local memory instead of a bit mask (chunk) three-case search, cap growing with the level (merge) turns on absolute coordinates (merge) a warp for every merge (merge) warps for capped walks only (merge)

noise increasing noise cliff noise sawtooth noise sawtooth noise sawtooth

2.10 3.37 1.89 1.05 1.02 0.84 2.37 0.72 1.15 0.71

1.58 1.25 2.89 1.08 0.93 0.76 3.43 0.57 0.95 0.52

1.77 1.05 1.91 1.03 0.90 0.80 1.37 0.39 1.10 0.60

1.81 1.17 1.99 1.03 0.88 0.82 1.29 0.36 1.06 0.55

bit-mask lane hulls are faster than a stack in local memory with a pop loop on every GPU, by 1.6 to 2.1 times on noise. On increasing data, where the stack is deepest, the difference is 3.4 times on the T4, which has the least local-memory bandwidth, and 1.05 to 1.25 times on the others. The three-case search of Theorem A.2, with a walk cap that grows with the level, is two to three times slower on noisy data and about as fast on cliffs, where bridges are long. Plain double predicates on absolute coordinates, without the coordinate hierarchy and its filter, are up to 24% faster on sawtooth data, where merges walk far, and up to 12% on noise; the accuracy of Section 7.6 depends on the hierarchy. The last two variants change the schedule. A warp for every merge makes the merge levels 1.4 to 2.8 times faster on sawtooth data, where many walks reach their cap and a single thread then searches, and 1.3 to 3.4 times slower on noise, where many short merges leave most of a warp idle. In the last variant a thread walks first and a second kernel finishes the merges whose walk reached its cap, with a warp each. It makes sawtooth merges 1.4 to 1.9 times faster on all four GPUs, but is up to 15% slower on noise on three of them, which is why the default keeps the single rule.

7.5

Performance breakdown

On one row of 108 float32 elements the fit takes 20.0, 8.4, 4.3 and 2.7 ms on the T4, L4, A100 and H100 (medians). A device-to-device copy of the input takes 3.4, 3.5, 0.6 and 0.27 ms, so the fit takes 5.9, 2.4, 7 and 10 times as long as the copy. On noise the chunk stage takes 63 to 73% of the kernel time on the four GPUs, the merge levels 10 to 15%, the reconstruction 7 to 19% and the prefix sums the rest. The chunk stage is limited by the instructions of the lane

5.68 63.6 61.5 0.84 5.88

3.18 34.9 33.8 0.54 3.65

hulls and their merges. The fit is therefore furthest from the copy on GPUs with the most bandwidth per unit of compute, and closest on the L4, whose memory is slowest relative to its arithmetic. With every point a vertex (increasing data) the same row takes 42.8, 19.4, 11.9 and 7.3 ms of kernel time, 2.4 to 2.8 times as long as noise, since the chunk stage writes every vertex and the fit reads them all. On the T4 the reconstruction’s share grows from 10 to 42%. On shorter rows the shares change. At 107 noisy elements the merge levels take 25, 37, 45 and 48% of the kernel time on the four GPUs, and on the H100 as long as the chunk stage (0.23 against 0.20 ms). Each of the 14 levels costs about 0.016 ms however few merges it has, which is the latency of a launch.

7.6

Numerical accuracy

The reference is a float64 solution on the CPU of the same cast input, so that the comparison measures the solver and not the cast; it is the library’s own float64 pool adjacent violators, which the test suite checks against SciPy. Ten data families (increasing, decreasing, constant, piecewise constant with noise, Gaussian, random walk, sawtooth, a cliff, nearly collinear points, and data offset by 106 ) are crossed with 𝑛 ∈ {103, 105, 107 }, both precisions, and three weightings: none, uniform in [0.1, 1.1], and 𝑒𝑈 [0,8] , a dynamic range of 3 · 103 . We report four quantities: • the largest step down of the output, relative to max |𝑦|; • the largest deviation from the reference, max𝑖 |𝑧𝑖 −𝑧𝑖★ |/(1 + |𝑧𝑖★ |); • the excess of the weighted least-squares objective over the Í reference’s, relative to the reference’s, or to 𝜀 𝑖 𝑤𝑖 𝑦𝑖2 where that is near zero; • the largest violation of the KKT conditions, computed in Í float64 from the output alone, relative to 2 𝑖 𝑤𝑖 |𝑦𝑖 |. The Í multipliers 𝜆𝑖 = −2 𝑗 ≤𝑖 𝑤 𝑗 (𝑧 𝑗 − 𝑦 𝑗 ) must be non-negative, zero at the end, and zero where the output steps up. This measure shrinks with 𝑛 for a fixed error per element, which is why the deviation is reported too. We also compare the number of blocks with the reference’s. In float64 no output steps down on these families. The deviation is at most 1.3 · 10−10 (extreme weights, 𝑛 = 107 ), the objective is within 7·10−10 , the KKT residual is at most 8·10−14 , and the number of blocks is the same everywhere. This includes 107 increasing points with extreme weights, whose cumulative sums reach 1010 . A version that forms slopes from differences of absolute sums, without the coordinate hierarchy, pools 681 of these points that are vertices of the reference and steps down by up to 9 · 10−7 .

Merge Hulls Without Copying Vertices: GPU Isotonic Regression

In float32, over all families, sizes and weightings, the deviation is at most 1.2 · 10−7 , the objective within 7 · 10−9 and the KKT residual below 6·10−8 . No output steps down, here or in the test suite, which checks monotonicity on its own families with three weightings. On 107 consecutive floats, each one ulp above the previous, every point stays its own block without weights, as in the reference. With weights, about 8% of neighbouring outputs coincide: the products 𝑤𝑖 𝑦𝑖 are rounded to float32 before the hull is built, and neighbouring slopes may then round to the same value. A further family stresses cancellation. Each 64 elements hold 32 values of +𝐵, one of −32𝐵 and 31 small, nearly collinear values, for 𝐵 from 103 to 1012 and 𝑛 from 500 to 105 . The cumulative sums rise and fall back within every chunk, so a coordinate is the difference of terms far larger than itself. In float32 the output does not step down and deviates from the reference by at most 4 · 10−9 of max |𝑦|. In float64 it deviates by at most 1.4·10−15 of max |𝑦|, and neighbouring outputs step down by up to 5 · 10−16 of max |𝑦|, a few units in the last place of the data’s scale. Monotone output is therefore not guaranteed in float64 either. A cumulative maximum over the output would enforce it at the cost of one more pass over the data; the operator does not apply one. These results use the default settings. The option exact=False drops the float32 pairs within the lanes and the filter of the chunk level; float64 is unaffected, as its sums use no pairs. On a T4 it takes 1.30 against 2.30 ms for 104 rows of 103 elements and 13.7 against 20.0 ms for one row of 108 (medians). Without the pairs, a lane’s cumulative sum of 32 terms can be 32 times a single term, or 32𝑒 8 times with extreme weights. Its float32 rounding then loses the last bits of each increment and decides which points are vertices. On the 107 unweighted consecutive floats, exact=False keeps 1.17 million of the 107 blocks. It pools 88% of the points and so has less to merge and fit, which is part of its speed. With extreme weights it deviates from the reference by up to 4.4 · 10−3 and steps down by up to 8.5 · 10−4 of the data’s scale.

7.7

Case study: calibrating a click model

Isotonic regression is a standard non-parametric calibration of classifier scores [18, 29]. We reproduce a common industrial setting on the Criteo display-advertising data [7], 45.8 million impressions over seven days in time order. A logistic regression with hashed features is trained on the first five days with a tenth of the non-clicks, a negative downsampling that leaves the predicted click probabilities far too high [12]. The sixth day calibrates the model and the seventh evaluates it, with one curve for all traffic and, separately, one curve for each of the 280 frequent values of a categorical feature. Isotonic calibration brings the expected calibration error from 0.41 to 0.0015; Appendix D compares it with the analytic correction and Platt scaling. The fits agree with scikit-learn’s to 6 · 10−16 for one curve and to 7 · 10−10 per segment, where the keys (segment plus score) are float64. We time the fits alone on both sides. On a T4, one curve over the 6.5 million scores of a day takes 0.042 s, and scikit-learn takes 2.0 s on one core of the host CPU. Over all 45.8 million scores the times are 0.155 and 17.2 s, including the sort of the scores and the pooling of equal ones. With the scores in host memory, as scikitlearn’s are, the copy to the GPU doubles the time to 0.31 s. The 280

Table 5: Fit and fit plus gradient on the T4, ms, decreasing unweighted float32; median of ten calls timed with CUDA events, and the interquartile range of the fit plus gradient.

batch × 𝑛 106 × 16 104 × 103 1 × 107

fit

ours +grad (IQR)

2.19 5.52 (4.77–5.69) 2.30 2.90 (2.88–2.90) 2.12 4.19 (4.18–4.22)

torchsort fit +grad 74.0 100 5 493

84.7 107 6 917

per-segment curves take 0.033 against 0.77 s. They are fitted as one problem: the scores are keyed by segment and the targets of each segment offset by twice its index, so that no block can extend across a segment boundary. The operator takes a rectangular batch, so segments of unequal length are either padded, exactly, with values above the row’s maximum, or laid end to end in one row as here; the merge tree then spreads that row across the whole GPU.

7.8

Gradients

Table 5 gives the fit alone and with its gradient with respect to 𝑦, for a random 𝑔, on the T4. The backward pass costs a fraction of the forward pass on rows of 103 elements and about as much on long rows. On 106 rows of 16 it costs more than the fit, as it pays a fixed cost per tile while the forward pass packs several rows to a warp. The lead over torchsort is therefore smaller with the gradient than for the fit alone.

7.9

Differentiable ranking and a training step

The soft rank of Blondel et al. [5] with 𝐿2 regularization of strength 𝜀 projects 𝜃 = 𝑥/𝜀 onto the permutahedron: it sorts 𝜃 in decreasing order into 𝑠, solves a decreasing isotonic regression of 𝑧 = 𝑠 −𝑤 with 𝑤 = (𝑛, . . . , 1), and returns 𝑠 minus that fit in the original order. We time torchsort’s soft_rank against a copy of its implementation in which only the isotonic solver, forward and backward, is replaced by TorchIsotonic; the sort, the inverse permutation and the backward pass around the solver are torchsort’s code. The values are standard normal and the GPU is the T4 (Table 6). We use two strengths: 𝜀 = 1, torchsort’s default, at which the ranks of these values are pooled into one block per row, and 𝜀 = 1/𝑛, at which about 70% of the elements start a block and the soft ranks are close to the hard ones. With 𝜀 = 1/𝑛, forward and backward together are 2.4 to 231 times faster with TorchIsotonic, the gain growing with the length of the rows; with 𝜀 = 1 they are 2.3 to 141 times faster. The training step is a synthetic benchmark of a compact scoring model: a twolayer network (32 → 128 → 1) scores each item, the soft rank of the scores is taken per query, and the loss is one minus the Spearman correlation with the target ranks. Its gain depends on the soft rank’s share of the step. For 4 queries of 105 items with 𝜀 = 1/𝑛, the soft rank takes 93% of torchsort’s step and 12% of ours, and the step is 9.1 times faster. For 256 queries of 103 items it takes half of torchsort’s step, and the step is 1.8 times faster. With 𝜀 = 1 the steps are 1.6 to 4.4 times faster. The model is cheap, so the soft rank’s share is large; with a heavier model the gain shrinks with that share. We checked both solvers against a float64 solution of the same float32 problem 𝑧. The values and gradients of TorchIsotonic agree

Elena Raikova

Table 6: Soft rank (forward and backward) and a training step with a Spearman loss on the T4, ms, median of 20 calls, 𝜀 = 1/𝑛. The two differ only in the isotonic solver. The training step is a forward, backward and Adam step of a two-layer scoring model on 𝑛 items per query.

batch × 𝑛 1024 × 100 256 × 103 32 × 104 64 × 104 4 × 105 8 × 105 1 × 106

soft rank torchsort ours 2.46 6.78

1.04 1.00

41.9

1.75

208 877

2.25 3.80

training step torchsort ours

14.8 46.3

8.38 10.9

127

13.9

with it closely, except at one boundary between two blocks of exactly equal means, where the fit is not differentiable. torchsort’s values deviate by up to 2.5 · 10−4 with 𝜀 = 1, where it accumulates long blocks in float32, and its gradients by up to 4% in the ℓ2 norm with 𝜀 = 1/𝑛 and 𝑛 ≥ 104 . Appendix C gives the measures and their analysis.

8

Related Work

Isotonic regression. The equivalence of isotonic regression with the greatest convex minorant of the cumulative sum diagram is classical [1, 22]. Stout [27] surveys the sequential algorithms, from pool adjacent violators to the 𝑂 (𝑛) bounds for the 𝐿2 case. Kearsley, Tapia and Trosset [15] and Kearsley [14] parallelize by partitioning the row, solving the parts and merging their solutions, the scheme behind Spark’s implementation. The merge of the parts’ solutions is sequential there; segment lists are a way to do this merge on a GPU. Henzi, Mösching and Dümbgen [13] accelerate the algorithm for isotonic distributional regression, where one data set yields a batch of weighted problems, the setting of our operator. On accelerators, Sander et al. [23] solve isotonic regression approximately by Dykstra’s projections inside a differentiable top-𝑘. torchsort [17], built for the soft sorting and ranking of Blondel et al. [5], runs the exact pool adjacent violators with one thread per row, forward and backward, and is our main baseline. Parallel hulls of sorted points. The merge of two 𝑥-separated hulls along a bridge is due to Preparata and Hong [21], and its 𝑂 (log 𝑛) search to Overmars and van Leeuwen [19]. Their dynamic hull keeps chains in concatenable queues, balanced trees that are split at a bridge and concatenated in 𝑂 (log 𝑛), so a merge relinks nodes instead of copying vertices. Segment lists apply the same idea in a flat form: at most one descriptor per chunk, cut at the bridge, concatenated by moving descriptors, and read by indexing. Gæde et al. [8] reduced the case analysis of Overmars and van Leeuwen to three cases; we use their form, mirrored to lower hulls, in Theorem A.2. For presorted points, Goodrich [9] reaches 𝑂 (log 𝑛) depth on a CREW PRAM. Berkman, Schieber and Vishkin [3] reach 𝑂 (log log 𝑛) on a CRCW PRAM with an implicit representation of the hull, since writing it into a contiguous array, a compaction, would take Ω(log 𝑛/log log 𝑛) time. Our layout, one slot per chunk,

serves the same purpose on a GPU. Browne et al. [6] give an 𝑂 (𝑛)work, 𝑂 (log 𝑛)-span algorithm in the binary forking model. Hulls on GPUs. Work on two-dimensional hulls on GPUs treats unsorted point clouds, with QuickHull and parallel partitioning [26, 28] or by filtering most points with extreme witnesses before building the hull of the rest, as Keith, Ferrada and Navarro do [16]. Presorted input needs no such filter. Closest in design is the warpcooperative three-dimensional hull in CuACD [25]. A warp there partitions its points among lane groups, merges the sub-hulls pairwise in logarithmically many rounds, and finds bridges by giftwrapping with warp shuffles. We also use the warp as the unit of design and merge in pairwise rounds. The differences are the input (sorted, two-dimensional, up to 108 points in one problem), the per-lane bit-mask stack, the walk whose cost is amortized against removed vertices, the hierarchy with segment references, and floating-point predicates with error bounds in place of integer determinants. Stream compaction [4] is the usual GPU method to keep the survivors of a step; the grid-copying control of Section 7.3 uses it. Choosing a thread, a warp or more per unit of work by its size is common, as in the CSR-Adaptive sparse matrix-vector product [11], and our per-level rule is of this kind. Our hi–lo pairs use float-float arithmetic on GPUs [10]. Floating-point filters for robust predicates go back to Shewchuk [24], and Bartels, Fisikopoulos and Weiser [2] derive static and semi-static filters. Our filter compares its bound by exponents only and falls back to double arithmetic on local coordinates (Appendix B).

9

Conclusion

The fit of isotonic regression can be computed in parallel within a row on a GPU by merging lower hulls. Segment lists lower the cost of data movement in that merging: they remove the copying of surviving vertices, which dominates on dense hulls and matters little on sparse ones. The resulting operator is fast enough to change the cost of a training step that uses a differentiable soft rank. Limitations. The operator takes a rectangular batch; rows of unequal length are padded or laid end to end. The merge levels follow a fixed schedule, which is far from the best on some data, and on mid-sized rows they cost more in launches than in work. The executed code adds sequential rounds of descriptor copying to the depth of the model (Proposition 5.3). The predicates are filtered and their second evaluation is double arithmetic, so accuracy and monotone output are measured properties, and near ties the gradient follows the solver’s partition.

References [1] R. E. Barlow, D. J. Bartholomew, J. M. Bremner, and H. D. Brunk. 1972. Statistical Inference under Order Restrictions. Wiley, London and New York. [2] Tinko Bartels, Vissarion Fisikopoulos, and Martin Weiser. 2023. Fast FloatingPoint Filters for Robust Predicates. BIT Numerical Mathematics 63, 2 (2023), 31. arXiv:2208.00497. doi:10.1007/s10543-023-00975-x [3] Omer Berkman, Baruch Schieber, and Uzi Vishkin. 1996. A Fast Parallel Algorithm for Finding the Convex Hull of a Sorted Point Set. International Journal of Computational Geometry & Applications 6, 2 (1996), 231–241. doi:10.1142/ S0218195996000162 [4] Markus Billeter, Ola Olsson, and Ulf Assarsson. 2009. Efficient Stream Compaction on Wide SIMD Many-Core Architectures. In Proceedings of the Conference on High Performance Graphics (HPG). ACM, 159–166. doi:10.1145/1572769. 1572795

Merge Hulls Without Copying Vertices: GPU Isotonic Regression

[5] Mathieu Blondel, Olivier Teboul, Quentin Berthet, and Josip Djolonga. 2020. Fast Differentiable Sorting and Ranking. In Proceedings of the 37th International Conference on Machine Learning (ICML) (Proceedings of Machine Learning Research, Vol. 119), Hal Daumé III and Aarti Singh (Eds.). PMLR, 950–959. arXiv:2002.08871. [6] Reilly Browne, Rezaul Chowdhury, Shih-Yu Tsai, and Yimin Zhu. 2023. CacheOblivious Parallel Convex Hull in the Binary Forking Model. arXiv preprint arXiv:2305.10389 (2023). [7] Criteo Labs. 2014. Display Advertising Challenge dataset. https://www.kaggle. com/c/criteo-display-ad-challenge. [8] Emil Toftegaard Gæde, Inge Li Gørtz, Ivor van der Hoog, Christoffer Krogh, and Eva Rotenberg. 2024. Simple and Robust Dynamic Two-Dimensional Convex Hull. In Proceedings of the Symposium on Algorithm Engineering and Experiments (ALENEX). SIAM, 144–156. arXiv:2310.18068. doi:10.1137/1.9781611977929.11 [9] Michael T. Goodrich. 1987. Finding the Convex Hull of a Sorted Point Set in Parallel. Inform. Process. Lett. 26, 4 (1987), 173–179. doi:10.1016/0020-0190(87) 90002-0 [10] Guillaume Da Graça and David Defour. 2006. Implementation of Float-Float Operators on Graphics Hardware. In 7th Conference on Real Numbers and Computers (RNC7). Nancy, France, 23–32. arXiv:cs/0603115. [11] Joseph L. Greathouse and Mayank Daga. 2014. Efficient Sparse Matrix-Vector Multiplication on GPUs Using the CSR Storage Format. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC14). IEEE, 769–780. doi:10.1109/SC.2014.68 [12] Xinran He, Junfeng Pan, Ou Jin, Tianbing Xu, Bo Liu, Tao Xu, Yanxin Shi, Antoine Atallah, Ralf Herbrich, Stuart Bowers, and Joaquin Quiñonero Candela. 2014. Practical Lessons from Predicting Clicks on Ads at Facebook. In Proceedings of the Eighth International Workshop on Data Mining for Online Advertising (ADKDD). ACM, 1–9. doi:10.1145/2648584.2648589 [13] Alexander Henzi, Alexandre Mösching, and Lutz Dümbgen. 2022. Accelerating the Pool-Adjacent-Violators Algorithm for Isotonic Distributional Regression. Methodology and Computing in Applied Probability 24, 4 (2022), 2633–2645. arXiv:2006.05527. doi:10.1007/s11009-022-09937-2 [14] Anthony J. Kearsley. 2006. Projections onto Order Simplexes and Isotonic Regression. Journal of Research of the National Institute of Standards and Technology 111, 2 (2006), 121–125. doi:10.6028/jres.111.011 [15] Anthony J. Kearsley, Richard A. Tapia, and Michael W. Trosset. 1996. An Approach to Parallelizing Isotonic Regression. In Applied Mathematics and Parallel Computing: Festschrift for Klaus Ritter, Herbert Fischer, Bruno Riedmüller, and Stefan Schäffler (Eds.). Physica-Verlag, Heidelberg, 141–147. doi:10.1007/978-3642-99789-1_10 [16] Alan Keith, Héctor Ferrada, and Cristóbal A. Navarro. 2022. Accelerating the Convex Hull Computation with a Parallel GPU Algorithm. In 2022 41st International Conference of the Chilean Computer Science Society (SCCC). IEEE, 1–7. arXiv:2209.12310. doi:10.1109/SCCC57464.2022.10000307 [17] Teddy Koker. 2021. torchsort: Fast, differentiable sorting and ranking in PyTorch. https://github.com/teddykoker/torchsort. [18] Alexandru Niculescu-Mizil and Rich Caruana. 2005. Predicting Good Probabilities with Supervised Learning. In Proceedings of the 22nd International Conference on Machine Learning. 625–632. doi:10.1145/1102351.1102430 [19] Mark H. Overmars and Jan van Leeuwen. 1981. Maintenance of Configurations in the Plane. J. Comput. System Sci. 23, 2 (1981), 166–204. doi:10.1016/00220000(81)90012-X [20] John C. Platt. 2000. Probabilities for SV Machines. In Advances in Large Margin Classifiers, Alexander J. Smola, Peter Bartlett, Bernhard Schölkopf, and Dale Schuurmans (Eds.). MIT Press, Cambridge, MA, 61–74. doi:10.7551/mitpress/ 1113.003.0008 [21] Franco P. Preparata and Se June Hong. 1977. Convex Hulls of Finite Sets of Points in Two and Three Dimensions. Commun. ACM 20, 2 (1977), 87–93. doi:10. 1145/359423.359430 [22] Tim Robertson, F. T. Wright, and R. L. Dykstra. 1988. Order Restricted Statistical Inference. Wiley, Chichester. [23] Michael E. Sander, Joan Puigcerver, Josip Djolonga, Gabriel Peyré, and Mathieu Blondel. 2023. Fast, Differentiable and Sparse Top-k: a Convex Analysis Perspective. In Proceedings of the 40th International Conference on Machine Learning (ICML) (Proceedings of Machine Learning Research, Vol. 202), Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (Eds.). PMLR, 29919–29936. arXiv:2302.01425. [24] Jonathan Richard Shewchuk. 1997. Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates. Discrete & Computational Geometry 18, 3 (1997), 305–363. doi:10.1007/PL00009321 [25] Ruoxi Shi, Xinyue Wei, Fanbo Xiang, Zexiang Xu, and Hao Su. 2026. CuACD: A Fully GPU-Resident Approximate Convex Decomposition. arXiv preprint arXiv:2609.28731 (2026). [26] Srikanth Srungarapu, Durga Prasad Reddy, Kishore Kothapalli, and P. J. Narayanan. 2011. Fast Two Dimensional Convex Hull on the GPU. In IEEE Workshops of the International Conference on Advanced Information Networking and Applications (WAINA). IEEE, 7–12. doi:10.1109/WAINA.2011.64

[27] Quentin F. Stout. 2022. Fastest Known Isotonic Regression Algorithms. https: //web.eecs.umich.edu/~qstout/IsoRegAlg.pdf. Status report. [28] Stanley Tzeng and John D. Owens. 2012. Finding Convex Hulls Using Quickhull on the GPU. arXiv:1201.2936 doi:10.48550/arXiv.1201.2936 [29] Bianca Zadrozny and Charles Elkan. 2002. Transforming Classifier Scores into Accurate Multiclass Probability Estimates. In Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 694–699. doi:10.1145/775047.775151

A

The three-case bridge search

The nested search of Section 5.2 tests a candidate vertex of one chain against its tangent to the other, a bisection inside a bisection. The search below halves a candidate range on one of the two chains at every step instead. Lemma A.1 (three cases, lower hulls). Let 𝐴 and 𝐵 be the hulls of two neighbouring nodes without their shared point 𝑠, so that 𝑥 (𝑎) < 𝑥 (𝑠) < 𝑥 (𝑏) for 𝑎 ∈ 𝐴, 𝑏 ∈ 𝐵, and suppose 𝑠 is not a vertex of 𝐻 (𝐴 ∪ 𝐵). Let (𝑎 ∗, 𝑏 ∗ ) be the bridge, with 𝑎 ∗ the leftmost and 𝑏 ∗ the rightmost point of contact of its line. For an edge 𝛼 of 𝐴, an edge 𝛽 of 𝐵, and points 𝑙 ∈ 𝛼, 𝑟 ∈ 𝛽: (1) if slope(𝛼) ≥ slope(𝑙𝑟 ) then 𝑎 ∗ is not right of 𝛼’s left end point; (2) if slope(𝑙𝑟 ) ≥ slope(𝛽) then 𝑏 ∗ is not left of 𝛽’s right end point; (3) otherwise the lines of 𝛼 and 𝛽 meet at a point 𝛾 with 𝑥 (𝑙) < 𝑥 (𝛾) < 𝑥 (𝑟 ), and if 𝑥 (𝛾) ≤ 𝑥 (𝑠) then 𝑎 ∗ is not left of 𝛼’s right end point, while if 𝑥 (𝛾) ≥ 𝑥 (𝑠) then 𝑏 ∗ is not right of 𝛽’s left end point. Proof. Let ℎ be the convex function whose graph is 𝐻 (𝐴 ∪ 𝐵) and ℓ the line of the bridge, which supports 𝐴 and 𝐵 from below. (1) Suppose 𝑎 ∗ lies right of 𝛼’s left end point; as 𝑎 ∗ is a vertex, 𝛼 is then an edge of ℎ left of 𝑎 ∗ , so the line of 𝛼 supports ℎ and 𝑟 𝑦 ≥ ℎ(𝑥 (𝑟 )) ≥ 𝑙 𝑦 + slope(𝛼)(𝑥 (𝑟 ) −𝑥 (𝑙)), i.e. slope(𝑙𝑟 ) ≥ slope(𝛼). With the hypothesis this is an equality, 𝑟 lies on the line of 𝛼, ℎ is linear on [𝑥 (𝑙), 𝑥 (𝑟 )] ∋ 𝑥 (𝑎 ∗ ), and 𝑎 ∗ is not the leftmost point of contact, a contradiction. (2) is symmetric. (3) The slopes differ, so the lines meet. From slope(𝑙𝑟 ) > slope(𝛼) the point 𝑟 lies strictly above the line of 𝛼 and on that of 𝛽, which is below the line of 𝛼 left of 𝛾; so 𝑥 (𝑟 ) > 𝑥 (𝛾), and symmetrically 𝑥 (𝑙) < 𝑥 (𝛾). Let 𝑥 (𝛾) ≤ 𝑥 (𝑠) and suppose 𝑎 ∗ is not right of 𝛼’s left end point. Then 𝛼 is not before the edge out of 𝑎 ∗ in 𝐴, so slope(ℓ) ≤ slope(𝛼), and ℓ is not above the line of 𝛼 at 𝛼’s left end point, hence nowhere to its right. At 𝑥 (𝑏 ∗ ) > 𝑥 (𝑠) ≥ 𝑥 (𝛾) the line of 𝛼 is strictly below the line of 𝛽, which is not above 𝐵: 𝑏 𝑦∗ = ℓ (𝑥 (𝑏 ∗ )) < 𝛽-line(𝑥 (𝑏 ∗ )) ≤ 𝑏 𝑦∗ , a contradiction. The case 𝑥 (𝛾) ≥ 𝑥 (𝑠) is symmetric. □ Each case discards a prefix or suffix of one chain, so the search keeps a contiguous candidate range on each side and halves one of them per step. On segment lists the halving is done in two levels: while a range spans several segments, 𝛼 is the edge between the middle segment’s last vertex and the next segment’s first (an edge of the chain, as consecutive vertices of the hull), and a case keeps either the segments up to the middle one or those after it; within one segment the vertices are consecutive slots and 𝛼 is the middle edge. A side of 𝑘 segments is thus reduced to a vertex in log2 𝑘 + log2 𝐶 + 𝑂 (1) steps, each reading two vertices and two descriptors, and the remaining tangent from that vertex to the other side is the same bisection with a turn predicate. The predicates are signs of polynomials of degree two (slope comparisons) and three (𝑥 (𝛾) against 𝑥 (𝑠)) in the coordinates; any interior points 𝑙, 𝑟 serve, so the edge midpoints are used.

Elena Raikova

Theorem A.2 (merge phase, three-case search). With every walk capped at a constant 𝑇 steps and the search above beyond it, and with a parallel copy for the descriptor moves, the merge phase over 𝑚 = ⌈𝑛/𝐶⌉ chunk hulls does 𝑂 (𝑛 + 𝑚 log 𝑚) work at depth 𝑂 (𝑇 log 𝑛 + log2 𝑛); the work is 𝑂 (𝑛) whenever 𝐶 ≥ log2 𝑚. Proof. The walks and the moves cost as in Theorem 5.2. A search at level 𝑟 takes at most 2(𝑟 + log2 𝐶) + 𝑂 (1) steps, so the Í searches cost 𝑟 (𝑚/2𝑟 +1 + 1) 𝑂 (𝑟 + log 𝐶) = 𝑂 (𝑚 log 𝐶 + log2 𝑚) = 𝑂 (𝑛). The moves cost 𝑂 (𝑚 log 𝑚), which is 𝑂 (𝑛) when 𝐶 ≥ log2 𝑚. The depth of level 𝑟 is at most 𝑇 + 2(𝑟 + log2 𝐶) +𝑂 (1), and the sum over 𝑟 ≤ log2 𝑚 is 𝑂 (𝑇 log 𝑚 + log2 𝑚 + log 𝑚 log 𝐶) = 𝑂 (𝑇 log 𝑛 + log2 𝑛). □ Theorem A.2 bounds a variant of the code. We implemented the three-case search for the levels run with a thread per merge, with a walk capped at 8(𝑟 + log2 𝐶) steps. The levels run with a warp per merge keep the warp-parallel nested search, so the variant does not reach the bound at those levels, and its times show the cost of the search where it is used. On an A100 the variant and the default take the same time when bridges reach across the whole row, with a last point far below the rest: 1.41 against 1.37 ms for the merges of 108 elements. On noisy data the variant’s merges are twice as slow (1.16 against 0.58 ms): few walks reach the cap there, and the variant’s larger kernel lowers occupancy.

B

The filters of the predicates

Merge levels. The merge levels decide turns turn(𝑎, 𝑏, 𝑑) = (𝑏𝑥 − 𝑎𝑥 ) (𝑑 𝑦 − 𝑎 𝑦 ) − (𝑏 𝑦 − 𝑎 𝑦 )(𝑑𝑥 − 𝑎𝑥 ) of vertices addressed by segment and slot. A coordinate of a vertex of chunk 𝑗 is represented as 𝑋 = 𝑐 𝑗 +𝑐 ′𝑗 +ℎ. Here 𝑐 𝑗 and 𝑐 ′𝑗 are the hi and lo parts of the chunk’s start, a renormalized double-double pair with |𝑐 ′𝑗 | ≤ 𝑢|𝑐 𝑗 |, ℎ is the vertex’s coordinate relative to the chunk, a double, and 𝑢 = 2−53 . The slow path computes the turn from differences of these parts (Section 6). The fast path uses absolute coordinates 𝑋ˆ = fl(𝑐 𝑗 + ℎ), one addition each. It computes 𝑏𝑥′ = fl(𝑏ˆ𝑥 − 𝑎ˆ𝑥 ) and likewise 𝑏 𝑦′ , 𝑑𝑥′ , 𝑑 𝑦′ , then 𝑝 1 = fl(𝑏𝑥′ 𝑑 𝑦′ ), 𝑝 2 = fl(𝑏 𝑦′ 𝑑𝑥′ ) and 𝑐 = fl(𝑝 1 − 𝑝 2 ). Error. Let 𝑀 be at least the magnitude of each of the six absolute coordinates and of each chunk start 𝑐 𝑗 of the three points. Then |𝑋ˆ − 𝑋 | ≤ |𝑐 ′𝑗 | + 𝑢 |𝑋ˆ | ≤ 2𝑢𝑀 to first order. The chunk starts matter for the ordinates: a cumulative sum can be far from zero at a chunk’s start and return close to zero within the chunk, so that |𝑐 𝑗 | exceeds the coordinate. For the abscissae, sums of positive weights, the chunk start never exceeds the coordinate. With Σ = |𝑏𝑥′ | + |𝑏 𝑦′ | + |𝑑𝑥′ | + |𝑑 𝑦′ |, the errors of the differences are at most 4𝑢𝑀 +𝑢 |𝑏𝑥′ | and so on. Those of the products are at most 4𝑢𝑀 (|𝑏𝑥′ | + |𝑑 𝑦′ |) + 3𝑢 |𝑝 1 | and 4𝑢𝑀 (|𝑏 𝑦′ | + |𝑑𝑥′ |) + 3𝑢 |𝑝 2 |, and that of 𝑐 is at most 4𝑢𝑀 Σ + 3𝑢 (|𝑝 1 | + |𝑝 2 |) + 𝑢 |𝑐 | ≤ 4𝑢 (𝑀Σ + |𝑝 1 | + |𝑝 2 |)

(4)

to first order, as |𝑐 | ≤ (1 + 𝑢)(|𝑝 1 | + |𝑝 2 |). The code uses twice this, 𝐸 = 8𝑢 (𝑀 Σ + |𝑝 1 | + |𝑝 2 |), which absorbs the second-order terms. If |𝑐 | > 𝐸, then 𝑐 has the sign of the turn of the represented points. The test by exponents. For a normal 𝑣 ≠ 0 let 𝑒 (𝑣) = ⌊log2 |𝑣 |⌋, so ˆ two that 2𝑒 (𝑣) ≤ |𝑣 | < 2𝑒 (𝑣)+1 . Let 𝑒𝑎 be the larger exponent of 𝑎’s coordinates, 𝑒𝑠 the largest exponent of the ordinates’ chunk starts of the three points, 𝑒𝑑 the largest exponent of the four differences,

and 𝑒𝑚 = max(𝑒𝑎 , 𝑒𝑠 , 𝑒𝑑 ) + 1. Every coordinate of 𝑏ˆ and 𝑑ˆ is at most ˆ + |difference| < 2𝑒𝑎 +1 + 2𝑒𝑑 +1 ≤ 2𝑒𝑚 +1 , and every chunk start is |𝑎| below 2𝑒𝑠 +1 ≤ 2𝑒𝑚 , so 𝑀 = 2𝑒𝑚 +1 qualifies. Further, Σ < 4 · 2𝑒𝑑 +1 = 2𝑒𝑑 +3 and |𝑝 1 | + |𝑝 2 | ≤ max(|𝑑𝑥′ |, |𝑑 𝑦′ |)(|𝑏𝑥′ | + |𝑏 𝑦′ |) < 2𝑒𝑑 +1 Σ ≤ 2𝑒𝑚 Σ. Hence 𝐸 < 8𝑢 (2𝑒𝑚 +1 + 2𝑒𝑚 ) Σ < 2−50 · 2𝑒𝑚 +2 · 2𝑒𝑑 +3 = 2𝑒𝑚 +𝑒𝑑 −45 .

(5)

Since |𝑐 | ≥ 2𝑒 (𝑐 ) , the condition 𝑒 (𝑐) ≥ 𝑒𝑚 + 𝑒𝑑 − 44

(6)

implies |𝑐 | ≥ 2𝑒𝑚 +𝑒𝑑 −44 > 𝐸, and the fast result is returned. The code reads the biased exponents 𝑒˜ = 𝑒 + 1023 from the high word of each double with an integer shift and mask. In these, (6) reads 𝑒˜ (𝑐) ≥ 𝑒˜𝑚 + 𝑒˜𝑑 − 1067 with 𝑒˜𝑚 = max(𝑒˜𝑎 , 𝑒˜𝑠 , 𝑒˜𝑑 ) + 1, which takes integer maxima and one comparison. The derivation assumes normal numbers, so a zero or subnormal 𝑐 or difference, whose biased exponent is 0, always takes the slow path. Replacing each quantity by a power of two makes the test coarser than |𝑐 | > 𝐸 by a constant factor, which sends a few more nearly collinear triples to the slow path. What it certifies. A turn that passes the test has the sign that exact arithmetic gives on the represented points, the inputs of the slow path. A turn that fails it is computed by the slow path, in double arithmetic on differences at the scale of the chunks’ sums, which can still err for nearly collinear points. The represented points are themselves not the exact cumulative sums: they carry the rounding of the chunks’ sums and, in float32, of the products 𝑤𝑖 𝑦𝑖 . The predicates are therefore not exact, and the accuracy of Section 7.6 is a measurement. Lane and chunk levels. Here the coordinates are pairs of the working type, with unit roundoff 𝑢, and the fast path uses their hi parts. Within a lane, a coordinate is a renormalized pair whose hi part is its value rounded once, so it is within 𝑢 times the largest magnitude of the lane’s sums, the bound used by the lane chain. Across the lanes of a chunk, a coordinate is the sum of a lane’s start and a sum within the lane, both pairs, and its hi part is the rounded sum of their hi parts. It is within 𝑢 (|𝑦ℎ | + |𝑠 | + |𝑡 |) ≤ 2𝑢 𝑚 of the pair’s value, where 𝑠 and 𝑡 are the two terms and 𝑚 = |𝑠 | + |𝑡 |; 𝑚 can far exceed |𝑦ℎ | when the terms cancel, and the code carries it with each point. For the abscissae, sums of positive weights, |𝑥ℎ | bounds the terms. With 𝑀𝑦 the largest 𝑚 of the three points and 𝑀𝑥 their largest |𝑥ℎ |, the cross product is within 4𝑢 (𝑀𝑦 (|𝑥 1 | + |𝑥 2 |) + 𝑀𝑥 (|𝑦1 | + |𝑦2 |)) + 4𝑢 (|𝑝 1 | + |𝑝 2 |) to first order. The code uses 12𝑢 and 8𝑢 for the two coefficients. The bound is evaluated in the working type, as these levels run at the single-precision rate, and a turn within it is recomputed from the pairs with two-sum and two-product.

C

Gradient validation

Tests. The test suite checks the gradients against finite differences on the CPU, in 𝑦 and in 𝑤, in float64. On the GPU it checks them against the CPU’s on data without near-ties, and against block sums over the GPU’s own block starts on every data family, sawtooth and constant data included. Separately, it checks blocks whose means differ but round to the same float32 value: they stay separate blocks.

Merge Hulls Without Copying Vertices: GPU Isotonic Regression

Soft rank. For the soft rank of Section 7.9 we compare the solvers with a float64 solution of the same float32 problem 𝑧 = 𝑠 − 𝑤. The values of TorchIsotonic agree with it to 5 · 10−8 of their scale. Those of torchsort agree to 6 · 10−8 with 𝜀 = 1/𝑛; with 𝜀 = 1, where each row is one long block whose mean torchsort accumulates in float32, they deviate by 8 · 10−7 at 𝑛 = 100 and by up to 2.5 · 10−4 at 𝑛 = 106 . For the gradients we compare the vector-Jacobian product of the isotonic step, 𝐴𝑔 for a standard normal 𝑔, where 𝐴 averages over the blocks. We report its largest error relative to the reference’s largest entry (ℓ∞ ; the reference’s entries reach 4.6 to 5.0 with 𝜀 = 1/𝑛) and its error in the ℓ2 norm relative to the reference’s. Since 𝐴𝑔 is constant on each block and differs between blocks, its changes mark the blocks each solver used. The gradients of TorchIsotonic agree with the reference to 1.6 · 10−7 in ℓ∞ , with the same blocks, except on the row of 106 with 𝜀 = 1/𝑛. There one of the reference’s 712 735 block boundaries differs, and the means on its two sides are exactly equal: the fit is not differentiable at that point, and each partition gives a valid one-sided derivative. That boundary accounts for an ℓ∞ difference of 0.21 and an ℓ2 difference of 1.4 · 10−3 . torchsort’s gradients agree to 1.1·10−5 with 𝜀 = 1 and for 𝑛 ≤ 103 . With 𝜀 = 1/𝑛 and 𝑛 ≥ 104 , 11, 104 and 1 322 of its block boundaries differ from the reference’s. At each of them the means on the two sides differ, by up to 5·10−8 of the data’s scale, which is below float32’s resolution. These differences

are therefore not explained by equal means of neighbouring blocks; we compared only those means, not the stability of the partition within the blocks. torchsort’s gradient then differs by 0.43 to 0.60 in ℓ∞ and by 0.5 to 4% in ℓ2 .

D

Calibration metrics

Table 7: Calibration of the click model of Section 7.7 on day 7: expected calibration error (ECE) and the change in logloss, with 95% intervals paired over the rows of day 7. One isotonic curve is compared with Platt scaling, one curve per segment with one isotonic curve on the segments’ rows. The model’s AUC is 0.767; uncalibrated, it predicts a mean click probability of 0.67 for a click rate of 0.26. calibration none downsampling correction Platt scaling [20] isotonic, one curve

ECE

Δ log-loss

0.41 0.032 0.0080 0.0015

−5.0 (−5.3 to −4.7) ·10−4

segments (60% of rows), one curve 0.0023 one curve per segment 0.0026 −1.24 (−1.42 to −1.07) ·10−3

Related documents

• Merge Hulls Without Copying Vertices: GPU Isotonic Regression#1388604
• Correlation-Based Distillation Yields More Mergeable Speech-Music Encoders#1187150
• GPU-Accelerated Computation of Persistent Homology for Topological Analysis of Image Data#1247886
• An Iterative Method with Asymptotic Orthogonality for Simultaneous Eigenpair Computation#925887
• Same Flow, Different Paths: Variance Reduction in Flow Matching#1123663