diff --git a/docs/source/api/cuml.linear_model.rst b/docs/source/api/cuml.linear_model.rst index 00f9ab2bcc..dfd455a449 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 552d66b7bb..96698fd683 100644 --- a/docs/source/api/index.rst +++ b/docs/source/api/index.rst @@ -294,6 +294,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/compatibility.rst b/docs/source/cuml-accel/compatibility.rst index ed3c267ea7..5082811f73 100644 --- a/docs/source/cuml-accel/compatibility.rst +++ b/docs/source/cuml-accel/compatibility.rst @@ -301,6 +301,16 @@ scores like ``sklearn.metrics.r2_score`` (for regression) or - If ``positive=True`` or ``solver="lbfgs"``. +.. dropdown:: ``RidgeCV`` + :name: 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. + + .. dropdown:: ``Lasso`` :name: lasso diff --git a/python/cuml/cuml/__init__.py b/python/cuml/cuml/__init__.py index 40e201e969..bbd6ad65de 100644 --- a/python/cuml/cuml/__init__.py +++ b/python/cuml/cuml/__init__.py @@ -50,6 +50,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 @@ -155,6 +156,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 75db6f6733..f63cc90967 100644 --- a/python/cuml/cuml/accel/_overrides/sklearn/linear_model.py +++ b/python/cuml/cuml/accel/_overrides/sklearn/linear_model.py @@ -18,6 +18,7 @@ "LogisticRegression", "ElasticNet", "Ridge", + "RidgeCV", "Lasso", ) @@ -63,6 +64,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 a8fddd0a94..a5607c809f 100644 --- a/python/cuml/cuml/benchmark/algorithms.py +++ b/python/cuml/cuml/benchmark/algorithms.py @@ -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 # import functools @@ -389,6 +389,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 @@ -423,6 +424,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 @@ -575,6 +577,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 702272a1e2..bba052777d 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 eb320c7342..0317c70f72 100644 --- a/python/cuml/cuml/linear_model/__init__.py +++ b/python/cuml/cuml/linear_model/__init__.py @@ -11,3 +11,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..a92a629007 --- /dev/null +++ b/python/cuml/cuml/linear_model/ridge_cv.py @@ -0,0 +1,547 @@ +# +# 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.doc_utils import generate_docstring +from cuml.internals.base import Base +from cuml.internals.interop import InteropMixin, UnsupportedOnGPU +from cuml.internals.mixins import FMajorInputTagMixin, RegressorMixin +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 +from cuml.model_selection import KFold + + +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) + + +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. + + 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_ = ReflectedAttr() + intercept_ = ReflectedAttr() + cv_results_ = ReflectedAttr() + alpha_ = ReflectedAttr() + best_score_ = ReflectedAttr() + + _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: + # 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 { + "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_": 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_"] = cp.asarray(model.cv_results_) + return out + + def _attrs_to_cpu(self, model): + out = { + "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_"] = cp.asnumpy(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``). Computation happens in the input dtype. + # ------------------------------------------------------------------ + 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): + """Fit via Generalized (leave-one-out) cross-validation.""" + n_samples, n_features = X.shape + # 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) + 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) + + # "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) + 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 + + self.coef_ = coef + if cp.isscalar(intercept) or intercept.ndim == 0: + self.intercept_ = float(intercept) + else: + 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_ = best_alpha + self.best_score_ = best_score + 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_ = cv_results.reshape(shape) + + def _fit_cv(self, X, y, sample_weight): + """Fit via an explicit k-fold search over :class:`~cuml.Ridge`.""" + + 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 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) + test = cp.asarray(test) + model = Ridge( + alpha=float(alpha), + fit_intercept=self.fit_intercept, + output_type="cupy", + ) + 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( + alpha=best_alpha, + fit_intercept=self.fit_intercept, + output_type="cupy", + ) + model.fit(X, y, sample_weight=sample_weight) + + self.coef_ = model.coef_ + intercept = model.intercept_ + if cp.isscalar(intercept) or intercept.ndim == 0: + self.intercept_ = float(intercept) + else: + self.intercept_ = intercept + self.alpha_ = best_alpha + self.best_score_ = float(mean_scores[best_idx]) + + @generate_docstring() + @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: + 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, + ) + if self.cv is None: + self._fit_gcv(X, y, sample_weight) + else: + self._fit_cv(X, y, sample_weight) + + return self 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..20a4082532 --- /dev/null +++ b/python/cuml/cuml_accel_tests/integration/test_ridge_cv.py @@ -0,0 +1,79 @@ +# +# 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_best_alpha(regression_data): + X, y = regression_data + model = RidgeCV(alphas=ALPHAS).fit(X, y) + assert model._gpu is not None + + # 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): + 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) diff --git a/python/cuml/tests/test_sklearn_compatibility.py b/python/cuml/tests/test_sklearn_compatibility.py index 06aaa8c613..c2c09e489a 100644 --- a/python/cuml/tests/test_sklearn_compatibility.py +++ b/python/cuml/tests/test_sklearn_compatibility.py @@ -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 # @@ -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(),