Table of Contents

モチベーション

関数の種類

目的関数の最小値を求めたい機会は多々ある，目的関数を微分することで最小値を求めることがあるが，関数の形によっては，部分的に微分できない場合がある．例えば， $y=\|x\|$ は $x=0$ で微分不可能である．

劣微分・劣勾配

凸関数を最小化するために勾配情報が欲しいが，凸関数はいつでも微分可能であるとは限らない．例えば， $y=|x|$ は凸関数だが， $x=0$ で微分可能でない．このような関数に対して勾配（劣勾配）を定義してみる．

真凸関数 $f:\mathbb{R}^n \to \mathbb{R} \cup \{\infty\}$ かつ $x \in dom(f)$ の元で， $\forall y \in \mathbb{R}^n$ が

$$ f(\boldsymbol{y}) \geq f(\boldsymbol{x}) + \boldsymbol{g}^{\mathrm{T}}(\boldsymbol{y}-\boldsymbol{x}) $$

を満たすとき， $\boldsymbol{g}$ を $\boldsymbol{x}$ における劣勾配という．劣勾配の集合 $\partial f(\boldsymbol{x})$ を劣微分という．

$f(x)=|x|$ の $x=0$ における劣微分は，以下より， $\partial f(0) = [-1,1]$ と求まる．

もっと複雑な関数を微分したい

座標降下法(Coordinate Descent)または勾配降下法(Gradient Descent)を用いる．

座標降下法(Coodinate Descent)によるL1正則化(Lasso)の最適化

LassoとはLeast absolute shrinkage and selection operatorの略であり，以下の式で表される．以下を最小にするような $\boldsymbol{w}$ を探したいが，微分できない項を含んでいる．ここで， $\boldsymbol{x_n}^{\mathrm{T}} = (1, x_{n1} , x_{n2}, …, x_{nD})$ ， $\boldsymbol{w}^{\mathrm{T}} = (w_0, w_1, …, w_D)$ である． $\boldsymbol{w}_{-j}$ は列ベクトル $\boldsymbol{w}$ からj列目を除いたものを表す．

$$ E(\boldsymbol{w}) = \frac{1}{2} \sum_{n=1}^{N} \{\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x_n} - t_n\}^2 + \lambda \| \boldsymbol{w}_{-0} \| = \frac{1}{2} (\boldsymbol{Xw-t}) ^{\mathrm{T}} (\boldsymbol{Xw-t}) + \lambda \| \boldsymbol{w}_{-0} \| $$

$w_0$ は正則化項に含まれないので， $E(\boldsymbol{w})$ を $w_0$ で微分することで $w_0$ の最小値が求まる．

$$ \frac{ \partial E(\boldsymbol{w}) }{ \partial w_0 } = \sum_{n=1}^{N} \{\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x_n} - t_n\} = \sum_{n=1}^{N} \{w_0 + \boldsymbol{w_{-0}}^{\mathrm{T}} \boldsymbol{{x_n}_{-0}} - t_n\} = 0 \\ \therefore w_0 = \frac{1}{N} \sum_{n=1}^{N} \{t_n - \boldsymbol{w_{-0}}^{\mathrm{T}} \boldsymbol{{x_n}_{-0}}\} $$

次に， $E(\boldsymbol{w})$ を $w_k ( \neq 0, 0 \lt k \leq D )$ で偏微分した値を0とする．

$$ \frac{ \partial E(\boldsymbol{w}) }{ \partial w_k } = \boldsymbol{X_{(:,k)}^{\mathrm{T}}} \{ \boldsymbol{X_{(:,k)}} w_k + \boldsymbol{X_{(:,-k)}} \boldsymbol{w_{-k}} - \boldsymbol{t} \} + \lambda sign(w_k) = 0 $$

$\boldsymbol{X_{(:,k)}}$ は行列Xのk列ベクトルを表し， $\boldsymbol{X_{(:,k)}^{\mathrm{T}}}$ はそれの転置ベクトルを表す．また， $\boldsymbol{X_{(:,-k)}}$ は行列Xからk列ベクトルを除いたものを表す．

