【LightGBM】カスタム目的関数(custom objective)で非対称な損失を使う

はじめに

LightGBMはobjective="regression"のように、あらかじめ用意された目的関数(損失関数)を選んで使うのが普通です。 ですが実は、目的関数を自分で書いて渡すこともできます。これが**カスタム目的関数(custom objective)**です。

例えばこんな場面で効いてきます。

  • 需要予測で、在庫切れ(下振れ)は機会損失が大きいので、上振れより重く罰したい
  • 外れ値に引っ張られにくい、頑健な損失を使いたい
  • コンペの評価指標に合わせて、独自の損失を最適化したい

今回はこのカスタム目的関数の仕組みを理解し、「下振れ予測を重く罰する非対称な損失」を実装して、予測が実際にどう変わるかを見ていきます。

仕組み:勾配(grad)とヘシアン(hess)を返すだけ

LightGBMのような勾配ブースティングは、損失関数を予測値で微分した情報を使って木を成長させます。 そのため、カスタム目的関数がやることはシンプルで、予測値に対する損失の「1階微分(grad)」と「2階微分(hess)」を返すだけです。

  • grad(gradient)… 損失の1階微分。「どちらにどれだけ動かせば損失が減るか」
  • hess(hessian)… 損失の2階微分。「その曲がり具合」

具体的には、(preds, dataset)の2つを引数に取る関数を用意します。predsはその時点でのモデルの予測値、datasetは学習データ(正解ラベルはdataset.get_label()で取り出せます)です。 この2つから予測1件ごとのgradhessを計算し、予測と同じ長さの配列2つとして返します。

def my_objective(preds, dataset):
    y = dataset.get_label()
    grad = ...  # 損失をpredsで1階微分したもの
    hess = ...  # 損失をpredsで2階微分したもの
    return grad, hess

注意:渡し方はバージョンで変わっている

少し前まではlgb.train(..., fobj=my_objective)のようにfobj引数で渡していましたが、LightGBM 4.x ではfobj引数は廃止され、paramsobjectiveに関数を渡す形になりました。今回のバージョン(4.6.0)ではこう書きます。

params = dict(objective=my_objective, ...)
model = lgb.train(params, dtrain, num_boost_round=300)

まずは自作でL2損失を再現して検証する

いきなり複雑な損失を書く前に、組み込みのL2損失(regression)を自作で再現してみます。 これが組み込みと一致すれば、grad/hessの実装が正しいと確認できます。

L2損失は loss = 0.5 * (pred - y)^2 なので、微分すると

  • grad = pred - y
  • hess = 1

です。

import numpy as np
import lightgbm as lgb

rng = np.random.RandomState(42)
n = 4000
X = rng.uniform(-3, 3, (n, 3))
y = X[:,0]**2 + 2*X[:,1] + rng.randn(n)*1.0

Xtr, Xte = X[:3000], X[3000:]
ytr, yte = y[:3000], y[3000:]
dtr = lgb.Dataset(Xtr, ytr)

base = dict(num_leaves=31, learning_rate=0.05, min_child_samples=20, verbose=-1)

# 組み込みL2
m_builtin = lgb.train(dict(base, objective="regression"), dtr, num_boost_round=300)
p_builtin = m_builtin.predict(Xte)

# 自作L2: grad = pred - y, hess = 1
def l2_obj(preds, dataset):
    y = dataset.get_label()
    grad = preds - y
    hess = np.ones_like(preds)
    return grad, hess

m_custom = lgb.train(dict(base, objective=l2_obj), dtr, num_boost_round=300)
p_custom = m_custom.predict(Xte)

print("予測の最大絶対差:", np.max(np.abs(p_builtin - p_custom)).round(4))
print("組み込みRMSE:", np.sqrt(np.mean((p_builtin-yte)**2)).round(4))
print("自作RMSE   :", np.sqrt(np.mean((p_custom-yte)**2)).round(4))
予測の最大絶対差: 0.0
組み込みRMSE: 1.0766
自作RMSE   : 1.0766

予測が完全に一致しました(最大差 0.0)。grad/hessの実装が正しいことが確認できたので、安心して独自の損失に進めます。

非対称な損失を実装する

では本題です。「下振れ(pred < y)を上振れより重く罰する」非対称な損失を作ります。 残差 r = pred - y として、下振れ(r < 0)のときだけ重みalphaを掛けたL2にします。

  • loss = 0.5 * w * (pred - y)^2w = alpha(下振れ時) or 1(上振れ時))
  • grad = w * (pred - y)
  • hess = w
def asym_obj(alpha):
    def obj(preds, dataset):
        y = dataset.get_label()
        residual = preds - y                       # >0: 上振れ, <0: 下振れ
        weight = np.where(residual < 0, alpha, 1.0) # 下振れだけ重く
        grad = weight * residual
        hess = weight
        return grad, hess
    return obj

# 下振れを10倍重く罰する
m_asym = lgb.train(dict(base, objective=asym_obj(10.0)), dtr, num_boost_round=300)
p_asym = m_asym.predict(Xte)

for nm, p in [("対称(L2)", p_custom), ("非対称", p_asym)]:
    r = p - yte
    print(f"{nm:8s} 平均残差={r.mean(): .3f}  下振れ率(pred<y)={(r<0).mean()*100:4.1f}%")
対称(L2)   平均残差=-0.024  下振れ率(pred<y)=51.0%
非対称      平均残差= 0.587  下振れ率(pred<y)=30.0%

