正則化（L1, L2, ElasticNet, GroupLasso）とProximal Operator（近接勾配法）の実装

# はじめに

最近は今までライブラリに頼ってブラックボックス化していたものを、一つずつでもなくせるように色々と実装しています。前回はSGDによるロジスティック回帰を、今回は正則化・近接勾配法を実装したので備忘録程度にまとめました。

- [ロジスティック回帰で各種勾配降下法のアニメーション・可視化（python実装）](https://qiita.com/hokuto_HIRANO/items/ff845d2385798d1aed44)

# 近接勾配法

機械学習などでは目的関数を最小化する作業があります。そして、関数の最小化と言えば「微分」です。もし、目的関数が**微分可能**であれば勾配降下法（SGDなど）が使用できます。しかし、（部分的に）**微分できない**目的関数も存在します。そのような時は**近接勾配法**が利用できます。近接勾配法は式の導入こそややこしくはありますが、実際にプログラムするところは簡単で、勾配降下法の更新式に**Proximal Operator**という項を掛け合わせるだけです。

再度、ロジスティック回帰を例にあげながら説明します。

- **目的関数L=f(w)が微分可能であるとき**

ロジスティック回帰で言うところのF(w)は、以下のようになります。

```math
L = f(w) = - \sum_{n=1}^{N}\ln\{s(x_n)^{y_n}(1 - s(x_n))^{1-y_n}\}\\
s() : シグモイド関数\\
x : 説明変数\\
y : 正解ラベル（0, 1）
```

この場合の重みwの更新式は勾配降下法により以下のようになります。

```math
w \leftarrow w - \eta\frac{\partial L}{\partial w}
```

- **目的関数L=f(w)+g(w)が（部分的に）微分できないとき**

例えば、以下に示すようなL１正則化ロジスティック回帰の目的関数です。

```math
L = f(w) + g(w)\\
  = - \sum_{n=1}^{N}\ln\{s(x_n)^{y_n}(1 - s(x_n))^{1-y_n}\} + \lambda \| w \|_1\\

※ g(w) = \lambda \| w \|_1\\
```

**L１正則化の罰則項であるg(w)はw=0のところで連続ではないので微分できません**

このように、目的関数の中に（部分的に）微分できないg(w)のような項を含む場合は近接勾配法を利用します。

近接勾配法による重みwの更新式は以下のようになります。

```math
w \leftarrow prox_{\gamma g}(w - \eta\frac{\partial L}{\partial w})\\
※ prox_{\gamma g} : Proximal Operator
```

勾配降下法との違いは勾配に**Proximal Operatorという変数が掛け合わされているだけ**です。Proximal Operatorについては次章で紹介します。

ここまでで一旦、近接勾配法によるロジスティック回帰の枠組みをプログラムします。

```.py
# 近接勾配法
def proximal_gradient(grad_f, prox, gamma, objective, init_w, tol = 1e-6):
    w = init_w
    result = objective(w)
        
    while 1:
        
        w_new = prox(w - gamma * grad_f(w), gamma)
        result_new = objective(w_new) 
        
        # 停止条件を説明変数差分のノルムに
        if np.linalg.norm(w_new - w, 2) < tol:
            break;
                    
        w = w_new
        result = result_new
    return w_new, result_new


# X:学習用データベクトルを列ベクトルとして並べた行列
# y:ラベルを並べたベクトル
def LogisticRegression(X, y, lam):
    #　ロジスティック回帰内部の計算
    sigma = lambda a : 1.0/(1.0+np.exp(-a))
    p=lambda Z, w:sigma(np.dot(Z.T, w))
    X=np.concatenate((np.ones((1,np.size(X,1))), X), 0) #バイアス項の追加

    # 勾配
    grad_E = lambda w: np.dot(X, p(X, w) - y)
    (u, l, v) = np.linalg.svd(X)
    gamma = 1.0/max(l.real*l.real)
    
    # 目的関数
    objective = lambda w: - np.sum( y * np.log(p(X, w)) + \
                          (1 - y) *np.log(( 1 - p(X, w)))) + \
                          l1_norm(lam, w) # <= 正則化手法に応じて変更（次の章で紹介）
                                          # l1_norm(lam, w)
                                          # l2_norm(lam, w)
                                          # elastic_norm(lam, w)
                                          # group_norm(lam, w)
                        
    # proximal operator
    prox = lambda w, gamma:prox_l1(w, gamma, lam) # <= 正則化手法に応じて変更（次の章で紹介）
                                                       # prox_l1(w, gamma, lam)
                                                       # prox_l2(w, gamma, lam)
                                                       # prox_elastic(w, gamma, lam)
                                                       # prox_group(w, gamma, lam)
    
    # 近接勾配法
    w_init = np.zeros(X.shape[0])
    (w_hat, result) = proximal_gradient(grad_E, prox, gamma, objective, w_init, 1e-6)
    return w_hat, result
```


# 各種Proximal Operator

Proximal Operatorは正則化項のようなg(x)を使って計算されるパラメータで、式による導入が可能です（[詳細](http://yamagensakam.hatenablog.com/entry/2018/02/14/075106)）。以下に、一般的に使用される正則化項g(w)とそのProximal Operatorを紹介します。加えて、それらを計算する関数も載せています。

## L１正則化

L１正則化は意味のない説明変数の重みwがゼロになりやすいという特徴を持ち、過学習を抑える。よって、「不要なパラメータを削りたい」・「変数選択したい」という時によく使われる。なお、L1正則化はLassoとも呼ばれる。

```math
g(w) = \lambda \| w \|_1\\

prox_{\gamma g}(w) = S_{\gamma \lambda}(w) = \left\{
\begin{array}{ll}
w_i - \gamma \lambda & (w_i \geq \gamma \lambda) \\
0 & (- \gamma \lambda \lt 0 \lt \gamma \lambda) \\
w_i + \gamma \lambda & (w_i \leq \gamma \lambda)
\end{array}
\right.
```

```.py
def l1_norm(lam, w):
    penalty = lam * np.sum(np.abs(w))
    return penalty

def prox_l1(w, gamma, lam):
    prox = np.zeros(w.shape[0])
    thresh = gamma * lam
    prox[w >= thresh] = w[w >= thresh] - thresh
    prox[w <= -thresh] = w[w <= -thresh] + thresh
    
    return prox
```

## L２正則化

L2正則化はL1正則化ほど極端ではないが（重みwがゼロになりやすいわけではないが）、過学習を抑えて汎化された滑らかなモデルを得やすくなる。なお、L2正則化はRidgeとも呼ばれる。

```math
g(w) = \lambda \| w \|_2 \\

prox_{\gamma g}(w) = \frac{1}{1 + 2 \gamma \lambda}w
```

```.py
def l2_norm(lam, w):
    penalty = lam * np.sum(np.abs(w**2))
    return penalty

def prox_l2(w, gamma, lam):
    prox = w/(1+2*gamma*lam)
    return prox
```

## Elastic Net

Elastic NetはL1正則化とL2正則化を組み合わせた正則化手法です。

```math
g(w) = \lambda_1 \| w \|_1 + \lambda_2 \| w \|_2 \\

prox_{\gamma g}(w) = \frac{1}{1 + 2 \gamma \lambda_2}S_{\gamma \lambda_1}(w)
```

```.py
def elastic_norm(lam, w):
    lam1 = lam[0]
    lam2 = lam[1]
    penalty = l1_norm(lam1, w) + l2_norm(lam2, w)
    return penalty

def prox_elastic(w, gamma, lam):
    lam1 = lam[0]
    lam2 = lam[1]
    Sv = prox_l1(w, gamma, lam1)
    prox = Sv/(1+2*gamma*lam2)
    return prox
```

## Group Lasso（重複あり）

Group Lassoは説明変数がいくつかのグループに分けられる時に、グループごとに正則化を行う手法です（重みベクトルwのうちグループgに属する変数を並べたのがw<sub>g</sub>）。これにより、グループごとにスパースになる・ならないが決まります。そして、以下は**グループ間で説明変数の重複がない**場合です。

```math
g(w) = \lambda \sum_{g \in G} \| w_g \|_2 \\

prox_{\gamma g}(w) = \left\{
\begin{array}{ll}
w_g - \gamma \lambda \frac{w_g}{\| w_g \|_2} & (\| w_g \|_2 \geq \gamma \lambda) \\
0 & (\| w_g \|_2 \leq \gamma \lambda)
\end{array}
\right.
```

```.py
def group_norm(lam, w):
    penalty = 0

    for i in group_list:
        group_index = np.where(group_ids == i)
        penalty += np.sum(np.abs(w[group_index]**2))

    return lam * penalty

def prox_group(w, gamma, lam):
    # group_ids と group_list は事前に定義して、
    # グローバル変数として使っています。
    prox = np.zeros(w.shape[0])
    thresh = gamma * lam
    
    for i in group_list:
        group_index = np.where(group_ids == i)
        group_l2 = np.sum(np.abs(w[group_index]**2))
        
        if group_l2 > thresh:
            prox[group_index] = w[group_index] - gamma*lam*w[group_index]/group_l2
    
    return prox
```

ちなみに、グループ間で重複がある場合のGroup Lassoも提案されているようですがまだ追えていません。だれか、Proximal Operatorを求めていないでしょうか？？？

# 実験結果

糖尿病患者のデータセットを用いて糖尿病かどうかを学習する際の重みWの挙動を可視化した。

- [Regularization-ProximalOperator](https://github.com/hokuto-HIRANO/Proximal-Gradient-Method) : 一連のプログラム

---

**L１正則化**

重みが完全に０になる。

![l1-150.gif](https://qiita-image-store.s3.amazonaws.com/0/247972/ae0dc36d-bf43-f093-cfcd-d2627666c9ec.gif)

---

**L２正則化**

L１ほど重みが０になりにくい。

![l2-150.gif](https://qiita-image-store.s3.amazonaws.com/0/247972/3a073091-967b-a713-33b6-4d7b3b3efbeb.gif)

---

**Elastic Net**

- L１ラムダとL２ラムダが同じ割合　

![elastic-same-50.gif](https://qiita-image-store.s3.amazonaws.com/0/247972/d6da66a3-f2bc-e92b-7e64-8283c2762a75.gif)

- L１ラムダがL２ラムダの２倍　

L１の方が強く効くので重みがゼロになりやすい。

![elastic-l1x2-50.gif](https://qiita-image-store.s3.amazonaws.com/0/247972/df9f870f-59ef-74d9-3c99-87a8de21d067.gif)

- L２ラムダがL１ラムダの２倍

L２の方が強く効く。

![elastic-l2x2-50.gif](https://qiita-image-store.s3.amazonaws.com/0/247972/1435b81c-0ca6-bff6-a0a8-48d14206ccd0.gif)

---

**Group Lasso**

グループごとに重みがゼロになる。

- グループの分け方　

1. [bias]
2. [Pregnancies, Glucose] 
3. [Blood~, Skin~] 
4. [Insulin, BMI] 
5. [Diabetes~, Age]

![group-22-240.gif](https://qiita-image-store.s3.amazonaws.com/0/247972/b7a30ea0-cf9a-e47d-d91f-a8af7f0297d6.gif)

- グループの分け方　

1. [bias]
2. [Pregnancies, Glucose, Blood~, Skin~] 
3. [Insulin, BMI, Diabetes~, Age]

![group-44-240.gif](https://qiita-image-store.s3.amazonaws.com/0/247972/6375c342-79b0-2080-3891-090669039c1e.gif)



# 参考

- [近接勾配法とproximal operator](http://yamagensakam.hatenablog.com/entry/2018/02/14/075106)