MLエンジニアへの道 #4 - ソフトマックス回帰

この記事はソフトマックス回帰(多項ロジスティック回帰)を紹介します。

分類モデルの評価方法に入る前に、ソフトマックス回帰について紹介します。ロジスティック回帰は二項分類のみを扱うものであるのに対してソフトマックス回帰は多クラス分類を行うことができます。

ソフトマックス関数

複数のクラスを分類するための直感的な方法の一つは、すべてのクラスに対して複数のロジスティック関数を適合させ、最も高い確率を持つクラスを選ぶことです。しかし、これは多クラス分類で望むことではありません。なぜなら、私たちはすべてのクラスの分布を得たいのであり、各クラスの異なる分布のargmaxを求めたいわけではないからです。それならば、なぜ次のようにクラスごとの確率を総確率で割らないのか、疑問に思うかもしれません。

p_k = \frac{p_k(x)}{\sum p_i(x)} \\ = \frac{\frac{1}{1+e^{-\beta_k x}}}{\sum \frac{1}{1+e^{-\beta_i x}}}

ここで、 $p_k$ は $x$ が与えられたときにカテゴリが $k$ である確率を示し、それは $k$ のロジスティック関数によって計算された確率をすべてのカテゴリのロジスティック関数の総確率で割ることで求められます。しかし、見ての通り、この式は……厄介です。

では、ロジットはどうでしょう？ロジスティック回帰のためにロジットと説明変数の線形回帰を取りました。ここでもロジットを使うことで簡素化できるでしょうか？見てみましょう。

p_k = \frac{\beta_k x}{\sum \beta_i x}

一見するとこれで良いように思えますが、これは良くありません。なぜなら、ロジットは $-\infty$ から $\infty$ の範囲にあり、負になることもあるからです。確率は0から1の範囲でなければならないため、残念ながらロジットは使用できません。確率とロジットの代わりに、オッズを使用することができます。ロジットの指数を取ると、オッズにたどり着きます。( $e^{\log(odds)} = odds$ )

p_k = \frac{e^{\beta_k x}}{\sum e^{\beta_i x}}

指数化は $-\infty$ から $\infty$ の範囲の数を $0$ から $\infty$ の範囲にマッピングできるため、 $p$ が0から1の範囲になることを保証できます。さらに、これは複数のロジスティック関数を使用するよりもはるかに簡単です。この関数は ソフトマックス関数 と呼ばれ、ロジスティック回帰の一般化されたバージョンであり、多クラス分類にとって非常に便利な関数です。

ソフトマックス関数は冗長なパラメータを持つ

ソフトマックス関数のユニークな特徴の1つは、過剰パラメータ化されていること、つまり冗長なパラメータを持っていることです。その理由を見てみましょう。パラメータベクトル $\psi$ をすべてのクラスのパラメータベクトル $\beta$ から引くと仮定します。すると、ソフトマックス関数は次のようになります。

p_k = \frac{e^{(\beta_k-\psi) x}}{\sum e^{(\beta_i-\psi) x}}

上記の式を整理しましょう。

p_k = \frac{e^{\beta_k x}e^{-\psi x}}{\sum e^{\beta_i x}e^{-\psi x}} \\ = \frac{e^{\beta_k x}}{\sum e^{\beta_i x}}

$e^{-\psi \cdot x}$ はキャンセルできるため、元のソフトマックス関数に戻ります。これがソフトマックス関数について何を意味しているのでしょうか？それは、同じ予測を生み出すパラメータの組み合わせが複数存在することを意味します。これは直感的にも理解することができます。なぜなら、すべてのクラスに対して、10倍のオッズを導出するパラメータは、正規化によって元のパラメータと全く同じ確率予測を行うからです。

しかし、より重要なのは、 $\psi$ を $\beta_k$ の逆ベクトル（ $\beta_k - \psi = \overrightharpoon{0}$ となるような）に設定でき、 $\beta_k$ を実質的に排除して同じ結果にたどり着くことができるということです。これも直感的に理解することができます。なぜなら、一つのクラス以外のすべてのクラスの確率分布がわかっている場合、その欠けているクラスの確率を他のすべてのクラスの合計を1から引くことで推測できるからです。

ロジスティック関数との関係性

クラスの数が2の場合、ソフトマックス関数がロジスティック関数と等価であることを、ソフトマックス関数が過剰パラメータ化されているという事実を使って示すことができます。クラスの数が2（クラスが1と2として表される）場合、クラス1の確率をソフトマックス関数で次のように表すことができます：

