[BUG] Batched silhouette score can differ from non-batched score under process concurrency

[csadorf]

Describe the bug

On an A100, concurrent processes can intermittently produce a batched silhouette score that differs by more than 1e-4 from the non-batched cuVS score for the same deterministic input.

In one run, 14 mismatches occurred across 49,152 comparisons from 512 deterministic seeds executed as eight concurrent processes. Every affected seed passed when rerun serially. The observed discrepancies were approximately 1e-3 to 1e-2.

Steps/Code to reproduce bug

Build and run:

$CONDA_PREFIX/bin/nvcc -ccbin $CONDA_PREFIX/bin/x86_64-conda-linux-gnu-g++ \
  -std=c++20 -I$CONDA_PREFIX/include -I$CONDA_PREFIX/include/rapids \
  -I$CONDA_PREFIX/targets/x86_64-linux/include -L$CONDA_PREFIX/lib \
  -Xlinker -rpath -Xlinker $CONDA_PREFIX/lib repro.cu -lcuvs -lrmm -o repro
seq 0 511 | xargs -P8 -n1 ./repro

repro.cu

#include <cuvs/stats/silhouette_score.hpp>

#include <raft/core/device_mdarray.hpp>
#include <raft/core/resource/cuda_stream.hpp>
#include <raft/core/resource/cuda_stream_pool.hpp>
#include <raft/util/cudart_utils.hpp>

#include <rmm/cuda_stream_pool.hpp>

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <memory>
#include <optional>
#include <random>
#include <string_view>
#include <utility>
#include <vector>

namespace {

struct metric_case {
  std::string_view name;
  cuvs::distance::DistanceType metric;
};

constexpr std::array<metric_case, 4> metrics{{
  {"cosine", cuvs::distance::DistanceType::CosineExpanded},
  {"euclidean", cuvs::distance::DistanceType::L2SqrtUnexpanded},
  {"sqeuclidean", cuvs::distance::DistanceType::L2Expanded},
  {"l1", cuvs::distance::DistanceType::L1},
}};

}  // namespace

int main(int argc, char** argv)
{
  constexpr int64_t rows       = 1000;
  constexpr int64_t cols       = 2;
  constexpr int labels         = 2;
  constexpr int repetitions    = 4;
  constexpr float tolerance    = 1e-4f;
  constexpr std::array<int64_t, 3> chunks{rows, rows / 3, rows / 5};

  auto seed = argc > 1 ? std::strtoul(argv[1], nullptr, 10) : 0;

  std::mt19937 rng(seed);
  std::uniform_real_distribution<float> centers(-1.0f, 1.0f);
  std::normal_distribution<float> noise(0.0f, 1.5f);

  std::array<std::array<float, cols>, labels> center{};
  for (auto& cluster : center) {
    for (auto& value : cluster) {
      value = centers(rng);
    }
  }

  std::vector<int64_t> order(rows);
  for (int64_t i = 0; i < rows; ++i) {
    order[i] = i;
  }
  std::shuffle(order.begin(), order.end(), rng);

  std::vector<float> X(rows * cols);
  std::vector<int> y(rows);
  for (int64_t row = 0; row < rows; ++row) {
    auto label = static_cast<int>(order[row] / (rows / labels));
    y[row]     = label;
    for (int64_t col = 0; col < cols; ++col) {
      X[row * cols + col] = center[label][col] + noise(rng);
    }
  }

  raft::resources default_handle;
  raft::resources pool_handle;
  raft::resource::set_cuda_stream_pool(pool_handle, std::make_shared<rmm::cuda_stream_pool>(4));

  auto stream = raft::resource::get_cuda_stream(default_handle);
  auto d_X     = raft::make_device_matrix<float, int64_t>(default_handle, rows, cols);
  auto d_y     = raft::make_device_vector<int, int64_t>(default_handle, rows);

  raft::update_device(d_X.data_handle(), X.data(), X.size(), stream);
  raft::update_device(d_y.data_handle(), y.data(), y.size(), stream);
  raft::resource::sync_stream(default_handle);

  auto X_view =
    raft::make_device_matrix_view<const float, int64_t>(d_X.data_handle(), rows, cols);
  auto y_view = raft::make_device_vector_view<const int, int64_t>(d_y.data_handle(), rows);

  bool failed = false;
  for (auto const& metric : metrics) {
    auto non_batched = cuvs::stats::silhouette_score(
      default_handle, X_view, y_view, std::nullopt, labels, metric.metric);

    for (auto const& handle :
         {std::pair{"default", &default_handle}, std::pair{"pool", &pool_handle}}) {
      for (auto chunk : chunks) {
        for (int repetition = 0; repetition < repetitions; ++repetition) {
          auto batched = cuvs::stats::silhouette_score_batched(
            *handle.second, X_view, y_view, std::nullopt, labels, chunk, metric.metric);

          auto difference = std::abs(batched - non_batched);
          if (difference > tolerance) {
            failed = true;
            std::cerr << std::setprecision(10) << "seed=" << seed
                      << " handle=" << handle.first << " metric=" << metric.name
                      << " chunk=" << chunk << " repetition=" << repetition
                      << " non-batched=" << non_batched << " batched=" << batched
                      << " difference=" << difference << '\n';
          }
        }
      }
    }
  }

  return failed ? EXIT_FAILURE : EXIT_SUCCESS;
}

Expected behavior

Batched and non-batched silhouette scores should agree within 1e-4 for the same input.

The tolerance intentionally allows ordinary floating-point reduction-order variation.

Environment details (please complete the following information):

• GPU: A100 / SM80
• NVIDIA driver: 560.35.05
• CUDA: 12.2.140
• libcuvs and headers: 26.08.00a117 cuda12_260731_2db820e7
• RAFT headers: 26.08.00a61 cuda12_260804_97d22c81
• RMM: 26.08.00a62 cuda12_260804_1a39f9e8
• Method of cuVS install: conda

Additional context

The failure depends on process concurrency: every affected deterministic seed passed when rerun serially.