LightGBMのカスタムobjectiveで組み込みを再現する

この記事は MicroAd Advent Calendar 2024 の21日目の記事です。

https://qiita.com/advent-calendar/2024/microad

新卒で機械学習エンジニアをしている Sakishita です。
今回は普段使っている LightGBM について、数式通りに実装したら数値が合わず色々と躓いたので仕様を深堀りした内容をご紹介します。

## 概要

LightGBMでは、問題に合わせてブースティングに用いられる目的関数をユーザーが設計することで、モデル性能の向上を図ることができます。
この記事ではその準備段階として、カスタムで定義した関数を目的関数として指定する方法を紹介し、組み込みの`"binary"`（二値分類）と`"multiclass"`（多クラス分類）が再現してカスタム目的関数が正しく動作していることを確認します。

## 分類問題の目的関数

機械学習モデルでは、正解ラベル $\boldsymbol{y}$ とモデル出力 $\boldsymbol{z}$ の間の"距離"を表す目的関数 $\mathcal{L}( \boldsymbol{z}, \boldsymbol{y})$ を最小化するような予測関数 $\boldsymbol{z}=f(X)$ を構築することが目標です。この目的関数が小さくなるようにモデル $f$ の内部パラメータを調整するため、$\boldsymbol{z}$ に対する微分などの情報を用いられます。
LightGBM は Gradient Boosted Decision Tree という名前でありながら実は Gradient Boosting ではなく Newton Boosting[^nb]という手法を使っている[^newton]ため、gradient（1階偏微分）と hessian（2階偏微分）を求める必要があります。
hessian（ヘッセ行列）とは通常第 $i$ 成分と第 $j$ 成分でそれぞれ偏微分した $N\times N$ の行列を指しますが、LightGBM ではサンプル間の独立性を仮定して **対角成分のみ** を計算に利用しています。
以下は、grad と hess の定義式です。

$$
\mathrm{grad}_i=\frac{\partial\mathcal{L}(\boldsymbol{z}, \boldsymbol{y})}{\partial z_i},\ \mathrm{hess}_i=\frac{\partial^2\mathcal{L}(\boldsymbol{z}, \boldsymbol{y})}{\partial z_i^2}
$$

（ $i$：サンプルのインデックス、$\mathcal{L}$：目的関数、$\boldsymbol{y}$：正解ラベル、$\boldsymbol{z}$：予測値）

