jiku log

データサイエンスの核心を掴む : 学びと発見の記録

統計検定1級対策に関する記事一覧 #統計検定

はじめに

本記事について

当ブログでは,統計検定1級・統計検定準1級に関連する数理統計学のTips記事や,私が日々学んでいることを紹介している。

統計検定に合格するためには,参考図書を読み問題演習をしっかりこなすことが重要だと考えるが,当ブログでは私が過去に参考書や問題集を読んでいたときに,つまずいたことや理解を深めるために考えたことを紹介する。そのため,「このサイトの内容さえ読めば,受験対策は完璧!」という性質のブログではないのだが,今後統計検定を受験する方にとって,少しでも参考になればと考える。

これらのTips記事は随時更新されるが,記事の数が増えてきたので,「増訂版 日本統計学会公式認定 統計検定1級対応 統計学」の章立てと,「統計検定1級出題範囲表」に沿ってインデックスを作った。本記事はこのインデックスを示すものである。

自己紹介

私は2023年に統計検定1級の統計数理に,2021年に統計検定1級の統計応用(理工学)に合格した。詳細な自己紹介は,下記サイトに記載した。

第6章 統計応用共通手法 ※記事なし

(記事無し)

「確率的機械学習 入門編I」を読む ~第11章 線形回帰 ③リッジ回帰~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第11章 線形回帰」における,リッジ回帰に関する読書メモである。

11.3 リッジ回帰

最尤推定法は過剰適合を起こしやすい。これに対する単純な対策は,重みについて平均0の正規分布による事前分布 p(\boldsymbol{w}) = \mathcal{N}(\boldsymbol{w} \mid \boldsymbol{0}, \lambda^{-1} \boldsymbol{I})を用いたMAP推定を使うことである。これはリッジ回帰と呼ばれる。

MAP推定値は以下のように計算される。


 \begin{align}
\hat{\boldsymbol{w}}_{map} 
&= \operatorname*{argmin}_{\boldsymbol{w}} \mathrm{RSS} (\boldsymbol{w}) +  \lambda \lVert \boldsymbol{w} \rVert _2^2 \\ \\
&= \operatorname*{argmin}_{\boldsymbol{w}} \frac{1}{2 \sigma^2} (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{w})^T (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{w}) + \frac{1}{2 \tau^2} \boldsymbol{w}^T \boldsymbol{w} \tag{1} \\ \\
\end{align}
ただし, \lambda = \sigma^2 / \tau^2である。

11.3.1 MAP推定値の計算

MAP推定値は,次のような罰則付きの目的関数の最小化に対応している。


 \begin{align}
J(\boldsymbol{w}) = (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{w})^T (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{w}) +  \lambda \lVert \boldsymbol{w} \rVert _2^2 \\ \\
\end{align}

この微分を0と置くとMAP推定値が得られる。


 \begin{align}
&\nabla_{\boldsymbol{w}} J(\boldsymbol{w}) = 2(\boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} - \boldsymbol{X}^T \boldsymbol{y} + \lambda \boldsymbol{w}) = \boldsymbol{0} \\ \\
\therefore \: & \hat{\boldsymbol{w}}_{map} = (\boldsymbol{X}^T \boldsymbol{X} + \lambda \boldsymbol{I}_D)^{-1} \boldsymbol{X}^T \boldsymbol{y} \\ \\
\end{align}

QRを用いた求め方

MAP推定値は,式(1)から,


 \begin{align}
\hat{\boldsymbol{w}}_{map} 
= \operatorname*{argmin}_{\boldsymbol{w}} \frac{1}{\sigma^2} (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{w})^T (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{w}) + \frac{1}{\tau^2} \boldsymbol{w}^T \boldsymbol{w} \\ \\
\end{align}
となる。

この式の右辺を f(\boldsymbol{w})とすると,


 \begin{align}
f(\boldsymbol{w}) &= (\tilde{\boldsymbol{y}} - \tilde{\boldsymbol{X}} \boldsymbol{w})^T (\tilde{\boldsymbol{y}} - \tilde{\boldsymbol{X}} \boldsymbol{w}) \\ \\

\tilde{\boldsymbol{X}} &= 
  \begin{pmatrix}
    \boldsymbol{X}/ \sigma \\ 
    \sqrt{\boldsymbol{\Lambda}} \\  
  \end{pmatrix}
 , \quad
\tilde{\boldsymbol{y}} = 
  \begin{pmatrix}
    \boldsymbol{y}/\sigma \\ 
    \boldsymbol{0}_{D \times 1} \\  
  \end{pmatrix}
\\ \\
\boldsymbol{\Lambda} &= (1/\tau^2)\boldsymbol{I} \\ \\

\end{align}
のように通常の最小二乗法における誤差項のようになる。すなわち,

 \begin{align}
\hat{\boldsymbol{w}}_{map} = (\tilde{\boldsymbol{X}}^T \tilde{\boldsymbol{X}})^{-1} \tilde{\boldsymbol{X}}^T \boldsymbol{y} \\ \\
\end{align}
になるので,通常の最小二乗法のようにQR分解を用いて解けばよい。

SVDを用いた求め方

特徴量の方がサンプル数よりも多い場合,すなわち D \gt Nである場合を考える。この場合,QR分解よりもSVDを用いたほうが速い。

 \boldsymbol{X} \in \mathbb{R}^{N \times D}のSVDを,


 \begin{align}
\boldsymbol{X} = \boldsymbol{U} \boldsymbol{S} \boldsymbol{V}^T,\quad \boldsymbol{U} \in \mathbb{R}^{N \times N}, \boldsymbol{S} \in \mathbb{R}^{N \times N}, \boldsymbol{V} \in \mathbb{R}^{N \times D} \\ \\
\end{align}
とする。また, \boldsymbol{R} = \boldsymbol{U} \boldsymbol{S}すなわち \boldsymbol{X} = \boldsymbol{R} \boldsymbol{V}^T とする。

MAP推定値は,


 \begin{align}
\hat{\boldsymbol{w}}_{map} = (\boldsymbol{X}^T \boldsymbol{X} + \lambda \boldsymbol{I}_D)^{-1} \boldsymbol{X}^T \boldsymbol{y} \\ \\
\end{align}
であるが,

 \begin{align}
\hat{\boldsymbol{w}}_{map} = \boldsymbol{V} (\boldsymbol{R}^T \boldsymbol{R} + \lambda \boldsymbol{I}_N)^{-1} \boldsymbol{R}^T \boldsymbol{y} \\ \\
\end{align}
と書き換えられるので,計算量を減らすことができる。

