🤥

Sparse Estimation with Python 勉強ノート 1 疎な線形回帰

2024/03/06に公開

鈴木讓著
スパース推定100問 with Python
2021-01-29
https://www.kyoritsu-pub.co.jp/book/b10003298.html

Christopher M. Bishop 著, 元田浩栗田多喜夫樋口知之松本裕治村田昇監訳
パターン認識と機械学習ベイズ理論による統計的予測 (Pattern Recognition and Machine Learning)
2012-01-20
https://www.maruzen-publishing.co.jp/item/?book_no=294524

研究でスパース行列分解や変数選択を扱うことがあったので、興味本位でスパース推定について勉強しようと思いました。
RedgeやLassoのベイズ的な観点からの導出は知っているものの、真っ当な方向からアプローチする方法も知っておきたいです。この記事ではこの本(Sparse Estimation with Python; SEwP)の内容を勉強しつつ不親切な部分を補完したメモと、加えてベイズ的な方法による結果について書きます。
式などについては怪しい部分が多いので、データを基底関数で置き換えた場合(PRMLのやりかた)を確認した方を書きました。

なんで「正則化項」が「正則化」と呼ばれるのか初めて知りました。線形回帰っておもしろいですね。

この記事でわかること

(1) 線形回帰について

回帰関数 $y(\bold{x},\bold{w}) = \bold{x}^\top \bold{w}$ の単純な線形回帰は、「ターゲット(目的変数) $t$ が、尤度(データ $\bold{X}$ を観測したときのパラメータ $\bold{w}$ が、ターゲットにどれだけ従っているか)である $p(\bold{t}|\bold{X},\bold{w})$ 、具体的には $\prod_{n=1}^N \mathcal{N}(t_n|y(\bold{x}_n,\bold{w}), β^{-1})$ を最大化する」という 最尤推定 の考えから導出される。
最尤推定により尤度の最大値を求めると、パラメータ $\bold{w}$ が関わる項が $-\frac{β}{2}||\bold{t} - y(\bold{X},\bold{w})||^2$ として現れ、これを最大化する、つまり 二乗和誤差 を最小化するという問題に帰着する。
二乗和誤差の最小値を求める条件を計算することで、 $\bold{w}$ の解が $\bold{X}$ の擬似逆行列を含む $\bold{w}_{\text{ML}} = (\bold{X}^\top \bold{X})^{-1} \bold{X}^\top \bold{t}$ であることが導かれ、これが線形回帰のパラメータの最尤解である。
しかし、最尤推定では擬似逆行列の逆行列部分 $(\bold{X}^\top \bold{X})^{-1}$ が正則でない場合(例えばデータ $\bold{X}$ が線形従属に近い、つまり相関の高い特徴次元の列を持つ時など)に計算が不安定になり、パラメータ $\bold{w}$ が大きな値になってしまう。
パラメータが大きいということは回帰関数 $y(\bold{x},\bold{w})$ が極端な形になっているということであり、これは 過学習 になっている。

(2) 線形回帰を安定させる正則化について

尤度 $p(\bold{t}|\bold{X},\bold{w}, β)$ だけでなく、パラメータ $\bold{w}$ が $\mathcal{N}(\bold{w}|\bold{0}, α^{-1}\bold{I})$ 、つまり0付近にあるだろうという 事前分布 $p(\bold{w}|α)$ を考えることで、ベイズの定理から 事後分布 $p(\bold{w}|\bold{X},\bold{w}, α, β)$ を導出でき、これを最大化する MAP推定 が可能になる。
事後分布の最大化条件を計算することで、尤度部分の二乗誤差 $-\frac{β}{2}||\bold{t} - y(\bold{X},\bold{w})||^2$ に加えて、事後分布の部分が $-\frac{α}{2}||\bold{w}||^2$ という形で出てくるため、これを合わせて $λ=α/β$ とした正則化二乗和誤差 $\frac{1}{2}||\bold{t} - y(\bold{X},\bold{w})||^2 + \frac{λ}{2}||\bold{w}||^2$ の最小化問題に帰着する。
正則化二乗和誤差の最小化条件をパラメータについて解くと $\bold{w}_{\text{MAP}} = (\bold{X}^\top \bold{X} - λ\bold{I})^{-1} \bold{X}^\top \bold{t}$ が得られ、これが正則化線形回帰のパラメータである。
逆行列の中に単位行列 $I$ が入っているため、 $\bold{X}^\top \bold{X}$ が非正則に近くても逆行列の中全体は正則になり、極端に大きなパラメータになりにくくなる。
この結果は事前分布によりパラメータ $\bold{w}$ が極端な値を取りにくいという情報が織り込まれている所に起因しており、正則化項であるL2ノルム $||\bold{w}||^2$ の ロバスト性 を使った線形回帰を Redge回帰 という。

(3) スパース推定について

