【1時間で分かる】実務で役に立たなかった正則化・LASSOトリビアの供養😇

先日、LASSOに関する雑学をtweetしたところ、意外なことに複数名の方からレスポンスをいただきました。私の職場(ビジネス支援寄りのデータサイエンティスト集団。笑)では全く使ったことがない知見だったので、こんなトリビアに需要があったのがとても意外に思えました。

しかし、思い返して見ると、かくいう私も、LASSOをはじめて勉強した際、数学的な面白さに魅了されて周辺領域も少しだけ深く調ベていたことがありました。そこで、この記事では、私の実務では全く役に立たなかった正則化、特にLASSOに関する雑学達に触れて、供養してやりたいと思います。(来世で役に立ちますように!)

文量と数式が多く、前提知識なくまともに読むと1時間くらいかかるかもしれません。適宜掻い摘んで読んでいただければと思います。
※「全く役に立たない」と表現しているのは、あくまで私の職場での話です。きっと役に立っているシーンもたくさんあるはずので、あしからず。です。

忙しい人向けのまとめ

  • 正則化とは過学習を防ぐための手法。下記観点で捉えると理解が深まるかも
    • ベイズ観点:回帰係数がゼロに近いという事前確率の前提をおいた手法
    • バイアス・分散観点:バイアスを許容し分散を減らす手法
    • 情報量規準観点:情報量規準の最小化問題における近似解を得る手法
  • LASSOは、スパースな解が得られ(=変数選択ができる)、数学的な性質も良い(=解きやすい)素敵な手法
  • 数学的な性質が良いと言っても解析解は得られないので、数理最適化手法を用いて数値解を求める。でも、そこには偏微分地獄が待っている。。
  • 実務では意識することはないが、LASSOの弱点を補う発展手法が生み出されている

問題設定

登場する記号達

本題に入る前に問題設定を整理しておきます。よくある教師ありの線形回帰問題で、目的変数$y$と説明変数$X$の間に次のような関係を仮定します。
\begin{eqnarray}
y&=&X\beta+\epsilon \\
\epsilon&\sim&\mathcal{N}(0, \sigma^2I)
\end{eqnarray}

あるいは上記をひとまとめにして、
$$y\sim\mathcal{N}(X\beta, \sigma^2I)$$
と書くこともできます。

ここで、$y, X, \beta, \epsilon, \sigma^2, I$はそれぞれ、次のようなものを想定しています。

  • $y$は目的変数を並べたベクトル(目的変数ベクトル)
  • $X$は説明変数ベクトルを並べた行列(説明変数行列)
  • $\beta$は回帰係数を並べたベクトル(回帰係数ベクトル)
  • $\epsilon$は残差を並べたベクトル(残差ベクトル)
  • $\sigma^2$は残差が従う正規分布の標準偏差
  • $I$は単位行列
  • $N$はサンプル数(上式には登場しませんが、あとで登場するのでここで紹介しておきます)

※問題を簡単化するために、$y$は中央化済み($\sum_i y_i = 0$)、$X$は標準化済み($j$番目の説明変数それぞれに対して、$\sum_i X_{ij}=0$, $\sum_i X_{ij}^2=0$)とします。
※本当は、定数項$\beta_0$も加味して、$y=X\beta+\beta_0\vec{1}+\epsilon$という問題を解くのですが、$(X, y)$を中央化すると、$\beta_0$が相殺されて消えます。
※$\beta_0$を含んだ「正しい」定式化については、文末で紹介する参考書籍で丁寧にまとめられていますので、気になる方はそちらを当たって見ると良いかもしれません。

解きたい問題

線形回帰問題では、($X$と)$y$が分かった時に、その状態が最も発生しやすい回帰係数ベクトル$\beta$を推定したい、というモチベーションがあります。数式で表すとこうです。

$$\hat{\beta}=\mathop{\rm arg~max}\limits_\beta P(\beta|y) \tag{1}\label{1}$$

$P(\beta|y)$は事後確率と呼ばれることがあります。$P(\beta|y)$の関数形が不明なので、ベイズの定理を使って次のように変形します。

$$P(\beta|y)=\frac{P(y|\beta)P(\beta)}{P(y)} $$

上記のうち、$P(y)$は、$\beta$の値によって変化しないので、式\eqref{1}は次のような問題に変形できます。

$$\hat{\beta}=\mathop{\rm arg~max}\limits_\beta P(y|\beta)P(\beta) \tag{2}\label{2}$$

こいつを解くことを今回の問題設定としたいと思います。

ちなみに、ブログ出所で恐縮ですが、Suyama Atsushi氏のブログによりますと、

  • ベイズ推定:$P(\beta|y)$の確率分布を推定する問題
  • MAP推定:式\eqref{2}の最適解を求める問題
  • 最尤推定:式\eqref{2}の$P(y|\beta)$のみを最適化する問題

