sklearn/linear_model/_glm/glm.py - main - githubmirror/scikit-learn

sklearn/linear_model/_glm/glm.py
1 217 строк45 KB
Christian Lorentzen
ENH Newton-CD part 7: add NewtonCDSolver (#34561)
07 авг 2026, 19:08
07 авг 2026, 19:08c74a434
# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

"""
Generalized Linear Models with Exponential Dispersion Family
"""

from numbers import Integral, Real

import numpy as np
import scipy.optimize

from sklearn._loss.loss import (
    HalfGammaLoss,
    HalfPoissonLoss,
    HalfPoissonLossArrayAPI,
    HalfSquaredError,
    HalfTweedieLoss,
    HalfTweedieLossIdentity,
)
from sklearn.base import BaseEstimator, RegressorMixin, _fit_context
from sklearn.linear_model._glm._newton_solver import (
    NewtonCDGramSolver,
    NewtonCDSolver,
    NewtonCholeskySolver,
    NewtonSolver,
)
from sklearn.linear_model._linear_loss import LinearModelLoss
from sklearn.utils import check_array
from sklearn.utils._array_api import (
    _average,
    _is_numpy_namespace,
    _matching_numpy_dtype,
    get_namespace,
    get_namespace_and_device,
    move_to,
)
from sklearn.utils._openmp_helpers import _openmp_effective_n_threads
from sklearn.utils._param_validation import Hidden, Interval, StrOptions
from sklearn.utils.fixes import _get_additional_lbfgs_options_dict
from sklearn.utils.optimize import _check_optimize_result, _newton_cg
from sklearn.utils.validation import (
    _check_sample_weight,
    check_is_fitted,
    validate_data,
)


class _GeneralizedLinearRegressor(RegressorMixin, BaseEstimator):
    """Regression via a penalized Generalized Linear Model (GLM).

    GLMs based on a reproductive Exponential Dispersion Model (EDM) aim at fitting and
    predicting the mean of the target y as y_pred=h(X*w) with coefficients w.
    Therefore, the fit minimizes the following objective function with L2 priors as
    regularizer::

        1/(2*sum(s_i)) * sum(s_i * deviance(y_i, h(x_i*w))
        + alpha * l1_ratio ||w||_1 + 1/2 * alpha * (1 - l1_ratio) * ||w||_2^2

    with inverse link function h, s=sample_weight and per observation (unit) deviance
    deviance(y_i, h(x_i*w)). Note that for an EDM, 1/2 * deviance is the negative
    log-likelihood up to a constant (in w) term.
    The parameter ``alpha`` corresponds to the lambda parameter in glmnet.

    Instead of implementing the EDM family and a link function separately, we directly
    use the loss functions `from sklearn._loss` which have the link functions included
    in them for performance reasons. We pick the loss functions that implement
    (1/2 times) EDM deviances.

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

    .. versionadded:: 0.23

    Parameters
    ----------
    alpha : float, default=1
        Constant that multiplies the penalty terms and determines the
        regularization strength. ``alpha = 0`` is equivalent to unpenalized
        GLMs. In this case, the design matrix `X` must have full column rank
        (no collinearities).
        Values of `alpha` must be in the range `[0.0, inf)`.

    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, do not work with some
           solvers. See the parameter `solver` below, to know the compatibility between
           the penalty and solver.

        .. versionadded:: 1.10

    fit_intercept : bool, default=True
        Specifies if a constant (a.k.a. bias or intercept) should be
        added to the linear predictor (`X @ coef + intercept`).

    solver : {'lbfgs', 'newton-cg', 'newton-cholesky'}, default='lbfgs'
        Algorithm to use in the optimization problem:

        'lbfgs'
            Calls scipy's L-BFGS-B optimizer.

        'newton-cd'
            Uses Newton-Raphson steps in an iterated reweighted least squares fashion:
            The normal equations are cast as a weighted least squares problem with
            elastic-net penalty. The inner solver then uses a coordinate descent based
            solver. This way the full Hessian is used but never explicitly constructed.
            It can solve for all values of `l1_ratio`.
            This solver is a good choice for `n_features` > `n_samples`.

            .. versionadded:: 1.10

        'newton-cd-gram'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner coordinate descent based
            solver that uses the full Hessian/Gram matrix. It can solve for all
            values of `l1_ratio`.
            This solver is a good choice for `n_samples` >> `n_features`. Be aware
            that the memory usage of this solver has a quadratic dependency on
            `n_features` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.10

        'newton-cg'
            Uses a slightly adapted version of scipy's Newton-CG optimizer. This is
            sometimes called the truncated Newton method. Due to the fact that it
            does not construct the Hessian matrix but only uses gradients and
            vector products of the Hessian, it is a good solver when `X` is sparse
            or when `X` has many features.

            .. versionadded:: 1.10

        'newton-cholesky'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner Cholesky based solver.
            This solver is a good choice for `n_samples` >> `n_features`, 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` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.2

        .. 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):

           ================= ========================
           solver            l1_ratio
           ================= ========================
           'lbfgs'           l1_ratio=0
           'newton-cd-gram'  0<=l1_ratio<=1
           'newton-cg'       l1_ratio=0
           'newton-cholesky' l1_ratio=0
           ================= ========================

    max_iter : int, default=100
        The maximal number of iterations for the solver.
        Values must be in the range `[1, inf)`.

    tol : float, default=1e-4
        Stopping criterion. For the lbfgs solver,
        the iteration will stop when ``max{|g_j|, j = 1, ..., d} <= tol``
        where ``g_j`` is the j-th component of the gradient (derivative) of
        the objective function.
        Values must be in the range `(0.0, inf)`.

    warm_start : bool, default=False
        If set to ``True``, reuse the solution of the previous call to ``fit``
        as initialization for ``coef_`` and ``intercept_``.

    verbose : int, default=0
        For the lbfgs solver set verbose to any positive number for verbosity.
        Values must be in the range `[0, inf)`.

    Attributes
    ----------
    coef_ : array of shape (n_features,)
        Estimated coefficients for the linear predictor (`X @ coef_ +
        intercept_`) in the GLM.

    intercept_ : float
        Intercept (a.k.a. bias) added to linear predictor.

    n_iter_ : int
        Actual number of iterations used in the solver.

    _base_loss : BaseLoss, default=HalfSquaredError()
        This is set during fit via `self._get_loss()`.
        A `_base_loss` contains a specific loss function as well as the link
        function. The loss to be minimized specifies the distributional assumption of
        the GLM, i.e. the distribution from the EDM. Here are some examples:

        =======================  ========  ==========================
        _base_loss               Link      Target Domain
        =======================  ========  ==========================
        HalfSquaredError         identity  y any real number
        HalfPoissonLoss          log       0 <= y
        HalfGammaLoss            log       0 < y
        HalfTweedieLoss          log       dependent on tweedie power
        HalfTweedieLossIdentity  identity  dependent on tweedie power
        =======================  ========  ==========================

        The link function of the GLM, i.e. mapping from linear predictor
        `X @ coeff + intercept` to prediction `y_pred`. For instance, with a log link,
        we have `y_pred = exp(X @ coeff + intercept)`.
    """

    # We allow for NewtonSolver classes for the "solver" parameter but do not
    # make them public in the docstrings. This facilitates testing and
    # benchmarking.
    _parameter_constraints: dict = {
        "alpha": [Interval(Real, 0.0, None, closed="left")],
        "l1_ratio": [Interval(Real, 0, 1, closed="both")],
        "fit_intercept": ["boolean"],
        "solver": [
            StrOptions(
                {"lbfgs", "newton-cd", "newton-cd-gram", "newton-cg", "newton-cholesky"}
            ),
            Hidden(type),
        ],
        "max_iter": [Interval(Integral, 1, None, closed="left")],
        "tol": [Interval(Real, 0.0, None, closed="neither")],
        "warm_start": ["boolean"],
        "verbose": ["verbose"],
    }

    def __init__(
        self,
        *,
        alpha=1.0,
        l1_ratio=0.0,
        fit_intercept=True,
        solver="lbfgs",
        max_iter=100,
        tol=1e-4,
        warm_start=False,
        verbose=0,
    ):
        self.alpha = alpha
        self.l1_ratio = l1_ratio
        self.fit_intercept = fit_intercept
        self.solver = solver
        self.max_iter = max_iter
        self.tol = tol
        self.warm_start = warm_start
        self.verbose = verbose

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y, sample_weight=None):
        """Fit a Generalized Linear Model.

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

        y : array-like of shape (n_samples,)
            Target values.

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

        Returns
        -------
        self : object
            Fitted model.
        """
        if self.l1_ratio > 0 and self.solver not in ("newton-cd", "newton-cd-gram"):
            msg = (
                f"The solver '{self.solver}' does not support l1_ratio > 0; got "
                f"l1_ratio={self.l1_ratio}."
            )
            raise ValueError(msg)
        xp, _, device = get_namespace_and_device(X)
        X, y = validate_data(
            self,
            X,
            y,
            accept_sparse="csc" if self.solver == "newton-cd" else ["csc", "csr"],
            dtype=[xp.float64, xp.float32],
            order="F" if self.solver == "newton-cd" else None,
            y_numeric=True,
            multi_output=False,
        )
        y, sample_weight = move_to(y, sample_weight, xp=xp, device=device)

        loss_dtype = X.dtype
        y = check_array(y, dtype=loss_dtype, order="C", ensure_2d=False)

        if sample_weight is not None:
            # Note that _check_sample_weight calls check_array(order="C") required by
            # losses.
            sample_weight = _check_sample_weight(sample_weight, X, dtype=loss_dtype)

        n_samples, n_features = X.shape
        self._base_loss = self._get_loss(xp=xp, device=device)

        linear_loss = LinearModelLoss(
            base_loss=self._base_loss,
            fit_intercept=self.fit_intercept,
        )

        if not linear_loss.base_loss.in_y_true_range(y):
            raise ValueError(
                "Some value(s) of y are out of the valid range of the loss"
                f" {self._base_loss.__class__.__name__!r}."
            )

        # TODO: if alpha=0 check that X is not rank deficient

        # NOTE: Rescaling of sample_weight:
        # We want to minimize
        #     obj = 1/(2 * sum(sample_weight)) * sum(sample_weight * deviance)
        #         + a * L1 + 1/2 * b * L2
        # with
        #     deviance = 2 * loss, a = alpha * l1_ratio, b = alpha * (1 - l1_ratio).
        # The objective is invariant to multiplying sample_weight by a constant. We
        # could choose this constant such that sum(sample_weight) = 1 in order to end
        # up with
        #     obj = sum(sample_weight * loss) + a * L1 + 1/2 * b * L2.
        # But LinearModelLoss.loss() already computes
        #     average(loss, weights=sample_weight)
        # Thus, without rescaling, we have
        #     obj = LinearModelLoss.loss(...)

        loss_dtype_np = _matching_numpy_dtype(X, xp=xp)
        if self.warm_start and hasattr(self, "coef_"):
            coef_xp, _ = get_namespace(self.coef_)
            coef = move_to(self.coef_, xp=np, device="cpu")
            if self.fit_intercept:
                # LinearModelLoss needs intercept at the end of coefficient array.
                intercept = move_to(self.intercept_, xp=np, device="cpu")
                coef = np.concatenate((coef, np.array([intercept])))
            coef = coef.astype(loss_dtype_np, copy=False)
        else:
            coef = linear_loss.init_zero_coef(X, dtype=loss_dtype_np)
            if self.fit_intercept:
                coef[-1] = linear_loss.base_loss.link.link(
                    _average(y, weights=sample_weight)
                )

        l1_reg_strength = self.l1_ratio * self.alpha
        l2_reg_strength = (1 - self.l1_ratio) * self.alpha
        n_threads = _openmp_effective_n_threads()

        # Algorithms for optimization:
        # Note again that our losses implement 1/2 * deviance.
        if self.solver == "lbfgs":
            func = linear_loss.loss_gradient

            opt_res = scipy.optimize.minimize(
                func,
                coef,
                method="L-BFGS-B",
                jac=True,
                options={
                    "maxiter": self.max_iter,
                    "maxls": 50,  # default is 20
                    "gtol": self.tol,
                    # The constant 64 was found empirically to pass the test suite.
                    # The point is that ftol is very small, but a bit larger than
                    # machine precision for float64, which is the dtype used by lbfgs.
                    "ftol": 64 * np.finfo(float).eps,
                    **_get_additional_lbfgs_options_dict("iprint", self.verbose - 1),
                },
                args=(X, y, sample_weight, 0, l2_reg_strength, n_threads),
            )
            self.n_iter_ = _check_optimize_result(
                "lbfgs", opt_res, max_iter=self.max_iter
            )
            coef = opt_res.x
            coef = xp.asarray(
                coef.copy(order="C" if not _is_numpy_namespace(xp) else "K"),
                dtype=X.dtype,
                device=device,
            )
        elif self.solver == "newton-cg":
            func = linear_loss.loss
            grad = linear_loss.gradient
            hess = linear_loss.gradient_hessian_product  # hess = [gradient, hessp]

            coef, self.n_iter_ = _newton_cg(
                grad_hess=hess,
                func=func,
                grad=grad,
                x0=coef,
                args=(X, y, sample_weight, 0, l2_reg_strength, n_threads),
                maxiter=self.max_iter,
                tol=self.tol,
                verbose=self.verbose,
            )
        elif self.solver in ("newton-cd", "newton-cd-gram", "newton-cholesky"):
            if self.solver == "newton-cholesky":
                sol = NewtonCholeskySolver
                params = dict()
            elif self.solver == "newton-cd":
                sol = NewtonCDSolver
                params = dict(l1_reg_strength=l1_reg_strength)
            else:
                sol = NewtonCDGramSolver
                params = dict(l1_reg_strength=l1_reg_strength)
            sol = sol(
                coef=coef,
                linear_loss=linear_loss,
                **params,
                l2_reg_strength=l2_reg_strength,
                tol=self.tol,
                max_iter=self.max_iter,
                n_threads=n_threads,
                verbose=self.verbose,
            )
            coef = sol.solve(X, y, sample_weight)
            self.n_iter_ = sol.iteration
        elif issubclass(self.solver, NewtonSolver):
            sol = self.solver(
                coef=coef,
                linear_loss=linear_loss,
                l2_reg_strength=l2_reg_strength,
                tol=self.tol,
                max_iter=self.max_iter,
                n_threads=n_threads,
            )
            coef = sol.solve(X, y, sample_weight)
            self.n_iter_ = sol.iteration
        else:
            raise ValueError(f"Invalid solver={self.solver}.")

        if self.fit_intercept:
            self.intercept_ = coef[-1]
            self.coef_ = coef[:-1]
        else:
            # set intercept to zero as the other linear models do
            self.intercept_ = 0.0
            self.coef_ = coef

        return self

    def _linear_predictor(self, X):
        """Compute the linear_predictor = `X @ coef_ + intercept_`.

        Note that we often use the term raw_prediction instead of linear predictor.

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

        Returns
        -------
        y_pred : array of shape (n_samples,)
            Returns predicted values of linear predictor.
        """
        xp, _ = get_namespace(X)
        check_is_fitted(self)
        X = validate_data(
            self,
            X,
            accept_sparse=["csr", "csc", "coo"],
            dtype=[xp.float64, xp.float32],
            ensure_2d=True,
            allow_nd=False,
            reset=False,
        )
        return X @ self.coef_ + self.intercept_

    def predict(self, X):
        """Predict using GLM with feature matrix X.

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

        Returns
        -------
        y_pred : array of shape (n_samples,)
            Returns predicted values.
        """
        # check_array is done in _linear_predictor
        raw_prediction = self._linear_predictor(X)
        y_pred = self._base_loss.link.inverse(raw_prediction)
        return y_pred

    def score(self, X, y, sample_weight=None):
        """Compute D^2, the percentage of deviance explained.

        D^2 is a generalization of the coefficient of determination R^2.
        R^2 uses squared error and D^2 uses the deviance of this GLM, see the
        :ref:`User Guide <regression_metrics>`.

        D^2 is defined as
        :math:`D^2 = 1-\\frac{D(y_{true},y_{pred})}{D_{null}}`,
        :math:`D_{null}` is the null deviance, i.e. the deviance of a model
        with intercept alone, which corresponds to :math:`y_{pred} = \\bar{y}`.
        The mean :math:`\\bar{y}` is averaged by sample_weight.
        Best possible score is 1.0 and it can be negative (because the model
        can be arbitrarily worse).

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

        y : array-like of shape (n_samples,)
            True values of target.

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

        Returns
        -------
        score : float
            D^2 of self.predict(X) w.r.t. y.
        """
        # TODO: Adapt link to User Guide in the docstring, once
        # https://github.com/scikit-learn/scikit-learn/pull/22118 is merged.
        #
        # Note, default score defined in RegressorMixin is R^2 score.
        # TODO: make D^2 a score function in module metrics (and thereby get
        #       input validation and so on)
        raw_prediction = self._linear_predictor(X)  # validates X

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

        # required by losses
        y = check_array(y, dtype=raw_prediction.dtype, order="C", ensure_2d=False)

        if sample_weight is not None:
            # Note that _check_sample_weight calls check_array(order="C") required by
            # losses.
            sample_weight = _check_sample_weight(sample_weight, X, dtype=y.dtype)

        base_loss = self._base_loss

        if not base_loss.in_y_true_range(y):
            raise ValueError(
                "Some value(s) of y are out of the valid range of the loss"
                f" {base_loss.__name__}."
            )

        constant = _average(
            base_loss.constant_to_optimal_zero(y_true=y, sample_weight=None),
            weights=sample_weight,
        )

        # Missing factor of 2 in deviance cancels out.
        deviance = base_loss(
            y_true=y,
            raw_prediction=raw_prediction,
            sample_weight=sample_weight,
            n_threads=1,
        )
        y_mean = base_loss.link.link(_average(y, weights=sample_weight))
        deviance_null = base_loss(
            y_true=y,
            raw_prediction=xp.tile(y_mean, (y.shape[0],)),
            sample_weight=sample_weight,
            n_threads=1,
        )
        return float(1 - (deviance + constant) / (deviance_null + constant))

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        tags.input_tags.sparse = True
        try:
            # Create instance of BaseLoss if fit wasn't called yet. This is necessary as
            # TweedieRegressor might set the used loss during fit different from
            # self._base_loss.
            base_loss = self._get_loss()
            tags.target_tags.positive_only = not base_loss.in_y_true_range(-1.0)
        except (ValueError, AttributeError, TypeError):
            # This happens when the link or power parameter of TweedieRegressor is
            # invalid. We fallback on the default tags in that case.
            pass  # pragma: no cover
        return tags

    def _get_loss(self, xp=None, device=None):
        """This is only necessary because of the link and power arguments of the
        TweedieRegressor.

        Note that we do not need to pass sample_weight to the loss class as this is
        only needed to set loss.constant_hessian on which GLMs do not rely.
        """
        return HalfSquaredError()


