【改】PRML第4章のロジスティック回帰をPythonで実装

この記事では、PRML第4章で述べられている、(多クラス)ロジスティック回帰を実装します。
(NOTE: この記事は、「PRML第4章のロジスティック回帰をPythonで実装」の一部を書き直したものです。具体的には、最小化アルゴリズムの部分を書き直しました：古い方の記事ではコスト函数の最小化にscipy.optimize.minimizeを用いていましたが、この記事では自分で書いたソルバーを用いています。)

対応するjupyter notebookは筆者のgithubリポジトリにあります。連載の方針や、PRMLの他のアルゴリズムについては、まとめページ( https://qiita.com/amber_kshz/items/e912a59ffc1365cb6adf )をご覧ください。

前回と同様、PRMLの本に載っている理論(+$\alpha$)のおさらいをし、数式をコードに翻訳して、最後に実際に動かした結果を見ることにします。また、最小化を行う箇所では、複数の最小化アルゴリズムを試し、それらの挙動の違いを見ることにします。

1 理論のおさらい

まず、PRMLの4章(正確には4.3.4節)に書かれていることのおさらいです。すこーしだけ本と記号の定義が異なります¹が、それ以外は同じなので、内容を知っている方は軽く読み飛ばしていただいて大丈夫です。

1.1 設定

以下のように記号を定義します。

• $N \in \mathbb{N}$ : を訓練データの個数
• $d \in \mathbb{N}$ : 入力(特徴量)の次元
• $\mathcal{C} = \left\{ 0, 1, \dots, C-1 \right\}$ : ありうるラベルの集合
• $x_0, x_1, \dots , x_{N-1} \in \mathbb{R}^d$を訓練データの入力(特徴量)、$y_0, y_1, \dots, y_{N-1} \in \mathcal{C}$を訓練データのラベル、$y := {}^t (y_0, \dots, y_{N-1}) \in \mathcal{C}^N$とします。
• $t_n \in \left\{ 0,1\right\}^{C}$ を$y_n$の1-of-$C$符号化、($y_n = c$なる$n$に対して、$t_n$は第$c$成分だけが1で他は0のベクトル)とします。

1.2 モデル

多クラスロジスティック回帰では、次のモデルを仮定します(PRMLの(4.104)式参照。)。
$$
\begin{align}
p(y|x,\theta) = \frac{\exp\left[{\theta^{(y)}}^T \phi(x)\right] }{ \sum_{y'=0}^{C-1} \exp\left[{\theta^{(y')}}^T \phi(x)\right] }
\end{align}
$$
ただし、

• $\phi : \mathbb{R}^d \rightarrow \mathbb{R}^M $, $\phi(x) = (\phi_0(x), \dots, \phi_{M-1}(x))^T$は基底関数、$M$は基底関数の本数、
• $\theta^{(y)} \in \mathbb{R}^M$は重みを表すパラメーター(学習時にはこのパラメーターを調整します。)、
• $\theta = \left( {\theta^{(0)}}^T, \dots, {\theta^{(C-1)}}^T \right)^T \in \mathbb{R}^{CM}$　です。

後の実装の際には、以下の"線型"基底関数を用います：
$$
\begin{align}
\phi_j(x) = \begin{cases}
1 & (j=0)\
x^{(j-1)} & (j = 1.\dots, d)
\end{cases}
\end{align}
$$
ただし、$x^{(j)}$は$x \in \mathbb{R}^d$の第$j$成分を表します。

また、$\Phi$を以下で定められる$N \times M$行列とします:
$$
\begin{align}
& \Phi = (\Phi_{i,j})_{i=0,\dots, N-1, j=0,\dots, M-1}, \
& \Phi_{i,j} = \phi_j(x_i)
\end{align}
$$

1.3 損失函数とその微分

1.3.1 損失函数

ロジスティック回帰の学習では、次の損失函数を最小化するような$\theta$を選びます(PRMLの(4.107), (4.108)参照。)。第1項は尤度の対数に$-1$をかけたもの、第2項は$l^2$正則化項です。
$$
\begin{align}
J(\theta) := - \frac{1}{N} \sum_{n=0}^{N-1} \sum_{c=0}^{C-1} t_{n,c} \log p_{n,c} + \frac{\lambda}{2N} \| \theta \|^2,
\end{align}
$$
ただし、
$$
\begin{align}
\lambda & > 0 \
p_{n,c} &:= p(c|x_n,\theta) \
t_{n,c} &:= \delta_{y_n, c}
\end{align}
$$
です。

1.3.2 勾配

後で損失函数の最小化を行うために、損失函数の$\theta$に関する微分(勾配)を用います。
上記の損失函数を$\theta^{(c)}_{j}$で微分することにより、次の勾配が得られます(PRMLの(4.109)参照。)。
$$
\begin{align}
\frac{\partial J}{\partial \theta^{(c)}_{j}} = \frac{1}{N} \sum_{n=0}^{N-1} ( p_{n,c} - t_{n,c}) \phi_j(x_n) +\frac{\lambda}{N} \theta^{(c)}_{j}
\end{align}
$$

1.3.3 Hessian

