1次元ロジスティック回帰¶

1次元のロジスティック回帰は，

1. $d$次元ベクトル${\bf x}$を$f({\bf x}) = {\bf w}^{\rm T}{\bf x} + b$によって1次元に写像

2. それをシグモイド関数で$(0, 1)$区間へ写像：$s({\bf x}) \equiv \sigma(f({\bf x}))$

3. これを用いてラベル$c(0\ or\ 1)$に対する尤度を定義：

   $p(y=c \mid {\bf x};{\bf w},b) = s({\bf x})^c(1-s({\bf x}))^{1-c}$

   （${\bf w},b$は確率変数ではなくロジスティックモデルのパラメータなので，$;$によって区切られた形で示されている）

4. 与えられた$N$個のデータ点${\bf x}_n$とそのラベル値$c_n$について，各データ点が全て独立であるとの仮定のもとで以下の尤度関数を最大化する${\bf w},b$を求める

   $$\mathcal{l}({\bf w},b) = \prod_{n=1}^N p(y \mid {\bf x};{\bf w},b) = \prod_{n=1}^N s({\bf x}_n)^{c_n}(1-s({\bf x}_n))^{1-c_n}$$

という手順によって${\bf x}$から2値ラベル$c$を予測する写像を得る（すなわち2クラス分類を行う）手法である．以下に書く手順の詳細を述べる．

シグモイド関数は$\sigma(a)=\frac{1}{1+\exp(-a)}$のような関数で，以下のような形をしている．

これは$(-\infty, \infty)$の区間の実数を$(0, 1)$に写像する関数なので，あるデータ点${\bf x}$を$f({\bf x}) = {\bf w}^{\rm T}{\bf x} + b$によって1次元に写像したあと，シグモイド関数に通した値$\sigma(f({\bf x}))$は，$(0,1)$の値をとる確率値であると考えることができるようになる．

今，データ点${\bf x}$にはラベルとして0か1のどちらかが対応しているとする．このとき，${\bf x}$のラベルが1である確率と0である確率は，それぞれ

$$p(y=1 \mid {\bf x};{\bf w},b) = \sigma(f({\bf x})) \equiv s({\bf x})$$ $$p(y=0 \mid {\bf x};{\bf w},b) = 1 - \sigma(f({\bf x})) \equiv 1 - s({\bf x})$$

と表すことができる．ラベルを$c(0\ or\ 1)$で表せば，データ点${\bf x}$が与えられたときのラベルの尤度は

$$p(y=c \mid {\bf x};{\bf w},b) = s({\bf x})^c(1 - s({\bf x}))^{1-c}$$

とまとめて書くことができる．今，$N$個のデータ点があるとして，それぞれが全て独立であると考えると，データ点の集合${\bf x}_n(n=1,\dots,N)$とラベルの集合$c_n(n=1,\dots,N)$に対する，このモデルのパラメータである${\bf x},b$の尤度は各データの尤度の積で表せるので，

$$\mathcal{l}({\bf w},b) = \prod_{n=1}^N p(y \mid {\bf x};{\bf w},b) = \prod_{n=1}^N s({\bf x}_n)^{c_n}(1-s({\bf x}_n))^{1-c_n}$$

と定義できる．この最大化を考えるのだが，まず$(0,1)$区間の値の積をデータ数分だけ計算しようとすると，浮動小数点のオーバーフローが起こりやすいため，対数をとって和の形で表すことにする．対数関数は単調増加であるので，対数をとったとしても大小関係は変わらない．また，実装上の容易さの観点から，符号を反転することで最大化ではなく最小化問題に置き換えることにする．すると，目的関数は以下のように書き下せる．

$$\begin{eqnarray} \def\sx{s({\bf x}_n)} \mathcal{L}({\bf w},b) &=& -\sum_{n=1}^N \ln \left\{ \sx^{c_n}(1 - \sx)^{1 - c_n} \right\} \\ &=& -\sum_{n=1}^N \left\{ c_n \ln \sx + (1 - c_n) \ln (1 - \sx) \right\} \end{eqnarray}$$

