Другое

Исправление ошибок трансляции в пользовательских функциях потерь LightGBM

Узнайте, как исправлять ошибки трансляции при вычислении производных для пользовательских функций потерь отрицательного биномиального распределения в LightGBM. Изучите аналитические и численные решения с полными примерами кода.

Как вычислить производные для пользовательской функции потерь отрицательного биномиального распределения в LightGBM, когда функция производной из scipy вызывает ошибки трансляции (broadcasting)? Я реализую пользовательскую функцию потерь отрицательного биномиального распределения для моей модели LGBM со следующим кодом:

python
def negative_binomial_objective(y_true, y_pred):
    p = 0.5

    def loss(x, t):
        loss = (
            gamma(x) + gamma(t + 1) - gamma(x + t) - x * np.log(p) - t * np.log(1 - p)
        )
        return loss

    def compute_grad(x):
        return derivative(f=lambda x: loss(x, y_true), x=x)

    grad = compute_grad(x=y_pred)
    hess = derivative(f=compute_grad, x=y_pred)  # вторая производная

    return grad, hess

Проблема в том, что функция производной из scipy расширяет x для соответствия параметру order (по умолчанию установлен на 8). Это вызывает ошибку, что мои операнды не могут быть совместно трансформированы (broadcasted) вместе (например, внутри gamma(x + t)), потому что они имеют формы (N,8) и (N,). Как можно найти производную этой функции, не сталкиваясь с проблемами трансляции?

Содержание

Понимание проблемы трансляции

Проблема возникает потому, что scipy.differentiate.derivative по умолчанию использует схему конечных разностей с 8 точками, которая вычисляет вашу функцию в нескольких точках вокруг каждого входного значения. Это расширяет размерность вашего входного массива с (N,) до (N,8), вызывая ошибки трансляции при попытке выполнить операции с y_true, который имеет форму (N,).

Ваша текущая реализация:

python
def negative_binomial_objective(y_true, y_pred):
    p = 0.5

    def loss(x, t):
        loss = (
            gamma(x) + gamma(t + 1) - gamma(x + t) - x * np.log(p) - t * np.log(1 - p)
        )
        return loss

    def compute_grad(x):
        return derivative(f=lambda x: loss(x, y_true), x=x)

    grad = compute_grad(x=y_pred)
    hess = derivative(f=compute_grad, x=y_pred)  # вторая производная

    return grad, hess

Аналитическое решение для функции потера отрицательного биномиального распределения

Для распределения отрицательного биномиального мы можем вывести аналитические градиенты и гессианы, что более эффективно и полностью избегает проблем численного дифференцирования.

Функция правдоподобия отрицательного биномиального распределения:

L=lnΓ(x)+lnΓ(t+1)lnΓ(x+t)xln(p)tln(1p)L = \ln \Gamma(x) + \ln \Gamma(t + 1) - \ln \Gamma(x + t) - x \ln(p) - t \ln(1 - p)

Где:

  • xx - предсказанный параметр (параметр дисперсии)
  • tt - целевое значение
  • pp - параметр вероятности

Градиент (первая производная):

Lx=ψ(x)ψ(x+t)ln(p)\frac{\partial L}{\partial x} = \psi(x) - \psi(x + t) - \ln(p)

Гессиан (вторая производная):

2Lx2=ψ(x)ψ(x+t)\frac{\partial^2 L}{\partial x^2} = \psi'(x) - \psi'(x + t)

Здесь ψ\psi - дигамма-функция, а ψ\psi' - тригамма-функция.

Альтернатива численному дифференцированию

Если вы предпочитаете использовать численные производные, вот решение, которое избегает проблем трансляции:

python
import numpy as np
from scipy.special import gamma, digamma, trigamma
from scipy.differentiate import derivative

def negative_binomial_objective(y_true, y_pred):
    p = 0.5
    
    # Векторизованная функция потерь
    def loss_vec(x, t):
        return (gamma(x) + gamma(t + 1) - gamma(x + t) - 
                x * np.log(p) - t * np.log(1 - p))
    
    # Вычисление градиента с использованием конечных разностей (ручная реализация)
    def grad_vec(x, t, eps=1e-6):
        return (loss_vec(x + eps, t) - loss_vec(x - eps, t)) / (2 * eps)
    
    # Вычисление гессиана с использованием конечных разностей
    def hess_vec(x, t, eps=1e-6):
        return (loss_vec(x + eps, t) - 2 * loss_vec(x, t) + loss_vec(x - eps, t)) / (eps ** 2)
    
    # Векторизуем функции для корректной обработки массивов
    grad = np.vectorize(grad_vec)(y_pred, y_true)
    hess = np.vectorize(hess_vec)(y_pred, y_true)
    
    return grad, hess

