Huber損失と二乗誤差の違いは？

Huber損失は誤差が小さい範囲では二乗誤差、大きい範囲では絶対誤差として振る舞います。これにより外れ値の影響を抑えつつ、小さな誤差では滑らかな最適化が行えます。

RANSAC や Theil-Sen との使い分けは？

外れ値が10%以下なら Huber 回帰で十分です。明確な外れ値が大量にある場合は RANSAC で候補を絞るか、中央値ベースの Theil-Sen 推定が頑健です。

ロバスト回帰でも前処理は必要？

外れ値耐性は上がりますが、特徴量のスケーリング（標準化）は依然重要です。スケールが大きく異なる特徴量があると Huber の閾値が効きにくくなります。

ロバスト回帰

中級

02.1.3

ロバスト回帰

最終更新 2020-02-26 読了時間 4 分

まとめ

ロバスト回帰は、残差の 2 乗ではなく Huber 損失のように「大きな残差ほど罰の増え方を抑える」損失で直線を当てはめ、外れ値 1 点に直線が引きずられるのを防ぐ手法です。Huber 損失はしきい値 $\delta$ より小さい残差には二乗誤差、大きい残差には線形の誤差を使うので、普通の点は最小二乗法と同じに扱い、外れ値だけ影響を頭打ちにします。scikit-learn では HuberRegressor の epsilon がこのしきい値、alpha が L2 正則化の強さです。このページでは 30 点のうち 1 点だけを 500 に置き換えたデータで、最小二乗法・リッジ回帰・Huber 回帰の 3 本の直線を重ね、Huber だけが残り 29 点の傾きを保つ様子を確かめます。

先に線形回帰で、残差を 2 乗するため残差 10 の点が残差 1 の点の 100 倍の重みを持つことを確認しておくと、このページで損失の形を変える理由が分かります。

直感 #

センサーの読み取りミスで 30 点のうち 1 点だけ $y = 500$ になったとします。最小二乗法はその 1 点の残差 400 を 2 乗して 160,000 と数えるので、他の 29 点の残差（それぞれ 10 以下）を全部足しても勝てず、直線は外れ値の方へ傾きます。Huber 損失は残差が $\delta$ を超えたら 2 乗をやめて 1 次で増やすので、外れ値の罰は 400 の定数倍で止まり、29 点の多数派が直線の位置を決めます。残差が小さい点では二乗損失と同じなので、外れ値のないデータでは最小二乗法とほぼ同じ直線が得られます。

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：誤差の大きい点の重みを下げるだけで、直線は多数派に戻る。HuberRegressor の epsilon は「残差がどこから大きいと見なすか」のしきい値で、下のコードでは 3.0 にしています。

式で書くとどこで損失の形が切り替わるか #

しきい値 $\delta$ を境に、内側は放物線、外側は直線になる 1 つの関数が Huber 損失です。残差 $r = y - \hat{y}$ とし、しきい値 $\delta > 0$ を用いると Huber損失は

$$ \ell_\delta(r) = \begin{cases} \dfrac{1}{2} r^2, & |r| \le \delta, \\\\ \delta \left(|r| - \dfrac{1}{2}\delta \right), & |r| > \delta. \end{cases} $$

$|r| \le \delta$ の内側では通常の二乗誤差と同じ挙動、外側では $|r|$ に比例して増えるため、残差 400 の外れ値でも損失は $\delta \times 400$ 程度で止まります。導関数（影響関数）は内側で $r$、外側で $\pm\delta$ の折れ線になり、1 点が直線を引っぱる力が $\delta$ で頭打ちになります。2 つ目の式の $-\frac{1}{2}\delta$ は、$|r| = \delta$ で 2 本の式がなめらかにつながるための補正項です。

Pythonを用いた実験や説明 #

まず Huber 損失を二乗損失・絶対値損失と重ねて形を確かめ、次に 30 点のうち 1 点を 500 に置き換えたデータで HuberRegressor・Ridge・LinearRegression の直線を比べます。最初のブロックは描画と数値計算のライブラリを読み込むだけです。

1
2
import matplotlib.pyplot as plt
import numpy as np

Huber損失と関連損失の比較 #

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
def huber_loss(r: np.ndarray, delta: float = 1.5):
    half_sq = 0.5 * np.square(r)
    lin = delta * (np.abs(r) - 0.5 * delta)
    return np.where(np.abs(r) <= delta, half_sq, lin)