と言います。したがって今回の問題は、「MAP推定」と呼ばれる分野に対応します。「最尤推定」については、次節で少し触れたいと思います。

最尤法との関係

$P(y|\beta)$って何でしょうか?
こいつは、正規分布$\mathcal{N}(X\beta, \sigma^2I)$から、$y$が得られる確率です。尤度とも言います。数式で表現すると以下のようになります。(以下の$N$はサンプル数です。目的変数ベクトル$y$の長さと言っても良いかもしれません)

$$P(y|\beta)=\frac{1}{\sqrt{2\pi\sigma^2}^{N}}\exp \bigl( -\frac{1}{2\sigma^2}(y-X\beta)^T(y-X\beta) \bigr)$$

実は、式\eqref{2}で、$P(\beta)$を無視して、尤度$P(y|\beta)$だけを最大化しようというのが、最尤法の考え方です。特に、今回の問題設定における最尤法は、モデルに線形関数を、残差に正規分布を仮定しているので、最小二乗法となります

数式を追って見ていきます。計算の都合上、尤度の対数とったもの(対数尤度)を最大化することを考えます。最後の式変形では、$\beta$が寄与する部分(残差の二乗和に等しいもの)だけを抽出して最小化問題として再定義しています。

\begin{eqnarray}
\hat{\beta}^{OLS}&=&\mathop{\rm arg~max}\limits_\beta P(y|\beta)\\
&=&\mathop{\rm arg~max}\limits_\beta \log \bigl( P(y|\beta)) \bigr) \\
&=&\mathop{\rm arg~max}\limits_\beta \bigl( -\frac{N}{2}\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}(y-X\beta)^T(y-X\beta) \bigr)\\
&=&\mathop{\rm arg~min}\limits_\beta \bigl( (y-X\beta)^T(y-X\beta) \bigr) \tag{3}\label{3}
\end{eqnarray}

確かに最小二乗法と同じ問題になりました。

ついでに$\beta$の推定値も求めてみます。式\eqref{3}が最小の時、式\eqref{3}を$\beta$で偏微分したものがゼロになるので、$\hat{\beta}^{OLS}$が次のように求められます。ここで、添字の$OLS$は、最小二乗法(Ordinary Least Squares)の意です。

\begin{eqnarray}
\frac{\partial}{\partial \beta} (y-X\beta)^T(y-X\beta)&=&2X^TX\beta-2X^Ty=0\\
\hat{\beta}^{OLS}&=&(X^TX)^{-1}X^Ty
\end{eqnarray}

$X^TX$の逆行列を使う都合上、$X^TX$の性質が悪い(固有値が小さいと言って良いかもしれません)と、良い推定値が得られないことがあります。説明変数に対してサンプル数$N$が少ない場合や(劣決定)、説明変数間の相関が大きい場合に(多重相関性)、$X^TX$の性質が悪くなりがちです。

ちなみに、、ここまであまり触れてきませんでしたが、モデルの残差$\epsilon$の標準偏差$\sigma^2$も未知なので、最尤法により推定することがあります。推定方法は$\beta$と同様で、次の問題を解くことで$\sigma^2$の推定値${\hat{\sigma^2}}^{OLS}$が求められます。

\begin{eqnarray}
\hat{\sigma^2}^{OLS}&=&\mathop{\rm arg~max}\limits_\hat{\sigma^2} P(y|\beta, \sigma^2)\\
&=&\mathop{\rm arg~max}\limits_\hat{\sigma^2} \log \bigl( P(y|\beta, \sigma^2)) \bigr)\\
&=&\mathop{\rm arg~max}\limits_\hat{\sigma^2} \bigl( -\frac{N}{2}\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}(y-X\beta)^T(y-X\beta) \bigr)
\end{eqnarray}

※ここまでの説明で尤度としてきた$P(y|\beta)$の表記を$P(y|\beta, \sigma^2)$と改めておきます。

こちらも$\beta$と同様で、$\mathop{\rm arg~max}$の中の数式を$\sigma^2$で偏微分したものがゼロになるようにすれば、$\sigma^2$の推定値${\hat{\sigma^2}}^{OLS}$が得られることになります。

\begin{eqnarray}
\frac{\partial}{\partial \sigma^2}\log \bigl( P(y|\beta, \sigma^2)) \bigr) &=& -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4}(y-X\beta)^T(y-X\beta)=0\\
\hat{\sigma^2}^{OLS}&=&\frac{1}{N}(y-X\beta)^T(y-X\beta)
\end{eqnarray}

正則化とはなんですの?

本題の前に、まず正則化の定義を確認します。

wikipediaからの引用:
数学・統計学・計算機科学において、特に機械学習と逆問題において、正則化(せいそくか、英: regularization)とは、不良設定問題を解いたり過学習を防いだりするために、情報を追加する手法である。モデルの複雑さに罰則を科すために導入され、なめらかでないことに罰則をかけたり、パラメータのノルムの大きさに罰則をかけたりする。