事前分布 $p(\bold{w}|α)$ にラプラス分布(平均付近に値が密集していて、裾野が厚い分布) $\text{Laplace}(\bold{w}|\bold{0},β^{2},α^{-1})$ を仮定すると、パラメータが0に集中するという信念が反映された事後分布が作れる。
このMAP推定の計算によって、L2ノルム正則化項の代わりにL1ノルムの項 $\frac{α}{β} ||\bold{w}||_{\text{L}1}$ が導かれ、 $λ=2αβ$ と置いたとき誤差関数 $\frac{1}{2}||\bold{t} - y(\bold{X},\bold{w})||^2 + λ||\bold{w}||_{\text{L}1}$ が導出される。
この誤差関数は非滑らかな凸関数であるL1ノルムを含んでいるので、最小化条件は普通の微分ではなく 劣微分 を使う。劣微分は凸関数の微分不可能な点を境界に3つの場合分けが発生し、微分不可能点の劣勾配に内点の閉区間を構成する。大雑把に言えば劣微分係数が一定に固定されるような内点の値の条件が発生する。
誤差関数の最小値の条件をパラメータ $w$ について解くと、劣微分の区間は各特徴次元について $φ_d = \sum_{n=1}^N (t_n - \sum_{d'≠d}^D x_{n,d'}w_{d'}) x_{n,d}$ と正則化係数の絶対値 $|λ|$ の大きさの関係によって分岐することが計算で確認できて、特に $φ_d \leq |λ|$ の条件が満たされる場合、その特徴次元のパラメータは0になる、つまり $\bold{w}_{\text{Lasso}}$ が スパース な解になる。
このスパース線形回帰を Lasso回帰 という。
Lasso回帰では、回帰に重要でない特徴次元が0になるため、回帰問題を解くことが自動的に変数選択になる。

正規分布(平均 $\bold{μ}$ 、共分散行列 $\bold{Σ}$ )、特に平均が $\bold{0}$ で共分散が全次元で等しい $σ\bold{I}$ 時

\begin{align*} \mathcal{N}(\bold{z} | \bold{μ}, \bold{Σ}) & = \frac{1}{(2π)^{\frac{D}{2}}} \frac{1}{|\bold{Σ}|^{\frac{1}{2}}} e^{-\frac{1}{2} (\bold{z}-\bold{μ})^\top \bold{Σ}^{-1} (\bold{z}-\bold{μ})} \\ \mathcal{N}(\bold{z} | \bold{0}, σ^{-1}\bold{I}) & =\frac{σ}{(2π)^{\frac{D}{2}}} e^{-\frac{σ}{2} ||\bold{z}||^2} \end{align*}

ベイズの定理

p(\bold{w}|\bold{X},\bold{w}, α, β) ∝ p(\bold{t}|\bold{X},\bold{w}, β) p(\bold{w}|α)

ノルム

\begin{align*} \text{L}_2 & : ||\bold{z}|| = (\sum_{d=1}^D z_d^2)^\frac{1}{2} \\ \text{L}_1 & : ||\bold{z}||_{\text{L}1} = \sum_{d=1}^D |z_d^2| \end{align*}

次のライブラリを使う。

import numpy as np
np.set_printoptions(precision=6, suppress=True)
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt
import arviz
arviz.style.use("arviz-darkgrid")

線形回帰

次元 $D \geq 1$ の特徴量ベクトルをサンプル数 $N \geq 1$ 個集めるとき、通常は $N>D$ を仮定して計算を考える。

\begin{align*} \bold{X} & = [ \bold{x}_n \in \mathbb{R}^{D} ] \in \mathbb{R}^{N×D} \\ & = \begin{bmatrix} x_{1,1} & \dots & x_{N,1} \\ \vdots & \ddots & \vdots \\ x_{D,1} & \dots & x_{N,D} \\ \end{bmatrix} \\ \bold{t} & = [ t_n ] \in \mathbb{R}^{D} \end{align*}

しかし $N<D$ の場合最小二乗の解が求まらない、情報量基準を使った変数選択が計算困難になるなどよくないことが起こる。
このような状況を改善するにはRidgeによる正則化線形回帰やLassoを使った変数選択について学ぶ。

1. 線形回帰

PRML 3.1.1

$\bold{t}$ を関数 $y(\bold{x}, \bold{w}) = w_0 + \bold{w}^\top \bold{x}$ で回帰する。
確率的な考え方では、ターゲットデータ $t$ は決定論的な関数 $y$ へ、平均0分散 $β^{-1}$ のガウスノイズを加えたものだと考えられる。

\begin{align*} t & = y(\bold{x}, \bold{w}) + ε \\ & \sim p(t|\bold{x}, \bold{w}, β) \\ & = \mathcal{N}(t|y(\bold{x}, \bold{w}), β^{-1}) \end{align*}

このとき、条件付き期待値(つまりtが従う分布の平均の点)は次のように簡単に表せる。

\mathbb{E}[t|\bold{x}] = \int t p(t|\bold{x}) dt = y(\bold{x}, \bold{w})

以上の議論を単一の $t, \bold{x}$ から $\bold{t}, \bold{X}$ へ拡張する。

\begin{align*} p(\bold{t} | \bold{X}, \bold{w}, β) & = \prod_{n=1}^N \mathcal{N}(t_n | y(\bold{x}_n, \bold{w}), β^{-1}) \\ & = \prod_{n=1}^N \mathcal{N}(t_n | \bold{w}^\top \bold{x}_n, β^{-1}) \\ ⇔ \ln p(\bold{t} | \bold{X}, \bold{w}, β) & = \sum_{n=1}^N \ln \mathcal{N}(t_n | \bold{w}^\top \bold{x}_n, β^{-1}) \\ & = \frac{N}{2} \ln β - \frac{N}{2} \ln 2π - βE(\bold{w}) \end{align*}

ただし $E$ は二乗和誤差 $E(\bold{w}) = \frac{1}{2} \sum_{n=1}^N (t_n - \bold{w}^\top \bold{x}_n)^2$ 。
これにより尤度 $p(\bold{t} | \bold{X}, \bold{w}, β)$ が書き下せたので、後は 最尤推定 で $\bold{w}, β$ を決定するだけになった。 $\bold{w}$ の最大化は $E$ を最小化することと等価(尤度に $\bold{w}$ が含まれるのは最後の $E$ の項だけだから)なので、この尤度の最大値の必要条件は $\frac{∂E}{∂\bold{w}} = 0$ である。

\begin{align*} 0 & = \frac{∂}{∂\bold{w}} p(\bold{t} | \bold{X}, \bold{w}, β) \\ ⇔ 0 & = \frac{∂}{∂\bold{w}} E(\bold{w}) \\ ⇔ 0 & = \frac{∂}{∂\bold{w}} \frac{β}{2} \sum_{n=1}^N (t_n - \bold{w}^\top \bold{x}_n)^2 \\ ⇔ 0 & = - β \sum_{n=1}^N (t_n - \bold{w}^\top \bold{x}_n) \bold{x}_n^\top \\ ⇔ 0 & = \sum_{n=1}^N t_n \bold{x}_n^\top - \bold{w}^\top \sum_{n=1}^N \bold{x}_n \bold{x}_n^\top \\ ⇔ \underset{\text{Max Likelihood}}{\bold{w}_{\text{ML}}} & = (\bold{X}^\top \bold{X})^{-1} \bold{X}^\top \bold{t} \end{align*}

特に $\bold{X}^{†} = (\bold{X}^\top \bold{X})^{-1} \bold{X}^\top$ は Moore-Penrose疑似逆行列 で、正方行列でない行列の逆行列の性質を持つものである。
この式は最小二乗法の 正規方程式(normal equation)と呼ばれていて、データ行列 $\bold{X}$ は計画行列(design matrix)と呼ばれる。

また、バイアスパラメータ $\bold{x}_0$ を明示的に扱うと以下のようになり、 $\bold{x}_0$ が求められる。

\begin{align*} 0 & = \frac{∂}{∂w_0} E(\bold{w}) \\ ⇔ 0 & = \frac{∂}{∂w_0} \frac{β}{2} \sum_{n=1}^N (t_n - (w_0 + \bold{w}^\top \bold{x}_n ))^2 \\ ⇔ 0 & = \sum_{n=1}^N t_n - \underbrace{\sum_{n=1}^N}_{↳N} w_0 - \sum_{d=1}^D w_d \sum_{n=1}^N x_{n,d} \\ ⇔ w_0 & = \underbrace{\frac{1}{N} \sum_{n=1}^N t_n}_{\bar{\bold{t}}} - \sum_{d=1}^D w_d \underbrace{\frac{1}{N} \sum_{n=1}^N x_{n,d}}_{\bar{\bold{x}}_n} \end{align*}

また、別の話題なので深く触れないが、分散 $β^{-1}$ についても最尤解が求められて、それは回帰関数 $y(\bold{x}, \bold{w})$ 周囲における目標値 $t_n$ の残差分散になる。

β_{\text{ML}}^{-1} = \frac{1}{N} \sum_{n=1}^N (t_n - \bold{w}_{\text{ML}^\top \bold{x}_n})^2

以上の議論で、データ点にガウスノイズが加わったと見てその期待値を最尤推定する導出の最小二乗法が理解できた。
これには先述の問題点があり、正規方程式に出てくる擬似逆行列が非正則に近い場合、求解で数値計算的な不安定が発生する。この問題は特異値分解(SVD) や、正則化項を使って無理やり正則な行列にしてしまうことで解決できる。

SEwPにおける導出

(PRMLに比べると少々天下り的でわかりにくい解説だったため、PRMLの方を理解していれば読み飛ばしてもいいと思う)

$\bold{X}, \bold{t}$ を $y(\bold{x}, \bold{w}) = w_0 + \bold{x}^\top \bold{w}$ でfitする問題を考える。このとき簡単のため $\bold{X}, \bold{t}$ を標準化し、中心を $\bold{0}$ にして切片が $w_0 = 0$ となるようにしておく。
シンプルな最小二乗法を使うと次の誤差関数を最小化する問題に帰着される。

\begin{align*} E(\bold{w}) & = \frac{1}{2} || \bold{t} - [y(\bold{x}, w)]_n ||^2 \\ & = \frac{1}{2} || \bold{t} - (\bold{w}_0 + \bold{X} \bold{w}) ||^2 \end{align*}

$||\bold{z}|| = \sqrt{\sum_{d=1}^D z_d^2}$ : L2ノルム
$\bold{w}_0$ は全要素 $w_0$

この設定で $E$ を最小化する。必要条件は $\frac{∂E}{∂w_0} = 0$ で、以下のように標準化が効いて $w_0$ が消えることがわかる。

\begin{align*} 0 & = \frac{∂E(\bold{w})}{∂w_0} \\ ⇔ 0 & = \frac{∂}{∂w} \sum_{n=1}^N (t_n - (w_0 + \sum_{d=1}^D x_{n,d} w_d))^2 \\ ⇔ 0 & = - N(\underset{↳0}{\bar{\bold{t}}} - w_0 - \sum_{d=1}^D \underset{↳0}{\bar{\bold{x}}_d} w_d) \\ ⇔ w_0 & = \frac{0}{N} \end{align*}

というわけでデータが標準化されている設定では切片 $w_0$ は無いものとして良い。
ここで、 $E$ のベクトル微分について次の式から $y(\bold{x}, \bold{w})$ の重み $\bold{w}$ を閉形式で得ることができる。

\begin{align*} 0 & = \frac{∂E(\bold{w})}{∂\bold{w}} \\ ⇔ 0 & = \begin{bmatrix} \frac{∂}{∂w_1} \sum_{n=1}^N (t_n - \sum_{d=1}^D x_{n,d} w_d)^2 \\ \vdots \\ \frac{∂}{∂w_D} \sum_{n=1}^N (t_n - \sum_{d=1}^D x_{n,d} w_d)^2 \\ \end{bmatrix} \\ ⇔ 0 & = - \begin{bmatrix} x_{1,1} & \dots & x_{N,1} \\ \vdots & \ddots & \vdots \\ x_{D,1} & \dots & x_{N,D} \\ \end{bmatrix} \begin{bmatrix} t_1 - \sum_{d=1}^D x_{1,d} w_d \\ \vdots \\ t_N - \sum_{d=1}^D x_{N,d} w_d \\ \end{bmatrix} \\ ⇔ 0 & = - \bold{X}^\top (\bold{t} - \bold{X} \bold{w}) \\ \end{align*}

これを解くことで、 $E$ を最小化する重み $\hat{\bold{w}}$ が求まる。

\bold{w}_{\text{ML}} = (\bold{X}^\top \bold{X})^{-1} \bold{X}^\top \bold{t}

また、標準化しない場合、切片 $w_0$ は $\bold{w}_{\text{ML}}$ を用いて表すと次になる。

w_0 = \bar{\bold{t}} - \sum_{d=1}^D \bar{\bold{x}}_d w_{d \text{ML}}

これにより、 $\bold{X}^\top \bold{X}$ が正則な場合は解が求まるが、 $N<D$ となる場合(SEwPでは「スパースな状況」と呼んでいるがデータ自体にスパース性はないため誤りだと思うより正確には「特徴次元が過大な状況」か)、及び特徴列に重複するコピーが含まれている場合では $\bold{X}$ が非正則になる。

\text{rank}(\bold{X}^\top \bold{X}) ≤ \text{rank}(\bold{X}) ≤ \min(N,D) = N < D

このような状況では、「上手く効く変数」を選択し、いらない変数を消すことで $N\geqD$ に整形することで困難を突破する。
例えば情報量基準や分散などを使う方法があるが、スパース推定の文脈では正則化項を加えた $E$ を解析することで非正則な逆行列の計算を回避する。

線形回帰の実装

SEwP 1.5 の通り、6次元×500サンプルのデータで、3次元ずつ強い相関があるデータを生成する。
実行結果は以下。

# データの準備
# data: 3列

N = 500
D = 6
data = np.zeros([N, D+1])
data_tmp = np.zeros([N, 2])
data_tmp[:,0] = np.random.normal(0, 1, N)
data_tmp[:,1] = np.random.normal(0.1, 1, N)
data[:,-1] = 4 * data_tmp[:,0] - 2 * data_tmp[:,1] + 3 * np.random.normal(0, 1, N)
data[:,0] = data_tmp[:,0] + 0.5 * np.random.normal(0, 1, N)
data[:,1] = data_tmp[:,0] + 0.5 * np.random.normal(0, 1, N)
data[:,2] = data_tmp[:,0] + 0.5 * np.random.normal(0, 1, N)
data[:,3] = data_tmp[:,1] + 0.5 * np.random.normal(0, 1, N)
data[:,4] = data_tmp[:,1] + 0.5 * np.random.normal(0, 1, N)
data[:,5] = data_tmp[:,1] + 0.5 * np.random.normal(0, 1, N)


data = (data - data.mean(keepdims=True))/data.std()  # 標準化により平均0に移動
data_train, data_val = train_test_split(data, test_size=1/5)
x_train = data_train[:, :D]     # (400, 6)
t_train = data_train[:, D:D+1]  # (400, 1)
x_val = data_val[:, :D]         # (100, 6)
t_val = data_val[:, D:D+1]      # (100, 1)


# 逆行列との積の計算 utility
# x: 逆を求めたい行列
# t: xを消した右辺
# w: 求めたいパラメータw
def solve(x, t):
    w = np.linalg.solve(x, t)
    return w


# 絶対誤差 utility
def error(p, t):
    return np.abs(p - t).mean()

max_likelihood.py

# モデル y(x,w_ML) 
# x: データ行列
# t: ターゲット

class ModelMl:
    def __init__(self):
        return None
    
    def fit(self, x, t):
        self.w = solve(x.T @ x, x.T @ t)
        return None
    
    def infer(self, x):
        pred = x @ self.w
        return pred


model = ModelMl()
model.fit(x_train, t_train)

pred = model.infer(x_train)
print("train error = ", error(pred, t_train))
# train error =  1.1098709629244718

pred = model.infer(x_val)
print("val error = ", error(pred, t_val))
# val error =  1.1842113090570163

2. Redge回帰

擬似逆行列がやな感じにならないためには、(答えを先にいうと) 擬似逆行列に含まれる逆部分 $\bold{X}^\top \bold{X}$ に単位行列 $\bold{I}$ を加算して正則にすればよい。

\bold{w}_{\text{MAP}} = (\bold{X}^\top \bold{X} + λ\bold{I})^{-1} \bold{X}^\top \bold{t}

$λ>0$ : 正則化の強さ

正則化項の導入

PRML 3.1.4

この結果は非常に自然な発想から導かれる。まず、非正則に近い(つまり線形従属に近い列を持つ)行列の数値計算では、値が不安定になりパラメータ $\bold{w}$ が極端に大きな値を取るようなことがある。
最尤推定に出てくる誤差 $E$ へこれを抑制する項を与えて、極端に大きなパラメータが出てきたら $E$ が大きくなってしまうようにすれば良いと考えられる。

E(\bold{w}) = \frac{1}{2} || \bold{t} - \bold{X}\bold{w} ||^2 + \frac{λ}{2} \underset{\text{Weight Decay}}{||\bold{w}||^2}

これを使って、同様に最尤推定で $\bold{w}$ について解く。

\begin{align*} 0 & = \frac{∂}{∂\bold{w}} E(\bold{w}) \\ ⇔ 0 & = \frac{∂}{∂\bold{w}} (\frac{β}{2} \sum_{n=1}^N (t_n - \bold{w}^\top \bold{x}_n)^2 + \frac{λ}{2} \bold{w}^\top \bold{w} ) \\ ⇔ 0 & = - β \sum_{n=1}^N (t_n - \bold{w}^\top \bold{x}_n) \bold{x}_n^\top + λ \bold{w} \\ ⇔ 0 & = \sum_{n=1}^N t_n \bold{x}_n^\top - \bold{w}^\top \sum_{n=1}^N \bold{x}_n \bold{x}_n^\top - λ \bold{w} \\ ⇔ \bold{w}_{\text{MAP}} & = (\bold{X}^\top \bold{X} + λ\bold{I})^{-1} \bold{X}^\top \bold{t} \end{align*}

逆行列部分に正則化が施され、擬似逆行列が正則になった。
$λ>0$ であれば、 $N<D$ となる場合でも $\bold{X}^\top \bold{X} + λ\bold{I}$ 正則であることが保証される。
なぜなら、 $\bold{X}^\top \bold{X}$ は半正定値になるので、固有値 $\{ A_d \}$ は全て非負になるため、 $\bold{X}^\top \bold{X} + λ\bold{I}$ の固有値は $\{ A_d + Nλ \}$ として導かれ、 $Nλ>0$ より全て正になる。則ち、 $λ>0 ⇒ \bold{X}^\top \bold{X} + λ\bold{I}$ が正則である。

この式から分かる通り、誤差 $E$ にパラメータの大きさに関して罰則を与えることは、正規方程式の求解において擬似逆行列の正則化することと等価である。この結果は統計学では パラメータ縮小推定(palameter shrinkage)と呼ばれ、誤差関数が $\bold{w}$ の二次形式で表現できていて解が閉じていることでも便利である。

MAP推定からの導出

ベイズ的な観点で言うと、このパラメータの大きさに関わる正則化は「回帰関数を平均としてデータの分散に正規分布を仮定したモデルに対し、回帰関数のパラメータも正規分布に従っている」という信念を反映した MAP推定(最大事後確率推定)から導ける。MAP推定は最尤推定と異なり、尤度を最大化するのではなく、尤度と事前分布の積(/エビデンス)から計算される事後分布を最大化する。
これから、尤度 $p(\bold{t}|\bold{X}, \bold{w}, β)$ に加えて、推定したいパラメータの事前分布 $p(t | \bold{x}, \bold{w}, β)$ をモデリングすることで、結果的に正則化が得られるという重要な事実を確認していく。

まず、ターゲット $\bold{t}$ が、平均 $y(\bold{x}, \bold{w}) = \bold{w}^\top \bold{x}$ 分散 $β^{-1}$ の正規分布に従うサンプルモデル(=尤度)を考える。

p(\bold{t} | \bold{X}, \bold{w}, β) = \prod_{n=1}^N \mathcal{N}(t_n | y(\bold{x}_n, \bold{w}), β^{-1})

この段階で先述の最尤推定を行った最尤パラメータ $\bold{w}_{\text{ML}}, β_{\text{ML}}$ を尤度に置けば、 予測分布(predictive distribution) $p(\bold{t} | \bold{X}, \bold{w}_{\text{ML}}, β_{\text{ML}})$ が得られる。

ここで $\bold{w}$ についての事前分布を考える。今回尤度が正規分布なので、共役事前分布である正規分布を仮定するのが筋だろう。回帰曲線のパラメータ $\bold{w}$ は、平均(ベクトル) $\bold{0}$ 、共分散行列 $α^{-1}\bold{I}$ の正規分布に従うという信念を与える。

\begin{align*} p(\bold{w}|α) & = \mathcal{N}(\bold{w}|\bold{0}, α^{-1}\bold{I}) \\ & = \left|\frac{α}{2π}\right|^{\frac{D+1}{2}} e^{-\frac{α}{2}||\bold{w}||^2} \end{align*}

ベイズの定理から、尤度と事前分布の積により「データを観測したあとの $\bold{w}$ についての信念」つまり事後分布が得られて、これを最大化する点を推定するのがMAP推定である。

\begin{align*} p(\bold{w} | \bold{X}, \bold{t}, α, β) & = \frac{p(\bold{t} | \bold{X}, \bold{w}, β) p(\bold{w}|α)}{\underbrace{p(\bold{X}, \bold{t}, \bold{w} | α, β)}_{\text{const}}} \\ & ∝ p(\bold{t} | \bold{X}, \bold{w}, β) p(\bold{w}|α) \\ & = \prod_{n=1}^N \mathcal{N}(t_n | y(\bold{x}_n, \bold{w}), β^{-1}) \mathcal{N}(\bold{w}|\bold{0}, α^{-1}\bold{I}) \end{align*}

ここから最小値の条件 $\frac{∂}{∂\bold{w}} p(\bold{w} | \bold{X}, \bold{t}, α, β) = 0$ 、簡単に処理するために対数を取ってこれを行うと次のようになる。

\begin{align*} 0 & = \frac{∂}{∂\bold{w}} \ln p(\bold{w} | \bold{X}, \bold{t}, α, β) \\ ⇔ 0 & = \frac{∂}{∂\bold{w}} ( \underbrace{ -\frac{β}{2}\sum_{n=1}^N (y(\bold{x}_n, \bold{w}) - t_n)^2 + \frac{N}{2}\ln β - \frac{N}{2}\ln 2π }_{\ln p(\bold{t} | \bold{X}, \bold{w}, β)} \underbrace{ -\frac{α}{2} \bold{w}^\top \bold{w} + \frac{D+1}{2} \ln α - \frac{D+1}{2} \ln 2π }_{\ln p(\bold{w}|α)} ) \\ ⇔ 0 & = \frac{∂}{∂\bold{w}} (-\frac{β}{2}\sum_{n=1}^N (y(\bold{x}_n, \bold{w}) - t_n)^2 - \frac{α}{2} \bold{w}^\top \bold{w}) \\ E(\bold{w}) & = \frac{1}{2}\sum_{n=1}^N (y(\bold{x}_n, \bold{w}) - t_n)^2 + \frac{\underset{↳α/β}{λ}}{2} \bold{w}^\top \bold{w}) \end{align*}

$\bold{w}$ を含まない定数部分だけ先に微分で消して、残った項を見てみれば一目瞭然。
これは正則化項を加えた最小二乗法に等しい。

以上より、線形回帰モデルの尤度を最尤推定することが2乗和誤差による最小二乗法に等しく、 $\bold{w}$ の事前分布を与えたMAP推定がこれに正則化項を加えたものとなる。
このL2ノルムによる正則化項 $||\bold{x}||^2$ を使う回帰を Redge回帰 という。

Redge回帰の実装

Redge回帰を実装する。

redge.py

class ModelRedge:
    def __init__(self):
        return None
    
    def fit(self, x, t, lam=1):
        self.w = solve(x.T @ x + lam * np.eye(D), x.T @ t)
        return None
    
    def infer(self, x):
        pred = x @ self.w
        return pred
    

model = ModelRedge()
lam = 1
model.fit(x_train, t_train, lam=lam)

pred = model.infer(x_train)
print("train error = ", error(pred, t_train))
# train error =  1.1100189132415847

pred = model.infer(x_val)
print("val error = ", error(pred, t_val))
# val error =  1.1829521749232386

正則化パラメータの大きさとモデルの精度、各回帰パラメータの大きさを図にした。
$λ$ が大きくなればなるほど $w_d$ は小さくなっているのがわかる。また、適切な $λ > 0$ を与えることで $λ = 0$ の単純な回帰より検証精度が高まっていることもわかる。

redge.py

lam_list = [10**i - 1 for i in np.linspace(0, 3, 50)]

model = ModelRedge()
w_list = []
error_train_list = []
error_val_list = []
for lam in lam_list:
    model.fit(x_train, t_train, lam=lam)
    w_list.append(model.w)
    pred = model.infer(x_train)
    error_train_list.append(error(pred, t_train))
    pred = model.infer(x_val)
    error_val_list.append(error(pred, t_val))

w = np.stack(w_list).squeeze()
error_train = np.stack(error_train_list).squeeze()
error_val = np.stack(error_val_list).squeeze()

fig, axes = plt.subplots(1, 2, figsize=(10,5))
axes[0].plot(lam_list, w)
axes[0].set_xlabel("λ")
axes[0].set_ylabel("w")
axes[1].plot(lam_list, error_train, c="b", label="train error")
axes[1].plot(lam_list, error_val, c="r", label="validation error")
axes[1].set_xlabel("λ")
axes[1].legend()

3. Lasso回帰

正則化項の導入により線形回帰のパラメータが大きな値を取ることを回避することに成功した。
この正則化項をL2ノルムからL1ノルム、つまりパラメータの絶対値に変えることで、驚くべきことに変数選択が可能な正則化が導出される。
Lassoによる変数選択とは、データ $\bold{X}$ 内のなかで、回帰に貢献しない特徴(次元)にかかるパラメータを0にする、つまり特徴を無視することであり、これによりスパース(疎)な解が得られるようになる。

また、Lassoをベイズ的な観点で解釈すれば「パラメータの事前分布に正規分布ではなくラプラス分布を与えた場合のMAP推定解」となる。

劣微分による絶対値の微分

E(\bold{w}) = \frac{1}{2} \sum_{n=1}^D (t_n - \bold{w}^\top \bold{x}_n)^2 + λ ||\bold{w}||_{\text{L1}}

L1ノルム $||\bold{w}||_{\text{L1}} = \sum_{i=0}^I |w_i|$ は $w=0$ で微分不可能な関数であるため、単純な微分により最小化の条件をつけることができない。ここで凸関数の微分を拡張した劣微分を考え、これを突破する。

凸関数: 任意の $0 < s < 1$ と $x_1, x_2 \in \mathbb{R}$ について $f(sx_1 + (1-s)x_2) \leq sf(x_1) + (1-s)f(x_2)$ が成り立つ関数のこと。例えば $f(x) = |x|$ はこれを満たす。また、凸関数同士の和も線形和も凸を保つ。

劣微分: 凸関数 $f: \mathbb{R} \mapsto \mathbb{R}$ と任意の内点 $x_0 \in \mathbb{R}$ に対して $f(x) \geq f(x_0) + c (x-x_0)$ を満たす $\mathcal{C} = \{c \in \mathbb{R} \}$ の集合。劣微分 $∂f(x_0) = \mathcal{C}$ は次の変形で分かる通り、閉区間の中に収まる。

\begin{align*} f(x) & ≥ f(x_0) + c (x-x_0) \\ ⇔ c & ≤ \frac{f(x) - f(x_0)}{x-x_0} \\ ⇔ c & \in [\lim_{x\nearrow x_0} \frac{f(x) - f(x_0)}{x-x_0}, \lim_{x\searrow x_0} \frac{f(x) - f(x_0)}{x-x_0}] \end{align*}

また、劣微分が1点に定まることは、その関数は微分可能であることの必要十分条件である。

今回扱う絶対値関数 $f(x) = |x|$ について考えると、内点 $x_0$ を条件分けして劣勾配を取ることができて、 $f$ の劣微分は次のように表せる。

∂f(x_0) = \begin{cases} \{1\} & x_0 > 0 \\ [-1,1] & x_0 = 0 \\ \{-1\} & x_0 < 0 \\ \end{cases}

劣微分については次のサイトがわかりやすかった。

この結果を使うことでLassoの誤差関数 $E(\bold{w})$ の微分を求めることができ、最小化の条件が作れる。特徴1次元ずつに注目しているので、和の中が微分で $w_d$ 1つ消えることに気をつける。

\begin{align*} 0 & = \frac{∂}{∂w_d} E(w_d) \\ ⇔ 0 & = - \sum_{n=1}^N (t_n - \sum_{d'≠d}^D x_{n,d'}w_{d'} - x_{n,d}w_d) x_{n,d} + λ \begin{cases} 1 & w_d > 0 \\ [-1,1] & w_d = 0 \\ -1 & w_d < 0 \\ \end{cases} \\ ⇔ 0 & = - w_d \sum_{n=1}^N x_{n,d}^2 + \underbrace{\sum_{n=1}^N (t_n - \sum_{d'≠d}^D x_{n,d'}w_{d'}) x_{n,d}}_{:= φ} + \begin{cases} λ & w_d > 0 \\ [-λ,λ] & w_d = 0 \\ -λ & w_d < 0 \\ \end{cases} \\ ⇔ w_d & = \begin{cases} \frac{φ - λ}{\sum_{n=1}^N x_{n,d}^2} & λ < φ \\ 0 & -λ ≤ φ ≤ λ \\ \frac{φ + λ}{\sum_{n=1}^N x_{n,d}^2} & φ < -λ \\ \end{cases} \end{align*}