これを負の対数尤度（negative log likelihood）関数と呼ぶ．また，モデルの予測のはずれ具合を表すので，損失関数（loss function）とも呼ばれる．これを最小化するパラメータ${\bf w},b$は，この関数をパラメータそれぞれで偏微分して0とおけば求まるような気がする．シグモイド関数の導関数が

$$\sigma(a) = \frac{1}{1 + \exp(-a)}$$ $$\frac{d \sigma}{d a} = \frac{\exp(-a)}{(1 + \exp(-a))^2} = \frac{1}{1 + \exp(-a)}\frac{\exp(-a)}{1 + \exp(-a)} = \sigma(a)\frac{1 + \exp(-a) - 1}{1 + \exp(-a)} = \sigma(a)\left( 1 - \sigma(a) \right)$$

となることに注意して，パラメータ${\bf w},b$それぞれにおける偏導関数を求めてみる．

${\bf w}$における偏導関数¶

$$\frac{\partial \mathcal{L}}{\partial {\bf w}} = -\sum_{n=1}^N \left\{ c_n \frac{\partial \ln\sigma(f({\bf x}_n))}{\partial f({\bf x}_n)}\frac{\partial f({\bf x}_n)}{\partial {\bf w}} + (1 - c_n) \frac{\partial \ln(1 - \sigma(f({\bf x}_n)))}{\partial f({\bf x}_n)} \frac{\partial f({\bf x}_n)}{{\bf w}} \right\}$$ $$\begin{eqnarray} \def\xn{{\bf x}_n} \frac{\partial \mathcal{L}}{\partial {\bf w}} = -\sum_{n=1}^N \xn (c_n - \sx) \end{eqnarray}$$

$b$における偏導関数¶

${\bf w}$と同様に

$$\begin{eqnarray} \def\xn{{\bf x}_n} \frac{\partial \mathcal{L}}{\partial b} = -\sum_{n=1}^N (c_n - \sx) \end{eqnarray}$$

パラメータの学習¶

以上の2つの導関数

$$\frac{\partial \mathcal{L}}{\partial {\bf w}} = -\sum_{n=1}^N \xn (c_n - \sx)$$ $$\frac{\partial \mathcal{L}}{\partial b} = -\sum_{n=1}^N (c_n - \sx)$$

は，シグモイド関数$\sx$を含んでいるため，$=0$とおいて解析的に${\bf w}$と$b$の値を決定することができない．そこで，この偏導関数を使って，勾配降下法（gradient descent）によるパラメータ${\bf w},b$の更新を行うことで，損失関数$\mathcal{L}({\bf w},b)$の最小化を行う．

このための更新式は以下である．

$${\bf w} \leftarrow {\bf w} - \eta\frac{\partial \mathcal{L}}{\partial {\bf w}}$$ $$b \leftarrow b - \eta\frac{\partial \mathcal{L}}{\partial b}$$

ただし，$\eta$は学習率とよばれ，一回の更新でどのくらいパラメータを変化させるかという度合いを表す．これが大きいと振動しながら最適解へ近づくことになり，小さいと収束に時間を要することになる．

確率的勾配降下法（stochastic gradient descent: SGD）¶

上記の更新式によるパラメータの更新では，勾配を計算するために毎回，全てのデータ点に渡って和をとる必要がある．これは，データセットが大きいな場合は非常に計算コストが大きくなってしまう．

そこで，データセットの中から部分的にいくつかのデータ点を選びとって，それらに対してだけ上記の和を計算して勾配とし，パラメータを更新する方法がよくとられる．これを確率的勾配降下法（stochastic gradient descent: SGD）と呼ぶ．どのデータ点部分集合を選択するかを確率的に選ぶためこう呼ばれる．このとき，一回の更新に用いられるデータセットの部分集合のことをミニバッチ（minibatch）と呼ぶことがある．

