ガウスカーネルでl1正則化付き回帰をやってみる

$l2$正則化であれば正則化項が微分可能なので特に難しいことはなく、ラグランジュの未定乗数法で解くことができました。

今回は基本的なモデルはガウスカーネル回帰で固定して、正則化項の部分だけ$l1$正則化に変えてみましょう。

$l2$正則化に比べて$l1$正則化は得られる解がスパースになりやすいという特徴があります。

解がスパースというのは、成分が0になっている要素が多い、という意味です。
解がスパースであれば、重みが非零の部分の特徴だけ見ればいいので効率的ですし、そもそも特徴のなかであまり予測に必要のないものをあぶり出すなどの分析にも役立ったりします。

しかし、$l1$正則化を用いた場合は、原点で正則化項が微分不可能なのでラグランジュの未定乗数法がそのまま適用することができません。

今回は拡張ラグランジュ関数を使って、交互方向乗数法によって$l1$正則化をつけたガウスカーネルの最適化を行ってみます。

### 問題設定

$$
\min_{\theta} \frac{1}{2}|K\theta - y|^2 + \lambda |\theta|
$$

をときます。
ただし$\bm K$はグラム行列とし、$\lambda$は正則化係数とします。

### 拡張ラグランジュ関数と交互方向乗数法

$$
\min_{\bm{\theta,z}} f(\bm \theta) + g(\bm z)
$$

$$
\mathrm{s.t.}\ \bm{A\theta} + \bm{Bz} -\bm c = \bm 0 
$$

を解くためのアルゴリズムです。
この時、

拡張ラグランジュ関数を

$$L(\bm{\theta,z,u})= f(\bm \theta)+g(\bm z)+ \bm u^T(\bm{A\theta} + \bm{Bz} - \bm c) + \frac{1}{2} |\bm{A\theta} + \bm{Bz}-\bm c|^2$$

と定義します。($\bm u$は未定乗数)

このとき、以下のように更新式を作ります。

$$
  \bm\theta^{(k+1)}  = \mathrm {argmin}_{\bm \theta}\  L(\bm\theta,\bm z^{(k)},\bm u^{(k)})$$

$$\begin{aligned}
  \bm z^{(k+1)} &= \mathrm{argmin}_{\bm z}\  L(\bm\theta^{(k+1)},\bm z,\bm u^{(k)}),\cr
  \bm u^{(k+1)} &= \bm u^{(k)} + \bm\theta^{(k+1)} - \bm z^{(k+1)}
\end{aligned}$$

今回は、これを用いて$l1$正則化項つきガウスカーネルモデルを解きます。

### 更新式の導出

ガウスカーネルモデルのスパース回帰に対して交互方向乗数法を適用します。

$$\begin{aligned}
    f(\theta)&= \frac{1}{2}|\bm{K\theta} - \bm{y}|^2,\cr
    g(z) &= \lambda|\bm z|,\cr
    \bm A &= -\bm B = \bm I ,\cr
    c &= \bm{0}
\end{aligned}$$

交互方向乗数法を適用すると、各パラメータの更新式は以下のように導出できます。

拡張ラグランジュ関数は

$$L(\bm{\theta,z,u})= f(\bm \theta)+g(\bm z)+ \bm u^T(\bm{A\theta} + \bm{Bz} - \bm c) + \frac{1}{2} |\bm{A\theta} + \bm{Bz}-\bm c|^2$$

です。
$\theta$に関する更新式は、$\frac{\partial L}{\partial \bm \theta} =\bm 0$となる$\bm\theta$を求めればよく、これは普通に微分できるので

$$
\bm \theta^{(k+1)}= (\bm K^2 + \bm I)^{-1}(\bm{Ky}+\bm z^{(k)} - \bm u^{(k)})
$$

となります。

$\bm z$に関する更新式は、$L(\bm \theta^{(k+1)}, \bm z, \bm u)$を最小にする$z$を求めればよいですが、普通に$\bm z$について微分はできません。

$$
L(\bm \theta^{(k+1)}, \bm z , \bm u^{(k)}) = \lambda |\bm z| +  \bm u^{(k)}\cdot(\bm \theta^{(k+1)} -\bm z )+ \frac{1}{2}|\bm \theta^{(k+1)}- \bm z|^2
$$

これは成分ごとの関数の和としてかけるので、成分ごとに独立にソフト閾値処理を行うことで最小にする$z$を決めることができます。

