13か月前公開・13か月前更新・4 min read

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

機械学習Python

https://picsum.photos/seed/a5b61d6a9dad42bfa9723d8f9e3a9111/600/800

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

参考文献

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

https://bookclub.kodansha.co.jp/product?item=0000148211

Huber損失最小化学習

概要

ロバスト学習を行う方法の1つであるHuber損失最小化学習について紹介する。その方法は、以下に定義するHuber損失を最小化することである。これは、いわゆる

l^1

損失と

l^2

損失の両方を加味した結果提案されたものとなっている。

\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損失と、学習に用いた訓練標本の数

n

と各訓練標本の目的変数とその予測値の差

r_i

を用いて目的関数

J(\theta)

を以下のように構成する。

J(\theta) = \sum^{n}_{i = 1} \rho_{Huber}(r_i, \eta)

この目的関数をそのまま最小化するように推定するパラメータを更新していくというわけではなく、今回は代理損失を用いる。

代理損失とは、Huber損失の上界をとるような別の最適化しやすい関数を指し、本ドキュメントでは以下で紹介するパラメータ

W

を用いて表現される。

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

数式中の記号	コード中の変数	概要
$y$	`y`	目的変数で、実数値をとる。
$y_i$	`y[i-1]`	$i$ 番目の標本が持つ目的変数の値。
$\Phi$	`X`	適当な関数 $\Phi$ を説明変数 $X$ に適用して作ったデザイン行列。確率変数ではない。
$\phi(x_j)$	`X[:, j-1]`	デザイン行列から $j$ 列目をとったもの。
$\theta$	`coef`	求めたい回帰係数。推定値だけどハットは省略した。
$f_{\theta}(x_i)$	`yhat`	訓練標本 $i$ に対する推測値。
$r_i$	`resid[i-1]`	$y_i - f_{\theta}(x_i)$ で、推測値と実値の誤差。
$\eta$	`eta`	学習の際に用いる閾値。値は適当にセットした。
$w_i$	`w_diag[i-1]`	表が崩れてしまったので別途下に記載
$n$	`len(y)`	訓練標本の数となっている。
$W$	`W`	対角行列 $diag(w_1,\dots,w_n)$ を指す。
$\dagger$	`-`	行列 $A$ に対して $A^{\dagger}$ は $A$ の一般逆行列を指す。

w_i =

\mathbb{I}(|r| \le \eta) + \frac{\eta}{|r_i|} \mathbb{I}(|r| > \eta)

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

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

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

# インポート
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`を更新する。

def makew(ndarr, eta):
  w_diag = np.where(ndarr <= eta, ndarr, eta) / ndarr
  W = np.diag(w_diag)
  return W

def make_abs_resid(y, yhat):
  resid = np.abs(y - yhat)
  return resid

yhat = X @ coef
resid = make_abs_resid(y, yhat)
eta = 0.5
W = makew(resid, eta)

# 疑似コードに収束するまで、と書いてあるが今回は収束を確認するための閾値を設けずに1万回だけ更新する。

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