また，勾配の計算の際，上式のようにただ総和をとるのではなく，平均を計算することで，通常の勾配降下法の場合はデータセットの大きさに，確率的勾配降下法の場合はミニバッチの大きさに，勾配の大きさが影響されないようにすることができる．このようにすることで学習率$\eta$の決定も容易になる．

実装¶

まず，2つのクラスから生成されたデータ点群が与えられたとする．それぞれのデータ点は中心を$(0,0),(5,5)$にもつ分散1の正規分布から100個ずつサンプルされた点である．

これをひとつのデータ行列にまとめ，対応するラベルベクトルを作成する．ここでは，クラス1に属するデータ点にはラベル値$c=0$を，クラス2に属するデータ点にはラベル値$c=1$を対応させる．

シグモイド関数と，データ点群行列（$N \times d$）からクラス尤度を算出する関数p_y_given_x(x)を定義する．

データ点群と正解ラベルベクトルをとり，勾配を返す関数を定義する．

パラメータを適当に初期化したのち，データセットを10個ずつのミニバッチに分けて，確率的勾配降下法を用いてパラメータを最適化する．

上のコードではミニバッチ最適化をさらに10回繰り返している．繰り返しごとにエラーを記録し，最後にその変化をグラフ化している．パラメータを更新するたびに誤差が減っていっていることが分かる．

学習されたパラメータによる直線：${\bf w}^{\rm T}{\bf x} + b = 0$の図示¶

多次元ロジスティック回帰¶

上記の例では1次元のロジスティック回帰を考えた．これは2クラス分類問題を考えるのに利用できることが分かった．では，$K$クラス分類を考える場合はどうなるだろうか．この場合，出力が$K$次元のクラス所属度確率（class-membership probability）ベクトルとなるようにソフトマックス関数というのを用いる．ソフトマックス関数はシグモイド関数の多変量版であり，$K$次元ベクトル${\bf a}$の$k$番目の要素を$a_k(k=1,\dots,K)$と書くことにすると，

$$\def\softmax{\boldsymbol\pi} \pi(a_k) = \frac{\exp(a_i)}{\sum_{j=1}^K \exp(a_j)}$$

と定義される．以降，便宜上$(\pi(a_1),\dots,\pi(a_K))^{\rm T}$というベクトルを$\softmax({\bf a})$と書くことにする．

ソフトマックス関数にベクトルを入力すると，出力される値は同次元のベクトルとなるが，各要素の合計値が必ず1となるようなベクトルとなる．すなわち

$$\sum_{k=1}^K \pi(a_k) = 1$$

が常に成り立つ．よって，出力されるベクトルは離散確率分布であると考えることができる．

$K$クラス分類問題において，$D$次元ベクトル${\bf x}$と，それに対応する正解クラス$c_k$が与えられたとき，多次元ロジスティック回帰モデルは以下のような手順で計算を行う．

1. $D$次元ベクトルであるデータ点${\bf x}$を$K$次元に写像：$f({\bf x}) = {\bf W}^{\rm T}{\bf x} + {\bf b}$（${\bf W}$は$D \times K$の重み行列，${\bf b}$は$K$次元ベクトル）

2. ソフトマックス関数に通してクラス所属度確率ベクトルを計算：${\bf y} = \softmax(f({\bf x}))$

3. ${\bf x}$が所属するクラスが$c_k$のとき，正解データとして$k$個目の要素だけが1でその他の要素が0であるような教師ベクトル${\bf t}$（例：$(0,\dots,1,\dots,0)^{\rm T}$）が用意されているとする．これを用いて，データ点が与えられたときのパラメータ${\bf W},{\bf b}$の尤度を以下のように定義する：

   $$p({\bf y}={\bf t} \mid {\bf x};{\bf W},{\bf b}) = \prod_{k=1}^K y_k^{t_k}$$ このとき，${\bf t}$の形を念頭に，$y_k^{t_k}$に着目すると，これは$k$が${\bf x}$が所属しているクラスのインデックスに等しいときだけ$t_k=1$となって$y_k$に等しくなり，それ以外の$k$のときは$t_k=0$より$y_k^0=1$となる仕組みをとっていることが分かる．これを全ての要素$k=1,\dots,K$に渡って掛けあわせているので，結果的にクラス所属度確率ベクトルの要素のうち，正解クラスへの所属度確率だけを取り出す操作に等しい．