この式は正則化パラメータ $λ$ を大きくすればするほど $w_d$ が0の値を取る範囲が大きくなることを示している。つまり、正則化パラメータをL2ノルムからL1ノルムに変えれば、最適化時の条件に劣微分という区間条件が入り込んで、データのターゲットへの効きやすさ $φ$ が $λ$ を超えない特徴量次元 $d$ の回帰パラメータ $w_d$ は0になる。
これは驚くべきことだ。L1ノルム正則化の劣微分が生み出す区間という存在が、モデルの最適化と同時にデータを勝手に選択してくれることを意味している。

ただし、解を行列演算で閉じて書くことができない欠点があるため、全次元の解析計算を回して解を収束させるプログラムが要る。

ラプラス事前分布を使った導出

プロの書いたベイズ縮小推定の記事は必見。

詳しくは上のサイトに譲るが、MAP推定において $\bold{w}$ の事前分布として正規分布ではなくラプラス分布を用いる。ラプラス分布は「回帰パラメータは平均付近にまとまって存在しているだろう」という信念を表す 縮小事前分布 で、他にもStudent-T分布なども似た形をしている。

まず通常通り尤度を定義する。

\bold{t} \sim \mathcal{N}(\bold{t} | \bold{X}\bold{w}, β\bold{I})

