ホーム » 「ラッソ」タグがついた投稿

タグアーカイブ: ラッソ

Lassoの最適化アルゴリズムの考察

こんにちは。

前回の最後に、
– Coordinate descent(座標降下法)
– ISTA (メジャライザー最小化)
– FISTA (高速化)
の3つの最適化方法を紹介しました。今回は上の2つを解説します。そして次回実装へうつりましょう。

まずは簡単なLassoからです。

Coordinate descent

いきなり解説に入る前に何ですが、座標降下法ってメチャクチャシンプルなんですよ。
けど全く有名じゃないんですよね。(多分、理由は記事の最後)まあ、アルゴリズムみましょう


目的関数f(x) = g(x) + \sum_{i=1}^{n} h_i(x)の最適化を考えます。
ここでgはconvexかつdifferentiabl、h_i(x)はconvexとします。
(例えばconvexかつundifferentiablなf(x)を考えた場合、反例があります。詳しくはREADMORE)
k=1, 2, \cdots回目の更新を右肩の添え字で表すと、各変数に対して

    \[x_1^{(k)} = \texttt{argmin } f(x_1, x_2^{(k-1)}, x_3^{(k-1)}, \cdots, x_n^{(k-1)})\]

    \[x_2^{(k)} = \texttt{argmin } f(x_1^{(k-1)}, x_2, x_3^{(k-1)}, \cdots, x_n^{(k-1)})\]

    \[x_3^{(k)} = \texttt{argmin } f(x_1^{(k-1)}, x_2^{(k-1)}, x_3, \cdots, x_n^{(k-1)})\]

のように更新するアルゴリズムです。


Pointは1つ、変数を1つずつ更新するんです。その他の変数は定数とみるんです。なので

この図の感じですね。しかし次の問題を考えてみましょうx \in \mathbb{R}^2に対して

    \[\texttt{Minimize } x_1^2 + x_2^2 + 5|x_1 + x_2|\]

です。例えば初期値をx=(1,-1)で座標降下法をアプライしてみると
x_1, x_2のどちらから更新してもx=(1,-1)で停滞することがわかります。しかし
これは明らかに最適解はx=(0,0)ですよね。こういったことからスパースな問題にはあまり良くないようです。

次にメジャライザーの方を行きましょう。

Iterative Shrinkage Thresholding

ちなみにこのアルゴリズムは初めて聞きました。EMアルゴリズムでも似たことが行われますね。
こちらもアルゴリズムは単純で

テイラー展開を行い、目的関数を置き換えることで最適化問題を解きやすい形にしよう

というものです。


一次元の場合、fx=a周りでのテイラー展開

    \[f(x) \approx f(a) + f^{(1)}(a)(x-a) + \frac{1}{2!} f^{(2)}(a)(x-a)^2 + \frac{1}{3!} f^{(3)}(a)(x-a)^3 + \cdots\]

多次元の場合、

    \[f(\vc{x})= f(x_1,x_2, \ldots, x_n)\]

    \[f(\bf{x}) \approx f(\bf{a}) + Df(\bf{a}) (\bf{x}-\bf{a}) + \frac{1}{2} (\bf{x}-\bf{a})^T Hf(\bf{a}) (\bf{x}-\bf{a})\]


下の図がイメージになります。上界というのは「より大きいやつ」という意味です。

この図からアルゴリズムがどうやって動いて何で最適化できるのかは直感的にわかったかと思います。
しかし、このメジャライザーが本当に存在するのか?など疑問に思う人がいるでしょう。
では考察しましょう。

まず初期値をx^{(t)}とします。そして目的関数を\mathcal{L}(x),メジャライザーを\mathcal{T}(x)とします。
そして、メジャライザーは次のような性質を満たしてほしいです。
\mathcal{L}(x^{(t)}) = \mathcal{T}(x^{(t)})
\forall x\neq x^{(t)} に対して \mathcal{L}(x) \leq \mathcal{T}(x)
すると性質から明らかにTの最小化がLの最小化に繋がることがわかります。
では実際にそんな都合のいい\mathcal{T}(x)を探しましょう。の前にちょっと寄り道します。