簡単に言えば、過学習を防ぐための手法です。これまで説明してきた問題設定ですと、$X^TX$の性質が悪い問題をうまく解くための手法とも言えます。

線形回帰における正則化は、次のような定式化がされることがあります。$\hat{\beta}^{Reg}$の添字$Reg$は正則化(Regularization)の意です。

\begin{eqnarray}
\hat{\beta}^{Reg}=&\mathop{\rm arg~min}\limits_\beta (y-X\beta)^T(y-X\beta)+λp(\beta)
\end{eqnarray}

$λp(\beta)$は正則化項(または、ペナルティ項とも)と呼ばれます。正則化項には、モデルの複雑さを表現するような関数が用いられ、「尤度の改善」と、「モデルの複雑さに起因する汎化性能の悪化(過学習)」のトレードオフを考慮することになります。

正則化項をもう少し詳しく見ていきます。

  • $λ$はハイパーパラメータで、何らかの方法(後に示す情報量規準や、機械学習的でも用いられるクロスバリデーション等)で設定されます。このパラメータが大きいほど、「モデルの複雑さ」が重視されやすくなることが、感覚的にイメージできるかと思います。
  • $p(\beta)$には、モデルの複雑さを表現する関数が用いられます。しばしば、$\beta$の$L$-$q$ノルム${\|\beta\|_q}$が用いられます。定義を以下に示しておきます。

$$\|\beta\|_q={}^{q}\sqrt{\sum_i |\beta_i|^q}$$

端的に言えば、$L$-$q$ノルムはベクトルの大きさを表していますが、$q$の値によって、性質が異なります。以下で、代表的なものを紹介したいと思います。

  • $q=0$の場合は少し特殊で、$\beta$に含まれる非ゼロの要素の個数となります。
  • $q=1$の場合は要素の絶対値の和に対応します。この場合の正則化付き線形回帰問題をLASSO回帰と言います。
  • $q=2$の場合はベクトルの長さに対応します。この場合をRIDGE回帰と言います。
  • $q>2$の場合は、、、お目にかかったことがないので説明割愛させてください。。

このブログでは、$p(\beta)$に、$L$-$q$ノルム${\|\beta\|_q}$を用いた、以下の問題について議論を進めて行ければと思います。

\begin{eqnarray}
\hat{\beta}^{Reg}=&\mathop{\rm arg~min}\limits_\beta (y-X\beta)^T(y-X\beta)+λ{\|\beta\|_q}^q \tag{4}\label{4}
\end{eqnarray}

ベイズ的な観点

ベイズ的な観点に立つと、「回帰係数$\beta$がゼロに近いという事前確率の前提をおいた手法」と言えるかもしれません。

式\eqref{2}における$P(\beta)$のことを事前確率と呼びます。

実は式\eqref{4}は、$P(\beta) \propto \exp(-λ{\|\beta\|_q}^{q})$という事前確率の分布を仮定して、式\eqref{2}を解いているに他なりません。

$P(\beta)$は、$q$の値によって次の確率分布に対応します。

  • $q=0$:デルタ関数+一様分布(造語)
  • $q=1$:ラプラス分布
  • $q=2$:正規分布

いずれにしても、平均0の分布となっており、$P(\beta)$の事前確率として、$\beta$は0に近い値であることを仮定していることになります。

バイアスと分散の観点

バイアスと分散の観点では、「正則化とは、バイアスを許容し分散を減らす手法」と言えるかもしれません。

簡単のため、RIDGE回帰($q=2$の場合)を考えます。RIDGE回帰における最適化問題は\eqref{4}に$q=2$を代入して、次のようにかけます。

\begin{eqnarray}
\hat{\beta}^{RIDGE}=&\mathop{\rm arg~min}\limits_\beta (y-X\beta)^T(y-X\beta)+λ{\|\beta\|_2}^2
\end{eqnarray}

最小二乗法と同様に数式の偏微分がゼロになることから、次のように推定値が得られます。

\begin{eqnarray}
\frac{\partial}{\partial \beta} \bigl( (y-X\beta)^T(y-X\beta)+λ{\|\beta\|_2}^2 \bigr)&=&2X^TX\beta−2X^Ty+2λ\beta=0\\
\hat{\beta}^{RIDGE}&=&(X^TX+λI)^{-1}X^Ty\tag{5}\label{5}
\end{eqnarray}

$\hat{\beta}^{RIDGE}$のバイアスと分散はどうなってるんでしょう?本来は、真の$\beta$の値は分かりませんが、仮にこれが既知であると考えましょう。分かりやすく*マークをつけて表記します。

$$y=X\beta^{*}+\epsilon$$

上式を式\eqref{5}に代入すると次関係が得られます。