\begin{align*} \bold{\bold{w}} & \sim p(\bold{w}|β,α) \\ & = \text{Laplace}(\bold{w}|0,β^2,α) \\ & = \prod_{d=1}^D \frac{α}{2β} e^{-\frac{α|w_d|}{β}} \end{align*}

ベイズの定理から事後分布を導くと、最小化条件の微分でいらない項を消したときにLassoの誤差関数が現れる。

\begin{align*} p(\bold{w}|\bold{X}, \bold{t}, β, α) & ∝ \mathcal{N}(\bold{t} | \bold{X}\bold{w}, β\bold{I}) \text{Laplace}(\bold{w}|0,β^2,α) \\ & ∝ e^{-\frac{||\bold{t}-\bold{X}\bold{w}||^2}{2β^2} -\frac{α \sum_{d=1}^D |w_d|}{β}} \\ \frac{∂}{∂\bold{w}} \ln p(\bold{w}|\bold{X}, \bold{t}, β, α) & = \frac{∂}{∂\bold{w}} (\frac{1}{2}||\bold{t}-\bold{X}\bold{w}||^2 + \underset{↳2αβ}{λ} \sum_{d=1}^D |w_d|) \end{align*}

Lasso回帰の実装

Lasso回帰を実装する。プログラムは SEwP 1.3 のものを改造した。