どんな形のメジャライザーが都合いい?

メジャライザーT(x)を探すといえどこのままでは全くわかりません。
T(x) = x^2にすればいいのか?いや、それもわかりません。なので

都合のいいメジャライザーの形を予想します。

いきなりですが次の最適解xを考えましょう。次元は1です。(\lambda \geq 0)

    \[\texttt{argmin } \frac{1}{2} (y-x)^2 + \lambda |x|\]

前回の記事でやったようにこれは場合分けで簡単にできますね
\forall x \geq 0に対して

    \[\frac{1}{2} (x^2 - 2xy + y^2 + 2\lambda x) = \frac{1}{2} (x - (y-\lambda))^2 + \texttt{const}\]

同様にして\forall x < 0に対して

    \[\frac{1}{2} (x - (y+\lambda))^2 + \texttt{const}\]

では最後にyの範囲で場合分けすると前回のSoft thresholding functionが出てきましたね。

    \[x = S ( y, \lambda )\]

詳しく復習すると
y < - \lambdaの時

    \[x= y + \lambda\]

- \lambda \leq y \leq \lambdaの時

    \[x = 0\]

y > \lambdaの時

    \[x = y - \lambda\]

これが大切です。

 つまり平方完成して解が求まりました。

ではLassoの目的を再掲しましょう。

    \[\texttt{argmin } ||y - Xb||_2^2 + \lambda|b|_1^1\]

今の例と似てますよね。

 しかしここままでは平方完成できない

ここで

 平方完成できるやつを考えよう!!!

これがほしいメジャライザーを予測するということです。

では元に戻って都合のいい\mathcal{T}(x)を探しましょう。平方完成できるやるがいいんでしたね。
テイラーの2次近似の形を想像しながら

    \[\mathcal{T}(\mathbf{x}) = \mathcal{L}(\mathbf{x}^{(t)}) + \frac{d\mathcal{L}(\mathbf{x}^{(t)})}{d\mathbf{x}}(\mathbf{x} - \mathbf{x}^{(t)}) + \frac{\rho}{2}\bigl\vert\mathbf{x}-\mathbf{x}^{(t)}\bigr\vert^2\]

という形で考えましょう、(\frac{\rho}{2}は計算の都合上のもの)そうすれば

    \[= \frac{\rho}{2} \Bigl\vert \bigl( \mathbf{x}^{(t)} - \frac{1}{\rho} \frac{d\mathcal{L}(\mathbf{x}^{(t)})}{d\mathbf{x}}\bigr) - \mathbf{x} \Bigr\vert^2 + \texttt{const}\]

とできますね。よって

    \[\text{argmin } \Bigl(\mathcal{T}(\mathbf{x}) + \lambda |\mathbf{x}|\Bigr)\]

    \[x^{(t+1)} = S \Bigl( \frac{1}{\rho} \frac{d\mathcal{L}(\mathbf{x}^{(t)})}{d\mathbf{x}}, \lambda \Bigr)\]

となります。

 実はまだ終わっていない

そうなんです。先ほどメジャライザーを考えたときに\rhoを使いましたね。
この\rhoを決める必要があります。メジャライザーの性質(2)から

    \[\mathcal{T}(\mathbf{x}) - \mathcal{L}(\bf{x}) = \frac{1}{2} (\bf{x} - \bf{x}^{(t)})^\top (\rho I - X^TX) (\bf{x} - \bf{x}^{(t)})\]

が「0以上」を任意の\bf{x}について言える必要があるので

 (\rho I - X^TX)が半正定値

の必要があります。つまり、X^TXの最大の固有値が\rhoとして取りうる最小の値になります。