class PoissonRegressor(_GeneralizedLinearRegressor):
    """Generalized Linear Model with a Poisson distribution.

    This regressor uses the 'log' link function.

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

    .. versionadded:: 0.23

    Parameters
    ----------
    alpha : float, default=1
        Constant that multiplies the penalty terms and determines the
        regularization strength. ``alpha = 0`` is equivalent to unpenalized
        GLMs. In this case, the design matrix `X` must have full column rank
        (no collinearities).
        Values of `alpha` must be in the range `[0.0, inf)`.

    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, do not work with some
           solvers. See the parameter `solver` below, to know the compatibility between
           the penalty and solver.

        .. versionadded:: 1.10

    fit_intercept : bool, default=True
        Specifies if a constant (a.k.a. bias or intercept) should be
        added to the linear predictor (`X @ coef + intercept`).

    solver : {'lbfgs', 'newton-cg', 'newton-cholesky'}, default='lbfgs'
        Algorithm to use in the optimization problem:

        'lbfgs'
            Calls scipy's L-BFGS-B optimizer.

        'newton-cd'
            Uses Newton-Raphson steps in an iterated reweighted least squares fashion:
            The normal equations are cast as a weighted least squares problem with
            elastic-net penalty. The inner solver then uses a coordinate descent based
            solver. This way the full Hessian is used but never explicitly constructed.
            It can solve for all values of `l1_ratio`.
            This solver is a good choice for `n_features` > `n_samples`.

            .. versionadded:: 1.10

        'newton-cd-gram'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner coordinate descent based
            solver that uses the full Hessian/Gram matrix. It can solve for all
            values of `l1_ratio`.
            This solver is a good choice for `n_samples` >> `n_features`. Be aware
            that the memory usage of this solver has a quadratic dependency on
            `n_features` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.10

        'newton-cg'
            Uses a slightly adapted version of scipy's Newton-CG optimizer. This is
            sometimes called the truncated Newton method. Due to the fact that it
            does not construct the Hessian matrix but only uses gradients and
            vector products of the Hessian, it is a good solver when `X` is sparse
            or when `X` has many features.

            .. versionadded:: 1.10

        'newton-cholesky'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner Cholesky based solver.
            This solver is a good choice for `n_samples` >> `n_features`, 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` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.2

        .. 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):

           ================= ========================
           solver            l1_ratio
           ================= ========================
           'lbfgs'           l1_ratio=0
           'newton-cd-gram'  0<=l1_ratio<=1
           'newton-cg'       l1_ratio=0
           'newton-cholesky' l1_ratio=0
           ================= ========================

    max_iter : int, default=100
        The maximal number of iterations for the solver.
        Values must be in the range `[1, inf)`.

    tol : float, default=1e-4
        Stopping criterion. For the lbfgs solver,
        the iteration will stop when ``max{|g_j|, j = 1, ..., d} <= tol``
        where ``g_j`` is the j-th component of the gradient (derivative) of
        the objective function.
        Values must be in the range `(0.0, inf)`.

    warm_start : bool, default=False
        If set to ``True``, reuse the solution of the previous call to ``fit``
        as initialization for ``coef_`` and ``intercept_`` .

    verbose : int, default=0
        For the lbfgs solver set verbose to any positive number for verbosity.
        Values must be in the range `[0, inf)`.

    Attributes
    ----------
    coef_ : array of shape (n_features,)
        Estimated coefficients for the linear predictor (`X @ coef_ +
        intercept_`) in the GLM.

    intercept_ : float
        Intercept (a.k.a. bias) added to linear predictor.

    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_ : int
        Actual number of iterations used in the solver.

    See Also
    --------
    TweedieRegressor : Generalized Linear Model with a Tweedie distribution.

    Examples
    --------
    >>> from sklearn import linear_model
    >>> clf = linear_model.PoissonRegressor()
    >>> X = [[1, 2], [2, 3], [3, 4], [4, 3]]
    >>> y = [12, 17, 22, 21]
    >>> clf.fit(X, y)
    PoissonRegressor()
    >>> clf.score(X, y)
    np.float64(0.990)
    >>> clf.coef_
    array([0.121, 0.158])
    >>> clf.intercept_
    np.float64(2.088)
    >>> clf.predict([[1, 1], [3, 4]])
    array([10.676, 21.875])
    """

    _parameter_constraints: dict = {
        **_GeneralizedLinearRegressor._parameter_constraints
    }

    def __init__(
        self,
        *,
        alpha=1.0,
        l1_ratio=0.0,
        fit_intercept=True,
        solver="lbfgs",
        max_iter=100,
        tol=1e-4,
        warm_start=False,
        verbose=0,
    ):
        super().__init__(
            alpha=alpha,
            l1_ratio=l1_ratio,
            fit_intercept=fit_intercept,
            solver=solver,
            max_iter=max_iter,
            tol=tol,
            warm_start=warm_start,
            verbose=verbose,
        )

    def _get_loss(self, xp=None, device=None):
        if xp is None or _is_numpy_namespace(xp):
            return HalfPoissonLoss()
        else:
            return HalfPoissonLossArrayAPI(xp=xp, device=device)

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        tags.array_api_support = self.solver == "lbfgs"
        return tags