lasso.py

class ModelLasso:
    def __init__(self):
        return None
    
    def fit(self, x, t, lam=1):
        n, dim = x.shape
        eps = 1
        self.w = np.ones(dim)
        w_old = np.ones(dim)

        while eps > 0.00001:
            for d in range(dim):
                phi = t.squeeze()
                for dd in range(dim):
                    phi = (phi - x[:,dd] * self.w[dd]) if (d != dd) else phi
                z = ((x[:, dd] @ phi)/n) / (x[:,d] @ x[:,d]/n)
                self.w[d] = (np.sign(z) * np.maximum(np.abs(z) - lam, np.zeros(1)))[0]
            eps = np.linalg.norm(w_old - self.w)
            w_old = self.w.__copy__()
        return None
    
    def infer(self, x):
        pred = x @ self.w
        return pred


model = ModelLasso()
lam = 1
model.fit(x_train, t_train, lam=lam)

pred = model.infer(x_train)
print("train error = ", error(pred, t_train))
# train error =  1.9019276714302054

pred = model.infer(x_val)
print("val error = ", error(pred, t_val))
# val error =  1.8426633560643888

$λ$ の値に対する挙動を観察する。
確かに正則化を強めるほど0に張り付くパラメータが多くなっている。しかし最尤推定の回帰パラメータと比較すると正にあったはずのものまで $λ ≃ 0$ になっていて、これはおかしい... 実装ミスな気がするので気がついた方いれば教えて下さい。