p_1(x) = \frac{e^{\beta_1 x}}{e^{\beta_1 x}+e^{\beta_2 x}}

ここで、 $\psi = \beta_2$ と設定し、ソフトマックス関数が過剰パラメータ化されているという事実を用いて、パラメータベクトルからそれを引きます。

p_1(x) = \frac{e^{(\beta_1 - \beta_2) x}}{e^{(\beta_1 - \beta_2) x}+e^{0 x}} \\ = \frac{e^{(\beta_1 - \beta_2) x}}{e^{(\beta_1 - \beta_2) x}+1}

$\beta_1 - \beta_2$ を $\beta$ で置き換えると、次のようになります：

p_1(x) = = \frac{e^{\beta x}}{e^{\beta x}+1} \\ = \frac{1}{1+e^{-\beta x}}

これはロジスティック関数と同じ式です！したがって、ソフトマックス関数はロジスティック関数の一般化されたバージョンであると言えます。

ステップ 1. データ探索

ここでは、ロジスティック回帰で使用したのと同じデータセット、Irisデータセットを使用します。この部分は前回の記事で既に行ったので省略します。Setosa対その他を分類する代わりに、今回はSetosa、Versicolor、Virginicaの全ての種を分類することを目指します。

ステップ 2. データ前処理

種を0、1、2のような数値で表す代わりに、それぞれの種に対応するインデックスを持つベクトルを使用できます。このプロセスをOne-Hotエンコーディングとよび、keras.utilsのto_categoricalを使って実装できます。

from keras.utils import to_categorical 
iris_df.drop('species', axis= 1, inplace= True)
target_df = pd.DataFrame(columns= ['species'], data= iris.target)
 
X = np.array(iris_df)
y = to_categorical(np.array(target_df['species']))

説明変数にはX、分類したい結果変数にはyを使用するのが一般的です。さらに、データセットを訓練データとテストデータに分割する必要があります。

from sklearn.model_selection import train_test_split
 
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size= 0.33, random_state= 101)

上記のコードでは、データセットをランダムに分割して、訓練データとテストデータがそれぞれ元のデータセットの67%と33%になるようにしています。

ステップ 3. モデル

モデルの定義

この記事の最初の部分で述べたように、ソフトマックス関数を使用して多クラス分類を行うことができます。関数は次の通りです:

p_k = \frac{e^{\beta_k x}}{\sum e^{\beta_i x}}

上記をベクトル形式に書き直して、この分類タスクの全体の確率分布を表すことができます:

\overrightharpoon{p} = \begin{bmatrix} p_0(x) \\ p_1(x) \\ p_2(x) \end{bmatrix} = \frac{1}{\sum_{i=0}^{2} e^{\beta_i x}}\begin{bmatrix} e^{\beta_0 x} \\ e^{\beta_1 x} \\ e^{\beta_2 x} \end{bmatrix}

コスト関数

もう覚えていると思いますが、最尤推定（MLE）を使用して、データが生成されたモデルを最もよく近似するパラメータを見つけることができます。この状況における尤度関数を定義してみましょう。

ロジスティック回帰で使用したのと同じ論理がここでも適用されます。モデルが正しく分類するデータが多いほど、データに対するモデルの尤度は高くなるべきです。これは、正しいラベルに関連付けられる確率が高いほど、モデルの尤度が高くなることを意味します。このことから、尤度関数は次のようになります:

L(\beta) = \prod_{y_i=0}(p_0(x_i)) * \prod_{y_i=1}(p_1(x_i)) * \prod_{y_i=2}(p_2(x_i))

コスト関数のために負の対数尤度をとると、次のようになります:

J(\beta) = -\sum_{i=1}^{n}\sum_{k=0}^{2}1\{y_i = k\}log(p_k(x_i))

なじみがあるように見えますね？そうです、これは多クラス状況でのクロスエントロピーと全く同じです！ロジスティック回帰と同様に、負の対数尤度は、真の確率分布と予測された確率分布の間のクロスエントロピーと一致します。

学習メカニズム

ソフトマックス関数は、同じ予測分布に対して複数のパラメータを持つため、多くの最尤推定の方法を使用できません。しかし、勾配降下法を使用できます。そして勾配降下法は、コスト関数をパラメータに関して偏微分する必要があります。ではまず、パラメータに関してコスト関数を表現してみましょう:

J(\beta) = -\sum_{i=1}^{n}\sum_{k=0}^{2}1\{y_i = k\}log(\frac{e^{\beta_k x_i}}{\sum_{j=0}^{2}e^{\beta_j x_i}}) \\ = -\sum_{i=1}^{n}\sum_{k=0}^{2}1\{y_i = k\}(log(e^{\beta_k x_i})-log(\sum_{j=0}^{2}e^{\beta_j x_i})) \\ = -\sum_{i=1}^{n}\sum_{k=0}^{2}1\{y_i = k\}(\beta_k x_i-log(\sum_{j=0}^{2}e^{\beta_j x_i}))

ここで、 $\beta x_i$ はロジットの線形方程式の略です。このコスト関数を、クラス $0$ の傾きと切片に関してチェーンルールを使用して偏微分してみましょう。傾きに対しては：

\frac{\partial}{\partial\beta_0} = -\sum_{i=1}^{n}1\{y_i = 0\}(x_i-\frac{1}{\sum_{j=0}^{2}e^{\beta_j x_i}}e^{\beta_0 x_i}x_i) \\ = -\sum_{i=1}^{n}1\{y_i = 0\}(x_i-\frac{e^{\beta_0 x_i}}{\sum_{j=0}^{2}e^{\beta_j x_i}}x_i) \\ = -\sum_{i=1}^{n}1\{y_i = 0\}(1 - p_0(x_i))x_i

切片に対しては：

\frac{\partial}{\partial\beta_0} = -\sum_{i=1}^{n}1\{y_i = 0\}(1-\frac{1}{\sum_{j=0}^{2}e^{\beta_j x_i}}e^{\beta_0 x_i}) \\ = -\sum_{i=1}^{n}1\{y_i = 0\}(1-\frac{e^{\beta_0 x_i}}{\sum_{j=0}^{2}e^{\beta_j x_i}}) \\ = -\sum_{i=1}^{n}1\{y_i = 0\}(1-p_0(x_i))

上記はすべてのクラスに適用されます。これらの導関数は、クラス $k$ のパラメータはクラス $k$ のデータによってのみ更新されることを意味します。なぜならクラスが $k$ 以外のデータに対しては導関数は0になるからです。

さらに、勾配降下法では、勾配の逆を取り、パラメータに加えます。つまり、この場合予測確率と1との差をデータに掛けたものの合計を対応するパラメータに加えるだけで済むということがわかります。これにより、コードの実装は非常に簡単になります。

コードの実装

上記の情報を使用して、SoftmaxRegressionGDモデルを作成しましょう。

class SoftmaxRegressionGD():
  def __init__(self, lr=0.001):
    self.W = np.zeros((X.shape[1], y.shape[1]))
    self.b = np.zeros(y.shape[1])
    self.lr = lr # Learning rate
    self.history = [] # History of loss
 
  def predict(self, X):
    odds = np.exp(np.matmul(X, self.W.T)+self.b)
    total_odds = np.sum(odds, axis=1)
    return odds / total_odds[:, np.newaxis]
 
  def fit(self, X, y, epochs=100):
    for i in range(epochs):
      pred = self.predict(X)
 
      self.history.append(log_loss(y, pred))
 
      diff = 1 - pred
      grad_W = np.matmul(X.T, diff*y)
      grad_b = np.sum(diff*y, axis=0)
 
      self.W += self.lr * grad_W
      self.b += self.lr * grad_b
    return self.history

操作のベクトル化はロジスティック回帰と比較して直感的でないかもしれませんので、時間をかけて各ステップで何が行われているのかをしっかり理解することを強くお勧めします。

SoftmaxRegressionGDを初期化し、以下のように訓練データに適合させることができます。

sm = SoftmaxRegressionGD()
history = sm.fit(X_train, y_train, epochs=500)

各エポックにおける負の対数尤度損失を追跡することで、勾配降下法がどのように機能しているかをプロットできます。

plt.plot(history)
plt.title("Loss vs Epoch")
plt.ylabel("Loss")
plt.xlabel("Epoch")
plt.show()

上記のプロットからわかるように、損失は各エポックで徐々に減少しています。これは、勾配降下法が機能している良い兆候です。

そろそろ記事が長くなってきたので、次の記事で分類モデルの評価方法に関して話したいと思います。

リソース

UFLDL Tutorial. n d. Softmax Regression Stanford University.