delta = 1.5
r_vals = np.arange(-2, 2, 0.01)
h_vals = huber_loss(r_vals, delta=delta)

plt.figure(figsize=(8, 6))
plt.plot(r_vals, np.square(r_vals), "red", label=r"二乗損失 $r^2$")
plt.plot(r_vals, np.abs(r_vals), "orange", label=r"絶対値損失 $|r|$")
plt.plot(r_vals, h_vals, "green", label=fr"Huber損失 ($\delta={delta}$)")
plt.axhline(0, color="k", linewidth=0.8)
plt.grid(True, alpha=0.3)
plt.legend()
plt.xlabel("残差 $r$")
plt.ylabel("損失値")
plt.title("二乗・絶対値・Huber損失の比較")
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-4 行目数式をそのまま 2 本の式と条件分岐で書く
half_sq が $|r| \le \delta$ 側の $\frac{1}{2}r^2$、lin が外側の $\delta(|r| - \frac{1}{2}\delta)$ で、np.where がしきい値で切り替えています。両方を全要素で計算してから選ぶ書き方なので、r が配列でも 1 行で済みます。lin の - 0.5 * delta を消すと $|r| = \delta$ で曲線に段差ができるので、実際に消して描いてみると補正項の役割が分かります。
6-8 行目δ = 1.5 と残差の範囲 −2〜2
しきい値を 1.5 にしているのは、描画範囲 $[-2, 2]$ の中に放物線から直線への切り替わりが入るようにするためです。delta を 3.0 にすると範囲内では二乗損失と完全に重なり、0.5 にするとほぼ絶対値損失に見えます。np.arange(-2, 2, 0.01) は 400 点の細かい格子で、折れ目をなめらかに描くための刻みです。
11-13 行目3 つの損失を同じ軸に重ねる
赤が二乗損失 $r^2$、橙が絶対値損失 $|r|$、緑が Huber 損失です。原点付近では緑が赤に沿い、$|r| > 1.5$ で橙と同じ傾きの直線に変わるのが見どころで、この「途中で乗り換える」形が外れ値の罰を頭打ちにします。$r = 2$ で赤は 4、緑は 2.6 程度と差がつき始めます。

外れ値を含むデータでの挙動 #

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
np.random.seed(42)

N = 30
x1 = np.arange(N)
x2 = np.arange(N)
X = np.c_[x1, x2]
epsilon = np.random.rand(N)
y = 5 * x1 + 10 * x2 + epsilon * 10

y[5] = 500  # 外れ値を1点だけ挿入

plt.figure(figsize=(8, 6))
plt.plot(x1, y, "ko", label="data")
plt.xlabel("$x_1$")
plt.ylabel("$y$")
plt.legend()
plt.title("外れ値を含むデータ")
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

3-8 行目傾き 15 の直線にノイズを足した 30 点
x1 と x2 は同じ 0〜29 の値なので、$y = 5x_1 + 10x_2$ は実質 $y = 15x_1$ の 1 変数問題です。2 列は完全に相関していて係数の分け方は一意に決まりませんが、図に描くのは予測値 predict(X) なので直線は問題なく引けます。ノイズは rand の一様乱数を 10 倍した 0〜10 の範囲で、外れ値と比べればごく小さい幅です。
10 行目1 点だけを 500 に書き換える
インデックス 5（本来の値は 85 前後）を 500 に置き換え、残差 400 以上の外れ値を 1 つだけ作ります。この 1 行が次のブロックで 3 本の直線に差を生む原因で、消せば 3 本はほぼ重なります。値を 150 程度に下げると最小二乗法の傾きのずれも小さくなり、外れ値の大きさと直線の歪みの関係を試せます。
12-13 行目横軸を x1 だけにして 2 次元で描く
説明変数は 2 列ありますが同じ値なので、x1 を横軸に取れば 1 枚の散布図で全体が見えます。黒丸 30 個のうち 1 個だけが $y = 500$ に飛び出しているのを、この時点で目視しておくと次の図の読み方が決まります。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
from sklearn.linear_model import HuberRegressor, Ridge, LinearRegression

plt.figure(figsize=(8, 6))

huber = HuberRegressor(alpha=0.0, epsilon=3.0)
huber.fit(X, y)
plt.plot(x1, huber.predict(X), "green", label="Huber回帰")