Векторизованная реализация

Более эффективная векторизованная реализация с использованием аналитических производных:

python
import numpy as np
from scipy.special import digamma, trigamma

def negative_binomial_objective(y_true, y_pred):
    """
    Пользовательская функция потерь отрицательного биномиального распределения для LightGBM.
    
    Параметры:
    -----------
    y_true : array-like
        Истинные целевые значения
    y_pred : array-like  
        Предсказанные значения (параметр дисперсии)
        
    Возвращает:
    --------
    grad : array-like
        Значения градиента
    hess : array-like
        Значения гессиана
    """
    p = 0.5  # параметр вероятности
    eps = 1e-6  # малая константа для численной стабильности
    
    # Обеспечиваем положительные предсказания для гамма-функций
    y_pred = np.maximum(y_pred, eps)
    
    # Аналитический градиент с использованием дигамма-функции
    grad = digamma(y_pred) - digamma(y_pred + y_true) - np.log(p)
    
    # Аналитический гессиан с использованием тригамма-функции  
    hess = trigamma(y_pred) - trigamma(y_pred + y_true)
    
    # Обеспечиваем положительный гессиан для численной стабильности
    hess = np.maximum(hess, eps)
    
    return grad, hess

Полный пользовательский объект LightGBM

Вот полная реализация, которую можно использовать непосредственно с LightGBM:

python
import numpy as np
import lightgbm as lgb
from scipy.special import digamma, trigamma

def negative_binomial_objective(y_true, y_pred):
    """
    Пользовательская функция потерь отрицательного биномиального распределения для LightGBM.
    
    Эта функция реализует логарифм правдоподобия распределения отрицательного биномиального
    с аналитическими вычислениями градиента и гессиана.
    
    Параметры:
    -----------
    y_true : array-like
        Истинные целевые значения (счета)
    y_pred : array-like  
        Предсказанные значения (должны быть положительными для гамма-функций)
        
    Возвращает:
    --------
    grad : array-like
        Первая производная (градиент) функции потерь
    hess : array-like
        Вторая производная (гессиан) функции потерь
        
    Математическая формулировка:
    -----------------------
    Логарифм правдоподобия: L = ln Γ(x) + ln Γ(t + 1) - ln Γ(x + t) - x ln(p) - t ln(1 - p)
    
    Градиент: ∂L/∂x = ψ(x) - ψ(x + t) - ln(p)
    Гессиан: ∂²L/∂x² = ψ'(x) - ψ'(x + t)
    
    Где:
    - x: предсказанный параметр дисперсии
    - t: истинное целевое значение  
    - p: параметр вероятности (фиксирован на 0.5)
    - ψ: дигамма-функция
    - ψ': тригамма-функция
    """
    p = 0.5  # параметр вероятности
    eps = 1e-6  # константа для численной стабильности
    
    # Обеспечиваем положительные предсказания для гамма-функций
    y_pred = np.maximum(y_pred, eps)
    
    # Аналитический градиент с использованием дигамма-функции
    grad = digamma(y_pred) - digamma(y_pred + y_true) - np.log(p)
    
    # Аналитический гессиан с использованием тригамма-функции
    hess = trigamma(y_pred) - trigamma(y_pred + y_true)
    
    # Обеспечиваем положительный гессиан для численной стабильности
    hess = np.maximum(hess, eps)
    
    return grad, hess

def negative_binomial_metric(y_true, y_pred):
    """
    Пользовательская метрика оценки для модели отрицательного биномиального распределения.
    
    Параметры:
    -----------
    y_true : array-like
        Истинные целевые значения
    y_pred : array-like
        Предсказанные значения
        
    Возвращает:
    --------
    metric_name : str
        Название метрики
    metric_value : float
        Вычисленное значение метрики
    is_higher_better : bool
        Являются ли более высокие значения показателем лучшей производительности
    """
    p = 0.5
    eps = 1e-6
    
    y_pred = np.maximum(y_pred, eps)
    
    # Отрицательное логарифмическое правдоподобие в качестве метрики
    nll = -(gamma(y_pred) + gamma(y_true + 1) - gamma(y_pred + y_true) - 
            y_pred * np.log(p) - y_true * np.log(1 - p))
    
    # Возвращаем среднее отрицательное логарифмическое правдоподобие
    return 'neg_log_likelihood', np.mean(nll), False