\begin{eqnarray}
\hat{\beta}^{RIDGE}&=&(X^TX+λI)^{-1}X^T(X\beta^{*}+\epsilon)\\
&=&(X^TX+λI)^{-1}X^TX\beta^{*}+(X^TX+λI)^{-1}\epsilon\\
&=&\bigl( I – λ(X^TX+λI)^{-1} \bigr) \beta^{*}+(X^TX+λI)^{-1}X^T\epsilon\\
\end{eqnarray}

さらに、ベクトルの確率変数$x$(上式の残差ベクトル$\epsilon$に対応しています)に対して、行列$A$, ベクトル$b$による一時変換$Ax+b$を行った場合の平均、分散共分散行列の関係式:
\begin{eqnarray}
\mathop{\rm E}(AX+b)&=&A\mathop{\rm E}(X)+b\\
\mathop{\rm cov}(AX+b)&=&A\mathop{\rm cov}(X)A^T
\end{eqnarray}
を踏まえると、$\hat{\beta}^{RIDGE}$の平均、バイアス、分散共分散行列が次のように求められます。

\begin{eqnarray}
\mathop{\rm E}(\hat{\beta}^{RIDGE}) &=& \bigl( I – λ(X^TX+λI)^{-1} \bigr) \beta^{*}\\
\mathop{\rm bias} ( \hat{\beta}^{RIDGE} )&=&\mathop{\rm E}(\hat{\beta}^{RIDGE})-b^{*}=- λ(X^TX+λI)^{-1} \beta^{*}\\
\mathop{\rm cov}(\hat{\beta}^{RIDGE}) &=& \sigma^2 \bigl( (X^TX+λI)^{-1}X^TX(X^TX+λI)^{-1}\bigr)
\end{eqnarray}

λ→0の時と、λ→∞の時を考えると、λの調整によってバイアスと分散がチューニングされることがイメージしやすいと思います。

λ→0の時(最小二乗法)

  • bias→$0$、cov→$σ^2(X^TX)^{-1}$

λ→∞の時(正則化を利かせまくったRIDGE回帰)

  • bias→$-b^{*}$、cov→$O$

厳密な証明は省きますが、λ→∞の場合は、$X^TX+λI$→$λI$とみなせるので、上記のような変形になります。

λ→0の時にバイアス(bias)がゼロに、λ→∞の時に分散(cov)がゼロ行列に収束すると言うのがポイントだと思います。

上記は極端な例※でしたが、正則化手法では、λ2を0と∞の間の適当な値にチューニングすることで、バイアスと分散のバランスを良い感じに調整しているという感覚が掴めたのではないかと思います。

※実際、λ→∞の場合、$\mathop{\rm E}(\hat{\beta}^{RIDGE})$がゼロになってしまうので、意味ある推定ができているとは言いにくいのですが、雰囲気が伝わればと思い、敢えて極端なケースを考えています笑。

情報量規準の観点

情報量規準の観点では、「正則化とは、情報量基準の最小化問題における近似解を得る手法」と言えるかもしれません。

情報量規準とは、推定した確率分布が真の確率分布とどれらい近いのかを測る尺度で、AIC、BICなどが代表的です。これらは小さければ小さいほど真の確率分布に近いこと(≒汎化性能が良いこと)を示します。

  • $\mathop{\rm AIC}=-2\log P(y|\beta)+2k$
  • $\mathop{\rm BIC}=-2\log P(y|\beta)+\log(N)k$

ここで、$P(y|\beta)$は先ほども登場した尤度、$k$はモデルに含まれる非ゼロのパラメータの個数です。$\beta$のみをパラメータとみなすなら(後で補足します)、$k=\|\beta\|_0$となります。

最尤法では尤度を最大化することを良しとしていましたが、情報量規準では、説明変数の追加によって尤度が改善されても、説明変数の個数増加が尤度の改善度合いよりも大きい場合、汎化性能が悪化したと判定します。つまり、「説明変数が多くなればなるほど、尤度は良くなるが、汎化性能が悪くなる」、という、「過学習問題に関するトレードオフ」が考慮されています。

$\sigma^2$が$\beta$に依存しない定数だとすると、上記のAIC、BICは、今回の問題設定では、次のように書き下せます。

  • $\mathop{\rm AIC}=\frac{1}{\sigma^2}(y-X\beta)^T(y-X\beta)+2\|\beta\|_0+const. $
  • $\mathop{\rm BIC}=\frac{1}{\sigma^2}(y-X\beta)^T(y-X\beta)+\log(N)\|\beta\|_0+const.$

$L$-$0$ノルムにおける正則化問題(式\eqref{4}に$q=0$を代入したもの)と同様の形式となります。つまり、$L$-$0$ノルムを用いた正則化では、過学習も加味した指標(すなわち情報量規準)を最適化しているとみなすことができる、ということです。

$L$-$0$ノルムを使った正則化手法は、とても理想的に思えるのですが、2つ問題があるように思われます。