$w_k$ について整理すると， $w_k > 0$ の場合と $w_k <0$ の場合ができる．

$$ \begin{eqnarray} w_k = \begin{cases} \frac{1}{\boldsymbol{X_{(:,k)}^{\mathrm{T}}} \boldsymbol{X_{(:,k)}}} \{ \boldsymbol{X_{(:,k)}^{\mathrm{T}}} \{ \boldsymbol{t} - \boldsymbol{X_{(:,-k)}} \boldsymbol{w_{-k}} \} - \lambda\} & ( w_k \gt 0) & (1)\\ \frac{1}{\boldsymbol{X_{(:,k)}^{\mathrm{T}}} \boldsymbol{X_{(:,k)}}} \{ \boldsymbol{X_{(:,k)}^{\mathrm{T}}} \{ \boldsymbol{t} - \boldsymbol{X_{(:,-k)}} \boldsymbol{w_{-k}} \} + \lambda\} & ( w_k \lt 0) & (2) \end{cases} \end{eqnarray} $$

ここで劣微分を導入する．具体的には，(1)式または(2)式を満たさない場合は問答無用で $w_k=0$ とするだけ．λの大きさによって $w_k=0$ とする範囲を制御できる．これらをまとめてsoft thresholding operatorと呼ぶ．

で，以下の手順に沿って $\boldsymbol{w}$ を求める．以下のように，目的変数（最小化したい変数）を一つずつ触って，目的関数（今回はLasso）を最小化する方法を，座標降下法(Coodinate Descent)という．

  1. $\boldsymbol{w}$ を適当な値で初期化
  2. $w_0$ を計算
  3. while $\boldsymbol{w}$ が収束条件を満たすまで
    • for k in 1,2,…,D
      • $w_k$ を計算
    • $w_0$ を計算
import numpy as np
import copy

class MyLasso:
  def __init__(self, alpha=1.0, max_iter=100):
    self.alpha = alpha
    self.max_iter = max_iter
    self.coef_ = None
    self.intercept_ = None

  def fit(self, X, t):
    X = np.column_stack((np.ones(len(X)),X)) # (1,x1,x2,...)
    w = np.zeros(X.shape[1])
    
    w[0] = np.sum(t - np.dot(X[:,1:], w[1:]))/(X.shape[0])
    
    for iter in range(self.max_iter):
        for i in range(1, len(w)):
            _w = copy.deepcopy(w)
            _w[i] = 0.0

            w[i] = 0.0
            a = (np.dot(X[:, i], t - np.dot(X, _w)) - self.alpha)/(X[:, i]**2).sum()
            if a > 0.0 : w[i] = a
            a = (np.dot(X[:, i], t - np.dot(X, _w)) + self.alpha)/(X[:, i]**2).sum()
            if a < 0.0 : w[i] = a

        w[0] = np.sum(t - np.dot(X[:,1:], w[1:]))/(X.shape[0])

        self.intercept_ = w[0]
        self.coef_ = w[1:]

    return self

  def predict(self, X):
    y = np.dot(X, self.coef_)
    y += self.intercept_*np.ones(len(y))
    return y

$y=x^2+2x+3$ にσ=1.0のノイズを加えた50点を教師データとして，Lassoで予測した結果が以下である．

from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures

def f(x):
    return  x**2 + 2*x + 3

N = 50
X = np.zeros((N,1))
np.random.seed(seed=2)
X[:,0] = np.random.uniform(-3, 3, N)
noise = 1.0*np.random.randn(N)
t = f(X[:,0])+ noise

X_train, X_test, t_train, t_test = train_test_split(X, t, random_state=0)

fig, ax = plt.subplots(1, 3, figsize = (15, 5))

for degree, ax in zip([1,3,9],ax):
    pipeline = Pipeline([
        ('polynominal_features', PolynomialFeatures(degree=degree)),
        ('linear_regression', MyLasso(alpha=0.01))
    ])

    lr = pipeline.fit(X_train, t_train)
    print(pipeline.steps[1][1].intercept_)
    print(pipeline.steps[1][1].coef_)
    
    Xcont = np.linspace(np.min(X), np.max(X), 200).reshape((200, 1))
    ax.plot(Xcont, f(Xcont), label='original')
    ax.plot(X, t,'.', label='training data')
    ax.plot(Xcont, lr.predict(Xcont), label='predictive')
    ax.set_title('degree : {}'.format(degree))
    ax.legend()
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$t$')