4. $N$個のデータ点がそれぞれ独立に観測されるとして，データセット全体に渡るパラメータの尤度を以下のように定義する：

   $$\mathcal{l}({\bf W},{\bf b}) = p({\bf Y}={\bf T} \mid {\bf X};{\bf W},{\bf b}) = \prod_{n=1}^N \prod_{k=1}^K y_{nk}^{t_{nk}}$$ ただし

   • ${\bf Y}$は$N$個のデータのクラス所属度ベクトル（列ベクトル）を行方向に並べたもの$({\bf y}_1,\dots,{\bf y}_N)$であり，

     $${\bf Y} = \left( \begin{array}{ccc} y_{11} & & y_{N1} \\ \vdots & & \vdots \\ y_{1k} & \cdots & y_{Nk} \\ \vdots & & \vdots \\ y_{1K} & & y_{NK} \end{array} \right) = \softmax({\bf W}^{\rm T}{\bf X} + {\bf B})$$

     という$K \times N$行列

   • ${\bf T}$は$N$個のデータの教師ベクトル（列ベクトル）を行方向に並べたもの$({\bf t}_1,\dots,{\bf t}_N)$であり，

     $${\bf T} = \left( \begin{array}{ccc} t_{11} & & t_{N1} \\ \vdots & & \vdots \\ t_{1k} & \cdots & t_{Nk} \\ \vdots & & \vdots \\ t_{1K} & & t_{NK} \end{array} \right)$$

     という$K \times N$行列

   • ${\bf X}$はデータ行列（デザイン行列とも呼ばれる）であり，データベクトル（列ベクトル）を行方向に並べたもの$({\bf x}_1,\dots,{\bf x}_N)$で，

     $${\bf X} = \left( \begin{array}{ccc} x_{11} & & x_{N1} \\ \vdots & & \vdots \\ x_{1d} & \cdots & x_{Nd} \\ \vdots & & \vdots \\ x_{1D} & & x_{ND} \end{array} \right)$$

     という$D \times N$行列

   • ${\bf B}$は全く同じ$K$次元バイアスベクトル${\bf b}$を行方向に$N$個コピーした$K \times N$行列

     $${\bf B} = \left( \begin{array}{ccc} b_{1} & & b_{1} \\ \vdots & & \vdots \\ b_{k} & \cdots & b_{k} \\ \vdots & & \vdots \\ b_{K} & & b_{K} \end{array} \right)$$

     という$K \times N$行列になる．

5. 上記の尤度を最大にするように${\bf W},{\bf b}$を求めることで学習を行う

勾配計算の準備¶

尤度$\mathcal{l}({\bf W},{\bf b})$を例によって負の対数尤度に変換すると，

$$\mathcal{L}({\bf W},{\bf b}) = -\sum_{n=1}^N \sum_{k=1}^K t_{nk} \ln y_{nk}$$

となる．更新式を求めるために勾配を算出しよう．

パラメータ${\bf W},{\bf b}$に関する勾配は，

$$\frac{\partial \mathcal{L}}{\partial {\bf W}} = \left( \frac{\partial \mathcal{L}}{\partial {\bf w}_1},\dots,\frac{\partial \mathcal{L}}{\partial {\bf w}_j},\dots,\frac{\partial \mathcal{L}}{\partial {\bf w}_K} \right)$$ $$\frac{\partial \mathcal{L}}{\partial {\bf b}} = \left( \frac{\partial \mathcal{L}}{\partial b_1},\dots,\frac{\partial \mathcal{L}}{\partial b_j},\dots,\frac{\partial \mathcal{L}}{\partial b_K} \right)$$

を考えればよいから，問題となるのは$\frac{\partial \mathcal{L}}{\partial {\bf w}_j},\frac{\partial \mathcal{L}}{\partial b_j}$である．