1つ目の問題:
「$\sigma^2$が$\beta$に依存しない定数」とする条件が合理的だと思えないことです。実際、「最尤法との関係」の節で$\sigma^2$の推定値を導出しましたが、推定値は$\beta$の値に依存していました。

AICを例に具体的に見ていきます。導出過程は省略しますが、AICの定義において$\sigma^2$を$\beta$の関数で置き換えると、AICは次のような形になります。

$$\mathop{\rm AIC}=N\log\bigl(\frac{2\pi}{N}(y-X\beta)^T(y-X\beta)\bigr)+N+2(|\beta|_0+2)$$

$L$-$0$の正則化問題とは少し形状が変わってしまいましたね。したがって、$L$-$0$の正則化は、情報量規準(例えばAIC、BIC)の最適化問題を解いているとは厳密には言えず、あくまで、”近似解のようなもの”を求めている、と解釈するのが良いと考えています。

なお、AICの定義と上式の各項目の対応関係(AIC定義→上式の項目)は以下の通りです。

  • $-2\log P(y|\beta)$→$N\log\bigl(\frac{2\pi}{N}(y-X\beta)^T(y-X\beta)\bigr)+N$
  • $2k$→$2(|\beta|_0+2)$

$2k$が、$2(\|\beta\|_0)$ではなく、$2(\|\beta\|_0+2)$に書き換えられた理由は、$\sigma^2$を推定すべきパラメータと見做していることと、今回触れていない定数項$\beta_0$の存在にあります。導出過程や、定数項の扱いについて知りたい方は文末で紹介している参考書籍に当たってみると良いかもしれません。

2つ目の問題:
最適解を求める計算が面倒なことです。$L$-$0$の代替として、$L$-$1$や$L$-$2$を用いることにより計算上の問題を緩和することができます。この点については、次章「なんでLASSOなのか?」で触れたいと思います。

以上が、正則化が情報量規準の最適化問題の近似解っぽいものを求めている説明となります。

なんでLASSOを使うの?

スパース性と数学的性質の間で

正則化では、$L$-$1$ノルム(LASSO)または$L$-$1$ノルムと$L$-$2$ノルム(RIDGE)のハイブリッドなど、$L$-$1$ノルムを中心とした構成がなされることが多い気がします。それは、$L$-$1$ノルムがスパース性(≒変数選択ができる)と数学的性質の良さ(≒計算しやすい)の両方を兼ね備えた「できる子」だからだと思われます。

スパース性:
$\beta$の推定値にゼロが多く含まれる性質のことを、スパース性と言います。スパース性を持った正則化手法の代表格は$L$-$0$ノルムを用いたもので、非ゼロの要素の個数を含む損失関数を最適化するため、得られる解にはゼロが多く含まれることになります。$L$-$1$ノルム(LASSO)にも同様の性質があります。一方、$L$-$2$( RIDGE)では駄目です。このことを次節「LASSOの幾何学的解釈」で説明します。

数学的性質:
数理最適化において、凸関数の最適化は扱いやすいとされます。凸関数であれば、極地は一つに定まるので、極所解にハマることを考えずに済む、ということが大きい気がします。ここでは詳細を説明しませんが、$L$-$1$、$L$-$2$は凸関数ですが、$L$-$0$は凸関数ではありません。

したがって、スパース性と数学的な性質の良さ(凸関数)の両方の性質を兼ね備えたのは、$L$-$1$(LASSO)のみとなります。

今更ですが、LASSOとは、Least Absolute Shrinkage and Selection Operatorの略称です。つまり、正則化(Shrinkage)と同時に変数選択(Selection)までやってしまう「できる子」なのです。

LASSOの幾何学的解釈

前述した、「LASSOではスパースな解が得られるが、RIDGEではだめな理由」を、図形的な解釈により説明します。図形的な解釈に進む前に、元の問題設定(\eqref{4})を少しいじります。

\begin{eqnarray}
\hat{\beta}^{Reg}&=&\mathop{\rm arg~min}\limits_\beta (y-X\beta)^T(y-X\beta)+λ{|\beta|_q}^q\\
&=&\mathop{\rm arg~min}\limits_\beta (y-X\beta)^T(y-X\beta), s.t. {|\beta|_q}^q \leq t
\end{eqnarray}

※二番目の式への変形には、ラグランジュの未定定数法、及び、KKT条件を用いますが、ここでは説明を割愛します。詳細はこちらのブログが参考になるかもしれません。私は上部しか理解してないですが、こいつ(ラグランジュの未定定数法)は、サポートベクターマシンでも登場するので、アルゴリズムを深く理解するなら、ちゃんと理解しておかないといけないんでしょうね。。

つまり、${|\beta|_q}^q \leq t$で示される領域内で、$(y-X\beta)^T(y-X\beta)$を最小化すれば良いということになります。もっと言えば、$(y-X\beta)^T(y-X\beta)$の等高線と、${|\beta|_q}^q \leq t$で示される領域とが接するポイントが最適解が得られるポイントとなります。