11.3.2 リッジ回帰とPCAの関連

リッジ回帰によって得られる予測値は, \hat{\boldsymbol{y}} = \boldsymbol{X} \hat{\boldsymbol{w}}_{map}であるが,これを \boldsymbol{X}のSVD \boldsymbol{X} = \boldsymbol{U} \boldsymbol{S} \boldsymbol{V}^Tを用いて表現すると,


 \begin{align}
\hat{\boldsymbol{y}} = \boldsymbol{X} \hat{\boldsymbol{w}}_{map} = \sum_{j=1}^D \boldsymbol{u}_j \frac{\sigma^2 }{\sigma^2 + \lambda} \boldsymbol{u}_j^T \boldsymbol{y} \\ \\
\end{align}
となる。

 \boldsymbol{X}の下位の特異値が小さい場合, \sigma^2 / (\sigma^2 + \lambda)における \lambda の影響が大きくなり,結果として予測値における \boldsymbol{u}_j方向の成分が小さくなる。この現象は縮小と呼ばれる。


関連した方法に,主成分回帰と呼ばれる方法もある。これは,まず主成分分析を用いて次元数を K次元まで減らしてから回帰を行なうという方法である。
ただし主成分回帰はリッジ回帰と異なり,下位主成分の情報を除くことになるので,リッジ回帰よりも性能は良くない。

11.3.3 正則化の強さの選び方

正則化のパラメータ \lambdaの選び方は,

  • 交差検証を用いる
  • 経験ベイズを用いて, \hat{\lambda} = \mathrm{argmax}_\lambda \log p(\mathcal{D} \mid \lambda)とする

といった方法が挙げられる。

まとめと感想

「第11章 線形回帰」における,リッジ回帰についてまとめた。

リッジ回帰は,安定的に最小二乗法を適用するための方法であり,実務上でもよく用いられる。特徴量の数がデータの数よりも多い場合( D \gt N)には,特に有用である。

リッジ回帰のパラメータの求め方として,QR分解を用いる方法とSVDを用いる方法があるが,これは特徴量の数とデータの数によって良し悪しが変わるので,状況に応じて使い分けられるようにしたい。

本記事を最後まで読んでくださり,どうもありがとうございました。

「確率的機械学習 入門編I」を読む ~第11章 線形回帰 ②最小二乗法・その2~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第11章 線形回帰」における,最小二乗法に関する読書メモである。

11.2 最小二乗法による線形回帰

11.2.3 最尤推定値を計算する他の手法

本項では,最尤推定値を計算する方法を説明している。

オフセットと傾きを別々に解く

線形回帰の式は, p(y \mid \boldsymbol{x}, \boldsymbol{\theta}) = \mathcal{N}(y \mid w_0 + \boldsymbol{w}^T \boldsymbol{x}, \sigma^2)という形をしている。

オフセット(またはバイアス)は w_0であり,傾きは \boldsymbol{w}である。これを求めるには,

  • 中心化された入力,すなわち \boldsymbol{x}_n^c = \boldsymbol{x}_n - \bar{\boldsymbol{x}},またはこれをを並べた入力行列 \boldsymbol{X}_c
  • 中心化された出力 \boldsymbol{y}_c = \boldsymbol{y} - \bar{\boldsymbol{y}}

を用いて,以下のように推定する。


 \begin{align}
\hat{\boldsymbol{w}} = (\boldsymbol{X}_c^T \boldsymbol{X}_c)^{-1} \boldsymbol{X}_c^T \boldsymbol{y}_c \\ \\
\hat{w}_0 = \bar{y} - \bar{\boldsymbol{x}}^T \hat{\boldsymbol{w}} \\ \\
\end{align}

単回帰(1次元の入力)

1次元では, C_{xy} = \mathrm{Cov}[X, Y]および C_{xx} = \mathrm{Cov}[X, X] = V[X]を用いて,以下のように表される。


 \begin{align}
\hat{w}_1 = \frac{C_{xy}}{C_{xx}} \\ \\
\hat{w}_0 = E[y] - w_1 E[x] \\ \\
\end{align}

偏回帰

入力が2次元のとき,すなわち Y = w_0 + w_1 X_1 + w_2 X_2 + \epsilonとなる場合を考える。

 X_2を一定にした X_1に対する Y偏回帰係数 R_{YX_1 \cdot X_2}は,以下で与えられる。


 \begin{align}
w_1 = R_{YX_1 \cdot X_2} = \frac{ \partial }{ \partial x } E[ Y \mid X_1=x, X_2 ] \\ \\
\end{align}

3変数以上の場合でも,着目している変数 x_j以外の変数を全て一定にしたときの x_jの単位変化に対する出力 yの変化量である。

生成モデルの視点からの最尤推定値の導出

線形回帰モデルは,生成モデルを用いて E[y \mid \boldsymbol{x} ] = f(\boldsymbol{x})を用いて表すことができる。

多変量正規分布を用いて p(\boldsymbol{x}, y)(同時分布)の当てはめを行なうとする。
このとき,同時分布の最尤推定量は,経験平均・経験共分散になる。


 \begin{align}
\boldsymbol{\mu}_x &= \frac{1}{N} \sum_n \boldsymbol{x} \\ \\
\mu_y &= \frac{1}{N} \sum_n y_n \\ \\
\boldsymbol{\Sigma}_{xx} &= \frac{1}{N} \boldsymbol{X}_c^T \boldsymbol{X}_c \\ \\
\boldsymbol{\Sigma}_{xy} &= \frac{1}{N} \boldsymbol{X}_c^T \boldsymbol{y}_c \\ \\
\end{align}

多変量正規分布の条件付き期待値を用いると,次式を得る。


 \begin{align}
E[ y \mid \boldsymbol{x} ] = \mu_y + \boldsymbol{\Sigma}_{xy}^T \boldsymbol{\Sigma}_{xx}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}_x) \\ \\
\end{align}
この式と, w_0, \boldsymbol{w}の対応を見ると,以下のようになる。

 \begin{align}
w_0 &= \mu_y - \boldsymbol{w}^T \boldsymbol{x} = \bar{y} - \boldsymbol{w}^T \boldsymbol{x} \\ \\
\boldsymbol{w} &= \boldsymbol{\Sigma}_{xx}^{-1} \boldsymbol{\Sigma}_{xy} = (\boldsymbol{X}_c^T \boldsymbol{X}_c)^{-1} \boldsymbol{X}_c^T \boldsymbol{y}_c \\ \\
\end{align}

