scipyで任意の目的関数を最小化する

この記事はMicroAd Advent Calendar 2018の9日目の記事です。

注意事項

私は数理最適化の専門家ではないため、最適化問題のアルゴリズムには全く詳しくありません。
ただある程度機械学習をやっていると、scikit-learnなどの機械学習ライブラリでは対処できない方法に出会うときがあり、そういう場面での何らかのヒントになれば良いと思って記事を書きます。
という予防線を張っておきます。

はじめに

機械学習（特に教師有り学習）は、基本的に任意の特徴空間における関数近似の問題です。
近似する関数のパラメータを求めるために、何らかの目的関数¹を最小化（最大化も考えられるけど一般的には最小化）することを考えます。
例えばロジスティック回帰では、以下のような負の対数尤度 (Negative log-likelihood) が目的関数です。

E(w) = - \sum_{i=1}^{N} y_i log(h_w (x_i)) + (1-y_i) log (1 - h_w(x_i))

ここで、$w$はパラメータベクトル、$N$はサンプルサイズ、$x_i$は特徴量ベクトル、 $y_i$ は教師信号、$h_w (x_i)$は以下のようなシグモイド関数です。

h_w (x_i) = \frac{1}{1 + e^{-w^T x_i}}

こういった一般的な方法は普通にscikit-learnなどに実装されているので、何の問題もありません。
しかし、論文などで提案されている方法を実装しようとする際、論文では独自のオレオレ目的関数が提案されていることがほとんどです。
こういった場合、自分で目的関数の最小化を実装しなければなりません。
目的関数の勾配を計算できれば、最急降下法や確率的勾配降下法なら普通に実装できるかもしれません。
ただ、このような数値計算には誤差等の問題がつきもので、出来るなら信頼のあるライブラリに任せてしまうほうが安全です。
ゆとり世代ど真ん中の私としては、自分で実装する選択肢はありませんでした。

scipyによる目的関数最小化

Pythonのscipy.optimizeモジュールに、最適化問題を解くアルゴリズムの実装があります。
順を追って使い方の説明をしていきます。

普通の関数の最小化

まずは一番簡単な例として、目的関数として以下の二次関数を考えます。

E(\theta) = (\theta - 2) ^ 2

目的関数の$\theta$に関する微分は以下のとおりです。

\frac{dE(\theta)}{d\theta} = 2(\theta - 2)

この最適解は$\theta=2$ですね。
では、scipyを用いて解いてみます。BFGS法 (scipy.optimize.fmin_bfgs()) というアルゴリズムで解いてみました。

import scipy.optimize

def objective_function(theta):
    '''目的関数''' 
    return (theta - 2) ** 2

def gradient(theta):
    '''勾配'''
    return 2 * (theta - 2)

theta_opt = scipy.optimize.fmin_bfgs(f=objective_function, x0=[0], fprime=gradient)
print(theta_opt)

scipy.optimize.fmin_bfgsの引数として、fに目的関数を、x0にパラメータの初期値を、fprimeに勾配の関数（導関数）を与えます。
結果は、最適解であるtheta_opt = 2.0となりました。完全に理論解と一致していますね。
ここからがすごい部分ですが、何と勾配を与えなくても解けてしまいます。fprimeは任意の引数なのです。
この場合、内部で自動微分が行われています。
この辺りは全く詳しくないので説明は出来ませんが、何と勾配を導出しなくても目的関数だけ与えれば計算が出来てしまうのです。
一応プログラムを記載しておきますが、以下のような感じです。

theta_opt = scipy.optimize.fmin_bfgs(f=objective_function, x0=[0])
print(theta_opt)

最適解はtheta_opt = 1.99999998となりました。微妙に誤差は生じてしまいますね。

ロジスティック回帰の目的関数の最小化

次に実践的な例として、ロジスティック回帰の目的関数の最小化を行います。
先程の例との違いは、パラメータの最適解を求めるために学習データが存在することです。
これらは、scipy.optimize.fmin_bfgsの引数としてargsに渡すことができます。
まず説明のため、scikit-learnのmake_classificationを使用して分類用の人工データを作成します。
ちなみにこれ以降、乱数のシードは公開日の1209に統一しています。機械学習でもこういった願掛けは大事ですよね。違うか。

from sklearn.datasets import make_classification
X, y = make_classification(n_samples=1000, n_features=10, n_informative=8, n_classes=2, random_state=1209)

後にscikit-learnでロジスティック回帰を学習した結果と比較したいので、学習データとテストデータに分割します。

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1209, stratify=y)

目的関数を以下のように定義しておきます。
scikit-learnのロジスティック回帰の実装と合わせるため、目的関数に$w$のL2ノルムを加えています。

import numpy as np

def sigmoid(x):
    '''シグモイド関数'''
    return 1 / (1 + np.exp(-x))

def negative_log_likelihood(w, X, y):
    '''負の対数尤度'''
    z = np.dot(X, w)
    return -np.sum((y * np.log(sigmoid(z)) + (1 - y) * np.log(1 - sigmoid(z))))

def objective_function(w, X, y):
    '''目的関数'''
    return negative_log_likelihood(w, X, y) + np.linalg.norm(w)

後は、scipy.optimize.fmin_bfgsに渡すだけです。
ただし、切片を考慮して工夫を加えてあります。