class GammaRegressor(_GeneralizedLinearRegressor):
    """Generalized Linear Model with a Gamma distribution.

    This regressor uses the 'log' link function.

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

    .. versionadded:: 0.23

    Parameters
    ----------
    alpha : float, default=1
        Constant that multiplies the penalty terms and determines the
        regularization strength. ``alpha = 0`` is equivalent to unpenalized
        GLMs. In this case, the design matrix `X` must have full column rank
        (no collinearities).
        Values of `alpha` must be in the range `[0.0, inf)`.

    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, do not work with some
           solvers. See the parameter `solver` below, to know the compatibility between
           the penalty and solver.

        .. versionadded:: 1.10

    fit_intercept : bool, default=True
        Specifies if a constant (a.k.a. bias or intercept) should be
        added to the linear predictor (`X @ coef + intercept`).

    solver : {'lbfgs', 'newton-cg', 'newton-cholesky'}, default='lbfgs'
        Algorithm to use in the optimization problem:

        'lbfgs'
            Calls scipy's L-BFGS-B optimizer.

        'newton-cd'
            Uses Newton-Raphson steps in an iterated reweighted least squares fashion:
            The normal equations are cast as a weighted least squares problem with
            elastic-net penalty. The inner solver then uses a coordinate descent based
            solver. This way the full Hessian is used but never explicitly constructed.
            It can solve for all values of `l1_ratio`.
            This solver is a good choice for `n_features` > `n_samples`.

            .. versionadded:: 1.10

        'newton-cd-gram'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner coordinate descent based
            solver that uses the full Hessian/Gram matrix. It can solve for all
            values of `l1_ratio`.
            This solver is a good choice for `n_samples` >> `n_features`. Be aware
            that the memory usage of this solver has a quadratic dependency on
            `n_features` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.10

        'newton-cg'
            Uses a slightly adapted version of scipy's Newton-CG optimizer. This is
            sometimes called the truncated Newton method. Due to the fact that it
            does not construct the Hessian matrix but only uses gradients and
            vector products of the Hessian, it is a good solver when `X` is sparse
            or when `X` has many features.

            .. versionadded:: 1.10

        'newton-cholesky'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner Cholesky based solver.
            This solver is a good choice for `n_samples` >> `n_features`, 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` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.2

        .. 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):

           ================= ========================
           solver            l1_ratio
           ================= ========================
           'lbfgs'           l1_ratio=0
           'newton-cd-gram'  0<=l1_ratio<=1
           'newton-cg'       l1_ratio=0
           'newton-cholesky' l1_ratio=0
           ================= ========================

    max_iter : int, default=100
        The maximal number of iterations for the solver.
        Values must be in the range `[1, inf)`.

    tol : float, default=1e-4
        Stopping criterion. For the lbfgs solver,
        the iteration will stop when ``max{|g_j|, j = 1, ..., d} <= tol``
        where ``g_j`` is the j-th component of the gradient (derivative) of
        the objective function.
        Values must be in the range `(0.0, inf)`.

    warm_start : bool, default=False
        If set to ``True``, reuse the solution of the previous call to ``fit``
        as initialization for `coef_` and `intercept_`.

    verbose : int, default=0
        For the lbfgs solver set verbose to any positive number for verbosity.
        Values must be in the range `[0, inf)`.

    Attributes
    ----------
    coef_ : array of shape (n_features,)
        Estimated coefficients for the linear predictor (`X @ coef_ +
        intercept_`) in the GLM.

    intercept_ : float
        Intercept (a.k.a. bias) added to linear predictor.

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

        .. versionadded:: 0.24

    n_iter_ : int
        Actual number of iterations used in the solver.

    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
    --------
    PoissonRegressor : Generalized Linear Model with a Poisson distribution.
    TweedieRegressor : Generalized Linear Model with a Tweedie distribution.

    Examples
    --------
    >>> from sklearn import linear_model
    >>> clf = linear_model.GammaRegressor()
    >>> X = [[1, 2], [2, 3], [3, 4], [4, 3]]
    >>> y = [19, 26, 33, 30]
    >>> clf.fit(X, y)
    GammaRegressor()
    >>> clf.score(X, y)
    np.float64(0.773)
    >>> clf.coef_
    array([0.073, 0.067])
    >>> clf.intercept_
    np.float64(2.896)
    >>> clf.predict([[1, 0], [2, 8]])
    array([19.483, 35.795])
    """

    _parameter_constraints: dict = {
        **_GeneralizedLinearRegressor._parameter_constraints
    }

    def __init__(
        self,
        *,
        alpha=1.0,
        l1_ratio=0.0,
        fit_intercept=True,
        solver="lbfgs",
        max_iter=100,
        tol=1e-4,
        warm_start=False,
        verbose=0,
    ):
        super().__init__(
            alpha=alpha,
            l1_ratio=l1_ratio,
            fit_intercept=fit_intercept,
            solver=solver,
            max_iter=max_iter,
            tol=tol,
            warm_start=warm_start,
            verbose=verbose,
        )

    def _get_loss(self, xp=None, device=None):
        return HalfGammaLoss()