●ブログ筆者註 :
多変量正規分布の条件付き分布については,以下の過去記事に解説している。
stern-bow.hatenablog.com

11.2.4 適合度の指標

本項では,回帰モデルがデータにどれくらい適合しているか,という適合度(goodness of fit)を評価する簡単な方法を説明している。

残差プロット

残差プロットは,1次元の入力 x_nに対する残差 r_n = y_n - \hat{y}_nの散布図である。

残差は \mathcal{N}(0, \sigma^2)にしたがうことを仮定しているため,残差プロットは r=0の直線の周りに均等にばらつくことが期待される。

下図に,本書図1.7に示したデータに1次式・2次式を当てはめたときの残差プロットを示す。
左図が,1次式を当てはめたときの残差プロットである。これは残差が曲線上に分布しているので,モデルが適合していないと言える。
右図が,2次式を当てはめたときの残差プロットである。これは残差が r=0の周りにばらついているので,このモデルの方がより適合していると言える。

多項式回帰における残差プロット

入力が多次元の場合は,真の出力を横軸,予測値を縦軸とするようなプロットを作成する。当てはまりが良い場合は,対角線上に点が並ぶことになる。

多項式回帰における実測値・予測値のプロット

まとめと感想

「第11章 線形回帰」における,最小二乗法についてまとめた。

最小二乗法はシンプルな方法ではあるが,実際に最尤推定量を求めようとする際には,様々な工夫が必要になることがある。本項の内容は,最尤推定量を様々な観点で考察しており,最尤推定量が求めづらいときの解決手段を与えているようで参考になった。

また,モデルの適合度に関する議論も,実務において線形モデルを用いたときの妥当性チェックを行なうときに有用であると感じた。
残差プロットは簡単なチェック方法ではあるが,その理論的な裏付けも説明されており,分析結果を他者に説明する際にも参考にできる。

本記事を最後まで読んでくださり,どうもありがとうございました。

「確率的機械学習 入門編I」を読む ~第11章 線形回帰 ①最小二乗法・その1~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第11章 線形回帰」における,最小二乗法に関する読書メモである。

第11章 線形回帰

11.1 はじめに

本章では,線形回帰について説明している。線形回帰は

  • 入力 \boldsymbol{x} (独立変数,説明変数,共変量,ともいう)
  • 出力 y (従属変数,ターゲット,ともいう)

について,出力の期待値は入力の線形関数 E[y \mid \boldsymbol{x} ]になることを仮定している。

11.2 最小二乗法による線形回帰

11.2.1 用語

本項では,用語やモデルについて説明している。

線形回帰

線形回帰モデルは,通常,以下のモデルのことを指す。


 \begin{align}
p(y \mid \boldsymbol{x}, \boldsymbol{\theta}) = \mathcal{N}(y \mid w_0 + \boldsymbol{w}^T \boldsymbol{x}, \sigma^2) \\ \\
\end{align}
入力が1次元のときは単回帰,入力が多次元のときは重回帰と呼ぶ。

多変量線形回帰

出力が多次元になる場合は,多変量線形回帰と呼ばれ,以下の式で表される。


 \begin{align}
p(\boldsymbol{y} \mid \boldsymbol{x}, \boldsymbol{W}) = \prod_{j=1}^J \mathcal{N}(y_j \mid \boldsymbol{w}_j^T \boldsymbol{x}, \sigma_j^2) \\ \\
\end{align}

特徴量抽出器

非線形なデータに対して線形回帰を適用する場合,特徴量抽出器と呼ばれる関数 \boldsymbol{\phi}によって非線形変換 \boldsymbol{\phi} (\boldsymbol{x})を用いれば,


 \begin{align}
p(y \mid \boldsymbol{x}, \boldsymbol{\theta}) = \mathcal{N}(y \mid \boldsymbol{w}^T \boldsymbol{\phi} (\boldsymbol{x}), \sigma^2) \\ \\
\end{align}
というモデルになる。これは,モデルは入力について線形でなくても,パラメータに対しては線形である。
具体例としては,多項式回帰が挙げられる。

11.2.2 最小二乗推定

最小二乗推定では,訓練データにおける負の対数尤度を最小化する。

負の対数尤度の最小化は,以下で表される残差二乗和


 \begin{align}
\mathrm{RSS}(\boldsymbol{w}) = \frac{1}{2} (\boldsymbol{X} \boldsymbol{w} - \boldsymbol{y})^T (\boldsymbol{X} \boldsymbol{w} - \boldsymbol{y}) \\ \\
\end{align}
の最小化と同等である。

最小二乗法

RSSを最小にする \boldsymbol{w}は,RSSの勾配を0と置いた正規方程式を解けばよい。これで得られる解は最小二乗(ordinary least squares, OLS)解であり,


 \begin{align}
&\nabla_{\boldsymbol{w}}\mathrm{RSS}(\boldsymbol{w}) = \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} - \boldsymbol{X}^T \boldsymbol{y} = \boldsymbol{0} \\ \\
\therefore &\hat{\boldsymbol{w}} = (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y} \\ \\
\end{align}
で得られる。

なおRSSのヘッセ行列が \boldsymbol{H}(\boldsymbol{w}) = \boldsymbol{X}^T \boldsymbol{X}で与えられ,これは正定値行列なので,最小二乗法の目的関数は一意な大域的最適解を持つ。
そのため最小二乗解は一意に定まる。

RSSの等高線

最小二乗法の幾何学的解釈

●ブログ筆者註 :
本項の内容は,本ブログの過去記事の内容に近いものであった。
stern-bow.hatenablog.com

計算上の課題

OLS解 \hat{\boldsymbol{w}} = (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y}の問題点は,逆行列 (\boldsymbol{X}^T \boldsymbol{X})^{-1}が悪条件になっていたり,非正則になっていたりする可能性がある点である。

回避策としては,

  • 特異値分解を用いて擬似逆行列を計算する
  • 正規方程式をQR分解を用いて解くことにより,逆行列計算を行なわずに最小二乗解を得る

といったことが挙げられる。

●ブログ筆者註 :
QR分解については,過去記事にまとめた。
stern-bow.hatenablog.com

重み付き最小二乗法

最小二乗法の式は,各事例 \boldsymbol{x}_nで出力を予測した時の分散 \sigma^2が同じものであることを仮定している。
しかし不均一分散回帰では,各事例の分散が異なるため,モデルは


 \begin{align}