そこで固有値の大きさを見積れる定理があるのでそれを使いましょう。
Gershgorin circle theorem
Gershgorin’s Theorem for Estimating Eigenvalues
この定理からmaxをとることで

    \[\texttt{eigenvalue } \le \max_j \sum_i\bigl\vert{\mathrm{X^\top X}}_{ij}\bigr\vert\]

と固有値の大きさが見積もれます。

長かったですがとりあえずLassoはおしまいです。

でわ。

READMORE

Lassoの進化、Group Lassoとは

こんにちは。

Lassoでスパース、スパースと言ってましたが実際はスパース推定という言葉がよく使われます。これについての説明を軽くしてからGroup lassoの紹介をします。

スパース推定とは

スパースは”疎”を意味します。スパース推定とはどのパラメータが0になるかを推定すること。つまり、データの本質がわかります。直感的には次の図から

もう少し実用的な例は?


右の入力から左の出力を考えるLassoモデルを考えると

のようにアレルゲン反応のスパース推定ができます。

つまり、スパース推定により関与しないパラメータがわかります

Group Lassoとは

グループラッソとは

「潰れる変数がグループになったLassoモデル」

先ほどのアレルゲンを例にとると例えば、ヒノキなどの個体単位ではなく「花粉」というグループ単位でスパース性が検証できるモデルのことです。Lassoはの正規化項は次のものでした。

    \[\Omega_{Lasso}(\beta)  =  \sum_{i=1}^{n} |\beta_i|\]

Group Lassoの正規化項は次で定義されます

    \[\Omega_{Group}(\beta)  =  \sum_{g=1}^{G}  \sqrt[]{p_g} |\beta_{g}|_2\]

ここでgはg番目のグループを表すindex。(ただし、g = 1,\cdots, Gp_gはグループgの大きさ)
前述通りGroup Lassoでは特徴をグループ化します。よって、事前に類似の傾向がありそうな特徴の情報を考慮します。

なるほど、と思った方と、ん?、と思った方がいると思います。

グループ単位で本当に潰れるのか?
仮にそうならどうやって変数を扱うのか?

これは図を用いてチェックしましょう。

Group Lassoの解

  • 参考文献の論文に従い、解説します。

まずJ個の変数からなるもっとも一般的な回帰問題を考えます

    \[Y = \sum_{j=1}^{J} X_j \beta_j + \epsilon\]

Yはn\times 1ベクトル、\epsilon \sim \mathcal{N}(0, \sigma^2I)X_jj番目のデータに対応したn\times p_j行列で\beta_jはサイズp_j \times 1: j=1,\cdots,Jの係数からなるベクトルとする。さらに各X_jは直交行列であると仮定する。
すなはちX_j^t X_j = I_{p_j}; j=1,\cdots ,Jとする。さらに、X = (X_1, X_2, \cdots, X_J)\beta = (\beta_1^t, \beta_2^t, \cdots, \beta_J^t)^tとすると上式はY=X\beta + \epsilonかける。


長ったらしく書きましたが要は「グループの大きさが各p_iのサイズ」です。なのでp_1 = p_2 = \cdots = p_J = 1を考えるとこれはLassoそのものとなります。


\eta \in \mathbb{R}^d, d \geq 1d \times dの正定値対称行列Kに対して次を定める。

    \[|\eta|_K = (\eta^t K \eta)^{1/2}\]

ただし、|\eta| = |\eta|_{I_d}とする。正定値行列K_1,\cdots,K_Jが与えられた時、Group Lasso回帰では次の解を考える。ただし\lambda \geq 0

    \[\frac{1}{2} | Y - \sum_{j=1}^{J} X_j \beta_j  |^2 + \lambda \sum_{j=1}^{J} |\beta_j|_{K_j}\]

Bakin(1999)はこれをグループ変数によるLassoの拡張版として提案しました。


ここも特に気にする必要はなく、大切なのは正定値対称行列Kにより変数に「重み」が掛かっているところです。機械学習ではこのように変数に重みを加える動作をよくします。一例として
マハラノビス距離(Mahalanobis distance)
があります。僕たちが無意識に距離として使っているものはユークリッドノルムでK=I、つまり単位行列の時です。

