はじめに #
等式制約あり最適化問題を扱う拡張ラグランジュ法 (augmented Lagrangian method) をPythonで実装しました。
拡張ラグランジュ法 #
以下の等式制約つき最適化問題を考えます。
$$ \min \ f(x) \\\ \mathrm{s.t.} \ g_i (x) = 0 \ (i=1, 2, ..., N) $$\(f(x), g_i(x)\)は凸とします。この最適化問題に対し、次式の拡張ラグランジュ関数\(L_A(x)\)を定義します。
$$ L_A(x) = f(x) + \sum_{i=1}^{N} \lambda_i g_i(x) + \frac{\rho}{2} \sum_{i=1}^{N}g_i(x)^2 $$ここで、\(\lambda_i \ (i=1, 2, ..., N)\)はラグランジュ乗数、\(\rho\)はペナルティ項の係数です。
拡張ラグランジュ法のアルゴリズムは以下のようになります。
Step 1. \(\lambda_i, \rho\)の初期値を設定する。
Step 2. \(L_A(x)\)を最小化する\(x\)の値\(x^*\)を求める。最小化する手法は何でもよい。
Step 3. \(x^*\)がKKT条件を満たせば終了。そうでなければStep 4に進む。
Step 4. \(\lambda_i, \rho\)を次式で更新し、Step 2に戻る。
$$\lambda_i \leftarrow \lambda_i + \rho g_i(x) \\\ \rho \leftarrow \alpha \rho \ (\alpha > 1)$$Step 2で拡張ラグランジュ関数\(L_A(x)\)を最小化する手法は、制約なし最適化問題を解ける手法であれば何でも構いません（例：再急降下法、ニュートン法）。 また、Step 3のKKT (Karush-Kuhn-Tucker) 条件とは、最適解が満たすべき条件のことです。ここで、通常の（ペナルティ項のない）ラグランジュ関数は次式となります。
$$ L(x, \lambda) = f(x) + \sum_{i=1}^{N} \lambda_i g_i(x)$$\(L(x, \lambda)\)を用いると、等式制約つき最小化問題のKKT条件は以下になります。
$$\nabla_x L(x, \lambda) = \nabla f(x) + \sum_{i=1}^{N} \lambda_i \nabla g_i(x) = 0 \\\ \nabla_\lambda L(x, \lambda) = g_i(x) = 0 \ (i=1, 2, ..., N)$$すなわち、ラグランジュ関数の\(x\)と\(\lambda\)に関するそれぞれの勾配が0になると収束したことになります。実装上は、許容誤差以下になれば収束したと判定します。
背景 #
拡張ラグランジュ法の背景を解説します。制約つき最適化問題を解く方法として、以下のペナルティ法があります。
$$ \min \ f(x) + \frac{\rho}{2} \sum_{i=1}^{N}g_i(x)^2 $$係数\(\rho\)の値を無限大に近づけることによって、\(g_i(x)=0\)の制約を満たすことができますが、\(\rho\)の値が大きくなると数値的に不安定になる課題があります。
一方、拡張ラグランジュ法では、係数\(\rho\)の値が無限でなくとも、ある有限の大きな値であれば、制約を満たす最適解が得られます。詳細は以下の書籍に記載されています。
Pythonによる実装 #
拡張ラグランジュ法をPythonで実装します。拡張ラグランジュ関数の最小化には再急降下法を使用します。バージョンは以下の通りです。
| ソフトウェア | バージョン |
|---|---|
| python | 3.9.7 |
| numpy | 1.20.3 |
| pandas | 1.3.4 |
例題 #
以下の最小化問題を例題とします。
$$\min \ x_0^2 + 2x_1^2 + x_2^2 \\\ \mathrm{s.t.} \ g_0(x) = 2x_0 + x_1 - 5 = 0 \\\ \ g_1(x) = x_1 + 2x_2-3=0$$評価関数の勾配と、等式制約のヤコビアンはそれぞれ以下となります。
$$ \nabla f(x) = \begin{bmatrix} 2x_0 & 4x_1 & 2x_2 \end{bmatrix}^\top \\\ \nabla g(x)^\top = \begin{bmatrix} \nabla g_0(x)^\top \\\ \nabla g_1(x)^\top \end{bmatrix} = \begin{bmatrix} 2 & 1 & 0 \\\ 0 & 1 & 2 \end{bmatrix} $$拡張ラグランジュ関数\(L_A(x)\)を再急降下法で最小化するため、\(L_A(x)\)の勾配も計算できるようにします。
$$L_A(x) = f(x) + \sum_{i=1}^{N} \lambda_i g_i(x) + \frac{\rho}{2} \sum_{i=1}^{N}g_i(x)^2$$であるため、勾配は次式となります。
$$\nabla L_A(x) = \nabla f(x) + \sum_{i=1}^{N} \lambda_i \nabla g_i(x) + \rho \sum_{i=1}^{N}g_i(x) \nabla g_i(x) \\\ = \nabla f(x) + \lambda \nabla g(x)^\top + \rho g(x) \nabla g(x)^\top $$関数とクラスの定義 #
まず、以下をクラスや関数で定義します。
- 再急降下法
- 評価関数とその勾配
- 等式制約とそのヤコビアン
- 拡張ラグランジュ関数とその勾配
import numpy as np
import pandas as pd
class GradientDescent:
"""再急降下法"""
def __init__(self, fun, der, xi=0.3, tau=0.9, tol=1e-4, ite_max=20):
self.fun = fun # 目的関数
self.der = der # 関数の勾配
self.xi = xi # Armijo条件の定数
self.tau = tau # 方向微係数の学習率
self.tol = tol # 勾配ベクトルのL2ノルムがこの値より小さくなると計算を停止
self.path = None # 解の点列
self.ite_max = ite_max # 最大反復回数
def minimize(self, x):
path = [x]
for i in range(self.ite_max):
grad = self.der(x)
if np.linalg.norm(grad, ord=2)<self.tol:
break
else:
beta = 1
while self.fun(x-beta*grad) > (self.fun(x)-self.xi*beta*np.dot(grad,grad)):
# Armijo条件を満たすまでループする
beta *= self.tau
x -= beta*grad
path.append(x)
self.opt_x = x # 最適解
self.opt_result = self.fun(x) # 関数の最小値
self.path = np.array(path) # 探索解の推移
return x
def objective(x):
"""評価関数"""
return x[0]**2 + 2*x[1]**2 + x[2]**2
def grad_objective(x):
"""評価関数の勾配"""
return np.array([2*x[0], 4*x[1], 2*x[2]], dtype=float)
def constraint(x):
"""等式制約"""
return np.array([2*x[0]+x[1]-5, x[1]+2*x[2]-3])
def grad_constraint(x):
"""等式制約のヤコビアン"""
return np.array([[2, 1, 0],
[0, 1, 2]], dtype=float)
class AugLag:
"""拡張ラグランジュ関数"""
def __init__(self, obj, const, lambd, rho):
self.obj = obj # 評価関数を求める関数
self.const = const # 等式制約を求める関数
self.lambd = lambd # ラグランジュ乗数
self.rho = rho # ペナルティ係数
def __call__(self, x):
c = self.const(x)
return self.obj(x) + sum(self.lambd*c) + self.rho*sum(c**2)/2
class GradAugLag:
"""拡張ラグランジュ関数の勾配"""
def __init__(self, grad_obj, const, grad_const, lambd, rho):
self.grad_obj = grad_obj # 評価関数の勾配を求める関数
self.const = const # 等式制約を求める関数
self.grad_const = grad_const # 等式制約のヤコビアンを求める関数
self.lambd = lambd # ラグランジュ乗数
self.rho = rho # ペナルティ係数
def __call__(self, x):
c = self.const(x)
dc = self.grad_const(x)
return self.grad_obj(x) + np.dot(self.lambd,dc) + self.rho*np.dot(c,dc)
Pandasは計算結果の格納にのみ使用しており、実質的にはNumPyのみ最適化に使用しています。
再急降下法は関数の勾配の方向に解を更新する最適化手法です。詳細は以下の記事を参考にしてください。 直線探索を使った最急降下法をPythonで実装 - Helve Tech Blog
拡張ラグランジュ関数とその勾配については、GradientDescentクラスの内部でfun(x), der(x)というように関数の形式で値を計算したいため、__call__()メソッドを実装しています。
最適化計算の実行 #
上記のクラス・関数を使った拡張ラグランジュ法による最適化のスクリプトは以下になります。ただし、コードの分かりやすさを優先しているため、同じ条件で評価関数や等式制約を何度も計算しているなど、計算効率は悪いです。
# 拡張ラグランジュ法に関するパラメータ
rho0 = 5 # ペナルティ係数の初期値
aplha = 1.5 # ペナルティ係数を更新する係数
tol = 1e-4 # KKT条件の許容誤差
lambda0 = np.array([1, 1], dtype=float) # ラグランジュ乗数（制約の数と同じ次元のベクトル）
max_iter = 100 # 最大反復回数
x_init = np.array([0.0, 0.0, 0.0], dtype=float) # 解の初期値
x = x_init
lambd = lambda0
rho = rho0
gal = GradAugLag(grad_objective, constraint, grad_constraint, lambd, rho)
columns = ["obj","x0","x1","x2","dldx0","dldx1","dldx2","g0","g1","lambda0","lambda1","rho"]
result = pd.DataFrame([[objective(x)]+x.tolist()+gal(x).tolist()+constraint(x).tolist()\
+lambd.tolist()+[rho]], columns=columns, index=[0])
for i in range(max_iter):
al = AugLag(objective, constraint, lambd, rho)
gal = GradAugLag(grad_objective, constraint, grad_constraint, lambd, rho)
gd = GradientDescent(al, gal)
x = gd.minimize(x) # 再急降下法による解の更新
grad_lag = GradAugLag(grad_objective, constraint, grad_constraint, lambd, rho=0)
if max(abs(grad_lag(x)))<tol and max(abs(constraint(x)))<tol:
# KKT条件による収束判定
break
# ラグランジュ乗数・ペナルティ係数の更新
lambd += rho*constraint(x)
rho *= aplha
temp = pd.DataFrame([[objective(x)]+x.tolist()+gal(x).tolist()+constraint(x).tolist()\
+lambd.tolist()+[rho]], columns=columns, index=[i+1])
result = pd.concat([result, temp])
print(f"評価関数：{objective(x)}")
print("最適解：")
for i in range(len(x)):
print(f"x{i}: {x[i]}")
収束判定として、ラグランジュ関数のヤコビアンgrad_lag()と等式制約constraint()の全ての要素の絶対値が、許容誤差tol以下となることを条件にしています。
実行結果
最適解と評価関数は以下のようになります（実行時間は1秒程度です）。また、resultには各ステップの評価関数の値や解の値などを格納しています。
評価関数：6.900124790320017
最適解：
x0: 2.099936900209656
x1: 0.8001650843826812
x2: 1.0999370349831807
Pyomoによる検算 #
同じ問題をPyomoと非線形ソルバIPOPTで解いて検算を行います。
import pyomo.environ as pe
m = pe.ConcreteModel()
m.x = pe.Var(range(3), domain=pe.Reals)
m.obj = pe.Objective(expr = m.x[0]**2 + 2*m.x[1]**2 + m.x[2]**2)
m.const0 = pe.Constraint(expr = 2*m.x[0] + m.x[1] == 5)
m.const1 = pe.Constraint(expr = m.x[1] + 2*m.x[2] == 3)
opt = pe.SolverFactory('ipopt')
res = opt.solve(m)
print(f"評価関数：{m.obj()}")
for i in range(3):
print(f"x{i}: {m.x[i]()}")
実行結果
評価関数：6.9
x0: 2.1
x1: 0.8
x2: 1.0999999999999999
拡張ラグランジュ法と同じ評価関数と最適解が得られました。
参考文献 #
実装にあたって、以下のサイト・書籍を参考にさせていただきました。