$\frac{\partial \mathcal{L}}{\partial {\bf w}_j}$の計算¶

$$\begin{eqnarray} \frac{\partial \mathcal{L}}{\partial {\bf w}_j} &=& -\frac{\partial}{\partial {\bf w}_j} \sum_{n=1}^N \sum_{k=1}^K t_{nk} \ln y_{nk} \\ &=& -\sum_{n=1}^N \sum_{k=1}^K t_{nk} \frac{\partial \ln y_{nk}}{\partial y_{nk}} \frac{\partial y_{nk}}{\partial a_{nj}} \frac{\partial a_{nj}}{\partial {\bf w}_j} \\ &=& -\sum_{n=1}^N \sum_{k=1}^K t_{nk} \frac{1}{y_{nk}} \frac{\partial y_{nk}}{\partial a_{nj}} \frac{\partial a_{nj}}{\partial {\bf w}_j} \end{eqnarray}$$

ここで，${\bf a}_n = {\bf W}^{\rm T}{\bf x}_n + {\bf b}$とおき，この$k$個目の要素を$a_{nk}$と書いた．また，

$${\bf y}_n = \softmax({\bf a}_n) = (\pi(a_{n1}),\dots,\pi(a_{nK}))^{\rm T}$$

であるから，$y_{nk} = \pi(a_{nk})$である．微分の連鎖律の部分については，

$$\begin{eqnarray} y_{nk} &=& \pi(a_{nk}) \\ a_{nk} &=& {\bf w}_k^{\rm T}{\bf x}_n + b_k \end{eqnarray}$$

に注意する．

では，$\frac{\partial y_{nk}}{\partial a_{nj}}$および$\frac{\partial a_{nj}}{\partial {\bf w}_j}$について個別に考えていく．

$\frac{\partial y_{nk}}{\partial a_{nj}}$の計算¶

ここで添字の都合上$\frac{\partial y_{nk}}{\partial a_{nj}}$を一旦$\frac{\partial y_{nk}}{\partial a_{ni}}$と書くと，

• $i = k$のとき

  $$\begin{eqnarray} \frac{\partial y_{nk}}{\partial a_{ni}} = \frac{\partial \pi(a_{nk})}{\partial a_{ni}} &=& \frac{\exp(a_{nk})}{\sum_{j=1}^K \exp(a_{nj})} - \frac{\exp(a_{nk})^2}{\left( \sum_{j=1}^K \exp(a_{nj}) \right)^2} \\ &=& \pi(a_{nk})(1 - \pi(a_{nk})) \end{eqnarray}$$

• $i \neq k$のとき

  $$\begin{eqnarray} \frac{\partial y_{nk}}{\partial a_{ni}} = \frac{\partial \pi(a_{nk})}{\partial a_{ni}} &=& -\frac{\exp(a_{ni})\exp(a_{nk})}{\left( \sum_{j=1}^K \exp(a_{nj}) \right)^2} \\ &=& -\pi(a_{ni})\pi(a_{nk}) \\ &=& \pi(a_{nk})(0 - \pi(a_{ni})) \\ \end{eqnarray}$$

となるので，再び添字$i$を$j$に置き換えると，

$$\frac{\partial y_{nk}}{\partial a_{nj}} = y_{nk}({\bf I}_{jk} - y_{nj})$$

と書ける．ここで${\bf I}_{jk}$は$K \times K$単位行列${\bf I}$の$(j,k)$要素である．これは$j=k$のとき1であり，それ以外の場合は0をとる．

$\frac{\partial a_{nj}}{\partial {\bf w}_j}$の計算¶

$$\frac{\partial a_{nj}}{\partial {\bf w}_j} = \frac{\partial}{\partial {\bf w}_j}({\bf w}_j^{\rm T}{\bf x}_n + b_j) = {\bf x}_n$$

である．

${\bf w}_j$に関する勾配¶

上記の結果を用いると，