汚らしい図で恐縮ですが、下記をイメージしていただければと思います。

ポイントは以下の通りです。

  • i) $q=1$(LASSO)の場合:
    制約条件が示す領域が四角なので、角($\beta_1$, $\beta_2$軸上の点)で接しやすく、ゼロを含んだ解が得られやすい
  • ii) $q=2$(RIDGE)の場合:
    制約条件が示す領域が円なので軸上の点で接することはほとんどなく、ゼロを含んだ解は得られにくい

ここでは$q=1, 2$の場合のみを扱いましたが、一般的な$L$-$q$ノルムの示す領域の形状に関しては、こちらのブログをご覧いただくとイメージしやすいかもしれません。

どうやってLASSOの解を求めるの?

LASSOで用いる$L$-$1$は数学的性質が良いとは言ったものの、角張ったところで微分が不連続になるため、解析的に解を求めることができません。ここでは、よくお見かけする近接勾配法の概要に触れらればと思います。

近接勾配法を用いた解法

近接勾配法は、滑らかな凸関数$f$と滑らかとは限らない凸関数$g$の和($f(x)+g(x)$)を最小化する$x$を求めるの手法です。

LASSOの文脈ですと、以下のような対応を考えれば良いです。

  • $f(\beta)=(y-X\beta)^T(y-X\beta)$
  • $g(\beta)=λ{\|\beta\|_0}$

アルゴリズムとしては、次の式で$t$ステップ目の$x_t$を更新していき、解が収束すればそれを最適解の近似値とします。

\begin{eqnarray}
x_{t+1}={\rm prox}_{\eta g}(x_t – \eta \nabla f(x_t))\tag{6}\label{6}
\end{eqnarray}

ここで、式中に登場する記号について補足します。

  • ${\rm prox}_{g}(y)$はproximal operator(近接写像)と呼ばれる関数で、次式で定義されます

\begin{eqnarray}
\mathop{\rm prox_{g}} (y) =\mathop{\rm arg~min}_{x}\bigl( g(x)+\frac{1}{2}\|x-y\|_2^2 \bigr) \tag{7}\label{7}
\end{eqnarray}

  • $\eta$は更新の幅を決めるパラメータで何らかの方法で設定します

LASSOにおける近接勾配法に関しては、こちらのQiitaの記事に詳細にまとめてくださっているので、詳細が気になる方、実装してみたい方、はぜひご覧ください。

あえて、この記事で概要を紹介したのは、私自身が、LASSOやその拡張手法を理解・実装しようと思った際に、偏微分地獄に陥り、偏微分の大事さを思い知ったためです。

次節「偏微分地獄」で、LASSOにおいて、偏微分がどのように使われているのか、雰囲気をお伝えできればと思います。

偏微分地獄

LASSOにおける近接勾配法で偏微分が登場するのは、式\eqref{6}における$\nabla f(x_t)$と式\eqref{7}における$\mathop{\rm arg~min}$です。

LASSOの文脈で、それぞれ書き下して見ましょう。

式\eqref{6}における$\nabla f(x_t)$:

$$\nabla f(\beta)=\frac{\partial}{\partial \beta} \bigl( (y-X\beta)^T(y-X\beta) \bigr)=2X^TX\beta-2X^Ty$$

※これまで、$\beta$を変数と見てきたので、$\nabla f(x_t)$→$\nabla f(\beta)$と読み替えています。
※ナブラ記号を$\nabla = (\frac{\partial}{\partial \beta_1}, \frac{\partial}{\partial \beta_2}, \cdots) = \frac{\partial}{\partial \beta}$と表記しています。

式\eqref{7}における$\rm arg~min$:

最小化問題ですので、例のごとく、偏微分をとってゼロとなるものを最適解とします。近接勾配法では、偏微分=ゼロの解が簡単に書き下せることを前提としているので、それに該当しない場合には別の手法を検討する必要がありそうです。もちろん、LASSOの場合には、偏微分=ゼロの解が簡単に書き下せます!

天下り式ですが、次のように展開されます。

まずは、偏微分ですが、$g(\beta)=λ{\|\beta\|_0}$を踏まえると、次の形式となります。

$$
\frac{\partial}{\partial \beta} \bigl( g(\beta)+\frac{1}{2}\|\beta-y\|_2^2 \bigr) =
\begin{cases}
\lambda + \beta_i – y_i & \beta_i>0\ \text{のとき}\\
-\lambda – \beta_i – y_i & \beta_i<0\ \text{のとき}\\
\end{cases}
$$

上式=0を$\beta$について解いて、prox関数を下記のように得ることができます。