&p(\boldsymbol{y} \mid \boldsymbol{x}; \boldsymbol{\theta}) = \mathcal{N}(\boldsymbol{y} \mid \boldsymbol{X} \boldsymbol{w}, \boldsymbol{\Lambda}^{-1}) \\ \\
&\boldsymbol{\Lambda} = \mathrm{diag}(1/\sigma^2(\boldsymbol{x}_n)) \\ \\
\end{align}
のようになる。これは重み付き線形回帰の式であり,最尤推定値は,

 \begin{align}
\hat{\boldsymbol{w}} = (\boldsymbol{X}^T \boldsymbol{\Lambda} \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{\Lambda} \boldsymbol{y} \\ \\
\end{align}
となる。

まとめと感想

「第11章 線形回帰」における,最小二乗法についてまとめた。

最小二乗法は,実務のデータ分析において最も基本的な方法であり,復習的な内容も多かった。

ただ本章で特筆するべき点は,紹介されているモデルの種類の多さである。単回帰・重回帰だけでなく,多変量線形回帰や非線形化,不均一分散回帰などさまざまな派生形を紹介していた。
このようにモデルの種類とその違いを理解しておくと,実務で利用して精度が出なかったときに,モデルの修正案として利用することができるためである。

本記事を最後まで読んでくださり,どうもありがとうございました。

「確率的機械学習 入門編I」を読む ~第10章 ロジスティック回帰 ⑧ベイズロジスティック回帰~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第10章 ロジスティック回帰」における,ベイズロジスティック回帰に関する読書メモである。

10.5 ベイズ的ロジスティック回帰

これまでの節では,パラメータの点推定に着目していた。本節では,データの不確実性を捉えるために,事後確率 p(\boldsymbol{x} \mid \mathcal{D})を計算する,ベイズロジスティック回帰について説明している。

ロジスティック回帰の場合,事後確率を正確に求めることは不可能なので,近似を用いることになるが,本節ではラプラス近似を用いた方法について説明している。

10.5.1 ラプラス近似

ラプラス近似は,MAP推定量 \hat{\boldsymbol{w}}と,MAP推定量で計算されるヘッセ行列 \boldsymbol{H}を用いて,事後分布を


 \begin{align}
p(\boldsymbol{w} \mid \mathcal{D}) \approx \mathcal{N}(\boldsymbol{w} \mid \hat{\boldsymbol{w}}, \boldsymbol{H}^{-1}) \\ \\
\end{align}
で近似する。

下のFigure 10.13にベイズロジスティック回帰におけるパラメータ \boldsymbol{w}の尤度やラプラス近似を示している。詳細は後述する。

ベイズロジスティック回帰の例

データと決定境界(ブログ筆者補足)

Figure 10.13の(a)(左上)は,データと決定境界を示している。

2クラス分類問題のデータであり,(a)図では決定境界の例として4つの直線を示している。4つの直線の式は, \boldsymbol{w}^T \boldsymbol{x} = 0を示しており,重み \boldsymbol{w}について,

  •  \boldsymbol{w}_1 = (3, 1)
  •  \boldsymbol{w}_2 = (4, 2)
  •  \boldsymbol{w}_3 = (5, 3)
  •  \boldsymbol{w}_4 = (7, 3)

に相当している。

なおこのデータにおいて,最尤推定量は \hat{\boldsymbol{w}}_{mle} = (8.0, 3.4)である。

尤度関数(ブログ筆者補足)

Figure 10.13の(b)(右上)は,尤度関数の等高線を示している。

尤度関数 \mathcal{L}(\boldsymbol{w}) \boldsymbol{w} \in \mathbb{R}^2を入力にしているので,2次元上に表示される。各点は, \boldsymbol{w}=(\boldsymbol{w}(1), \boldsymbol{w}(2))のような成分を持つ重みに該当する。

図中の1,2,3,4は,重み \boldsymbol{w}_1, \boldsymbol{w}_2, \boldsymbol{w}_3, \boldsymbol{w}_4に相当している。
また原点から右上方向に延びている直線は,傾きが \hat{\boldsymbol{w}}_{mle}の直線である。

 \boldsymbol{w}_1, \boldsymbol{w}_2, \boldsymbol{w}_3, \boldsymbol{w}_4の各点はいずれも,ある定数をかけると \hat{\boldsymbol{w}}_{mle}に近づくため,各点はこの直線に近い場所に位置する。

正規化されていない事後分布(ブログ筆者補足)

Figure 10.13の(c)(左下)は,正規化されていない事後分布(の対数)を示している。

事後分布は,


 \begin{align}
p(\boldsymbol{w} \mid \mathcal{D}) \propto f(\mathcal{D} \mid \boldsymbol{w}) p(\boldsymbol{w}) \\ \\
\end{align}
のように,尤度 f(\mathcal{D} \mid \boldsymbol{w})と事前分布の積に比例する。単に掛け算しただけでは,事後分布の積分が1になる保証がないので,「正規化されていない事後分布」と言える。

事前分布には,ガウス事前分布 \mathcal{N}(\boldsymbol{w} \mid \boldsymbol{0}, \sigma^2 \boldsymbol{I})を用いる。この図では, \sigma^2=100を用いている。

また,MAP推定量は赤点で示されている。

ラプラス近似(ブログ筆者補足)

Figure 10.13の(d)(右下)は,事後分布をラプラス近似した結果を示している。

ラプラス近似は事後分布を正規近似する。そのため図は,この正規分布の集中楕円の等高線を示している。

近似前の(c)と近似後の(d)は形状が異なるが,最頻値(MAP推定量)を示す赤点は一致している。

またこの集中楕円において,長軸方向は(b)で示したMAP推定量の直線と同じ方向になっており,これは \boldsymbol{w}の大きさに関する不確実性を示している。

また単軸方向は長軸方向と直交しており,決定境界の向き(法線方向)に関する不確実性を捉えている。

事後予測分布

下のFigure 10.14に,事後予測分布の例を示した。詳細は後述する。

ロジスティック回帰モデルの事後予測分布

MAP推定量の等高線(ブログ筆者補足)

Figure 10.14の(a)(左上)は,MAP推定量をパラメータに代入した事後予測分布 p(y=1 \mid \boldsymbol{x}, \hat{\boldsymbol{w}}_{map})を示している。

MAP推定量を代入しただけでは,決定境界の不確実性がほとんどないことが分かる。

事後予測分布から得られたサンプル(ブログ筆者補足)

Figure 10.14の(b)(右上)は,事後予測分布から得られたサンプルによる決定境界を示している。

すなわち,確率分布 \mathcal{N}(\boldsymbol{w} \mid \hat{\boldsymbol{w}}, \boldsymbol{H}^{-1})にしたがう乱数 \boldsymbol{w}を複数個生成し,これを用いて直線 \boldsymbol{w}^T \boldsymbol{x} = 0を作成している。

事後予測分布から得られたサンプルでは,決定境界に関する不確実性は大きいことが分かる。

サンプル平均(ブログ筆者補足)

Figure 10.14の(c)(左下)は,(b)で得られた直線を平均化し,等高線で表したものである。いわばモンテカルロ近似と同じことをしていることになる。

プロビット近似(ブログ筆者補足)

Figure 10.14の(d)(右下)は,後述するプロビット近似の出力を図示している。

10.5.2 事後予測分布の近似

本項では,新たなデータ \boldsymbol{x}が得られたときに出力 yを予測する事後予測分布


 \begin{align}
p(y \mid \boldsymbol{x}, \mathcal{D}) = \int p(y \mid \boldsymbol{x}, \boldsymbol{w}) p(\boldsymbol{w} \mid \mathcal{D}) d \boldsymbol{w} \\ \\
\end{align}
の計算方法について説明している。

プロビット近似

事後予測分布について,単純な方法は点推定量 \hat{\boldsymbol{w}}を用いたプラグイン近似 p(y \mid \boldsymbol{x}, \hat{\boldsymbol{w}})を計算することだが,予測の不確実性は得られない。

モンテカルロ近似も考えられるが,これは複数のサンプルを得る必要があるため,計算が高速ではない。

他の近似方法として,プロビット近似が挙げられる。プロビット近似は,事後予測分布の積分の中身について,

ことによって,積分計算を標準正規分布の累積分布関数で表すという方法である。


上記の近似を用いると,


 \begin{align}
p(y \mid \boldsymbol{x}, \mathcal{D}) &\approx \int \Phi(\lambda a) \mathcal{N}(a \mid m, v) da
= \Phi \left( \frac{ \lambda m }{ (1+\lambda^2 v)^{1/2} }   \right) \approx \sigma(\kappa(v) m) \\ \\
\kappa &\equiv (1+\pi v /8)^{-1/2} \\ \\
\end{align}
となる。


したがって,ロジットを a = \boldsymbol{x}^T \boldsymbol{w}とすれば,


 \begin{align}
p(y=1 \mid  \boldsymbol{x}, \mathcal{D}) &\approx \sigma(\kappa (v) m) \\ \\
m &= E[a] = \boldsymbol{x}^T \boldsymbol{\mu} \\ \\
v &= V[a] = V[ \boldsymbol{x}^T \boldsymbol{w} ] =  \boldsymbol{x}^T \boldsymbol{\Sigma} \boldsymbol{x} \\ \\
\end{align}
という近似式が得られる。

まとめと感想

「第10章 ロジスティック回帰」における,ベイズなロジスティック回帰についてまとめた。

ロジスティック回帰モデルをベイズ化することにより,予測の不確実性を扱えるようになるが,事後分布は解析的に求められないことが多いので,近似を用いることになる。

本節ではラプラス近似を紹介していたが,図の内容を細かく見ていくことで,具体的にどのような計算をしているのか,理解が深まった。

またプロビット近似は,シグモイド関数とプロビット関数の形が似ていることを利用した近似であり,興味深かった。

本記事を最後まで読んでくださり,どうもありがとうございました。

「確率的機械学習 入門編I」を読む ~第10章 ロジスティック回帰 ⑦頑健なロジスティック回帰~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第10章 ロジスティック回帰」における,頑健なロジスティック回帰に関する読書メモである。

10.4 頑健なロジスティック回帰

本節では,データに外れ値(outlier)があっても影響をされにくい,頑健なロジスティック回帰(robust logistic regression)について説明している。

なお外れ値として,主にラベルの付け間違いによるラベル雑音(label noise)を扱っている。

10.4.1 尤度の混合分布モデル

頑健なロジスティック回帰を定義する最も単純な方法は,一般的な条件付モデルを用いて,

  • 確率 \piでランダムなラベル y
  • 確率 1 - \piで通常なラベル

が予測されるように尤度を修正することである。

二値分類の場合は,


 \begin{align}
p(y \mid \boldsymbol{x}) = \pi \mathrm{Ber}(y \mid 0.5) + (1 - \pi) \mathrm{Ber}(y \mid \sigma(\boldsymbol{w}^T \boldsymbol{x})) \\ \\
\end{align}
となる。

●ブログ筆者註 :
上式の意味だが,前半は確率 \piで,「0または1がランダムに出力される」ということを示している。

下図に,外れ値(「×」と記載)を含むロジスティック回帰の例を示す。
左が通常のロジスティック回帰を当てはめた場合である。左図に示すように,決定境界の信用区間の幅が広がっている。
右が頑健なロジスティック回帰を当てはめた場合である。右図では,決定境界の幅は広くなっておらず,外れ値に近いと言える。

外れ値(×)を含むロジスティック回帰の例

10.4.2 二重温度付き損失

本節では,「二重温度付き損失」を用いた頑健なロジスティック回帰を説明している。

温度付き交差エントロピー

基本的なアイディアは,

損失関数が凸な場合は,決定境界から遠く,ラベル付けが間違っているデータ点が悪影響を及ぼす

という問題に対応すること,である。


この問題に対応するために,通常の交差エントロピー損失の代わりに,「温度付け」された損失関数を用いる。
通常の相対エントロピー損失は,真のラベル分布 \boldsymbol{y}と,予測分布 \hat{\boldsymbol{y}}を用いて,


 \begin{align}
\mathcal{L}(\boldsymbol{y}, \hat{\boldsymbol{y}}) = H(\boldsymbol{y}, \hat{\boldsymbol{y}}) = \sum_c y_c \log \hat{y}_c \\ \\
\end{align}

ここで,温度付き交差エントロピー損失を次のように定義する。


 \begin{align}
\mathcal{L}(\boldsymbol{y}, \hat{\boldsymbol{y}}) = \sum_c \left[ y_c (\log _{t_1} y_c - \log _{t_1} \hat{y}_c) - 
\frac{1}{2 - t_1} (y_c^{2-t_1} - \hat{y}_c^{2-t_1} )     \right] \\ \\
\end{align}
ただし, 0 \leq t_1 \leq 1である。

真の分布 \boldsymbol{y}が,クラス cだけが1となる1ホットベクトルである場合は,


 \begin{align}
\mathcal{L}(c, \hat{\boldsymbol{y}}) = - \log _{t_1} \hat{y}_c - \frac{1}{2 - t_1} \left( 1 - \sum_{c'=1}^C \hat{y}_{c'}^{2-t_1}  \right) \\ \\
\end{align}
となる。

なお \log _tは次で与えられる,対数関数の温度付きの一般化である。


 \begin{align}
\log_t(x) \equiv \frac{1}{1 - t} (x^{1-t} - 1) \\ \\
\end{align}
これは単調増加な凹関数(上に凸の関数)であり, t=1のとき通常の対数関数に戻る。

●ブログ筆者註 :
 0 \leq t \leq 1であるため,たとえば t=0.5を考えると,この関数が凹関数であることが分かる。

また関数 g(t) = x^{1-t}について, g(1) = 0 g'(t) = -x^{1-t} \log xであるので,


 \begin{align}
\lim_{t \to 1} \frac{x^{1-t} - 1}{1 - t} = - \lim_{t \to 1} \frac{g(t) - g(1)}{t - 1} = -g'(1) = \log x \\ \\
\end{align}
であることが分かる。

温度付き対数関数は, -1/(1-t)の下界を持つため,温度付き交差エントロピー関数は上界を持つ。
すなわち,下に示すFigure 10.11に示すように,通常のロジスティック損失では外れ値の損失の値は大きくなるが,温度付きロジスティック損失は損失の大きさが頭打ちになる。

通常・温度付きの損失関数と逆リンク関数の比較

温度付きソフトマックス関数

決定境界に近いラベル付けが間違っているデータ点に関しては,ソフトマックス関数よりも裾の重い逆リンク関数を用いる。

通常のソフトマックス関数は,


 \begin{align}
\hat{y}_c = \frac{ a_c }{ \sum_{c'=1}^C \exp (a_{c'}) } = \exp \left[ a_c - \log \sum_{c'=1}^C \exp (a_{c'})  \right] \\ \\
\end{align}
である。