ridge = Ridge(alpha=1.0, random_state=0)
ridge.fit(X, y)
plt.plot(x1, ridge.predict(X), "orange", label="リッジ回帰 (α=1.0)")

ols = LinearRegression()
ols.fit(X, y)
plt.plot(x1, ols.predict(X), "r-", label="最小二乗法 (OLS)")

plt.plot(x1, y, "kx", alpha=0.7)
plt.xlabel("$x_1$")
plt.ylabel("$y$")
plt.legend()
plt.title("外れ値を含むときの推定曲線の違い")
plt.grid(alpha=0.3)
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

5-7 行目epsilon=3.0 がしきい値、alpha=0.0 で正則化なし
HuberRegressor の epsilon は本文の $\delta$ に相当し、残差をスケール推定値で割った値がこれを超えると線形の罰に切り替わります（1.0 以上が必須、既定は 1.35）。alpha=0.0 で L2 正則化を切っているので、違いは損失の形だけです。epsilon を 100 のように大きくすると全点が二乗損失の範囲に入り、緑の直線は赤の最小二乗法に重なります。
9-11 行目リッジ回帰は外れ値対策にならないことを示す
Ridge(alpha=1.0) は係数の大きさに罰則をかけるだけで、損失は二乗誤差のままです。外れ値の残差 400 は 2 乗して数えられるので、橙の直線は赤の最小二乗法とほぼ同じ位置に来ます。「正則化」と「ロバスト化」は別の問題への対処だと、この 1 本を並べることで確認できます。
13-15 行目比較の基準になる最小二乗法
LinearRegression は残差二乗和を最小にするので、1 点の残差 400 を減らすために直線全体を外れ値側へ持ち上げます。赤の直線が黒い×印の多数派から上にずれ、傾きも真の 15 から変わっているのが外れ値の影響そのものです。
17 行目観測点を最後に薄く重ねる
3 本の直線を描いた後に "kx"（黒い×）で 30 点を重ね、alpha=0.7 で直線を隠さないようにしています。緑の Huber 回帰が×印の列の上を通り、赤と橙がその上に浮いているかを見るための配置です。

実行結果の読み方 #

1 枚目の図では、緑の Huber 損失が原点付近で赤の二乗損失に重なり、$|r| > 1.5$ から橙の絶対値損失と同じ傾きに変わります。3 枚目の図では、緑（Huber 回帰）が 29 点の×印の列をなぞる位置に来るのに対し、赤（最小二乗法）と橙（リッジ回帰）は $y = 500$ の 1 点に引かれて上にずれ、2 本はほぼ重なります。リッジ回帰が最小二乗法と並ぶのは、L2 正則化が係数の大きさを抑えるだけで損失の形を変えないためです。epsilon を 1.35（既定値）まで下げると Huber の直線はさらに多数派に寄り、100 まで上げると赤に一致します。

使いどころと落とし穴 #

センサー値や手入力データのように、全体の 1〜10% に異常値が混ざるデータでは HuberRegressor を最小二乗法の代わりに最初から使う価値があります。外れ値が 10% を超えるなら、分位点回帰の中央値回帰や RANSACRegressor のように外れ値を明示的に除く手法を検討します。落とし穴は 2 つあります。1 つ目はスケールで、epsilon は残差をスケール推定値で割った値に対するしきい値なので、特徴量の単位がばらばらだと意図したしきい値になりません。StandardScaler を前に置いてください。2 つ目は「ロバスト＝万能」という誤解で、Huber 損失は $y$ 方向の外れ値には強いものの、$x$ 方向に離れた点（てこ点）には最小二乗法と同じく引かれます。

外れ値の影響 #

外れ値の割合を変えた場合、ロバスト回帰と OLS がどう異なるか比較できます。

参考文献 #

Huber, P. J. (1964). Robust Estimation of a Location Parameter. The Annals of Mathematical Statistics, 35(1), 73–101.
Hampel, F. R. et al. (1986). Robust Statistics: The Approach Based on Influence Functions. Wiley.
Huber, P. J., & Ronchetti, E. M. (2009). Robust Statistics (2nd ed.). Wiley.

線形回帰 — 標準的な最小二乗
分位点回帰 — もう一つの頑健な手法
中央値絶対誤差 — ロバストな評価指標
Huber Loss — MAEとMSEを統合したロバスト損失関数