$$\frac{\partial \mathcal{L}}{\partial {\bf w}_j} = -\sum_{n=1}^N \sum_{k=1}^K t_{nk} {\bf x}_n ({\bf I}_{jk} - y_{nj}) = -\sum_{n=1}^N {\bf x}_n \left\{ \sum_{k=1}^K t_{nk}{\bf I}_{kj} - \sum_{k=1}^K t_{nk}y_{nj} \right\}$$

となる．ここで，${\bf I}_{kj}$は，$j=k$のときだけ1で，それ以外のときは0だから，$\sum_{k=1}^K t_{nk}{\bf I}_{kj} = t_{nj}$と書ける．また，$\sum_{k=1}^N t_{nk}y_{nj}$における$y_{nj}$は，この和をとる間は不変であり，またこのとき$t_{nk}$は$k=1,\dots,K$のいずれか1つでだけ1をとり，それ以外のときは0である教師信号だから，$\sum_{k=1}^N t_{nk}y_{nj} = y_{nj}$として問題ない．

すると，結局

$$\frac{\partial \mathcal{L}}{\partial {\bf w}_j} = \sum_{n=1}^N (y_{nj} - t_{nj}){\bf x}_n$$

というシンプルな形に変形できる．

$b_j$に関する勾配¶

$b_j$で偏微分することで変化するのは，連鎖律の最後が$\frac{\partial a_{nj}}{\partial b_j} = 1$になるだけなので，

$$\frac{\partial \mathcal{L}}{\partial b_j} = \sum_{n=1}^N (y_{nj} - t_{nj})$$

である．

更新式¶

上記の勾配を用いて更新式を定義する．

$${\bf w}_j \leftarrow {\bf w}_j - \eta \frac{\partial \mathcal{L}}{\partial {\bf w}_j} = {\bf w}_j - \eta \sum_{n=1}^N (y_{nj} - t_{nj}){\bf x}_n$$ $$b_j \leftarrow b_j - \eta \frac{\partial \mathcal{L}}{\partial b_j} = b_j - \eta \sum_{n=1}^N (y_{nj} - t_{nj})$$

パラメータの扱いについて¶

1次元／多次元ロジスティック回帰双方においてモデルのパラメータである重みとバイアスを${\bf w}, b(or\ {\bf b})$のように別々に扱っていたが，これはデータ点ベクトルを$\tilde{\bf x} = (1, x_0, \dots, x_D)$のように最初の次元がつねに1となる拡張ベクトルとして扱い，重みを$\tilde{\bf w} = (w_0, w_1, \dots, w_K)$のようにしてバイアスを1つめの次元に追加して扱えば，$\tilde{\bf w}^{\rm T}\tilde{\bf x}$としてこれまでの${\bf w}^{\rm T}{\bf x} + b$を扱える．多次元の場合も同様で，$\tilde{\bf W}^{\rm T}\tilde{\bf x}$と扱えるようになる．するとこれまでの更新式の導出などにおけるバイアスの考慮は（数式上）必要がなくなる．自動的に重みの中の要素のひとつとして扱われるようになるからである．

実装¶

まずはデータ点とラベルのセットを用意する．

上で定義したソフトマックス関数は複数の多次元ベクトルに対してベクトルごとにsoftmaxを計算する仕様になっていないので，複数ベクトルを行列として渡した場合にも正しくベクトルごとにsoftmaxが適用されたものを返す関数として再定義する．

次に，$\softmax({\bf W}^{\rm T}{\bf X} + {\bf b})$を計算する関数p_y_given_x(x)を用意する．これは与えられたデータ点が列方向に並んだ$N \times D$次元の行列を与えると，行ごと（データ点ごと）に各クラスへの所属度確率ベクトルを計算し，これを列方向へ並べた$N \times K$行列を返す関数である．

grad関数は上で定義した勾配を求める関数である．定義では総和をとっていたが，ここでは平均をとることにする．

実際の学習を行う前に，学習率$\eta$とミニバッチサイズを設定する．学習ループではミニバッチごとにデータセット全体を学習することを100回繰り返している．結果として各クラスへ所属度確率の誤差の平均の推移をプロットしている．

さらに，識別境界を以下に図示する．