LightGBM で強引に Multi-Task(は???) Regression を行う

!!注意!!：この記事の内容は率直に言って全く褒められた行為ではありません。それをご了承の上でネタとしてご覧ください。

はじめに

前回*1に引き続き今回も LightGBM のお話です*2。前回はちょい使いどころがある話でしたが、今回に関しては本当に誰得ネタだと思います。

僕自身は DSB2019 の反省*3 *4をしている際に Regression 用の Custom Objective の作り方を学び、前回の記事を書いている中で Binary Classification についても作り方を学びました*5 。その流れで「んじゃ今度は Multi-Class でもやってみるか～」となり、更にその流れで「Multi-Task はどうなんや？」となってこの記事を書くに至りました。タイトルで「Multi-Task(は???) Regression」となっているのはこの表現が合っているのか自信が無い & 実装上の理由*6があって、これについては後程説明します。

先におことわりしておくと、僕は GBDT の仕組みを何となくしか理解しておらずあくまで外部仕様を把握して使っているという状態です。なので、この記事が誤った内容や不正確な内容を含んでいる可能性は多分にあり、またそれを非常に恐れています。読む方は疑いの目線を持って読んで頂きつつ、誤りを発見してもどうか...どうか優しくしてください😇

はじめに
準備：Multi-Class Classification 用に Custom Objective を作る
本題：じゃあ Multi-Task はどう？
これって意味あるの？
本当に... 本当に意味は無い？
おわりに
参考： GBDT の multi-class classification に関する情報
おまけ：他に方法はあるの？

今回使用しているライブラリの import

from typing import List, Tuple, Optional

import numpy as np
import pandas as pd
import category_encoders as ce
from sklearn.model_selection import StratifiedKFold

import lightgbm as lgb

準備：Multi-Class Classification 用に Custom Objective を作る

実装の説明をする前に公式 docs の Training API *7 の fobj の項でポイントを確認します。

For binary task, the preds is margin. For multi-class task, the preds is group by class_id first, then group by row_id. If you want to get i-th row preds in j-th class, the access way is score[j * num_data + i] and you should group grad and hess in this way as well.

一つ目のポイントは、Custom Objective を用いる場合は帰ってくるのは $[0, 1$ ] の範囲の値では無く $(- \infty, +\infty )$ の範囲の margin であるということです。よって渡す自作関数の中で確率として扱う値に変換する必要があり、微分もこれを踏まえた上で行なわければなりません*8。
因みに objective に組み込みのものを使用して metric に自作関数を使う場合は [0, 1] の範囲の値が返ってくる想定で実装する必要があります。実にややこしい😇

二つ目に関しては自作関数が受け取るのは 1D-array だよという話で、要するに自作関数の中で以下のような変形が求められます*9。

f:id:Tawara:20200511191202p:plain — Custom Objective において必要な preds への変形処理

transpose しなくても処理自体は可能ですが、わかりやすいのでやることにしています。

以上の二つを踏まえて Custom Objective を実装します。今回は組み込みとして実装されている Multi-Class LogLoss を練習がてら実装しました。DSB のときと同様に class として定義し、loss 及び grad & hess を返すメソッドをそれぞれ用意していますが、全体がめちゃ長いのでここでは grad & hess に関する部分の一部を抜粋します。

class MultiLoglossForLGBM:
    """Self-made multi-class logloss for LightGBM."""

    def __init__(self, n_class: int=3, use_softmax: bool=True, epsilon: float=1e-32) -> None:
        """Initialize."""
        self.name = "my_mlnloss"
        self.n_class = n_class
        self.prob_func = self._get_prob_value if use_softmax else lambda x: x
        self.epsilon = epsilon
    
    # -------------- ( 省略 ) -------------- #

    def _calc_grad_and_hess(
        self, preds: np.ndarray, labels: np.ndarray, weight: Optional[np.ndarray]=None
    ) -> Tuple[np.ndarray]:
        """Calc Grad and Hess"""
        # # get prob value by softmax
        prob = self.prob_func(preds)           # <= margin を確率値に直す
        # # convert labels to 1-hot
        labels = self._get_1hot_label(labels)  # <= labels (1D-array) を 1hot (2D-array) に変換

        grad = prob - labels
        hess = prob * (1 - prob)
        if weight is not None:
            grad = grad * weight[:, None]
            hess = hess * weight[:, None]
        return grad, hess
    
    # -------------- ( 省略 ) -------------- #
    
    def return_grad_and_hess(self, preds: np.ndarray, data: lgb.Dataset) -> Tuple[np.ndarray]:
        """Return Grad and Hess for lightgbm"""
        labels = data.get_label()
        weight = data.get_weight()
        n_example = len(labels)
        
        # # reshape preds: (n_class * n_example,) => (n_class, n_example) =>  (n_example, n_class)
        preds = preds.reshape(self.n_class, n_example).T  # <= preds (1D-array) を 2D-array に直す
        # # calc grad and hess.
        grad, hess =  self._calc_grad_and_hess(preds, labels, weight)

        # # reshape grad, hess: (n_example, n_class) => (n_class, n_example) => (n_class * n_example,) 
        grad = grad.T.reshape(n_example * self.n_class)   # <= 1D-array に戻す
        hess = hess.T.reshape(n_example * self.n_class)   # <= 1D-array に戻す
        
        return grad, hess
    
    # -------------- ( 省略 ) -------------- #

