sklearn/_loss/link.py - main - githubmirror/scikit-learn

sklearn/_loss/link.py
268 строк7 KB
Christian Lorentzen
ENH make `_loss.link.py` array API compatible (#33345)
05 мар 2026, 13:33
05 мар 2026, 13:33ba3ea9b
"""
Module contains classes for invertible (and differentiable) link functions.
"""

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

from abc import ABC, abstractmethod
from dataclasses import dataclass
from math import ulp

from scipy.stats import gmean

from sklearn.utils._array_api import _expit, _logit, get_namespace
from sklearn.utils.extmath import softmax


@dataclass
class Interval:
    low: float
    high: float
    low_inclusive: bool
    high_inclusive: bool

    def __post_init__(self):
        """Check that low <= high"""
        if self.low > self.high:
            raise ValueError(
                f"One must have low <= high; got low={self.low}, high={self.high}."
            )

    def includes(self, x):
        """Test whether all values of x are in interval range.

        Parameters
        ----------
        x : ndarray
            Array whose elements are tested to be in interval range.

        Returns
        -------
        result : bool
        """
        xp, _ = get_namespace(x)
        if self.low_inclusive:
            low = xp.greater_equal(x, self.low)
        else:
            low = xp.greater(x, self.low)

        if not xp.all(low):
            return False

        if self.high_inclusive:
            high = xp.less_equal(x, self.high)
        else:
            high = xp.less(x, self.high)

        # Note: np.all returns numpy.bool_
        return bool(xp.all(high))


def _inclusive_low_high(interval):
    """Generate values low and high to be within the interval range.

    This is used in tests only.

    Returns
    -------
    low, high : tuple of floats
        The returned values low and high lie within the interval.
    """
    eps = 10 * ulp(1)
    if interval.low == -float("inf"):
        low = -1e10
    elif interval.low < 0:
        low = interval.low * (1 - eps) + eps
    else:
        low = interval.low * (1 + eps) + eps

    if interval.high == float("inf"):
        high = 1e10
    elif interval.high < 0:
        high = interval.high * (1 + eps) - eps
    else:
        high = interval.high * (1 - eps) - eps

    return float(low), float(high)


class BaseLink(ABC):
    """Abstract base class for differentiable, invertible link functions.

    Convention:
        - link function g: raw_prediction = g(y_pred)
        - inverse link h: y_pred = h(raw_prediction)

    For (generalized) linear models, `raw_prediction = X @ coef` is the so
    called linear predictor, and `y_pred = h(raw_prediction)` is the predicted
    conditional (on X) expected value of the target `y_true`.

    The methods are not implemented as staticmethods in case a link function needs
    parameters.
    """

    is_multiclass = False  # used for testing only

    # Usually, raw_prediction may be any real number and y_pred is an open
    # interval.
    # interval_raw_prediction = Interval(-float("inf"), float("inf"), False, False)
    interval_y_pred = Interval(-float("inf"), float("inf"), False, False)

    @abstractmethod
    def link(self, y_pred):
        """Compute the link function g(y_pred).

        The link function maps (predicted) target values to raw predictions,
        i.e. `g(y_pred) = raw_prediction`.

        Parameters
        ----------
        y_pred : array
            Predicted target values.

        Returns
        -------
        array
            Output array, element-wise link function.
        """

    @abstractmethod
    def inverse(self, raw_prediction):
        """Compute the inverse link function h(raw_prediction).

        The inverse link function maps raw predictions to predicted target
        values, i.e. `h(raw_prediction) = y_pred`.

        Parameters
        ----------
        raw_prediction : array
            Raw prediction values (in link space).

        Returns
        -------
        array
            Output array, element-wise inverse link function.
        """


class IdentityLink(BaseLink):
    """The identity link function g(x)=x."""

    def link(self, y_pred):
        return y_pred  # TODO: Should we copy?

    inverse = link


class LogLink(BaseLink):
    """The log link function g(x)=log(x)."""

    interval_y_pred = Interval(0, float("inf"), False, False)

    def link(self, y_pred):
        xp, _ = get_namespace(y_pred)
        return xp.log(y_pred)

    def inverse(self, raw_prediction):
        xp, _ = get_namespace(raw_prediction)
        return xp.exp(raw_prediction)


class LogitLink(BaseLink):
    """The logit link function g(x)=logit(x)."""

    interval_y_pred = Interval(0, 1, False, False)

    def link(self, y_pred):
        return _logit(y_pred)

    def inverse(self, raw_prediction):
        return _expit(raw_prediction)


class HalfLogitLink(BaseLink):
    """Half the logit link function g(x)=1/2 * logit(x).

    Used for the exponential loss.
    """

    interval_y_pred = Interval(0, 1, False, False)

    def link(self, y_pred):
        return 0.5 * _logit(y_pred)

    def inverse(self, raw_prediction):
        return _expit(2 * raw_prediction)


class MultinomialLogit(BaseLink):
    """The symmetric multinomial logit function.

    Convention:
        - y_pred.shape = raw_prediction.shape = (n_samples, n_classes)

    Notes:
        - The inverse link h is the softmax function.
        - The sum is over the second axis, i.e. axis=1 (n_classes).

    We have to choose additional constraints in order to make

        y_pred[k] = exp(raw_pred[k]) / sum(exp(raw_pred[k]), k=0..n_classes-1)

    for n_classes classes identifiable and invertible.
    We choose the symmetric side constraint where the geometric mean response
    is set as reference category, see [2]:

    The symmetric multinomial logit link function for a single data point is
    then defined as

        raw_prediction[k] = g(y_pred[k]) = log(y_pred[k]/gmean(y_pred))
        = log(y_pred[k]) - mean(log(y_pred)).

    Note that this is equivalent to the definition in [1] and implies mean
    centered raw predictions:

        sum(raw_prediction[k], k=0..n_classes-1) = 0.

    For linear models with raw_prediction = X @ coef, this corresponds to
    sum(coef[k], k=0..n_classes-1) = 0, i.e. the sum over classes for every
    feature is zero.

    Reference
    ---------
    .. [1] Friedman, Jerome; Hastie, Trevor; Tibshirani, Robert. "Additive
        logistic regression: a statistical view of boosting" Ann. Statist.
        28 (2000), no. 2, 337--407. doi:10.1214/aos/1016218223.
        https://projecteuclid.org/euclid.aos/1016218223

    .. [2] Zahid, Faisal Maqbool and Gerhard Tutz. "Ridge estimation for
        multinomial logit models with symmetric side constraints."
        Computational Statistics 28 (2013): 1017-1034.
        http://epub.ub.uni-muenchen.de/11001/1/tr067.pdf
    """

    is_multiclass = True
    interval_y_pred = Interval(0, 1, False, False)

    def symmetrize_raw_prediction(self, raw_prediction):
        xp, _ = get_namespace(raw_prediction)
        return raw_prediction - xp.mean(raw_prediction, axis=1)[:, None]

    def link(self, y_pred):
        xp, _ = get_namespace(y_pred)
        # geometric mean as reference category
        gm = gmean(y_pred, axis=1)
        return xp.log(y_pred / gm[:, None])

    def inverse(self, raw_prediction):
        return softmax(raw_prediction)


_LINKS = {
    "identity": IdentityLink,
    "log": LogLink,
    "logit": LogitLink,
    "half_logit": HalfLogitLink,
    "multinomial_logit": MultinomialLogit,
}