[^nb]: 根を求める[ニュートン法](https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95)ではなく[最適化におけるニュートン法](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization)を決定木学習に応用したもの。詳しくは [Gradient Boosting と XGBoost｜Zenn](https://zenn.dev/zaburo_ch/articles/71b1af9bba1571) を参照。

[^newton]: XGBoost と 同様に -grad/hess について木を最適化しているとの記述がある https://github.com/microsoft/LightGBM/issues/5233 

### 二値分類

二値分類でよく使われる目的関数（Log loss、Cross Entropy loss などと呼ばれる）は、正解ラベル $\boldsymbol{y}$、予測確率 $\boldsymbol{p}$ に対して以下のように表されます。

$$
\mathcal{L}( \boldsymbol{p}, \boldsymbol{y}) = -\sum_{i=1}^N( y_i\log p_i+(1-y_i)\log(1-p_i))
$$

ただし、$0\leq p \leq 1$ の制約などの都合でこのままでは計算しづらいため

```math
p=\operatorname{sigmoid}(z)=\frac{1}{1+\exp(-z)}
```

という変数変換を用いて実数全体を動く変数 $z$ に対する関数として計算することが多いです。

$y_i$, $z_i$（$i=1, \dots, N$）の関数である $\mathcal{L}( \boldsymbol{z}, \boldsymbol{y})$ を $z_i$ で偏微分して grad、hess を求めると

```math
\begin{align*}
\mathrm{grad}_i=\frac{\partial\mathcal{L}}{\partial z_i}& = p_i-y_i=\frac{-\sigma_i}{1+e^{\sigma_i z_i}}\\
\mathrm{hess}_i=\frac{\partial^2\mathcal{L}}{\partial z_i^2}& = p_i(1-p_i) =|\mathrm{grad}_i|(1-|\mathrm{grad}_i|)
\end{align*}
```

（ $\sigma\equiv 2 y-1$：$y=0/1$ を $\sigma=-/+$ に変換したもの）

となります。
$z$ を $p$ に変換した表現と、$\sigma$ や絶対値を使った表現の2通りありますが同じものです。どちらにせよ $e^x$ の計算をしなければいけないので計算量や誤差は大差ないはずですが、LightGBM の組み込み実装では後者を使っています。

Cross Entropy loss の目的関数の意味づけや勾配の詳しい導出を知りたい方は記事末尾の補足ノートをご参照ください。

### 多クラス分類

二値分類で $y=0\ \text{or}\ 1$ だったのを $y=0, \dots, K-1$ のように多値に拡張した場合、Cross Entropy loss 関数は $y_i=k$ のとき $t_{ik}=1$、それ以外では $t_{ik}=0$ の one-hot エンコーディングした変数を使って、まとめて

$$
\mathcal{L}(\{p_{ik}\}, \{t_{ik}\})=-\sum_i\sum_k\{t_{ik}\log(p_{ik})\}
$$

と表せます。
ここでも $0\leq p_{ik}\leq 1$ かつ $\sum_k p_{ik}=1$ という制約条件を扱いやすくするために

```math
p_{ik} = \operatorname{softmax}(z_{ik})=\frac{e^{z_{ik}}}{\sum_{k'=0}^{K-1}e^{z_{ik'}}}
```

の変数変換を導入します（$K=2$、$k=1$ の場合を考えると二値のときの sigmoid 関数と一致します）。
grad と hess をサンプルのインデックス $i$ とクラスのインデックス $k$ について各々計算して

```math
\begin{align*}
\mathrm{grad}_{ik}=\frac{\partial \mathcal{L}}{\partial z_{ik}}&=-t_{ik}+p_{ik}\\
\mathrm{hess}_{ik}=\frac{\partial^2 \mathcal{L}}{\partial z_{ik}^2}&=p_{ik}(1-p_{ik})
\end{align*}
```

が得られます。
ただし、LightGBM の実装では **hess が $K/(K-1)$ 倍されています** [^hess]。これはクラス数に対して確率の総和=1の制約で自由度が1減っていること・サンプル間の相関を無視していることの補正や、hessian が小さすぎて収束が遅くなることの対策らしいです[^hesss]。

[^hess]: [mutliclass_objective の実装](https://github.com/microsoft/LightGBM/blob/master/src/objective/multiclass_objective.hpp#L105)

[^hesss]: [Gradient Boost Modelにおける多クラス分類の目的関数について#Hessianの乗数について](https://zenn.dev/bilzard/articles/obj-func-of-multi-class-on-gbm)

## 実装

LightGBM で目的関数を設定し、その挙動を確認します。
sklearn API では、`objective` パラメータに`"binary"`, `"multiclass"`, ...などの組み込みで用意されている目的関数を指す文字列か、以下の形式の関数をカスタム目的関数として渡します。
- `func(y_true: np.ndarray, y_pred: np.ndarray) -> (grad: np.ndarray, hess: np.ndarray)`

`y_pred` の形状は正解ラベルから自動で決まるようになっており、0,1だけの場合は `shape=(n_samples,)`、0から$K-1$までの整数である場合は `shape=(n_samples, n_classes)` になります。
`grad`, `hess` はそれぞれ `y_pred` と同じ形状のndarrayです。
カスタム関数を渡した場合は、$\boldsymbol{z}$ から $\boldsymbol{p}$ への変数変換が自動で行われなくなるので、 **後処理を自分でやる必要がある** ことに注意が必要です。

train API を使用する場合は、`train` 関数の `params` パラメータに辞書として `{"objective": func}` のように渡します。
- `func(pred: np.ndarray, train_data: lgb.Dataset) -> (grad, hess)`

関数の引数形式が微妙に異なるため注意。正解ラベルは `train_data.labels` に格納されているので関数内で取り出す必要があります。
正解ラベル以外のデータを参照するなど、複雑な目的関数を設計する場合はこちらが便利です。

以下で`"binary"`（二値分類）と`"multiclass"`（多クラス分類）をカスタム関数で再現できることを確かめます。

:::note warn
LightGBM 組み込みの目的関数を使用する場合はデフォルトで `boost_from_average=True` となっており「初期スコアをラベルの平均に調整する」機能が適用されていますが、カスタム関数ではこの機能が適用されません[^boost]。このため学習結果が組み込みの目的関数を使ったときと異なる場合があります。
:::

`boost_from_average` の実装をまだ追い切れていないため、本記事では `boost_from_average=False` として比較を行います。（詳しくは [Focal loss implementation for LightGBM • Max Halford](https://maxhalford.github.io/blog/lightgbm-focal-loss/) を参照）

[^boost]: "used only in `regression`, `binary`, `multiclassova` and `cross-entropy` applications" https://lightgbm.readthedocs.io/en/latest/Parameters.html#boost_from_average

動作環境は

```
Python 3.12.4
lightgbm                  4.5.0
numpy                     2.0.0
scikit-learn              1.5.1
```

です。LightGBM v3系では API が異なるため注意してください。

### `"binary"`の再現

#### サンプルデータ生成

ベンチマーク用に適当な二値分類問題を生成します。

```python
import lightgbm as lgb
import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn import metrics


# 適当な問題を生成
X, y = make_classification(n_samples=1000, n_features=10, n_informative=5, n_redundant=3, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
```

#### 組み込みの`"binary"`

```python
lgb_clf = lgb.LGBMClassifier(
    n_estimators=100,
    objective='binary',  # 'cross_entropy' でも同じ
    random_state=42,
    boost_from_average=False,
    verbose=-1,
)
lgb_clf.fit(X_train, y_train)
metrics.roc_auc_score(y_test, lgb_clf.predict_proba(X_test)[:, 1])
# => 0.9872987298729874
```

`boost_from_average=False` を入れないとカスタム関数と一致しなくなる。

#### 二値 Cross Entropy loss のカスタム関数

```python
def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def binary_cross_entropy_loss(y_true: np.ndarray, y_pred: np.ndarray):
    label = (y_true > 0) * 2 - 1  # pos -> 1, neg -> -1
    response = -label / (1.0 + np.exp(label * y_pred))
    abs_response = np.abs(response)
    grad = response
    hess = abs_response * (1.0 - abs_response)
    # pred = sigmoid(y_pred)
    # grad = pred - y_true
    # hess = pred * (1 - pred)
    return grad, hess

custom_clf = lgb.LGBMModel(  # lgb.LGBMClassifier でも良い（警告抑制のため）
    n_estimators=100,
    objective=binary_cross_entropy_loss,
    random_state=42,
    # "boost_from_average": False,
    verbose=-1,
)
custom_clf.fit(X_train, y_train)
custom_pred = sigmoid(custom_clf.predict(X_test))  # 後処理
metrics.roc_auc_score(y_test, custom_pred)
# => 0.9872987298729874
assert (lgb_clf.predict_proba(X_test)[:, 1] == sigmoid(custom_clf.predict(X_test))).all()
# OK
```

予測データが一致していることが確認できました。

Cross Entropy loss の実装はコメントアウトしてあるほうでも同様の結果が得られます[^fast]。
また、予測モデルのクラスは `LGBMClassifier` でも大丈夫ですが、出力の変数変換がされていないことの注意喚起のためwarningが出るので `LGBMModel` の方がベターです。

[^fast]: Pythonの実装では若干こちらのほうが早い

<details><summary>warningの内容</summary>

```
>>> custom_clf.predict(X_test)
UserWarning: Cannot compute class probabilities or labels due to the usage of customized objective function.
Returning raw scores instead.
```

</details>


### `"multiclass"` の再現

#### サンプルデータ生成

適当ですが、広告データセットを想定して「Clickもコンバージョンもなし」=0、「Clickのみでコンバージョンはなし」=1、「Clickもコンバージョンもあり」=2という3値分類問題を生成しました。

```python
from sklearn.datasets import make_multilabel_classification
import pandas as pd

X, y = make_multilabel_classification(
    n_samples=1000,
    n_features=10,
    n_classes=2,
    n_labels=2,
    random_state=0,
)
y_df = pd.DataFrame(y, columns=["click", "click_to_cv"])
y_df["cv"] = y_df["click"] * y_df["click_to_cv"]
y_df["click_cv"] = y_df["click"] + y_df["cv"]
X_train, X_test, y_train, y_test = train_test_split(X, y_df, test_size=0.2, random_state=42)
```

#### 組み込みの`"multiclass"`

```python
multiclass_clf = lgb.LGBMClassifier(
    n_estimators=100,
    objective='multiclass',
    n_classes=3,
    random_state=42,
    boost_from_average=False,
    verbose=-1,
)
multiclass_clf.fit(X_train, y_train["click_cv"])
multiclass_pred = multiclass_clf.predict_proba(X_test)
print(
    "multiclass:",
    metrics.roc_auc_score(y_test["click"], multiclass_pred[:, 1:].sum(axis=1)),
    metrics.roc_auc_score(y_test["cv"], multiclass_pred[:, 2]),
)
# multiclass: 0.9569 0.9181325367850426
```

#### 多クラス Cross Entropy loss のカスタム関数

```python
def softmax(logit: np.ndarray):
    vmax = logit.max(axis=1)
    p: np.ndarray = np.exp(logit - vmax[:, np.newaxis])
    p /= p.sum(axis=1)[:, np.newaxis]
    return p

def multi_cross_entropy_loss(y_true: np.ndarray, y_pred: np.ndarray):
    num_class = y_pred.shape[1]
    label = np.zeros(y_pred.shape)
    for k in range(num_class):
        label[y_true == k, k] = 1
    pred = softmax(y_pred)
    grad = pred - label
    hess = pred * (1 - pred) * num_class / (num_class - 1)  # 補正項
    return grad, hess

custom_clf = lgb.LGBMModel(
    n_estimators=100,
    objective=multi_cross_entropy_loss,
    n_classes=3,
    random_state=42,
    verbose=-1,
)
custom_clf.fit(X_train, y_train["click_cv"])
custom_pred = softmax(custom_clf.predict(X_test))  # 後処理
print(
    "custom:",
    metrics.roc_auc_score(y_test["click"], custom_pred[:, 1:].sum(axis=1)),
    metrics.roc_auc_score(y_test["cv"], custom_pred[:, 2]),
)
# custom: 0.9569 0.9181325367850426
assert (multiclass_pred == custom_pred).all()
# OK
```

予測データが一致していることが確認できました。

### train API での実装

<details><summary>二値分類</summary>

```python:組み込みbinary
lgb_train = lgb.Dataset(X_train, y_train)
binary_booster = lgb.train(
    params={
        "objective": "binary",
        "random_state": 42,
        "boost_from_average": False,
        "verbose": -1,
    },
    train_set=lgb_train,
    num_boost_round=100,
)

print(
    """train binary:  """,
    metrics.roc_auc_score(y_test, binary_booster.predict(X_test)),
)
```

```python:カスタムbinary
def train_binary_cross_entropy_loss(preds: np.ndarray, train_data: lgb.Dataset):
    return binary_cross_entropy_loss(train_data.label, preds)


lgb_train = lgb.Dataset(X_train, y_train)
custom_booster = lgb.train(
    params={
        "objective": train_binary_cross_entropy_loss,
        "random_state": 42,
        # "boost_from_average": False,
        "verbose": -1,
    },
    train_set=lgb_train,
    num_boost_round=100,
)

print(
    """train custom:  """,
    metrics.roc_auc_score(y_test, sigmoid(custom_booster.predict(X_test))),
)

assert (binary_booster.predict(X_test) == sigmoid(custom_booster.predict(X_test))).all()
```

</details>


<details><summary>多クラス分類</summary>

```python:組み込みmulticlass
lgb_train = lgb.Dataset(X_train, y_train)
multiclass_booster = lgb.train(
    params={
        "objective": "multiclass",
        "num_class": 3,
        "random_state": 42,
        "boost_from_average": False,
        "verbose": -1,
    },
    train_set=lgb_train,
    num_boost_round=100,
)

multiclass_pred = multiclass_booster.predict(X_test)
```

```python:カスタムmulticlass
def train_multi_cross_entropy_loss(preds: np.ndarray, train_data: lgb.Dataset):
    return multi_cross_entropy_loss(train_data.label, preds)


lgb_train = lgb.Dataset(X_train, y_train)
custom_booster = lgb.train(
    params={
        "objective": train_multi_cross_entropy_loss,
        "num_class": 3,
        "random_state": 42,
        # "boost_from_average": False,
        "verbose": -1,
    },
    train_set=lgb_train,
    num_boost_round=100,
)
custom_pred = softmax(custom_booster.predict(X_test))

assert (multiclass_pred == custom_pred).all()
```

</details>


## まとめ

- 二値分類、多クラス分類の目的関数 Cross Entropy loss の定義とその gradient、hessian の計算結果を紹介
- LightGBM での目的関数の指定方法を紹介
- **実装上の注意点**
  - カスタム目的関数を与えると変数変換の後処理が自動で行われないので、忘れずに自分で実装する
  - 組み込みの目的関数ではデフォルトで `boost_from_average=True` となっているが、カスタム関数ではこれが適用されないため学習結果に差異が生じる
  - 組み込みの多クラス分類`"multiclass"` は Log loss の hessian に $K/(K-1)$ 倍の補正ファクターがかけられている
  - sklearn API と train API は目的関数に与えられる引数が異なる

おそらく、これで躓きポイントは一通りまとめられたと思います。
実装の参考となれば幸いです。

## 参考

- [Focal loss implementation for LightGBM • Max Halford](https://maxhalford.github.io/blog/lightgbm-focal-loss/)
    - カスタムobjective実装について最も詳しい解説。LightGBM v3 系なので一部APIが異なることに注意が必要だが、`boost_from_average`の扱いやmetricsの実装を含め詳しく解説されている
- 補足：目的関数の導出ノート
    - [log lossとその勾配の導出｜Notion](https://intriguing-kitchen-267.notion.site/log-loss-14647086bf81802ca6bcc418597bebd9)
    - [多値Cross Entropyの勾配の導出｜Notion](https://intriguing-kitchen-267.notion.site/Cross-Entropy-15047086bf8180968c40f5a7bb4e469a)
- custom objective の公式説明
    - [lightgbm.LGBMClassifier — LightGBM 4.5.0 documentationl](https://lightgbm.readthedocs.io/en/stable/pythonapi/lightgbm.LGBMClassifier.html)
    - [lightgbm.train — LightGBM 4.5.0 documentation](https://lightgbm.readthedocs.io/en/stable/pythonapi/lightgbm.train.html#lightgbm.train:~:text=for%20more%20information.-,Note,-A%20custom%20objective)
- 先行研究（LightGBM v4 以前のためAPIが異なることに注意）
    - [LightGBM で強引に Multi-Task(は???) Regression を行う - 俵言](https://tawara.hatenablog.com/entry/2020/05/14/120016)
    - [LightGBMのObjectiveをpythonで再現する #Python - Qiita](https://qiita.com/nadare/items/245ebf265bc64c5d7996)
    - [LighGBMでcustom objectを使う｜Zenn](https://zenn.dev/ymd/articles/b446e4b04814d3)
    - [How to implement custom loss function in scikit-learn? | Interviews Vector](https://www.interviewsvector.com/blog/how-implement-custom-loss-function-in-scikit-learn)
- 理論解説
    - Gradient Boost や Newton Boost について [Gradient Boosting と XGBoost｜Zenn](https://zenn.dev/zaburo_ch/articles/71b1af9bba1571)
    - [Gradient Boost Modelにおける多クラス分類の目的関数について](https://zenn.dev/bilzard/articles/obj-func-of-multi-class-on-gbm)
    - [Custom Objective for LightGBM | Hippocampus's Garden](https://hippocampus-garden.com/lgbm_custom/)
- eval に使われる metric 関数についてもカスタム関数を設定することができる： [Python: LightGBM でカスタムメトリックを扱う - CUBE SUGAR CONTAINER](https://blog.amedama.jp/entry/lightgbm-custom-metric)
- XGBoost でも同様のことができる [Custom Objective and Evaluation Metric — xgboost 2.1.1 documentation](https://xgboost.readthedocs.io/en/stable/tutorials/custom_metric_obj.html)