lightgbm.train へ fobj として渡すのは return_grad_and_hess という名前のメソッドです。こいつの中で、受け取った preds を 2D-array に直す操作と、計算した grad と hess を 1D-array に戻す操作を行っているのが分かると思います。内部的なメソッドである _calc_grad_and_hess の中では、preds を確率値 (prob)に変換(softmax を適用)する処理と labels を 1-hot に変換する処理を行った後に grad と hess を計算しています。

実装の一部しか載せていないので説明していない部分がありますが、Notebook を最後に載せておくので興味があればそちらをご確認ください。Gradient と Hessian の計算についてもそちらに載せてあります。

というわけで早速試してみましょう。今回の検証は Multi-Class Classification の例題として有名な Iris Dataset*10 を使って行います。以下のような形のデータ(train_all は csv を pandas で読み込んだもの)で、Label(Target) である Speicies には３種類の Class ('Iris-setosa', 'Iris-versicolor', 'Iris-virginica') があります。

train_all.head()

	Id	SepalLengthCm	SepalWidthCm	PetalLengthCm	PetalWidthCm	Species
0	1	5.1	3.5	1.4	0.2	Iris-setosa
1	2	4.9	3	1.4	0.2	Iris-setosa
2	3	4.7	3.2	1.3	0.2	Iris-setosa
3	4	4.6	3.1	1.5	0.2	Iris-setosa
4	5	5	3.6	1.4	0.2	Iris-setosa

Speicies に対して Label Encoding を行ったのち、適当に分割して先程の自作関数を適用しましょう。

ord_enc = ce.OrdinalEncoder(
     mapping=[{"col": TARGET, "mapping": {c: i for i, c in enumerate(train_all.Species.unique())}}],
    cols=[TARGET])
train_all_multi_class = ord_enc.fit_transform(train_all)
train_all_multi_class.head()

	Id	SepalLengthCm	SepalWidthCm	PetalLengthCm	PetalWidthCm
0	1	5.1	3.5	1.4	0.2
1	2	4.9	3	1.4	0.2
2	3	4.7	3.2	1.3	0.2
3	4	4.6	3.1	1.5	0.2
4	5	5	3.6	1.4	0.2

scikit-learn の StratifiedKFold を用いて分割します。

X_0 = train_all_multi_class.iloc[:,1:5].values
y_0 = train_all_multi_class.iloc[:,-1].values

kf = StratifiedKFold(n_splits=N_FOLD, random_state=RANDOM_SEED, shuffle=True)
train_val_splits = list(kf.split(X_0, y_0))

# # use fold 0
train_index, valid_index = train_val_splits[0]
X_0_tr, y_0_tr = X_0[train_index], y_0[train_index]
X_0_val, y_0_val =  X_0[valid_index], y_0[valid_index]

データの準備が出来たので、自作した MultiLoglossForLGBM のインスタンスを生成し fobj と feval にそれぞれ return_grad_and_hess と return_loss を渡して学習を行います。lightgbm.train に渡す params 引数については、 num_class を指定することさえ注意すれば行けるっぽいです。

MODEL_PARAMS_LGB = {
    'num_class': 3,  # <= class 数を指定
    "metric": "None",
    "first_metric_only": True,
    "eta": 0.01,
    "max_depth": -1,
    "seed": RANDOM_SEED,
    "num_threads": NUM_THREADS,
    "verbose": -1
}
FIT_PARAMS_LGB = {
    "num_boost_round": 10000,
    "early_stopping_rounds": 100,
    "verbose_eval":100}

my_mlnloss = MultiLoglossForLGBM(n_class=3, use_softmax=True)  # <= 自作 class を初期化 

lgb_tr = lgb.Dataset(X_0_tr, y_0_tr)
lgb_val = lgb.Dataset(X_0_val, y_0_val)

model_mymlnloss = lgb.train(
    params=MODEL_PARAMS_LGB, train_set=lgb_tr, **FIT_PARAMS_LGB,
    valid_names=['train', 'valid'], valid_sets=[lgb_tr, lgb_val],
    fobj=my_mlnloss.return_grad_and_hess,           # <= gradient と hessian を返す関数
    feval=lambda preds, data: [
        my_mlnloss.return_loss(preds, data),        # <= loss を返す関数
        multi_class_accuracy_for_lgbm(preds, data)  # <= 評価用の multi-class accuracy (自作) 
    ]
)

以下のようにちゃんと学習されます。

Training until validation scores don't improve for 100 rounds
[100]   train's my_mlnloss: 0.208129    train's my_macc: 0.983333   valid's my_mlnloss: 0.214356    valid's my_macc: 0.966667
[200]   train's my_mlnloss: 0.0847299   train's my_macc: 0.983333   valid's my_mlnloss: 0.105429    valid's my_macc: 0.966667
[300]   train's my_mlnloss: 0.047182    train's my_macc: 0.991667   valid's my_mlnloss: 0.0967384   valid's my_macc: 0.966667
Early stopping, best iteration is:
[280]   train's my_mlnloss: 0.0527208   train's my_macc: 0.991667   valid's my_mlnloss: 0.0962209   valid's my_macc: 0.966667
Evaluated only: my_mlnloss

因みに組み込みの multiclass を使用した場合の結果は以下。

MODEL_PARAMS_LGB = {
    "objective": "multiclass",  # <= set implemented multi logloss
    'num_class': 3,
    "first_metric_only": True,
    "eta": 0.01,
    "max_depth": -1,
    "seed": RANDOM_SEED,
    "num_threads": NUM_THREADS,
    "verbose": -1
}
FIT_PARAMS_LGB = {
    "num_boost_round": 10000,
    "early_stopping_rounds": 100,
    "verbose_eval":100}

my_mlnloss = MultiLoglossForLGBM(n_class=3, use_softmax=False)   # <= objective が組み込みなので softmax は使用しない.
lgb_tr = lgb.Dataset(X_0_tr, y_0_tr)
lgb_val = lgb.Dataset(X_0_val, y_0_val)

model = lgb.train(
    params=MODEL_PARAMS_LGB, train_set=lgb_tr, **FIT_PARAMS_LGB,
    valid_names=['train', 'valid'], valid_sets=[lgb_tr, lgb_val],
    feval=lambda preds, data: [
        my_mlnloss.return_loss(preds, data),        # <= multi-class logloss  (自作)
        multi_class_accuracy_for_lgbm(preds, data)  # <= multi-class accuracy (自作) 
    ]
)

比較のために自作の multi-class logloss も metric として feval に設定しています。

Training until validation scores don't improve for 100 rounds
[100]   train's multi_logloss: 0.420851     train's my_mlnloss: 0.420851    train's my_macc: 0.966667   valid's multi_logloss: 0.419136     valid's my_mlnloss: 0.419136    valid's my_macc: 0.966667
[200]   train's multi_logloss: 0.208767     train's my_mlnloss: 0.208767    train's my_macc: 0.983333   valid's multi_logloss: 0.214871     valid's my_mlnloss: 0.214871    valid's my_macc: 0.966667
[300]   train's multi_logloss: 0.123661     train's my_mlnloss: 0.123661    train's my_macc: 0.975      valid's multi_logloss: 0.136866     valid's my_mlnloss: 0.136866    valid's my_macc: 0.966667
[400]   train's multi_logloss: 0.0850451    train's my_mlnloss: 0.0850451   train's my_macc: 0.983333   valid's multi_logloss: 0.105641     valid's my_mlnloss: 0.105641    valid's my_macc: 0.966667
[500]   train's multi_logloss: 0.0631742    train's my_mlnloss: 0.0631742   train's my_macc: 0.991667   valid's multi_logloss: 0.0983892    valid's my_mlnloss: 0.0983892   valid's my_macc: 0.966667
[600]   train's multi_logloss: 0.0474029    train's my_mlnloss: 0.0474029   train's my_macc: 0.991667   valid's multi_logloss: 0.0969201    valid's my_mlnloss: 0.0969201   valid's my_macc: 0.966667
Early stopping, best iteration is:
[553]   train's multi_logloss: 0.0539797    train's my_mlnloss: 0.0539797   train's my_macc: 0.991667   valid's multi_logloss: 0.0964299    valid's my_mlnloss: 0.0964299   valid's my_macc: 0.966667
Evaluated only: multi_logloss

metric については自作と組み込みで値が一致していますね。ただ、学習の挙動が自前実装と微妙に違うような？🤔
少し調べましたが「組み込みの方が C++ で実装してるから違いが出るんじゃない？」みたいな issue *11とか、

github.com

また別の Issue があったりしました。

github.com

一応学習は出来てるし実装は間違っていないと信じたいところですが...

とりあえず学習は出来たということにして次に進みます。

本題：じゃあ Multi-Task はどう？

さて、前回のような小ネタ記事であれば準備の時点で(文字数的に)終了なので既にお腹一杯かもしれませんが、残念ながらここからが本番です。

前述したように、num_class さえ指定すれば Class 数分の予測値(margin) が得られることが分かっています。ならば、Custom Objective さえうまく作ってやれば Multi-Task の回帰だって形式的には出来るのではないでしょうか？

一つ問題になるのは Label の与え方です。複数の予測対象に対する実数値を事例ごとに渡してあげる必要がある(なので 2D-array で渡したい)のですが、 lightgbm.Dataset に 1D-array 以外を渡すと(そのときは怒られ無いのですが)学習を行うために lightgbm.Booster を初期化する時点でお叱りを受けます。

# -------------------- (前略)  -------------------- #

TypeError: Wrong type(ndarray) for label.
It should be list, numpy 1-D array or pandas Series

「んじゃ 2D-array を flatten して渡したらいいのでは？」とやったら「長さ揃えろや」って怒られます。(それはそう)

# -------------------- (前略) -------------------- #

LightGBMError: Length of label is not same with #data

これについての回避は可能で、自作関数の内部が無法地帯なので色々案がありました。ただ、よくよく考えると最もシンプルな方法は lightgbm.Dataset にメンバ変数として持たせることです。

lgb_tr = lgb.Dataset(X_1_tr, np.arange(X_1_tr))  # <= dummy の label を与えておく (もしくは何も与えない)
lgb_tr.multi_label = y_1_tr                      # <= 本当のラベルを持たせる

これであとは自作関数の中で data.multi_label にアクセスすれば解決出来ますが、流石にお行儀が悪い気がするので lgb.Dataset を継承した class を定義することにします。とはいっても multi_label を持たせることだけが目的なのでめちゃくちゃ簡易なものであり、またあくまで python 側で動く Custom Objective のためだけの class であることにご注意ください。

class MultiLabelDatasetForLGBM(lgb.Dataset):
    """
    Makeshift Class for storing multi label.
    
    label: numpy.ndarray (n_example, n_target)
    """

    def __init__(
        self, data, label=None, reference=None, weight=None, group=None, init_score=None, silent=False,
        feature_name='auto', categorical_feature='auto', params=None,  free_raw_data=True
    ):
        """Initialize."""
        if label is not None:
            # # make dummy 1D-array
            dummy_label = np.arange(len(data))

        super(MultiLabelDatasetForLGBM, self).__init__(
            data, dummy_label, reference, weight, group, init_score, silent,
            feature_name, categorical_feature, params, free_raw_data)
        
        self.mult_label = label
        
    def get_multi_label(self):
        """Get 2D-array label"""
        return self.mult_label
    
    def set_multi_label(self, multi_label: np.ndarray):
        """Set 2D-array label"""
        self.mult_label = multi_label
        return self

渡した label は multi_label として持ち、lightgbm 側が label として認識する奴には怒られないように dummy を渡しています。あとは自作関数の中で data.get_multi_label() を呼べばよいだけです。

さて、肝心の Custom Objective/Metric ですが、今回はシンプルにそれぞれの Task の二乗誤差を足したものを考えます。書くと長そうなので式から察してください。

$N: \verb|number of examples| \\ M: \verb|number of targets (tasks)| \\ Y: \verb|matrix of predicted values| \\ T: \verb|matrix representing target values| \\ \displaystyle Y = \left ( \vec{y}_1, \vec{y}_2, ... , \vec{y}_N\right )^\mathsf{T} \verb|, where | \ \vec{y}_i = \left( y_{i 1}, y_{i 2}, ... , y_{i M} \right)^\mathsf{T} \\ \displaystyle T = \left ( \vec{t}_1, \vec{t}_2, ... , \vec{t}_N\right )^\mathsf{T} \verb|, where | \displaystyle \vec{t}_i = \left( t_{i 1}, t_{i 2}, ... , t_{i M}\right)^\mathsf{T}$

$\displaystyle \mathrm{MultiTaskMSE} \left(Y, T \right) = \frac{1}{N} \sum_{i=1}^N 2 \ \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) \\ \displaystyle \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = \sum_{j=1}^M \frac{1}{2} \left( y_{ij} - t_{ij} \right)^2 \\ \displaystyle \frac{\partial}{\partial y_{ij}} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = y_{ij} - t_{ij} \\ \displaystyle \frac{\partial^2}{\partial y_{ij}^2} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = 1$

$\mathrm{loss}$ が事例ごとの損失を計算する関数で、これを $y_{ij}$ で偏微分しています。言い換えると、事例 $i$ の target $j$ に対する予測を用いて事例 $i$ に対する損失を偏微分しています。自作関数が返す grad, hess はその $\left(i, j \right)$ 成分がそれぞれ $\displaystyle \frac{\partial}{\partial y_{ij}} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right)$ , $\displaystyle \frac{\partial^2}{\partial y_{ij}^2} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right)$ であるような行列(を lightgbm が所望する形式の 1D-array に変換したもの) となります。

というわけで実装(の一部を抜粋したもの)は以下の通りです。

class MultiMSEForLGBM:
    """Self-made multi-task(?) mse for LightGBM."""
    
    def __init__(self, n_target: int=3) -> None:
        """Initialize."""
        self.name = "my_mmse"
        self.n_target = n_target

    # -------------- ( 省略 ) -------------- #

    def _calc_grad_and_hess(
        self, preds: np.ndarray, labels: np.ndarray, weight: Optional[np.ndarray]=None
    ) -> Tuple[np.ndarray]:
        """Calc Grad and Hess"""
        grad = preds - labels
        hess = np.ones_like(preds)     
        if weight is not None:
            grad = grad * weight[:, None]
            hess = hess * weight[:, None]

        return grad, hess

    # -------------- ( 省略 ) -------------- #

    def return_grad_and_hess(self, preds: np.ndarray, data: lgb.Dataset) -> Tuple[np.ndarray]:
        """Return Grad and Hess for lightgbm"""
        labels = data.get_multi_label()  # <= 改造した Dataset から multi-label を受け取る
        weight = data.get_weight()
        n_example = len(labels)
        
        # # reshape preds: (n_target * n_example,) => (n_target, n_example) =>  (n_example, n_target)
        preds = preds.reshape(self.n_target, n_example).T  # <= preds (1D-array) を 2D-array に直す
        # # calc grad and hess.
        grad, hess =  self._calc_grad_and_hess(preds, labels, weight)

        # # reshape grad, hess: (n_example, n_target) => (n_class, n_target) => (n_target * n_example,) 
        grad = grad.T.reshape(n_example * self.n_target)   # <= 1D-array に戻す
        hess = hess.T.reshape(n_example * self.n_target)   # <= 1D-array に戻す
        
        return grad, hess

準備で紹介した MultiLoglossForLGBM と grad や hess の計算はもちろん異なりますが、.get_multi_label() で label を受け取る以外の流れはほぼほぼ一緒です。

これを用いて、非常に強引ですが先程の Iris Dataset を Multi-Task Regression として解いてみます。３つの Class それぞれについて回帰で予測を行うという設定です。
先程は Label-Encoding をしていましたが、今度は OneHot-Encoding を行って 2D-array の label を作ります。

ohot_enc = ce.OneHotEncoder(cols=[TARGET], use_cat_names=True)
train_all_multi_reg = ohot_enc.fit_transform(train_all)
train_all_multi_reg.head()

	Id	SepalLengthCm	SepalWidthCm	PetalLengthCm	PetalWidthCm	Species_Iris-setosa
0	1	5.1	3.5	1.4	0.2	1
1	2	4.9	3	1.4	0.2	1
2	3	4.7	3.2	1.3	0.2	1
3	4	4.6	3.1	1.5	0.2	1
4	5	5	3.6	1.4	0.2	1

分割については Multi-Class Classification として解いた場合と同様の分割を使用します。

X_1 = train_all_multi_reg.iloc[:, 1:5].values
y_1 = train_all_multi_reg.iloc[:, 5:8].values

# # use the same split as multi-class classification
train_index, valid_index = train_val_splits[1]
X_1_tr, y_1_tr = X_1[train_index], y_1[train_index]
X_1_val, y_1_val =  X_1[valid_index], y_1[valid_index]

いざ、学習です。loss だけだとよくわからないので、metric で accuracy も測るようにしておきます。

MODEL_PARAMS_LGB = {
    'num_class': 3,  # <= class 数を指定
    "eta": 0.01,
    "metric": "None",
    "first_metric_only": True,
    "max_depth": -1,
    "seed": RANDOM_SEED,
    "num_threads": NUM_THREADS,
    "verbose": -1
}
FIT_PARAMS_LGB = {
    "num_boost_round": 10000,
    "early_stopping_rounds": 100,
    "verbose_eval":50}

my_mmse = MultiMSEForLGBM(n_target=3)  # <= 自作 class を初期化 

lgb_tr = MultiLabelDatasetForLGBM(X_1_tr, y_1_tr)     # <= 改造 Dataset
lgb_val = MultiLabelDatasetForLGBM(X_1_val, y_1_val)  # <= 改造 Dataset

model_my_mmse = lgb.train(
    MODEL_PARAMS_LGB, lgb_tr, **FIT_PARAMS_LGB,
    valid_names=['train', 'valid'], valid_sets=[lgb_tr, lgb_val],
    fobj=my_mmse.return_grad_and_hess,                       # <= gradient と hessian を返す関数
    feval=lambda preds, data: [
        my_mmse.return_loss(preds, data),                    # <= loss を返す関数
        multi_class_accuracy_for_lgbm_altered(preds, data),  # <= multi-class accuracy (自作) 
    ]
)

出力は以下。

Training until validation scores don't improve for 100 rounds
[50]    train's my_mmse: 0.40525    train's my_macc: 0.958333   valid's my_mmse: 0.412759   valid's my_macc: 0.966667
[100]   train's my_mmse: 0.179346   train's my_macc: 0.983333   valid's my_mmse: 0.19676    valid's my_macc: 0.966667
[150]   train's my_mmse: 0.0925749  train's my_macc: 0.983333   valid's my_mmse: 0.11856    valid's my_macc: 0.966667
[200]   train's my_mmse: 0.0589131  train's my_macc: 0.983333   valid's my_mmse: 0.0899017  valid's my_macc: 0.966667
[250]   train's my_mmse: 0.0456372  train's my_macc: 0.983333   valid's my_mmse: 0.079527   valid's my_macc: 0.966667
[300]   train's my_mmse: 0.0401923  train's my_macc: 0.983333   valid's my_mmse: 0.0766546  valid's my_macc: 0.966667
[350]   train's my_mmse: 0.0373631  train's my_macc: 0.983333   valid's my_mmse: 0.0767275  valid's my_macc: 0.966667
[400]   train's my_mmse: 0.0354808  train's my_macc: 0.983333   valid's my_mmse: 0.0779402  valid's my_macc: 0.966667
Early stopping, best iteration is:
[330]   train's my_mmse: 0.0384971  train's my_macc: 0.983333   valid's my_mmse: 0.0763893  valid's my_macc: 0.966667
Evaluated only: my_mmse

一応学習は出来ているようです。ただ少し困った点は accuracy の面では Multi-Class Classification として解いたときと差が出ていないことでしょうか。まあ全体で事例数が150しかない(validation が 30しかない) & 簡単なデータなのでので仕方ないっちゃ仕方ないですね。

よし、これで Multi-Task な Regression が... 出来たのでしょうか？

これって意味あるの？

このブログのタイトルを見た段階で、「は??? (何言ってんだコイツ...)」となった方は沢山いらっしゃるはずです。

「複数の Target があって Multi-Task の回帰問題として解く」と聞けば、おそらくほとんどの人は「Multi-Task Learning をするのであろう」と考えると思います。「Multi-Task Learning の正確な定義を教えてください」って言われると正直苦しいですが、例えば NN を用いた Multi-Task Learning*12 *13であれば、ネットワーク構造の一部を Task 間で共有したりするのが分かりやすいと思います。あるいは出力値間での関係性を制約に入れたり。要は Task 間に持たせた何らかの関係性がお互いの学習に(良い)影響を及ぼすのを目的としています。

では、上で無理矢理やった、一つの LightGBM モデルで複数の Task に対する回帰を行ったものはどうでしょうか？

まず、上記の学習は結局 LightGBM の Multi-Class Classification の仕組みに乗っかっています。そしてこれが非常に重要なポイントですが、GBDTの性質上それぞれの木は一つの Class に対して作られます。このため、学習の際には iteration 数 × Class 数分の木が生成されることになり、つまるところ NN の Mulit-Task Learning の例で出した「モデルの一部を Task 間で共有する」といったことは行われていません。 (参考情報は下の方にまとめたのでそちらをご覧ください。教えて頂いた皆様ありがとうございました！)。

モデルを共有していないことはわかりました。それでは Task に対する出力値間に制約があったりするのでしょうか？上述した計算式のうち、事例ごとの部分だけ再確認してみます。

$\displaystyle \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = \sum_{j=1}^M \frac{1}{2} \left( y_{ij} - t_{ij} \right)^2 \\ \displaystyle \frac{\partial}{\partial y_{ij}} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = y_{ij} - t_{ij} \\ \displaystyle \frac{\partial^2}{\partial y_{ij}^2} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = 1$

LightGBM の学習で用いるのは Gradient と Hessian です。Hessian は定数になっているので特に Gradient に注目しますが、よくよく考えると $\frac{\partial}{\partial y_{ij}} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right)$ (Task $j$ の学習のために計算して返す Gradient) には Task $j$ についての予測値( $y_{ij}$ )と正解( $t_{ij}$ )しか出て来てません。よって Task 間でお互いの予測値が学習に影響を及ぼしてはいないことになります。なんてこったい...😱

僕が行ったことは、外側から見ると一つのモデルで複数の回帰問題を解いているのですが、ただ並列に学習を行っているだけで Multi-Task Learning 的なうまみは全くないんですよね。悲しい。メタ的には Hyper Parameter (木の数など)を共有というか揃えることになりますが、この Loss だと別々にやっちゃっていいよなあ...

因みに「じゃあ Multi-Class Classification はどうなの？」って思う方もいらっしゃるかもしれません。一つの木が一つの Class にしか対応せずモデル構造を共有していない点では同じですからね。ただ、明確に違う点は Custom Objective を作る際に必要だと言っていた確率に直す操作( softmax )の存在です。 Multi-Class Logloss の計算を確認してみます。

$N: \verb|number of examples| \\ K: \verb|number of classes| \\ Y: \verb|matrix of predicted values| \\ T: \verb|matrix representing labels| \\ \displaystyle Y = \left ( \vec{y}_1, \vec{y}_2, ... , \vec{y}_N\right )^\mathsf{T} \verb|, where | \ \vec{y}_i = \left( y_{i 1}, y_{i 2}, ... , y_{i K} \right)^\mathsf{T} \\ \displaystyle T = \left ( \vec{t}_1, \vec{t}_2, ... , \vec{t}_N\right )^\mathsf{T} \verb|, where | \displaystyle \vec{t}_i = \left( t_{i 1}, t_{i 2}, ... , t_{i K}\right)^\mathsf{T} \verb| (1-hot vector)|$

$\displaystyle \mathrm{MultiLogLoss}\left(Y, T \right) = \frac{1}{N} \sum_{i = 1}^N \mathrm{loss} \left( \vec{y}_i, \vec{t}_i \right) \\ \displaystyle \mathrm{loss} \left( \vec{y}_i, \vec{t}_i \right) = - \sum_{j = 1}^K t_{ij} \log \left(p_{ij} \right) \\ \displaystyle \frac{\partial}{\partial y_{ij}} \mathrm{loss} \left( \vec{y}_i, \vec{t}_i \right) = p_{ij} - t_{ij} \\ \displaystyle \frac{\partial^2}{\partial y_{ij}^2} \mathrm{loss} \left( \vec{y}_i, \vec{t}_i \right) = p_{ij} \left( 1 - p_{ij} \right) \\ \displaystyle p_{ij} = \mathrm{softmax}_j \left( \vec{y}_{i} \right) = \frac{\exp \left(y_{ij} \right)}{\sum_{j'=1}^K \exp \left(y_{ij'} \right)}$

Gradient と Hessian に出てくる $p_{ij}$ は softmax を通して計算されています。他の Class に対する予測値も学習に影響を与えているわけで、モデル構造を共有していなくても Multi-Class として学習させる意味がここにあります。まあそうでなければ個別に Binary-Classification のモデル学習させればよいだけですもんね。そして僕が苦労してやったやつは全く意味が無かったのであった...。

本当に... 本当に意味は無い？

実は、やりようによってはうまく使えるケースがあるんんじゃないかと思っています。モデルを共有していなくても Task 間で予測値が学習に影響を及ぼしあえばいいわけです。

上で Iris Dataset を Multi-Task(?) Regression で強引に解いたとき実質的には個別に学習が行われちゃってますが、本当は 3つの Class には関係性があります。本来の問題設定では各事例はいずれか一つの Class に分類されるので、当然ながら「予測値を足したら 1 になる」という制約が存在するわけです。これを損失関数に追加しましょう。

$\displaystyle \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = \sum_{j=1}^M \frac{1}{2} \left( y_{ij} - t_{ij} \right)^2 + \frac{1}{2} \left( 1 - \sum_{j=1}^M y_{ij} \right)^2 \\ \displaystyle\frac{d}{d y_{ij}} \mathrm{loss}(\vec{y}_i, \vec{t}_i) = \left( y_{ij} - t_{ij} \right) - \left( 1 - \sum_{k=1}^M y_{ik} \right) = 2 y_{ij} - t_{ij} + \sum_{k \neq j} y_{ik} - 1 \\ \displaystyle \frac{\partial^2}{\partial y_{ij}^2} \mathrm{loss} \left(\vec{y}_i, \vec{t}_i \right) = 2$

「えーでも第二項の制約のせいで変なことになったりしない？」という声が聞こえるので、ハイパラで重み付けするのもありかもしれないですね。

見ればわかる通り、追加した第二項の制約を通して他の Task の予測値が Gradient に影響を及ぼすことになります。これであれば、Multi-Task(?) Regression をすることに意味が出るはずです。
というわけでこのパターンもやってみました(実装に関しては略)。出力だけ載せます。

Training until validation scores don't improve for 100 rounds
[50]    train's my_mmse_2: 0.607067     train's my_mmse: 0.474448   train's my_macc: 0.958333   valid's my_mmse_2: 0.609939     valid's my_mmse: 0.477332   valid's my_macc: 0.966667
[100]   train's my_mmse_2: 0.307231     train's my_mmse: 0.289642   train's my_macc: 0.958333   valid's my_mmse_2: 0.314787     valid's my_mmse: 0.297199   valid's my_macc: 0.966667
[150]   train's my_mmse_2: 0.19553      train's my_mmse: 0.193196   train's my_macc: 0.983333   valid's my_mmse_2: 0.207681     valid's my_mmse: 0.205337   valid's my_macc: 0.966667
[200]   train's my_mmse_2: 0.135525     train's my_mmse: 0.135206   train's my_macc: 0.983333   valid's my_mmse_2: 0.153247     valid's my_mmse: 0.152911   valid's my_macc: 0.966667
[250]   train's my_mmse_2: 0.0992235    train's my_mmse: 0.0991635  train's my_macc: 0.983333   valid's my_mmse_2: 0.121746     valid's my_mmse: 0.121674   valid's my_macc: 0.966667
[300]   train's my_mmse_2: 0.0765443    train's my_mmse: 0.0765146  train's my_macc: 0.983333   valid's my_mmse_2: 0.102669     valid's my_mmse: 0.102635   valid's my_macc: 0.966667
[350]   train's my_mmse_2: 0.0622024    train's my_mmse: 0.0621767  train's my_macc: 0.983333   valid's my_mmse_2: 0.0907846    valid's my_mmse: 0.0907574  valid's my_macc: 0.966667
[400]   train's my_mmse_2: 0.0531537    train's my_mmse: 0.0531241  train's my_macc: 0.983333   valid's my_mmse_2: 0.0838569    valid's my_mmse: 0.0838261  valid's my_macc: 0.966667
[450]   train's my_mmse_2: 0.0473215    train's my_mmse: 0.0472925  train's my_macc: 0.983333   valid's my_mmse_2: 0.0796961    valid's my_mmse: 0.0796606  valid's my_macc: 0.966667
[500]   train's my_mmse_2: 0.0435603    train's my_mmse: 0.0435309  train's my_macc: 0.983333   valid's my_mmse_2: 0.0772807    valid's my_mmse: 0.0772427  valid's my_macc: 0.966667
[550]   train's my_mmse_2: 0.0410891    train's my_mmse: 0.0410623  train's my_macc: 0.983333   valid's my_mmse_2: 0.0761265    valid's my_mmse: 0.0760921  valid's my_macc: 0.966667
[600]   train's my_mmse_2: 0.0394512    train's my_mmse: 0.0394266  train's my_macc: 0.983333   valid's my_mmse_2: 0.0759097    valid's my_mmse: 0.0758785  valid's my_macc: 0.966667
[650]   train's my_mmse_2: 0.0383218    train's my_mmse: 0.0382973  train's my_macc: 0.983333   valid's my_mmse_2: 0.0760149    valid's my_mmse: 0.0759848  valid's my_macc: 0.966667
Early stopping, best iteration is:
[594]   train's my_mmse_2: 0.0396163    train's my_mmse: 0.0395913  train's my_macc: 0.983333   valid's my_mmse_2: 0.075896 valid's my_mmse: 0.0758644  valid's my_macc: 0.966667
Evaluated only: my_mmse_2

my_mmse_2 が制約を追加したもの、my_mmse が追加前のものです。数字が色々出て来て分かりにくいですが伝えたいこととしては以下。

とりあえず同等の validation accuracy は出ている(もちろん少数データなのであくまで目安)
当初は my_mmse_2 の方が数字が大きいが、学習が進むにつれて my_mmse に近づく
軽微な差であるが、制約なしで学習させた時よりも validation の my_mmse の値が低い