では先ほどの

「どうやってグループ単位で変数が潰れるのか?」

をみましょう。

下の図はグループが2つ(各係数はベクトルとスカラー)、つまり\beta_1 =(\beta_{11}, \beta_{12})^t, \beta_2の場合を考え、K_jが単位行列の場合の正規化項を表しています(ラッソはダイヤ、リッジは円だったやつの三次元バージョン)

つまり、

  • Figure(a)は   ||(\beta_1^t,\beta_2) ||_1^1   = 1
  • Figure(e)は   ||\beta_1||_K +  ||\beta_2||_1^1 = 1
  • Figure(i)は   ||(\beta_1^t,\beta_2) ||_2^2   = 1

もっと噛み砕くと

  • Figure(a)は   LassoのL_1-norm
  • Figure(e)は   GroupLassoのnorm
  • Figure(i)は   RidgeのL_2-norm

です。Figure(e)において\beta_{11}- \beta_{12}平面を見ると\beta_1というグループ単位でゼロに潰れることがわかります。

だからGroup Lassoはグループ単位ででスパース性を与えているんです。

でわ。

READMORE

Lassoを数式から実装まで(理論編)

こんにちは。

その1でラッソの概要と大きな特徴であるスパース性を確認しました。
今回からはラッソ実装に向けて数学を頑張りましょう。

Notationのについて

  • Obj (objective function)
  • OLS (Ordinary least squares) (回帰の残差平方和)
  • m: データ数
  • n: 次元(特徴)
  • \theta_j: j番目の特徴\theta
  • x_j^{(i)}: i番目の観測データxのうちjに対応するもの

目的関数について

今回、つまりラッソの目的関数は以下の通り

    \[Obj(\theta) & = OLS(\theta) + \lambda | \theta |_1^1\]

    \[= \frac{1}{2} \sum_{i=1}^m \left[y^{(i)} - \sum_{j=0}^n \theta_j  x_j^{(i)}\right]^2 + \lambda \sum_{j=0}^n |\theta_j|\]

(OLS)の勾配について

上の目的関数は微分できません。なのでひとまず置いておいて
OLSをいつも通り微分して微分を計算します。
後々のため、今回は勾配でなく、\theta_jでの微分を計算します。

    \[\frac{\partial }{\partial \theta_j} OLS(\theta) & = -  \sum_{i=1}^m x_j^{(i)}  \left[y^{(i)} - \sum_{j=0}^n \theta_j x_j^{(i)}\right]\]

    \[= -  \sum_{i=1}^m x_j^{(i)}  \left[y^{(i)} - \sum_{k \neq j}^n \theta_k x_k^{(i)} - \theta_j x_j^{(i)}\right]\]

    \[= -  \sum_{i=1}^m x_j^{(i)} \left[y^{(i)} - \sum_{k \neq j}^n \theta_k x_k^{(i)} \right] +  \theta_j \sum_{i=1}^m (x_j^{(i)})^2\]

簡単のためこれを次のように表します。

    \[:= - \rho_j + \theta_j z_j\]

L1-normの微分について

L_1-normは微分できないとわかりました。なので
微分できるように、微分という概念を拡張します。そこで出てくるのが

 劣勾配(Subgradient) 

です。wikiのリンク


定義(劣微分)
 Convex f: I \rightarrow \mathbb{R} の点 x_0 における劣微分は次の条件を満たす c \in \mathbb{R} で定義する。

(簡単のため\mathbb{R}^1についてはなし、微分の代わりに勾配という言葉を使っている)

    \[\frac{\partial f }{\partial  x}\bigg |<em>{x=x_0} = { c \in \mathbb{R} : f(x) \geq f(x_0) + c(x - x_0) , \forall x \in I  }\]

