From 23ef3b68e5119e49aace58dc055c624703861dbd Mon Sep 17 00:00:00 2001 From: Max Buckley Date: Mon, 31 Aug 2026 21:50:01 +0200 Subject: [PATCH 1/3] Keep the sparse Jensen-Shannon post-processing in the input precision The final step of the sparse Jensen-Shannon distance is [=] __device__(value_t input) { return raft::sqrt(0.5 * input); } `0.5` is a `double`, so for `value_t == float` the multiply promotes the expression and `raft::sqrt` resolves to the fp64 overload. Unlike a stray `double` in an add or multiply, an fp64 sqrt is not one instruction: it is a Newton-Raphson refinement sequence. The generated SASS for the float instantiation of the map kernel was DMUL R14, R4, 0.5 MUFU.RSQ64H R13, R15 DMUL R16, R12, R12 DFMA R16, R14, -R16, 1 ... 6 more DMUL/DFMA ten fp64 instructions per element, which run at 1/64 the fp32 rate on consumer parts. Writing the constant as `value_t(0.5)` picks the float overload and reduces this to `FMUL` + `MUFU.RSQ` + four float ops. `0.5` is exactly representable in binary, so the `double` instantiation is bit-for-bit unchanged. The `float` instantiation now rounds once instead of twice, so results move by at most an ulp. Measured on an RTX 5090 (sm_120a, CUDA 13.2), random CSR inputs, median of 5: rows x cols, nnz/row map kernel full pairwise_distance 4096 x 4096, 32 263 -> 25 us 6.44 -> 6.20 ms 1.04x 8192 x 4096, 32 1046 -> 319 us 25.31 -> 24.57 ms 1.03x 16384 x 4096, 8 4167 -> 1395 us 29.21 -> 26.42 ms 1.11x 16384 x 16384, 32 4169 -> 1394 us 130.30 -> 127.53 ms 1.02x The map kernel itself is 3.0-10.4x faster. At 16384 rows it now moves 2.1 GB in 1.39 ms, i.e. it has gone from fp64-throughput-bound to sitting at the memory roofline; the 4096 case is faster still because its output fits in L2. End-to-end gains are smaller because the balanced COO SpMV that produces the input dominates the call. Output checksums are unchanged to six decimal places across all four shapes. --- cpp/src/distance/detail/sparse/lp_distance.cuh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/distance/detail/sparse/lp_distance.cuh b/cpp/src/distance/detail/sparse/lp_distance.cuh index 38025329b9..f396b647ef 100644 --- a/cpp/src/distance/detail/sparse/lp_distance.cuh +++ b/cpp/src/distance/detail/sparse/lp_distance.cuh @@ -266,7 +266,7 @@ class jensen_shannon_unexpanded_distances_t : public distances_t { raft::linalg::map( config_->handle, raft::make_device_vector_view(out_dists, n), - [=] __device__(value_t input) { return raft::sqrt(0.5 * input); }, + [=] __device__(value_t input) { return raft::sqrt(value_t(0.5) * input); }, raft::make_const_mdspan(raft::make_device_vector_view(out_dists, n))); } From df07908f24b0c2789e0a41cc790629cd36b81d09 Mon Sep 17 00:00:00 2001 From: Max Buckley Date: Mon, 31 Aug 2026 22:01:21 +0200 Subject: [PATCH 2/3] Keep the polynomial and RBF Gram epilogues in the input precision `polynomial_kernel`, `polynomial_kernel_nopad` and `rbf_kernel_expanded` all evaluate their epilogue in fp64 even when instantiated on `float`: - `pow(gain * x + offset, exponent)` has `exp_t == int`, so overload resolution picks `double pow(double, double)` and promotes the whole expression. - `exp(-1.0 * gain * (...))` promotes on the `-1.0` literal, and then resolves to the fp64 `exp`. `pow` and `exp` are not single instructions; both expand to a polynomial evaluation, so the promotion multiplies out. Before this change the `float` instantiations contained 92 and 19 fp64 instructions respectively; after, they contain none. The `double` instantiations are unchanged (98 and 22). Casting the exponent to `math_t` and dropping the `-1.0` selects the float overloads without changing the `double` path. ## Cost of the epilogue Measured with `KernelFactory::create(...)->evaluate()` on random column-major inputs, RTX 5090 (sm_120a), CUDA 13.2, median of 5. `linear` is the same GEMM with no epilogue; `tanh` has the same epilogue structure but no `double` literal, so both act as controls. float: shape kernel before after speedup 8192^2, d=128 linear 0.316 ms 0.315 ms 1.00x (control) 8192^2, d=128 poly 8.669 ms 0.652 ms 13.3x 8192^2, d=128 rbf 2.189 ms 0.746 ms 2.93x 8192^2, d=128 tanh 0.662 ms 0.664 ms 1.00x (control) 16384^2, d=32 linear 0.696 ms 0.695 ms 1.00x (control) 16384^2, d=32 poly 33.952 ms 2.112 ms 16.1x 16384^2, d=32 rbf 7.831 ms 2.198 ms 3.56x 16384^2, d=32 tanh 2.113 ms 2.112 ms 1.00x (control) 4096^2, d=1024 poly 2.725 ms 0.703 ms 3.88x 4096^2, d=1024 rbf 1.160 ms 0.752 ms 1.54x double is unchanged to within 0.05% on every shape. Before this change the polynomial epilogue cost 24-43x the GEMM it decorates, and the `float` path (33.95 ms) was within 1.4x of the `double` path (47.56 ms). Afterwards poly, tanh and linear+epilogue all land within a few percent of each other, i.e. the epilogue is memory-bound as it should be. ## Effect on cuML's SVM cuML's SVC calls these through its kernel cache. Built cuML 26.10 against this branch (`CPM_cuvs_SOURCE`), blobs data as in cuML's own `bench/sg/svc.cu`, median of 3, run to convergence. Solver iteration counts are reported so that a changed convergence path is not mistaken for a speedup. float, converged before after speedup n_iter 50000x1000 poly 9.43 ms 6.26 ms 1.51x 500 -> 500 50000x1000 rbf 12.04 ms 10.44 ms 1.15x 1319 -> 1319 50000x2 poly 808.60 ms 686.71 ms 1.18x 338727 -> 323526 50000x2 rbf 89.42 ms 62.36 ms 1.43x 10436 -> 9219 50000x1000 linear 4.50 ms 4.38 ms 1.03x 192 -> 192 (control) 50000x1000 tanh 12.89 ms 12.59 ms 1.02x 771 -> 771 (control) The two 50000x1000 rows converge in an identical number of iterations, so 1.51x and 1.15x are like-for-like. The 50000x2 rows changed iteration count; per iteration they are 1.12x and 1.27x. The 2048x100000 shape is ~1.01x, since at d=100000 the GEMM dominates. All `double` SVM configurations are unchanged. ## Numerics The `float` path now rounds once instead of going through fp64. Measured against the previous fp64-then-round result over the value range these kernels actually see (4M samples): polynomial degree=2 max relative error 1.19e-07 (1 ulp) polynomial degree=3 max relative error 1.19e-07 (1 ulp) polynomial degree=4 max relative error 1.19e-07 (1 ulp) rbf max relative error 1.67e-07 (2 ulp) `double` results are bit-identical. --- cpp/src/distance/detail/kernels/kernel_matrices.cu | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/src/distance/detail/kernels/kernel_matrices.cu b/cpp/src/distance/detail/kernels/kernel_matrices.cu index 9ed25f959c..d0d5251671 100644 --- a/cpp/src/distance/detail/kernels/kernel_matrices.cu +++ b/cpp/src/distance/detail/kernels/kernel_matrices.cu @@ -29,7 +29,7 @@ RAFT_KERNEL polynomial_kernel_nopad( { for (size_t tid = threadIdx.x + blockIdx.x * blockDim.x; tid < len; tid += blockDim.x * gridDim.x) { - inout[tid] = pow(gain * inout[tid] + offset, exponent); + inout[tid] = pow(gain * inout[tid] + offset, (math_t)exponent); } } @@ -51,7 +51,7 @@ RAFT_KERNEL polynomial_kernel( tidy += blockDim.y * gridDim.y) for (size_t tidx = threadIdx.x + blockIdx.x * blockDim.x; tidx < rows; tidx += blockDim.x * gridDim.x) { - inout[tidx + tidy * ld] = pow(gain * inout[tidx + tidy * ld] + offset, exponent); + inout[tidx + tidy * ld] = pow(gain * inout[tidx + tidy * ld] + offset, (math_t)exponent); } } @@ -118,7 +118,7 @@ RAFT_KERNEL rbf_kernel_expanded( for (size_t tidx = threadIdx.x + blockIdx.x * blockDim.x; tidx < rows; tidx += blockDim.x * gridDim.x) { inout[tidx + tidy * ld] = - exp(-1.0 * gain * (norm_x[tidx] + norm_y_val - inout[tidx + tidy * ld] * 2)); + exp(-gain * (norm_x[tidx] + norm_y_val - inout[tidx + tidy * ld] * 2)); } } } From ba762fd363d19a8ac68fb824ce70fbc7f415bacc Mon Sep 17 00:00:00 2001 From: Max Buckley Date: Thu, 3 Sep 2026 15:55:53 +0200 Subject: [PATCH 3/3] Use the canonical NVIDIA copyright notice --- cpp/src/distance/detail/kernels/kernel_matrices.cu | 2 +- cpp/src/distance/detail/sparse/lp_distance.cuh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/distance/detail/kernels/kernel_matrices.cu b/cpp/src/distance/detail/kernels/kernel_matrices.cu index d0d5251671..8b631d795c 100644 --- a/cpp/src/distance/detail/kernels/kernel_matrices.cu +++ b/cpp/src/distance/detail/kernels/kernel_matrices.cu @@ -1,5 +1,5 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2019-2026, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2019-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ diff --git a/cpp/src/distance/detail/sparse/lp_distance.cuh b/cpp/src/distance/detail/sparse/lp_distance.cuh index f396b647ef..8533ed83c8 100644 --- a/cpp/src/distance/detail/sparse/lp_distance.cuh +++ b/cpp/src/distance/detail/sparse/lp_distance.cuh @@ -1,5 +1,5 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */