最近、Twitter のタイムラインで QWK (Quadratic Weighted Kappa: 二次の重み付きカッパ係数) の最適化が話題になっていたので個人的に調べていた。 QWK は順序つきの多値分類問題を評価するための指標で、予測を大きく外すほど大きなペナルティが与えられるようになっている。 また、予測値の分布が結果に影響する点も特徴的で、この点が今回取り扱う最適化にも関係してくる。
QWK の最適化については、Kaggle 本と、その著者 @Maxwell_110 さんによる次のブログエントリが詳しい。 ようするに、真のラベルの分布に沿った形で予測しないと最適な結果が得られない、ということ。
Kaggle 本や、Web をざっと調べた限りでは scipy.optimize.minimize() 関数を使って Nelder Mead 法で最適化するのが一般的なようだった。
今回は、せっかくなので Optuna を使って TPE (Tree-structured Parzen Estimator) で最適化してみた。
使った環境は次のとおり。
$ sw_vers ProductName: Mac OS X ProductVersion: 10.14.6 BuildVersion: 18G3020 $ python -V Python 3.7.6
下準備
はじめに、利用するパッケージをインストールしておく。
$ pip install scikit-learn optuna lightgbm pandas
その上で、Kaggle で開催された次のコンペのデータセットをダウンロードしておく。
これは、前述した Kaggle 本で QWK の最適化が有効に働く例として挙げられているもの。 カレントワーキングディレクトリに、次のファイルがある状態にしておく。
$ ls | grep .zip sample_submission.csv.zip test.csv.zip train.csv.zip
多値分類問題として解いてみる
そもそも、順序つきのラベルを予測するタスクの場合、多値分類問題として解く方法と回帰問題として解く方法がある。 まずは、ベーシックな考え方として多値分類問題として解く方法を試してみる。
早速だけど、サンプルコードは次のとおり。 LightGBM を使って多値分類問題として解いている。 目的関数は Multi LogLoss で、5-Fold CV で QWK について評価している。
#!/usr/bin/env python3 import numbers import pandas as pd import numpy as np import lightgbm as lgb from sklearn.model_selection import KFold from sklearn.metrics import cohen_kappa_score def lgb_custom_metric_qwk_multiclass(preds, data): """LightGBM のカスタムメトリックを計算する関数 多値分類問題として解いた予測から QWK を計算する""" # 正解ラベル y_true = data.get_label() # ラベルの数 num_labels = data.params['num_class'] # 多値分類問題では本来は二次元配列が一次元配列になって提供される reshaped_preds = preds.reshape(num_labels, len(preds) // num_labels) # 最尤と判断したクラスを選ぶ y_pred = np.argmax(reshaped_preds, axis=0) # 最尤と判断したクラスを選ぶ # QWK を計算する return 'qwk', qwk(y_true, y_pred), True def qwk(y_true, y_pred): """QWK (Quadratic Weighted Kappa) を計算する関数""" return cohen_kappa_score(y_true, y_pred, weights='quadratic') def _numerical_only_df(df_): """数値のカラムだけ取り出す関数""" number_cols = [name for name, dtype in df_.dtypes.items() if issubclass(dtype.type, numbers.Number)] numerical_df = df_[number_cols] return numerical_df def main(): # データセットを読み込む train_df = pd.read_csv('train.csv.zip') # 説明変数 x_train = train_df.drop('Response', axis=1) # 簡単にするため、とりあえずすぐに使えるものだけ取り出して使う x_train = _numerical_only_df(x_train) # 目的変数 y_train = train_df.Response.astype(float) # 目的変数はゼロからスタートにする y_train -= 1 # LightGBM のデータセット表現にする lgb_train = lgb.Dataset(x_train, y_train) # 多値分類問題として解く lgbm_params = { 'objective': 'multiclass', 'metric': 'multi_logloss', 'num_class': len(y_train.unique()), 'first_metric_only': True, 'verbose': -1, } # データ分割の方法 folds = KFold(n_splits=5, shuffle=True, random_state=42) # 5-Fold CV で検証する result = lgb.cv(lgbm_params, lgb_train, num_boost_round=1000, early_stopping_rounds=100, folds=folds, verbose_eval=10, feval=lgb_custom_metric_qwk_multiclass, seed=42, ) # early stop したラウンドでの QWK を出力する print(f'CV Mean QWK: {result["qwk-mean"][-1]}') if __name__ == '__main__': main()
上記を実行した結果が次のとおり。
$ python qwkmc.py [10] cv_agg's multi_logloss: 1.37765 + 0.0076621 cv_agg's qwk: 0.442683 + 0.00532416 [20] cv_agg's multi_logloss: 1.25563 + 0.0095355 cv_agg's qwk: 0.507798 + 0.00703211 [30] cv_agg's multi_logloss: 1.20349 + 0.0100597 cv_agg's qwk: 0.519766 + 0.00748977 ...(省略)... [220] cv_agg's multi_logloss: 1.14153 + 0.0126548 cv_agg's qwk: 0.560446 + 0.00516118 [230] cv_agg's multi_logloss: 1.14198 + 0.0126144 cv_agg's qwk: 0.56074 + 0.00444787 [240] cv_agg's multi_logloss: 1.14242 + 0.0130142 cv_agg's qwk: 0.561591 + 0.0046819 CV Mean QWK: 0.5574903141807788
交差検証では QWK の平均として約 0.5574 が得られている。
回帰問題として解いてみる
続いては、同じデータを回帰問題として解くパターンについて。 ラベルが順序つきであれば、回帰問題として解いた上で結果を clip + 四捨五入するような方法もある。
サンプルコードは次のとおり。 やっていることは同じで、解き方が違うだけ。
#!/usr/bin/env python3 import numbers import numpy as np import pandas as pd import lightgbm as lgb from sklearn.model_selection import StratifiedKFold from sklearn.metrics import cohen_kappa_score def lgb_custom_metric_qwk_regression(preds, data): """LightGBM のカスタムメトリックを計算する関数 回帰問題として解いた予測から QWK を計算する""" # 正解ラベル y_true = data.get_label() # 予測ラベル y_pred = np.clip(preds, 0, 7).round() # 単純に予測値を clip して四捨五入する # QWK を計算する return 'qwk', qwk(y_true, y_pred), True def qwk(y_true, y_pred): """QWK (Quadratic Weighted Kappa) を計算する関数""" return cohen_kappa_score(y_true, y_pred, weights='quadratic') def _numerical_only_df(df_): """数値のカラムだけ取り出す関数""" number_cols = [name for name, dtype in df_.dtypes.items() if issubclass(dtype.type, numbers.Number)] numerical_df = df_[number_cols] return numerical_df def main(): # データセットを読み込む train_df = pd.read_csv('train.csv.zip') # 説明変数 x_train = train_df.drop('Response', axis=1) # 簡単にするため、とりあえずすぐに使えるものだけ取り出して使う x_train = _numerical_only_df(x_train) # 目的変数 y_train = train_df.Response.astype(float) # 目的変数はゼロからスタートにする y_train -= 1 # LightGBM のデータセット表現にする lgb_train = lgb.Dataset(x_train, y_train) # 回帰問題として解く lgbm_params = { 'objective': 'regression', 'metric': 'rmse', 'first_metric_only': True, 'verbose': -1, } # データ分割の方法 folds = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) # 5-Fold CV で検証する result = lgb.cv(lgbm_params, lgb_train, num_boost_round=1000, early_stopping_rounds=100, folds=folds, verbose_eval=10, feval=lgb_custom_metric_qwk_regression, seed=42, ) # early stop したラウンドでの QWK を出力する print(f'CV Mean QWK: {result["qwk-mean"][-1]}') if __name__ == '__main__': main()
上記の実行結果は次のとおり。
$ python qwkreg.py [10] cv_agg's rmse: 2.02046 + 0.0085424 cv_agg's qwk: 0.409782 + 0.0042456 [20] cv_agg's rmse: 1.91647 + 0.0125972 cv_agg's qwk: 0.506518 + 0.00716283 [30] cv_agg's rmse: 1.87532 + 0.0140473 cv_agg's qwk: 0.553375 + 0.00638559 ...(省略)... [220] cv_agg's rmse: 1.83657 + 0.0177681 cv_agg's qwk: 0.607652 + 0.00776567 [230] cv_agg's rmse: 1.83694 + 0.0177138 cv_agg's qwk: 0.608083 + 0.00775286 [240] cv_agg's rmse: 1.83728 + 0.0176133 cv_agg's qwk: 0.608366 + 0.00780103 CV Mean QWK: 0.6064040204540543
交差検証では QWK の平均として約 0.6064 が得られている。 先ほどよりもスコアが改善している。
回帰問題として解いた上で OOF Prediction を最適化する
続いては、今回の本題である QWK の閾値の最適化について。 OOF で予測した値を、閾値を最適化することでスコアが改善するかどうか確認してみる。
サンプルコードは次のとおり。
OptunaRounder というクラスを導入して、OOF Prediction を最適化している。
#!/usr/bin/env python3 import numbers import optuna import pandas as pd import numpy as np import lightgbm as lgb from sklearn.model_selection import KFold from sklearn.metrics import cohen_kappa_score class OptunaRounder: """Optuna を使って QWK の最適な閾値を探索するクラス""" def __init__(self, y_true, y_pred): # 真のラベル self.y_true = y_true # 予測したラベル self.y_pred = y_pred # ラベルの種類 self.labels = np.unique(y_true) def __call__(self, trial): """最大化したい目的関数""" # 閾値を Define by run で追加していく thresholds = [] # ラベルの数 - 1 が必要な閾値の数になる for i in range(len(self.labels) - 1): # 閾値の下限 (既存の最大 or ラベルの最小値) low = max(thresholds) if i > 0 else min(self.labels) # 閾値の上限 (ラベルの最大値) high = max(self.labels) # 閾値の候補を追加する t = trial.suggest_uniform(f't{i}', low, high) thresholds.append(t) # 閾値の候補を元に QWK を計算する opt_y_pred = self.adjust(self.y_pred, thresholds) return qwk(self.y_true, opt_y_pred) def adjust(self, y_pred, thresholds): """閾値にもとづいて予測を補正するメソッド""" opt_y_pred = pd.cut(y_pred, [-np.inf] + thresholds + [np.inf], labels=self.labels) return opt_y_pred def lgb_custom_metric_qwk_regression(preds, data): """LightGBM のカスタムメトリックを計算する関数 回帰問題として解いた予測から QWK を計算する""" # 正解ラベル y_true = data.get_label() # 予測ラベル y_pred = np.clip(preds, 0, 7).round() # 単純に clip して四捨五入する # QWK を計算する return 'qwk', qwk(y_true, y_pred), True def qwk(y_true, y_pred): """QWK (Quadratic Weighted Kappa) を計算する関数""" return cohen_kappa_score(y_true, y_pred, weights='quadratic') class ModelExtractionCallback(object): """lightgbm.cv() から学習済みモデルを取り出すためのコールバックに使うクラス""" def __init__(self): self._model = None def __call__(self, env): # _CVBooster の参照を保持する self._model = env.model def _assert_called_cb(self): if self._model is None: # コールバックが呼ばれていないときは例外にする raise RuntimeError('callback has not called yet') @property def boosters_proxy(self): self._assert_called_cb() # Booster へのプロキシオブジェクトを返す return self._model @property def raw_boosters(self): self._assert_called_cb() # Booster のリストを返す return self._model.boosters @property def best_iteration(self): self._assert_called_cb() # Early stop したときの boosting round を返す return self._model.best_iteration def _numerical_only_df(df_): """数値のカラムだけ取り出す関数""" number_cols = [name for name, dtype in df_.dtypes.items() if issubclass(dtype.type, numbers.Number)] numerical_df = df_[number_cols] return numerical_df def main(): # データセットを読み込む train_df = pd.read_csv('train.csv.zip') # 説明変数 x_train = train_df.drop('Response', axis=1) # 簡単にするため、とりあえずすぐに使えるものだけ取り出して使う x_train = _numerical_only_df(x_train) # Id 列をインデックスに設定する x_train = x_train.set_index('Id') # 目的変数 y_train = train_df.Response.astype(float) # 目的変数はゼロからスタートにする y_train -= 1 # LightGBM のデータセット表現にする lgb_train = lgb.Dataset(x_train, y_train) # 回帰問題として解く lgbm_params = { 'objective': 'regression', 'metric': 'rmse', 'first_metric_only': True, 'verbose': -1, } # データ分割の方法 folds = KFold(n_splits=5, shuffle=True, random_state=42) # 学習済みモデルを取り出すためのコールバックを用意する extraction_cb = ModelExtractionCallback() callbacks = [ extraction_cb, ] # 5-Fold CV で検証する result = lgb.cv(lgbm_params, lgb_train, num_boost_round=10000, early_stopping_rounds=100, folds=folds, verbose_eval=10, feval=lgb_custom_metric_qwk_regression, callbacks=callbacks, seed=42, ) # early stop したラウンドでの QWK を出力する print(f'CV Mean QWK: {result["qwk-mean"][-1]}') # 学習データを学習済みモデルを使って OOF で予測する boosters = extraction_cb.raw_boosters best_iteration = extraction_cb.best_iteration oof_y_pred = np.zeros_like(y_train) for (_, val_index), booster in zip(folds.split(x_train, y_train), boosters): x_val = x_train.iloc[val_index] oof_y_pred[val_index] = booster.predict(list(x_val.to_numpy()), num_iteration=best_iteration) # ひとまず clipping + 四捨五入する oof_y_pred = np.clip(oof_y_pred, 0, 7).round() # 最適化していない状態での OOF Prediction の QWK raw_oof_qwk = qwk(y_train, oof_y_pred) print(f'Raw OOF QWK: {raw_oof_qwk}') # Optuna を使って QWK の閾値を最適化する objective = OptunaRounder(y_train, oof_y_pred) # 探索に 100 秒かける study = optuna.create_study(direction='maximize') study.optimize(objective, timeout=100) # 見つけた閾値 best_thresholds = sorted(study.best_params.values()) print(f'Optimized thresholds: {best_thresholds}') # 閾値を最適化した場合の QWK optimized_oof_y_pred = objective.adjust(oof_y_pred, best_thresholds) optimized_oof_qwk = qwk(y_train, optimized_oof_y_pred) print(f'Optimized OOF QWK: {optimized_oof_qwk}') # テスト用データを読み込む test_df = pd.read_csv('test.csv.zip') test_df = test_df.set_index('Id') # 数値のカラムだけ取り出す numerical_test_df = _numerical_only_df(test_df) # 学習済みモデルの Averaging で予測する boosters_proxy = extraction_cb.boosters_proxy test_y_preds = boosters_proxy.predict(numerical_test_df, num_iteration=best_iteration) test_y_pred_avg = np.mean(test_y_preds, axis=0) # サンプル提出ファイルを読み込む submission_df = pd.read_csv('sample_submission.csv.zip') submission_df = submission_df.set_index('Id') # 最適化していない予測値 (単純なクリッピングと四捨五入) submission_df.Response = np.clip(test_y_pred_avg, y_train.min(), y_train.max() + 1).round().astype(int) + 1 submission_df.to_csv('unoptimized_submission.csv') # 最適化した予測値 optimized_test_y_pred_avg = objective.adjust(test_y_pred_avg, best_thresholds) submission_df.Response = optimized_test_y_pred_avg.astype(int) + 1 submission_df.to_csv('optimized_submission.csv') if __name__ == '__main__': main()
上記を実行した結果は次のとおり。
$ python qwkregopt.py [10] cv_agg's rmse: 2.01959 + 0.0101767 cv_agg's qwk: 0.408957 + 0.00513284 [20] cv_agg's rmse: 1.9165 + 0.008155 cv_agg's qwk: 0.507252 + 0.00192446 [30] cv_agg's rmse: 1.87492 + 0.00797763 cv_agg's qwk: 0.553388 + 0.00129979 ...(省略)... [200] cv_agg's rmse: 1.8374 + 0.00729844 cv_agg's qwk: 0.60726 + 0.00206979 [210] cv_agg's rmse: 1.83764 + 0.00757892 cv_agg's qwk: 0.607434 + 0.00210639 [220] cv_agg's rmse: 1.83771 + 0.00759351 cv_agg's qwk: 0.607882 + 0.00205622 CV Mean QWK: 0.6054069736940326 Raw OOF QWK: 0.605407593931218 [I 2020-03-01 12:39:28,349] Finished trial#0 resulted in value: 0.14760909846772208. Current best value is 0.14760909846772208 with parameters: {'t0': 0.6069999513322623, 't1': 6.071871351250824, 't2': 6.421127597260041, 't3': 6.6499083619322334, 't4': 6.825786239070231, 't5': 6.94896643994053, 't6': 6.988042174746154}. [I 2020-03-01 12:39:28,461] Finished trial#1 resulted in value: 0.3273196067078864. Current best value is 0.3273196067078864 with parameters: {'t0': 1.3022593301640533, 't1': 3.3523582269835885, 't2': 5.775809715879871, 't3': 6.748066158840195, 't4': 6.813600862422318, 't5': 6.952492269479133, 't6': 6.986909664863027}. [I 2020-03-01 12:39:28,580] Finished trial#2 resulted in value: 0.1845825834695184. Current best value is 0.3273196067078864 with parameters: {'t0': 1.3022593301640533, 't1': 3.3523582269835885, 't2': 5.775809715879871, 't3': 6.748066158840195, 't4': 6.813600862422318, 't5': 6.952492269479133, 't6': 6.986909664863027}. ...(省略)... [I 2020-03-01 12:41:08,371] Finished trial#660 resulted in value: 0.3366380459371958. Current best value is 0.6461578531751463 with parameters: {'t0': 1.9366086554521658, 't1': 2.765496885821397, 't2': 3.711341767628701, 't3': 3.8616368179727982, 't4': 4.309923993071266, 't5': 5.4467521001848525, 't6': 5.753790447645578}. Optimized thresholds: [1.9366086554521658, 2.765496885821397, 3.711341767628701, 3.8616368179727982, 4.309923993071266, 5.4467521001848525, 5.753790447645578] Optimized OOF QWK: 0.6461578531751463
最適化する前の OOF Prediction の QWK が約 0.6054 なのに対して、閾値を最適化すると約 0.6461 と、大きくスコアが向上している。
手元の交差検証だけでは何か間違っている恐れがあるので、テストデータに対するサブミッションまでやってみた。
先ほどのモジュールを実行すると unoptimized_submission.csv という最適化前のファイルと、optimized_submission.csv という最適化後のファイルができる。
すると、最適化する前の Public LB / Private LB が 0.61931 / 0.62190 だったのに対して、最適化した後は 0.65801 / 0.66005
と、ちゃんと効果があることが確認できた。
所感など
QWK の最適化は、真のラベルの分布に依存している。 そのため、学習済みモデルを使って予測する対象が学習時の分布と異なる場合、最適化の効果が得られない。 この点から、Private データの分布が異なる場合には最適化が有効ではないと考えられる。
だとすると、Private データの分布が学習時と同じであるという確証が得られない限りは使うのは賭けになるのでは、という疑問が浮かんだ。 とはいえ、うまくいった場合の効果は絶大なので、ある程度の仮定のもとに使う、という選択肢は現実的なのかもしれない。
いじょう。
参考