はっきり効果が出ました。

  • 対称なL2では、平均残差はほぼ0、下振れ率も約50%(上振れ・下振れが半々)
  • 非対称にすると、平均残差が+0.59(全体的に高めに予測するようになった)、下振れ率は30%まで低下

下振れを重く罰した結果、モデルが「罰の大きい下振れを避けるために、あえて高めに予測する」ようになったわけです。

残差の分布をヒストグラムにするとさらにわかりやすいです。

import matplotlib.pyplot as plt

plt.hist(p_custom-yte, bins=40, alpha=0.55, label="symmetric L2")
plt.hist(p_asym-yte, bins=40, alpha=0.55, label="asymmetric (alpha=10)")
plt.axvline(0, color="black", lw=1)
plt.xlabel("prediction - actual (residual)"); plt.ylabel("count")
plt.legend(); plt.show()

lightgbm-custom-objective

オレンジ(非対称)の山が右(上振れ側)にずれていて、0より左の下振れがぐっと減っているのが見て取れます。

2値分類で使うときの落とし穴:predictが確率にならない

ここまでは回帰でしたが、2値分類でカスタム目的関数を使うとpredictの挙動が変わるので、これも実際に確かめておきましょう。

2値分類の基本の損失はロジスティック損失(logloss)です。生スコアrawをシグモイドで確率 p = 1/(1+exp(-raw)) に変換したうえで、微分すると

  • grad = p - y
  • hess = p * (1 - p)

になります。組み込みのbinaryと、この自作版を並べてみます。

rng = np.random.RandomState(0)
n = 4000
X = rng.randn(n, 5)
logit = 1.5*X[:,0] - 2.0*X[:,1] + 0.5*X[:,2]
y = (rng.rand(n) < 1/(1+np.exp(-logit))).astype(float)

Xtr, Xte, ytr, yte = X[:3000], X[3000:], y[:3000], y[3000:]
dtr = lgb.Dataset(Xtr, ytr)
base = dict(num_leaves=31, learning_rate=0.05, min_child_samples=20, verbose=-1)

# 組み込み binary
m_builtin = lgb.train(dict(base, objective="binary"), dtr, num_boost_round=200)
pred_builtin = m_builtin.predict(Xte)

# 自作 logloss: p=sigmoid(raw), grad = p - y, hess = p*(1-p)
def logloss_obj(preds, dataset):
    y = dataset.get_label()
    p = 1/(1+np.exp(-preds))
    grad = p - y
    hess = p*(1-p)
    return grad, hess

m_custom = lgb.train(dict(base, objective=logloss_obj), dtr, num_boost_round=200)
raw_custom = m_custom.predict(Xte)            # そのまま
prob_custom = 1/(1+np.exp(-raw_custom))       # sigmoidをかける

print("組み込みpredict      :", round(pred_builtin.min(),3), "〜", round(pred_builtin.max(),3))
print("自作predict(そのまま) :", round(raw_custom.min(),3), "〜", round(raw_custom.max(),3))
print("自作predict + sigmoid :", round(prob_custom.min(),3), "〜", round(prob_custom.max(),3))
組み込みpredict      : 0.001 〜 1.0
自作predict(そのまま) : -7.577 〜 7.765
自作predict + sigmoid : 0.001 〜 1.0

ここがまさに落とし穴です。

  • 組み込みのbinarypredict0〜1の確率を返す
  • ところが自作版のpredictは**-7.5〜7.8のような生スコア(対数オッズ)**で、確率になっていない
  • 自分でシグモイド1/(1+exp(-x))をかけると、ちゃんと0〜1の確率に戻る

組み込みの分類は、この最後のシグモイド変換を裏でやってくれていたわけです。カスタム目的関数ではそこまで面倒を見てくれないので、確率が欲しければ自分で変換する必要があります。 (なお、0.5を閾値にクラスを判定したいだけなら、確率に直さず生スコアの符号=raw > 0を見るだけでもOKです。)

使うときの注意点

カスタム目的関数を使うときは、いくつか気をつける点があります。

  • 評価指標も自前で用意する必要があります。組み込みのmetricは効かなくなるので、feval引数にカスタム評価関数を渡します。
  • hessは正の値にしておくのが基本です。0や負になると木の分割がうまくいきません(今回のように重みを掛ける形なら安全です)。
  • 組み込みのregressionは「平均から学習を始める」最適化(boost_from_average)が効きますが、カスタム目的関数では0から始まります。今回は300ラウンド回したので差は出ませんでしたが、ラウンド数が少ないと序盤の立ち上がりが遅くなることがあります。

まとめ

  • LightGBMは目的関数を自作でき、やることはgrad(1階微分)とhess(2階微分)を返す関数を書くだけ
  • LightGBM 4.x ではfobj引数は廃止。paramsobjectiveに関数を渡す
  • まず自作でL2を再現して組み込みと一致させれば、実装の正しさを検証できる
  • 非対称な損失にすると、下振れ/上振れのコストの違いを予測に反映できる(今回は下振れを重く罰して予測が上方シフト)
  • 2値分類でカスタム目的関数を使うとpredictは確率でなく生スコア(対数オッズ)を返すので、確率が欲しければ自分でシグモイドをかける

「評価指標に合わせて損失をカスタムしたい」というのはKaggleでもよくある要求です。grad/hessさえ書ければ実現できるので、引き出しに入れておくと強力です。

参考リンク