class TweedieRegressor(_GeneralizedLinearRegressor):
    """Generalized Linear Model with a Tweedie distribution.

    This estimator can be used to model different GLMs depending on the
    ``power`` parameter, which determines the underlying distribution.

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

    .. versionadded:: 0.23

    Parameters
    ----------
    power : float, default=0
            The power determines the underlying target distribution according
            to the following table:

            +-------+------------------------+
            | Power | Distribution           |
            +=======+========================+
            | 0     | Normal                 |
            +-------+------------------------+
            | 1     | Poisson                |
            +-------+------------------------+
            | (1,2) | Compound Poisson Gamma |
            +-------+------------------------+
            | 2     | Gamma                  |
            +-------+------------------------+
            | 3     | Inverse Gaussian       |
            +-------+------------------------+

            For ``0 < power < 1``, no distribution exists.

    alpha : float, default=1
        Constant that multiplies the penalty terms and determines the
        regularization strength. ``alpha = 0`` is equivalent to unpenalized
        GLMs. In this case, the design matrix `X` must have full column rank
        (no collinearities).
        Values of `alpha` must be in the range `[0.0, inf)`.

    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, do not work with some
           solvers. See the parameter `solver` below, to know the compatibility between
           the penalty and solver.

        .. versionadded:: 1.10

    fit_intercept : bool, default=True
        Specifies if a constant (a.k.a. bias or intercept) should be
        added to the linear predictor (`X @ coef + intercept`).

    link : {'auto', 'identity', 'log'}, default='auto'
        The link function of the GLM, i.e. mapping from linear predictor
        `X @ coeff + intercept` to prediction `y_pred`. Option 'auto' sets
        the link depending on the chosen `power` parameter as follows:

        - 'identity' for ``power <= 0``, e.g. for the Normal distribution
        - 'log' for ``power > 0``, e.g. for Poisson, Gamma and Inverse Gaussian
          distributions

    solver : {'lbfgs', 'newton-cg', 'newton-cholesky'}, default='lbfgs'
        Algorithm to use in the optimization problem:

        'lbfgs'
            Calls scipy's L-BFGS-B optimizer.

        'newton-cd'
            Uses Newton-Raphson steps in an iterated reweighted least squares fashion:
            The normal equations are cast as a weighted least squares problem with
            elastic-net penalty. The inner solver then uses a coordinate descent based
            solver. This way the full Hessian is used but never explicitly constructed.
            It can solve for all values of `l1_ratio`.
            This solver is a good choice for `n_features` > `n_samples`.

            .. versionadded:: 1.10

        'newton-cd-gram'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner coordinate descent based
            solver that uses the full Hessian/Gram matrix. It can solve for all
            values of `l1_ratio`.
            This solver is a good choice for `n_samples` >> `n_features`. Be aware
            that the memory usage of this solver has a quadratic dependency on
            `n_features` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.10

        'newton-cg'
            Uses a slightly adapted version of scipy's Newton-CG optimizer. This is
            sometimes called the truncated Newton method. Due to the fact that it
            does not construct the Hessian matrix but only uses gradients and
            vector products of the Hessian, it is a good solver when `X` is sparse
            or when `X` has many features.

            .. versionadded:: 1.10

        'newton-cholesky'
            Uses Newton-Raphson steps (in arbitrary precision arithmetic equivalent to
            iterated reweighted least squares) with an inner Cholesky based solver.
            This solver is a good choice for `n_samples` >> `n_features`, 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` because it explicitly computes the Hessian matrix.

            .. versionadded:: 1.2

        .. 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):

           ================= ========================
           solver            l1_ratio
           ================= ========================
           'lbfgs'           l1_ratio=0
           'newton-cd-gram'  0<=l1_ratio<=1
           'newton-cg'       l1_ratio=0
           'newton-cholesky' l1_ratio=0
           ================= ========================

    max_iter : int, default=100
        The maximal number of iterations for the solver.
        Values must be in the range `[1, inf)`.

    tol : float, default=1e-4
        Stopping criterion. For the lbfgs solver,
        the iteration will stop when ``max{|g_j|, j = 1, ..., d} <= tol``
        where ``g_j`` is the j-th component of the gradient (derivative) of
        the objective function.
        Values must be in the range `(0.0, inf)`.

    warm_start : bool, default=False
        If set to ``True``, reuse the solution of the previous call to ``fit``
        as initialization for ``coef_`` and ``intercept_`` .

    verbose : int, default=0
        For the lbfgs solver set verbose to any positive number for verbosity.
        Values must be in the range `[0, inf)`.

    Attributes
    ----------
    coef_ : array of shape (n_features,)
        Estimated coefficients for the linear predictor (`X @ coef_ +
        intercept_`) in the GLM.

    intercept_ : float
        Intercept (a.k.a. bias) added to linear predictor.

    n_iter_ : int
        Actual number of iterations used in the solver.

    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
    --------
    PoissonRegressor : Generalized Linear Model with a Poisson distribution.
    GammaRegressor : Generalized Linear Model with a Gamma distribution.

    Examples
    --------
    >>> from sklearn import linear_model
    >>> clf = linear_model.TweedieRegressor()
    >>> X = [[1, 2], [2, 3], [3, 4], [4, 3]]
    >>> y = [2, 3.5, 5, 5.5]
    >>> clf.fit(X, y)
    TweedieRegressor()
    >>> clf.score(X, y)
    np.float64(0.839)
    >>> clf.coef_
    array([0.599, 0.299])
    >>> clf.intercept_
    np.float64(1.600)
    >>> clf.predict([[1, 1], [3, 4]])
    array([2.500, 4.599])
    """

    _parameter_constraints: dict = {
        **_GeneralizedLinearRegressor._parameter_constraints,
        "power": [Interval(Real, None, None, closed="neither")],
        "link": [StrOptions({"auto", "identity", "log"})],
    }

    def __init__(
        self,
        *,
        power=0.0,
        alpha=1.0,
        l1_ratio=0.0,
        fit_intercept=True,
        link="auto",
        solver="lbfgs",
        max_iter=100,
        tol=1e-4,
        warm_start=False,
        verbose=0,
    ):
        super().__init__(
            alpha=alpha,
            l1_ratio=l1_ratio,
            fit_intercept=fit_intercept,
            solver=solver,
            max_iter=max_iter,
            tol=tol,
            warm_start=warm_start,
            verbose=verbose,
        )
        self.link = link
        self.power = power

    def _get_loss(self, xp=None, device=None):
        if self.link == "auto":
            if self.power <= 0:
                # identity link
                return HalfTweedieLossIdentity(power=self.power)
            else:
                # log link
                return HalfTweedieLoss(power=self.power)

        if self.link == "log":
            return HalfTweedieLoss(power=self.power)

        if self.link == "identity":
            return HalfTweedieLossIdentity(power=self.power)