一番言いたいのは 3 です。といっても「validation が良くなりました！」と主張したいわけでは無く(差が微妙過ぎて何とも言えない)、追加した制約によって単純に並列で学習させるのとは挙動が変わったというところがポイントになります。
Multi-Class Classification を強引に解いているので、もうちょっと適した問題でやってみると面白いかもしれないですね。

おわりに

はい、実装とかしてる間は楽しかったのですが、書いてる途中から「これGBDTガチ勢に怒られるのでは...?」と恐怖しながら書いていた記事でした。正直未だに内部実装を理解してないので、間違ったことを言っていないかはかなり不安が残っています。

num_class * num_iteration の木を作ることから、 GBDT で Multi-Class Classification を解くのは微妙かもしれません(特にClass 数が膨大なとき)。今回やってみた Multi-Task(は???) Regression とかは別々に学習させればいいのではと言われると思いますし、僕もそう思っています。

一方で、最後にちょろっとやってみたように Regression Task 間で何らかの制約がある場合は面白いことが出来るかもしれません。勿論 Task 数がそんなに多くないという前提の下ではありますが、複数の Regression Task について「 $\mathrm{Task} \ A + \mathrm{Task} \ B = \mathrm{Task} \ C$ 」の関係( Task の階層構造)とか「 $|| \mathrm{Task} \ D - \mathrm{Task} \ E || \leq \mathrm{const.}$ 」の関係( Task の類似性)があった場合に、別個に Regression Task として解くよりも今回みたいなアプローチは面白い気がします*14。

以下が今回の実装などを載せた Notebook です。Custom Objective の実装は一通りやった(はず)と思うので GitHub とかにまとめても良い気がしますが、面倒なので気が向いたらやります。

www.kaggle.com

息抜きも兼ねて最近こまめにブログ更新してましたが今回のは流石に疲れました...。頃合いなのでそろそろ Kaggle に復帰しようと思います。
今回のネタもどこかで使える気がしますし、また頑張って行きましょう！

参考： GBDT の multi-class classification に関する情報

僕が丁度良いやつを探し当てられず Twitter で help を求めたら色んな方がリプで教えてくれました。

GBDTに詳しい人助けて... (これがないとブログが書き終わらない)https://t.co/iaP3zsAmqh
— イ表 (@tawatawara) May 12, 2020

教えて下さった threecourse(id:threecourse)さん、杏仁まぜそば(id:aotamasaki)さん、Nomi(id:nyanp)さん、まますさん、nyker_goto(id:dette) さん、marugari(id:marugari2)さん、ありがとうございました！

クラス数分の木を作ることに言及していた貴重な日本語記事

marugari2.hatenablog.jp

上記の Tweet も「なんか日本語で言及してる記事があった気がするけど思い出せね～」とつぶやいたのでした。ありがとうございました。

GitHub Issue「なんで multi-class classification のときは num_class * num_iterations の木を作るんですか？」

github.com

回答としては以下。

Because there is not a good multi-output tree, thus, it uses multi-trees to output multi values.

XGBoost の実装(iteration ごとの木を作る部分)

github.com

ngroup = model_.param.num_output_group がクラス数に相当。if 文で分岐しており、for (int gid = 0; gid < ngroup; ++gid) の中で gid ごとに木を作って追加してますね。
gpair_h へのアクセスの仕方を見て気付きましたが、XGBoost と LightGBM だと Multi-Class classification における preds, grad, hess の持ち方が異なる...?

調べてみた所以下の example code に誘導されるのですが、

github.com

受け取る preds は (n_example, n_class) の 2D-array で、grad と hess は 2D-array (preds と同形式のものを (n_example * n_class, 1) に reshape したもの ) として返しているように見えます。公式の docs*15で明言してないのもあって実際に自分でやってみた方が良さそう...。
何か変な感じですが、とりあえずLightGBM とは少しだけ実装を変える必要がありますね。いい気付きを得ました。

因みに LightGBM のそれっぽい部分だと、

github.com

num_tree_per_iteration_ が Class 数の入っている変数なのですが、for (int cur_tree_id = 0; cur_tree_id < num_tree_per_iteration_; ++cur_tree_id) の中で const size_t offset = static_cast<size_t>(cur_tree_id) * num_data_ を計算し、Class 毎に offset 分 index をずらしてアクセスしているようです(多分)。

シナモン先生：「本質的にはMARTの解説を読めば十分」

web.stanford.edu

まさかこんなきっかけでこの本を読むことになるとは思いませんでした。
教えて頂いた p.387 の Algorithm 10.4 と p.361 の Algorithm 10.3 を比較すると iteration ごとにクラス数(K) 分の木を作っていることが分かりやすいですね。

Gradient boosted decision trees for high dimensional sparse output (ICML2017)

https://dl.acm.org/doi/10.5555/3305890.3306010

主旨的にはクラス数が膨大なときに GBDT を適用するのが厳しいのでそこにうまく対処するお話なのですが、文中で GBDT で Multi-Label を解く際の仕組みを簡潔にまとめてくれているということで紹介いただきました。
Extreme Multi-Labels 問題を (GBDT で)解くためにどうするかという問題は大事だし、Kaggle でもあり得る設定ですね。これ教えてもらったときに以下の発表を思い出しました。

YJTC18 D-4 AnnexML: 近似最近傍検索を⽤いたextreme multi-label分類の⾼速化