多項式モデルで予測し，具体的には１次，３次，５次のモデルで予測してみた．

linear_regression_with_lasso

各モデルが予測した式は以下（係数を小数点以下第一位に丸めている）であり，３次以上の項の係数は0か非常に小さな値になっている．これはL1正則化によりスパースなモデルを求めることができたことを意味する．

１次式モデル： $1.8x+5.5$

３次式モデル： $0.0x^3 + 0.9x^2 + 2.1x + 3.1$

５次式モデル： $0.0x^5 + 0.0x^4 - 0.1x^3 + 0.9x^2 + 2.2x + 3.1$

Lassoを考慮したモデルがなぜスパースになるのかを知るには，ラグランジュの未定乗数法を知ると良い

ラグランジュの未定乗数法とは

制約条件 $g(x,y)=0$ が与えられた状態で，関数 $f(x,y)$ の最小値または最大値を求めるために便利なのがラグランジュの未定乗数法である．具体的には，ラグランジュ関数 $L(x,y)$ とラグランジュ乗数 $\lambda$ を用いて，以下を解けば良い．

$$ L(x,y,\lambda) = f(x,y) - \lambda g(x,y) \\ \frac{ \partial L }{ \partial x } = \frac{ \partial L }{ \partial y } = \frac{ \partial L }{ \partial \lambda } = 0 $$

なぜこれで停留点を求めることができるか

$$ \frac{ \partial L }{ \partial x } = \frac{ \partial L }{ \partial y } = 0 $$

を行列で書き直すと，

$$ \begin{pmatrix} \frac{ \partial f(x,y) }{ \partial x } \\ \frac{ \partial f(x,y) }{ \partial y } \end{pmatrix} = \lambda \begin{pmatrix} \frac{ \partial g(x,y) }{ \partial x } \\ \frac{ \partial g(x,y) }{ \partial y } \end{pmatrix} $$

となり，ベクトル $\partial g$ をλ倍したものがベクトル $\partial f$ であることを表している．では，ベクトル $\partial g$ とは何か．

制約面 $g(x,y)=0$ 上の点 $(x_0,y_0)$ と $(x_0+\delta x, y_0+\delta y)$ を考えると， $(x_0,y_0)$ の周りでのテイラー展開より以下が得られる．

$$ g(x_0+\delta x, y_0+\delta y) \simeq g(x_0,y_0) + \left. \frac{ \partial g(x,y) }{ \partial x } \right|_{ (x_0,y_0) } \delta x + \left. \frac{ \partial g(x,y) }{ \partial y }\right|_{ (x_0,y_0) } \delta y $$

いま， $(x_0,y_0)$ と $(x_0+\delta x, y_0+\delta y)$ は制約面上に存在しているので， $g(x_0,y_0)=g(x_0+\delta x, y_0+\delta y)$ が成立するためには，以下が成り立つ必要がある．

$$ \left. \frac{ \partial g(x,y) }{ \partial x } \right|_{ (x_0,y_0) } \delta x + \left. \frac{ \partial g(x,y) }{ \partial y } \right|_{ (x_0,y_0) } \delta y = 0 $$

内積で書き直すと，

$$ \begin{pmatrix} \frac{ \partial g(x,y) }{ \partial x } \\ \frac{ \partial g(x,y) }{ \partial y } \end{pmatrix}_{ (x_0,y_0) } \cdot \begin{pmatrix} \delta x \\ \delta y \end{pmatrix} = 0 $$

$\delta x$ と $\delta y$ が十分に小さければ， $(\delta x, \delta y)$ は点 $(x_0,y_0)$ における制約面 $g$ の接線を表すので，ベクトル $(\frac{ \partial g(x,y) }{ \partial x }, \frac{ \partial g(x,y) }{ \partial y })$ は接線の法線ベクトルを表す．

先程の

