Ccmmutty logo
Commutty IT
・4 min read

ロバスト学習についてその1 Huber損失最小化学習

機械学習Python
https://picsum.photos/seed/a5b61d6a9dad42bfa9723d8f9e3a9111/600/800

ロバスト学習その1 Huber損失最小化学習

参考文献

文献1.のChapter6を参考としている。

Huber損失最小化学習

概要

ロバスト学習を行う方法の1つであるHuber損失最小化学習について紹介する。 その方法は、以下に定義するHuber損失を最小化することである。 これは、いわゆるl1l^1損失とl2l^2損失の両方を加味した結果提案されたものとなっている。
ρHuber(r,η)=r22I(∣r∣≤η)+(η∣r∣−η22)I(∣r∣>η)\rho_{Huber}(r, \eta) = \frac{r^2}{2}\mathbb{I}(|r| \le \eta) + \left( \eta|r| - \frac{\eta^2}{2} \right)\mathbb{I}(|r| > \eta)
このHuber損失と、学習に用いた訓練標本の数nnと各訓練標本の目的変数とその予測値の差rir_iを用いて目的関数J(θ)J(\theta)を以下のように構成する。
J(θ)=∑i=1nρHuber(ri,η)J(\theta) = \sum^{n}_{i = 1} \rho_{Huber}(r_i, \eta)
この目的関数をそのまま最小化するように推定するパラメータを更新していくというわけではなく、今回は代理損失を用いる。
代理損失とは、Huber損失の上界をとるような別の最適化しやすい関数を指し、本ドキュメントでは以下で紹介するパラメータWWを用いて表現される。

Huber損失最小化学習のための記号とパラメータ

数式中の記号コード中の変数概要
yyy目的変数で、実数値をとる。
yiy_iy[i-1]ii番目の標本が持つ目的変数の値。
Φ\PhiX適当な関数Φ\Phiを説明変数XXに適用して作ったデザイン行列。確率変数ではない。
ϕ(xj)\phi(x_j)X[:, j-1]デザイン行列からjj列目をとったもの。
θ\thetacoef求めたい回帰係数。推定値だけどハットは省略した。
fθ(xi)f_{\theta}(x_i)yhat訓練標本iiに対する推測値。
rir_iresid[i-1]yi−fθ(xi)y_i - f_{\theta}(x_i)で、推測値と実値の誤差。
η\etaeta学習の際に用いる閾値。値は適当にセットした。
wiw_iw_diag[i-1]表が崩れてしまったので別途下に記載
nnlen(y)訓練標本の数となっている。
WWW対角行列diag(w1,…,wn)diag(w_1,\dots,w_n)を指す。
†\dagger-行列AAに対してA†A^{\dagger}はAAの一般逆行列を指す。
wi=w_i = I(∣r∣≤η)+η∣ri∣I(∣r∣>η)\mathbb{I}(|r| \le \eta) + \frac{\eta}{|r_i|} \mathbb{I}(|r| > \eta)

Huber損失最小化学習のための疑似コード

  1. θ\thetaの初期値を(Φ⊤Φ)†Φ⊤y\left( \Phi^{\top}\Phi \right)^{\dagger}\Phi^{\top}yとする。つまり、普通に線形回帰係数を得る。
  2. 推定値θ\thetaの値が収束するまで、以下の2つのステップを繰り替えす。
    1. 現在の推定値θ\thetaを使って、WWを計算する。
    2. θ\thetaを(Φ⊤WΦ)†Φ⊤Wy\left( \Phi^{\top}W\Phi \right)^{\dagger}\Phi^{\top}Wyで更新する。

データセットirisの呼び出しと、線形回帰係数を得るまで。

python
# インポート
import numpy as np
import statsmodels.api as sm
# データirisをインポートする
iris = sm.datasets.get_rdataset("iris", "datasets").data
# 説明変数の行列Xと被説明変数のベクトルyを作る
# sm.add_constant(X)と同じことをするためにiris["const"]を追加
iris["const"] = 1
X = iris[["const", "Sepal.Length", "Petal.Length"]].values
y = iris["Sepal.Width"].values
# 演算子@で内積をとり、np.linalg.solveで引数1の逆行列を引数2にかけたものを計算できる。
# 下の書き方は速度面に配慮した記法で、次の記事を参考にした。
# https://qiita.com/fujiisoup/items/e7f703fc57e2dfc441ad
coef = np.linalg.solve(X.T @ X, X.T @ y)
coef
/usr/local/lib/python3.7/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead. import pandas.util.testing as tm
array([ 1.03806906, 0.56118597, -0.33526674])

実際にcoefを更新する。

python
def makew(ndarr, eta):
  w_diag = np.where(ndarr <= eta, ndarr, eta) / ndarr
  W = np.diag(w_diag)
  return W
python
def make_abs_resid(y, yhat):
  resid = np.abs(y - yhat)
  return resid
python
yhat = X @ coef
resid = make_abs_resid(y, yhat)
eta = 0.5
W = makew(resid, eta)
python
# 疑似コードに収束するまで、と書いてあるが今回は収束を確認するための閾値を設けずに1万回だけ更新する。
python
for _ in range(10000):
    yhat = X @ coef
    resid = make_abs_resid(y, yhat)
    W = makew(resid, eta)
    coef = np.linalg.inv(X.T @ W @ X) @ X.T @ W @ y
python
coef
array([ 1.10332467, 0.54918993, -0.33358833])

Discussion

コメントにはログインが必要です。