lasso.py

lam_list = [10**i - 1 for i in np.linspace(0, 0.8, 50)]

model = ModelLasso()
w_list = []
error_train_list = []
error_val_list = []
for lam in lam_list:
    model.fit(x_train, t_train, lam=lam)
    w_list.append(model.w)
    pred = model.infer(x_train)
    error_train_list.append(error(pred, t_train))
    pred = model.infer(x_val)
    error_val_list.append(error(pred, t_val))

w = np.stack(w_list).squeeze()
error_train = np.stack(error_train_list).squeeze()
error_val = np.stack(error_val_list).squeeze()

fig, axes = plt.subplots(1, 2, figsize=(10,5))
axes[0].plot(lam_list, w)
axes[0].set_xlabel("λ")
axes[0].set_ylabel("w")
axes[1].plot(lam_list, error_train, c="b", label="train error")
axes[1].plot(lam_list, error_val, c="r", label="validation error")
axes[1].set_xlabel("λ")
axes[1].legend()

おまけ

高次変動正則化((k, q)-variation regularization): 重みを抑えることで関数形を縮小するのではなく、関数形の変動(微分)にペナルティを与えることで「滑らかな」回帰曲線を得る方法。NNのような非線形モデルの正則化として注目されているらしい。

y^{[k]}(\bold{x},\bold{w}) = \frac{∂^{k}}{∂\bold{x}^{k}} y(\bold{x},\bold{w}) \\ C(y(\bold{x},\bold{w}),k,q) = \int |y^{[k]}(\bold{x},\bold{w})|^q d\bold{x}