$$ \begin{pmatrix} \frac{ \partial f(x,y) }{ \partial x } \\ \frac{ \partial f(x,y) }{ \partial y } \end{pmatrix} = \lambda \begin{pmatrix} \frac{ \partial g(x,y) }{ \partial x } \\ \frac{ \partial g(x,y) }{ \partial y } \end{pmatrix} $$

より，ベクトル $\partial g$ とベクトル $\partial f$ が点 $(x_0,y_0)$ で並行であれば，その点で $f(x,y)$ が停留点を持つ．

「ベクトル $\partial g$ とベクトル $\partial f$ が点 $(x_0,y_0)$ で並行である」という部分が肝である．例えば，下図において，制約条件 $g(x,y)=0$ は赤の等高線， $f(x,y)$ はc1,c2,c3の等高線で表し，c1 < c2 < c3 とすると，勾配 $\nabla f$ は図のようになる(わざと $g(x,y)=0$ と $f(x,y)=c2$ が接するように描いている)．ラグランジュ乗数法により勾配 $\nabla g$ と勾配 $\nabla f$ が並行な点，すなわち(x0,y0)において， $g(x,y)=0$ を満たす $f(x,y)$ の停留点（下図の場合は最小値)となることがわかる．

lagrange_multiplier

$\nabla g$ に対する1次独立性が必要

制約条件 $g$ が複数ある場合，ラグランジュの未定乗数法によって最適解を求めるためには， $\nabla g$ がそれぞれ1次独立である必要がある．1次独立な $\nabla g$ と $\lambda$ の線形和によって， $\nabla f$ を表現する必要があることが理由．

最適性の２次の十分条件を満たせば局所最適解

関数 $f,g_1,…,g_p$ は二階微分可能とし，局所最適解を $x^*$ とする． $\nabla g_1(x^*),…,\nabla g_p(x^*)$ は１次独立であり，ラグランジュ乗数を $\lambda^*_1,…,\lambda^*_p \in \mathbb{ R }$ とする．

ここで，部分空間 $S$ を以下のように定義する．

$$ S = \{y \in \mathbb{ R } | \nabla g_i(x^*)^T y =0, i=1,...,p \} $$

更に，行列 $L \in \mathbb{ R }^{n \times n}$ を以下のように定義する．

$$ L = \nabla^2 f(x^*) + \displaystyle \sum_{ i = 1 }^{ m } \lambda_i^* \nabla^2 g_i(x^*) $$

これら $S,L$ について， $y \in S, y \neq 0$ を用いて，以下を満たせば $x^*$ は局所最適解である．

$$ y^T L y \gt 0 $$

つまり，ヘッセ行列$L$が正定値であれば， $x^*$ は局所最適解である．

ちなみに，２次の必要条件は，上記ヘッセ行列が半正定値であるという，十分条件を緩くしたものになる．

ラグランジュ乗数法のKKT条件

これまでは制約条件 $g(x,y)$ が0である元で議論してきたが， $g(x,y)$ に範囲がある場合を考える．

下図において， $g(x,y)=0$ を赤の等高線で表し，内側に向かって $g(x,y)$ が大きくなっているとすると， $g(x,y) \gt 0$ は楕円の内部に該当し，勾配 $\nabla g$ は図のような向きになる．また， $f(x,y)$ を緑の等高線で表し，c1 < c2 < c3 であるとすると，勾配 $\nabla f$ は図のような向きになる．

lagrange_multiplier_kkt1

この場合，「制約条件 $g(x,y) \geq 0$ の元で $f(x,y)$ を最大にする点」をラグランジュ乗数法により，点(x0,y0)と求めることができる．

また，「制約条件 $g(x,y) \geq 0$ の元で $f(x,y)$ を最小にする点」は，ラグランジュ乗数法により，下図の点(x1,y1)と求めることができる．

lagrange_multiplier_kkt2

更に，下図の場合，「制約条件 $g(x,y) \geq 0$ の元で $f(x,y)$ を最大にする点」は，c1 < c2 < c3 であるとするなら，ラグランジュ乗数法のλを0とし，単に $f(x,y)$ の最大点を求めれば良いので，点(x2,y2)と求まる．

