tl;dr
- 共役関数を定義した後，等式制約付き最適化問題に対するラグランジュ双対問題を導出する
- 双対問題を解く "dual ascent method" について概説する
- toy problem に対して dual ascent method を適用した結果について述べる
共役関数
上の定義では $f$ の共役関数を $\sup$ によって定義したが，次のように $\inf$ によって定義することもできる：
\begin{align} f^\star (\mathbf{\xi}) := -\inf_{x\in\mathbb{R}^n} \, \left\{ f(\mathbf{x}) - \left< \mathbf{\xi}, \mathbf{x} \right> \right\}. \end{align}
実際，以下のように確かめることができる：
\begin{align} f^\star (\mathbf{\xi}) &:= \sup_{x\in\mathbb{R}^n} \, \left\{ \left<\mathbf{\xi}, \mathbf{x} \right> - f(\mathbf{x}) \right\} \\ &= -\inf_{x\in\mathbb{R}^n} \, -\left\{ \left<\mathbf{\xi}, \mathbf{x} \right> - f(\mathbf{x}) \right\} \\ &= -\inf_{x\in\mathbb{R}^n} \, \left\{ f(\mathbf{x}) - \left<\mathbf{\xi}, \mathbf{x} \right> \right\}. \end{align}
ここで，一般に $\sup\, f(\mathbf{x}) = - \inf\, \{-f(\mathbf{x})\}$ であることを用いた．
共役関数の例
例えば，分離可能な凸関数 $f(x_1, x_2) = g(x_1) + h(x_2)$ の共役関数は，
$$ f^\star ( \xi_1, \xi_2 ) = g^\ast(\xi_1) + h^\ast(\xi_2) $$
で与えられる．
共役関数を用いた双対問題の導出
次の等式制約付き最適化問題を考える:
\begin{align} \min_{x\in\mathbb{R}^n} \quad &f(x) \\ \mathrm{s.t.} \quad &Ax = b. \end{align}
ここで，$A\in\mathbb{R}^{m\times n}, \, b\in\mathbb{R}^{m}$ はそれぞれ適当な大きさの定数行列，ベクトルである．また，$f$ は凸関数とする．
$$ L(x, y) = f(x) + y^\top (Ax - b) $$
で定義される．ここで，ベクトル $y\in\mathbb{R}^{m}$ はラグランジュ乗数を表す．双対関数は，
\begin{align} g(y) &= \inf_{x\in\mathbb{R}^n} L(x, y) \\ &= \inf \left\{ f(x) + y^\top (Ax - b) \right\} \\ &= \inf \left\{ f(x) - \left< -Ax+b, y \right>\right\} \\ &= \inf \left\{ f(x) - \left< -Ax, y \right> \right\} - \left<b, y \right> \\ &= \inf \left\{ f(x) - \left< -A^\top y, x \right> \right\} - \left<b, y \right> \\ &= -f^\star (-A^\top y) - \left<b, y \right> \end{align}
で与えられる．これより，主問題に対する双対問題は
\begin{align} \max \quad &g(y) \\ \mathrm{s.t.} \quad &y\in\mathbb{R}^m \end{align}
で与えられる．
いま，双対問題の最適解を $y^\star$ とする．このとき，主問題の最適解は $y^\star$ の情報を用いて次のように求めることができる：
\begin{align} x^\star = \arg\min_x L(x, y^\star). \end{align}
ただし，ラグランジュ関数 $L(x, y^\star)$ の最適解は一意であるときに限る（例えば，主問題の目的関数 $f$ が強凸関数であるとき）．
dual ascent method
双対問題を解くことを考える．ここでは，双対問題の目的関数 $g$ は微分可能であるとする．
まず，$x^+$ をラグランジュ関数を最小にする点として求める，すなわち， $$ x^+ = \arg\min_x L(x, y). $$
すると，関数 $g$ の点 $y$ における勾配ベクトルは $\nabla g(y)=Ax^+-b$ で与えられる．実際，$x^+ = \arg\min_x L(x, y)$ であるとき，
$$ g(y) = \inf_x L(x, y) = L(x^+, y) $$
であり，両辺 $y$ について微分すると，
\begin{align} \nabla_y g(y) &= \nabla_y L(x^+, y) \\ &= \nabla_y \left\{ f(x^+) + y^\top (Ax^+ - b) \right\} \\ &= Ax^+ - b \end{align}
を得る．これは，主問題の等式制約の残差に相当する．
dual ascent method は次のような更新則で構成される：
\begin{align} x^{k+1} &:= \arg\min_x L(x, y^k) \\ y^{k+1} &:= y^k + \alpha^k (Ax^{k} - b). \end{align}
ここで，$\alpha^k > 0$ は第 $k$ 反復目におけるステップサイズを表す．この手法が "dual ascent" と呼ばれる所以は，ステップサイズ $\alpha^k$ を適切に選ぶことによって，$g(y^{k+1}) > g(y^k)$ が成り立つことに拠る．
数値実験
折角なので dual ascent method を実装して数値実験してみる．解くべき最適化問題として次のような二次計画問題を考える：
\begin{align} \min_{x\in\mathbb{R}^n} \quad &\frac{1}{2}x^\top Px + q^\top x \\ \mathrm{s.t.} \quad &Ax = b, \end{align}
ここで，$P$ は $n$ 次対称定数行列，$q\in\mathbb{R}^n,\, A\in\mathbb{R}^{m\times n},\, b\in\mathbb{R}^{m}$ はそれぞれ適当な大きさの定数ベクトルとする．
以下で numpy と matplotlib を用いるため import しておく．
import numpy as np
import matplotlib.pyplot as plt
それぞれの行列，ベクトルは次のように生成しておく．
m = 5
n = 20
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)
また，必要な関数を次のように定義しておく．
def f(x):
return 0.5 * np.dot(x, P@x) + q.T@x
def L(x, y):
return f(x) + y.T@(A@x-b)
def argminL(y):
return np.linalg.solve(P, -q-A.T@y)
def g(y):
x = argminL(y)
return L(x, y)
さて，dual ascent method の反復は以下のように記述できる．
yk = np.random.randn(m)
maxiter = 1000
alpha = 0.001
f_log = []
g_log = []
for k in range(maxiter):
xk = argminL(yk)
y = yk + (A@xk - b) * alpha
f_log.append(f(xk))
g_log.append(g(y))
yk = y
ここでは，ステップサイズとして定数ステップ $\alpha^k = \frac{1}{1000}$ を用いた．また，f_log, g_log にはそれぞれ暫定解 $x^k$ における主問題の目的関数値，$y^k$ における双対問題の目的関数値を格納する．
f_log, g_log を図示すると次のようになる．
なお，可視化コードは以下．
fig = plt.figure()
ax = fig.add_subplot(111)
cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
ax.plot(f_log, label="f", color=cycle[2])
ax.plot(g_log, label="g", color=cycle[3])
ax.legend()
ax.set_xlabel("iter")
ax.set_ylabel("the value of objective function")
ax.grid(True)
他にも $Ax^k - b$ の値の推移など興味あるが，また追記したいと思う．