ここで,温度パラメータ t_2 \gt 1を用いて,温度付きソフトマックス関数を以下のように定義する。


 \begin{align}
\hat{y}_c &= \exp_{t_2} (a_c -  \lambda_{t_2} (\boldsymbol{a})) \\ \\
\exp_t(\boldsymbol{x}) &\equiv [ 1 + (1-t)x  ]_+^{1/(1-t)} \\ \\
\boldsymbol{a} &= [ \boldsymbol{w}_1^T \boldsymbol{x}, ..., \boldsymbol{w}_C^T \boldsymbol{x} ]^T \\ \\
\end{align}

2クラスの場合,ソフトマックス関数はシグモイド関数になる。温度付きのシグモイド関数の挙動はFigure10.11の右図にある通りだが,温度付きシグモイド関数は,値が大きく変化せず,外れ値に対しても頑健であることが分かる。

二重温度付きロジスティック回帰

を組合わせることによって,二重温度付きロジスティック回帰と呼ばれる手法が得られる。


下図に,ロジスティック回帰(上段)と頑健なロジスティック回帰(下段)の比較図を示した。データセット

  • (a)ノイズなし
  • (b)決定境界付近にノイズ
  • (c)決定境界から離れたところにノイズ
  • (d)ランダムノイズ

である。