# Пример использования
if __name__ == "__main__":
    # Генерируем тестовые данные
    np.random.seed(42)
    n_samples = 1000
    X = np.random.randn(n_samples, 5)
    
    # Генерируем целевые значения, распределенные по отрицательному биномиальному закону
    true_dispersion = np.random.gamma(2, 1, n_samples)
    y_true = np.random.negative_binomial(true_dispersion, 0.5)
    
    # Создаем набор данных LightGBM
    train_data = lgb.Dataset(X, label=y_true)
    
    # Настраиваем параметры
    params = {
        'objective': 'custom',
        'learning_rate': 0.1,
        'num_leaves': 31,
        'verbose': -1,
        'seed': 42
    }
    
    # Обучаем модель с пользовательской функцией потерь
    model = lgb.train(
        params,
        train_data,
        num_boost_round=100,
        fobj=negative_binomial_objective,
        feval=negative_binomial_metric,
        valid_sets=[train_data],
        valid_names=['train'],
        callbacks=[lgb.log_evaluation(10)]
    )

Тестирование и валидация

Чтобы убедиться, что ваша пользовательская функция потерь работает правильно, вот тестовый скрипт:

python
import numpy as np
from scipy.special import gamma, digamma, trigamma

def test_negative_binomial_objective():
    """Тест пользовательской функции потерь отрицательного биномиального распределения."""
    
    # Тестовые данные
    y_true = np.array([1, 2, 3, 4, 5])
    y_pred = np.array([1.5, 2.0, 2.5, 3.0, 3.5])
    
    # Ожидаемое вычисление градиента
    p = 0.5
    eps = 1e-6
    
    # Численный градиент для проверки
    def loss(x, t):
        return (gamma(x) + gamma(t + 1) - gamma(x + t) - 
                x * np.log(p) - t * np.log(1 - p))
    
    grad_numeric = np.array([(loss(x + eps, t) - loss(x - eps, t)) / (2 * eps) 
                            for x, t in zip(y_pred, y_true)])
    
    # Аналитический градиент
    grad_analytical = digamma(y_pred) - digamma(y_pred + y_true) - np.log(p)
    
    # Численный гессиан для проверки  
    hess_numeric = np.array([(loss(x + eps, t) - 2 * loss(x, t) + loss(x - eps, t)) / (eps ** 2) 
                            for x, t in zip(y_pred, y_true)])
    
    # Аналитический гессиан
    hess_analytical = trigamma(y_pred) - trigamma(y_pred + y_true)
    
    print("Сравнение градиентов:")
    print("Численный:", grad_numeric)
    print("Аналитический:", grad_analytical)
    print("Разница:", np.abs(grad_numeric - grad_analytical))
    print("Максимальная разница:", np.max(np.abs(grad_numeric - grad_analytical)))
    
    print("\nСравнение гессианов:")
    print("Численный:", hess_numeric)
    print("Аналитический:", hess_analytical)
    print("Разница:", np.abs(hess_numeric - hess_analytical))
    print("Максимальная разница:", np.max(np.abs(hess_numeric - hess_analytical)))
    
    # Проверка форм
    assert grad_analytical.shape == y_true.shape
    assert hess_analytical.shape == y_true.shape
    
    # Проверка, что гессиан положителен (для выпуклости)
    assert np.all(hess_analytical > 0), "Гессиан должен быть положительным"
    
    print("\nВсе тесты пройдены!")

if __name__ == "__main__":
    test_negative_binomial_objective()

Источники

  1. Документация по пользовательским объектам LightGBM
  2. Распределение отрицательного биномиального - Википедия
  3. Дигамма и тригамма функции - Документация SciPy
  4. Пользовательские функции потерь для LightGBM – Nima Sarang
  5. Пользовательский объект для LightGBM | Hippocampus’s Garden

Заключение

Чтобы решить проблему ошибки трансляции функции derivative scipy в вашем пользовательском объекте LightGBM:

  1. Используйте аналитические производные - Предпочтительный подход - вывести аналитические градиент и гессиан с использованием дигамма и тригамма функций, что более эффективно и точно.

  2. Реализуйте ручные конечные разности - Если необходимо использовать численные производные, реализуйте их вручную с правильной обработкой массивов для избежания проблем трансляции.

  3. Обеспечьте численную стабильность - Добавляйте малые константы (eps) для предотвращения деления на ноль и обеспечения положительных значений гессиана.

  4. Тщательно тестируйте - Проверяйте вашу реализацию путем сравнения численных и аналитических производных.

Аналитический подход с использованием дигамма и тригамма функций рекомендуется, так как он обеспечивает наилучшее сочетание точности, эффективности и численной стабильности для функций потерь отрицательного биномиального распределения в градиентном бустинге.

Авторы
Проверено модерацией
Модерация