lagrange_multiplier_kkt3

これらの条件をまとめたものが以下のKKT条件(Karush-Kuhn-Tucker condition)である．

$$ L(x,y,\lambda) = f(x,y) + \lambda g(x,y) \\ \frac{ \partial L }{ \partial x } = \frac{ \partial L }{ \partial y } = \frac{ \partial L }{ \partial \lambda } = 0 \\ g(x,y) \geq 0 \\ \lambda \geq 0 \\ \lambda g(x,y) = 0 $$
$$ L(x,y,\lambda) = f(x,y) - \lambda g(x,y) \\ \frac{ \partial L }{ \partial x } = \frac{ \partial L }{ \partial y } = \frac{ \partial L }{ \partial \lambda } = 0 \\ g(x,y) \geq 0 \\ \lambda \geq 0 \\ \lambda g(x,y) = 0 $$

条件式がややこしく見えるが，言っていることは単純で，「最大化したいのか最小化したいのかに注意し，勾配 $\nabla f$ と $\nabla g$ の向きをλの正負で制御して，ラグランジュ乗数法を解け」ということである．

$\lambda g(x) = 0$ は相補性条件と呼ばれ， $\lambda$ または $g(x)$ のどちらかが0になることを表している．

正則化とラグランジュ乗数法

以下の二乗和誤差と正則化項を合わせた関数 $E(\boldsymbol{w})$ を最小化したいとする．ξは正則化項の係数であり，正であるとする．

$$ E(\boldsymbol{w}) = \frac{1}{2} \sum_{n=1}^{N} \{\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x_n} - t_n\}^2 + \xi \sum_{i=1}^{D} | w_i |^q $$

いま，以下の二乗和誤差の最小値を

$$ f(\boldsymbol w) = \frac{1}{2} \sum_{n=1}^{N} \{\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x_n} - t_n\}^2 $$

以下の制約条件で求めるとする．

$$ g(\boldsymbol w) = \eta - \sum_{i=1}^{D} | w_i |^q $$

$f(\boldsymbol w)$ の最小化には，KKT条件 $g(\boldsymbol w) \geq 0$ と $\lambda \geq 0$ の元で，以下のラグランジュ関数を解く必要がある．

$$ L(\boldsymbol {w},\lambda) = f(\boldsymbol w) - \lambda g(\boldsymbol w) \\ \frac{ \partial L }{ \partial \boldsymbol{w} } = 0 $$

ここで，

$$ \frac{ \partial L }{ \partial \boldsymbol{w} } = \frac{ \partial}{ \partial \boldsymbol{w} } \{ \frac{1}{2} \sum_{n=1}^{N} \{\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x_n} - t_n\}^2 + \lambda \sum_{i=1}^{D} | w_i |^q \} = 0 $$

であるので， $\lambda = \xi$ とすると，内部的には，二乗和誤差+正則化の最小値を求めていることになる． $\xi > 0$ であるので， $\lambda > 0$ である必要が出てくる．KKT条件 $\lambda g(\boldsymbol w) = 0$ より， $g(\boldsymbol w) = 0$ なので以下が成り立つ．

$$ \eta = \sum_{i=1}^{D} | w_i |^q $$

q=1のときの $f(\boldsymbol w)$ と $g(\boldsymbol w)$ を以下のように図示してみた．

visualize_of_lasso

更に，断面を簡易的に図示してみた．ラグランジュ乗数法の仕組みにより，制約条件$g$の元で，関数$f$の停留点を求めるには，勾配ベクトル $\nabla g$ のスカラーλ倍したものが，勾配ベクトル $\nabla f$ と並行である必要があるので，停留点はもちろん，ηもλに依存する．

fragment_of_lasso

上図の停留点においてw1=0となっている．疎（スパース）な解が得られている例であり，lassoの特徴である．この特徴を用いることで，モデルの複雑さを抑えることができるが，必ずスパースな解が得られるとは限らない．例えば，下図のような場合はスパースでない解が得られ，関数f(w)が制約条件g(w)>0内に入っている場合は，関数f(w)そのものの最小値となる．

visualize_of_lasso_without_sparse

近接勾配法

たぶんつづく