diff --git a/src/treelearner/cuda/cuda_best_split_finder.cu b/src/treelearner/cuda/cuda_best_split_finder.cu
index 24b03fb5df37..ff1c6bea9422 100644
--- a/src/treelearner/cuda/cuda_best_split_finder.cu
+++ b/src/treelearner/cuda/cuda_best_split_finder.cu
@@ -30,6 +30,30 @@ __device__ __forceinline__ int CUDARoundInt(double x) {
   return static_cast<int>(x + 0.5);
 }
 
+// Relative tolerance for gain-tie detection in the best-split reductions.
+// Two gains within this relative band are treated as a plateau; the
+// reduction tie-breaks to the lower thread index, which corresponds to
+// the lower bin index — matching CPU's "first bin in scan order wins"
+// behaviour. Sized at ~5000x fp64 epsilon, well above the ~1e-15
+// reduction-order noise we observe in CUDA histograms but well below
+// any genuine gain difference seen in practice. Without this, ULP-level
+// FP noise causes CPU and CUDA to pick different bins from a true gain
+// plateau, which seeds the round-3+ structural divergence in cases like
+// reg_bagging and multi_dense.
+#define BEST_GAIN_TIE_RELATIVE_TOL 1e-12
+
+__device__ __forceinline__ bool OtherIsBetterWithTieBreak(
+    double other_gain, double gain, uint32_t other_thread_index, uint32_t thread_index) {
+  const double tol = fmax(fabs(gain), fabs(other_gain)) * BEST_GAIN_TIE_RELATIVE_TOL;
+  if (other_gain > gain + tol) {
+    return true;
+  }
+  if (fabs(other_gain - gain) <= tol && other_thread_index < thread_index) {
+    return true;
+  }
+  return false;
+}
+
 __device__ void ReduceBestGainWarp(double gain, bool found, uint32_t thread_index, double* out_gain, bool* out_found, uint32_t* out_thread_index) {
   const uint32_t mask = 0xffffffff;
   const uint32_t warpLane = threadIdx.x % warpSize;
@@ -37,7 +61,8 @@ __device__ void ReduceBestGainWarp(double gain, bool found, uint32_t thread_inde
     const bool other_found = __shfl_down_sync(mask, found, offset);
     const double other_gain = __shfl_down_sync(mask, gain, offset);
     const uint32_t other_thread_index = __shfl_down_sync(mask, thread_index, offset);
-    if ((other_found && found && other_gain > gain) || (!found && other_found)) {
+    const bool other_better = OtherIsBetterWithTieBreak(other_gain, gain, other_thread_index, thread_index);
+    if ((other_found && found && other_better) || (!found && other_found)) {
       found = other_found;
       gain = other_gain;
       thread_index = other_thread_index;
@@ -56,7 +81,8 @@ __device__ uint32_t ReduceBestGainBlock(double gain, bool found, uint32_t thread
     const bool other_found = __shfl_down_sync(mask, found, offset);
     const double other_gain = __shfl_down_sync(mask, gain, offset);
     const uint32_t other_thread_index = __shfl_down_sync(mask, thread_index, offset);
-    if ((other_found && found && other_gain > gain) || (!found && other_found)) {
+    const bool other_better = OtherIsBetterWithTieBreak(other_gain, gain, other_thread_index, thread_index);
+    if ((other_found && found && other_better) || (!found && other_found)) {
       found = other_found;
       gain = other_gain;
       thread_index = other_thread_index;
diff --git a/tests/python_package_test/test_dual.py b/tests/python_package_test/test_dual.py
index c910a4318260..e66e8e24f184 100644
--- a/tests/python_package_test/test_dual.py
+++ b/tests/python_package_test/test_dual.py
@@ -968,3 +968,84 @@ def test_cuda_linear_tree_handles_nan_like_cpu():
         preds[device_type] = lgb.train(params, ds, num_boost_round=30).predict(X)
     max_diff = float(np.max(np.abs(preds["cpu"] - preds["cuda"])))
     assert max_diff < 1e-6, f"CUDA linear tree (NaN) diverges from CPU: max|diff|={max_diff:.3e}"
+
+
+def _make_regression_for_parity(n=200, d=8, seed=0):
+    rng = np.random.default_rng(seed)
+    X = rng.standard_normal((n, d)).astype(np.float64)
+    coef = rng.standard_normal(d)
+    y = (X @ coef + 0.1 * rng.standard_normal(n)).astype(np.float64)
+    return X, y
+
+
+def _train_cpu_and_cuda(params_overrides, X, y, num_round):
+    out = {}
+    for device_type in ("cpu", "cuda"):
+        params = {
+            "objective": "regression",
+            "verbose": -1,
+            "deterministic": True,
+            "num_threads": 1,
+            "seed": 42,
+            "feature_pre_filter": False,
+            "device_type": device_type,
+            "gpu_use_dp": True,
+            "force_col_wise": True,
+            "num_leaves": 7,
+            "learning_rate": 0.1,
+            "min_data_in_leaf": 5,
+            "min_sum_hessian_in_leaf": 1e-3,
+            **params_overrides,
+        }
+        ds = lgb.Dataset(X, label=y, params={"verbose": -1, "feature_pre_filter": False})
+        out[device_type] = lgb.train(params, ds, num_boost_round=num_round)
+    return out
+
+
+@_REQUIRES_CUDA
+@pytest.mark.parametrize(
+    ("name", "params_overrides", "seed", "num_round"),
+    [
+        # Regression test for the gain-plateau argmax bug. Bagging configuration
+        # was the original failure: max|Δ|=0.39 at round 3, structurally
+        # divergent trees from round 3 onward. With the fix, all 5 rounds
+        # match at fp64 epsilon and trees are bit-identical.
+        (
+            "bagging",
+            {"bagging_fraction": 0.7, "bagging_freq": 1, "bagging_seed": 1},
+            11,
+            5,
+        ),
+        # Plain dense regression: predictions matched at fp64 epsilon prior
+        # to the fix but the encoded tree thresholds differed cosmetically
+        # (CPU and CUDA picked different bins from a true gain plateau).
+        # After the fix, trees are bit-identical.
+        ("dense", {}, 1, 5),
+        # max_depth regression: another configuration where round-1 trees
+        # had cosmetic threshold-encoding differences before the fix.
+        ("max_depth", {"max_depth": 3}, 12, 5),
+        # L2 regularisation: same family.
+        ("l2", {"lambda_l2": 1.0}, 7, 5),
+    ],
+)
+def test_cuda_split_gain_tie_break_matches_cpu(name, params_overrides, seed, num_round):
+    """CUDA must match CPU at fp64 epsilon when the best-split argmax has a
+    gain plateau (multiple bins with truly equal gain).
+
+    Prior to the tolerance-based tie-break in cuda_best_split_finder.cu's
+    ReduceBestGain* helpers, ULP-level FP noise in the parallel histogram
+    flipped which bin from the plateau had the slightly-higher numerical
+    gain on CUDA. CPU's exact computation picked a different bin (its
+    sequential scan + strict ``>`` comparison resolves to the lowest-index
+    bin on the plateau). The threshold-encoding mismatch was at round 1
+    cosmetic for predictions (data routed identically), but compounded
+    through score updates and surfaced as structural tree divergence by
+    round 3 in cases like reg_bagging.
+    """
+    X, y = _make_regression_for_parity(seed=seed)
+    pair = _train_cpu_and_cuda(params_overrides, X, y, num_round=num_round)
+    pred_cpu = pair["cpu"].predict(X, raw_score=True)
+    pred_cuda = pair["cuda"].predict(X, raw_score=True)
+    # fp64 epsilon ≈ 2.2e-16; allow a generous 1e-10 to absorb any
+    # remaining round-by-round drift from sources unrelated to this fix.
+    np.testing.assert_allclose(pred_cuda, pred_cpu, atol=1e-10)