post Image
ラグランジュの未定乗数法(言語処理のための機械学習入門)


はじめに

始めまして、たなしょうです。皆さんは以下のようなワードを聞いたとき、何も見ずにその手法を解説することができますか?

  • ラグランジュの未定乗数法
  • EM アルゴリズム
  • 隠れマルコフモデル (HMM)
  • サポートベクトルマシン (SVM)

僕の場合何となくで理解しているものが多く、正直何も見ずに解説することはできません。色々な教科書、授業で触れているにも関わらず、です。これを克服ためには、アルゴリズムを実装するのが一番です。そこで今回は文章力及び実装力向上も兼ねて、コロナ社から出版されている自然言語処理シリーズの「言語処理のための機械学習入門」に登場する上記のようなアルゴリズムを実装し、コードの解説を Qiita に投稿していこうと思います。記事中でp.xx 例題○○ などと書いてある場合は特筆しない限り前述のページ数、例題を意味します。

実装は Jupyter Notebook 上で Python を用いて行い、実装したコードはすべて僕の GitHub レポジトリに公開します。コードレビューや誤りの指摘、疑問質問など、積極的にコメントしてくださると幸いです。


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

さて、第一回となる今回は皆さんおなじみ、ラグランジュの未定乗数法を実装したいと思います。この記事の目的はあくまで実装なので詳しい解説は上記の本に譲りますが、簡単に言うとある凸計画問題に制約が加わった場合、これを解くための手法です。

今回は p.17 例題1.7 に載っている以下の最大化問題を解いてみます。

\begin{align}

\max . & - x^2 - y^2 \\
{\rm s.t.} & \ x + y - 1 = 0
\end{align}


関数のプロット

まずはこの目的関数と制約式をプロットしてみたいと思います。numpy と matplotlib を使い、コードを実装すると以下のようになります。

# coding: utf-8

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

def plot_plane(vrange, inc):
"""Plot the mesh plane and the constraint equation.

Args:
vrange: A (min, max) tuple.
inc: The incremental value.

"""
xmin, xmax = vrange
ymin, ymax = vrange
x = np.arange(xmin+inc, xmax+inc, inc)
y = np.arange(ymin+inc, ymax+inc, inc)
xs, ys = np.meshgrid(x, y)
z = -(xs ** 2) -(ys ** 2)
const_eq = -x + 1

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.imshow(z, cmap=plt.cm.gray)
ax.autoscale(False)
plt.plot(x * x.size, const_eq * y.size + y.size / 2)
plt.colorbar()
plt.title("Image plot of $- x^2 - y^2$ and $x + y -1 = 0$")
plt.gca().invert_yaxis()
plt.show()

plot_plane((-1, 1), 0.1)

目的関数をメッシュで描画した後、その上に制約式をプロットしています。ここで躓いたのが、目的関数をプロットした後に制約式をプロットしようとすると、直線の位置がずれたりしてなかなか狙い通りにプロットされないんですよね。試行錯誤してなんとかそれっぽくプロットされましたが。。。

上記のコードを実行すると次のようなグラフが描画されます。$x$ 軸、$y$ 軸がそれぞれ配列の添え字になっているので少し分かりにくいですが、$-1 < x, y < 1$ の範囲で目的関数の値及び制約式(青線)が描画されているのが分かります。グラフを見れば(もしくは目的関数を $x, y$ でそれぞれ偏微分すれば)わかる通り、目的関数は $x, y = 0$ で最大値 $0$ をとりますが、制約式のために $x, y$ は同時に $0$ になることはできません。ここで登場するのが、ラグランジュの未定乗数法です。

plot.png


ラグランジュ関数

目的関数を $f(\boldsymbol{x})$, 制約式を $g(\boldsymbol{x})$ とするとラグランジュ関数は次のように表せられます。ここで $\lambda$ は未定乗数を意味します。

L(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) + \lambda g(\boldsymbol{x})

ラグランジュ関数を実装すると以下のようになります。

def func(X):

"""Calculate a value of the method of Lagrange multipliers using given X.

Args:
X: It consists of x, y, and L(it means lambda, the Lagurange multiplier.)
Returns:
f: The calculated value.

"""


x, y, L = X
f = -(x ** 2) - (y ** 2) - L * (x + y - 1)
return f


ラグランジュ関数の偏微分

上記のラグランジュ関数を $\boldsymbol{x}$ に関して偏微分したものが $0$ になり、かつ与えられた制約を満たすとき、最適な解が得られます。すなわちラグランジュ関数を $x, y, \lambda$ でそれぞれ偏微分したものが全て $0$ になる $x, y, \lambda$ の組み合わせが制約のもとで目的関数を最大化します。

実際の偏微分は以下のプログラムのように近似したものとなります。

def dfunc(X):

"""Calculate partial derivatives using given X.

Args:
X: It consists of x, y, and L(it means lambda.)
Returns:
dLambda: The calculated values of partial derivatives.

"""


dLambda = np.zeros(len(X))
h = 1e-3
for i in range(len(X)):
dX = np.zeros(len(X))
dX[i] = h
dLambda[i] = (func(X + dX) - func(X - dX)) / (2 * h)
return dLambda

さて、あとは上記の $\boldsymbol{dLambda} = \boldsymbol{0}$ となる $x, y, \lambda$ の組み合わせを求めるだけです。これは x in range(...) のようにループを回して求めることもできるのですが、幸い scipy にはこういった問題を解くための fsolve という素晴らしい関数があります。これは第一引数に $f(\boldsymbol{x}) = 0$ としたい関数 $f$, 第二引数に関数 $f$ の引数 $\boldsymbol{x}$ の初期値を与えてやることで、$f(\boldsymbol{x}) = 0$ となる $\boldsymbol{x}$ を返してくれます。すなわち、以下のようなプログラムで、目的関数を最大化する $x, y$ 及び最大値を求めることができます。

from scipy.optimize import fsolve

max_arg = fsolve(dfunc, [-1, -1, 0])
print("max arguments: x = {0:.2}, y = {1:.2}".format(max_arg[0], max_arg[1]))
print("max value: {0:.2}".format(func(max_arg)))

これを実行すると以下のような出力が得られます。すなわち、$x, y = 0.5$ のとき目的関数は制約 $x + y – 1 = 0$ のもと最大値 $-0.5$ をとることが分かりました。グラフを見てもちょうどそのあたり((15, 15) のあたり)で目的関数が最大値をとっていることが分かります。

max arguments: x = 0.5, y = 0.5, lambda = 0.5

max value: -0.5


まとめ

今回は初めての投稿ということで簡単にラグランジュの未定乗数法を実装してみましたが、いざ実装してみるとなかなか大変でした(特に matplot が。。。)。ラグランジュの未定乗数法は EM アルゴリズムやサポートベクトルマシンなど機械学習において重要なアルゴリズム中でも使われている手法であり、ウォーミングアップにちょうどよいセレクトだったかと思います。

プログラムの実装には Python によるデータ分析入門およびこのページを参考にしました(パクったとも言う)。

今後は EM アルゴリズムやサポートベクトルマシン、HMM などを実装していきたいと思います。


宣伝

僕のプロフィール、作品は Wantedly に公開しています。


『 Python 』Article List