$$
\bigl( \mathop{\rm prox_{\lambda\| \cdot \|_1 }}(y) \bigr)_i = \begin{cases}
y_i – \lambda & y_i > \lambda \ \text{のとき} \\
0 & |y_i| \leq \lambda \ \text{のとき} \\
y_i +\lambda & y_i < -\lambda \ \text{のとき}
\end{cases}
$$

次章「LASSOの弱点と発展モデル」でも一部ご紹介しますが、LASSOやLASSOの発展手法、あるいは、正則化手法において、正則化項の設定方法にかなりのバリエーションがあり、例えば新しい手法(今回のLASSOのようにQiitaで優しく解説されてないやつ笑)や、自分が考えた手法を試そうと思うと、どうしても自分で偏微分を頑張らねばなりません。

偏微分のための数学Tips

随所に出てくる偏微分を理解するために必要になった数学Tipsを紹介します。
偏微分マスターの方は読み飛ばしていただいてかまいせん。

例えば、LASSO界隈を探索していると、しばしばこんな微分に出会します。ここではこんな微分をやっつけるためのTipsを紹介しようと思います。

  • $X\beta$を$\beta$で微分する
  • $\beta^TX^TX\beta$を$\beta$で微分する
  • $\mathop{\rm Tr}(B^TX^TXB)$を$B$で偏微分する
    ※$B$は目的変数が多次元の場合に対応した回帰係数行列です。

最低限のTipsたち

ここでは偏微分はこなすために必要と思われる最低限のTipsを3つご紹介します。

  • アインシュタインの縮約記法:
    • なにこれ?:
      ①和の記号(Σ)を省略します(例:$Σ_i a_ib_i$ → $a_ib_i$)
      ②偏微分をカンマで表現します(例:$\partial b_i / \partial b_l$ = $b_{i,l}$)
    • なんで必要?:
      ①行列・ベクトル計算をシンプルに表現できます
      例:${AB}_{ij} = Σ_k a_{ik}b_{kj}$ → $a_{ik}b_{kj}$
      ②偏微分記号もシンプルに表記できます
      例:$\partial Ab / \partial b_l $ →$ A_{ij}b_{j, l}$
  • の微分:
    • なにこれ?
      掛け算を微分するための公式です。
      公式:$[g(x)f(x)]’=g'(x)f(x)+g(x)f'(x)$
    • なんで必要?
      行列積やベクトルとの内積を微分するのに役立ちます。
      例:$\partial b_iA_{ij}b_{j} / \partial b_l = b_{i,l}A_{ij}b_{j} + b_{i}A_{ij}b_{j,l}$
  • クロネッカーのデルタ$δ_{ij}$:
    • なにこれ?
      $i$と$j$が一致する場合に1を返し、それ以外は0を返す記号です。
      定義:$δ_{ij}=1(i=j), 0(i \neq j)$
    • なんで必要?
      随所に現れる$b_{j,l}$を$δ_{ij}$に置き換えて式の見通しを良くします。
      例:$A_{ij}b_{j,l} = A_{ij}δ_{jl}=A_{il}$
       ※$A_{il}$成分以外は$δ_{ij}=0(i \neq j)$の影響で消滅

Tips活用例

ご紹介したTipsを使って「$\mathop{\rm Tr}(B^TX^TXB)$を$B$で偏微分する」をやっつけたいと思います。

STEP1:
まず、$B^TX^TXB$の$i,j$成分をアインシュタインの縮約記法で表現します。

$$\{B^TX^TXB\}_{ij} = B_{ik}X_{lk}X_{lm}B_{mj}$$

ちょっと複雑なんですが、行列やベクトルの内積をアインシュタインの縮約記法で表記する際には、

  • $\{ABCD$…$Z\}_{ij}=A_{ia}B_{ab}C_{bc}D_{cd}$…$Y_{xy}Z_{yj}$のように添字がしりとりのように連関すること
  • 転置行列では、${A^T}_{ij}=A_{ji}$のように添字が前後で入れ替わること

をイメージすると、簡単に書き下せるようになります。

STEP2:
$\mathop{\rm Tr}(B^TX^TXB)$をアインシュタインの縮約記法で表現します。

一般に、行列A(正方行列)のトレースはアインシュタインの縮約記法を使うと、$\mathop{\rm Tr}(A)=\sum_i A_{ii}=A_{ii}$とかけます。このことと、STEP1の結果を使うと、

$$\mathop{\rm Tr}(B^TX^TXB)={ B^TX^TXB }_{ii} = B_{ki}X_{lk}X_{lm}B_{mi}$$

となります。$j=i$としただけです。

STEP3:
積の微分公式を使って$\mathop{\rm Tr}(B^TX^TXB)$を$B_{pq}$で偏微分する。

$$\frac{\partial}{\partial B_{pq}} \mathop{\rm Tr}(B^TX^TXB)=(B_{ki}X_{lk}X_{lm}B_{mi})_{,pq}=B_{ki,pq}X_{lk}X_{lm}B_{mi}+B_{ki}X_{lk}X_{lm}B_{mi,pq}$$

