はじめに
シリーズ第3弾で今回は制約付き非線形計画問題を解く主双対内点法を実装した*1。前回の準ニュートン法の記事はこちら。
yamagensakam.hatenablog.com
制約付き非線形計画問題を式で表すと、以下のような最適化問題になる。
なお今回取り上げる主双対内点法は、すぐ下に記すKKT条件が最適性の必要条件となることが保証されなければ成立しない。これが保証されるために制約条件が満たすべき条件が存在する。そのような制約条件に関する条件を制約想定と言い、1次独立条件やSlater条件などの種類がある。制約想定に関してはこちらが参考になる*2。
qiita.com
手法
内点法は最適性の1次の必要条件であるKKT条件を満たす点を実行可能領域の内部を通りながら探す手法である。まずKKT条件について簡単に触れた後、主双対内点法の考え方、アルゴリズムを述べる。
式(1)の最適化問題に対して、以下のような関数を考える。
これはラグランジュ関数と呼ばれる関数で、式(1)の最適化問題が上述の制約想定を満たす場合、この問題における1次の最適性の必要条件はラグランジュ関数を使って、
と表すことができる。これはKKT条件と呼ばれ、主双対内点法ではこのKKT条件を満たす(に十分に近い点)を探索することになる。
スラック変数とバリア関数の導入
ここで式(1)にスラック変数とバリア関数というものを導入した以下のような最適化問題を考える。
がスラック変数でこの変数の導入により不等式制約は等式制約に置き換わる。また、がバリア関数で実行可能領域の内部では有効な値を取り、であれば制約の境界であるに近づくにつれ無限に発散していく。従って、バリア関数は不等式制約を満たさないことに対するペナルティとして働く。で式(4)は元の最適化問題式(1)と等価になる。そのため、主双対内点法では式(4)の停留点方向に移動しながら徐々にをに近づけていくという操作を繰り返すことで式(1)の停留点を見つける。
主双対内点法では現時点がだとして、と更新する。以下、更新方向、ステップ幅、をどのように更新するか述べる。
更新方向の計算
式(4)の停留点ではKKT条件より、
が満たされる。そのため更新方向は式(5)を近似的に満たすように決める。つまり、
を満たすように(を決めたいのだが、式(6)の方程式を解くのは困難なので1次近似して解く。1次近似なので、
を式(6)の等式に代入し、2次項を無視すると、以下の1次連立方程式が得られる。
ただし、
であり、はアダマール積を表す。この式(8)を解くことで更新方向が得られる。
ステップ幅の決め方
次にステップ幅だが、まずとなる必要があるため、これを満たすを見つける*3。更にその条件を満たす中でもなるべく適切なを選びたいので、以下の目的関数と制約のペナルティを組み合わせたメリット関数を定義し、
を最小にするを直線探索により見つける*4。はペナルティの重みを表すパラメータで、今回の実装では固定値をあらかじめ与えている。
の更新
最後にの更新だが、これはが式(6)を満たすのならなので*5、となっているのが理想的である。このことから、適切なパラメータをあらかじめ与えて、
として更新する。
主双対内点法のアルゴリズムは以下の通り。
- 初期点とバリア関数のパラメータの初期値を定める。
- が十分に小さければ終了。
- 式(8)を解いて探索方向を決める。
- ステップ幅を決める。
- として更新する。
- 式(11)でを求める。
- 2に戻る。
その他補足
実装する上で困ったこと、まだよく分かってないこと、その他細かい話など。
初期点について
内点法は基本的に実行可能領域の内部を通って解を見つけるアルゴリズムなので初期点も内部に置かないといけないと思ったのだが、実は今回述べたアルゴリズムだと初期点を実行可能領域外の点にしても試した範囲では実行可能領域の停留点に収束した。
今回の方法が実行可能領域外を初期点にしても収束するという理論的保証があるかはよく分かっていないが、雰囲気的にはが十分に大きければ式(4)の最適解は実行可能領域内の境界から遠く離れた点であり、式(8)を満たす方向に点を動かすということはその式(4)の最適解に近づけることを意味するので、たとえ初期点が実行可能領域外だったとしても内部に入る方向に動き最終的には内部の停留点に収束するという感じだろう*6。ということで実装では、適当な値を初期値としている。
ラグランジュ関数のヘッセ行列(式(9)参照)が正定値でない場合、式(8)を満たすが下降方向とは限らない。一方、ステップ幅を探索するときは、なるべくメリット関数が小さくなるを探すので、が上昇方向であれば限りなくに近いが得られなかなか収束しないということが起こりうる。そのため、今回の実装では上述の通りが小さくなりすぎないよう探索の範囲を絞っている。また、そのほかの対策方法として、ヘッセ行列の代わりに準ニュートン法同様にヘッセ行列を近似する正値対称行列を用いる場合もある模様。
実装
以下が実装した主双対内点法のpythonコード。
import numpy as np
import copy
class problem():
def __init__(self, objective,
inequality_constraints_list,
equality_constraints_list):
for constraint in inequality_constraints_list:
if objective.dim != constraint.dim:
raise Exception('dimension mismatch')
for constraint in equality_constraints_list:
if objective.dim != constraint.dim:
raise Exception('dimension mismatch')
self.dim = objective.dim
self.objective = objective
self.inequality_constraints_list = inequality_constraints_list
self.equality_constraints_list = equality_constraints_list
def merit_function(problem, rho, eta, x, s):
res = problem.objective.func(x)
res -= rho * np.sum(np.log(s))
g_sum = 0
for i, constraint in enumerate(problem.inequality_constraints_list):
g_sum += np.abs(constraint.func(x) + s[i])
res += eta * g_sum
h_sum = 0
for i, constraint in enumerate(problem.equality_constraints_list):
h_sum += np.abs(constraint.func(x))
res += eta * h_sum
return res
def backtracking_line_search(problem, x, s, u, dx, ds, du, rho, eta, beta):
phi = lambda x, s:merit_function(problem, rho, eta, x, s)
alpha1 = 1.0
alpha2 = 1.0
if(np.min(ds) < 0):
alpha1 = np.min(-s[ds < 0]/ds[ds < 0])
if(np.min(du) < 0):
alpha2 = np.min(-u[du < 0]/du[du < 0])
alpha = min(1.0, alpha1 * beta, alpha2 * beta)
min_merit = np.inf
res_alpha = alpha
init_alpha = alpha
while(alpha > init_alpha * 0.1):
new_x = x + alpha * dx
new_s = s + alpha * ds
m = phi(new_x, new_s)
if (m < min_merit):
min_merit = m
res_alpha = alpha
alpha *= beta
return res_alpha
def interior_point(problem, eps = 1e-9, eta =0.1, beta = 0.9, t = 0.5):
l = len(problem.inequality_constraints_list)
m = len(problem.equality_constraints_list)
n = problem.dim
rho = 1
x = np.ones(n)
s = np.ones(l)
u = rho / s
v = np.zeros(m)
while(1):
if rho < eps:
break
H = problem.objective.hessian(x)
d = problem.objective.grad(x)
Jg = np.zeros((l, n))
g = np.zeros(l)
for i, constraint in enumerate(problem.inequality_constraints_list):
H += constraint.hessian(x) * u[i]
Jg[i, :] = constraint.grad(x)
d += Jg[i, :] * u[i]
g[i] = constraint.func(x)
Jh = np.zeros((m, n))
h = np.zeros(m)
for i, constraint in enumerate(problem.equality_constraints_list):
H += constraint.hessian(x) * v[i]
Jh[i, :] = constraint.grad(x)
d += Jh[i, :] * v[i]
h[i] = constraint.func(x)
Du = np.diag(u)
Ds = np.diag(s)
P = np.zeros((n + 2 * l + m, n + 2 * l + m))
P[:n, :n] = H
P[:n, n + l: n + 2 * l] = Jg.T
P[:n, n + 2 * l:] = Jh.T
P[n:n + l, n:n + l] = Du
P[n:n + l, n + l:n + 2 * l] = Ds
P[n + l:n + 2 * l, :n] = Jg
P[n + l:n + 2 * l, n:n + l] = np.eye(l)
P[n + 2 * l:, :n] = Jh
r = np.concatenate([-d, rho - s * u, -(g + s), -h])
delta = np.linalg.solve(P, r)
dx = delta[:n]
ds = delta[n:n + l]
du = delta[n + l:n + 2 * l]
dv = delta[n + 2 * l:]
alpha = backtracking_line_search(problem, x, s, u,
dx, ds, du, rho, eta, beta)
x += alpha * dx
s += alpha * ds
u += alpha * du
v += alpha * dv
rho = t * u @ s / l
return x, s, u, v
式(1)の最適化問題はproblem
というクラスに記述し、そのインスタンスを関数interior_point
に渡すと停留点を探して出力する。problem
のobjective
に目的関数を表すクラスのインスタンスをセットし、inequality_constraints_list
には不等式制約の関数(式(1)の)を表すクラスのインスタンスのリスト、inequality_constraints_list
には等式制約の関数(式(1)の)を表すクラスのインスタンスのリストをそれぞれセットする。各クラスはの値を返すメソッドfunc
、を返すメソッドgrad
、を返すメソッドhessian
を実装してる必要がある。
なお、式(9)の連立方程式はnp.linalg.solve
で解いているが、零要素が多いためもう少し効率的に解く方法があるかもしれない。また、この実装だと不等式制約が存在しない問題は解けない(等式制約が存在しない問題は解ける)。
SVMは分類問題を解く教師あり学習の一種であり、その「学習」は最適化問題を解くことに対応する。通常はSMOなどの専用の効率が良い最適化アルゴリズムなどで解くが、今回はこれを主双対内点法で解いてみる。
SVMの学習で解くことになる最適化問題を簡単に説明する。SVMの学習では教師データ集合が得られてるとして、それをうまいこと分離できる関数とバイアス項を求める。具体的にはヒンジ損失関数に正則化項を加えた以下の最小化問題を解く。
この問題は関数のに関する最適化問題であり、加えて第2項(ヒンジ損失関数)は微分できないので、このまま解くのは難しい。そこでリプレゼンター定理などを使ってゴニョゴニョしたりヒンジ損失関数の制約条件への変換を行ったりすると、以下の凸2次計画問題に変形できる。
ここではカーネル関数と呼ばれ適切な正定値対称関数を選ぶ(例えば、など)。または行列目にを持つグラム行列である。
これらの理屈などについては、以下の資料や参考文献として挙げた『カーネル多変量解析』という書籍を参考にするとよい*7。
https://www.ism.ac.jp/~fukumizu/ISM_lecture_2004/Lecture2004_kernel_method.pdf
この最適化問題式(13)の解とを用いて、新たに与えられた入力は
としてクラスを分類できる。
以下が実装したSVMとそれを用いて検証用データの学習と分類を行うpythonコード。上述の主双対内点法の関数を利用している。
import numpy as np
import interior_point_method as ipm
import matplotlib.pyplot as plt
class objective_function:
def __init__(self, K, C, var_dim, constraint_num):
self.K = K
self.C = C
self.a_dim = var_dim
self.z_dim = constraint_num
self.dim = var_dim + constraint_num + 1
def func(self, a_z_b):
return 0.5 * a_z_b[:self.a_dim] @ self.K @ a_z_b[:self.a_dim] + \
self.C * np.sum(a_z_b[self.a_dim:self.a_dim + self.z_dim])
def grad(self, a_z_b):
return np.concatenate([self.K @ a_z_b[:self.a_dim],
self.C * np.ones(self.z_dim),
[0]])
def hessian(self, a_z_b):
H = np.zeros((self.dim, self.dim))
H[:self.a_dim, :self.a_dim] = self.K
return H
class train_data_constraint_creator:
def __init__(self, K, y):
self.K = K
self.y = y
class train_data_constraint:
def __init__(self, K, y, idx):
self.idx = idx
self.K = K
self.y = y
self.a_dim = K.shape[0]
self.z_dim = y.shape[0]
self.dim = self.a_dim + self.z_dim + 1
def func(self, a_z_b):
return -a_z_b[self.a_dim + self.idx] + 1 - \
self.y[self.idx] * \
(a_z_b[:self.a_dim] @ self.K[self.idx, :] + a_z_b[-1])
def grad(self, a_z_b):
return np.concatenate([-self.K[self.idx, :] * self.y[self.idx],
np.ones(self.z_dim) * -1,
[-self.y[self.idx]]])
def hessian(self, a_z_b):
return np.zeros((self.dim, self.dim))
def create(self, idx):
return self.train_data_constraint(self.K, self.y, idx)
class nonnegative_constraint_creator:
def __init__(self, K, y):
self.constraint_num = y.shape[0]
self.var_dim = K.shape[0]
self.constraint_num = y.shape[0]
self.dim = self.var_dim + self.constraint_num + 1
class function_nonnegative_constraint:
def __init__(self, var_dim, constraint_num, idx):
self.idx = idx
self.a_dim = var_dim
self.z_dim = constraint_num
self.dim = self.a_dim + self.z_dim + 1
def func(self, a_z_b):
return -a_z_b[self.a_dim + self.idx]
def grad(self, a_z_b):
res = np.zeros(self.dim)
res[self.a_dim + self.idx] = -1
return res
def hessian(self, a_z_b):
return np.zeros((self.dim, self.dim))
def create(self, idx):
return self.function_nonnegative_constraint(self.var_dim,
self.constraint_num, idx)
class svm:
def __init__(self, C, kernel_func):
self.C = C
self.K = None
self.X = None
self.kernel_func = kernel_func
def fit(self, X, y):
self.X = X
self.K = np.zeros((X.shape[0], X.shape[0]))
for i in range(X.shape[0]):
for j in range(i, X.shape[0]):
self.K[i, j] = self.kernel_func(X[i, :], X[j, :])
self.K[j, i] = self.K[i, j]
tdcc = train_data_constraint_creator(self.K, y)
ncc = nonnegative_constraint_creator(self.K, y)
Pro = ipm.problem(objective_function(self.K, self.C,
self.K.shape[1], y.shape[0]),
[tdcc.create(i) for i in range(y.shape[0])] +
[ncc.create(i) for i in range(y.shape[0])],
[])
a_z_b_opt, s_opt, u_opt, v_opt = ipm.interior_point(Pro, eps = 1e-15)
self.a = a_z_b_opt[:self.K.shape[0]]
self.b = a_z_b_opt[-1]
self.u = u_opt
self.support_vectors = np.where(self.u[:X.shape[0]] > 1e-8)
def predict(self, X_):
if X_ is None:
raise Exception('not learned')
K_ = np.zeros((self.K.shape[0], X_.shape[0]))
for i in range(self.K.shape[0]):
for j in range(X_.shape[0]):
K_[i, j] = self.kernel_func(self.X[i], X_[j])
return np.sign(self.a @ K_ + self.b)
def make_data(xmin, xmax, n):
np.random.seed(1111)
cnt = 0
t = 0
N = n * 9
X_train = np.zeros((N, 2))
y_train = np.zeros(N)
center1 = np.linspace(-(-xmin[0] + xmax[0])/3.5,
(-xmin[0] + xmax[0])/3.5, 3)
center2 = np.linspace(-(-xmin[1] + xmax[1])/3.5,
(-xmin[1] + xmax[1])/3.5, 3)
for i in range(3):
for j in range(3):
y_train[t:t + n] = 1 if cnt % 2 == 0 else -1
X_train[t:t + n, :] = 1.5 * np.random.randn(n, 2) + \
np.array([center1[i], center2[j]])
t += n
cnt += 1
T1 = np.linspace(xmin[0], xmax[0], 60)
T2 = np.linspace(xmin[1], xmax[1], 60)
X_test = np.zeros((T1.shape[0] * T2.shape[0], 2))
i = 0
for t1 in T1:
for t2 in T2:
X_test[i, :] = np.array([t1, t2])
i += 1
return X_train, y_train, X_test
def main():
xmin = np.array([-9, -9])
xmax = np.array([9, 9])
n = 20
X_train, y_train, X_test = make_data(xmin, xmax, n)
gamma = 0.5
kernel_func = lambda x_i, x_j : np.exp(-gamma * (x_i - x_j) @ (x_i - x_j))
sv = svm(0.1, kernel_func)
sv.fit(X_train, y_train)
y_pred = sv.predict(X_test)
y_pred[y_pred == -1] = 0
size = 100
sc = plt.scatter(X_test[:, 0], X_test[:, 1],
size, y_pred, marker = ',', cmap='bwr')
plt.scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1] ,
edgecolors="black", marker="o" , color='r',
label= "train data (y = 1)")
plt.scatter(X_train[y_train == -1, 0], X_train[y_train == -1, 1] ,
edgecolors="black", marker="o" , color='b',
label= "train data (y = -1)")
plt.axis([xmin[0], xmax[0], xmin[1], xmax[1]])
plt.legend(loc='upper left')
plt.show()
if __name__ == '__main__':
main()
先述のproblem
クラスとinterior_point
関数を使ってSVMを実装している。基本的にはこれまで述べてきたことをそのまま書いてるだけなので特筆すべきことはないと思う。これを実行すると以下の結果が得られる。
ガウスカーネルを使っているので非線形な分離できている。
ちなみに、グラム行列が正則でない場合(例えば、を小さくしすぎた場合や全く同じ学習サンプルが2つ以上ある場合などにあり得る)、主双対内点法のコードにおけるP
も正則でなくなり計算が破綻する。このあたりの弱点はSMO法などで克服される。
ここに挙げた主双対内点法のコードやSVMのコードに加えて、線形のSVMやその他の最適化問題で確認したコードもあるので、それらは全て以下のgithubリポジトリに上げている。
github.com
最後に
今回は制約付き非線形計画問題を解く主双対内点法を実装した。制約付き非線形計画問題を解く手法としては他にも逐次2次凸計画法や乗数法(拡張ラグランジュ関数法)などがある。ちなみに乗数法については、以前その改良の1つである交互方向乗数法(ADMM)を実装した記事を書いている。
yamagensakam.hatenablog.com
ということで、非線形計画問題については一旦ここまでにして次から整数計画問題に移る予定。
参考文献
- いつも通り今回の内容もこちらの書籍を主に参考にして書いた。
- SVMやカーネル法、リプレゼンター定理などについてはこちらの書籍がわかりやすい。