損失函数の最小化は勾配があれば実行可能ですが、2階微分(Hessian)の情報を追加で用いることもできます。それが、PRMLの.4.3.3節で説明されているNewton法による最小化です²。

多クラスロジスティック回帰のためのHessianは、以下のようになります(PRML(4.110)式)
$$
\begin{align}
H_{(c, j), (c', j')} &:= \frac{\partial^2 J}{\partial \theta^{(c)}_{j} \partial \theta^{(c')}_{j'}} \
&= \frac{1}{N} \sum_{n=0}^{N-1} p_{n, c} \left( \delta_{c, c'} - p_{n, c'}\right) \Phi_{n, j} \Phi_{n, j'}
+ \frac{\lambda}{N} \delta_{c, c'} \delta_{j, j'}
\end{align}
$$

1.4 損失函数の最小化

さて、上に与えた損失函数の形を見れば分かるとおり、ロジスティック回帰の損失函数の最小化は解析的に手に負えません³。
ということで、手に負えないときは数値計算! 数値的になんとか損失函数を最小化したいと思います。ここでは2つの手法を用いることにします。

この節では一般的な記法を用いるために、

• $f$を最小化したい実数値函数
• $x$をその変数
• $H$を$f$のHessian

とします。

また、ここでは損失函数は下に凸⁴なので、

• 最小値が存在すれば一意である
• 数値的な最適化で$x$を更新する手法の場合、初期値にこだわる必要がない

ことに注意します。

1.4.1 最急降下法

最急降下法⁵では、変数$x$を、函数$f$が一番大きく減少する方向、つまり勾配の逆方向に更新します。
数式で表すと、$x$を次のように更新していきます($x_n$を$n$ステップ目での$x$の値とします。)。
$$
\begin{align}
x_{n+1} = x_n - \alpha \nabla f(x_n)
\end{align}
$$
ただし、$\alpha > 0$は学習率と呼ばれるパラメーターであり、人が手で決めることになります⁶。

1.4.2 Newton法

Newton法(または多次元の場合はNewton-Raphson法)では、勾配に加えて、2階微分の情報を活用します。Newton法では、
$$
\begin{align}
x_{n+1} = x_{n} - \left[H(x_n)\right]^{-1} \nabla f(x_n)
\end{align}
$$
のように変数を更新していきます。

ここで、第2項は$\left[H(x_n)\right]^{-1} \nabla f(x_n)$いくつかの手法で計算することができます。

• 1つは、数値的に「厳密に」求める手法です。具体的には、$H(x)$の逆行列をとって$\nabla f(x)$にかける、または線型方程式$H(x)y = \nabla f(x)$をLU分解などで$y$について解くなどの手法です⁷。
• もう1つは、線型方程式$H(x)y = \nabla f(x)$を共役勾配法などの繰り返し手法で解くことです⁸。

前者の方が素直に見えますが、後者にも利点があります: 共役勾配法を用いる場合、係数行列そのものは必要でなく、「任意のベクトルが与えられたとき、係数行列とそのベクトルの積を与える函数」さえあれば十分です。この性質のありがたみについては、後ほど詳しく見たいと思います。

1.4.3 最急降下法とNewton法、どちらを用いるか?

最急降下法では勾配の情報があれば十分だったのに対し、Newton法ではHessianの情報も必要になり、$x$の更新にかかる時間も大きくなります。一方で、Newton法は更新の回数という意味では勾配降下法より速く収束することが知られています。

この記事では、両方の手法(特にNewton法については上にも指摘した2つの手法)を実装し、その振る舞いを比較することにします。

また、Andrew Ng先生のこちらの動画も参考になるかと思います　https://www.youtube.com/watch?v=iwO0JPt59YQ

2 数式からコードへ

ではいよいよ、ここまでの数式をコードに翻訳していきましょう！次のようなステップで進んで行こうと思います。

• まず、ここまで出てきた変数を行列(2-D array)で表すための函数を書きます。
• 損失函数やその微分を表す函数を定義します。
• 与えられた函数を最小化する函数も定義しておきます(学習の際に用います。)。
• 以上の函数を組み合わせて、Logistic回帰分類器を書きます。

2.1 行列表記

まず、ここまで出てきた$p$や$t$などの2つの添え字を伴う変数を行列(正確には2-D numpy array)で表しておきます。こうしておくと、後の計算が行列の積などで書け、いろいろと便利になります。
具体的には、arrayを以下のように定義します。

• Phi : (N, M) array, Phi[n, m] := $\phi_m(x_n)$
• Theta : (C, M) array, Theta[c, j] := $\theta^{(c)}_{j}$
• T : (N,C) array,　T[n, c] := $t_{n, c}$
• P : (N, C) array, P[n, c] := $p_{n, c}$

それぞれのarrayを生成する函数を定義しておきましょう。

まず、Phiから。これは(今の設定では)1を含む列を加えるだけなので容易です。

def calcPhimat(X):
    '''
    This function generates the design matrix Phi from X, the input data
        
    Parameters
    ----------
    X : 2-D numpy array
        (N,d) numpy array, with X[n, i] = i-th element of x_n
            
    Returns
    ----------
    Phi : 2-D numpy array
        The design matrix
    '''
    N = len(X)
    if len(np.shape(X)) == 1:
        d = 1
    else:
        d = np.shape(X)[1]
    Phi = np.zeros((N,d+1))
    Phi[:,0] = np.ones(N)
    Phi[:,1:] = np.reshape(X, (N,d))
    return Phi

次に、Tを生成する函数。「bool型をfloat型のarrayに代入するとTrueが1に、Falseが0に変換される」ことを利用します。

def calcTmat(y, C):
    '''
    This function generates the matrix T from the training label y and the number of classes C
    
    Parameters
    ----------
    y : 1-D numpy array
        The elements of y should be integers in [0, C-1]
    C : int
        The number of classes
        
    Returns
    ----------
    T : (len(y), C) numpy array
        T[n, c] = 1 if y[n] == c else 0
    '''
    N = len(y)
    T = np.zeros((N, C))
    for c in range(C):
        T[:, c] = (y == c)
    return T

最後に、Pです。
Pを$P$、Phiを$\Phi$、Thetaを$\Theta$と行列で表すと、
$$
\begin{align}
P_{n,c} &= \frac{\exp\left[\sum_{j=0}^{M-1} \theta^{(c)}_{j} \phi_{j}(x_n)\right] }{ \sum_{c'=0}^{C-1} \exp\left[ \sum_{j=0}^{M-1} \theta^{(c')}_{j} \phi_{j}(x_n) \right]} \
&= \frac{ \exp \left[ \left( \Phi \Theta^T \right)_{n, c'} \right] }{ \sum_{c'=0}^{C-1} \exp\left[ \left( \Phi \Theta^T \right)_{n, c'} \right] }
\end{align}
$$
となります。このように行列の積を利用すると、対応するコードは以下のように書けます。

def calcPmat(Theta, Phi):
    '''
    This function generates the matrix P from the weight Theta and the design matrix Phi
    
    Parameters
    ----------
    Theta : 2-D numpy array
        Matrix representing the weight parameter
    Phi : 2-D numpy array
        The design matrix
        
    Returns
    ----------
    P : 2-D numpy array
    '''
    P = np.exp( Phi @ (Theta.T) )
    P = P/np.reshape( np.sum(P, axis= 1), (len(Phi),1)  ) 
    return P

2.2 損失函数とその微分

2.2.1 損失函数

上記の行列表記を用いると、
$$
\begin{align}
J(\Theta)
&= - \frac{1}{N} \sum_{n=0}^{N-1} \sum_{c=0}^{C-1} t_{n,c} \log p_{n,c} + \frac{\lambda}{2N} \| \theta \|^2 \
&= -\frac{1}{N} sum(T \ast \log P) + \frac{\lambda}{2N} \| \Theta \|_{2}^{2}
\end{align}
$$
となります。ただし、$\ast$は行列の要素ごとの積、$\log$は行列の要素ごとに対数をとることを表し、$sum$は行列の全要素の和を取ることを表します。

def cost_function(thtvec, Phi, T, lam):
    '''
    This function calculate the loss function and its gradient
    
    Parameters
    ----------
    thtvec : 1-D numpy array 
        (M*C,) array, which represents the weight parameter Theta, in flattened form
    Phi : 2-D numpy array
        (N, M) array, design matrix
    T : 2-D numpy array
        (N, C) array, where T[n, c] = 1 if y[n] == c else 0
    lam : float
        The regularization constant
    
    Returns
    ----------
    J : float
        The value of cost function        
    '''
    N, M = np.shape(Phi)
    C = np.shape(T)[1]
    Theta = np.reshape(thtvec, (C, M))
    P = calcPmat(Theta, Phi)
    J = -1.0/N*np.sum(T*np.log(P)) + lam/(2.0*N)*np.linalg.norm(thtvec)**2
    return J

2.2.2 勾配

また、$\frac{\partial J(\Theta)}{\partial \Theta}$を(C,M)行列とみなすと、
$$
\begin{align}
\frac{\partial J(\Theta)}{\partial \Theta}
&= \left( \frac{1}{N} \sum_{n=0}^{N-1} ( p_{n,c} - t_{n,c}) \phi_j(x_n) + \frac{\lambda}{N} \theta^{(c)}_{j} \right)_{c,j} \
&= \frac{1}{N} (P-T)^T \Phi + \frac{\lambda}{N} \Theta
\end{align}
$$
と書けます。

この記事での最小化のソルバー(と通常の多くのソルバー)は勾配がベクトル(1-D array)で与えられることを仮定しているので、この「勾配行列」を最後いnp.reshapeでベクトルに変換します。

def grad_cost_function(thtvec, Phi, T, lam):
    '''
    This function calculates the gradient of the cost function 
    
    Parameters
    ----------
    See the help for cost_function
        
    Returns
    ----------
    grad_vec : 2-D array
        The gradient of the cost function with respect to the weight parameter Theta, in flattened form
    '''
    N, M = np.shape(Phi)
    C = np.shape(T)[1]
    Theta = np.reshape(thtvec, (C, M))
    P = calcPmat(Theta, Phi)
    grad_mat = 1.0/N*((P - T).T) @ Phi + lam/N*Theta
    grad_vec = np.reshape(grad_mat, len(thtvec))
    return grad_vec

2.2.3 Hessian

Hessianの計算はもう少しややこしくなります。
念のため、ここでもう一度Hessianの表式を見ておきましょう(先ほどとは添え字を変更しています)
$$
\begin{align}
H_{(i, j), (k, l)} &= \frac{\partial^2 J}{\partial \theta^{(i)}_{j} \partial \theta^{(k)}_{l}} \
&= \frac{1}{N} \sum_{n=0}^{N-1} p_{n, i} \left( \delta_{i, k} - p_{n, k}\right) \Phi_{n, j} \Phi_{n, l}
{} + \frac{\lambda}{N} \delta_{(i,j), (k,l)}
\end{align}
$$

最初に、第1項に対応する$(C, M, C, M)$ arrayを計算することにします(正則化に由来する項はただの定数×単位行列なので、$(MC, MC)$arrayに直した後で足します。)。
計算にあたってはnumpy broadcastingとnp.einsumを利用します。

• まず、$p_{n, i} \left( \delta_{i, k} - p_{n, k}\right)$をtmparr_Aとして計算しておきます。具体的にはtmparr_A[n, i, k]$ = p_{n, i} \left( \delta_{i, k} - p_{n, k} \right)$とします。
• 次に、$\Phi_{n, j} \Phi_{n, l}$をtmparr_Bとして計算します。具体的には、tmparr_B[n, j, l] $= \Phi_{n, j} \Phi_{n, l}$とします。
• 最後に、この2つのarrayをnp.einsumで組合せて結果を得ます。

なお、Hessianの計算には、$\mathcal{O}(NM^2C^2)$の計算量がかかります。というのも、$(CM)^2$個ある行列の各要素の計算に$\mathcal{O}(N)$かかるからです。

def hess_cost_function(thtvec, Phi, T, lam):
    '''
    This function calculates the Hessian of the cost function 
    
    Parameters
    ----------
    See the help for cost_function
    
    Returns
    ----------
    H : 2-D array
        The hessian
    '''
    N, M = np.shape(Phi)
    C = np.shape(T)[1]
    Theta = np.reshape(thtvec, (C, M))
    P = calcPmat(Theta, Phi)
    tmparr_A = np.reshape(P, (N, C, 1)) * ( np.reshape(np.identity(C), (1, C, C)) - np.reshape(P, (N, 1, C)) )
    tmparr_B = np.reshape(Phi, (N, M, 1)) * np.reshape(Phi, (N, 1, M)) 
    H_tmp = np.einsum('nik,njl->ijkl', tmparr_A, tmparr_B)/N
    H = np.reshape(H_tmp, (C*M, C*M)) + lam/N*np.identity(C*M)
    return H

2.2.4 Hessianとベクトルの積

1.4.2節で述べたように、共役勾配法(CG法)を用いて$H^{-1}(x)\nabla f(x)$を計算する場合、Hessianそのものは必要なく、「与えられたベクトルとHessianの積を計算する函数」があれば十分になります(詳しくは https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.cg.html 参照。)。
後で用いるために、ここではそのような「積を与える函数」を与える函数を定義しておきます。この「積を与える函数」における計算は、Hessianそのものの計算に比べて大幅に短い時間ですみます(この節の最後に評価します。)。

まず、「与えられたベクトル」を(C*M,) array vとします。また、$(V_{i,j})_{i= 0, \dots, C-1, j = 0, \dots, M-1}$をvの行列(2-D array)表記とします。入力のベクトルとHessianの積は、次で与えられます：
$$
\begin{align}
&{} \sum_{k,l} H_{(i, j), (k, l)} V_{k, l} \
&= \sum_{k,l} \frac{\partial^2 J}{\partial \theta^{(i)}_{j} \partial \theta^{(k)}_{l}} V_{k, l} \
&= \frac{1}{N} \sum_{n=0}^{N-1} \sum_{k,l} p_{n, i} \left( \delta_{i, k} - p_{n, k}\right) \Phi_{n, j} \Phi_{n, l} V_{k, l} + \frac{\lambda}{N} \sum_{k,l} \delta_{(i,j), (k,l)} V_{k,l} \
&= \frac{1}{N} \sum_{n=0}^{N-1} \left( \sum_{l} p_{n, i} \Phi_{n, j} \Phi_{n, l} V_{i, l} - \sum_{k,l} p_{n, i} p_{n, k} \Phi_{n, j} \Phi_{n, l} V_{k, l} \right) + \frac{\lambda}{N} V_{i,j} \
&= \frac{1}{N} \sum_{n=0}^{N-1} p_{n, i} \Phi_{n, j} \left( \sum_{l} \Phi_{n, l} V_{i, l}
{} - \sum_{k,l} p_{n, k} \Phi_{n, l} V_{k, l} \right) + \frac{\lambda}{N} V_{i,j}
\end{align}
$$

複雑に見えますが、step by stepで計算していきましょう。計算の途中の結果をtmpに保存するようにしておきます。

• まず、tmp[n, k] $= \sum_{l=0}^{M-1} \Phi_{n, l} V_{k, l}$ = (Phi @ V.T) [n, k]とします。 添え字を外して考えると、tmp = Phi @ V.Tで求まります。PhiとVの形はそれぞれ(N, M)と(C, M)なので、この計算には$\mathcal{O}(NMC)$の時間がかかります.
• 次に、$ \sum_{l} \Phi_{n, l} V_{i, l}- \sum_{k,l} p_{n, k} \Phi_{n, l} V_{k, l}$ = tmp[n, i] $- \sum_{k=0}^{C-1}$ P[n, k] * tmp[n, k]を計算し、これをtmp[n, i]とおきます。この計算には$\mathcal{O}(NC^2)$の時間がかかります(先に積をとってから和を取る場合)。
• これを用いて、次の計算を行います：$\sum_{n=0}^{N-1} p_{n, i} \Phi_{n, j} \left( \sum_{l} \Phi_{n, l} V_{i, l}- \sum_{k,l} p_{n, k} \Phi_{n, l} V_{k, l} \right)$ = $\sum_{n=0}^{N-1} $P[n, i] * Phi[n, j] * tmp[n, i] = $\sum_{n=0}^{N-1}$ (P*tmp)[n, i] * Phi[n, j] = ((P*tmp).T @ Phi)[i,j]。この計算には、$\mathcal{O}(NMC)$の時間がかかります。

まとめると、Hessianと与えられたベクトルの積を計算するのにかかる時間は$\mathcal{O}(NC(M+C))$となり、Hessianそのものを求めるための計算量より大幅に小さくなります⁹。

また、コードは以下のようになります。後の線型方程式ソルバーの仕様上、scipy.sparse.linalg.LinearOperatorオブジェクトを返すようにしておきます。

from scipy.sparse.linalg import LinearOperator

def hessp_cost_function(thtvec, Phi, T, lam):
    '''
    This function returns the function representing the product of Hessian and a given vector
    
    Parameters
    ----------
    See the help for cost_function
    
    Returns
    ----------
    A : scipy.sparse,linalg.LinearOperator object
        A linear operator representing the Hessian of the cost function
    '''
    N, M = np.shape(Phi)
    C = np.shape(T)[1]
    Theta = np.reshape(thtvec, (C, M))
    P = calcPmat(Theta, Phi)
    def hp(v):
        V = np.reshape(v, (C, M))
        tmp = Phi @ V.T
        tmp = tmp - np.reshape(np.sum(P*tmp, axis=1), (N, -1))
        product = ((P*tmp).T @ Phi)/N + lam/N*V
        product = product.flatten()
        return product
    A = LinearOperator((C*M, C*M), matvec=hp)
    return A

2.3 最小化

ここまでで損失函数は用意できたので、これを最小化する函数を書きます。

2.3.1 最急降下法

最急降下法では、次の更新を繰り返して(近似的)最適解を得ます。
$$
\begin{align}
x_{n+1} = x_n - \alpha \nabla f(x_n).
\end{align}
$$
繰り返しを停止する基準としては、「函数$f$の値の変化が事前に与えた値ftolを下回れば停止する」をとります¹⁰。

def minimize_GD(func, x0, grad, alpha=0.01, maxiter=1e4, ftol=1e-5):
    '''
    This function minimizes the given function using gradient descent method
    
    Parameters
    ----------
    func : callable
        Function to be minimized (real-valued)
    x0 : 1-D array
        Initial value of the variable
    grad : callable
        The gradient of func (returns 1-D array)
    alpha : float
        Learning rate
    maxiter : int
        Maximum number of iteration
    ftol : float
        The threshold for stopping criterion. If the change of the value of function is smaller than this value, the iteration stops.
    
    Returns
    ----------
    result : dictionary
        result['x'] ... variable, result['nit'] ... the number of iteration, result['func']...the value of the function, result['success']... whether the minimization is successful or not
    '''
    x = x0
    nit = 0
    while nit < maxiter:
        xold = x
        x = x - alpha*grad(x)
        nit += 1
        if abs(func(x) - func(xold)) < ftol:
            break
    success = (nit < maxiter)
    return {'x': x, 'nit':nit, 'func':func(x), 'success':success}

2.3.2 Newton法

1.4.2節でも述べたように、Newton法の更新式は
$$
\begin{align}
x_{n+1} = x_{n} - \left[H(x_n)\right]^{-1} \nabla f(x_n)
\end{align}
$$
で与えられます。

そして、これまた1.4.2節でも述べたように、第2項$\left[H(x_n)\right]^{-1} \nabla f(x_n)$を計算する方法は大きく分けて2つあります：

• 1つは、numpy.linalg.solve( https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html )を用いて線型方程式$H(x_n) y = \nabla f(x_n)$を$y$について(数値的に)「厳密に」解く手法です。
• もう1つは、先ほどの線型方程式を繰り返しを用いて近似的に解く手法です。ここでは、共役勾配法scipy.sparse.linalg.cg( https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.cg.html )を用いることにします。

既に述べたように、後者の手法を用いる場合はHessianそのものは必要ではなく、「Hessianとベクトルの積を計算する函数」を与えれば十分です。2.2.3節と2.2.4節で見たように、後者の計算量はHessianを求める計算量より低いため、計算が速くなることが期待されます¹¹。比較のため、両方の手法を実装します。

まず、前者の手法で更新を行う函数です。

def minimize_newton(func, x0, grad, hess, maxiter=1e4, ftol=1e-5):
    '''
    This function minimizes the given function using Newton method
    
    Parameters
    ----------
    func : callable
        Function to be minimized (real-valued)
    x0 : 1-D array
        Initial value of the variable
    grad : callable
        The gradient of func (returns 1-D array)
    hess: callable
        The hessian of func (returns 2-D array)
    maxiter : int
        Maximum number of iteration
    ftol : float
        The threshold for stopping criterion. If the change of the value of function is smaller than this value, the iteration stops.
    
    Returns
    ----------
    result : dictionary
        result['x'] ... variable, result['nit'] ... the number of iteration, result['func']...the value of the function, result['success']... whether the minimization is successful or not
    '''
    x = x0
    nit = 0
    while nit < maxiter:
        xold = x
        x = x - np.linalg.solve(hess(x), grad(x))
        nit += 1
        if abs(func(x) - func(xold)) < ftol:
            break
    success = (nit < maxiter)
    return {'x': x, 'nit':nit, 'func':func(x), 'success':success}

次に、共役勾配法を用いたものです。これをNewton-CG法と呼ぶことにします。

from scipy.sparse.linalg import cg

def minimize_ncg(func, x0, grad, hessp, maxiter=1e4, ftol=1e-5):
    '''
    This function minimizes the given function using Newton-CG method
    
    Parameters
    ----------
    func : callable
        Function to be minimized (real-valued)
    x0 : 1-D array
        Initial value of the variable
    grad : callable
        The gradient of func (returns 1-D array)
    hessp: callable
        A function which returns scipy.sparse.linalg.LinearOperator object representing Hessian
    maxiter : int
        Maximum number of iteration
    ftol : float
        The threshold for stopping criterion. If the change of the value of function is smaller than this value, the iteration stops.
    
    Returns
    ----------
    result : dictionary
        result['x'] ... variable, result['nit'] ... the number of iteration, result['func']...the value of the function, result['success']... whether the minimization is successful or not
    '''
    x = x0
    nit = 0
    while nit < maxiter:
        xold = x
        x = x - cg(A=hessp(x), b=grad(x))[0]
        nit += 1
        if abs(func(x) - func(xold)) < ftol:
            break
    success = (nit < maxiter)
    return {'x': x, 'nit':nit, 'func':func(x), 'success':success} 

2.3.4 まとめたもの

以上をまとめた最小化の函数を書いておきます(methodから具体的なアルゴリズムを選択できるようにしてあります。)。

def minimize(func, x0, method, grad, hess=None, hessp=None, alpha=None, maxiter=1e4, ftol=1e-5):
    '''
    This function minimizes the given function using the method speicified by user.
    
    Parameters
    ----------
    func : callable
        Function to be minimized (real-valued)
    x0 : 1-D array
        Initial value of the variable
    method : string
        One of 'gd', 'newton', 'newton-cg'.
    grad : callable
        The gradient of func (returns 1-D array)
    hess : callable
        The hessian of func (returns 2-D array). Required when method='newton'
    hessp: callable
        A function which returns scipy.sparse.linalg.LinearOperator object representing Hessian. Required when method='newton-cg'
    alpha : float
        Learning rate
    maxiter : int
        Maximum number of iteration
    ftol : float
        The threshold for stopping criterion. If the change of the value of function is smaller than this value, the iteration stops.
    
    Returns
    ----------
    result : dictionary
        result['x'] ... variable, result['nit'] ... the number of iteration, result['func']...the value of the function, result['success']... whether the minimization is successful or not
    '''
    if method == 'gd':
        if alpha is None:
            print("Gradient descent needs a learning rate parameter alpha.")
            return
        else:
            return minimize_GD(func, x0, grad, alpha, maxiter, ftol)
    elif method == 'newton':
        if hess is None:
            print("Newton-Raphson method requires the hessian.")
            return
        else:
            return minimize_newton(func, x0, grad, hess, maxiter, ftol)
    elif method == 'newton-cg':
        if hessp is None:
            print("Newton-CG method requires the hessian-product function.")
            return
        else:
            return minimize_ncg(func, x0, grad, hessp, maxiter, ftol)
    else:
        print("method should be one of gd, newton, newton-cg")

2.4 ロジスティック回帰分類器

準備が長かったですが、ようやく分類器にたどりつきました。
ロジスティック回帰を行うためのクラスLogisticClfを作ります。
次の属性を持たせることにします:

• 分類するラベルの種類の個数$C$ (self.C)
• 正則化パラメーター$\lambda$ (self.lam)
• 重みパラメーター$\theta$(self.Theta)

また、次のメソッドを持たせます。

• コンストラクタ (__init__)
• 学習を行うメソッド(fit さきほどの損失函数と勾配を用いて、fittingを行います。)
• 入力に対して、各ラベルの確率を返すメソッド(predict_proba)
• 入力に対して、予測されるラベルを返すメソッド(predict)

また、minimizeのためのx0(変数の初期値)は、ここでは全てが0の配列にしておきます(凸最適化なので、初期値はあまり気にする必要はないかと思います。)

class LogisticClf:
    def __init__(self, C, lam):
        self.C = C  # the number of labels
        self.lam = lam #regularization parameter
        self.Theta = None
            
    def fit(self, X, y, method='NR', alpha=None, maxiter=1e4, ftol=1e-5, show_message=False):
        '''
        Parameters
        ----------
        X : 1-D or 2-D numpy array
            (N,) or (N, d) array, representing the training input data
        y : 1-D numpy arra
            (N,) array, representing training labels
        '''
        Phi = calcPhimat(X)
        T = calcTmat(y, self.C)
        N, M = np.shape(Phi)
        
        tht0 = np.zeros(M*self.C)
        time_start = time.time()
        result = minimize(func=lambda x : cost_function(x, Phi, T, self.lam), 
                          x0=tht0, 
                          method=method,
                          grad=lambda x : grad_cost_function(x, Phi, T, self.lam),
                          hess=lambda x : hess_cost_function(x, Phi, T, self.lam),
                          hessp = lambda x : hessp_cost_function(x, Phi, T, self.lam),
                          alpha=alpha,
                          maxiter=maxiter,
                          ftol=ftol
                         )
        time_end = time.time()
        if show_message:
            print(result['success'])
            print(f"nit: {result['nit']}")
            print(f"calcualtion time : {time_end - time_start}seconds")
        self.Theta = np.reshape(result['x'], (self.C, M))
        
    def predict_proba(self, X):
        '''
        Parameters
        ----------
        X : 1-D or 2-D numpy array
            (N,) or (N, d) array, representing the training input data
        Returns
        ----------
        proba : 2-D numpy arra
            (len(X), self.C) array, where proba[n, c] represents the probability that the n-th instance belongs to c-th class
        '''
        return calcPmat(self.Theta, calcPhimat(X))
        
    def predict(self, X):
        '''
        Parameters
        ----------
        X : 1-D or 2-D numpy array
            (N,) or (N, d) array, representing the training input data
        Returns
        ----------
        classes : 1-D numpy arra
            (len(X), ) array, where classes[n] represents the predicted class to which the n-th instance belongs
        '''
        tmp = self.predict_proba(X)
        return np.argmax(tmp, axis=1 )

3 実験

ということで、ようやく実験です。
今回は簡単な2次元のデータと、MNISTの手書き文字データを用います。

3.1 Toy data

まずは、簡単な2次元のデータを用います。$N=200, d=2, C=3$とします。

X, y = datasets.make_blobs(n_samples=200, n_features=2, centers=3)

このデータで、先ほどのLogisticClfを学習させます。
最急降下法、Newton法、Newton-CG法のそれぞれを用いて、学習の結果を比べます。
以下では、カラーバーの数字の0, 1, 2が予測されたラベルを表します。

まず、勾配降下法の結果を示します。

success : True
nit: 693
calcualtion time : 0.241013765335083seconds

次に、通常のNewton法です

success : True
nit: 8
calcualtion time : 0.009000778198242188seconds

最後に、Newton-CG法の結果です。

success : True
nit: 8
calcualtion time : 0.01500082015991211seconds

まず、いずれも、綺麗に分類できているのが見て取れるかと思います。なお、この図では決定境界が直線になっていますが、これは我々が基底関数$\phi$として、"線型"なものを用いたことに由来します。非線形な基底函数を用いた場合、一般に境界は直線になるとは限りません。

また、繰り返し回数という意味でも実際の計算時間という意味でも、Newton法とNewton-CG法が最急降下法に比べて速く収束していることが分かります。
これは、今のデータでは$\theta$の次元が$MC=9$と小さいことによるものと考えられます(次元が低いため、Hessianを扱うコストは大きくなりません。)。データ数$N$を1000程度まで増やしても同様の振る舞いが見られました。

3.2 手書き文字データ

次に、より次元の高いデータで試してみましょう。ここでは、MNISTの手書き文字データを用います。

digits = datasets.load_digits()
dat_train, dat_test, label_train, label_test = train_test_split(digits.data, digits.target, test_size=0.25)

今度はデータ数と次元があがります。

Training data : 1347
Test data : 450
Input dimension : 64

特にこの場合、$MC = (d+1)C = 650$となり、先ほどの例より大幅に大きくなります。
まずは正則化パラメーターを固定($\lambda=1.0$)して、solverを変えたときの挙動を見てみましょう。

最急降下法

success : True
nit: 2104
calcualtion time : 3.5881898403167725seconds
train accuracy score: 0.9977728285077951
test accuracy score: 0.9711111111111111

Newton法

success : True
nit: 9
calcualtion time : 28.27373433113098seconds
train accuracy score: 1.0
test accuracy score: 0.9577777777777777

Newton-CG法

success : True
nit: 9
calcualtion time : 1.4760844707489014seconds
train accuracy score: 1.0
test accuracy score: 0.9577777777777777

Newton法とNewton-CG法で繰り返し回数が少ないのは先ほどと同様です。しかし、今度は先ほどとは異なり、Newton法では学習に時間がかかってしまっています。これは、$MC$が大きくなってHessianを求めるコストが高くなったためと考えられます。一方、Newton-CG法ではHessianを直接扱う必要がないためか、非常に速く収束しています¹²。

興味深いことに、Newton法とNewton-CG法では訓練精度が1.0に達しているのに対し、最急降下法の場合は訓練精度が1.0に達していません。一方で、テスト精度は最急降下法で学習したときのほうが高くなっています。これは、early stoppingの一例ではないかと思います。

次に、正則化パラメーター$\lambda$の値も変えて試してみましょう(以下では計算が速いNewton-CG法を用います。)。結果は以下の通りでした。

success : True
nit : 11
calcualtion time : 2.740156650543213seconds
lambda = 0.1
train accuracy score: 1.0
test accuracy score: 0.9533333333333334

success : True
nit : 9
calcualtion time : 0.9400537014007568seconds
lambda = 1.0
train accuracy score: 1.0
test accuracy score: 0.9577777777777777

success : True
nit : 8
calcualtion time : 0.4060232639312744seconds
lambda = 10.0
train accuracy score: 0.9992576095025983
test accuracy score: 0.9666666666666667

だいたい正しく分類できていることに加え、$\lambda$が大きくなるにつれて汎化誤差が下がっていることも見て取れます。正則化パラメーターを大きくすることは、biasを上げる代わりにvarianceを下げることに対応するので、この振る舞いは直感的にも納得のいくものかと思います。

4 まとめ

今回は、分類問題に対するロジスティック回帰を実装しました。
特に、損失函数の最小化を数値的に行うにあたっていくつかの手法があることを概観し、3つの最適化手法を実装しました。
数値実験を通じて、

• パラメーターの次元が大きいときでもNewton-CG法を用いると高速に学習を行えることを見ました。
• 正則化パラメーターが結果に与える影響についても観察しました。

旧版の記事では何も考えずに、2階微分を用いずにscipy.optimize.minimizeを用いたため計算が遅くなっていました。それをHessianを用いたNewton-CG法に書き換えることにより、学習の大幅な高速化ができました¹³。

1. 本の$K$を$C$、$k$を$c$、$y_{n,k}$を$p_{n,c}$, $w_k$を$\theta^{(c)}$と表しています。 ↩

2. 肝心の「なぜ勾配だけでできる勾配降下法ではなく、わざわざHessianを用いてNewton法を用いて最小化を行うのか」についてはPRMLの中ではあまり説明されていないかと思います。この点については、この後この記事で軽く触れたいと思います。 ↩

3. 線型回帰の場合は手で簡単に求められたのとは対照的ですね。 ↩

4. 損失函数が凸であることは、Hessianが正定値であることから分かります。 ↩

5. 英語ではgradient descentの方がよく聞くので、素直に訳すなら「勾配降下法」な気がします。ちなみにstochastic gradient descentの和訳としては「確率的勾配降下法」が一般的なようなので、やや不思議ですね。 ↩

6. 学習率をステップ毎に自動で決める手法もありますが、ここでは取り扱いません。気になる方は https://www.slideshare.net/pecorarista/steepest-descent などをご覧ください。 ↩

7. 実際は逆行列を求めずにLU分解を経由して線型方程式を解く方が速くなります。詳細は、伊理 正夫、藤野 和建(1985)「数値計算の常識」共立出版の第5章「逆行列よさようなら」などをご参照ください。 ↩

8. 線型方程式に対する共役勾配法は係数行列が正定値の場合に使えます。ここで考えているロジスティック回帰の損失函数のHessianは正定値行列なので、問題なく利用可能です。 ↩

9. ただし、Hessianを計算してから線型方程式を解く場合はHessianを計算するのは一度ですみますが、共役勾配法で線型方程式を解く場合はHessianとベクトルの積を複数回計算する必要があります。そのため、実際の計算時間は共役勾配法の繰り返し回数にも依存することになると思います。 ↩

10. 他にも、「$x$の変化が事前に与えた値より小さくなれば停止する」「$\nabla f(x)$が小さくなれば停止する」などの基準がありえます。ここではひとまず、$f$の値の変化を基準にとってみます。 ↩

11. 本当は、共役勾配法でどの程度の回数の繰り返しが行われるかを評価しなければ、全体で見た計算時間の見積もりになりませんが。。。 ↩

12. 正確には、「Hessianを直接扱う必要がない」だけでなく、「numpy.linalg.solveのalgorithmよりもscipy.linalg.sparse.cgのアルゴリズムが速い」による寄与が大きいという可能性もあります。これを比較するために、あえて直接Hessianを扱ってNewton-CG法を使ってみたところ、通常のNewton法とほぼ同程度の学習時間がかかりました。従って、ここでは線型方程式を解くアルゴリズムそのものの速さよりも、Hessianを直接扱わなくても済むことの寄与が大きいと考えられます。 ↩

13. やっぱり横着はよくないですね。今回の記事の執筆を通じて、「何が原因かをきちんと分析し、自分で考える/書くこと」の大切さも認識しなおしたのでした ↩