# 切片を含めて行列演算を行うため全サンプルに1つ要素を追加（値は1）
intercept = np.ones(X_train.shape[0]).reshape(X_train.shape[0], 1)
X_train = np.concatenate((X_train, intercept), axis = 1)
intercept = np.ones(X_test.shape[0]).reshape(X_test.shape[0], 1)
X_test = np.concatenate((X_test, intercept), axis = 1)

# パラメータの初期値を設定（平均0、標準偏差1の正規乱数）
np.random.seed(1209)
w = np.random.randn(X_train.shape[1])

w_opt = scipy.optimize.fmin_bfgs(f=objective_function, x0=w, args=(X_train, y_train))

print(f'係数:{w_opt[:-1]}\n切片:{w_opt[-1]}')

またテストデータに対して予測を行い、confusion matrixを計算しておきます。

from sklearn.metrics import confusion_matrix

y_proba_scipy = sigmoid(np.dot(X_test, w_opt))
y_pred_scipy = (y_proba_scipy >= 0.5) * 1
print(confusion_matrix(y_test, y_pred_scipy))

比較のため、scikit-learnのロジスティックでも同じことをやってみます。

from sklearn.linear_model import LogisticRegression

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1209, stratify=y)
model = LogisticRegression(random_state=1209)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
print(f'係数:{model.coef_}\n切片:{model.intercept_}')
print(confusion_matrix(y_test, y_pred))

結果を見てみましょう。
confusion_matrixに関しては完全に一致しています。よって分類性能は全く同じです。
学習後のパラメータについては、多少誤差はあるものの大体は一致しています。

# scipy
係数:[ 0.1895906   1.19134589 -0.76077283 -0.49307632  0.45103504 -0.61238509
  0.15386609 -0.57944113  0.36421485 -0.16949388]
切片:0.9138857938228944
[[138  13]
 [ 21 128]]

# scikit-learn
係数:[[ 0.18475573  1.17797849 -0.74865278 -0.48591727  0.4412075  -0.60217896
   0.15192604 -0.57577281  0.36166131 -0.16599469]]
切片:[0.88543542]
[[138  13]
 [ 21 128]]

任意の目的関数の最小化

最後に発展的な例として、オレオレ目的関数を最小化します。
以下のように、異なるデータについて同時に目的関数を最小化したい場合を考えます。

E(w_1, w_2) = - \sum_{i=1}^{N} y_{1i} log(h_{w_1} (x_{1i})) + (1-y_{1i}) log (1 - h_{w_1}(x_{1i})) - \sum_{i=1}^{N} y_{2i} log(h_{w_2} (x_{2i})) + (1-y_{2i}) log (1 - h_{w_2}(x_{2i})) + \frac{1}{2} (\|w_1\|^2_2 + \|w_2\|^2_2)

これはもうscikit-learnでは無理ですね。
これまでと同様に、scipy.optimize.fmin_bfgsを用いて最適解を求めてみましょう。
まず、目的関数を定義します。

import numpy as np

def sigmoid(x):
    '''シグモイド関数'''
    return 1 / (1 + np.exp(-x))

def negative_log_likelihood(w, X, y):
    '''負の対数尤度'''
    z = np.dot(X, w)
    return -np.sum((y * np.log(sigmoid(z)) + (1 - y) * np.log(1 - sigmoid(z))))

def objective_function(w, X1, y1, X2, y2):
    '''目的関数'''
    dim1 = X1.shape[1]

    w1 = w[:dim1]
    w2 = w[dim1:]

    return  negative_log_likelihood(w1, X1, y1) + negative_log_likelihood(w2, X2, y2) + (np.linalg.norm(w1) + np.linalg.norm(w2)) / 2

後は同様に、分類用の人工データを作成して最適解を解いてみます。

from sklearn.datasets import make_classification
import numpy as np

# X1, y1の作成
X1, y1 = make_classification(n_samples=1000, n_features=10, n_informative=8, n_classes=2, random_state=1209)
intercept = np.ones(X1.shape[0]).reshape(X1.shape[0], 1)
X1 = np.concatenate((X1, intercept), axis = 1)

# X2, y2の作成
X2, y2 = make_classification(n_samples=1000, n_features=8, n_informative=5, n_classes=2, random_state=1209)
intercept = np.ones(X2.shape[0]).reshape(X2.shape[0], 1)
X2 = np.concatenate((X2, intercept), axis = 1)

# パラメータの初期化
w1 = np.random.randn(X1.shape[1])
w2 = np.random.randn(X2.shape[1])
w = np.concatenate((w1, w2))

w_opt = scipy.optimize.fmin_bfgs(f=objective_function, x0=w, args=(X1, y1, X2, y2))
print(f'w1={w_opt[:w1.shape[0]]}\nw2={w_opt[w1.shape[0]:]}')

求められた最適解は以下のようになりました。

w1=[ 0.24380479  1.22542828 -0.77760999 -0.4975671   0.47881832 -0.55993407
  0.16043769 -0.54981939  0.36073229 -0.19703398  1.08620108]
w2=[-0.09020432  0.21506923  0.19820872  0.00251996  1.08074929 -0.0553308
 -0.54615085 -0.21130013 -0.32719551]

まとめ

scipyを用いて任意の目的関数を最小化する方法を紹介しました。
勾配を求めなくても計算出来るのが非常に良いですね。
TensorFlowやKerasでも同じことが出来ると思うので、試してみようと思います。

1. 画像処理の世界ではエネルギー関数とも呼ばれますね。 ↩