Lassoがスパース推定になる理由について

はじめに

会社の輪読会でふと「lassoがスパース推定になる理由って結局何なの？」という話が出ました. よくそれを説明する際には以下のひし形を使った図が使われますよね.
[https://ja.wikipedia.org/wiki/ラッソ回帰#/media/ファイル:L1_and_L2_balls.svg:embed:cite]
輪読会でもそれに触れられたんですが， どうもしっくりきませんでした.

そこでふと思ったのですが「$L^2$正則化であるRidgeと比べてLassoは数式レベルで理解したことないなあ．なんか悔しいなあ」と帰り道に思っていました.

そこで, 週末にいくつか参考書を読んで, 自分なりにまとめたので公開しようと思います. 読んで参考になった本は以下の2冊です.
[https://www.amazon.co.jp/統計的機械学習の数理100問-Python-機械学習の数理100問シリーズ-鈴木-讓/dp/432012507X:title]

[https://www.amazon.co.jp/機械学習のための連続最適化-機械学習プロフェッショナルシリーズ-金森-敬文/dp/406152920X/ref=sr_1_1?__mk_ja_JP=カタカナ&keywords=連続最適化&qid=1637485102&s=books&sr=1-1:title]

本記事の議論は鈴木先生の参考書の5章を大いに参考にしています. 丁寧にまとめられていたので, 興味のある方は是非本屋さんで立ち読みして購入を検討してみて下さい. さて, まずはLasso推定がスパース推定になっている理由(本質的な概念)について先に触れておきます.

そもそも何が影響でスパース推定になっているのか？

Lassoがスパース推定になっている理由は「微分不可能な点の存在が本質的である」ということだと数式を追っていて思いました. 実際, (線形回帰)lasso推定は以下の式を$\bm{\beta}$について最小化する問題として定式化されます.

$$L(\bm{\beta}) = \| \bm{y} - X\bm{\beta}) \|_2^2 + \lambda \|{\bm{\beta}}\|_1 \ (\lambda > 0) \ \ \ \ \ \ (1)$$

最小化するので, Ridge回帰と同様に(1)式を$\bm{\beta}$で微分します. しかし, $\|\bm{\beta} \|_1 = \sum_{i=j}^p | \beta_j |$は各$\beta_j = 0$で微分不可能です. これでは解析する上では不便なので, 微分の概念を少し広げたものを代用するのが最適化の分野ではよく行われます. この概念は劣微分(subdifferential)と呼ばれるものです. 本記事では, まず劣微分の概念を導入します. その後, 劣微分を用いてlasso推定の数学的な解析を行い, lasso推定がスパース推定になっていることを確かめるという手順で進んでいきます.

劣微分(Subdifferential)

以下では簡単なため$x \in \mathbb{R}$の場合を考えることにします.

定義(凸関数(convex function))

$S \subset \mathbb{R}, f: S \to \mathbb{R}$とする. このとき, 関数$f$が集合$S$上の凸関数(convex function)
$\overset{\text{def}}{\Longleftrightarrow}$

$$\forall x, \forall y \in S, 0 < \alpha < 1, \ \ f(\alpha x + (1-\alpha) y) \leq \alpha f(x) + (1 - \alpha ) f(y).$$

凸関数の定義の主張は図に書くと非常に分かりやすいです. 2点$x, y \in S \subset \mathbb{R}$を結ぶ直線が$S$上で関数$f(x)$よりも上に書けることを要請するのが凸関数の定義です.

例1

$\forall x \in \mathbb{R}, f(x) = |x|$は凸関数.

例2

$\forall x \in \mathbb{R}, f(x) = x^2$は凸関数.
ここで, 劣微分の概念を導入します.

定義 (劣微分(subdifferential))

$f: \mathbb{R} \to \mathbb{R}$な凸関数とする. このとき, $\forall x, x_0 \in \mathbb{R}$について

$$f(x) \geq f(x_0) + \xi (x - x_0) \ \ \ \ \ \ (2)$$

を満たすような集合$\xi$を$f$の$x = x_0$の劣微分(subdifferential)といい, 記号$\partial f(x_0)$で表すことがある.

lasso推定を解析する上では, $f(x) = |x|$の$x = 0$に関する劣微分が重要である. そこで, 式(2)において$x_0 = 0$としたときの$|x| \geq \xi x$を満たす集合$\xi$について考えます.

補題

$$\forall x \in \mathbb{R}, |x| \geq \xi x \Longleftrightarrow |\xi| \leq 1.$$

つまり, 関数$f(x) = |x|$の$x = 0$の劣微分は閉区間$[-1, 1]$である.
Proof.
($\Longrightarrow$) $x > 0$のとき, 仮定より$x \geq \xi x \Longleftrightarrow 1 \geq \xi$. また, $x < 0$のときは $-x \geq \xi x \Longleftrightarrow z \geq -1$. よって, $|\xi| \leq 1$.
($\Longleftarrow$)

$$\xi x \leq |\xi x| \leq |\xi| \cdot |x| \leq |x|.$$

(証明終了)

Lasso推定

$X \in \mathbb{R}^{n \times p}$の行列であり, 標準化が行われているものとします. つまり各$\bm{x}_j = (x_{1j} x_{2j} \cdots x_{nj})^T$に対して,

$$\bm{x}_j = \frac{\bm{x}_j - \mathbb{E}[\bm{x}_j]}{\sqrt{\mathbb{V}[\bm{x}_j}]}$$

という処理がなされているものとする. また$\bm{y} \in \mathbb{R}^n$についても$\bm{y} = \bm{y} - \bm{1}\cdot \mathbb{E}[\bm{y}]$の中心化が行われているとします. それでは以上の用意のもとで

$$L(\bm{\beta}) = \frac{1}{2n} \| \bm{y} - X\bm{\beta}) \|_2^2 + \lambda \|{\bm{\beta}}\|_1 \ \ \ (\lambda \geq 0) \ \ \ \ \ \ (3)$$

の最小化問題を考えていきます.
最初は議論の見通しを良くするために

$$\frac{1}{n}\sum_{i = 1}^n x_{ij}x_{im} = \left\{ \begin{array}{cc} 1 & (m = j) \\ 0 & (m \neq j) \end{array} \right. \ \ \ \ \ \ (4)$$

を仮定します. また, $s_j := \frac{1}{n}\sum_{i=1}^n x_{ij}y_i$とおいておきます. まずは$L(\bm{\beta})$の$\beta_j$に関する微分を考えてみます. ここで, $\|\bm{\beta} \| = \sum_{j=1}^p |\beta_j|$より, $\beta_j = 0$に関して劣微分となるので

$$\frac{\partial L}{\partial \beta_j} = - \frac{1}{n}\sum_{i = 1}^n x_{ij}(y_i - \sum_{m=1}^p\beta_m x_{im}) + \lambda \left\{ \begin{array}{cc} 1 & (\beta_j > 0) \\ -1 & (\beta_j < 0) \\ \left[-1, 1\right] & (\beta_j = 0) \end{array} \right.$$

が得られる. また, 式(4)より$x_{ij} \cdot x_{ij}$の項しか残らないので

$$\begin{align*} \frac{\partial L}{\partial \beta_j} &= - \frac{1}{n}\sum_{i=1}^n y_i + \beta_j + \lambda \left\{ \begin{array}{cc} 1 & (\beta_j > 0) \\ -1 & (\beta_j < 0) \\ \left[-1, 1\right] & (\beta_j = 0) \end{array} \right. \\ &= \left \{ \begin{array}{cc} - s_j + \beta_j + \lambda & (\beta_j > 0) \\ - s_j + \beta_j - \lambda & (\beta_j < 0) \\ -s_j + \beta_j + \lambda \left[ -1, 1 \right] & (\beta_j = 0) \end{array} \right. \end{align*}$$

よって$\frac{\partial L}{\partial \beta_j} = 0$とおくと

$$\begin{equation*} \beta_j = \left\{ \begin{array}{cc} s_j - \lambda & (s_j > \lambda) \\ -s_j + \lambda & (s_j < - \lambda) \\ 0 & (- \lambda \leq s_j \leq \lambda) \end{array} \right. \end{equation*}$$

と書ける. よって $s_j = \sum_{i=1}^n x_{ij} y_i \in [-\lambda, \lambda]$であるとき, $\beta_j = 0$となる. これが, lasso推定がスパース推定になる理由の概要です. 今回は議論の見通しを良くするために，　(4)式を仮定しました. 以下では, この仮定を外した一般の場合についてのlasso推定を見ていきます．

Lasso推定:一般の場合

次は式$(4)$の仮定を外した一般の場合を考えてみる. この場合は$\beta_m(m \neq j)$を定数と思って各$\beta_j$について個別に更新するようなアルゴリズム更新を行う. この更新は連続凸最適化の文脈では近接勾配法(proxiaml gradient method)と呼ばれる手法の1つの特殊例になっているらしいです[金森 16].

さて, 実際に更新アルゴリズムを導いてみよう. $\beta_m(m \neq j)$を固定すると

$$\begin{align*} \frac{\partial L}{\partial \beta_j} &= - \frac{1}{n}\sum_{i=1}^n x_{ij}\Big( y_i - \sum_{m\neq j}^p \beta_m x_{im} - \beta_j x_{ij} \Big) + \lambda \left\{ \begin{array}{cc} 1 & (\beta_j > 0) \\ -1 & (\beta_j < 0) \\ \left[-1, 1\right] & (\beta_j = 0) \end{array} \right. \\ &= - \frac{1}{n}\sum_{i=1}^n x_{ij} (r_{ij} - \beta_j x_{ij}) + \lambda \left\{ \begin{array}{cc} 1 & (\beta_j > 0) \\ -1 & (\beta_j < 0) \\ \left[-1, 1\right] & (\beta_j = 0) \end{array} \right. \ (\because r_{ij} := y_i - \sum_{m \neq j}^p \beta_m x_{im} \text{とおいた}.) \end{align*}$$

となり, ここで$s_j := \frac{1}{n}\sum_{i=1}^n r_{ij}x_{ij}$とおくと

$$\begin{align*} \frac{\partial L}{\partial \beta_j} &= \frac{\beta_j}{n}\sum_{i=1}^n x_{ij}^2 - s_j + \lambda \left\{ \begin{array}{cc} 1 & (\beta_j > 0) \\ -1 & (\beta_j < 0) \\ \left[-1, 1\right] & (\beta_j = 0) \end{array} \right. \\ &= \frac{\beta_j}{n} - s_j + \lambda \left\{ \begin{array}{cc} 1 & (\beta_j > 0) \\ -1 & (\beta_j < 0) \\ \left[-1, 1\right] & (\beta_j = 0) \end{array} \right. \ \ (\because X\text{は標準化しているので} \sum_{i=1}^n x_{ij}^2 = 1) \end{align*}$$

が得られる. あとは

$$\begin{equation*} \beta_j = \left\{ \begin{array}{cc} s_j - \lambda & (s_j > \lambda) \\ -s_j + \lambda & (s_j < - \lambda) \\ 0 & (- \lambda \leq s_j \leq \lambda) \end{array} \right. \ \ \ \ \ \ (6) \end{equation*}$$

として$\beta_j$の更新を行う.

Lasso推定更新アルゴリズムの解釈

最後にlasso推定の更新則である式(6)の解釈を述べて終わりにしたいと思います. $\beta_j = 0$となるときは変数$\bm{x}_j$は$\bm{y}$を予測するのに重要でないことを表しています. 式(6)によると$s_j$が小さいとき$\beta_j = 0$となりやすいです. これは線形回帰でいう$\beta_j x_{ij}$項に対する残差$r_{ij} = y_i - \sum_{m \neq j}^p \beta_m x_{im}$が小さいことを表している. これは$y_i$を予測するには$\beta_j$以外の回帰係数を用いれば十分, つまり, $\bm{x}_j$以外の変数を用いれば予測には十分であることを示唆しています. よって, lasso推定更新アルゴリズム(6)式は妥当な理由でスパース化していることを表していることがわかります．

参考文献

[鈴木 20]鈴木譲, 統計的機械学習の数理100問 with Python,共立出版(2020)
[金森 16]金森敬文,鈴木大慈,竹内一郎,佐藤一誠, 機械学習のための連続最適化,講談社(2016)