From 9ad43ea37609f9ea521167692ca0acb3d990cd5a Mon Sep 17 00:00:00 2001 From: Tim Head Date: Wed, 8 Jul 2026 08:49:31 +0200 Subject: [PATCH 1/5] Add basic RidgeCV estimator --- docs/source/api/cuml.linear_model.rst | 1 + docs/source/api/index.rst | 3 + docs/source/cuml-accel/limitations.rst | 9 + python/cuml/cuml/__init__.py | 2 + .../accel/_overrides/sklearn/linear_model.py | 10 + python/cuml/cuml/benchmark/algorithms.py | 10 + .../benchmark/automated/bench_regression.py | 6 +- .../cuml/benchmark/configs/single_gpu.yaml | 24 + python/cuml/cuml/linear_model/__init__.py | 3 +- python/cuml/cuml/linear_model/ridge_cv.py | 574 ++++++++++++++++++ .../integration/test_ridge_cv.py | 74 +++ python/cuml/tests/test_ridge_cv.py | 191 ++++++ 12 files changed, 905 insertions(+), 2 deletions(-) create mode 100644 python/cuml/cuml/linear_model/ridge_cv.py create mode 100644 python/cuml/cuml_accel_tests/integration/test_ridge_cv.py create mode 100644 python/cuml/tests/test_ridge_cv.py diff --git a/docs/source/api/cuml.linear_model.rst b/docs/source/api/cuml.linear_model.rst index 063546c148..1f315ff8aa 100644 --- a/docs/source/api/cuml.linear_model.rst +++ b/docs/source/api/cuml.linear_model.rst @@ -13,6 +13,7 @@ cuml.linear_model LinearRegression LogisticRegression Ridge + RidgeCV Lasso ElasticNet MBSGDClassifier diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst index 94822a4a19..ef2c462b97 100644 --- a/docs/source/api/index.rst +++ b/docs/source/api/index.rst @@ -295,6 +295,9 @@ Regression and Classification * - :obj:`~cuml.linear_model.Ridge` - Ridge regression. - :mod:`cuml.linear_model` + * - :obj:`~cuml.linear_model.RidgeCV` + - Ridge regression with built-in cross-validation. + - :mod:`cuml.linear_model` * - :obj:`~cuml.solvers.SGD` - Stochastic Gradient Descent solver. - :mod:`cuml.solvers` diff --git a/docs/source/cuml-accel/limitations.rst b/docs/source/cuml-accel/limitations.rst index 9216b3007e..cdffaae8c1 100644 --- a/docs/source/cuml-accel/limitations.rst +++ b/docs/source/cuml-accel/limitations.rst @@ -301,6 +301,15 @@ Ridge - If ``positive=True`` or ``solver="lbfgs"``. +RidgeCV +^^^^^^^ + +``RidgeCV`` will fall back to CPU in the following cases: + +- If ``cv`` is not ``None`` (only the default Leave-One-Out/GCV path is + accelerated on GPU). +- If a custom ``scoring`` is provided. + Lasso ^^^^^ diff --git a/python/cuml/cuml/__init__.py b/python/cuml/cuml/__init__.py index c1d0a62d1e..11044e4ba2 100644 --- a/python/cuml/cuml/__init__.py +++ b/python/cuml/cuml/__init__.py @@ -52,6 +52,7 @@ from cuml.linear_model.mbsgd_classifier import MBSGDClassifier from cuml.linear_model.mbsgd_regressor import MBSGDRegressor from cuml.linear_model.ridge import Ridge +from cuml.linear_model.ridge_cv import RidgeCV from cuml.manifold.t_sne import TSNE from cuml.manifold.umap import UMAP from cuml.metrics import accuracy_score, adjusted_rand_score, r2_score @@ -137,6 +138,7 @@ def __getattr__(name): "RandomForestClassifier", "RandomForestRegressor", "Ridge", + "RidgeCV", "SGD", "SparseRandomProjection", "SVC", diff --git a/python/cuml/cuml/accel/_overrides/sklearn/linear_model.py b/python/cuml/cuml/accel/_overrides/sklearn/linear_model.py index f2ac0e1f05..a35a593f41 100644 --- a/python/cuml/cuml/accel/_overrides/sklearn/linear_model.py +++ b/python/cuml/cuml/accel/_overrides/sklearn/linear_model.py @@ -19,6 +19,7 @@ "LogisticRegression", "ElasticNet", "Ridge", + "RidgeCV", "Lasso", ) @@ -66,6 +67,15 @@ def _gpu_fit(self, X, y, sample_weight=None): return self +class RidgeCV(ProxyBase): + _gpu_class = cuml.linear_model.RidgeCV + + def _gpu_fit(self, X, y, sample_weight=None, **params): + # `**params` is only used by sklearn for metadata routing, which we + # don't support; it's absorbed here to match sklearn's signature. + return self._gpu.fit(X, y, sample_weight=sample_weight) + + class _ElasticNetMixin: def _gpu_fit(self, X, y, sample_weight=None, check_input=True): # check_input is ignored, only here to fix signature mismatch with sklearn diff --git a/python/cuml/cuml/benchmark/algorithms.py b/python/cuml/cuml/benchmark/algorithms.py index 1006ec9229..268fbfcdc6 100644 --- a/python/cuml/cuml/benchmark/algorithms.py +++ b/python/cuml/cuml/benchmark/algorithms.py @@ -330,6 +330,7 @@ def all_algorithms(): cuml_ElasticNet = cuml.linear_model.ElasticNet cuml_Lasso = cuml.linear_model.Lasso cuml_Ridge = cuml.linear_model.Ridge + cuml_RidgeCV = cuml.linear_model.RidgeCV cuml_KernelRidge = cuml.kernel_ridge.KernelRidge cuml_LogisticRegression = cuml.linear_model.LogisticRegression cuml_RandomForestClassifier = cuml.ensemble.RandomForestClassifier @@ -356,6 +357,7 @@ def all_algorithms(): cuml_LinearRegression = cuml_ElasticNet = cuml_Lasso = cuml_Ridge = ( None ) + cuml_RidgeCV = None cuml_KernelRidge = None cuml_LogisticRegression = None cuml_RandomForestClassifier = cuml_RandomForestRegressor = None @@ -496,6 +498,14 @@ def all_algorithms(): accepts_labels=True, accuracy_function=metrics.r2_score, ), + AlgorithmPair( + sklearn.linear_model.RidgeCV, + cuml_RidgeCV, + shared_args={"alphas": np.logspace(-3, 3, 50)}, + name="RidgeCV", + accepts_labels=True, + accuracy_function=metrics.r2_score, + ), AlgorithmPair( sklearn.kernel_ridge.KernelRidge, cuml_KernelRidge, diff --git a/python/cuml/cuml/benchmark/automated/bench_regression.py b/python/cuml/cuml/benchmark/automated/bench_regression.py index 99dd9e66bf..d40dcaa1dc 100644 --- a/python/cuml/cuml/benchmark/automated/bench_regression.py +++ b/python/cuml/cuml/benchmark/automated/bench_regression.py @@ -1,5 +1,5 @@ # -# SPDX-FileCopyrightText: Copyright (c) 2021-2025, NVIDIA CORPORATION. +# SPDX-FileCopyrightText: Copyright (c) 2021-2026, NVIDIA CORPORATION. # SPDX-License-Identifier: Apache-2.0 # @@ -64,6 +64,10 @@ def bench_ridge(gpubenchmark, bench_step, regression1): # noqa: F811 _benchmark_algo(gpubenchmark, "Ridge", bench_step, regression1) +def bench_ridge_cv(gpubenchmark, bench_step, regression1): # noqa: F811 + _benchmark_algo(gpubenchmark, "RidgeCV", bench_step, regression1) + + def bench_knnregressor(gpubenchmark, bench_step, regression1): # noqa: F811 _benchmark_algo( gpubenchmark, "KNeighborsRegressor", bench_step, regression1 diff --git a/python/cuml/cuml/benchmark/configs/single_gpu.yaml b/python/cuml/cuml/benchmark/configs/single_gpu.yaml index 8097e283f6..df14e79959 100644 --- a/python/cuml/cuml/benchmark/configs/single_gpu.yaml +++ b/python/cuml/cuml/benchmark/configs/single_gpu.yaml @@ -94,6 +94,30 @@ benchmarks: default: {rows: [2600000]} nightly: {rows: [10400000]} + - id: ridge_cv_fit + algorithm: RidgeCV + dataset: regression + operation: fit + tags: [linear, regression] + # RidgeCV runs a leave-one-out GCV over 50 alphas (O(n*p^2) per alpha), so + # it is much heavier than a single Ridge fit; row counts are scaled down. + variants: + narrow: + features: [16] + tiers: + default: {rows: [8000000]} + nightly: {rows: [16000000]} + medium: + features: [128] + tiers: + default: {rows: [1000000]} + nightly: {rows: [4000000]} + wide: + features: [512] + tiers: + default: {rows: [250000]} + nightly: {rows: [1000000]} + - id: elasticnet_fit algorithm: ElasticNet dataset: regression diff --git a/python/cuml/cuml/linear_model/__init__.py b/python/cuml/cuml/linear_model/__init__.py index a9e869be2e..9fc08295b8 100644 --- a/python/cuml/cuml/linear_model/__init__.py +++ b/python/cuml/cuml/linear_model/__init__.py @@ -1,5 +1,5 @@ # -# SPDX-FileCopyrightText: Copyright (c) 2019-2025, NVIDIA CORPORATION. +# SPDX-FileCopyrightText: Copyright (c) 2019-2026, NVIDIA CORPORATION. # SPDX-License-Identifier: Apache-2.0 # @@ -10,3 +10,4 @@ from cuml.linear_model.mbsgd_classifier import MBSGDClassifier from cuml.linear_model.mbsgd_regressor import MBSGDRegressor from cuml.linear_model.ridge import Ridge +from cuml.linear_model.ridge_cv import RidgeCV diff --git a/python/cuml/cuml/linear_model/ridge_cv.py b/python/cuml/cuml/linear_model/ridge_cv.py new file mode 100644 index 0000000000..0fc3e49a2d --- /dev/null +++ b/python/cuml/cuml/linear_model/ridge_cv.py @@ -0,0 +1,574 @@ +# +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. +# SPDX-License-Identifier: Apache-2.0 +# +import cupy as cp +import numpy as np + +from cuml.common.array_descriptor import CumlArrayDescriptor +from cuml.common.doc_utils import generate_docstring +from cuml.internals.array import CumlArray +from cuml.internals.base import Base +from cuml.internals.interop import ( + InteropMixin, + UnsupportedOnGPU, + to_cpu, + to_gpu, +) +from cuml.internals.mixins import FMajorInputTagMixin, RegressorMixin +from cuml.internals.outputs import reflect +from cuml.internals.validation import check_inputs +from cuml.linear_model.base import LinearPredictMixin +from cuml.linear_model.ridge import Ridge + +# The GCV decomposition is computed in the input dtype (float32 stays float32). +# scikit-learn always upcasts to float64 for numerical robustness, so float32 +# results may differ from sklearn by more than float64 results do. This dtype is +# only used for the cheap scalar scoring in the explicit-``cv`` k-fold path, +# where float64 keeps the R^2 computation stable at negligible cost. +_SCORE_DTYPE = cp.float64 + + +def _diag_dot(D, B): + """Compute ``dot(diag(D), B)`` for a 1d ``D`` and 1d/2d ``B``.""" + if B.ndim > 1: + D = D[(slice(None),) + (None,) * (B.ndim - 1)] + return D * B + + +def _decomp_diag(v_prime, Q): + """Compute the diagonal of ``dot(Q, dot(diag(v_prime), Q.T))``.""" + return cp.sum(v_prime * Q**2, axis=1) + + +def _check_gcv_mode(n_samples, n_features, gcv_mode): + """Resolve the ``gcv_mode`` parameter to a concrete strategy.""" + if gcv_mode == "svd": + return "svd" + # "auto" and "eigen" both fall back to a Gram/covariance eigendecomposition, + # picking whichever matrix is smaller. + return "gram" if n_samples <= n_features else "cov" + + +class RidgeCV( + InteropMixin, + RegressorMixin, + LinearPredictMixin, + FMajorInputTagMixin, + Base, +): + """Ridge regression with built-in cross-validation. + + By default, it performs efficient Leave-One-Out Cross-Validation (via + Generalized Cross-Validation) to select the best ``alpha`` from ``alphas``. + When an explicit ``cv`` is provided, a k-fold search over :class:`~cuml.Ridge` + is used instead. + + Parameters + ---------- + alphas : array-like of shape (n_alphas,), default=(0.1, 1.0, 10.0) + Array of alpha values to try. Regularization strength; must be a + positive float. When using Leave-One-Out cross-validation (``cv=None``), + alphas must be strictly positive. + fit_intercept : bool, default=True + Whether to calculate the intercept for this model. If set to false, no + intercept will be used in calculations (i.e. data is expected to be + centered). + scoring : str, callable, default=None + Only ``None`` (negative mean squared error for ``cv=None``, or + :math:`R^2` for an explicit ``cv``) is supported on GPU. A non-``None`` + value raises ``NotImplementedError``. + cv : int or None, default=None + Determines the cross-validation splitting strategy: + + - ``None``, to use efficient Leave-One-Out cross-validation. + - integer, to specify the number of folds. + + gcv_mode : {'auto', 'svd', 'eigen'}, default=None + Flag indicating which strategy to use when performing Leave-One-Out + cross-validation: + + - ``'auto'`` / ``None`` : same as ``'eigen'``. + - ``'svd'`` : use a singular value decomposition of X. + - ``'eigen'`` : use an eigendecomposition of ``X X'`` when + ``n_samples <= n_features`` or ``X' X`` otherwise. + + store_cv_results : bool, default=False + Whether to store the cross-validation values in the ``cv_results_`` + attribute. Only compatible with ``cv=None``. + alpha_per_target : bool, default=False + Whether to optimize the alpha value separately for each target (for + multi-output settings). Only compatible with ``cv=None``. + output_type : {'input', 'array', 'dataframe', 'series', 'df_obj', \ + 'numba', 'cupy', 'numpy', 'cudf', 'pandas'}, default=None + Return results and set estimator attributes to the indicated output + type. If None, the output type set at the module level + (`cuml.global_settings.output_type`) will be used. See + :ref:`output-data-type-configuration` for more info. + verbose : int or boolean, default=False + Sets logging level. It must be one of `cuml.common.logger.level_*`. + See :ref:`verbosity-levels` for more info. + + Attributes + ---------- + cv_results_ : array of shape (n_samples, n_alphas) or \ + (n_samples, n_targets, n_alphas), optional + Cross-validation values for each alpha (only available if + ``store_cv_results=True`` and ``cv=None``). + coef_ : array of shape (n_features,) or (n_targets, n_features) + Weight vector(s). + intercept_ : float or array of shape (n_targets,) + Independent term in the decision function. Set to 0.0 if + ``fit_intercept=False``. + alpha_ : float or array of shape (n_targets,) + Estimated regularization parameter, or, if ``alpha_per_target=True``, + the estimated regularization parameter for each target. + best_score_ : float or array of shape (n_targets,) + Score of the base estimator with the best alpha. + + Notes + ----- + For additional docs, see `scikit-learn's RidgeCV + `_. + + Examples + -------- + >>> from cuml.datasets import make_regression + >>> from cuml.linear_model import RidgeCV + >>> X, y = make_regression(n_samples=50, n_features=5, random_state=0) + >>> reg = RidgeCV(alphas=[0.1, 1.0, 10.0]).fit(X, y) + >>> reg.alpha_ # doctest: +SKIP + 0.1 + """ + + coef_ = CumlArrayDescriptor() + intercept_ = CumlArrayDescriptor() + cv_results_ = CumlArrayDescriptor() + + _cpu_class_path = "sklearn.linear_model.RidgeCV" + + @classmethod + def _get_param_names(cls): + return [ + *super()._get_param_names(), + "alphas", + "fit_intercept", + "scoring", + "cv", + "gcv_mode", + "store_cv_results", + "alpha_per_target", + ] + + @classmethod + def _params_from_cpu(cls, model): + if model.scoring is not None: + raise UnsupportedOnGPU("Custom `scoring` is not supported") + + if model.cv is not None: + # The k-fold path relies on matching scikit-learn's fold splitting + # and scoring exactly, which we don't attempt. Fall back to CPU. + raise UnsupportedOnGPU("`cv != None` is not supported") + + return { + "alphas": model.alphas, + "fit_intercept": model.fit_intercept, + "scoring": model.scoring, + "cv": model.cv, + "gcv_mode": model.gcv_mode, + "store_cv_results": model.store_cv_results, + "alpha_per_target": model.alpha_per_target, + } + + def _params_to_cpu(self): + return { + "alphas": self.alphas, + "fit_intercept": self.fit_intercept, + "scoring": self.scoring, + "cv": self.cv, + "gcv_mode": self.gcv_mode, + "store_cv_results": self.store_cv_results, + "alpha_per_target": self.alpha_per_target, + } + + def _attrs_from_cpu(self, model): + out = { + "coef_": to_gpu(model.coef_), + "intercept_": to_gpu(model.intercept_), + "alpha_": model.alpha_, + "best_score_": model.best_score_, + **super()._attrs_from_cpu(model), + } + if self.store_cv_results and hasattr(model, "cv_results_"): + out["cv_results_"] = to_gpu(model.cv_results_) + return out + + def _attrs_to_cpu(self, model): + out = { + "coef_": to_cpu(self.coef_), + "intercept_": to_cpu(self.intercept_), + "alpha_": self.alpha_, + "best_score_": self.best_score_, + **super()._attrs_to_cpu(model), + } + if self.store_cv_results and hasattr(self, "cv_results_"): + out["cv_results_"] = to_cpu(self.cv_results_) + return out + + def __init__( + self, + alphas=(0.1, 1.0, 10.0), + *, + fit_intercept=True, + scoring=None, + cv=None, + gcv_mode=None, + store_cv_results=False, + alpha_per_target=False, + output_type=None, + verbose=False, + ): + super().__init__(verbose=verbose, output_type=output_type) + self.alphas = alphas + self.fit_intercept = fit_intercept + self.scoring = scoring + self.cv = cv + self.gcv_mode = gcv_mode + self.store_cv_results = store_cv_results + self.alpha_per_target = alpha_per_target + + def __sklearn_tags__(self): + tags = super().__sklearn_tags__() + tags.target_tags.multi_output = True + return tags + + # ------------------------------------------------------------------ + # GCV decomposition/solve helpers (dense-only ports of sklearn's + # ``_RidgeGCV``). All computation happens in float64. + # ------------------------------------------------------------------ + def _solve_eigen_gram( + self, alpha, y, sqrt_sw, eigvals, Q, QT_y, QT_sqrt_sw, XT + ): + w = 1.0 / (eigvals + alpha) + c = Q @ _diag_dot(w, QT_y) + d = _decomp_diag(w, Q) + if self.fit_intercept: + sw_sum = sqrt_sw @ sqrt_sw + Ginv_sqrt_sw = Q @ _diag_dot(w, QT_sqrt_sw) + d -= Ginv_sqrt_sw * sqrt_sw / sw_sum + if y.ndim == 2: + d = d[:, None] + looe = c / d + coef = XT @ c + return looe, coef + + def _solve_eigen_covariance( + self, alpha, y, sqrt_sw, eigvals, V, X, XT_y, XT_sqrt_sw + ): + w = 1.0 / (eigvals + alpha) + Hinv = (V * w) @ V.T + Hinv_XT_y = Hinv @ XT_y + Hinv_XT_sqrt_sw = Hinv @ XT_sqrt_sw + X_Hinv_XT_y = X @ Hinv_XT_y + X_Hinv_XT_sqrt_sw = X @ Hinv_XT_sqrt_sw + alpha_c = y - X_Hinv_XT_y + alpha_d = 1 - cp.sum((X @ Hinv) * X, axis=1) + if self.fit_intercept: + sw_sum = sqrt_sw @ sqrt_sw + alpha_Ginv_sqrt_sw = sqrt_sw - X_Hinv_XT_sqrt_sw + alpha_d -= alpha_Ginv_sqrt_sw * sqrt_sw / sw_sum + if y.ndim == 2: + alpha_d = alpha_d[:, None] + looe = alpha_c / alpha_d + coef = Hinv_XT_y + return looe, coef + + def _solve_svd_design_matrix( + self, alpha, y, sqrt_sw, singvals, U, V, UT_y, UT_sqrt_sw + ): + n_samples, n_features = U.shape[0], V.shape[0] + if n_samples <= n_features: + # Wide X case (n_samples <= n_features). + w = alpha / (singvals**2 + alpha) + alpha_c = U @ _diag_dot(w, UT_y) + alpha_d = _decomp_diag(w, U) + else: + # Long X case (n_features < n_samples). + w = alpha / (singvals**2 + alpha) - 1 + alpha_c = U @ _diag_dot(w, UT_y) + y + alpha_d = _decomp_diag(w, U) + 1 + if self.fit_intercept: + sw_sum = sqrt_sw @ sqrt_sw + if n_samples <= n_features: + alpha_Ginv_sqrt_sw = U @ _diag_dot(w, UT_sqrt_sw) + else: + alpha_Ginv_sqrt_sw = U @ _diag_dot(w, UT_sqrt_sw) + sqrt_sw + alpha_d -= alpha_Ginv_sqrt_sw * sqrt_sw / sw_sum + if y.ndim == 2: + alpha_d = alpha_d[:, None] + looe = alpha_c / alpha_d + coef = V @ _diag_dot(singvals / (singvals**2 + alpha), UT_y) + return looe, coef + + def _fit_gcv(self, X, y, sample_weight, input_dtype): + """Fit via Generalized (leave-one-out) cross-validation.""" + n_samples, n_features = X.shape + # Compute in the input dtype (float32 stays float32). + work = X.dtype + + alphas = np.asarray(self.alphas, dtype=np.float64) + if alphas.ndim == 0: + alphas = alphas.reshape(1) + if (alphas <= 0).any(): + raise ValueError( + "alphas must be strictly positive when cv=None, got " + f"{self.alphas}" + ) + + # Center and rescale following sklearn's `_preprocess_data` + # (`rescale_with_sw=True`). + if self.fit_intercept: + if sample_weight is not None: + sw_sum = sample_weight.sum() + X_offset = (X * sample_weight[:, None]).sum(axis=0) / sw_sum + y_offset = ( + (y * sample_weight[:, None]).sum(axis=0) / sw_sum + if y.ndim == 2 + else (y * sample_weight).sum() / sw_sum + ) + else: + X_offset = X.mean(axis=0) + y_offset = y.mean(axis=0) + X = X - X_offset + y = y - y_offset + else: + X_offset = cp.zeros(n_features, dtype=work) + y_offset = ( + cp.zeros(y.shape[1], dtype=work) + if y.ndim == 2 + else work.type(0.0) + ) + + if sample_weight is not None: + sqrt_sw = cp.sqrt(sample_weight) + X = X * sqrt_sw[:, None] + y = y * (sqrt_sw[:, None] if y.ndim == 2 else sqrt_sw) + else: + sqrt_sw = cp.ones(n_samples, dtype=work) + + gcv_mode = _check_gcv_mode( + n_samples, n_features, self.gcv_mode or "auto" + ) + if gcv_mode == "gram": + K = X @ X.T + eigvals, Q = cp.linalg.eigh(K) + decomposition = (eigvals, Q, Q.T @ y, Q.T @ sqrt_sw, X.T) + solve = self._solve_eigen_gram + elif gcv_mode == "cov": + cov = X.T @ X + eigvals, V = cp.linalg.eigh(cov) + decomposition = (eigvals, V, X, X.T @ y, X.T @ sqrt_sw) + solve = self._solve_eigen_covariance + else: + U, singvals, VT = cp.linalg.svd(X, full_matrices=False) + decomposition = (singvals, U, VT.T, U.T @ y, U.T @ sqrt_sw) + solve = self._solve_svd_design_matrix + + n_y = 1 if y.ndim == 1 else y.shape[1] + n_alphas = len(alphas) + per_target = self.alpha_per_target and n_y > 1 + + if self.store_cv_results: + cv_results = cp.empty((n_samples * n_y, n_alphas), dtype=work) + + best_coef = best_score = best_alpha = None + for i, alpha in enumerate(alphas): + looe, coef = solve(float(alpha), y, sqrt_sw, *decomposition) + squared_errors = looe**2 + if self.store_cv_results: + cv_results[:, i] = squared_errors.reshape(-1) + + if per_target: + score = cp.mean(-squared_errors, axis=0) + else: + score = cp.mean(-squared_errors) + + if best_score is None: + best_coef = coef + if per_target: + best_score = score.reshape(-1) + best_alpha = cp.full(n_y, alpha, dtype=work) + else: + best_score = score + best_alpha = alpha + elif per_target: + to_update = score > best_score + best_coef[:, to_update] = coef[:, to_update] + best_score[to_update] = score[to_update] + best_alpha[to_update] = alpha + elif score > best_score: + best_coef, best_score, best_alpha = coef, score, alpha + + coef = best_coef + if y.ndim == 2: + coef = coef.T + if y.ndim == 1 or y.shape[1] == 1: + coef = coef.reshape(-1) + + # Set intercept, mirroring `LinearModel._set_intercept`. + if self.fit_intercept: + if coef.ndim == 1: + intercept = y_offset - X_offset @ coef + else: + intercept = y_offset - X_offset @ coef.T + else: + intercept = 0.0 + + # Cast fitted attributes back to the input dtype (like sklearn). + self.coef_ = CumlArray(coef.astype(input_dtype)) + if cp.isscalar(intercept) or getattr(intercept, "ndim", 1) == 0: + self.intercept_ = float(intercept) + else: + self.intercept_ = CumlArray( + cp.asarray(intercept, dtype=input_dtype) + ) + # `alpha_`/`best_score_` are small metadata; keep them on host to match + # scikit-learn (a float, or a numpy array when alpha_per_target=True). + if per_target: + self.alpha_ = cp.asnumpy(cp.asarray(best_alpha)).astype( + input_dtype + ) + self.best_score_ = cp.asnumpy(cp.asarray(best_score)).astype( + input_dtype + ) + else: + self.alpha_ = float(best_alpha) + self.best_score_ = float(best_score) + + if self.store_cv_results: + if y.ndim == 1: + shape = (n_samples, n_alphas) + else: + shape = (n_samples, n_y, n_alphas) + self.cv_results_ = CumlArray( + cv_results.reshape(shape).astype(input_dtype) + ) + + def _fit_cv(self, X, y, sample_weight, input_dtype): + """Fit via an explicit k-fold search over :class:`~cuml.Ridge`.""" + from cuml.model_selection import KFold + + if self.store_cv_results: + raise ValueError( + "cv!=None and store_cv_results=True are incompatible" + ) + if self.alpha_per_target: + raise ValueError( + "cv!=None and alpha_per_target=True are incompatible" + ) + + alphas = np.asarray(self.alphas, dtype=np.float64) + if alphas.ndim == 0: + alphas = alphas.reshape(1) + if (alphas < 0).any(): + raise ValueError(f"alphas must be non-negative, got {self.alphas}") + + if isinstance(self.cv, int): + splitter = KFold(n_splits=self.cv) + splits = list(splitter.split(X, y)) + elif hasattr(self.cv, "split"): + splits = list(self.cv.split(X, y)) + else: + splits = list(self.cv) + + mean_scores = cp.empty(len(alphas), dtype=_SCORE_DTYPE) + for i, alpha in enumerate(alphas): + fold_scores = [] + for train, test in splits: + train = cp.asarray(train) + test = cp.asarray(test) + model = Ridge( + alpha=float(alpha), + fit_intercept=self.fit_intercept, + output_type="cupy", + ) + if sample_weight is None: + model.fit(X[train], y[train]) + else: + model.fit( + X[train], y[train], sample_weight=sample_weight[train] + ) + y_pred = model.predict(X[test]) + fold_scores.append(_r2_score(y[test], y_pred)) + mean_scores[i] = sum(fold_scores) / len(fold_scores) + + best_idx = int(cp.argmax(mean_scores)) + best_alpha = float(alphas[best_idx]) + + model = Ridge( + alpha=best_alpha, + fit_intercept=self.fit_intercept, + output_type="cupy", + ) + if sample_weight is None: + model.fit(X, y) + else: + model.fit(X, y, sample_weight=sample_weight) + + self.coef_ = CumlArray(cp.asarray(model.coef_, dtype=input_dtype)) + intercept = model.intercept_ + if isinstance(intercept, CumlArray): + intercept = intercept.to_output("cupy") + if cp.isscalar(intercept) or getattr(intercept, "ndim", 1) == 0: + self.intercept_ = float(intercept) + else: + self.intercept_ = CumlArray( + cp.asarray(intercept, dtype=input_dtype) + ) + self.alpha_ = best_alpha + self.best_score_ = float(mean_scores[best_idx]) + + @generate_docstring() + @reflect(reset=True) + def fit(self, X, y, sample_weight=None) -> "RidgeCV": + """Fit the RidgeCV model with X and y.""" + if self.scoring is not None: + raise NotImplementedError( + "Only `scoring=None` is supported on GPU" + ) + + X, y, sample_weight = check_inputs( + self, + X, + y, + sample_weight, + dtype=("float32", "float64"), + accept_sparse=False, + accept_multi_output=True, + ensure_min_samples=2, + reset=True, + ) + input_dtype = X.dtype + + # Keep the input dtype (float32 stays float32) rather than upcasting. + X = cp.asarray(X) + y = cp.asarray(y) + if sample_weight is not None: + sample_weight = cp.asarray(sample_weight) + + if self.cv is None: + self._fit_gcv(X, y, sample_weight, input_dtype) + else: + self._fit_cv(X, y, sample_weight, input_dtype) + + return self + + +def _r2_score(y_true, y_pred): + """Uniform-average :math:`R^2` matching sklearn's default scorer.""" + y_true = cp.asarray(y_true, dtype=_SCORE_DTYPE) + y_pred = cp.asarray(y_pred, dtype=_SCORE_DTYPE) + ss_res = ((y_true - y_pred) ** 2).sum(axis=0) + ss_tot = ((y_true - y_true.mean(axis=0)) ** 2).sum(axis=0) + r2 = 1.0 - ss_res / ss_tot + return float(cp.mean(r2)) diff --git a/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py b/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py new file mode 100644 index 0000000000..d9cf241d48 --- /dev/null +++ b/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py @@ -0,0 +1,74 @@ +# +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. +# SPDX-License-Identifier: Apache-2.0 +# + +import numpy as np +import pytest +from sklearn.datasets import make_regression +from sklearn.linear_model import RidgeCV +from sklearn.metrics import r2_score +from sklearn.preprocessing import StandardScaler + +ALPHAS = [0.1, 1.0, 10.0, 100.0] + + +@pytest.fixture(scope="module") +def regression_data(): + X, y = make_regression( + n_samples=500, + n_features=20, + n_informative=10, + noise=0.1, + random_state=42, + ) + X = StandardScaler().fit_transform(X) + return X, y + + +@pytest.mark.parametrize("fit_intercept", [True, False]) +def test_ridge_cv_gcv_runs_on_gpu(regression_data, fit_intercept): + X, y = regression_data + model = RidgeCV(alphas=ALPHAS, fit_intercept=fit_intercept) + model.fit(X, y) + # cv=None (GCV) is accelerated on GPU + assert model._gpu is not None + y_pred = model.predict(X) + assert r2_score(y, y_pred) > 0.5 + + +def test_ridge_cv_selects_reasonable_alpha(regression_data): + X, y = regression_data + model = RidgeCV(alphas=ALPHAS).fit(X, y) + assert model._gpu is not None + assert model.alpha_ in ALPHAS + + +def test_ridge_cv_explicit_cv_falls_back_to_cpu(regression_data): + X, y = regression_data + model = RidgeCV(alphas=ALPHAS, cv=5) + model.fit(X, y) + # An explicit cv is not supported on GPU, should fall back to CPU + assert model._gpu is None + y_pred = model.predict(X) + assert r2_score(y, y_pred) > 0.5 + + +def test_ridge_cv_custom_scoring_falls_back_to_cpu(regression_data): + X, y = regression_data + model = RidgeCV(alphas=ALPHAS, scoring="neg_mean_squared_error") + model.fit(X, y) + # Custom scoring is not supported on GPU, should fall back to CPU + assert model._gpu is None + y_pred = model.predict(X) + assert r2_score(y, y_pred) > 0.5 + + +def test_ridge_cv_multi_output_runs_on_gpu(): + X, y = make_regression( + n_samples=400, n_features=15, n_targets=3, noise=0.1, random_state=0 + ) + model = RidgeCV(alphas=ALPHAS).fit(X, y) + assert model._gpu is not None + assert model.coef_.shape == (3, 15) + np.testing.assert_array_equal(np.asarray(model.predict(X)).shape, (400, 3)) diff --git a/python/cuml/tests/test_ridge_cv.py b/python/cuml/tests/test_ridge_cv.py new file mode 100644 index 0000000000..46dd13b392 --- /dev/null +++ b/python/cuml/tests/test_ridge_cv.py @@ -0,0 +1,191 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. +# SPDX-License-Identifier: Apache-2.0 + +import numpy as np +import pytest +import sklearn.linear_model +from sklearn.datasets import make_regression + +import cuml + +ALPHAS = [0.01, 0.1, 1.0, 10.0, 100.0] + + +def _make(n_samples=200, n_features=20, n_targets=1, dtype=np.float64): + X, y = make_regression( + n_samples=n_samples, + n_features=n_features, + n_informative=max(n_features // 2, 1), + n_targets=n_targets, + noise=1.0, + random_state=42, + ) + return X.astype(dtype), y.astype(dtype) + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("gcv_mode", [None, "auto", "svd", "eigen"]) +@pytest.mark.parametrize("fit_intercept", [True, False]) +@pytest.mark.parametrize( + "n_samples, n_features", + [(200, 20), (40, 80)], # n > p and n < p +) +def test_gcv_matches_sklearn( + dtype, gcv_mode, fit_intercept, n_samples, n_features +): + X, y = _make(n_samples=n_samples, n_features=n_features, dtype=dtype) + + cu = cuml.RidgeCV( + alphas=ALPHAS, fit_intercept=fit_intercept, gcv_mode=gcv_mode + ).fit(X, y) + sk = sklearn.linear_model.RidgeCV( + alphas=ALPHAS, fit_intercept=fit_intercept, gcv_mode=gcv_mode + ).fit(X, y) + + # cuML computes the GCV in the input dtype, while sklearn always upcasts to + # float64. In float32 the forward error is ~kappa_eff * eps_f32, so we use + # per-quantity tolerances (~10x over measured worst case). Note `coef_` is + # bound by `atol`, not `rtol`: make_regression leaves ~half the true + # coefficients near zero, where a normal float32 absolute error is a large + # relative error. + if dtype == np.float32: + coef_tol = dict(atol=1e-3, rtol=1e-3) + intercept_tol = dict(atol=1e-3) + predict_tol = dict(atol=5e-3, rtol=1e-3) + best_score_tol = dict(rtol=1e-2) + else: + coef_tol = dict(atol=1e-6, rtol=1e-5) + intercept_tol = dict(atol=1e-6, rtol=1e-5) + predict_tol = dict(atol=1e-5, rtol=1e-4) + best_score_tol = dict(rtol=1e-5) + + assert cu.alpha_ == pytest.approx(sk.alpha_) + np.testing.assert_allclose(cu.coef_, sk.coef_, **coef_tol) + np.testing.assert_allclose(cu.intercept_, sk.intercept_, **intercept_tol) + np.testing.assert_allclose( + cu.best_score_, sk.best_score_, **best_score_tol + ) + np.testing.assert_allclose(cu.predict(X), sk.predict(X), **predict_tol) + + +@pytest.mark.parametrize("fit_intercept", [True, False]) +def test_gcv_sample_weight_matches_sklearn(fit_intercept): + X, y = _make(n_samples=200, n_features=20) + sw = np.random.default_rng(0).uniform(0.5, 2.0, size=X.shape[0]) + + cu = cuml.RidgeCV(alphas=ALPHAS, fit_intercept=fit_intercept).fit( + X, y, sample_weight=sw + ) + sk = sklearn.linear_model.RidgeCV( + alphas=ALPHAS, fit_intercept=fit_intercept + ).fit(X, y, sample_weight=sw) + + assert cu.alpha_ == pytest.approx(sk.alpha_) + np.testing.assert_allclose(cu.coef_, sk.coef_, atol=1e-6, rtol=1e-5) + np.testing.assert_allclose( + cu.predict(X), sk.predict(X), atol=1e-5, rtol=1e-4 + ) + + +@pytest.mark.parametrize("alpha_per_target", [False, True]) +def test_gcv_multi_output_matches_sklearn(alpha_per_target): + X, y = _make(n_samples=300, n_features=20, n_targets=3) + + cu = cuml.RidgeCV(alphas=ALPHAS, alpha_per_target=alpha_per_target).fit( + X, y + ) + sk = sklearn.linear_model.RidgeCV( + alphas=ALPHAS, alpha_per_target=alpha_per_target + ).fit(X, y) + + np.testing.assert_allclose(cu.alpha_, sk.alpha_) + assert cu.coef_.shape == (3, 20) + np.testing.assert_allclose(cu.coef_, sk.coef_, atol=1e-6, rtol=1e-5) + np.testing.assert_allclose( + cu.predict(X), sk.predict(X), atol=1e-5, rtol=1e-4 + ) + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +def test_gcv_dtype_and_shapes(dtype): + X, y = _make(n_samples=120, n_features=10, n_targets=3, dtype=dtype) + + # 1D y + m = cuml.RidgeCV(alphas=ALPHAS).fit(X, y[:, 0]) + assert m.coef_.shape == (10,) + assert m.coef_.dtype == dtype + assert np.isscalar(m.intercept_) or m.intercept_.shape == () + + # multi-output y + m = cuml.RidgeCV(alphas=ALPHAS).fit(X, y) + assert m.coef_.shape == (3, 10) + assert m.coef_.dtype == dtype + assert m.intercept_.shape == (3,) + + # no intercept + m = cuml.RidgeCV(alphas=ALPHAS, fit_intercept=False).fit(X, y[:, 0]) + assert m.intercept_ == 0.0 + + +def test_store_cv_results_shape(): + X, y = _make(n_samples=150, n_features=12, n_targets=1) + m = cuml.RidgeCV(alphas=ALPHAS, store_cv_results=True).fit(X, y) + assert m.cv_results_.shape == (150, len(ALPHAS)) + + X, y = _make(n_samples=150, n_features=12, n_targets=3) + m = cuml.RidgeCV(alphas=ALPHAS, store_cv_results=True).fit(X, y) + assert m.cv_results_.shape == (150, 3, len(ALPHAS)) + + +@pytest.mark.parametrize("cv", [3, 5]) +@pytest.mark.parametrize("fit_intercept", [True, False]) +def test_kfold_path_matches_sklearn(cv, fit_intercept): + X, y = _make(n_samples=300, n_features=15) + + cu = cuml.RidgeCV(alphas=ALPHAS, cv=cv, fit_intercept=fit_intercept).fit( + X, y + ) + sk = sklearn.linear_model.RidgeCV( + alphas=ALPHAS, cv=cv, fit_intercept=fit_intercept + ).fit(X, y) + + assert cu.alpha_ == pytest.approx(sk.alpha_) + np.testing.assert_allclose( + cu.predict(X), sk.predict(X), atol=1e-2, rtol=1e-2 + ) + + +def test_invalid_alphas_gcv(): + X, y = _make() + + est = cuml.RidgeCV(alphas=[0.0, 1.0]) + with pytest.raises(ValueError, match="strictly positive"): + est.fit(X, y) + + +def test_incompatible_options_with_cv(): + X, y = _make() + with pytest.raises(ValueError, match="store_cv_results"): + cuml.RidgeCV(alphas=ALPHAS, cv=5, store_cv_results=True).fit(X, y) + with pytest.raises(ValueError, match="alpha_per_target"): + cuml.RidgeCV(alphas=ALPHAS, cv=5, alpha_per_target=True).fit(X, y) + + +def test_custom_scoring_not_implemented(): + X, y = _make() + + est = cuml.RidgeCV(alphas=ALPHAS, scoring="r2") + with pytest.raises(NotImplementedError): + est.fit(X, y) + + +def test_sklearn_roundtrip(): + X, y = _make(n_samples=200, n_features=20) + cu = cuml.RidgeCV(alphas=ALPHAS).fit(X, y) + + sk = cu.as_sklearn() + assert isinstance(sk, sklearn.linear_model.RidgeCV) + np.testing.assert_allclose(cu.predict(X), sk.predict(X), rtol=1e-5) + + cu2 = cuml.RidgeCV.from_sklearn(sk) + np.testing.assert_allclose(cu2.predict(X), sk.predict(X), rtol=1e-4) From 02bfdca825b443dd65559ddae90c0dc534a8f9bd Mon Sep 17 00:00:00 2001 From: Tim Head Date: Wed, 8 Jul 2026 13:01:22 +0200 Subject: [PATCH 2/5] Tidy up --- python/cuml/cuml/linear_model/ridge_cv.py | 148 ++++++++-------------- 1 file changed, 52 insertions(+), 96 deletions(-) diff --git a/python/cuml/cuml/linear_model/ridge_cv.py b/python/cuml/cuml/linear_model/ridge_cv.py index 0fc3e49a2d..ce37301e77 100644 --- a/python/cuml/cuml/linear_model/ridge_cv.py +++ b/python/cuml/cuml/linear_model/ridge_cv.py @@ -20,13 +20,7 @@ from cuml.internals.validation import check_inputs from cuml.linear_model.base import LinearPredictMixin from cuml.linear_model.ridge import Ridge - -# The GCV decomposition is computed in the input dtype (float32 stays float32). -# scikit-learn always upcasts to float64 for numerical robustness, so float32 -# results may differ from sklearn by more than float64 results do. This dtype is -# only used for the cheap scalar scoring in the explicit-``cv`` k-fold path, -# where float64 keeps the R^2 computation stable at negligible cost. -_SCORE_DTYPE = cp.float64 +from cuml.model_selection import KFold def _diag_dot(D, B): @@ -41,15 +35,6 @@ def _decomp_diag(v_prime, Q): return cp.sum(v_prime * Q**2, axis=1) -def _check_gcv_mode(n_samples, n_features, gcv_mode): - """Resolve the ``gcv_mode`` parameter to a concrete strategy.""" - if gcv_mode == "svd": - return "svd" - # "auto" and "eigen" both fall back to a Gram/covariance eigendecomposition, - # picking whichever matrix is smaller. - return "gram" if n_samples <= n_features else "cov" - - class RidgeCV( InteropMixin, RegressorMixin, @@ -126,11 +111,6 @@ class RidgeCV( best_score_ : float or array of shape (n_targets,) Score of the base estimator with the best alpha. - Notes - ----- - For additional docs, see `scikit-learn's RidgeCV - `_. - Examples -------- >>> from cuml.datasets import make_regression @@ -144,6 +124,8 @@ class RidgeCV( coef_ = CumlArrayDescriptor() intercept_ = CumlArrayDescriptor() cv_results_ = CumlArrayDescriptor() + alpha_ = CumlArrayDescriptor() + best_score_ = CumlArrayDescriptor() _cpu_class_path = "sklearn.linear_model.RidgeCV" @@ -166,8 +148,10 @@ def _params_from_cpu(cls, model): raise UnsupportedOnGPU("Custom `scoring` is not supported") if model.cv is not None: - # The k-fold path relies on matching scikit-learn's fold splitting - # and scoring exactly, which we don't attempt. Fall back to CPU. + # Only the GCV (cv=None) path is accelerated. cuML's k-fold path uses + # its own splitter/scoring and won't match sklearn's fold selection + # exactly, so under cuml.accel we fall back to CPU to preserve + # sklearn-identical results. raise UnsupportedOnGPU("`cv != None` is not supported") return { @@ -195,8 +179,8 @@ def _attrs_from_cpu(self, model): out = { "coef_": to_gpu(model.coef_), "intercept_": to_gpu(model.intercept_), - "alpha_": model.alpha_, - "best_score_": model.best_score_, + "alpha_": to_gpu(model.alpha_), + "best_score_": to_gpu(model.best_score_), **super()._attrs_from_cpu(model), } if self.store_cv_results and hasattr(model, "cv_results_"): @@ -207,8 +191,8 @@ def _attrs_to_cpu(self, model): out = { "coef_": to_cpu(self.coef_), "intercept_": to_cpu(self.intercept_), - "alpha_": self.alpha_, - "best_score_": self.best_score_, + "alpha_": to_cpu(self.alpha_), + "best_score_": to_cpu(self.best_score_), **super()._attrs_to_cpu(model), } if self.store_cv_results and hasattr(self, "cv_results_"): @@ -244,7 +228,7 @@ def __sklearn_tags__(self): # ------------------------------------------------------------------ # GCV decomposition/solve helpers (dense-only ports of sklearn's - # ``_RidgeGCV``). All computation happens in float64. + # ``_RidgeGCV``). Computation happens in the input dtype. # ------------------------------------------------------------------ def _solve_eigen_gram( self, alpha, y, sqrt_sw, eigvals, Q, QT_y, QT_sqrt_sw, XT @@ -310,10 +294,12 @@ def _solve_svd_design_matrix( coef = V @ _diag_dot(singvals / (singvals**2 + alpha), UT_y) return looe, coef - def _fit_gcv(self, X, y, sample_weight, input_dtype): + def _fit_gcv(self, X, y, sample_weight): """Fit via Generalized (leave-one-out) cross-validation.""" n_samples, n_features = X.shape - # Compute in the input dtype (float32 stays float32). + # Compute in the input dtype (float32 stays float32). scikit-learn + # always upcasts to float64 for numerical robustness, so float32 results + # may differ from sklearn by more than float64 results do. work = X.dtype alphas = np.asarray(self.alphas, dtype=np.float64) @@ -356,9 +342,13 @@ def _fit_gcv(self, X, y, sample_weight, input_dtype): else: sqrt_sw = cp.ones(n_samples, dtype=work) - gcv_mode = _check_gcv_mode( - n_samples, n_features, self.gcv_mode or "auto" - ) + # "auto"/"eigen" pick the smaller of the Gram (X X') or covariance + # (X' X) eigendecomposition; "svd" decomposes X directly. + if (self.gcv_mode or "auto") == "svd": + gcv_mode = "svd" + else: + gcv_mode = "gram" if n_samples <= n_features else "cov" + if gcv_mode == "gram": K = X @ X.T eigvals, Q = cp.linalg.eigh(K) @@ -424,23 +414,16 @@ def _fit_gcv(self, X, y, sample_weight, input_dtype): else: intercept = 0.0 - # Cast fitted attributes back to the input dtype (like sklearn). - self.coef_ = CumlArray(coef.astype(input_dtype)) - if cp.isscalar(intercept) or getattr(intercept, "ndim", 1) == 0: + self.coef_ = CumlArray(coef) + if cp.isscalar(intercept) or intercept.ndim == 0: self.intercept_ = float(intercept) else: - self.intercept_ = CumlArray( - cp.asarray(intercept, dtype=input_dtype) - ) - # `alpha_`/`best_score_` are small metadata; keep them on host to match - # scikit-learn (a float, or a numpy array when alpha_per_target=True). + self.intercept_ = CumlArray(intercept) + # A single float when selecting one alpha, or a per-target array when + # `alpha_per_target=True`. if per_target: - self.alpha_ = cp.asnumpy(cp.asarray(best_alpha)).astype( - input_dtype - ) - self.best_score_ = cp.asnumpy(cp.asarray(best_score)).astype( - input_dtype - ) + self.alpha_ = CumlArray(best_alpha) + self.best_score_ = CumlArray(best_score) else: self.alpha_ = float(best_alpha) self.best_score_ = float(best_score) @@ -450,13 +433,10 @@ def _fit_gcv(self, X, y, sample_weight, input_dtype): shape = (n_samples, n_alphas) else: shape = (n_samples, n_y, n_alphas) - self.cv_results_ = CumlArray( - cv_results.reshape(shape).astype(input_dtype) - ) + self.cv_results_ = CumlArray(cv_results.reshape(shape)) - def _fit_cv(self, X, y, sample_weight, input_dtype): + def _fit_cv(self, X, y, sample_weight): """Fit via an explicit k-fold search over :class:`~cuml.Ridge`.""" - from cuml.model_selection import KFold if self.store_cv_results: raise ValueError( @@ -481,8 +461,11 @@ def _fit_cv(self, X, y, sample_weight, input_dtype): else: splits = list(self.cv) - mean_scores = cp.empty(len(alphas), dtype=_SCORE_DTYPE) - for i, alpha in enumerate(alphas): + # Mean cross-validation score per alpha. With scoring=None, sklearn's + # RidgeCV scores each fold with the estimator's default (uniform-average + # R^2), which is exactly what `cuml.Ridge.score` computes. + mean_scores = [] + for alpha in alphas: fold_scores = [] for train, test in splits: train = cp.asarray(train) @@ -492,17 +475,15 @@ def _fit_cv(self, X, y, sample_weight, input_dtype): fit_intercept=self.fit_intercept, output_type="cupy", ) - if sample_weight is None: - model.fit(X[train], y[train]) - else: - model.fit( - X[train], y[train], sample_weight=sample_weight[train] - ) - y_pred = model.predict(X[test]) - fold_scores.append(_r2_score(y[test], y_pred)) - mean_scores[i] = sum(fold_scores) / len(fold_scores) - - best_idx = int(cp.argmax(mean_scores)) + sw_train = ( + None if sample_weight is None else sample_weight[train] + ) + model.fit(X[train], y[train], sample_weight=sw_train) + # Unweighted score to match scikit-learn's behaviour + fold_scores.append(model.score(X[test], y[test])) + mean_scores.append(sum(fold_scores) / len(fold_scores)) + + best_idx = int(np.argmax(mean_scores)) best_alpha = float(alphas[best_idx]) model = Ridge( @@ -510,21 +491,14 @@ def _fit_cv(self, X, y, sample_weight, input_dtype): fit_intercept=self.fit_intercept, output_type="cupy", ) - if sample_weight is None: - model.fit(X, y) - else: - model.fit(X, y, sample_weight=sample_weight) + model.fit(X, y, sample_weight=sample_weight) - self.coef_ = CumlArray(cp.asarray(model.coef_, dtype=input_dtype)) + self.coef_ = model.coef_ intercept = model.intercept_ - if isinstance(intercept, CumlArray): - intercept = intercept.to_output("cupy") - if cp.isscalar(intercept) or getattr(intercept, "ndim", 1) == 0: + if cp.isscalar(intercept) or intercept.ndim == 0: self.intercept_ = float(intercept) else: - self.intercept_ = CumlArray( - cp.asarray(intercept, dtype=input_dtype) - ) + self.intercept_ = intercept self.alpha_ = best_alpha self.best_score_ = float(mean_scores[best_idx]) @@ -548,27 +522,9 @@ def fit(self, X, y, sample_weight=None) -> "RidgeCV": ensure_min_samples=2, reset=True, ) - input_dtype = X.dtype - - # Keep the input dtype (float32 stays float32) rather than upcasting. - X = cp.asarray(X) - y = cp.asarray(y) - if sample_weight is not None: - sample_weight = cp.asarray(sample_weight) - if self.cv is None: - self._fit_gcv(X, y, sample_weight, input_dtype) + self._fit_gcv(X, y, sample_weight) else: - self._fit_cv(X, y, sample_weight, input_dtype) + self._fit_cv(X, y, sample_weight) return self - - -def _r2_score(y_true, y_pred): - """Uniform-average :math:`R^2` matching sklearn's default scorer.""" - y_true = cp.asarray(y_true, dtype=_SCORE_DTYPE) - y_pred = cp.asarray(y_pred, dtype=_SCORE_DTYPE) - ss_res = ((y_true - y_pred) ** 2).sum(axis=0) - ss_tot = ((y_true - y_true.mean(axis=0)) ** 2).sum(axis=0) - r2 = 1.0 - ss_res / ss_tot - return float(cp.mean(r2)) From be81c4843b9cbf5353d333636f7882975d04cd24 Mon Sep 17 00:00:00 2001 From: Tim Head Date: Wed, 8 Jul 2026 13:20:46 +0200 Subject: [PATCH 3/5] Improve test --- .../cuml/cuml_accel_tests/integration/test_ridge_cv.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py b/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py index d9cf241d48..20a4082532 100644 --- a/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py +++ b/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py @@ -37,11 +37,16 @@ def test_ridge_cv_gcv_runs_on_gpu(regression_data, fit_intercept): assert r2_score(y, y_pred) > 0.5 -def test_ridge_cv_selects_reasonable_alpha(regression_data): +def test_ridge_cv_selects_best_alpha(regression_data): X, y = regression_data model = RidgeCV(alphas=ALPHAS).fit(X, y) assert model._gpu is not None - assert model.alpha_ in ALPHAS + + # The grid fit should pick the alpha that maximizes the CV score, i.e. the + # same one we'd get by scoring each alpha on its own. + scores = {a: RidgeCV(alphas=[a]).fit(X, y).best_score_ for a in ALPHAS} + assert model.alpha_ == max(scores, key=scores.get) + assert model.best_score_ == pytest.approx(scores[model.alpha_]) def test_ridge_cv_explicit_cv_falls_back_to_cpu(regression_data): From 3d875b7b0639e0e82a06110561b74dfebccbf34d Mon Sep 17 00:00:00 2001 From: Tim Head Date: Thu, 9 Jul 2026 11:14:48 +0200 Subject: [PATCH 4/5] Add RidgeCV to common checks run --- python/cuml/tests/test_sklearn_compatibility.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/cuml/tests/test_sklearn_compatibility.py b/python/cuml/tests/test_sklearn_compatibility.py index 3070aeb0f9..07d683f687 100644 --- a/python/cuml/tests/test_sklearn_compatibility.py +++ b/python/cuml/tests/test_sklearn_compatibility.py @@ -28,6 +28,7 @@ MBSGDClassifier, MBSGDRegressor, Ridge, + RidgeCV, ) from cuml.manifold import TSNE, UMAP, SpectralEmbedding from cuml.multiclass import OneVsOneClassifier, OneVsRestClassifier @@ -106,6 +107,7 @@ LedoitWolf(), Lars(), Ridge(), + RidgeCV(), ElasticNet(), Lasso(), LinearRegression(), From bd5118e9244d5fd1049964b6fcb380e67e3147db Mon Sep 17 00:00:00 2001 From: Tim Head Date: Mon, 3 Aug 2026 09:54:33 +0200 Subject: [PATCH 5/5] Refresh --- python/cuml/cuml/linear_model/ridge_cv.py | 79 ++++++++++++++--------- 1 file changed, 48 insertions(+), 31 deletions(-) diff --git a/python/cuml/cuml/linear_model/ridge_cv.py b/python/cuml/cuml/linear_model/ridge_cv.py index ce37301e77..a92a629007 100644 --- a/python/cuml/cuml/linear_model/ridge_cv.py +++ b/python/cuml/cuml/linear_model/ridge_cv.py @@ -1,22 +1,15 @@ # -# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 # import cupy as cp import numpy as np -from cuml.common.array_descriptor import CumlArrayDescriptor from cuml.common.doc_utils import generate_docstring -from cuml.internals.array import CumlArray from cuml.internals.base import Base -from cuml.internals.interop import ( - InteropMixin, - UnsupportedOnGPU, - to_cpu, - to_gpu, -) +from cuml.internals.interop import InteropMixin, UnsupportedOnGPU from cuml.internals.mixins import FMajorInputTagMixin, RegressorMixin -from cuml.internals.outputs import reflect +from cuml.internals.outputs import ReflectedAttr, mlfunc from cuml.internals.validation import check_inputs from cuml.linear_model.base import LinearPredictMixin from cuml.linear_model.ridge import Ridge @@ -121,11 +114,11 @@ class RidgeCV( 0.1 """ - coef_ = CumlArrayDescriptor() - intercept_ = CumlArrayDescriptor() - cv_results_ = CumlArrayDescriptor() - alpha_ = CumlArrayDescriptor() - best_score_ = CumlArrayDescriptor() + coef_ = ReflectedAttr() + intercept_ = ReflectedAttr() + cv_results_ = ReflectedAttr() + alpha_ = ReflectedAttr() + best_score_ = ReflectedAttr() _cpu_class_path = "sklearn.linear_model.RidgeCV" @@ -177,26 +170,50 @@ def _params_to_cpu(self): def _attrs_from_cpu(self, model): out = { - "coef_": to_gpu(model.coef_), - "intercept_": to_gpu(model.intercept_), - "alpha_": to_gpu(model.alpha_), - "best_score_": to_gpu(model.best_score_), + "coef_": cp.asarray(model.coef_), + "intercept_": ( + model.intercept_ + if cp.isscalar(model.intercept_) + else cp.asarray(model.intercept_) + ), + "alpha_": ( + model.alpha_ + if cp.isscalar(model.alpha_) + else cp.asarray(model.alpha_) + ), + "best_score_": ( + model.best_score_ + if cp.isscalar(model.best_score_) + else cp.asarray(model.best_score_) + ), **super()._attrs_from_cpu(model), } if self.store_cv_results and hasattr(model, "cv_results_"): - out["cv_results_"] = to_gpu(model.cv_results_) + out["cv_results_"] = cp.asarray(model.cv_results_) return out def _attrs_to_cpu(self, model): out = { - "coef_": to_cpu(self.coef_), - "intercept_": to_cpu(self.intercept_), - "alpha_": to_cpu(self.alpha_), - "best_score_": to_cpu(self.best_score_), + "coef_": cp.asnumpy(self.coef_), + "intercept_": ( + self.intercept_ + if cp.isscalar(self.intercept_) + else cp.asnumpy(self.intercept_) + ), + "alpha_": ( + self.alpha_ + if cp.isscalar(self.alpha_) + else cp.asnumpy(self.alpha_) + ), + "best_score_": ( + self.best_score_ + if cp.isscalar(self.best_score_) + else cp.asnumpy(self.best_score_) + ), **super()._attrs_to_cpu(model), } if self.store_cv_results and hasattr(self, "cv_results_"): - out["cv_results_"] = to_cpu(self.cv_results_) + out["cv_results_"] = cp.asnumpy(self.cv_results_) return out def __init__( @@ -414,16 +431,16 @@ def _fit_gcv(self, X, y, sample_weight): else: intercept = 0.0 - self.coef_ = CumlArray(coef) + self.coef_ = coef if cp.isscalar(intercept) or intercept.ndim == 0: self.intercept_ = float(intercept) else: - self.intercept_ = CumlArray(intercept) + self.intercept_ = intercept # A single float when selecting one alpha, or a per-target array when # `alpha_per_target=True`. if per_target: - self.alpha_ = CumlArray(best_alpha) - self.best_score_ = CumlArray(best_score) + self.alpha_ = best_alpha + self.best_score_ = best_score else: self.alpha_ = float(best_alpha) self.best_score_ = float(best_score) @@ -433,7 +450,7 @@ def _fit_gcv(self, X, y, sample_weight): shape = (n_samples, n_alphas) else: shape = (n_samples, n_y, n_alphas) - self.cv_results_ = CumlArray(cv_results.reshape(shape)) + self.cv_results_ = cv_results.reshape(shape) def _fit_cv(self, X, y, sample_weight): """Fit via an explicit k-fold search over :class:`~cuml.Ridge`.""" @@ -503,7 +520,7 @@ def _fit_cv(self, X, y, sample_weight): self.best_score_ = float(mean_scores[best_idx]) @generate_docstring() - @reflect(reset=True) + @mlfunc(set_input_type=True) def fit(self, X, y, sample_weight=None) -> "RidgeCV": """Fit the RidgeCV model with X and y.""" if self.scoring is not None: