sklearn/linear_model/_logistic.py - main - githubmirror/scikit-learn

sklearn/linear_model/_logistic.py
2 685 строк102 KB
Christian Lorentzen
ENH Newton-CD part 7: add NewtonCDSolver (#34561)
07 авг 2026, 19:08
07 авг 2026, 19:08c74a434
"""
Logistic Regression
"""

# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

import inspect
import numbers
import warnings
from numbers import Integral, Real

import numpy as np
from scipy import optimize, sparse

from sklearn._loss.loss import (
    HalfBinomialLoss,
    HalfBinomialLossArrayAPI,
    HalfMultinomialLoss,
    HalfMultinomialLossArrayAPI,
)
from sklearn.base import _fit_context
from sklearn.callback import CallbackSupportMixin
from sklearn.linear_model._base import (
    BaseEstimator,
    LinearClassifierMixin,
    SparseCoefMixin,
)
from sklearn.linear_model._glm._newton_solver import (
    NewtonCDGramSolver,
    NewtonCDSolver,
    NewtonCholeskySolver,
)
from sklearn.linear_model._linear_loss import LinearModelLoss
from sklearn.linear_model._sag import sag_solver
from sklearn.metrics import get_scorer, get_scorer_names, make_scorer
from sklearn.model_selection import check_cv
from sklearn.preprocessing import LabelEncoder
from sklearn.svm._base import _fit_liblinear
from sklearn.utils import (
    check_array,
    check_consistent_length,
    check_random_state,
    compute_class_weight,
    metadata_routing,
)
from sklearn.utils._array_api import (
    _is_numpy_namespace,
    _matching_numpy_dtype,
    _ravel,
    _swapaxes,
    _unravel_index,
    check_same_namespace,
    get_namespace,
    get_namespace_and_device,
    move_to,
    size,
)
from sklearn.utils._indexing import _array_indexing
from sklearn.utils._param_validation import Hidden, Interval, StrOptions
from sklearn.utils.extmath import row_norms, softmax
from sklearn.utils.fixes import _get_additional_lbfgs_options_dict
from sklearn.utils.metadata_routing import (
    MetadataRouter,
    MethodMapping,
    _manual_routing,
    _raise_for_params,
    _routing_enabled,
    process_routing,
)
from sklearn.utils.multiclass import check_classification_targets
from sklearn.utils.optimize import _check_optimize_result, _newton_cg
from sklearn.utils.parallel import Parallel, delayed
from sklearn.utils.validation import (
    _check_method_params,
    _check_sample_weight,
    _deprecate_positional_args,
    check_is_fitted,
    validate_data,
)

_LOGISTIC_SOLVER_CONVERGENCE_MSG = (
    "Please also refer to the documentation for alternative solver options:\n"
    "    https://scikit-learn.org/stable/modules/linear_model.html"
    "#logistic-regression"
)


def _check_solver(solver, penalty, dual):
    if solver not in (
        "liblinear",
        "newton-cd",
        "newton-cd-gram",
        "saga",
    ) and penalty not in ("l2", None):
        raise ValueError(
            f"Solver '{solver}' supports only 'l2' or None penalties, got {penalty} "
            "penalty."
        )
    if solver != "liblinear" and dual:
        raise ValueError(f"Solver {solver} supports only dual=False, got dual={dual}")

    if penalty == "elasticnet" and solver not in (
        "saga",
        "newton-cd",
        "newton-cd-gram",
    ):
        raise ValueError(
            "Only solvers 'newton-cd', 'newton-cd-gram' and 'saga' support elasticnet "
            f"penalty, got solver={solver}."
        )

    if solver == "liblinear" and penalty is None:
        # TODO(1.10): update message to remove "as well as penalty=None".
        raise ValueError(
            "C=np.inf as well as penalty=None is not supported for the liblinear solver"
        )

    return solver


class _LbfgsCallbackBridge:
    """Forward sklearn callbacks to ``scipy.optimize.minimize`` for L-BFGS.

    ``scipy.optimize.minimize`` calls ``callback(xk)`` at the end of each L-BFGS
    iteration. We map that onto a pair of sklearn hooks per iteration:

    - end the iteration that just finished (``on_fit_task_end`` on the current
      subcontext),
    - start the next iteration's subcontext (``on_fit_task_begin``).

    The first iteration's subcontext is created at construction time using the
    initial weights ``w0``. After ``scipy.optimize.minimize`` returns, the last
    subcontext created by ``__call__`` corresponds to an iteration that never ran.
    ``close`` ends it as an empty task so every ``on_fit_task_begin`` call has a
    matching ``on_fit_task_end`` call.
    """

    def __init__(
        self,
        callback_ctx,
        estimator,
        X,
        y,
        metadata,
        w0,
        *,
        n_classes,
        is_binary,
        fit_intercept,
        coefs_order,
        xp,
        device,
    ):
        self._callback_ctx = callback_ctx
        self._estimator = estimator
        self._X = X
        self._y = y
        self._metadata = metadata
        self._coef_kwargs = dict(
            n_classes=n_classes,
            is_binary=is_binary,
            fit_intercept=fit_intercept,
            coefs_order=coefs_order,
            dtype=X.dtype,
            xp=xp,
            device=device,
        )
        # Start the first iteration's subcontext, using the initial weights.
        self._iter_ctx = self._start_iter(w0)

    def _start_iter(self, w):
        coef, intercept = _compute_coef_intercept(w, **self._coef_kwargs)
        ctx = self._callback_ctx.subcontext(task_name="lbfgs-iter")
        ctx.call_on_fit_task_begin(
            estimator=self._estimator,
            X=self._X,
            y=self._y,
            metadata=self._metadata,
            reconstruction_attributes={"coef_": coef, "intercept_": intercept},
        )
        return ctx

    def __call__(self, xk):
        coef, intercept = _compute_coef_intercept(xk, **self._coef_kwargs)
        # End hook for the iteration that scipy just finished.
        self._iter_ctx.call_on_fit_task_end(
            estimator=self._estimator,
            X=self._X,
            y=self._y,
            metadata=self._metadata,
            reconstruction_attributes={"coef_": coef, "intercept_": intercept},
        )
        # TODO(1.10): use the return value of ``call_on_fit_task_end`` (a bool
        # requesting early stopping) to ``raise StopIteration()``. scipy's
        # ``callback`` honours ``StopIteration`` only from scipy > 1.10, which
        # is the min supported scipy in scikit-learn 1.10.

        # Begin hook for the next iteration's subcontext.
        self._iter_ctx = self._start_iter(xk)

    def close(self):
        """Invoke ``on_fit_task_end`` for the left-over subcontext.

        The last subcontext created by ``__call__`` corresponds to an iteration
        that ``scipy.optimize.minimize`` did not actually run. It corresponds to
        an empty task so we don't forward anything to the ``on_fit_task_end``
        hooks.
        """
        self._iter_ctx.call_on_fit_task_end(estimator=self._estimator)


def _compute_coef_intercept(
    w, *, n_classes, is_binary, fit_intercept, coefs_order, dtype, xp, device
):
    """Helper to compute coef and intercept from a flattened array of parameters."""
    if is_binary:
        coef = xp.asarray(w.copy(order=coefs_order), dtype=dtype, device=device)
        if fit_intercept:
            intercept = coef[-1:]
            coef = coef[:-1][None, :]
        else:
            intercept = xp.zeros(1, dtype=dtype, device=device)
            coef = coef[None, :]
    else:
        multi_w = np.reshape(w, (n_classes, -1), order="F")
        coef = xp.asarray(multi_w.copy(order=coefs_order), dtype=dtype, device=device)
        if fit_intercept:
            intercept = coef[:, -1]
            coef = coef[:, :-1]
        else:
            intercept = xp.zeros(n_classes, dtype=dtype, device=device)

    return coef, intercept


def _logistic_regression_path(
    X,
    y,
    *,
    classes,
    Cs=10,
    fit_intercept=True,
    max_iter=100,
    tol=1e-4,
    verbose=0,
    solver="lbfgs",
    coef=None,
    dual=False,
    penalty="l2",
    intercept_scaling=1.0,
    random_state=None,
    check_input=True,
    max_squared_sum=None,
    sample_weight=None,
    l1_ratio=None,
    n_threads=1,
    callback_ctx=None,
    estimator=None,
):
    """Compute a Logistic Regression model for a list of regularization
    parameters.

    This is an implementation that uses the result of the previous model
    to speed up computations along the set of solutions, making it faster
    than sequentially calling LogisticRegression for the different parameters.
    Note that there will be no speedup with liblinear solver, since it does
    not handle warm-starting.

    Read more in the :ref:`User Guide <logistic_regression>`.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Input data.

    y : array-like of shape (n_samples,)
        Input data, label encoded target values, i.e. integers starting from 0.

    classes : ndarray
        A list of class labels known to the classifier.

    Cs : int or array-like of shape (n_cs,), default=10
        List of values for the regularization parameter or integer specifying
        the number of regularization parameters that should be used. In this
        case, the parameters will be chosen in a logarithmic scale between
        1e-4 and 1e4.

    fit_intercept : bool, default=True
        Whether to fit an intercept for the model. In this case the shape of
        the returned array is (n_cs, n_features + 1).

    max_iter : int, default=100
        Maximum number of iterations for the solver.

    tol : float, default=1e-4
        Stopping criterion. For the newton-cg and lbfgs solvers, the iteration
        will stop when ``max{|g_i | i = 1, ..., n} <= tol``
        where ``g_i`` is the i-th component of the gradient.

    verbose : int, default=0
        For the liblinear and lbfgs solvers set verbose to any positive
        number for verbosity.

    solver : {'lbfgs', 'liblinear', 'newton-cd', 'newton-cd-gram', 'newton-cg', \
            'newton-cholesky', 'sag', 'saga'}, default='lbfgs'
        Numerical solver to use.

    coef : array-like of shape (n_classes, features + int(fit_intercept)) or \
            (1, n_features + int(fit_intercept)) or \
            (n_features + int(fit_intercept)), default=None
        Initialization value for coefficients of logistic regression.
        Useless for liblinear solver.

    dual : bool, default=False
        Dual or primal formulation. Dual formulation is only implemented for
        l2 penalty with liblinear solver. Prefer dual=False when
        n_samples > n_features.

    penalty : {'l1', 'l2', 'elasticnet'}, default='l2'
        Used to specify the norm used in the penalization. The 'newton-cg',
        'sag' and 'lbfgs' solvers support only l2 penalties. 'elasticnet' is
        only supported by the 'saga' solver.

    intercept_scaling : float, default=1.
        Useful only when the solver `liblinear` is used
        and `self.fit_intercept` is set to `True`. In this case, `x` becomes
        `[x, self.intercept_scaling]`,
        i.e. a "synthetic" feature with constant value equal to
        `intercept_scaling` is appended to the instance vector.
        The intercept becomes
        ``intercept_scaling * synthetic_feature_weight``.

        .. note::
            The synthetic feature weight is subject to L1 or L2
            regularization as all other features.
            To lessen the effect of regularization on synthetic feature weight
            (and therefore on the intercept) `intercept_scaling` has to be increased.

    random_state : int, RandomState instance, default=None
        Used when ``solver`` == 'sag', 'saga' or 'liblinear' to shuffle the
        data. See :term:`Glossary <random_state>` for details.

    check_input : bool, default=True
        If False, the input arrays X and y (and sample_weight) will not be checked.

    max_squared_sum : float, default=None
        Maximum squared sum of X over samples. Used only in SAG solver.
        If None, it will be computed, going through all the samples.
        The value should be precomputed to speed up cross validation.

    sample_weight : array-like of shape (n_samples,), default=None
        Array of weights that are assigned to individual samples.
        If not provided, then each sample is given unit weight.

    l1_ratio : float, default=None
        The Elastic-Net mixing parameter, with ``0 <= l1_ratio <= 1``. Only
        used if ``penalty='elasticnet'``. Setting ``l1_ratio=0`` is equivalent
        to using ``penalty='l2'``, while setting ``l1_ratio=1`` is equivalent
        to using ``penalty='l1'``. For ``0 < l1_ratio <1``, the penalty is a
        combination of L1 and L2.

    n_threads : int, default=1
       Number of OpenMP threads to use.

    callback_ctx : CallbackContext or None, default=None
        The callback context of the fit task calling this function. If set to None, the
        callbacks will not be invoked.

    estimator : estimator instance or None, default=None
        The estimator instance in which fit this function is called, forwarded to the
        callbacks if callback_ctx is not None.

    Returns
    -------
    coefs : ndarray of shape (n_cs, n_classes, n_features + int(fit_intercept)) or \
            (n_cs, n_features + int(fit_intercept))
        List of coefficients for the Logistic Regression model. If fit_intercept is set
        to True, then the last dimension will be n_features + 1, where the last item
        represents the intercept.
        For binary problems the second dimension in n_classes is dropped, i.e. the shape
        will be `(n_cs, n_features + int(fit_intercept))`.

    Cs : ndarray
        Grid of Cs used for cross-validation.

    n_iter : array of shape (n_cs,)
        Actual number of iteration for each C in Cs.

    Notes
    -----
    You might get slightly different results with the solver liblinear than
    with the others since this uses LIBLINEAR which penalizes the intercept.

    .. versionchanged:: 0.19
        The "copy" parameter was removed.
    """
    if isinstance(Cs, numbers.Integral):
        Cs = np.logspace(-4, 4, Cs)

    solver = _check_solver(solver, penalty, dual)
    xp, _, device = get_namespace_and_device(X)
    # Only newton-cg has complete support of the array API, lbfgs still needs
    # coef / w0 as numpy arrays.
    coef_as_xp = solver == "newton-cg"

    # Preprocessing.
    if check_input:
        X = check_array(
            X,
            accept_sparse="csc" if solver == "newton-cd" else "csr",
            dtype=[xp.float64, xp.float32],
            order="F" if solver == "newton-cd" else None,
            accept_large_sparse=solver not in ("liblinear", "newton-cd", "sag", "saga"),
        )
        y = check_array(y, ensure_2d=False, dtype=None)
        check_consistent_length(X, y)

        if sample_weight is not None:
            sample_weight = _check_sample_weight(
                sample_weight, X, dtype=X.dtype, copy=True, ensure_same_device=True
            )

    n_samples, n_features = X.shape
    n_classes = size(classes)
    is_binary = n_classes == 2

    if solver == "liblinear" and not is_binary:
        raise ValueError(
            "The 'liblinear' solver does not support multiclass classification"
            " (n_classes >= 3). Either use another solver or wrap the "
            "estimator in a OneVsRestClassifier to keep applying a "
            "one-versus-rest scheme."
        )

    # TODO(newton-cd): add multiclass support
    if solver == "newton-cd" and not is_binary:
        raise ValueError(
            f"The '{solver}' solver does not support multiclass classification"
            f" (n_classes >= 3); got {n_classes=}."
        )

    random_state = check_random_state(random_state)

    if is_binary:
        # y is already encoded as values in {0, 1}.
        if coef_as_xp:
            w0 = xp.zeros(
                n_features + int(fit_intercept),
                dtype=X.dtype,
                device=device,
            )
        else:
            w0 = np.zeros(
                n_features + int(fit_intercept),
                dtype=_matching_numpy_dtype(X, xp=xp),
            )
        if solver == "liblinear":
            y = 2 * y - 1  # liblinear requires target values -1, +1
    else:
        if coef_as_xp:
            w0 = xp.zeros(
                (size(classes), n_features + int(fit_intercept)),
                dtype=X.dtype,
                device=device,
            )
        else:
            # It is important that w0 is F-contiguous.
            w0 = np.zeros(
                (size(classes), n_features + int(fit_intercept)),
                order="F",
                dtype=_matching_numpy_dtype(X, xp=xp),
            )

    # IMPORTANT NOTE:
    # All solvers relying on LinearModelLoss need to scale the penalty with n_samples
    # or the sum of sample weights because the implemented logistic regression
    # objective here is (unfortunately)
    #     C * sum(pointwise_loss) + penalty
    # instead of (as LinearModelLoss does)
    #     mean(pointwise_loss) + 1/C * penalty
    if solver in (
        "lbfgs",
        "newton-cd",
        "newton-cd-gram",
        "newton-cg",
        "newton-cholesky",
    ):
        # This needs to be calculated after sample_weight is multiplied by
        # class_weight. It is even tested that passing class_weight is equivalent to
        # passing sample_weights according to class_weight.
        sw_sum = n_samples if sample_weight is None else float(xp.sum(sample_weight))

    if coef is not None:
        if is_binary:
            if coef.ndim == 1 and coef.shape[0] == n_features + int(fit_intercept):
                w0[:] = move_to(coef, xp=xp, device=device) if coef_as_xp else coef
            elif (
                coef.ndim == 2
                and coef.shape[0] == 1
                and coef.shape[1] == n_features + int(fit_intercept)
            ):
                w0[:] = (
                    move_to(coef[0], xp=xp, device=device) if coef_as_xp else coef[0]
                )
            else:
                msg = (
                    f"Initialization coef is of shape {coef.shape}, expected shape "
                    f"{w0.shape} or (1, {w0.shape[0]})"
                )
                raise ValueError(msg)
        else:
            if (
                coef.ndim == 2
                and coef.shape[0] == n_classes
                and coef.shape[1] == n_features + int(fit_intercept)
            ):
                w0[:, : coef.shape[1]] = (
                    move_to(coef, xp=xp, device=device) if coef_as_xp else coef
                )
            else:
                msg = (
                    f"Initialization coef is of shape {coef.shape}, expected shape "
                    f"{w0.shape}"
                )
                raise ValueError(msg)

    if is_binary:
        loss = LinearModelLoss(
            base_loss=(
                HalfBinomialLoss()
                if _is_numpy_namespace(xp)
                else HalfBinomialLossArrayAPI(xp=xp, device=device)
            ),
            fit_intercept=fit_intercept,
        )
        if solver == "lbfgs":
            func = loss.loss_gradient
        elif solver == "newton-cg":
            func = loss.loss
            grad = loss.gradient
            hess = loss.gradient_hessian_product  # hess = [gradient, hessp]
        warm_start_sag = {"coef": w0[:, None]}
    else:  # multinomial
        loss = LinearModelLoss(
            base_loss=(
                HalfMultinomialLoss(n_classes=size(classes))
                if _is_numpy_namespace(xp)
                else HalfMultinomialLossArrayAPI(
                    n_classes=size(classes), xp=xp, device=device
                )
            ),
            fit_intercept=fit_intercept,
        )
        if solver in (
            "lbfgs",
            "newton-cd",
            "newton-cd-gram",
            "newton-cg",
            "newton-cholesky",
        ):
            # scipy.optimize.minimize and newton-cg accept only ravelled parameters,
            # i.e. 1d-arrays. LinearModelLoss expects classes to be contiguous and
            # reconstructs the 2d-array via w0.reshape((n_classes, -1), order="F").
            # As w0 is F-contiguous, ravel(order="F") also avoids a copy.
            if _is_numpy_namespace(xp) or not coef_as_xp:
                w0 = w0.ravel(order="F")
            else:
                w0 = _ravel(w0.T, xp=xp)
        if solver == "lbfgs":
            func = loss.loss_gradient
        elif solver == "newton-cg":
            func = loss.loss
            grad = loss.gradient
            hess = loss.gradient_hessian_product  # hess = [gradient, hessp]
        warm_start_sag = {"coef": w0.T} if w0.ndim > 1 else {"coef": w0}

    coefs = list()
    n_iter = xp.zeros(len(Cs), dtype=xp.int32, device=device)
    coefs_order = "C" if not _is_numpy_namespace(xp) else "K"
    callback_metadata = (
        {"sample_weight": sample_weight} if sample_weight is not None else None
    )
    for i, C in enumerate(Cs):
        if solver == "lbfgs":
            # In LogisticRegression.fit, Cs is always a one-element list, so we don't
            # add additional callback subcontexts inside this for-loop to avoid
            # introducing an unnecessary subtask level.
            # TODO(callbacks) When adding callback support to LogisticRegressionCV,
            # another subcontext level will be necessary, so we'll need to add a level
            # only in the LogisticRegressionCV, keeping the LogisticRegression version
            # without that extra level.
            lbfgs_bridge = None
            if callback_ctx is not None and hasattr(estimator, "_skl_callbacks"):
                lbfgs_bridge = _LbfgsCallbackBridge(
                    callback_ctx,
                    estimator,
                    X,
                    y,
                    callback_metadata,
                    w0,
                    n_classes=n_classes,
                    is_binary=is_binary,
                    fit_intercept=fit_intercept,
                    coefs_order=coefs_order,
                    xp=xp,
                    device=device,
                )
            l2_reg_strength = 1.0 / (C * sw_sum)
            iprint = [-1, 50, 1, 100, 101][
                np.searchsorted(np.array([0, 1, 2, 3]), verbose)
            ]
            opt_res = optimize.minimize(
                func,
                w0,
                method="L-BFGS-B",
                jac=True,
                args=(X, y, sample_weight, 0, l2_reg_strength, n_threads),
                options={
                    "maxiter": max_iter,
                    "maxls": 50,  # default is 20
                    "gtol": tol,
                    "ftol": 64 * np.finfo(float).eps,
                    **_get_additional_lbfgs_options_dict("iprint", iprint),
                },
                callback=lbfgs_bridge,
            )
            n_iter_i = _check_optimize_result(
                solver,
                opt_res,
                max_iter,
                extra_warning_msg=_LOGISTIC_SOLVER_CONVERGENCE_MSG,
            )
            w0, loss = opt_res.x, opt_res.fun
            if lbfgs_bridge is not None:
                lbfgs_bridge.close()
            if fit_intercept and not is_binary:
                # Use freedom to add the same constant to all classes in order to
                # achieve the symmetric parametrization with sum(intercept) = 0
                w0[-n_classes:] -= np.mean(w0[-n_classes:])
        elif solver == "newton-cg":
            l2_reg_strength = 1.0 / (C * sw_sum)
            args = (X, y, sample_weight, 0, l2_reg_strength, n_threads)
            w0, n_iter_i = _newton_cg(
                grad_hess=hess,
                func=func,
                grad=grad,
                x0=w0,
                args=args,
                maxiter=max_iter,
                tol=tol,
                verbose=verbose,
            )
            if fit_intercept and not is_binary:
                # Use freedom to add the same constant to all classes in order to
                # achieve the symmetric parametrization with sum(intercept) = 0
                w0[-n_classes:] -= xp.mean(w0[-n_classes:])
        elif solver == "newton-cholesky":
            l2_reg_strength = 1.0 / (C * sw_sum)
            sol = NewtonCholeskySolver(
                coef=w0,
                linear_loss=loss,
                l2_reg_strength=l2_reg_strength,
                tol=tol,
                max_iter=max_iter,
                n_threads=n_threads,
                verbose=verbose,
            )
            w0 = sol.solve(X=X, y=y, sample_weight=sample_weight)
            n_iter_i = sol.iteration
        elif solver in ("newton-cd", "newton-cd-gram"):
            if penalty == "l1":
                l1_ratio = 1.0
            elif penalty == "l2":
                l1_ratio = 0
            l1_reg_strength = l1_ratio / (C * sw_sum)
            l2_reg_strength = (1.0 - l1_ratio) / (C * sw_sum)
            sol = {"newton-cd": NewtonCDSolver, "newton-cd-gram": NewtonCDGramSolver}
            sol = sol[solver](
                coef=w0,
                linear_loss=loss,
                l1_reg_strength=l1_reg_strength,
                l2_reg_strength=l2_reg_strength,
                tol=tol,
                max_iter=max_iter,
                n_threads=n_threads,
                verbose=verbose,
            )
            w0 = sol.solve(X=X, y=y, sample_weight=sample_weight)
            n_iter_i = sol.iteration
        elif solver == "liblinear":
            coef_, intercept_, n_iter_i = _fit_liblinear(
                X,
                y,
                C,
                fit_intercept,
                intercept_scaling,
                None,
                penalty,
                dual,
                verbose,
                max_iter,
                tol,
                random_state,
                sample_weight=sample_weight,
            )
            if fit_intercept:
                w0 = np.concatenate([coef_.ravel(), intercept_])
            else:
                w0 = coef_.ravel()
            # n_iter_i is an array for each class. However, y is encoded
            # in {-1, 1}, so we only take the first element of n_iter_i.
            n_iter_i = n_iter_i.item()

        elif solver in ["sag", "saga"]:
            if is_binary:
                loss = "log"
            else:
                loss = "multinomial"
            # alpha is for L2-norm, beta is for L1-norm
            if penalty == "l1":
                alpha = 0.0
                beta = 1.0 / C
            elif penalty == "l2":
                alpha = 1.0 / C
                beta = 0.0
            else:  # Elastic-Net penalty
                alpha = (1.0 / C) * (1 - l1_ratio)
                beta = (1.0 / C) * l1_ratio

            w0, n_iter_i, warm_start_sag = sag_solver(
                X,
                y,
                sample_weight,
                loss,
                alpha,
                beta,
                max_iter,
                tol,
                verbose,
                random_state,
                False,
                max_squared_sum,
                warm_start_sag,
                is_saga=(solver == "saga"),
            )

        else:
            msg = (
                "solver must be one of {'lbfgs', 'liblinear', 'newton-cg', "
                "'newton-cholesky', 'sag', 'saga'}, "
                f"got '{solver}' instead."
            )
            raise ValueError(msg)

        if is_binary:
            if _is_numpy_namespace(xp):
                coefs.append(np.asarray(w0.copy(order=coefs_order), dtype=X.dtype))
            else:
                coefs.append(xp.asarray(w0, copy=True, dtype=X.dtype, device=device))
        else:
            if solver in (
                "lbfgs",
                "newton-cd",
                "newton-cd-gram",
                "newton-cg",
                "newton-cholesky",
            ):
                if _is_numpy_namespace(xp) or not coef_as_xp:
                    multi_w0 = np.reshape(w0, (n_classes, -1), order="F")
                else:
                    multi_w0 = xp.reshape(w0, (-1, n_classes)).T
            else:
                multi_w0 = w0
            coefs.append(xp.asarray(multi_w0, copy=True, dtype=X.dtype, device=device))

        n_iter[i] = n_iter_i

    return xp.stack(coefs), xp.asarray(Cs, device=device, dtype=X.dtype), n_iter


# helper function for LogisticCV
def _log_reg_scoring_path(
    X,
    y,
    train,
    test,
    *,
    classes,
    Cs,
    scoring,
    fit_intercept,
    max_iter,
    tol,
    verbose,
    solver,
    penalty,
    dual,
    intercept_scaling,
    random_state,
    max_squared_sum,
    sample_weight,
    l1_ratio,
    score_params,
):
    """Computes scores across logistic_regression_path

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Training data.

    y : array-like of shape (n_samples,)
        Input data, label encoded target values, i.e. integers starting from 0.

    train : list of indices
        The indices of the train set.

    test : list of indices
        The indices of the test set.

    classes : ndarray
        A list of class labels known to the classifier.

    Cs : int or list of floats
        Each of the values in Cs describes the inverse of
        regularization strength. If Cs is as an int, then a grid of Cs
        values are chosen in a logarithmic scale between 1e-4 and 1e4.

    scoring : str, callable or None
        The scoring method to use for cross-validation. Options:

        - str: see :ref:`scoring_string_names` for options.
        - callable: a scorer callable object (e.g., function) with signature
          ``scorer(estimator, X, y)``. See :ref:`scoring_callable` for details.
        - `None`: :ref:`accuracy <accuracy_score>` is used.

    fit_intercept : bool
        If False, then the bias term is set to zero. Else the last
        term of each coef_ gives us the intercept.

    max_iter : int
        Maximum number of iterations for the solver.

    tol : float
        Tolerance for stopping criteria.

    verbose : int
        For the liblinear and lbfgs solvers set verbose to any positive
        number for verbosity.

    solver : {'lbfgs', 'liblinear', 'newton-cg', 'newton-cholesky', 'sag', 'saga'}
        Decides which solver to use.

    penalty : {'l1', 'l2', 'elasticnet'}
        Used to specify the norm used in the penalization. The 'newton-cg',
        'sag' and 'lbfgs' solvers support only l2 penalties. 'elasticnet' is
        only supported by the 'saga' solver.

    dual : bool
        Dual or primal formulation. Dual formulation is only implemented for
        l2 penalty with liblinear solver. Prefer dual=False when
        n_samples > n_features.

    intercept_scaling : float
        Useful only when the solver `liblinear` is used
        and `self.fit_intercept` is set to `True`. In this case, `x` becomes
        `[x, self.intercept_scaling]`,
        i.e. a "synthetic" feature with constant value equal to
        `intercept_scaling` is appended to the instance vector.
        The intercept becomes
        ``intercept_scaling * synthetic_feature_weight``.

        .. note::
            The synthetic feature weight is subject to L1 or L2
            regularization as all other features.
            To lessen the effect of regularization on synthetic feature weight
            (and therefore on the intercept) `intercept_scaling` has to be increased.

    random_state : int, RandomState instance
        Used when ``solver`` == 'sag', 'saga' or 'liblinear' to shuffle the
        data. See :term:`Glossary <random_state>` for details.

    max_squared_sum : float
        Maximum squared sum of X over samples. Used only in SAG solver.
        If None, it will be computed, going through all the samples.
        The value should be precomputed to speed up cross validation.

    sample_weight : array-like of shape (n_samples,)
        Array of weights that are assigned to individual samples.
        If not provided, then each sample is given unit weight.

    l1_ratio : float
        The Elastic-Net mixing parameter, with ``0 <= l1_ratio <= 1``. Only
        used if ``penalty='elasticnet'``. Setting ``l1_ratio=0`` is equivalent
        to using ``penalty='l2'``, while setting ``l1_ratio=1`` is equivalent
        to using ``penalty='l1'``. For ``0 < l1_ratio <1``, the penalty is a
        combination of L1 and L2.

    score_params : dict
        Parameters to pass to the `score` method of the underlying scorer.

    Returns
    -------
    coefs : ndarray of shape (n_cs, n_classes, n_features + int(fit_intercept)) or \
            (n_cs, n_features + int(fit_intercept))
        List of coefficients for the Logistic Regression model. If fit_intercept is set
        to True, then the last dimension will be n_features + 1, where the last item
        represents the intercept.
        For binary problems the second dimension in n_classes is dropped, i.e. the shape
        will be `(n_cs, n_features + int(fit_intercept))`.

    Cs : ndarray of shape (n_cs,)
        Grid of Cs used for cross-validation.

    scores : ndarray of shape (n_cs,)
        Scores obtained for each Cs.

    n_iter : ndarray of shape (n_cs,)
        Actual number of iteration for each C in Cs.
    """
    xp, _, device = get_namespace_and_device(X)
    train_xp = xp.asarray(train, device=device)
    test_xp = xp.asarray(test, device=device)
    indices_dtype = train_xp.dtype
    X_train = _array_indexing(X, train_xp, indices_dtype, axis=0)
    X_test = _array_indexing(X, test_xp, indices_dtype, axis=0)
    y_train = _array_indexing(y, train_xp, indices_dtype, axis=0)
    y_test = _array_indexing(y, test_xp, indices_dtype, axis=0)

    sw_train, sw_test = None, None
    if sample_weight is not None:
        sample_weight = _check_sample_weight(sample_weight, X)
        sw_train = sample_weight[train_xp]
        sw_test = sample_weight[test_xp]

    if solver == "newton-cd":
        if sparse.issparse(X_train):
            X_train = X_train.tocsc()
        else:
            X_train = np.asfortranarray(X_train)

    # Note: We pass classes for the whole dataset to avoid inconsistencies,
    # i.e. different number of classes in different folds. This way, if a class
    # is not present in a fold, _logistic_regression_path will still return
    # coefficients associated to this class.
    coefs, Cs, n_iter = _logistic_regression_path(
        X_train,
        y_train,
        classes=classes,
        Cs=Cs,
        l1_ratio=l1_ratio,
        fit_intercept=fit_intercept,
        solver=solver,
        max_iter=max_iter,
        tol=tol,
        verbose=verbose,
        dual=dual,
        penalty=penalty,
        intercept_scaling=intercept_scaling,
        random_state=random_state,
        check_input=False,
        max_squared_sum=max_squared_sum,
        sample_weight=sw_train,
        callback_ctx=None,
        estimator=None,
    )

    log_reg = LogisticRegression(solver=solver)

    n_classes = size(classes)
    # The score method of Logistic Regression has a classes_ attribute.
    # As y is already expected to be label encoded, we don't set
    #     log_reg.classes_ = classes
    # but instead
    log_reg.classes_ = xp.arange(n_classes, dtype=X.dtype, device=device)

    scores = list()

    # Prepare the call to get the score per fold: calc_score
    scoring = get_scorer(scoring)
    if scoring is None:

        def calc_score(log_reg):
            return log_reg.score(X_test, y_test, sample_weight=sw_test)

    else:
        is_binary = n_classes <= 2
        score_params = score_params or {}
        score_params = _check_method_params(X=X, params=score_params, indices=test)
        # We need to pass the classes as "labels" argument to scorers that support
        # it, e.g. scoring = "neg_brier_score", because y_test may not contain all
        # class labels.
        # There are at least 2 possibilities:
        # 1. Metadata routing is enabled: A try except clause is possible with
        #   adding labels to score_params. We could then pass the already instantiated
        #   log_reg instance to scoring.
        # 2. We reconstruct the scorer and pass labels as kwargs explicitly.
        # We implement the 2nd option even if it seems a bit hacky because it works
        # with and without metadata routing.
        if hasattr(scoring, "_score_func"):
            sig = inspect.signature(scoring._score_func).parameters
        else:
            sig = []

        if "labels" in sig:
            pos_label_kwarg = {}
            if is_binary and "pos_label" in sig:
                # see _logistic_regression_path
                pos_label_kwarg["pos_label"] = n_classes - 1  # classes[-1]
            scoring = make_scorer(
                scoring._score_func,
                greater_is_better=True if scoring._sign == 1 else False,
                response_method=scoring._response_method,
                labels=log_reg.classes_,
                **pos_label_kwarg,
                **getattr(scoring, "_kwargs", {}),
            )

        def calc_score(log_reg):
            return scoring(log_reg, X_test, y_test, **score_params)

    for i in range(size(Cs)):
        w = coefs[i, ...]
        C = Cs[i]
        log_reg.C = C
        if fit_intercept:
            log_reg.coef_ = w[..., :-1]
            log_reg.intercept_ = w[..., -1]
        else:
            log_reg.coef_ = w
            log_reg.intercept_ = 0.0

        scores.append(calc_score(log_reg))

    scores = xp.asarray(scores, device=device, dtype=X.dtype)
    return coefs, Cs, scores, n_iter


class LogisticRegression(
    CallbackSupportMixin, LinearClassifierMixin, SparseCoefMixin, BaseEstimator
):
    """
    Logistic Regression (aka logit, MaxEnt) classifier.

    This class implements regularized logistic regression using a set of available
    solvers. **Note that regularization is applied by default**. It can handle both
    dense and sparse input `X`. Use C-ordered arrays or CSR matrices containing 64-bit
    floats for optimal performance; any other input format will be converted (and
    copied).

    The solvers 'lbfgs', 'newton-cg', 'newton-cholesky' and 'sag' support only L2
    regularization with primal formulation, or no regularization. The 'liblinear'
    solver supports both L1 and L2 regularization (but not both, i.e. elastic-net),
    with a dual formulation only for the L2 penalty. The Elastic-Net (combination of L1
    and L2) regularization is only supported by the 'saga' solver.

    For :term:`multiclass` problems (whenever `n_classes >= 3`), all solvers except
    'liblinear' optimize the (penalized) multinomial loss. 'liblinear' only handles
    binary classification but can be extended to handle multiclass by using
    :class:`~sklearn.multiclass.OneVsRestClassifier`.

    Read more in the :ref:`User Guide <logistic_regression>`.

    Parameters
    ----------
    penalty : {'l1', 'l2', 'elasticnet', None}, default='l2'
        Specify the norm of the penalty:

        - `None`: no penalty is added;
        - `'l2'`: add an L2 penalty term and it is the default choice;
        - `'l1'`: add an L1 penalty term;
        - `'elasticnet'`: both L1 and L2 penalty terms are added.

        .. warning::
           Some penalties may not work with some solvers. See the parameter
           `solver` below, to know the compatibility between the penalty and
           solver.

        .. versionadded:: 0.19
           l1 penalty with SAGA solver (allowing 'multinomial' + L1)

        .. deprecated:: 1.8
           `penalty` was deprecated in version 1.8 and will be removed in 1.10.
           Use `l1_ratio` and `C` instead. `l1_ratio=0` for `penalty='l2'`,
           `l1_ratio=1` for `penalty='l1'`, `l1_ratio` set to any float between 0 and 1
           for `penalty='elasticnet'`, and `C=np.inf` for `penalty=None`.

    C : float, default=1.0
        Inverse of regularization strength; must be a positive float.
        Like in support vector machines, smaller values specify stronger
        regularization. `C=np.inf` results in unpenalized logistic regression.
        For a visual example on the effect of tuning the `C` parameter
        with an L1 penalty, see:
        :ref:`sphx_glr_auto_examples_linear_model_plot_logistic_path.py`.

    l1_ratio : float, default=0.0
        The Elastic-Net mixing parameter, with `0 <= l1_ratio <= 1`. Setting
        `l1_ratio=1` gives a pure L1-penalty, setting `l1_ratio=0` gives a pure
        L2-penalty. Any value between 0 and 1 gives an Elastic-Net penalty of the form
        `l1_ratio * L1 + (1 - l1_ratio) * L2`.

        .. warning::
           Certain values of `l1_ratio`, i.e. some penalties, may not work with some
           solvers. See the parameter `solver` below, to know the compatibility between
           the penalty and solver.

        .. versionchanged:: 1.8
            Default value changed from None to 0.0.

        .. deprecated:: 1.8
            `None` is deprecated and will be removed in version 1.10. Always use
            `l1_ratio` to specify the penalty type.

    dual : bool, default=False
        Dual (constrained) or primal (regularized, see also
        :ref:`this equation <regularized-logistic-loss>`) formulation. Dual formulation
        is only implemented for l2 penalty with liblinear solver. Prefer `dual=False`
        when n_samples > n_features.

    tol : float, default=1e-4
        Tolerance for stopping criteria.

    fit_intercept : bool, default=True
        Specifies if a constant (a.k.a. bias or intercept) should be
        added to the decision function.

    intercept_scaling : float, default=1
        Useful only when the solver `liblinear` is used
        and `self.fit_intercept` is set to `True`. In this case, `x` becomes
        `[x, self.intercept_scaling]`,
        i.e. a "synthetic" feature with constant value equal to
        `intercept_scaling` is appended to the instance vector.
        The intercept becomes
        ``intercept_scaling * synthetic_feature_weight``.

        .. note::
            The synthetic feature weight is subject to L1 or L2
            regularization as all other features.
            To lessen the effect of regularization on synthetic feature weight
            (and therefore on the intercept) `intercept_scaling` has to be increased.

    class_weight : dict or 'balanced', default=None
        Weights associated with classes in the form ``{class_label: weight}``.
        If not given, all classes are supposed to have weight one.

        The "balanced" mode uses the values of y to automatically adjust
        weights inversely proportional to class frequencies in the input data
        as ``n_samples / (n_classes * np.bincount(y))``.

        Note that these weights will be multiplied with sample_weight (passed
        through the fit method) if sample_weight is specified.

        .. versionadded:: 0.17
           *class_weight='balanced'*

    random_state : int, RandomState instance, default=None
        Only used for `solver` == 'sag', 'saga' or 'liblinear' to shuffle the
        data. It has no effect on the other solvers.
        See :term:`Glossary <random_state>` for details.

    solver : {'lbfgs', 'liblinear', 'newton-cd', 'newton-cd-gram', 'newton-cg', \
            'newton-cholesky', 'sag', 'saga'}, default='lbfgs'

        Algorithm to use in the optimization problem. Default is 'lbfgs'.
        To choose a solver, you might want to consider the following aspects:

        - 'lbfgs' is a good default solver because it works reasonably well for a wide
          class of problems.
        - For :term:`multiclass` problems (`n_classes >= 3`), all solvers except
          'liblinear' minimize the full multinomial loss, 'liblinear' will raise an
          error.
        - 'newton-cholesky' is a good choice for
          `n_samples` >> `n_features * n_classes`, especially with one-hot encoded
          categorical features with rare categories. Be aware that the memory usage
          of this solver has a quadratic dependency on `n_features * n_classes`
          because it explicitly computes the full Hessian matrix.
        - For small datasets, 'liblinear' is a good choice, whereas 'sag'
          and 'saga' are faster for large ones;
        - 'liblinear' can only handle binary classification by default. To apply a
          one-versus-rest scheme for the multiclass setting one can wrap it with the
          :class:`~sklearn.multiclass.OneVsRestClassifier`.

        .. warning::
           The choice of the algorithm depends on the penalty chosen (`l1_ratio=0`
           for L2-penalty, `l1_ratio=1` for L1-penalty and `0 < l1_ratio < 1` for
           Elastic-Net) and on (multinomial) multiclass support:

           ================= ======================== ======================
           solver            l1_ratio                 multinomial multiclass
           ================= ======================== ======================
           'lbfgs'           l1_ratio=0               yes
           'liblinear'       l1_ratio=1 or l1_ratio=0 no
           'newton-cd'       0<=l1_ratio<=1           no
           'newton-cd-gram'  0<=l1_ratio<=1           yes
           'newton-cg'       l1_ratio=0               yes
           'newton-cholesky' l1_ratio=0               yes
           'sag'             l1_ratio=0               yes
           'saga'            0<=l1_ratio<=1           yes
           ================= ======================== ======================

        .. note::
           'sag' and 'saga' fast convergence is only guaranteed on features
           with approximately the same scale. You can preprocess the data with
           a scaler from :mod:`sklearn.preprocessing`.

        .. seealso::
           Refer to the :ref:`User Guide <Logistic_regression>` for more
           information regarding :class:`LogisticRegression` and more specifically the
           :ref:`Table <logistic_regression_solvers>`
           summarizing solver/penalty supports.

        .. versionadded:: 0.17
           Stochastic Average Gradient (SAG) descent solver. Multinomial support in
           version 0.18.
        .. versionadded:: 0.19
           SAGA solver.
        .. versionchanged:: 0.22
           The default solver changed from 'liblinear' to 'lbfgs' in 0.22.
        .. versionadded:: 1.2
           newton-cholesky solver. Multinomial support in version 1.6.

    max_iter : int, default=100
        Maximum number of iterations taken for the solvers to converge.

    verbose : int, default=0
        For the liblinear and lbfgs solvers set verbose to any positive
        number for verbosity.

    warm_start : bool, default=False
        When set to True, reuse the solution of the previous call to fit as
        initialization, otherwise, just erase the previous solution.
        Useless for liblinear solver. See :term:`the Glossary <warm_start>`.

        .. versionadded:: 0.17
           *warm_start* to support *lbfgs*, *newton-cg*, *sag*, *saga* solvers.

    n_jobs : int, default=None
        Does not have any effect.

        .. deprecated:: 1.8
           `n_jobs` is deprecated in version 1.8 and will be removed in 1.10.

    Attributes
    ----------

    classes_ : ndarray of shape (n_classes, )
        A list of class labels known to the classifier.

    coef_ : ndarray or CSR matrix of shape (1, n_features) or (n_classes, n_features)
        Coefficients of the features in the decision function.

        `coef_` is of shape (1, n_features) when the given problem is binary.

        By default, it will be created as a dense array, but can be turned to
        sparse (CSR format) through :meth:`sparsify` (which can be beneficial
        under L1 regularization when many coefficients are zero), and back to
        dense through :meth:`densify`.

    intercept_ : ndarray of shape (1,) or (n_classes,)
        Intercept (a.k.a. bias) added to the decision function.

        If `fit_intercept` is set to False, the intercept is set to zero.
        `intercept_` is of shape (1,) when the given problem is binary.

    n_features_in_ : int
        Number of features seen during :term:`fit`.

        .. versionadded:: 0.24

    feature_names_in_ : ndarray of shape (`n_features_in_`,)
        Names of features seen during :term:`fit`. Defined only when `X`
        has feature names that are all strings.

        .. versionadded:: 1.0

    n_iter_ : ndarray of shape (1, )
        Actual number of iterations for all classes.

        .. versionchanged:: 0.20

            In SciPy <= 1.0.0 the number of lbfgs iterations may exceed
            ``max_iter``. ``n_iter_`` will now report at most ``max_iter``.

    See Also
    --------
    SGDClassifier : Incrementally trained logistic regression (when given
        the parameter ``loss="log_loss"``).
    LogisticRegressionCV : Logistic regression with built-in cross validation.

    Notes
    -----
    For several reasons (floating point arithmetic, random number generators, etc.)
    the coefficients of a fitted model might differ slightly (among machines,
    scikit-learn versions, etc.) for the same input data. If that happens and you
    want to avoid it, you can try with a smaller `tol` parameter.

    Predict output may not match that of standalone liblinear in certain
    cases. See :ref:`differences from liblinear <liblinear_differences>`
    in the narrative documentation.

    Examples
    --------
    >>> from sklearn.datasets import load_iris
    >>> from sklearn.linear_model import LogisticRegression
    >>> X, y = load_iris(return_X_y=True)
    >>> clf = LogisticRegression().fit(X, y)
    >>> clf.predict(X[:2, :])
    array([0, 0])
    >>> clf.predict_proba(X[:2, :])
    array([[9.82e-01, 1.82e-02, 1.44e-08],
           [9.72e-01, 2.82e-02, 3.02e-08]])
    >>> clf.score(X, y)
    0.97

    For a comparison of the LogisticRegression with other classifiers see:
    :ref:`sphx_glr_auto_examples_classification_plot_classification_probability.py`.
    """

    _parameter_constraints: dict = {
        "penalty": [
            StrOptions({"l1", "l2", "elasticnet"}),
            None,
            Hidden(StrOptions({"deprecated"})),
        ],
        "C": [Interval(Real, 0, None, closed="right")],
        "l1_ratio": [Interval(Real, 0, 1, closed="both"), None],
        "dual": ["boolean"],
        "tol": [Interval(Real, 0, None, closed="left")],
        "fit_intercept": ["boolean"],
        "intercept_scaling": [Interval(Real, 0, None, closed="neither")],
        "class_weight": [dict, StrOptions({"balanced"}), None],
        "random_state": ["random_state"],
        "solver": [
            StrOptions(
                {
                    "lbfgs",
                    "liblinear",
                    "newton-cd",
                    "newton-cd-gram",
                    "newton-cg",
                    "newton-cholesky",
                    "sag",
                    "saga",
                }
            )
        ],
        "max_iter": [Interval(Integral, 0, None, closed="left")],
        "verbose": ["verbose"],
        "warm_start": ["boolean"],
        "n_jobs": [None, Integral],
    }

    def __init__(
        self,
        penalty="deprecated",
        *,
        C=1.0,
        l1_ratio=0.0,
        dual=False,
        tol=1e-4,
        fit_intercept=True,
        intercept_scaling=1,
        class_weight=None,
        random_state=None,
        solver="lbfgs",
        max_iter=100,
        verbose=0,
        warm_start=False,
        n_jobs=None,
    ):
        self.penalty = penalty
        self.C = C
        self.l1_ratio = l1_ratio
        self.dual = dual
        self.tol = tol
        self.fit_intercept = fit_intercept
        self.intercept_scaling = intercept_scaling
        self.class_weight = class_weight
        self.random_state = random_state
        self.solver = solver
        self.max_iter = max_iter
        self.verbose = verbose
        self.warm_start = warm_start
        self.n_jobs = n_jobs

    # TODO(callbacks): update/remove as more solvers get supported.
    def set_callbacks(self, *callbacks):
        """Set callbacks for the estimator.

        Parameters
        ----------
        *callbacks : callback instances
            The callbacks to set.

        Returns
        -------
        self : estimator instance
            The estimator instance itself.
        """
        if self.solver != "lbfgs":
            warnings.warn(
                f"Callbacks are only supported in {self.__class__.__name__} for"
                f" solver='lbfgs'. This estimator has solver='{self.solver}', callbacks"
                " will be ignored during the solving of the regression."
            )
        return super().set_callbacks(*callbacks)

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y, sample_weight=None):
        """
        Fit the model according to the given training data.

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape (n_samples, n_features)
            Training vector, where `n_samples` is the number of samples and
            `n_features` is the number of features.

        y : array-like of shape (n_samples,)
            Target vector relative to X.

        sample_weight : array-like of shape (n_samples,) default=None
            Array of weights that are assigned to individual samples.
            If not provided, then each sample is given unit weight.

            .. versionadded:: 0.17
               *sample_weight* support to LogisticRegression.

        Returns
        -------
        self
            Fitted estimator.

        Notes
        -----
        The SAGA solver supports both float64 and float32 bit arrays.
        """
        if self.penalty == "deprecated":
            if self.l1_ratio == 0 or self.l1_ratio is None:
                penalty = "l2"
                if self.l1_ratio is None:
                    warnings.warn(
                        (
                            "'l1_ratio=None' was deprecated in version 1.8 and will "
                            "trigger an error in 1.10. Use 0<=l1_ratio<=1 instead."
                        ),
                        FutureWarning,
                    )
            elif self.l1_ratio == 1:
                penalty = "l1"
            else:
                penalty = "elasticnet"
            if self.C == np.inf:
                penalty = None
        else:
            penalty = self.penalty
            warnings.warn(
                (
                    "'penalty' was deprecated in version 1.8 and will be removed in"
                    " 1.10. To avoid this warning, leave 'penalty' set to its default"
                    " value and use 'l1_ratio' or 'C' instead."
                    " Use l1_ratio=0 instead of penalty='l2',"
                    " l1_ratio=1 instead of penalty='l1',"
                    " l1_ratio set to a float between 0 and 1 instead of"
                    " penalty='elasticnet', and C=np.inf instead of penalty=None."
                ),
                FutureWarning,
            )

        solver = _check_solver(self.solver, penalty, self.dual)

        if penalty != "elasticnet" and (
            self.l1_ratio is not None and 0 < self.l1_ratio < 1
        ):
            warnings.warn(
                "l1_ratio parameter is only used when penalty is "
                "'elasticnet'. Got "
                "(penalty={})".format(penalty)
            )
        if (self.penalty == "l2" and self.l1_ratio != 0) or (
            self.penalty == "l1" and self.l1_ratio != 1
        ):
            warnings.warn(
                f"Inconsistent values: penalty={self.penalty} with "
                f"l1_ratio={self.l1_ratio}. penalty is deprecated. Please use "
                f"l1_ratio only."
            )
        if penalty == "elasticnet" and self.l1_ratio is None:
            raise ValueError("l1_ratio must be specified when penalty is elasticnet.")

        xp, _, device = get_namespace_and_device(X)
        sample_weight = move_to(sample_weight, xp=xp, device=device)
        xp_y, _ = get_namespace(y)

        if self.penalty is None:
            if self.C != 1.0:  # default value
                warnings.warn(
                    "Setting penalty=None will ignore the C and l1_ratio parameters"
                )
                # Note that check for l1_ratio is done right above
            C_ = xp.inf
            penalty = "l2"
        else:
            C_ = self.C

        msg = (
            "'n_jobs' has no effect since 1.8 and will be removed in 1.10. "
            f"You provided 'n_jobs={self.n_jobs}', please leave it unspecified."
        )
        if self.n_jobs is not None:
            warnings.warn(msg, category=FutureWarning)

        X, y = validate_data(
            self,
            X,
            y,
            accept_sparse="csc" if solver == "newton-cd" else "csr",
            dtype=[xp.float64, xp.float32],
            order="F" if solver == "newton-cd" else "C",
            accept_large_sparse=solver not in ("liblinear", "newton-cd", "sag", "saga"),
        )
        n_samples, n_features = X.shape
        check_classification_targets(y)
        le = LabelEncoder().fit(y)
        self.classes_ = le.classes_
        n_classes = size(self.classes_)
        if n_classes < 2:
            msg = (
                "This solver needs samples of at least 2 classes in the data, but the "
                f"data contains only one class: {self.classes_[0]}"
            )
            raise ValueError(msg)
        is_binary = n_classes == 2
        y_encoded = move_to(le.transform(y), xp=xp, device=device)
        y_encoded = xp.astype(y_encoded, X.dtype, copy=False)

        if sample_weight is not None or self.class_weight is not None:
            sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)

        if self.class_weight is not None:
            class_weight_ = compute_class_weight(
                self.class_weight,
                classes=self.classes_,
                y=y,
                sample_weight=sample_weight,
            )
            class_weight_ = xp.asarray(
                class_weight_[le.transform(y)], dtype=X.dtype, device=device
            )
            # Avoid overriding the input sample_weight.
            sample_weight = sample_weight * class_weight_

        # With lbfgs, the fit task will have a subtask even if max_iter is 0.
        # There's also always one extra empty subtask due to the scipy.optimize.minimize
        # callback logic.
        max_subtasks = max(self.max_iter, 1) + 1
        callback_ctx = self._init_callback_context(max_subtasks=max_subtasks)
        callback_metadata = (
            {"sample_weight": sample_weight} if sample_weight is not None else None
        )
        callback_ctx.call_on_fit_task_begin(
            estimator=self, X=X, y=y, metadata=callback_metadata
        )

        if solver == "liblinear":
            if not is_binary:
                raise ValueError(
                    "The 'liblinear' solver does not support multiclass classification"
                    " (n_classes >= 3). Either use another solver or wrap the "
                    "estimator in a OneVsRestClassifier to keep applying a "
                    "one-versus-rest scheme."
                )
            if np.max(X) > 1e30:
                raise ValueError(
                    "Using the 'liblinear' solver while X contains a maximum "
                    "value > 1e30 results in a frozen fit. Please choose another "
                    "solver or rescale the input X."
                )
            self.coef_, self.intercept_, self.n_iter_ = _fit_liblinear(
                X,
                2 * y_encoded - 1,  # liblinear requires target values -1, +1
                self.C,
                self.fit_intercept,
                self.intercept_scaling,
                None,  # class_weight
                penalty,
                self.dual,
                self.verbose,
                self.max_iter,
                self.tol,
                self.random_state,
                sample_weight=sample_weight,
            )
            return self

        if solver in ["sag", "saga"]:
            max_squared_sum = row_norms(X, squared=True).max()
        else:
            max_squared_sum = None

        if self.warm_start:
            warm_start_coef = getattr(self, "coef_", None)
        else:
            warm_start_coef = None
        if warm_start_coef is not None:
            warm_start_coef = move_to(warm_start_coef, xp=np, device="cpu")
            if self.fit_intercept:
                intercept_np = move_to(self.intercept_, xp=np, device="cpu")
                warm_start_coef = np.concatenate(
                    [warm_start_coef, intercept_np[:, None]], axis=1
                )

        # TODO: enable multi-threading if benchmarks show a positive effect,
        # see https://github.com/scikit-learn/scikit-learn/issues/32162
        n_threads = 1

        coefs, _, n_iter = _logistic_regression_path(
            X,
            y_encoded,
            classes=self.classes_,
            Cs=[C_],
            l1_ratio=self.l1_ratio,
            fit_intercept=self.fit_intercept,
            tol=self.tol,
            verbose=self.verbose,
            solver=solver,
            max_iter=self.max_iter,
            check_input=False,
            random_state=self.random_state,
            coef=warm_start_coef,
            penalty=penalty,
            max_squared_sum=max_squared_sum,
            sample_weight=sample_weight,
            n_threads=n_threads,
            callback_ctx=callback_ctx,
            estimator=self,
        )

        self.n_iter_ = xp.asarray(n_iter, dtype=xp.int32)

        self.coef_ = coefs[0, ...]
        if self.fit_intercept:
            if is_binary:
                self.intercept_ = self.coef_[-1:]
                self.coef_ = self.coef_[:-1][None, :]
            else:
                self.intercept_ = self.coef_[:, -1]
                self.coef_ = self.coef_[:, :-1]
        else:
            if is_binary:
                self.intercept_ = xp.zeros(1, dtype=X.dtype, device=device)
                self.coef_ = self.coef_[None, :]
            else:
                self.intercept_ = xp.zeros(n_classes, dtype=X.dtype, device=device)

        callback_ctx.call_on_fit_task_end(
            estimator=self,
            X=X,
            y=y,
            metadata=callback_metadata,
            reconstruction_attributes={},
        )

        return self

    def predict_proba(self, X):
        """
        Probability estimates.

        The returned estimates for all classes are ordered by the
        label of classes.

        For a multiclass / multinomial problem the softmax function is used to find
        the predicted probability of each class.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Vector to be scored, where `n_samples` is the number of samples and
            `n_features` is the number of features.

        Returns
        -------
        T : array-like of shape (n_samples, n_classes)
            Returns the probability of the sample for each class in the model,
            where classes are ordered as they are in ``self.classes_``.
        """
        check_is_fitted(self)
        check_same_namespace(X, self, attribute="coef_", method="predict_proba")

        is_binary = size(self.classes_) <= 2
        if is_binary:
            return super()._predict_proba_lr(X)
        else:
            decision_2d = self.decision_function(X)
            return softmax(decision_2d, copy=False)

    def predict_log_proba(self, X):
        """
        Predict logarithm of probability estimates.

        The returned estimates for all classes are ordered by the
        label of classes.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Vector to be scored, where `n_samples` is the number of samples and
            `n_features` is the number of features.

        Returns
        -------
        T : array-like of shape (n_samples, n_classes)
            Returns the log-probability of the sample for each class in the
            model, where classes are ordered as they are in ``self.classes_``.
        """
        check_same_namespace(X, self, attribute="coef_", method="predict_log_proba")
        xp, _ = get_namespace(X)
        return xp.log(self.predict_proba(X))

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        tags.input_tags.sparse = True
        tags.array_api_support = self.solver in ("lbfgs", "newton-cg")
        if self.solver == "liblinear":
            tags.classifier_tags.multi_class = False

        return tags


class LogisticRegressionCV(LogisticRegression, LinearClassifierMixin, BaseEstimator):
    """Logistic Regression CV (aka logit, MaxEnt) classifier.

    See glossary entry for :term:`cross-validation estimator`.

    This class implements regularized logistic regression with implicit cross
    validation for the penalty parameters `C` and `l1_ratio`, see
    :class:`LogisticRegression`, using a set of available solvers.

    The solvers 'lbfgs', 'newton-cg', 'newton-cholesky' and 'sag' support only L2
    regularization with primal formulation. The 'liblinear'
    solver supports both L1 and L2 regularization (but not both, i.e. elastic-net),
    with a dual formulation only for the L2 penalty. The Elastic-Net (combination of L1
    and L2) regularization is only supported by the 'saga' solver.

    For the grid of `Cs` values and `l1_ratios` values, the best hyperparameter
    is selected by the cross-validator
    :class:`~sklearn.model_selection.StratifiedKFold`, but it can be changed
    using the :term:`cv` parameter. All solvers except 'liblinear' can warm-start the
    coefficients (see :term:`Glossary<warm_start>`).

    Read more in the :ref:`User Guide <logistic_regression>`.

    Parameters
    ----------
    Cs : int or list of floats, default=10
        Each of the values in Cs describes the inverse of regularization
        strength. If Cs is as an int, then a grid of Cs values are chosen
        in a logarithmic scale between 1e-4 and 1e4.
        Like in support vector machines, smaller values specify stronger
        regularization.

    l1_ratios : array-like of shape (n_l1_ratios), default=None
        Floats between 0 and 1 passed as Elastic-Net mixing parameter (scaling between
        L1 and L2 penalties). For `l1_ratio = 0` the penalty is an L2 penalty. For
        `l1_ratio = 1` it is an L1 penalty. Any value between 0 and 1 gives
        an Elastic-Net penalty of the form `l1_ratio * L1 + (1 - l1_ratio) * L2`.
        All the values of the given array-like are tested by cross-validation and the
        one giving the best prediction score is used.

        .. warning::
           Certain values of `l1_ratios`, i.e. some penalties, may not work with some
           solvers. See the parameter `solver` below, to know the compatibility between
           the penalty and solver.

        .. deprecated:: 1.8
            `l1_ratios=None` is deprecated in 1.8 and will raise an error
            in version 1.10. Default value will change from `None` to `(0.0,)`
            in version 1.10.

    fit_intercept : bool, default=True
        Specifies if a constant (a.k.a. bias or intercept) should be
        added to the decision function.

    cv : int or cross-validation generator, default=None
        The default cross-validation generator used is Stratified K-Folds.
        If an integer is provided, it specifies the number of folds, `n_folds`, used.
        See the module :mod:`sklearn.model_selection` module for the
        list of possible cross-validation objects.

        .. versionchanged:: 0.22
            ``cv`` default value if None changed from 3-fold to 5-fold.

    dual : bool, default=False
        Dual (constrained) or primal (regularized, see also
        :ref:`this equation <regularized-logistic-loss>`) formulation. Dual formulation
        is only implemented for l2 penalty with liblinear solver. Prefer dual=False when
        n_samples > n_features.

    penalty : {'l1', 'l2', 'elasticnet'}, default='l2'
        Specify the norm of the penalty:

        - `'l2'`: add an L2 penalty term (used by default);
        - `'l1'`: add an L1 penalty term;
        - `'elasticnet'`: both L1 and L2 penalty terms are added.

        .. warning::
           Some penalties may not work with some solvers. See the parameter
           `solver` below, to know the compatibility between the penalty and
           solver.

        .. deprecated:: 1.8
           `penalty` was deprecated in version 1.8 and will be removed in 1.10.
           Use `l1_ratio` and `C` instead. `l1_ratio=0` for `penalty='l2'`,
           `l1_ratio=1` for `penalty='l1'`, `l1_ratio` set to any float between 0 and 1
           for `penalty='elasticnet'`, and `C=np.inf` for `penalty=None`.

    scoring : str or callable, default=None
        The scoring method to use for cross-validation. Options:

        - str: see :ref:`scoring_string_names` for options.
        - callable: a scorer callable object (e.g., function) with signature
          ``scorer(estimator, X, y)``. See :ref:`scoring_callable` for details.
        - `None`: :ref:`accuracy <accuracy_score>` is used.

        .. versionchanged:: 1.11
           The default will change from None, i.e. accuracy, to 'neg_log_loss' in
           version 1.11.

    solver : {'lbfgs', 'liblinear', 'newton-cd-gram', 'newton-cg', 'newton-cholesky', \
            'sag', 'saga'}, default='lbfgs'

        Algorithm to use in the optimization problem. Default is 'lbfgs'.
        To choose a solver, you might want to consider the following aspects:

        - 'lbfgs' is a good default solver because it works reasonably well for a wide
          class of problems.
        - For :term:`multiclass` problems (`n_classes >= 3`), all solvers except
          'liblinear' minimize the full multinomial loss, 'liblinear' will raise an
          error.
        - 'newton-cholesky' is a good choice for
          `n_samples` >> `n_features * n_classes`, especially with one-hot encoded
          categorical features with rare categories. Be aware that the memory usage
          of this solver has a quadratic dependency on `n_features * n_classes`
          because it explicitly computes the full Hessian matrix.
        - For small datasets, 'liblinear' is a good choice, whereas 'sag'
          and 'saga' are faster for large ones;
        - 'liblinear' might be slower in :class:`LogisticRegressionCV`
          because it does not handle warm-starting.
        - 'liblinear' can only handle binary classification by default. To apply a
          one-versus-rest scheme for the multiclass setting one can wrap it with the
          :class:`~sklearn.multiclass.OneVsRestClassifier`.

        .. warning::
           The choice of the algorithm depends on the penalty (`l1_ratio=0` for
           L2-penalty, `l1_ratio=1` for L1-penalty and `0 < l1_ratio < 1` for
           Elastic-Net) chosen and on (multinomial) multiclass support:

           ================= ======================== ======================
           solver            l1_ratio                 multinomial multiclass
           ================= ======================== ======================
           'lbfgs'           l1_ratio=0               yes
           'liblinear'       l1_ratio=1 or l1_ratio=0 no
           'newton-cd'       0<=l1_ratio<=1           no
           'newton-cd-gram'  0<=l1_ratio<=1           yes
           'newton-cg'       l1_ratio=0               yes
           'newton-cholesky' l1_ratio=0               yes
           'sag'             l1_ratio=0               yes
           'saga'            0<=l1_ratio<=1           yes
           ================= ======================== ======================

        .. note::
           'sag' and 'saga' fast convergence is only guaranteed on features
           with approximately the same scale. You can preprocess the data with
           a scaler from :mod:`sklearn.preprocessing`.

        .. versionadded:: 0.17
           Stochastic Average Gradient (SAG) descent solver. Multinomial support in
           version 0.18.
        .. versionadded:: 0.19
           SAGA solver.
        .. versionadded:: 1.2
           newton-cholesky solver. Multinomial support in version 1.6.

    tol : float, default=1e-4
        Tolerance for stopping criteria.

    max_iter : int, default=100
        Maximum number of iterations of the optimization algorithm.

    class_weight : dict or 'balanced', default=None
        Weights associated with classes in the form ``{class_label: weight}``.
        If not given, all classes are supposed to have weight one.

        The "balanced" mode uses the values of y to automatically adjust
        weights inversely proportional to class frequencies in the input data
        as ``n_samples / (n_classes * np.bincount(y))``.

        Note that these weights will be multiplied with sample_weight (passed
        through the fit method) if sample_weight is specified.

        .. versionadded:: 0.17
           class_weight == 'balanced'

    n_jobs : int, default=None
        Number of CPU cores used during the cross-validation loop.
        ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
        ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
        for more details.

    verbose : int, default=0
        For the 'liblinear', 'sag' and 'lbfgs' solvers set verbose to any
        positive number for verbosity.

    refit : bool, default=True
        If set to True, the scores are averaged across all folds, and the
        coefs and the C that corresponds to the best score is taken, and a
        final refit is done using these parameters.
        Otherwise the coefs, intercepts and C that correspond to the
        best scores across folds are averaged.

    intercept_scaling : float, default=1
        Useful only when the solver `liblinear` is used
        and `self.fit_intercept` is set to `True`. In this case, `x` becomes
        `[x, self.intercept_scaling]`,
        i.e. a "synthetic" feature with constant value equal to
        `intercept_scaling` is appended to the instance vector.
        The intercept becomes
        ``intercept_scaling * synthetic_feature_weight``.

        .. note::
            The synthetic feature weight is subject to L1 or L2
            regularization as all other features.
            To lessen the effect of regularization on synthetic feature weight
            (and therefore on the intercept) `intercept_scaling` has to be increased.

    random_state : int, RandomState instance, default=None
        Only used for `solver` == 'sag', 'saga' or 'liblinear' to shuffle the
        data. It has no effect on the other solvers.
        Note that this only applies to the solver and not the cross-validation
        generator. See :term:`Glossary <random_state>` for details.

    use_legacy_attributes : bool, default=True
        If True, use legacy values for attributes:

        - `C_` is an ndarray of shape (n_classes,) with the same value repeated
        - `l1_ratio_` is an ndarray of shape (n_classes,) with the same value repeated
        - `coefs_paths_` is a dict with class labels as keys and ndarrays as values
        - `scores_` is a dict with class labels as keys and ndarrays as values
        - `n_iter_` is an ndarray of shape (1, n_folds, n_cs) or similar

        If False, use new values for attributes:

        - `C_` is a float
        - `l1_ratio_` is a float
        - `coefs_paths_` is an ndarray of shape
          (n_folds, n_l1_ratios, n_cs, n_classes, n_features)
          For binary problems (n_classes=2), the 2nd last dimension is 1.
        - `scores_` is an ndarray of shape (n_folds, n_l1_ratios, n_cs)
        - `n_iter_` is an ndarray of shape (n_folds, n_l1_ratios, n_cs)

        .. versionchanged:: 1.10
           The default will change from True to False in version 1.10.
        .. deprecated:: 1.10
           `use_legacy_attributes` will be deprecated in version 1.10 and be removed in
           1.12.

    Attributes
    ----------
    classes_ : ndarray of shape (n_classes, )
        A list of class labels known to the classifier.

    coef_ : ndarray of shape (1, n_features) or (n_classes, n_features)
        Coefficient of the features in the decision function.

        `coef_` is of shape (1, n_features) when the given problem
        is binary.

    intercept_ : ndarray of shape (1,) or (n_classes,)
        Intercept (a.k.a. bias) added to the decision function.

        If `fit_intercept` is set to False, the intercept is set to zero.
        `intercept_` is of shape (1,) when the problem is binary.

    Cs_ : ndarray of shape (n_cs)
        Array of C i.e. inverse of regularization parameter values used
        for cross-validation.

    l1_ratios_ : ndarray of shape (n_l1_ratios)
        Array of l1_ratios used for cross-validation. If l1_ratios=None is used
        (i.e. penalty is not 'elasticnet'), this is set to ``[None]``

    coefs_paths_ : dict of ndarray of shape (n_folds, n_cs, n_dof) or \
            (n_folds, n_cs, n_l1_ratios, n_dof)
        A dict with classes as the keys, and the path of coefficients obtained
        during cross-validating across each fold (`n_folds`) and then across each Cs
        (`n_cs`).
        The size of the coefficients is the number of degrees of freedom (`n_dof`),
        i.e. without intercept `n_dof=n_features` and with intercept
        `n_dof=n_features+1`.
        If `penalty='elasticnet'`, there is an additional dimension for the number of
        l1_ratio values (`n_l1_ratios`), which gives a shape of
        ``(n_folds, n_cs, n_l1_ratios_, n_dof)``.
        See also parameter `use_legacy_attributes`.

    scores_ : dict
        A dict with classes as the keys, and the values as the
        grid of scores obtained during cross-validating each fold.
        The same score is repeated across all classes. Each dict value
        has shape ``(n_folds, n_cs)`` or ``(n_folds, n_cs, n_l1_ratios)`` if
        ``penalty='elasticnet'``.
        See also parameter `use_legacy_attributes`.

    C_ : ndarray of shape (n_classes,) or (1,)
        The value of C that maps to the best score, repeated n_classes times.
        If refit is set to False, the best C is the average of the
        C's that correspond to the best score for each fold.
        `C_` is of shape (1,) when the problem is binary.
        See also parameter `use_legacy_attributes`.

    l1_ratio_ : ndarray of shape (n_classes,) or (n_classes - 1,)
        The value of l1_ratio that maps to the best score, repeated n_classes times.
        If refit is set to False, the best l1_ratio is the average of the
        l1_ratio's that correspond to the best score for each fold.
        `l1_ratio_` is of shape (1,) when the problem is binary.
        See also parameter `use_legacy_attributes`.

    n_iter_ : ndarray of shape (1, n_folds, n_cs) or (1, n_folds, n_cs, n_l1_ratios)
        Actual number of iterations for all classes, folds and Cs.
        If `penalty='elasticnet'`, the shape is `(1, n_folds, n_cs, n_l1_ratios)`.
        See also parameter `use_legacy_attributes`.

    n_features_in_ : int
        Number of features seen during :term:`fit`.

        .. versionadded:: 0.24

    feature_names_in_ : ndarray of shape (`n_features_in_`,)
        Names of features seen during :term:`fit`. Defined only when `X`
        has feature names that are all strings.

        .. versionadded:: 1.0

    See Also
    --------
    LogisticRegression : Logistic regression without tuning the
        hyperparameter `C`.

    Examples
    --------
    >>> from sklearn.datasets import load_iris
    >>> from sklearn.linear_model import LogisticRegressionCV
    >>> X, y = load_iris(return_X_y=True)
    >>> clf = LogisticRegressionCV(
    ...     cv=5, random_state=0,
    ...     use_legacy_attributes=False,
    ...     l1_ratios=(0,),
    ...     scoring="neg_log_loss",
    ... ).fit(X, y)
    >>> clf.predict(X[:2, :])
    array([0, 0])
    >>> clf.predict_proba(X[:2, :]).shape
    (2, 3)
    >>> clf.score(X, y)
    -0.041...
    """

    # TODO(1.11): remove this when sample_weight as positional arg is removed
    # from the `score` signature
    __metadata_request__score = {"sample_weight": metadata_routing.UNUSED}
    _parameter_constraints: dict = {**LogisticRegression._parameter_constraints}

    for param in ["C", "warm_start", "l1_ratio"]:
        _parameter_constraints.pop(param)

    _parameter_constraints.update(
        {
            "Cs": [Interval(Integral, 1, None, closed="left"), "array-like"],
            "l1_ratios": ["array-like", None, Hidden(StrOptions({"warn"}))],
            "cv": ["cv_object"],
            "scoring": [
                StrOptions(set(get_scorer_names())),
                callable,
                None,
                Hidden(StrOptions({"warn"})),
            ],
            "refit": ["boolean"],
            "penalty": [
                StrOptions({"l1", "l2", "elasticnet"}),
                Hidden(StrOptions({"deprecated"})),
            ],
            "use_legacy_attributes": ["boolean", Hidden(StrOptions({"warn"}))],
        }
    )

    def __init__(
        self,
        *,
        Cs=10,
        l1_ratios="warn",
        fit_intercept=True,
        cv=None,
        dual=False,
        penalty="deprecated",
        scoring="warn",
        solver="lbfgs",
        tol=1e-4,
        max_iter=100,
        class_weight=None,
        n_jobs=None,
        verbose=0,
        refit=True,
        intercept_scaling=1.0,
        random_state=None,
        use_legacy_attributes="warn",
    ):
        self.Cs = Cs
        self.l1_ratios = l1_ratios
        self.fit_intercept = fit_intercept
        self.cv = cv
        self.dual = dual
        self.penalty = penalty
        self.scoring = scoring
        self.tol = tol
        self.max_iter = max_iter
        self.class_weight = class_weight
        self.n_jobs = n_jobs
        self.verbose = verbose
        self.solver = solver
        self.refit = refit
        self.intercept_scaling = intercept_scaling
        self.random_state = random_state
        self.use_legacy_attributes = use_legacy_attributes

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y, sample_weight=None, **params):
        """Fit the model according to the given training data.

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape (n_samples, n_features)
            Training vector, where `n_samples` is the number of samples and
            `n_features` is the number of features.

        y : array-like of shape (n_samples,)
            Target vector relative to X.

        sample_weight : array-like of shape (n_samples,) default=None
            Array of weights that are assigned to individual samples.
            If not provided, then each sample is given unit weight.

        **params : dict
            Parameters to pass to the underlying splitter and scorer.

            .. versionadded:: 1.4

        Returns
        -------
        self : object
            Fitted LogisticRegressionCV estimator.
        """
        _raise_for_params(params, self, "fit")

        if isinstance(self.l1_ratios, str) and self.l1_ratios == "warn":
            l1_ratios = None
            warnings.warn(
                (
                    "The default value for l1_ratios will change from None to (0.0,) "
                    "in version 1.10. From version 1.10 onwards, only array-like "
                    "with values in [0, 1] will be allowed, None will be forbidden. "
                    "To avoid this warning, explicitly set a value, "
                    "e.g. l1_ratios=(0,)."
                ),
                FutureWarning,
            )
        else:
            l1_ratios = self.l1_ratios

        if self.penalty == "deprecated":
            if self.l1_ratios is None:
                warnings.warn(
                    (
                        "'l1_ratios=None' was deprecated in version 1.8 and will "
                        "trigger an error in 1.10. Use an array-like with values"
                        "in [0, 1] instead."
                    ),
                    FutureWarning,
                )
            if np.all(np.asarray(l1_ratios) == 0) or l1_ratios is None:
                penalty = "l2"
            elif np.all(np.asarray(l1_ratios) == 1):
                penalty = "l1"
            else:
                penalty = "elasticnet"
        else:
            penalty = self.penalty
            warnings.warn(
                (
                    "'penalty' was deprecated in version 1.8 and will be removed in"
                    " 1.10. To avoid this warning, leave 'penalty' set to its default"
                    " value and use 'l1_ratios' and 'Cs' instead."
                    " Use l1_ratios=(0,) instead of penalty='l2',"
                    " l1_ratios=(1,) instead of penalty='l1',"
                    " l1_ratios set to floats between 0 and 1 instead of"
                    " penalty='elasticnet', and Cs=(np.inf,) instead of penalty=None."
                ),
                FutureWarning,
            )

        if self.scoring == "warn":
            warnings.warn(
                "The default value of the parameter 'scoring' will change from None, "
                "i.e. accuracy, to 'neg_log_loss' in version 1.11. To silence this "
                "warning, explicitly set the scoring parameter: "
                "scoring='neg_log_loss' for the new, scoring='accuracy' or "
                "scoring=None for the old default.",
                FutureWarning,
            )
            scoring = None
        else:
            scoring = self.scoring

        if self.use_legacy_attributes == "warn":
            warnings.warn(
                f"The fitted attributes of {self.__class__.__name__} will be "
                "simplified in scikit-learn 1.10 to remove redundancy. Set"
                "`use_legacy_attributes=False` to enable the new behavior now, or "
                "set it to `True` to silence this warning during the transition period "
                "while keeping the deprecated behavior for the time being. The default "
                "value of use_legacy_attributes will change from True to False in "
                f"scikit-learn 1.10. See the docstring of {self.__class__.__name__} "
                "for more details.",
                FutureWarning,
            )
            use_legacy_attributes = True
        else:
            use_legacy_attributes = self.use_legacy_attributes

        solver = _check_solver(self.solver, penalty, self.dual)

        if penalty == "elasticnet":
            if (
                l1_ratios is None
                or len(l1_ratios) == 0
                or any(
                    (
                        not isinstance(l1_ratio, numbers.Number)
                        or l1_ratio < 0
                        or l1_ratio > 1
                    )
                    for l1_ratio in l1_ratios
                )
            ):
                raise ValueError(
                    "l1_ratios must be an array-like of numbers between "
                    "0 and 1; got (l1_ratios=%r)" % l1_ratios
                )
            l1_ratios_ = l1_ratios
        else:
            if l1_ratios is not None and self.penalty != "deprecated":
                warnings.warn(
                    "l1_ratios parameter is only used when penalty "
                    "is 'elasticnet'. Got (penalty={})".format(penalty)
                )

            if l1_ratios is None:
                l1_ratios_ = [None]
            else:
                l1_ratios_ = l1_ratios

        xp, _, device = get_namespace_and_device(X)
        sample_weight = move_to(sample_weight, xp=xp, device=device)
        X, y = validate_data(
            self,
            X,
            y,
            accept_sparse="csr",  # CV will index data on first dimension
            dtype=[xp.float64, xp.float32],
            order="C",  # CV will index data on first dimension
            accept_large_sparse=solver not in ("liblinear", "newton-cd", "sag", "saga"),
        )
        n_samples, n_features = X.shape
        check_classification_targets(y)
        le = LabelEncoder().fit(y)
        self.classes_ = le.classes_
        n_classes = size(self.classes_)
        if n_classes < 2:
            msg = (
                "This solver needs samples of at least 2 classes in the data, but the "
                f"data contains only one class: {self.classes_[0]}"
            )
            raise ValueError(msg)
        is_binary = n_classes == 2
        y_encoded = move_to(le.transform(y), xp=xp, device=device)
        y_encoded = xp.astype(y_encoded, X.dtype, copy=False)

        if sample_weight is not None or self.class_weight is not None:
            sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)

        if self.class_weight is not None:
            class_weight_ = compute_class_weight(
                self.class_weight,
                classes=self.classes_,
                y=y,
                sample_weight=sample_weight,
            )
            class_weight_ = xp.asarray(
                class_weight_[le.transform(y)], dtype=X.dtype, device=device
            )
            # Avoid overriding the input sample_weight.
            sample_weight = sample_weight * class_weight_

        # The original class labels
        class_labels = self.classes_

        if is_binary:
            n_classes = 1
            class_labels = class_labels[1:]

        class_labels = move_to(class_labels, xp=np, device="cpu")

        if n_classes >= 3 and solver == "newton-cd" and sparse.issparse(X):
            raise ValueError(
                f"Solver 'newton-cd' does not support sparse X for multiclass settings"
                f" (n_classes >= 3); got {n_classes=}."
            )

        if solver in ["sag", "saga"]:
            max_squared_sum = row_norms(X, squared=True).max()
        else:
            max_squared_sum = None

        if _routing_enabled():
            routed_params = process_routing(
                self,
                "fit",
                sample_weight=sample_weight,
                **params,
            )
        else:
            score_kwargs = dict(params)
            if sample_weight is not None:
                score_kwargs["sample_weight"] = sample_weight
            routed_params = _manual_routing(
                {"splitter": {}, "scorer": {"score": score_kwargs}}
            )

        # init cross-validation generator
        cv = check_cv(self.cv, y, classifier=True)
        folds = list(cv.split(X, y_encoded, **routed_params.splitter.split))

        path_func = delayed(_log_reg_scoring_path)

        # The SAG solver releases the GIL so it's more efficient to use
        # threads for this solver.
        if self.solver in ["sag", "saga"]:
            prefer = "threads"
        else:
            prefer = "processes"

        fold_coefs_ = Parallel(n_jobs=self.n_jobs, verbose=self.verbose, prefer=prefer)(
            path_func(
                X,
                y_encoded,
                train,
                test,
                classes=self.classes_,
                Cs=self.Cs,
                fit_intercept=self.fit_intercept,
                penalty=penalty,
                dual=self.dual,
                solver=solver,
                tol=self.tol,
                max_iter=self.max_iter,
                verbose=self.verbose,
                scoring=scoring,
                intercept_scaling=self.intercept_scaling,
                random_state=self.random_state,
                max_squared_sum=max_squared_sum,
                sample_weight=sample_weight,
                l1_ratio=l1_ratio,
                score_params=routed_params.scorer.score,
            )
            for train, test in folds
            for l1_ratio in l1_ratios_
        )

        # fold_coefs_ is a list and would have shape (n_folds * n_l1_ratios, ..)
        # After reshaping,
        # - coefs_paths is of shape (n_classes, n_folds, n_Cs, n_l1_ratios, n_features)
        # - scores is of shape (n_classes, n_folds, n_Cs, n_l1_ratios)
        # - n_iter is of shape (1, n_folds, n_Cs, n_l1_ratios)
        coefs_paths, Cs, scores, n_iter_ = zip(*fold_coefs_)
        self.Cs_ = Cs[0]  # the same for all folds and l1_ratios
        n_folds = len(folds)
        n_cs = size(self.Cs_)
        n_l1_ratios = len(l1_ratios_)
        coefs_paths = xp.stack(coefs_paths)
        n_iter_ = xp.stack(n_iter_)
        if is_binary:
            coefs_paths = xp.reshape(coefs_paths, (n_folds, n_l1_ratios, n_cs, -1))
            # coefs_paths.shape = (n_folds, n_l1_ratios, n_Cs, n_features)
            coefs_paths = _swapaxes(coefs_paths, 1, 2, xp=xp)[None, ...]
        else:
            coefs_paths = xp.reshape(
                coefs_paths, (n_folds, n_l1_ratios, n_cs, n_classes, -1)
            )
            # coefs_paths.shape = (n_folds, n_l1_ratios, n_Cs, n_classes, n_features)
            coefs_paths = xp.moveaxis(coefs_paths, (0, 1, 3), (1, 3, 0))
        # n_iter_.shape = (n_folds, n_l1_ratios, n_Cs)
        n_iter_ = xp.reshape(n_iter_, (n_folds, n_l1_ratios, n_cs))
        self.n_iter_ = _swapaxes(n_iter_, 1, 2, xp=xp)[None, ...]
        # scores.shape = (n_folds, n_l1_ratios, n_Cs)
        scores = xp.stack(scores)
        scores = xp.reshape(scores, (n_folds, n_l1_ratios, n_cs))
        scores = _swapaxes(scores, 1, 2, xp=xp)[None, ...]
        # repeat same scores across all classes
        scores = xp.tile(scores, (n_classes, 1, 1, 1))
        self.scores_ = {class_labels[i]: scores[i, ...] for i in range(n_classes)}
        self.coefs_paths_ = {
            class_labels[i]: coefs_paths[i, ...] for i in range(n_classes)
        }

        self.C_ = list()
        self.l1_ratio_ = list()
        self.coef_ = xp.empty((n_classes, n_features), device=device, dtype=X.dtype)
        self.intercept_ = xp.zeros(n_classes, device=device, dtype=X.dtype)

        # All scores are the same across classes
        scores = self.scores_[class_labels[0]]

        if self.refit:
            # best_index over folds
            scores_sum = xp.sum(scores, axis=0)  # shape (n_cs, n_l1_ratios)
            best_index = _unravel_index(xp.argmax(scores_sum), scores_sum.shape)
            best_index_int = [int(idx) for idx in best_index]
            C_ = self.Cs_[best_index_int[0]]
            self.C_.append(C_)

            l1_ratio_ = l1_ratios_[best_index_int[1]]
            self.l1_ratio_.append(l1_ratio_)

            if is_binary:
                coef_init = xp.mean(coefs_paths[0, :, *best_index_int, :], axis=0)
            else:
                coef_init = xp.mean(coefs_paths[:, :, *best_index_int, :], axis=1)

            if solver == "lbfgs":
                coef_init = move_to(coef_init, xp=np, device="cpu")

            if solver == "newton-cd":
                if sparse.issparse(X):
                    X = X.tocsc()
                else:
                    X = np.asfortranarray(X)

            # Note that y is label encoded
            w, _, _ = _logistic_regression_path(
                X,
                y_encoded,
                classes=self.classes_,
                Cs=[float(C_)],
                solver=solver,
                fit_intercept=self.fit_intercept,
                coef=coef_init,
                max_iter=self.max_iter,
                tol=self.tol,
                penalty=penalty,
                verbose=max(0, self.verbose - 1),
                random_state=self.random_state,
                check_input=False,
                max_squared_sum=max_squared_sum,
                sample_weight=sample_weight,
                l1_ratio=l1_ratio_,
            )
            w = w[0, ...]

        else:
            # Take the best scores across every fold and the average of
            # all coefficients corresponding to the best scores.
            n_folds, n_cs, n_l1_ratios = scores.shape
            scores = xp.reshape(scores, (n_folds, -1))  # (n_folds, n_cs * n_l1_ratios)
            best_indices = xp.argmax(scores, axis=1)  # (n_folds,)
            best_indices = _unravel_index(best_indices, (n_cs, n_l1_ratios))
            best_indices = list(
                (int(i), int(j))
                for i, j in zip(*best_indices)  # (n_folds, 2)
            )
            # each row of best_indices has the 2 indices for Cs and l1_ratios
            if is_binary:
                fold_coef_paths = xp.stack(
                    [coefs_paths[0, i, *best_indices[i], :] for i in range(n_folds)]
                )
            else:
                fold_coef_paths = xp.stack(
                    [
                        coefs_paths[:, i, best_indices[i][0], best_indices[i][1], :]
                        for i in range(n_folds)
                    ]
                )
            w = xp.mean(fold_coef_paths, axis=0)

            best_indices = xp.asarray(best_indices, device=device)
            best_indices_C = best_indices[:, 0]
            self.C_.append(xp.mean(self.Cs_[best_indices_C]))

            if penalty == "elasticnet":
                best_indices_l1 = best_indices[:, 1]
                self.l1_ratio_.append(xp.mean(l1_ratios_[best_indices_l1]))
            else:
                self.l1_ratio_.append(0.0)

        self.C_ = xp.stack(self.C_)
        if is_binary:
            self.coef_ = w[:, :n_features] if w.ndim == 2 else w[:n_features][None, :]
            if self.fit_intercept:
                self.intercept_[0] = w[0, -1] if w.ndim == 2 else w[-1]
        else:
            self.C_ = xp.tile(self.C_, (n_classes,))
            self.l1_ratio_ = np.tile(self.l1_ratio_, n_classes)
            self.coef_ = w[:, :n_features]
            if self.fit_intercept:
                self.intercept_ = w[:, -1]

        self.C_ = xp.asarray(self.C_)
        self.l1_ratio_ = np.asarray(self.l1_ratio_)
        self.l1_ratios_ = np.asarray(l1_ratios_)
        if l1_ratios is None:
            # if elasticnet was not used, remove the l1_ratios dimension of some
            # attributes
            for cls, coefs_path in self.coefs_paths_.items():
                self.coefs_paths_[cls] = coefs_path[:, :, 0, :]
            for cls, score in self.scores_.items():
                self.scores_[cls] = score[:, :, 0]
            self.n_iter_ = self.n_iter_[:, :, :, 0]

        if not use_legacy_attributes:
            n_dof = X.shape[1] + int(self.fit_intercept)
            self.C_ = float(self.C_[0])
            newpaths = xp.concat(tuple(self.coefs_paths_.values()), axis=0)
            newscores = self.scores_[class_labels[0]]  # same for all classes
            newniter = self.n_iter_[0, ...]
            if l1_ratios is None:
                if n_classes <= 2:
                    newpaths = xp.reshape(newpaths, (1, n_folds, n_cs, 1, n_dof))
                else:
                    newpaths = xp.reshape(
                        newpaths, (n_classes, n_folds, n_cs, 1, n_dof)
                    )
                newscores = xp.reshape(newscores, (n_folds, n_cs, 1))
                newniter = xp.reshape(newniter, (n_folds, n_cs, 1))
                if self.penalty == "l1":
                    self.l1_ratio_ = 1.0
                else:
                    self.l1_ratio_ = 0.0
            else:
                n_l1_ratios = self.l1_ratios_.shape[0]
                self.l1_ratio_ = float(self.l1_ratio_[0])
                if n_classes <= 2:
                    newpaths = xp.reshape(
                        newpaths, (1, n_folds, n_cs, n_l1_ratios, n_dof)
                    )
                else:
                    newpaths = xp.reshape(
                        newpaths, (n_classes, n_folds, n_cs, n_l1_ratios, n_dof)
                    )
            # newpaths.shape = (n_classes, n_folds, n_cs, n_l1_ratios, n_dof)
            # self.coefs_paths_.shape should be
            # (n_folds, n_l1_ratios, n_cs, n_classes, n_dof)
            self.coefs_paths_ = xp.moveaxis(newpaths, (0, 1, 3), (3, 0, 1))
            # newscores.shape = (n_folds, n_cs, n_l1_ratios)
            # self.scores_.shape should be (n_folds, n_l1_ratios, n_cs)
            self.scores_ = xp.moveaxis(newscores, (1, 2), (2, 1))
            self.n_iter_ = xp.moveaxis(newniter, (1, 2), (2, 1))

        return self

    # TODO(1.11): remove this decorator along with `sample_weight` from the `score`
    #  signature
    @_deprecate_positional_args(version="1.11")
    def score(self, X, y, *, sample_weight=None, **score_params):
        """Score using the `scoring` option on the given test data and labels.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Test samples.

        y : array-like of shape (n_samples,)
            True labels for X.

        sample_weight : array-like of shape (n_samples,), default=None
            Sample weights.

            .. deprecated:: 1.9
              `sample_weight` needs to be passed as a keyword argument and not as a
              positional argument.

        **score_params : dict
            Parameters to pass to the `score` method of the underlying scorer.

            .. versionadded:: 1.4

        Returns
        -------
        score : float
            Score of self.predict(X) w.r.t. y.
        """
        # TODO(1.11): for backwards compatibility, when `sample_weight` becomes a part
        # of **score_params, it should be skipped in the following check so that people
        # can still pass it w/o metadata routing enabled.
        _raise_for_params(score_params, self, "score")
        # TODO(1.11): remove this when sample_weight is removed from the `score`
        # signature
        if sample_weight is not None:
            score_params["sample_weight"] = sample_weight

        scoring = self._get_scorer()
        if _routing_enabled():
            routed_params = process_routing(
                self,
                "score",
                **score_params,
            )
        else:
            sw = (
                {"sample_weight": score_params["sample_weight"]}
                if score_params.get("sample_weight") is not None
                else {}
            )
            routed_params = _manual_routing({"scorer": {"score": sw}})

        return scoring(
            self,
            X,
            y,
            **routed_params.scorer.score,
        )

    def get_metadata_routing(self):
        """Get metadata routing of this object.

        Please check :ref:`User Guide <metadata_routing>` on how the routing
        mechanism works.

        .. versionadded:: 1.4

        Returns
        -------
        routing : MetadataRouter
            A :class:`~sklearn.utils.metadata_routing.MetadataRouter` encapsulating
            routing information.
        """

        router = (
            MetadataRouter(owner=self)
            .add_self_request(self)
            .add(
                splitter=self.cv,
                method_mapping=MethodMapping().add(caller="fit", callee="split"),
            )
            .add(
                scorer=self._get_scorer(),
                method_mapping=MethodMapping()
                .add(caller="score", callee="score")
                .add(caller="fit", callee="score"),
            )
        )
        return router

    def _get_scorer(self):
        """Get the scorer based on the scoring method specified.
        The default scoring method is `accuracy`.
        """
        if self.scoring == "warn":  # TODO(1.11): remove
            return get_scorer("accuracy")
        scoring = self.scoring or "accuracy"
        return get_scorer(scoring)

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        tags.input_tags.sparse = True
        tags.array_api_support = self.solver in ("lbfgs", "newton-cg")
        return tags