非負制約を含む最適化問題に対する KKT 条件
表題の通りです．以下にまとめました．
『しっかり学ぶ数理最適化』演習問題（第 3 章）
『しっかり学ぶ数理最適化』第 3 章「非線形計画」の演習問題に関する個人的なメモ書きです．
3.1
二次正方行列の正定値性に関する等価な条件の導出の問題です．計算しましょう．
3.2-3.5
関数の凸性に関する問題です．3.2 は平均二乗誤差 $f(a, b) = \frac{1}{n} \sum_{i=1}^{n} (y_i - ax_i - b)^2$ の凸性を示す問題です．関数 $f$ の引数が $(a, b)$ であることに注意して，ヘッセ行列を計算し，その半正定値性を見てあげることによって示すことができます．具体的には，計算したヘッセ行列の二次形式を考えてあげて，その非負性を示せば十分です．3.3 は次の関数: \begin{align} f(x, y) = \max_{i\in [n]} \sqrt{(x_i - x)^2 + (y_i - y)^2} \label{1}\tag{1} \end{align} に関する凸性を示す問題です．ここで，$[n]$ は集合 $\{1,2,\dots,n\}$ を表すものとします．まず前提として，凸集合 $S$ 上で定義される凸関数 $f_1,\,f_2$ に対して，関数 $\max \{ f_1, \, f_2\}$ はまた凸関数であることが成り立ちます．この事実を繰り返し用いれば，$n$ 個の凸関数に対する $\max$ 値関数も凸関数であり，したがって，式 \eqref{1} で定義される関数の凸性を示すためには，$f_i (x,y) := \sqrt{(x_i - x)^2 + (y_i - y)^2}$ が凸関数であることを示せばよいです．書籍の解答では，3.2 と同様にヘッセ行列を計算して，その半正定値性を見ています．3.4 は $f(x) = \| x \|_2$ の凸性を示す問題で，書籍では勾配を用いた凸関数の特徴付けを用いて示しています．より簡単に，任意の $x,\,y\in\mathbb{R}^n$ と任意の実数 $\lambda\in [0,1]$ に対して，$\| \lambda x + (1-\lambda)y \|_2 \leqq \lambda\| x \|_2 + (1-\lambda)\| y \|_2$ が成立することを示してもよいと思います．これは，ノルムの性質より直ちに従います．3.5 はいわゆる log sum exp 関数と呼ばれる次の関数: \begin{align} f(x) = \log \left( \sum_{i=1}^{n} \mathrm{e}^{x_i} \right) \label{2}\tag{2} \end{align} の凸性を示す問題です．この関数は様々な分野で登場する関数で Wikipedia にもページがあります（LogSumExp）．特徴を一つ述べておくと，関数 $f$ の $x_i$ に関する偏微分は \begin{align} \frac{\partial f(x)}{\partial x_i} = \frac{\mathrm{e}^{x_i}}{\sum_{j}\mathrm{e}^{x_j}} \end{align} となり，よく知られた softmax 関数になります．また，本書でも述べられているように $\max\{x_1,x_2,\dots,x_n\}$ を近似する場面でもよく使われます．さて，この関数の凸性ですが，本書ではやはりヘッセ行列を計算し，その半正定値性を見ています．また，こちらのページ で議論されているように，Hölder の不等式に帰着させて示すのも綺麗だと思いました．
3.6
面白い問題です．与えられた非凸最適化問題を上手く変形することによって等価な凸計画問題を導出する問題です．
3.7
無制約二次計画問題を勾配降下法によって解くアルゴリズムにおいて，勾配降下方向にどれくらい進むのかについて決定する最適化問題を解く問題です．具体的には，アルゴリズムの第 $k$ 反復における近似解 $x^{(k)}$，勾配降下方向を $d(x^{(k)})$ としたとき，ステップサイズ $\alpha (>0)$ を変数とする関数 $g(\alpha) := f(x^{(k)} + \alpha d(x^{(k)}))$ を最小化する最適化問題を解く問題です．
3.8
無制約凸二次計画問題をニュートン法で実際に解く問題です．実際に解くといっても，実装する必要はなく，解析的に計算することができる問題です．実際に勾配ベクトル，ヘッセ行列を計算し，与えられた初期点から勾配降下方向に点を移動させてみると更新後の点における勾配ベクトルは $0$ ベクトルとなり，最適性条件を満たすことが分かります．
3.9
制約付き最適化問題を解く問題です．最適性条件を考えることによって解くことができます．
3.10
二次計画問題の双対問題を導出する問題です．本書では，行列 $Q$ が正則行列となっていますが，正定値対称行列が正しいです．対称行列でない場合は，本書 pp.91 と同様の議論によって一般性を失うことなく対称行列に変換することができます．改めて本問題で検討する二次計画問題は以下のように記述されます:
\begin{align} \min_{x} &\quad \frac{1}{2} x^\top Qx + c^\top x \\ \mathrm{s.t.} &\quad Ax = b, \, x\geqq 0. \end{align}
ここで，$Q$ は $n$ 次正定値対称行列，$c\in\mathbb{R}^{n},\, A\in\mathbb{R}^{m\times n},\, b\in\mathbb{R}^{m}$ は定数ベクトル，定数行列とします．
はじめに，双対問題の導出について簡単に述べておきます．与えられた最適化問題（"元問題"と呼ぶことにします，最小化問題とします）の一部の制約条件を緩和して，緩和した制約をペナルティ項として目的関数に組み込んだラグランジュ緩和問題を考えます．すると，ラグランジュ緩和問題の実行可能領域は元問題の実行可能領域よりも広がっており，またペナルティ項が非負であることから，ラグランジュ緩和問題の最適値は元問題の下界を達成することが確認されます．このラグランジュ緩和問題の最適値関数（決定変数はラグランジュ乗数となります）を最大化する問題を考えることによって，元問題の最適値とラグランジュ緩和問題の最適値が一致するようにすることが期待されます．この，ラグランジュ緩和問題の最適値関数を最大化する問題のことをラグランジュ双対問題と呼びます．[双対問題の導出の方法に関する記述おわり]
さて，3.10 の問題に戻りまして，本問題で考える二次計画問題に対するラグランジュ緩和問題を考えてみましょう．まず，制約条件は二つあり，等式制約 $Ax=b$ と決定変数 $x$ に関する非負条件です．等式制約に対するラグランジュ乗数を $u \in \mathbb{R}^{m}$，非負制約に対するラグランジュ乗数を $v\in\mathbb{R}^{n}_{\geqq 0}$ とします．すると，ラグランジュ関数は
\begin{align} L(x, u, v) &:= \frac{1}{2} x^\top Qx + c^\top x + u^\top (Ax-b) - v^\top x \\ &= \frac{1}{2} x^\top Qx + (c + A^\top u - v)^\top x - b^\top u \end{align}
となります．これを $x$ の非負条件の下で最小化する問題（ラグランジュ緩和問題）は次のように定式化されます:
\begin{align} \min_{x} &\quad \frac{1}{2} x^\top Qx + (c + A^\top u - v)^\top x - b^\top u \\ \mathrm{s.t.} &\quad x\in\mathbb{R}^n. \end{align}
いま，行列 $Q$ は正定値対称行列であるから，ラグランジュ緩和問題の目的関数は凸関数であり，その最適解を $\overline{x}$ とすれば，$\overline{x}$ は最適性条件:
\begin{align} \nabla_x L(\overline{x}, u, v) = Q\overline{x} + c + A^\top u - v = 0 \label{3}\tag{3} \end{align}
を満たします．なお，$\frac{1}{2} x^\top Qx$ の $x$ に関する勾配ベクトルが $Qx$ となるのは，$Q$ が対称行列だからです．式 \eqref{3} は $Q$ が正定値行列であるため，
\begin{align} \overline{x} = Q^{-1} (v - c - A^\top u) \end{align}
と表せます．この $\overline{x}$ を用いることによって，ラグランジュ緩和問題の最適値関数は陽に表すことができて，それを最大化する問題は
\begin{align} \max_{u,v} &\quad \frac{1}{2} (v - c - A^\top u)^{\top} Q^{-1} (v - c - A^\top u) - b^\top u \\ \mathrm{s.t.} &\quad v\geqq 0 \end{align}
であり，これが求めるラグランジュ双対問題です．
3.11
遂次二次計画法における Powell の修正 BFGS 公式を用いて近似行列 $B_k$ を更新したとき，$B_k$ が正定値ならば更新後の近似行列 $B_{k+1}$ も正定値となることを示す問題です．
参考
実行列のフロベニウスノルムのメモ
フロベニウスノルム
行列 $\mathbf{A}\in\mathbb{R}^{m\times n}$ に対して，その全成分の二乗和に対してルートをとったものをフロベニウスノルムといい，$\|\mathbf{A}\|_F$ で表す．
また，次に述べるように $\|\mathbf{A}\|_F$ はトレースを用いて次のように表すことができる：
$$ \|\mathbf{A}\|_F = \sqrt{\mathrm{tr}(\mathbf{A}^\top\mathbf{A})}. $$
証明
\begin{align} \|\mathbf{A}\|^2_F &\underset{\text{(1)}}{=} \sum_{i=1}^{m}\sum_{j=1}^{n}a_{ij}^2 \\ &\underset{\text{(2)}}{=} \sum_{j=1}^{n}\|\mathbf{a}_j\|^2_2 \\ &\underset{\text{(3)}}{=} \sum_{j=1}^{n}\left<\mathbf{a}_j,\mathbf{a}_j\right> \\ &\underset{\text{(4)}}{=}\sum_{j=1}^{n}(A^\top A)_{jj} \\ &\underset{\text{(5)}}{=} \mathrm{tr}(A^\top A). \end{align}
$\|\mathbf{x}\mathbf{x}^\top\|_F$ の導出
tl;dr
- $\mathbf{x}\mathbf{x}^\top$ のフロベニウスノルム（Frobenius norm）を導く
- 体 $K$ として，$K=\mathbb{R}^{m\times n}$ を考える
フロベニウスノルム
行列 $\mathbf{A}\in\mathbb{R}^{m\times n}$ に対して，その全成分の二乗和に対してルートをとったものをフロベニウスノルムといい，$\|\mathbf{A}\|_F$ で表す．
また，次に述べるように $\|\mathbf{A}\|_F$ はトレースを用いて次のように表すことができる：
$$ \|\mathbf{A}\|_F = \sqrt{\mathrm{tr}(\mathbf{A}^\top\mathbf{A})}. $$
証明は 行列のフロベニウスノルムとその性質 を参照． こちら に証明を記載しました．
導出
\begin{align} \|\mathbf{xx}^\top\|_F &= \sqrt{\mathrm{tr}((\mathbf{xx}^\top)^\top\mathbf{xx}^\top)} \\ &= \sqrt{\mathrm{tr}(\mathbf{x}\mathbf{x}^\top\mathbf{xx}^\top)} \\ &= \sqrt{\|\mathbf{x}\|^2_2\cdot\mathrm{tr}(\mathbf{x}\mathbf{x}^\top)} \\ &= \sqrt{\|\mathbf{x}\|^2_2\cdot \|\mathbf{x}\|^2_2} \\ &= \|\mathbf{x}\|^2_2. \end{align}
Augmented Lagrangians method（拡張ラグランジュ法）
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$ は凸関数とする． このような等式制約付きの凸計画問題に対して拡張ラグランジュ関数は次のように定義される．
定義
等式制約付き凸計画問題 (P) に対して，
\begin{align}
L_\rho(x, y) = f(x) + y^\top (Ax - b) + \frac{\rho}{2}\| Ax-b\|^2_2
\end{align}
で定義される関数 $L_\rho$ を拡張ラグランジュ関数（augmented Lagrangians function）という．ここで，$\rho > 0$ はペナルティである．
拡張ラグランジュ関数における $\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)$ を解く必要があるので毎回の反復で凸計画問題を扱う必要が出てきます．なんかそういうインスタンスを見つけて実験してみたいですね．