通常のロジスティック回帰では,ノイズが含まれると境界面が広くなることが分かる。一方で頑健なロジスティック回帰の場合,

  • (b)決定境界付近にノイズ : 温度付きソフトマックス関数のパラメータ t_2
  • (c)決定境界から離れたところにノイズ : 温度付き交差エントロピー損失のパラメータ t_1

の調整によって,ノイズにロバストになっていることが分かる。

ロジスティック回帰と頑健なロジスティック回帰の比較

まとめと感想

「第10章 ロジスティック回帰」における,頑健なロジスティック回帰についてまとめた。

外れ値がある場合,事前にデータを確認することによって外れ値を除外することは重要であるが,尤度関数の中に外れ値対策が盛り込まれていることは,二重の対策となるので重要だと感じた。

特にロジスティック回帰の場合,境界から遠い・近いの2パターンにおいて,損失関数やソフトマックス関数を修正することで対策が打てることは興味深かった。関数の形状は違えど, t \to 1の極限において通常の損失関数・ソフトマックス関数と一致するので,実用的なだけでなく理論的な裏付けがなされており,安心して利用できると感じた。

本記事を最後まで読んでくださり,どうもありがとうございました。

「確率的機械学習 入門編I」を読む ~第10章 ロジスティック回帰 ⑥多項ロジスティック回帰・その4~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第10章 ロジスティック回帰」における,多項ロジスティック回帰の各種の話題に関する読書メモである。

10.3 多項ロジスティック回帰

10.3.6 最大エントロピー分類器

多項ロジスティック回帰の式は,


 \begin{align}
p(y=c \mid \boldsymbol{x}, \boldsymbol{W}) = \frac{ \exp(\boldsymbol{w}_c^T \boldsymbol{x}) }{ Z(\boldsymbol{w}, \boldsymbol{x}) }
= \frac{ \exp(\boldsymbol{w}_c^T \boldsymbol{x}) }{ \sum_{c'=1}^C \exp(\boldsymbol{w}_{c'}^T \boldsymbol{x}) } \\ \\
\end{align}
であった。

この式は,各クラスに対して同じ特徴量 \boldsymbol{x})を用いる。一方で,以下のようにクラス依存な特徴量を使うようにモデルを修正できる。


 \begin{align}
p(y=c \mid \boldsymbol{x}, \boldsymbol{W}) = \frac{ 1 }{ Z(\boldsymbol{w}, \boldsymbol{x}) } \exp(\boldsymbol{w}^T \boldsymbol{\phi}( \boldsymbol{x}, c) ) \\ \\
\end{align}
このモデルは,最大エントロピー分類器(maxent分類器)と呼ばれる。


多項ロジスティック回帰は,最大エントロピー分類器に含まれる。すなわち,


 \begin{align}
&\boldsymbol{w} = [\boldsymbol{w}_1, ..., \boldsymbol{w}_C ] \\ \\
&\boldsymbol{\phi}( \boldsymbol{x}, c) = [\boldsymbol{0},..., \boldsymbol{x}, ..., \boldsymbol{0} ] \\ \\
\end{align}
のようにすると, \boldsymbol{w}^T \boldsymbol{\phi}(\boldsymbol{x}, c) = \boldsymbol{w}_c^T \boldsymbol{x}となる。

