From fb4e8edc95b2e00a18870f0e3f758ec44095f480 Mon Sep 17 00:00:00 2001 From: Max Buckley Date: Mon, 31 Aug 2026 22:01:21 +0200 Subject: [PATCH 1/2] 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 de7df75128fd1955af56ec7068ca5cf7eb360dad Mon Sep 17 00:00:00 2001 From: Max Buckley Date: Thu, 3 Sep 2026 10:18:09 +0200 Subject: [PATCH 2/2] Use the canonical NVIDIA copyright notice Applied by `pre-commit run --all-files`; the verify-copyright hook requires the "& AFFILIATES. All rights reserved." form. --- cpp/src/distance/detail/kernels/kernel_matrices.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 */