※ここでは、アインシュタインの縮約少し拡張して、$B_{pq}$による偏微分を$,_{pq}$と表現しています。

STEP4:
クロネッカーデルタを使って式の見通しをよくする。

\begin{eqnarray}
\frac{\partial}{\partial B_{pq}} \mathop{\rm Tr}(B^TX^TXB)&=&
B_{ki,pq}X_{lk}X_{lm}B_{mi}+B_{ki}X_{lk}X_{lm}B_{mi,pq}\\
&=&δ_{(k,i),(p,q)}X_{lk}X_{lm}B_{mi}+B_{ki}X_{lk}X_{lm}δ_{(m,i),(p,q)}\\
&=&X_{lp}X_{lm}B_{mq}+B_{kq}X_{lk}X_{lp}\\
&=&2X_{lp}X_{lm}B_{mq}
\end{eqnarray}

※ここでは、クロネッカーデルタを少し拡張して、次のような定義で用いています。

$$δ_{(k,i),(p,q)}=1((k,i)=(p,q)), 0((k,i)\neq(p,q))$$

平たく言うと、$(k,i)=(p,q)$以外の要素は消えると言うことなので、上述の式変形では$(k,i)=(p,q)$を代入し、$(k,i)$を消去してしまっています。

STEP5:
計算結果を行列の形式に戻します。

\begin{eqnarray}
\frac{\partial}{\partial B_{pq}} \mathop{\rm Tr}(B^TX^TXB)
&=&2X_{lp}X_{lm}B_{mq}\\
&=&\{2X^TXB\}_{pq}
\end{eqnarray}

よって、以下のように行列の微分が求められました。

$$\frac{\partial}{\partial B} \mathop{\rm Tr}(B^TX^TXB)
=2X^TXB$$

行列に纏わる各種微分に関しては、こちらのQiitaの記事も参考になりました。

LASSOの弱点を克服する発展モデルたち

「できる子LASSO」にも弱点があります。ここでは弱点とその対処法として生まれた発展モデルを紹介できればと思います。他にもっとあるだろとツッコミを受けそうですが、何卒ご勘弁を。

発展モデルとLASSOの弱点の対応関係:

  • elastic-net
    • この手法が解決するLASSOの弱点:
      多重共線性への弱点に対応するために生まれた手法っぽいです。弱点については、ホクソエムさんのブログにもう少し具体的に書かれています。
      • ホクソエムさんのブログからの引用:
        「相関の高い変数群がある場合、Lassoはその中の1つしか変数として選択できない。中略。Grouping効果を考慮できないということを意味しています。遺伝子データのような場合、事前知識から似ている特徴量が複数あり、しかもそれらを変数として選択したい状況もあるかと思います。このような場合、Lasso回帰では選択される変数が似ている特徴量の中のどれか一つとなってしまいます。そうすると、モデル作成者が達成したい目的はLasso回帰では達成できなくなります。」
    • この手法の特徴:
      中身はとっても単純で、L1とL2のハイブリッドの正則化項を追加したのみです。
    • elastic-netの論文はこちら
  • group-lasso:
    • この手法が解決するLASSOの弱点:
      one-hotエンコーディングしたカテゴリ変数など、特徴量をセットで選択したい場合があると思います。lassoだとセットでの選択には対応していません。group-lassoはセット選択を可能にする手法です。
    • この手法の特徴:
      正則化項として、特徴量セット毎のL2ノルム(二乗はしません)を用いています。
  • adaptive-lasso:
    • この手法が解決するLASSOの弱点:
      一致性が得られないという弱点に対応するための手法です。一致性とは、サンプル数が十分大きい時に、推定値が真の値に”近く”ことをいうそうです。
    • この手法の特徴:
      正則化項への重み付けを係数毎に使い分けることにより、一致性が担保できるそうです。係数の重みにはRIDGE回帰などで予め推定した係数の推定値の逆数(つまり小さい係数ほど正則化を強くする)を用いるようです。ホクソエムさんのブログにてRでの実装例が載っています。
    • adaptive-lassoの論文はこちら

最後に

最後まで、Myトリビアの供養にお付き合いいただきありがとうございます。
皆様の「ヘぇ〜」の数だけ大往生できる確率が高まったことでしょう。

各所で参考になりそうなWebや論文のリンクを紹介しましたが、今回触れた内容のベースには小西先生の多変量解析入門(下記書籍)がありますので、ご紹介させていただきます。だいぶ雑に説明した最尤法やAIC、BICの導出が丁寧に説明されていますので、気になった方はこちらにあたって見ると良いかもしれません。

その外、ロジスティック回帰、サポートベクターマシン、主成分分析、クラスター分析、にも丁寧に触れられていますので、私のトリビアとは違い、大変ありがたい内容になっているかと思います。

コメントを残す