正則化項 $C$ が小さくなることは、関数の高次変動が小さくなって回帰関数 $y$ が滑らかになることを意味する。

E(\bold{w}) = \frac{1}{2} ||\bold{t} - y(\bold{X},\bold{w}|| + λ \int |y^{[k]}(\bold{x},\bold{w})|^q d\bold{x}

ただし、積分を含む項は簡単に計算することができない。例えばモンテカルロ積分で数値計算しようとしても、NNのような誤差の反復評価の最適化に、M回の反復積分計算が混入しているのは計算量が厳しすぎる。(特に $D\geq3$ の時 $M$ が大きくなければ精度が怪しい)

\int_Ω |y^{[k]}(\bold{x},\bold{w})|^q d\bold{x} = p \lim_{M→∞} \frac{1}{M} \sum_{m=1}^M \frac{∂}{∂\bold{w}} |y^{[k]}(\bold{z_m \sim \text{U}(Ω)},\bold{w})|^q

驚くべきことに、モデルがNNの場合は無限の反復計算を行う場合、 $M=1$ のような小さな定数で良い。誤差の勾配に関してこれを適用すると、以下の式で反復 $t\to\infty$ の時 $E(\bold{w}^{(t)}) → \min_{\bold{w}} E(\bold{w})$ が言える。

\frac{∂E}{∂\bold{w}} ≃ - \sum_{n=1}^N (t_n - y(\bold{x}_n,\bold{w})) \frac{∂}{∂\bold{w}} y(\bold{x}_n,\bold{w}) + \frac{∂}{∂\bold{w}} |y^{[k]}(\bold{z_m \sim \text{U}(Ω)},\bold{w})|^q

これを勾配として使って最適化を行うことで、滑らかな回帰関数が求められる。
奥野彰文さんの資料によれば、実験では3次の変動を抑えた項を使うことでいい感じに滑らかな回帰曲線が得られている。

感想

正則化線形回帰について、これまでなんとなく使っていた正則化項、重み減衰について行列計算のレイヤーまで落として理解できた。線形回帰は噛めば噛むほど味がするガムのような非常に面白い題材で、基礎トレーニングに最適な問題としていろいろな人に布教していきたいと思った。

SEwPは計算ドリル代わりにと思って買った本だったが、想像の100倍難解で積ん読行きになってしまった。扱っている題材の知識ロードマップとしては非常に優れているように感じるが、説明が天下り的な感じがして、PRMLなどの名著に慣れた脳みそにはハードな書き方になっている...

有名な本で"Statistical Learning with Sparsity ~The Lasso and Generalizations"(Hastie/Tibshirani/Wainwright)があり、パッと見の数式の並び順で言えばこちらのほうが楽しく読めそうなので、英語の勉強も兼ねて頑張るべきかもしれない。 ~~でも英語苦手すぎて続けられる気がしないのでやっぱり邦訳がほしい~~