tl;dr
- Dual ascent method を振り返る
- Augmented Lagrangians method の概要について述べる
- toy problem に対して Augmented Lagrangians method を適用した結果について述べる
モチベ
前回の記事では，等式制約付き最適化問題に対するラグランジュ双対問題を導出し，双対問題を解くためのアルゴリズムである "dual ascent method" について概説した．
記事の中で数値実験として凸二次関数を線型の等式制約の下で最小化する問題を扱った．記事では主問題の目的関数と双対問題の目的関数が反復を重ねるごとに近くなっていく様子を提示したが，その結果はあくまでも上手くいったときの結果であり，初期点を変えたり，定数行列 $P,\, A$ や定数ベクトル $q,\, b$ を変えると途端に上手くいかない場合が多々ある．例えば，以下のような場合である．変更点は $m=5,\, n=20$ から $m=20,\, n=100$ にしただけである．
m = 10
n = 100
np.random.seed(1)
P = np.random.randn(n, n)
P = P.T @ P
q = np.random.randn(n)
A = np.random.randn(m, n)
b = np.random.randn(m)
最適化アルゴリズムの部分は 前回の記事 と同一であるため省略し，dual ascent method を適用して得られる各反復での主問題と双対問題の目的関数値の推移のみ以下に載せる．
グラフから明らかなように，反復を重ねると主問題の目的関数値（緑）は正の無限大に発散し，双対問題の目的関数値（赤）は負の無限大に発散してしまっている．
この記事では，dual ascent method を拡張した「拡張ラグランジュ法」と呼ばれる最適化アルゴリズムについて扱う．この手法は dual ascent method とほとんど同じ反復を行うアルゴリズムであるものの，より緩い仮定のもとで収束性を示すことができるなどの利点がある．収束解析については触れないが，数値実験として拡張ラグランジュ法が dual ascent method に比べて頑健に動作することを確認する．なお，本記事を通して，ベクトルは太字で記述しないので注意されたい．
拡張ラグランジュ関数
次の等式制約付き最適化問題 (P) を考える:
\begin{align} \mathrm{(P)}\qquad \begin{gathered} \min_{x\in\mathbb{R}^n} \quad &f(x) \\ \mathrm{s.t.} \quad &Ax = b. \end{gathered} \end{align}
ここで，$A\in\mathbb{R}^{m\times n}, \, b\in\mathbb{R}^{m}$ はそれぞれ適当な大きさの定数行列，ベクトルである．また，$f$ は凸関数とする． このような等式制約付きの凸計画問題に対して拡張ラグランジュ関数は次のように定義される．
拡張ラグランジュ関数における $\rho$ は正の値をとるが，$\rho=0$ としたときは通常のラグランジュ関数に対応する．また，拡張ラグランジュ関数 $L_\rho$ は次の最適化問題に対する（通常の）ラグランジュ関数と見ることもできる．
\begin{align} \min_{x\in\mathbb{R}^n} \quad &f(x) + \frac{\rho}{2}\|Ax-b\|^2_2 \\ \mathrm{s.t.} \quad &Ax = b. \end{align}
なお，この最適化問題は (P) と本質的に等価である．実際，任意の実行可能解 $\hat{x}$ に対し，目的関数に追加した項の値は $\frac{\rho}{2}\|A\hat{x}-b\|^2_2 = 0$ である．
拡張ラグランジュ法
前節で定義した拡張ラグランジュ関数を用いて，凸計画問題 (P) に対し，以下の反復を繰り返す手法を拡張ラグランジュ法（Augmented Lagrangian method; method of multipliers）という．
\begin{align} x^{k+1} &:= \arg\min_x L_\rho(x, y^k) \\ y^{k+1} &:= y^k + \rho (Ax^{k+1} - b). \end{align}
この反復は dual ascent method とほとんど同じであり，違いは $x$ の更新式において拡張ラグランジュ関数を用いている点と，$y$ の更新においてステップサイズ $\alpha^k$ の代わりに $\rho$ が用いられている点のみである．
例：凸二次計画問題に対する拡張ラグランジュ法
凸計画問題における目的関数 $f$ が $f(x)=\frac{1}{2}x^\top Px + q^\top x$ で与えられる場合を考える．ここで，$P$ は半正定値対称行列とする．この凸計画問題に対する拡張ラグランジュ関数 $L_\rho$ は
$$ L_\rho (x, y) = \frac{1}{2}x^\top Px + q^\top x + y^\top (Ax-b) + \frac{\rho}{2}\| Ax-b\|^2_2 $$
で与えられる．この $L_\rho$ は $x$ に関して凸であるから，拡張ラグランジュ法における $x$ の更新では，方程式 $\nabla_x L_\rho (x, y^k)$ を解いた解を $x^{k+1}$ とすればよい．拡張ラグランジュ関数の $x$ に関する勾配は
\begin{align} \nabla_x L_\rho (x, y^k) = Px + q + A^\top y^k + \rho (A^\top Ax - A^\top b) \end{align}
で与えられる．したがって，$\nabla_x L_\rho (x, y^k) = 0$，すなわち，正規方程式
\begin{align} (P + \rho A^\top A)x = -q + A^\top (\rho b - y^k) \end{align}
を解いて得られる解を $x^{k+1}$ とすればよい．
数値実験
前節の凸二次計画問題を扱う．以下で numpy と matplotlib を用いるため import しておく（省略）．
拡張ラグランジュ関数 $L_\rho$ と拡張ラグランジュ法の $x$ を更新する部分は次のように書ける．ただし，目的関数が凸二次関数の場合に限る．なお，f や他の定数については 前回の記事 を参照されたい．
def L_rho(x, y, rho):
return f(x) + y.T@(A@x-b) + 0.5 * rho * np.linalg.norm(A@x - b) ** 2
def argminLrho(y, rho):
return np.linalg.solve(P + rho * A.T@A, rho * A.T@b - A.T@y - q)
拡張ラグランジュ法は次のように書ける．引数として，初期点 y0，ペナルティパラメータ rho，反復回数 itermax を設定している．また，返り値として各反復における主問題の目的関数値 f_log，双対問題の目的関数値 g_log，残差 $\|Ax^k - b\|_2$ の値 residue を設定している．
def Augmented_Lagrangians_method(y0, rho, itermax=10):
yk = y0
f_log = np.zeros(itermax)
g_log = np.zeros(itermax)
residue = np.zeros(itermax)
for k in range(itermax):
xk = argminLrho(yk, rho)
y = yk + rho * (A@xk - b)
f_log[k] = f(xk)
g_log[k] = g(y)
residue[k] = np.linalg.norm(A@xk - b)
yk = y
return f_log, g_log, residue
まず，冒頭で述べた dual ascent method で上手くいかなかった例について実験してみる．適当な初期点の下で Augmented_Lagrangians_method() を呼び出せばよい．ここでは，初期点を標準正規分布から生成し，ペナルティパラメータとして $\rho=1$ と設定した．
f_log, g_log, residue =\
Augmented_Lagrangians_method(y0=np.random.randn(m), rho=1, itermax=10)
結果として，f_log と g_log をプロットした図を載せる．この図を見ると，僅か 4 反復程度で主問題の目的関数値と双対問題の目的関数値がほとんど一致していることが確認できる．
また，各反復における残差 $\|Ax^k - b\|_2$ は以下のようになる．縦軸は log スケールで表示しており，反復によって残差は指数的に減少していくことが確認できる．
蛇足
拡張ラグランジュ法について扱いました．今回は目的関数が凸二次関数というシンプルな形で記述される関数を扱ったので，$x$ の更新を最適化問題を解くことなく行うことができました．しかしながら，一般には $\arg\min_x L_\rho(x, y^k)$ を解く必要があるので毎回の反復で凸計画問題を扱う必要が出てきます．なんかそういうインスタンスを見つけて実験してみたいですね．