ソフト閾値処理とは

$$
  T(z)=\lambda |z| + u(\theta -z) + \frac{1}{2}(\theta -z)^2
$$
のとき、
$$
  \mathrm{argmin} \ T(z) = \max( 0,\theta+u-\lambda) + \min(0, \theta+u+\lambda)
$$

となることを指します。これは単純に絶対値が入った二次関数の最小なので、高校数学の知識で確かめられるので割愛です。

なので、これをつかうと、

$$
\bm z^{(k+1)}=\max(\bm 0,\bm \theta^{(k+1)}+\bm u^{(k)}-\lambda \bm 1) + \min(\bm 0,\bm \theta^{(k+1)}+\bm u^{(k)}+\lambda \bm 1)
$$
と更新式をかくことができます。
ただし$\bm 1$はすべての成分が1のベクトル、$\max$や$\min$は成分ごとの演算とします。

残りの$\bm u$に関する更新式は

$$
\bm u^{(k+1)} = \bm u^{(k)} + \bm \theta^{(k+1)} - \bm z^{(k+1)}
$$
となります。

### 実装
pythonで実装してみました。

基本的には前回と同じで、fit関数の部分だけ違います。

```python

class GaussianKernelNorm1(BaseEstimator, RegressorMixin):
    def __init__(self, h, c):
        """initialize
        :param h: ガウス幅
        :param c: 正則化係数
        """
        self.h = h
        self.c = c

def k(self, x_i, x_j):
        """kernel function"""
        return math.e**(-np.dot(x_i - x_j, x_i - x_j) / (2 * self.h**2))

def create_K(self, X):
        """create Kernel matrix

:param X: data matrix (data_n,feature_n)
        :return K: kernel matrix (data_n,data_n)
        """
        K = []
        for x_i in X:
            K_i = []
            for x_j in X:
                K_i.append(self.k(x_i, x_j))
            K.append(K_i)
        return np.array(K)

def fit(self, X, y, max_iter=10000, epsilon=0.0001):
        self.x_ = X
        K = self.create_K(X)
        n = K.shape[0]
        theta_next, z_next, u_next = self.init_params(n)
        for i in range(max_iter):
            u_prev = u_next
            z_prev = z_next
            theta_prev = theta_next
            theta_next, z_next, u_next = self.update_params(
                z_prev, u_prev, y, K, n)
            if np.linalg.norm(theta_next - theta_prev, ord=2) < epsilon:
                print(
                    f"iter {i+1} loss {self.loss(theta_next,y,K)} (Converged!)")
                break
            if (i + 1) % 100 == 0:
                print(f"iter {i+1} loss {self.loss(theta_next,y,K)}")
        self.coef_ = theta_next
        return self

def loss(self, theta, y, K):
        return np.linalg.norm(np.dot(K, theta) - y, ord=2) / 2 + self.c * np.linalg.norm(theta, ord=1)

def init_params(self, n):
        return np.random.rand(), np.random.rand(), np.random.rand()

def update_params(self, z_prev, u_prev, y, K, n):
        theta_next = np.linalg.solve(
            np.dot(K, K) + np.eye(n), np.dot(K, y) + z_prev - u_prev)
        z_next = np.maximum(np.zeros(n), theta_next, u_prev - self.c * np.ones(
            n)) + np.minimum(np.zeros(n), theta_next, u_prev + self.c * np.ones(n))
        u_next = u_prev + theta_next - z_next
        return theta_next, z_next, u_next

def predict(self, X):
        check_is_fitted(self, 'coef_')
        check_is_fitted(self, 'x_')
        pred_y = []
        for x_i in X:
            k_i = np.array([self.k(x_i, x_j) for x_j in self.x_])
            pred_y.append(np.dot(self.coef_, k_i))
        return np.array(pred_y)

def evaluate(self, X, y):
        pred_y = self.predict(X)
        return mean_squared_error(pred_y, y)

```
 ### おわりに
 
 前回の記事では$l2$正則化のモデルを作ったので、比較の実験などをしてみても面白いかもしれません。
ありがとうございました。

このブログを検索

ぐぐりらにっき

ガウスカーネルでl1正則化付き回帰をやってみる

コメント

コメントを投稿

ガウスカーネルでl1正則化付き回帰をやってみる

この記事をシェアする

コメント

コメントを投稿