10.3.7 階層的分類

ラベルの集合が,階層構造(hierarchy)や分類体系(taxonomy)をなしている場合がある。たとえば,

哺乳類━┳犬━┳ゴールデンレトリバー
    ┃  ┣ジャーマンシェパード
    ┃  ┣…
    ┗猫━┳キジトラ
       ┣シャム猫
       ┣…

のような状態である。

ラベルの階層構造の単純な例


多ラベル分類器 p(\boldsymbol{y} \mid \boldsymbol{x})で分類するとき,同じ階層にあるラベルは排他的である,という条件を加えるには,たとえば上図では根元が「Mammals(哺乳類)」と「Birds(鳥類)」なので,


 \begin{align}
p(\mathrm{Mammals} \mid \boldsymbol{x}) + p(\mathrm{Birds} \mid \boldsymbol{x}) = 1 \\ \\
\end{align}
とすればよい。

また次の階層では


 \begin{align}
p(\mathrm{Dogs} \mid \boldsymbol{x}) + p(\mathrm{Cats} \mid \boldsymbol{x}) = p(\mathrm{Mammals} \mid \boldsymbol{x}) \\ \\
\end{align}
とすればよい。

10.3.8 膨大な数のクラスの処理

階層的ソフトマックス関数

通常のソフトマックス分類器では,ラベル数 Cに対して正規化定数の計算は O(C)の計算量だが,ラベル数が多い*1と計算がボトルネックになる。

ラベル集合に木構造を導入すると,どのラベルの確率の計算も,根から葉までのパスに含まれる辺上の確率値のみを計算すればよいので,計算量は O(\log C)になる。
このような木構造を導入した分類器を階層的ソフトマックス関数と呼ぶ。

クラス不均衡とロングテール

クラス数が膨大な場合に現れる問題として,ほとんどのクラスについて非常に少ない事例しか得られないことが挙げられる。
このような場合,希少クラスは頻出クラスに比べて損失全体に与える影響が小さくなるため,モデルは頻出クラスばかりに着目することになる。

この対策として,クラスcのサンプル数 N_cを用いて,バイアス項 \boldsymbol{b} \mathcal{S}(\boldsymbol{b})_c = N_c / Nとなるように設定することである。

このモデルでは,重み \boldsymbol{w} = \boldsymbol{0}のとき,ラベルの経験分布に一致する。このモデルでは,入力に応じてこの事前分布からの乖離度合いを学習することができる。

まとめと感想

「第10章 ロジスティック回帰」における,多項ロジスティック回帰の各種の話題についてまとめた。

テキストデータを対象とした分類では,単語予測のように膨大なクラスを対象にすることがあったり,また階層的なラベル構造を扱う必要が出てくる。
本項で紹介していた手法は,多項ロジスティック回帰から発展して,より実務的な問題に直面したときの対策を紹介しており,とても興味深かった。

製造業においてもテキストデータを扱うことがある。例えば,入力された文章から想定されるリスクのクラスを分類する,といった事例である。
このような場合においては,本項で扱った各種知見が活かせると考えられる。

本記事を最後まで読んでくださり,どうもありがとうございました。

*1:例えば,単語予測などが挙げられる。

「確率的機械学習 入門編I」を読む ~第10章 ロジスティック回帰 ⑤多項ロジスティック回帰・その3~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第10章 ロジスティック回帰」における,多項ロジスティック回帰の目的関数の最適化に関する読書メモである。

10.3 多項ロジスティック回帰

10.3.3 勾配に基づく最適化

前項において,他行ロジスティック回帰の目的関数・勾配・ヘッセ行列を求めた。
勾配を用いると,確率的勾配降下法(SGD)のアルゴリズムを導出できる。
またヘッセ行列を用いると二次最適化法を導出できるが,ヘッセ行列の計算コストが高いので,一般には「制限メモリーBFGS」のような準ニュートン法を用いる。

10.3.4 バウンド最適化

本項では,8.7節で説明したバウンド最適化を用いて最適化する方法を紹介している。

目的関数の下界の設定

最大化したい目的関数 l(\boldsymbol{\theta})が凹関数であるとする(「凹関数」とは,上に凸の関数のことである)。

目的関数のヘッセ行列 \boldsymbol{H}(\boldsymbol{\theta})について, \boldsymbol{H}(\boldsymbol{\theta}) \succ \boldsymbol{B},すなわち \boldsymbol{H}(\boldsymbol{\theta}) - \boldsymbol{B}が正定値であるような,負定値行列 \boldsymbol{B}があるとすると,以下の不等式が成り立つ。


 \begin{align}
l(\boldsymbol{\theta}) &= -\mathrm{NLL}(\boldsymbol{\theta}) \\ \\
&=l(\boldsymbol{\theta}^t) + (\boldsymbol{\theta} - \boldsymbol{\theta}^t)^T \boldsymbol{g}(\boldsymbol{\theta}^t) + \frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}^t)^T \boldsymbol{B} (\boldsymbol{\theta} - \boldsymbol{\theta}^t) \\ \\
\boldsymbol{g}(\boldsymbol{\theta}^t) &= \nabla l (\boldsymbol{\theta}^t) \\ \\
\end{align}

この右辺をバウンド最適化の代理関数 Q(\boldsymbol{\theta}, \boldsymbol{\theta}^t)とすると,更新則は次のようになる。


 \begin{align}
\boldsymbol{\theta}^{t+1} = \boldsymbol{\theta}^t - \boldsymbol{B}^{-1} \boldsymbol{g} (\boldsymbol{\theta}^t) \\ \\
\end{align}

毎反復ごとに値が変わるヘッセ行列 \boldsymbol{H}(\boldsymbol{\theta}^t)を用いるニュートン法と異なり,固定された \boldsymbol{B}を用いて計算するので,計算コストが低い二次最適化法が得られる

多項ロジスティック回帰への適用

バウンド最適化をロジスティック回帰に適用する。

目標は,対数尤度 l(\boldsymbol{w}) = -N \times \mathrm{NLL}(\boldsymbol{w})を最大化することである。

本書P337の式(10.58)において,


 \begin{align}
\mu_{nc} = \mathcal{S}(f(\boldsymbol{x}_n ; \boldsymbol{\theta}))_c = \frac{ \exp(\boldsymbol{w}_c^T \boldsymbol{x}_n) }{ \sum_{c'}\exp(\boldsymbol{w}_{c'}^T \boldsymbol{x}_n)  } \\ \\
\end{align}
であり,

 \begin{align}