書き換えると
x_0における劣勾配は次の閉区間の集合である。a \leq b

    \[a = \lim</em>{x \rightarrow -x_0} \frac{f((x) - f(x_0)}{x - x_0}\]

    \[b = \lim_{x \rightarrow +x_0} \frac{f((x) - f(x_0)}{x - x_0}\]

    \[\frac{\partial f }{\partial  x}\bigg |_{x=x_0} =  [a,b]\]

これが勾配の拡張になっていることは微分可能な点においてその劣勾配が一点集合になることからわかる


実はこれは簡単で、場合分けして綺麗に微分できるところの間を1つの微分の集合と見るのです。
今回のL_1ノルム、つまり

    \[f(x) = |x|\]

を例にとると、

(1)   \begin{align*} \texttt{if }  x > 0  ~~ , ~~  \frac{\partial f}{\partial x} =  1 \end{align*}

(2)   \begin{align*} \texttt{if }  x < 0 ~~ , ~~  \frac{\partial f}{\partial x} =   -1 \end{align*}

そしてこの間の集合を0での微分とするのが劣勾配

(3)   \begin{align*} \texttt{if }  x = 0 ~~ , ~~  \frac{\partial f}{\partial x} =  [-1,1] \end{align*}

次のはイメージ図

ではこの劣勾配を使って今回の目的関数の第2項を微分してゼロとおきます。

(4)   \begin{align*} \frac{\partial }{\partial \theta_j}  Obj(\theta)  =     \frac{\partial }{\partial \theta_j} OLS(\theta) + \frac{\partial }{\partial \theta_j} \lambda | \theta |_1^1  \end{align*}

(5)   \begin{align*} 0 = -\rho_j + \theta_j z_j + \frac{\partial }{\partial \theta_j} \lambda | \theta_j |_1^1 \end{align*}

ここで先ほどと同様に場合分けを行って

(6)   \begin{align*} \texttt{if } \theta_j < 0  ~~ , ~~  0 = -\rho_j + \theta_j z_j  - \lambda  \end{align*}

(7)   \begin{align*} \texttt{if }  \theta_j > 0  ~~ , ~~  0 =  -\rho_j + \theta_j z_j +  \lambda  \end{align*}

そして、

(8)   \begin{align*} \texttt{if }  \theta_j = 0 ~~ , ~~  0 =   [-\rho_j  - \lambda ,-\rho_j + \lambda ]  \end{align*}

ですね。文字が多くなってしまったので再確認ですが、目標は\theta_jを求めることです。上の2つの式は簡単に求まりますが、3つ目の式については閉区間に0が含まれるという条件で考えましょう。なので

(9)   \begin{align*} 0 \in [-\rho_j  - \lambda ,-\rho_j + \lambda ] \end{align*}

(10)   \begin{align*} -\rho_j  - \lambda \leq 0 \end{align*}

(11)   \begin{align*} -\rho_j  + \lambda \geq 0 \end{align*}

なので

(12)   \begin{align*} - \lambda \leq \rho_j \leq \lambda \end{align*}

まとめると

(13)   \begin{align*} \texttt{for }  \rho_j < - \lambda  ~~ , ~~  \theta_j = \frac{\rho_j + \lambda}{z_j}  \end{align*}

(14)   \begin{align*} \texttt{for }  - \lambda \leq \rho_j \leq \lambda  ~~ , ~~  \theta_j = 0 \end{align*}

(15)   \begin{align*} \texttt{for }  \rho_j > \lambda   ~~ , ~~  \theta_j = \frac{\rho_j - \lambda}{z_j} \end{align*}

となります。ちなみにこの場合分けからSoft thresholding functionという関数が定義されます。

最後に劣微分と一緒にplotを見てみましょう。

微分できなくてよく使う関数をたまたま思い出しました。ニューラルネットワークで使うRELU fucntionですね。
活性化関数ReLUとその類似関数
他にも色々あるようです。

色々計算づくしでしたが微分を計算できました。あとは最適化のみですが今までは勾配法を使っていましたが「微分不可能な場合は勾配法は無理」なので他の最適化法が必要です。そこでLassoの問題ではいくつかのメジャーな方法があります。

  • Coordinate descent(座標降下法)
  • ISTA (メジャライザー最小化)
  • FISTA (高速化)

です。もっとも簡単なのは座標降下法です。メジャライザーはテイラー展開を用いて上界を扱うものですがこれらについては次回解説します。

でわ。

READMORE

Lassoはなんでスパース?

こんにちは。
素人にfish_shellは無理でした。kzです。

リッジが終わりましたね。ついに来ました

ラッソ

最近色々ラッソについて調べていたんですが、微分不可能な関数の最適化ってやっぱ難しいですね。機械学習において非常に重要な2つのキーワード
– ラグランジュの未定乗数法
– KKT条件
は別の記事でゆっくり解説します。では本題に入りましょう。

Lasso

過学習を考慮した回帰モデルの一つ
– L_1正則化項を使用した回帰model
スパース性を考えるときに用いる(これについては次の記事で詳しく説明します。)

(1)   \begin{align*} \beta = \texttt{argmin}_b { |y-Xb |^2_2 + \lambda|b|_1^1   } \end{align*}

リッジ回帰との唯一の違いは正規化項がL1(絶対値)であるということ。

微分できない?

ちょっと微分について復習しよう。おまけで複素解析も出てくるよ



L1とL2

これまでに何度か説明していると思いますがまずは
– 微分可能性
が大きな違いです。L1は尖っているので0で微分できませんね。

他の違いは?

リッジ回帰ではデカイパラメータbがでかくなりすぎないよう上限を設けましたね。
ラッソも上限はありますが、ゼロがKEYWORDです。なぜなら
– 無関係な特徴量はゼロで排除する
という特徴があります。ゆえに、モデルで初めに設定したものよりも少ない特徴量で済む可能性があります。これが先ほどのスパース性と言われる所以です。

パラメータが0を図から考察

では図を用いて直感的に理解しましょう。次の図はRidge, Lassoを議論するときに使われるとてもポピュラーなものです。

「スパース」とは「スカスカ」なイメージ

つまり、ゼロがいっぱいあるイメージでいい。
しかし重要なのは常にスパースなわけじゃないこと。

スパースと言われるのはこの「尖ったポイント」に最適解がある時。

これはL1だからできますよね。L2は円だったのでスパース性はありません。

パラメータが0を数式から考察

次のように超シンプルな回帰モデルを考えます。(以下転置使ってますが無くてもいいです)

(2)   \begin{align*} y = bx + \epsilon \end{align*}

次の最小化を考えます。

(3)   \begin{align*} \texttt{minimize  }  y^Ty - 2y^T xb + b x^Tx b + 2\lambda |b|_1^1 \end{align*}

(この2は計算の都合上のもの)

ここで解がb>0と仮定します。

すると、y = bx \Leftrightarrow <y,y> = b<y,x>より
y^Tx >0と仮定するのと同値であることがわかります。さて、目的関数をbに関して微分しましょう

仮定より今は微分できます

(4)   \begin{align*} \frac{\partial }{\partial b} \texttt{objective function} = -2y^Tx + 2x^Tx b + 2\lambda \end{align*}

であり、これをゼロと置くことでこの解は次のものだとわかります

(5)   \begin{align*} b = \frac{y^Tx - \lambda}{x^Tx } \end{align*}

ここで\lambdaを増やしていくと\lambda = y^Txでゼロになることがわかりますね。
この瞬間にLassoはスパース性を持ちます。

一方、リッジの場合は 

(6)   \begin{align*} \textrm{minimize  } &  y^Ty - 2y^T xb + b x^Tx b + 2\lambda |b|_2^2  \end{align*}

(7)   \begin{align*} b &= \frac{y^Tx}{x^Tx + \lambda} \end{align*}

となるので、b \neq 0よりパラメータはゼロにはならないですね。

でわ