\sum_c y_{nc} \log \mu_{nc} &= \sum_c y_{nc} \boldsymbol{w}_c^T \boldsymbol{x}_n - \sum_c y_{nc} \log (\sum_{c'}\exp(\boldsymbol{w}_{c'}^T \boldsymbol{x}_n) ) \\ \\
&= \sum_c y_{nc} \boldsymbol{w}_c^T \boldsymbol{x}_n - \log (\sum_{c}\exp(\boldsymbol{w}_{c}^T \boldsymbol{x}_n) ) \\ \\
\end{align}
となる。ここで, \boldsymbol{y}_nは1ホットベクトルであるため, \sum_{c} y_{nc} = 1であることを用いた。

よって目的関数は,


 \begin{align}
l(\boldsymbol{w}) = \sum_{n=1}^N \left[  \sum_{c=1}^C y_{nc} \boldsymbol{w}_c^T \boldsymbol{x}_n - \log \sum_{c=1}^C \exp(\boldsymbol{w}_{c}^T \boldsymbol{x}_n) \right] \\ \\
\end{align}
となる。

10.3.2項を参考にすると,勾配とヘッセ行列は,それぞれ以下のようになる。


 \begin{align}
\boldsymbol{g} (\boldsymbol{w}) &= \sum_{n=1}^N (\boldsymbol{y}_n - \boldsymbol{\mu}_n (\boldsymbol{w}) ) \otimes \boldsymbol{x}_n \in \mathbb{R}^{CD} \\ \\
\boldsymbol{H} (\boldsymbol{w}) &= -\sum_{n=1}^N (\mathrm{diag}(\boldsymbol{\mu}_n (\boldsymbol{w}) ) - \boldsymbol{\mu}_n (\boldsymbol{w}) \boldsymbol{\mu}_n (\boldsymbol{w}) ^T ) \otimes  (\boldsymbol{x}_n \boldsymbol{x}_n^T) \\ \\
\end{align}

ヘッセ行列の下界 \boldsymbol{B}は以下のように構成できることが知られている。


 \begin{align}
\boldsymbol{H} (\boldsymbol{w}) \succ -\frac{1}{2}[\boldsymbol{I} - \boldsymbol{1} \boldsymbol{1}^T /C  ] \otimes \left( \sum_{n=1}^N \boldsymbol{x}_n \boldsymbol{x}_n^T \right) \equiv \boldsymbol{B} \\ \\
\end{align}

これを用いると,更新則は以下のようになる。


 \begin{align}
\boldsymbol{w}^{t+1} = \boldsymbol{w}^t - \boldsymbol{B}^{-1} \boldsymbol{g} (\boldsymbol{w}^t) \\ \\
\end{align}

この \boldsymbol{B}^{-1}は,あらかじめ計算しておけばよいので,IRLSのようなヘッセ行列を反復ごとに計算する必要がある手法よりも高速な場合がある。

10.3.5 MAP 推定

本書P333の10.2.7項において,ロジスティック回帰におけるMAP推定が,L2正則化と同値であることを示していた。また10.2.7項では,L2正則化の利点(回帰係数が大きくなりすぎない)ことが紹介されていた。

多クラスの場合では,この利点に加えて,パラメータの識別可能性に関する利点も存在する。
尤度を最大化する解が唯一であるとき,パラメータは識別可能であるという。

●ブログ筆者註 :
多値ロジスティック回帰モデルは,識別可能ではない。すなわち,尤度を最大化するパラメータが一意に決まらない。

多値ロジスティック回帰モデルは,


 \begin{align}
p(y = c \mid \boldsymbol{x}, \boldsymbol{W}) = \frac{ \exp(\boldsymbol{w}_c^T \boldsymbol{x}) }{ \sum_{k=1}^C \exp(\boldsymbol{w}_k^T \boldsymbol{x}) } \\ \\
\end{align}
であるが,このパラメータ \boldsymbol{w}_cに定数 \boldsymbol{b}を加えたモデルもまた,

 \begin{align}
\frac{ \exp( (\boldsymbol{w}_c+\boldsymbol{b})^T \boldsymbol{x}) }{ \sum_{k=1}^C \exp( (\boldsymbol{w}_k+\boldsymbol{b})^T \boldsymbol{x}) } 
= \frac{ \exp(\boldsymbol{w}_c^T \boldsymbol{x}) }{ \sum_{k=1}^C \exp(\boldsymbol{w}_k^T \boldsymbol{x}) } \\ \\
\end{align}
となる。すなわちパラメータの平行移動に対する不変性が発生しており,パラメータが一意に定まらないことになる。

L2正則化を付けた目的関数は以下のようになる。


 \begin{align}
\mathrm{PNLL}(\boldsymbol{W}) = -\sum_{n=1}^N \log p(y_n \mid \boldsymbol{x}_n, \boldsymbol{W}) + \lambda \sum_{c=1}^C \lVert \boldsymbol{w}_c \rVert _2^2 \\ \\ 
\end{align}

最適解では,この勾配が0になるので,


 \begin{align}
&\mathrm{NLL}(\boldsymbol{w}) + 2 \lambda \boldsymbol{w} = \boldsymbol{0} \\ \\
\Leftarrow &\sum_n(\boldsymbol{y}_n - \boldsymbol{\mu}_n) \otimes \boldsymbol{x}_n = 2 \lambda \boldsymbol{w} \\ \\
\end{align}
そのため,本書P341にあるように, \sum_c w_{cj} = 0という,各特徴量次元について和を取ると重みは0になる,という制約が生じる。

この制約によって,パラメータが一意に定まるようになるので,L2正則化によって識別可能になる,と言える。

まとめと感想

「第10章 ロジスティック回帰」における,多項ロジスティック回帰の目的関数の最適化についてまとめた。

多項ロジスティック回帰モデルの最適化において,バウンド最適化を用いることで係数の更新が高速化されうる,ということは意外な発見であった。

バウンド最適化は,EMアルゴリズムなど隠れ変数が存在する場合に用いられるものであるという固定概念があったが,ヘッセ行列の下界を上手く設定することにより,固定の行列により勾配計算が可能になるという内容は興味深かった。

ヘッセ行列の下界に関する不等式が登場していた。この式について,

などにも目を通してみたが,よく理解できなかったので,改めて挑戦したいと思った。


本記事を最後まで読んでくださり,どうもありがとうございました。