[非線形計画問題トップに戻る](https://hackmd.io/@nagae/NLP)

## Frank-Wolfe法
**線形不等式制約** および非負制約つき非線形最適化問題：
$$\begin{align}
\min_{\boldsymbol{x}} & f(\boldsymbol{x}) \\
\text{s.t.}{} & \boldsymbol{A} \boldsymbol{x} \geq \boldsymbol{b}\\
& \boldsymbol{x} \geq \boldsymbol{0}
\end{align}$$
に対する解法の一つとして知られる **Frank-Wolfe法** (フランク-ヴォルフ法) は，ある解の周囲で目的関数を線形近似して得られる問題(**補助問題**)が**線形計画問題**となることを活用し，もとの解と補助解との**凸結合(線形結合)** として解を改訂する方法である．

暫定解$\boldsymbol{x}$の周辺で目的関数を Taylor 展開し，二次以上の項を無視すれば，ある解$\boldsymbol{y}$における目的関数の値は
$$\begin{align}
f(\boldsymbol{y}) \approx f(\boldsymbol{x}) + (\nabla f(\boldsymbol{x}))^{\top} (\boldsymbol{y}-\boldsymbol{x}) = \hat{f}(\boldsymbol{y})
\end{align}
$$と近似できる．$f(\boldsymbol{x})$および$(\nabla f(\boldsymbol{x}))^{\top}\boldsymbol{x}$は与件であることから，$\hat{f}(\boldsymbol{y})$を最小化する問題は，以下の線形計画問題：
$$\begin{align}
\min_{\boldsymbol{y}} & (\nabla f(\boldsymbol{x}))^{\top} \boldsymbol{y}\\
\text{s.t.}{} & \boldsymbol{A} \boldsymbol{y} \geq \boldsymbol{b}\\
& \boldsymbol{y} \geq \boldsymbol{0}
\end{align}$$
として定式化できる．この問題を**補助問題**(auxiliary problem)と呼び，その解$\boldsymbol{y}^{\ast}$を**補助解**(auxiliary solution)と呼ぶ．

Frank-Wolfe法では，こうして得られた補助解$\boldsymbol{y}$ともとの解$\boldsymbol{x}$との線形結合：
$$
(1-\alpha) \boldsymbol{x} + \alpha \boldsymbol{y}, \quad 0\leq \alpha \leq 1
$$あるいは，解の降下方向を$\boldsymbol{d} = \boldsymbol{y} - \boldsymbol{x}$として
$$
\boldsymbol{x} + \alpha \boldsymbol{d}, \quad 0\leq \alpha \leq 1
$$を新たな解として採用する．ここで，$\alpha$はステップ・サイズであり，新たな解における目的関数値:
$$
g(\alpha) = f\left(\boldsymbol{x} + \alpha \boldsymbol{d}\right)
$$が，もとの解における目的関数値$f(\boldsymbol{x})$よりも十分に小さくなるように選ばれる．$0 \leq \alpha \leq 1$の範囲で最適なステップ・サイズを決定する方法は後述する．

Frank-Wolfe法のアルゴリズムは以下のように整理される：
- [Step 0] 初期実効可能解$\boldsymbol{x}^{(0)} \in \Omega$を選ぶ．$n:=0$とする．
- [Step 1] $\boldsymbol{x}^{(n)}$が最適解と判断できるなら終了(収束判定)．
- [Step 2] 補助問題 $\displaystyle\min_{\boldsymbol{y}}\left\{\left.(\nabla f(\boldsymbol{x}^{(n)}))^{\top} \boldsymbol{y} \right| \boldsymbol{A} \boldsymbol{y} \geq \boldsymbol{b}, \boldsymbol{y} \geq \boldsymbol{0}\right\}$を解き，補助解$\boldsymbol{y}^{(n)}$を求める． 解の改訂方向を$\boldsymbol{d}^{(n)}:= \boldsymbol{y}^{(n)}-\boldsymbol{x}^{(n)}$とする．
- [Step 3] 範囲つき直線探索問題$\displaystyle\min_{0\leq \alpha \leq 1} g(\alpha) = f\left( \boldsymbol{x}^{(n)} + \alpha \boldsymbol{d}^{(n)}\right)$ を問いてステップ・サイズ$\alpha^{(n)}$を決定する．
- [Step 4] 解を$\boldsymbol{x}^{(n+1)} := \boldsymbol{x}^{(n)} + \alpha^{(n)} \boldsymbol{d}^{(n)}$と改訂する．
- [Step 5] $n:=n+1$として [Step 1]へ．

### Frank-Wolfe 法の具体例：
$\boldsymbol{x}=(x_1, x_2)$として以下の問題を考えてみよう：
$$\begin{align}
\max_{\boldsymbol{x}} \quad f(\boldsymbol{x}) &= \beta \ln (x_{1}+1) + (1-\beta) \ln (x_{2}+1)\\
\text{s.t.} \quad & x_{1} + x_{2} \leq 8\\
& -2x_{1} + x_{2} \leq 2\\
& 2x_{1} + 3x_{2} \leq 18\\
& x_{1}, x_{2}, x_{3} \geq 0
\end{align}$$
ただし，$\beta$は$0 \leq \beta \leq 1$なる定数である． この問題は，線形計画問題で取り扱ったリア充問題の目的関数（リア充度)として [Cobb-Douglas 型関数](https://wiis.info/economics/microeconomics/consumer-theory/cobb-douglas-utility-function/)を採用したものである． **Cobb-Douglas 型関数**は，時間をバイトまたはデートの一つに「全振り/極振り」するよりも，両方の活動に「バランスよく振り分ける」方がより高いリア充度が実現できることを反映させた関数である． $\beta$はバイト時間がリア充度に対する貢献の大きさを表すパラメタで，$\beta$が大きいほどバイト時間あたりの充実度が高いことを意味する． 以下では，$\beta=0.75$としよう． このとき，リア充度$f(\boldsymbol{x})$の等高線は，下図のような原点に対して凸な関数で表される．

<div style="text-align:center; padding:20pt">
<img src="https://hackmd.io/_uploads/BJre2xd_2.png" width=300><br>
非線形リア充問題の許容領域と目的関数の等高線
</div>

許容領域は緑色で塗り潰された範囲であるため，この中で原点から最も遠い等高線が通る点が最適解となる． これを Frank-Wolfe法で探してみよう．

#### 0) 初期解 $\boldsymbol{x}^{(0)} = (0, 0)$における Frank-Wolfe法の各ステップ
初期許容解を$\boldsymbol{x}^{(0)} = (0, 0)$としよう． 目的関数の勾配
$$
\nabla f(\boldsymbol{x}) = \begin{bmatrix}
\frac{\beta}{1+x_{1}} \\ \frac{1-\beta}{1+x_{2}}
\end{bmatrix}
$$であるから，$\nabla f(\boldsymbol{x}^{0}) = (\beta, 1-\beta)$である（下図左の青い矢印）．従って，補助問題の目的関数の等高線は，下図左のように与えられる．補助解は図から$\boldsymbol{y}^{(0)} = (8,0)$と求められる． 解の改訂方向（この場合，目的関数の増加方向）を$\boldsymbol{d}^{(0)} = (8, 0)$として$\alpha$を$[0, 1]$の範囲で動かした時の目的関数は下図の右図のように$\alpha$について単調増加となる（これは，もとの$f(\boldsymbol{x})$の等高線を見ても明らか）．そこで，最適ステップ・サイズを$\alpha^{(0)} = 1$とし， 解を
$$
\boldsymbol{x}^{(1)} = \boldsymbol{y}^{(1)} = (8, 0)
$$と改訂する．

<div style="text-align:center; padding:20pt">
<img src="https://hackmd.io/_uploads/rJoTqZd_2.png" width=600><br>

$\boldsymbol{x} = (0, 0)$ における補助問題とその解（左）と，
ステップサイズに応じた目的関数値(右)．
</div>

#### 1) $\boldsymbol{x}^{(1)}=(8,0)$における各ステップ
$\boldsymbol{x}^{(1)}=(8,0)$における目的関数の勾配は$\nabla f(\boldsymbol{x}^{(1)}) = (\frac{\beta}{9}, 1-\beta)=(0.083, 0.250)$ （下図左における青い矢印）で表される．

補助問題の目的関数の等高線はこの青い矢印を法線とする平行線で表され，補助解は $\boldsymbol{y}^{(1)} = (1.5, 2)$ (下図左の橙丸)で表される． 現在の解$\boldsymbol{x}^{(1)}$と補助解$\boldsymbol{y}^{(1)}$の凸結合
$$(1-\alpha) \boldsymbol{x}^{(1)} + \alpha \boldsymbol{y}^{(1)}$$で表される解は，図中の黒い点線で表され，それらの解に対する目的関数値はステップ・サイズ$\alpha$の関数として下図右に表されている． 目的関数が最大となるステップ・サイズは$\alpha^{(1)} \approx 0.196$と求められる． これを用いて，解を
$$
\boldsymbol{x}^{(2)} = \boldsymbol{x}^{(1)} + \alpha^{(1)} \boldsymbol{d}^{(1)} = (6.725, 0.981)
$$と改訂する．改訂後の解は下図左の赤丸で表示されている．

<div style="text-align:center; padding:20pt">
<img src="https://hackmd.io/_uploads/SycphzdOh.png" width=600><br>

左図：暫定解$\boldsymbol{x}^{(1)} = (8, 0)$ (青丸)における補助問題の解（橙丸)，および最適ステップサイズにおける次の解（赤丸）<br>右図：ステップサイズに応じた目的関数値
</div>

#### 2) $\boldsymbol{x}^{(2)} = (6.725, 0.981)$ における各ステップ
$\boldsymbol{x}^{(2)}=(6.725,0.981)$における目的関数の勾配は$\nabla f(\boldsymbol{x}^{(2)}) = (\frac{\beta}{7.725}, \frac{1-\beta}{0.981}) = (0.097, 0.126)$ （下図左における青い矢印）で表される．

補助問題の目的関数の等高線はこの青い矢印を法線とする平行線で表され，補助解は$\boldsymbol{y}^{(2)} = (6, 2)$ (下図左の赤丸)で表される． このとき， 下図右のグラフより凸結合$(1-\alpha)\boldsymbol{x}^{(2)} + \alpha \boldsymbol{y}^{(2)}$における目的関数を最大とするステップ・サイズは$\alpha=1$となるため，
$$
\boldsymbol{x}^{(3)} = \boldsymbol{y}^{(2)} = (6, 2)
$$と改訂する．改訂後の解は下図左の赤丸で表示されている．
<div style="text-align:center; padding:20pt">
<img src="https://hackmd.io/_uploads/r1HdT-Yd2.png" width=600><br>

左図：暫定解$\boldsymbol{x}^{(2)} = (6.725, 0.981)$ (青丸)における補助問題の解は， 最適ステップサイズが$\alpha=1$であるため，次の解（赤丸）と重なってい表示されている．<br>右図：ステップサイズに応じた目的関数値
</div>

#### 3) $\boldsymbol{x}^{(3)} = (6, 2)$における各ステップ
$\boldsymbol{x}^{(3)}=(6,2)$における目的関数の勾配は$\nabla f(\boldsymbol{x}^{(3)}) = (\frac{\beta}{7}, \frac{1-\beta}{3})=(0.107, 0.083)$ （下図左における青い矢印）で表される．補助問題の目的関数の等高線はこの青い矢印を法線とする平行線で表され，補助解は$\boldsymbol{y}^{(3)} = (8, 0)$ (下図左の橙丸)で表される． 現在の解$\boldsymbol{x}^{(3)}$と補助解$\boldsymbol{y}^{(3)}$の凸結合
$$(1-\alpha) \boldsymbol{x}^{(3)} + \alpha \boldsymbol{y}^{(3)}$$で表される解は図中の黒い点線で表され，それらの解に対する目的関数値はステップ・サイズ$\alpha$の関数として下図右に表されている． 目的関数が最大となるステップ・サイズは$\alpha^{(3)} \approx 0.250$と求められる． これを用いて，解を
$$
\boldsymbol{x}^{(4)} = \boldsymbol{x}^{(3)} + \alpha^{(3)} \boldsymbol{d}^{(3)} = (6.5, 1.5)
$$
と改訂する． 改訂後の解は下図左の赤丸で表示されている．

<div style="text-align:center; padding:20pt">
<img src="https://hackmd.io/_uploads/B1jQkGtu2.png" width=600><br>

左図：暫定解$\boldsymbol{x}^{(3)} = (6, 2)$ (青丸)における補助問題の解（橙丸)，および最適ステップサイズにおける次の解（赤丸）<br>右図：ステップサイズに応じた目的関数値
</div>

#### 4) $\boldsymbol{x}^{(4)} = (6.5, 1.5)$ の最適性
$\boldsymbol{x}^{(4)}$における勾配は$\nabla f(\boldsymbol{x}^{4}) = (\frac{\beta}{7.5}, \frac{1-\beta}{2.5}) = (0.1, 0.1)$である． この解において束縛的(binding)な制約条件は時間制約$g_{1}(\boldsymbol{x}) = x_{1} + x_{2} - 8 \geq 0$ だけであり， その法線ベクトル$\nabla g_{1}(\boldsymbol{x}^{(4)}) = (1, 1)$ は目的関数の勾配$\nabla f(\boldsymbol{x}^{(4)})$の非負定数倍に等しい． すなわち，
$$
\nabla f(\boldsymbol{x}^{(4)}) = \lambda_{1} \nabla g_{1}(\boldsymbol{x}^{(4)})
$$が$\lambda_{1} = 0.1$において成り立つ． 以上のことから， $\boldsymbol{x}^{(4)}$ は**最適解**である．

#### Frank-Wolfe法による解軌道
以上の手続きによって得られた解$(\boldsymbol{x}^{(0)}, \boldsymbol{x}^{(1)}, \cdots, \boldsymbol{x}^{(4)})$をプロットしたものが下図である．この図において， 初期解は白い四角，最適解は赤丸で表されており，途中の暫定解は青丸で表されている． この図より， Frank-Wolfe法により，最適解へと漸近的に収束していることが判る．
<div style="text-align:center; padding:20pt">
<img src="https://hackmd.io/_uploads/HknP4GKO3.png" width=600><br>

非線形リア充問題に対して Frank-Wolfe 法を適用した場合の解軌道
</div>

### 範囲つき直線探索の方法
本節では，$0 \leq \alpha \leq 1$ の範囲内で目的関数を最小化するステップ・サイズ$\alpha$を決定する問題：
$$
\min_{0 \leq \alpha \leq 1} g(\alpha) =
f\left((1-\alpha)\boldsymbol{x}^{(n)} + \alpha \boldsymbol{y}^{(n)}\right)
$$を解く代表的な方法の1つとして**黄金分割法**(golden section method)を紹介しておこう．

黄金分割法は，単峰関数$g(\alpha)$の極小値$\alpha^{\ast}=\displaystyle\arg\min_{a\leq \alpha\leq b} g(\alpha)$を求める方法の一つで「極小値が含まれない範囲」を除外していくことで，極小値が存在する範囲を逐次的に狭めていく方法である．
いま，$n$回目の繰り返しにおける極小値が含まれる区間を$\left[a^{(n)}, b^{(n)}\right]$としよう．この区間内に，2つの**内分点**$l^{(n)}<r^{(n)}$を設け，それぞれにおける目的関数値を比較する．

目的関数が単峰であれば, $g(l^{(n)})$と$g(r^{(n)})$の大小関係から，以下のいずれかが判明する：
1. $g(l^{(n)}) \leq g(r^{(n)})$ならば，$r^{(n)}$から$b^{(n)}$の間に極小値は含まれない．
2. $g(l^{(n)}) > g(r^{(n)})$ならば，$a^{(n)}$から$l^{(n)}$の間に極小値は含まれない．

そこで，$g(l^{(n)}) \leq g(r^{(n)})$の場合は区間$[r^{(n)}, b^{(n)}]$を候補から除外し，新しい区間を$[a^{(n+1)}, b^{(n+1)}] := [a^{(n)},r^{(n)}]$とする．その上で，**内分点**について$r^{(n+1)}:=l^{(n)}$とし，$l^{(n+1)}$を新たに選ぶ（下図左参照）．

一方，$vg(l^{(n)}) > g(r^{(n)})$の場合は区間$[a^{(n)}, l^{(n)}]$を候補から除外し，新しい区間を$[a^{(n+1)}, b^{(n+1)}]:=[l^{(n)}, b^{(n)}]$とする．その上で，**内分点**について$l^{(n+1)}:=r^{(n)}$とし，$r^{(n+1)}$を新たに選ぶ（下図右参照）．

<table><tr><td>


</td><td>


</td></tr></table>

上述の手続きの特徴は，各ステップにおいて新たに目的関数を評価すべき点が$l^{(n+1)}$もしくは$r^{(n+1)}$の**いずれか一方のみ**で良い点である．これは，各ステップにおける**内分点** $l^{(n)}, r^{(n)}$を**黄金比**:
$$
R = \frac{\sqrt{5}-1}{2} \approx 0.618
$$
に従って選ぶことで実現できる．具体的には，$l^{(n)}$および$r^{(n)}$を，候補区間$[a^{(n)}, b^{(n)}]$を$[1-R,R]$および$[R, 1-R]$に分割する点として選ぶ．これにより「$n$回目に除外された区間が元の候補区間$[a^{(n)}, b^{(n)}]$に占める比」と「$n+1$回目に新たに求められる内分点が新しい候補区間$[a^{(n+1), b^{(n+1)}}]$に占める比」が一致する．

なお，黄金分割法においては$n$ステップ目の候補区間の長さを$I^{(n)}=b^{(n)}-a^{(n)}$とするとき，以下の関係が成り立つ：
$$
I^{(n)} = I^{(n+1)} + I^{(n+2)}
$$

黄金分割法のアルゴリズムは以下のように整理される：
- [Step 0] 初期候補区間$[a^{(0)}, b^{(0)}]$を選び，その長さを$I^{(0)}=b^{(0)}-a^{(0)}$とする．黄金比$R=0.618$を用いて，内分点を$(l^{(n)}, r^{(n)}) = (a^{(0)}+(1-R)I^{(0)}, a^{(0)}+RI^{(0)})$とする．$n:=0$とする．
- [Step 1] 内分点における目的関数値$g(l^{(n)}), g(r^{(n)})$の大小関係を比較する：
- $g(l^{(n)})\leq g(r^{(n)})$の場合，候補区間および内分点を以下のように改訂する：
$$\begin{align}(a^{(n+1)}, r^{(n+1)}, b^{(n+1)}) &:= (a^{(n)}, l^{(n)}, r^{(n)})\\l^{(n+1)}&:= a^{(n)}+(1-R)I^{(n)}\end{align}$$
- $g(l^{(n)}) > g(r^{(n)})$の場合，候補区間および内分点を以下のように改訂する：
$$\begin{align}(a^{(n+1)}, l^{(n+1)}, b^{(n+1)}) &:= (l^{(n)}, r^{(n)}, b^{(n)})\\r^{(n+1)}&:=a^{(n)}+RI^{(n)}\end{align}$$
- [Step 2] $I^{(n+1)}:=b^{(n+1)}-a^{(n)}$とする．$I^{(n+1)}$が十分に小さければ，候補区間の中点$\alpha^{\ast}:= \frac{b^{(n+1)}+a^{(n+1)}}{2}$を最適解として終了．そうでなければ，$n:=n+1$として[Step 1]に戻る．

## 逐次二次計画法
上述の Frank-Wolfe法は，線形制約つき最適化問題を**線形計画問題として近似**する手法であった．しかし，制約条件が線形で無い場合には適用できず，収束もそれほど早くない（最適解の近傍で"ジグザグ運動"することが知られている）．逐次二次計画法は，最適化問題を**線形制約つき二次計画問題として近似**する手法である．

### 等式制約条件のみの問題
まず，等式制約条件のみの以下の問題を考える：
$$\begin{align}
\min_{\boldsymbol{x}} &f(\boldsymbol{x})\\
\text{s.t.}{} & g_{j}(\boldsymbol{x}) = 0, \quad j = 1, \cdots, M
\end{align}$$

この問題の Lagrangian は，Lagrange乗数を$\boldsymbol{\lambda} = (\lambda_{j} : j = 1, \cdots, M)$とすれば，
$$
L(\boldsymbol{x}, \boldsymbol{\lambda}) = f(\boldsymbol{x}) - \sum_{j=1}^{M} \lambda_{j} g_{j}(\boldsymbol{x})
$$
で表されるから，最適性条件(KKT条件)は，以下の連立方程式で表される：
$$\begin{align}
\nabla L(\boldsymbol{x}, \boldsymbol{\lambda}) &= \begin{bmatrix}
\nabla_{\boldsymbol{x}} L(\boldsymbol{x}, \boldsymbol{\lambda}) \\
\nabla_{\boldsymbol{\lambda}} L(\boldsymbol{x}, \boldsymbol{\lambda})
\end{bmatrix}\\
&= \begin{bmatrix}
\nabla f(\boldsymbol{x}) - (\nabla \boldsymbol{g}(\boldsymbol{x}))^{\top} \boldsymbol{\lambda}\\
\boldsymbol{g}(\boldsymbol{x})
\end{bmatrix} = \boldsymbol{0}
\end{align}$$

つまり，等式制約条件のみの最適化問題は，明示的な未知変数を$\boldsymbol{x}, \boldsymbol{\lambda}$とし，Lagrangian $L(\boldsymbol{x}, \boldsymbol{\lambda})$を目的関数とする，制約なしの最適化問題に帰着する．そこで，その解法に Newton 法を用いることを考えよう．

$k$回目における暫定解を$(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)})$とすると，Newton方向は
$$
\nabla L(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)}) + \nabla^{2} L(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)}) \begin{bmatrix}
\boldsymbol{x} - \boldsymbol{x}^{(n)} \\
\boldsymbol{\lambda} - \boldsymbol{\lambda}^{(n)}
\end{bmatrix} = \boldsymbol{0}
$$
の解として得られる．ここで，Lagrangianの Hessian は，
$$\begin{align}
\nabla^{2} L(\boldsymbol{x}, \boldsymbol{\lambda}) &= \begin{bmatrix}
\nabla_{\boldsymbol{x}}^{2} L(\boldsymbol{x}, \boldsymbol{\lambda}) & -(\nabla \boldsymbol{g}(\boldsymbol{x}))^{\top}\\
-\nabla \boldsymbol{g}(\boldsymbol{x}) & \boldsymbol{0}
\end{bmatrix}
\end{align}$$
で表される．ただし，
$$
\nabla_{\boldsymbol{x}}^{2} L(\boldsymbol{x}, \boldsymbol{\lambda}) = \nabla^{2} f(\boldsymbol{x}) - \sum_{j=1}^{M} \lambda_{j} \nabla^{2} g_{j}(\boldsymbol{x}).
$$
従って，Newton方向が従う方程式は
$$
\begin{align}
&\nabla L^{(n)}(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)}) + \nabla^{2} L(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)}) \begin{bmatrix}
\boldsymbol{x} - \boldsymbol{x}^{(n)} \\
\boldsymbol{\lambda} - \boldsymbol{\lambda}^{(n)}
\end{bmatrix}\\
&= \begin{bmatrix}
\nabla f^{(n)} + \nabla_{\boldsymbol{x}}^{2} L^{(n)} - (\nabla \boldsymbol{g}^{(n)})^{\top}\boldsymbol{\lambda} \\
\boldsymbol{g}^{(n)} - (\nabla \boldsymbol{g}^{(n)})(\boldsymbol{x}-\boldsymbol{x}^{(n)})
\end{bmatrix} = \boldsymbol{0}
\end{align}
$$
となる．ただし，記述の簡便化のために$$\begin{align}
\nabla f^{(n)} &= \nabla f(\boldsymbol{x}^{(n)}),&
\nabla_{\boldsymbol{x}}^{2} L^{(n)} &= \nabla_{\boldsymbol{x}}^{2} L(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)})\\
\boldsymbol{g}^{(n)} &= \begin{bmatrix}
g_{1}(\boldsymbol{x}^{(n)}) \\ \vdots \\ g_{M}(\boldsymbol{x}^{(n)})\end{bmatrix},&
\nabla \boldsymbol{g}^{(n)} &= \begin{bmatrix}(\nabla g_{1}(\boldsymbol{x}^{(n)}))^{\top}\\\vdots\\(\nabla g_{M}(\boldsymbol{x}^{(n)}))^{\top}\end{bmatrix}
\end{align}$$
とした．

等式制約のみであれば，この方程式を直接解くことでNewton方向$(\boldsymbol{x}-\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}-\boldsymbol{\lambda}^{(n)})$が求められるが，ここでは不等式制約つき問題への拡張も考えて，次の二次計画問題を解くことを考える．
$$\begin{align}\text{[QP-Aux]}\quad
\min_{\boldsymbol{d}, \boldsymbol{v}} & \frac{1}{2} \boldsymbol{d}^{\top} \nabla_{\boldsymbol{x}}^{2} L^{(n)} \boldsymbol{d} + (\nabla f^{(n)})^{\top} \boldsymbol{d} \\
\text{s.t.}{} & \boldsymbol{g}^{(n)} - (\nabla \boldsymbol{g}^{(n)}) \boldsymbol{d} = 0, \quad j = 1, \cdots, M
\end{align}$$
この問題の最適性条件は，
$$\begin{align}
&\nabla_{\boldsymbol{x}}^{2} L^{(n)} \boldsymbol{d} + \nabla f^{(n)}
-(\nabla \boldsymbol{g}^{(n)})^{\top} \boldsymbol{\lambda} = \boldsymbol{0}\\
&\boldsymbol{g}^{(n)} - (\nabla \boldsymbol{g}^{(n)}) \boldsymbol{d} = \boldsymbol{0}
\end{align}
$$
と表せる．これは Newton 方向が従う方程式に他ならない．従って，二次計画問題$\text{[QP-Aux]}$を解けば Newton 方向 $(\boldsymbol{d}, \boldsymbol{\lambda}-\boldsymbol{\lambda}^{(n)})$ が求められる．

補助問題$\text{[QP-Aux]}$を解いて$\boldsymbol{d}=\boldsymbol{0}$が得られる場合，補助問題の最適性条件は
$$\begin{align}
&\nabla f^{(n)}
-(\nabla \boldsymbol{g}^{(n)})^{\top} \boldsymbol{\lambda} = \boldsymbol{0}\\
&\boldsymbol{g}^{(n)} = \boldsymbol{0}
\end{align}
$$
となる．これは元の問題の最適性条件に他ならない．従って，探索方向が$\boldsymbol{d}=\boldsymbol{0}$となった時，$(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)})$は元の問題の最適解である．

なお，$\text{[QP-Aux]}$の最適性条件は，Lagrangian を
$$\begin{align}
J(\boldsymbol{d}, \boldsymbol{\lambda})
&= \frac{1}{2} \boldsymbol{d}^{\top} \nabla_{\boldsymbol{x}}^{2} L^{(n)} \boldsymbol{d} + \boldsymbol{d}^{\top}(\nabla f^{(n)}) - \boldsymbol{\lambda}^{\top}\left\{\boldsymbol{g}^{(n)} - (\nabla \boldsymbol{g}^{(n)}) \boldsymbol{d}\right\}
\end{align}$$
と定義すれば，
$$\begin{align}
\nabla J(\boldsymbol{d}, \boldsymbol{\lambda}) =
\begin{bmatrix}
\nabla_{\boldsymbol{d}} J(\boldsymbol{d}, \boldsymbol{\lambda})\\
\nabla_{\boldsymbol{\lambda}} J(\boldsymbol{d}, \boldsymbol{\lambda})
\end{bmatrix}
= \begin{bmatrix}
\nabla f^{(n)} + \nabla_{\boldsymbol{x}}^{2} L^{(n)} \boldsymbol{d} - (\nabla \boldsymbol{g}^{(n)})^{\top} \boldsymbol{\lambda}\\
\boldsymbol{g}^{(n)} - (\nabla \boldsymbol{g}^{(n)}) \boldsymbol{d}
\end{bmatrix}
\end{align}$$
であることから容易に求められる．

### 等式制約と不等式制約が混在する場合
等式制約と不等式制約の両方が混在する以下の最適化問題を考える：
$$\begin{align}
\min_{\boldsymbol{x}} &f(\boldsymbol{x})\\
\text{s.t.}{} & g_{j}(\boldsymbol{x}) \geq 0, \quad j = 1, \cdots, M\quad (\lambda_{j})\\
& h_{k}(\boldsymbol{x}) = 0, \quad k = 1, \cdots, L\quad (\mu_{k})
\end{align}$$
ただし，$\boldsymbol{\lambda}=(\lambda_{j}: j = 1, \cdots, M)$および$\boldsymbol{\mu}=(\mu_{k}: k = 1, \cdots, L)$は，それぞれ，不等式制約および等式制約に関する Lagrange 乗数である．この問題の最適性条件は，
$$\begin{align}
&\nabla f(\boldsymbol{x}) - (\nabla \boldsymbol{g}(\boldsymbol{x}))^{\top}\boldsymbol{\lambda} - (\nabla \boldsymbol{h}(\boldsymbol{x}))^{\top} \boldsymbol{\mu} = \boldsymbol{0}\\
&\begin{cases}
\boldsymbol{\lambda}^{\top} \boldsymbol{g}(\boldsymbol{x}) = 0\\
\boldsymbol{\lambda} \geq \boldsymbol{0}, \, \boldsymbol{g}(\boldsymbol{x}) \geq \boldsymbol{0}
\end{cases}\\
& \boldsymbol{h}(\boldsymbol{x}) = \boldsymbol{0}
\end{align}$$
で表される．

$k$回目における暫定解を$(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)}, \boldsymbol{\mu}^{(n)})$とした時，Newton方向を求める補助問題を以下のように定式化する：
$$\begin{align}\text{[QP-Aux]}\quad
\min_{\boldsymbol{d}, \boldsymbol{u}, \boldsymbol{v}} & \frac{1}{2} \boldsymbol{d}^{\top} \nabla_{\boldsymbol{x}}^{2} L^{(n)}\boldsymbol{d} + (\nabla f^{(n)})^{\top} \boldsymbol{d} \\
\text{s.t.}{} & g_{j}^{(n)} - (\nabla g_{j}^{(n)})^{\top} \boldsymbol{d} \geq 0, \quad j = 1, \cdots, M\\
& h_{k}^{(n)} - (\nabla h_{k}^{(n)})^{\top} \boldsymbol{d} = 0,
\quad k = 1, \cdots, L
\end{align}$$
この補助問題の最適性条件は，
$$\begin{align}
&\nabla_{\boldsymbol{x}}^{2} L^{(n)}\boldsymbol{d} + (\nabla f^{(n)})- (\nabla \boldsymbol{g}^{(n)})^{\top} \boldsymbol{\lambda} - (\nabla \boldsymbol{h}^{(n)})^{\top} \boldsymbol{\mu} = \boldsymbol{0}\\
&\begin{cases}\boldsymbol{\lambda}^{\top}
\left\{\boldsymbol{g}^{(n)} - (\nabla \boldsymbol{g}^{(n)})\boldsymbol{d}\right\} = 0\\
\boldsymbol{\lambda} \geq \boldsymbol{0}, \boldsymbol{g}^{(n)} - (\nabla \boldsymbol{g}^{(n)})\boldsymbol{d} \geq \boldsymbol{0}\\
\end{cases}\\
&\boldsymbol{h}^{(n)} - (\nabla \boldsymbol{h}^{(n)}) \boldsymbol{d} = \boldsymbol{0}
\end{align}$$
であり，やはり$\boldsymbol{d}=\boldsymbol{0}$であるときに$(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}, \boldsymbol{\mu})$は元の問題の最適性条件を満たす．

上記の$\text{[QP-Aux]}$は，$\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)}$が改訂されるたびに Lagrangian の Hessian
$$
\nabla_{\boldsymbol{x}}^{2} L^{(n)} = \nabla_{\boldsymbol{x}}^{2} L(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n)}) = \nabla^{2} f(\boldsymbol{x}^{(n)}) + \sum_{j=1}^{M}\lambda_{j}^{(n)} \nabla^{2} g_{j}(\boldsymbol{x}^{(n)})
$$
を計算し直す必要があり，大規模な問題に対しては容易ではない．また，$\nabla_{\boldsymbol{x}}^{2}L^{(n)}$が正定値でない場合には$\text{[QP-Aux]}$の解を求めるのが困難になる．これらは Newton 法の問題点と共通である．そこで，準Newton法と同様の考え方を用いて，$L^{(n)}$を正定値行列$\boldsymbol{B}^{(n)}$で近似した下記の問題を解くことを考える．
$$\begin{align}\text{[QP-Aux]}\quad
\min_{\boldsymbol{d}, \boldsymbol{u}, \boldsymbol{v}} & \frac{1}{2} \boldsymbol{d}^{\top} \boldsymbol{B}^{(n)}\boldsymbol{d} + (\nabla f^{(n)})^{\top} \boldsymbol{d} \\
\text{s.t.}{} & g_{j}^{(n)} - (\nabla g_{j}^{(n)})^{\top} \boldsymbol{d} \geq 0, \quad j = 1, \cdots, M\\
& h_{k}^{(n)} - (\nabla h_{k}^{(n)})^{\top} \boldsymbol{d} = 0,
\quad k = 1, \cdots, L
\end{align}$$

修正 Newton 法で用いた BFGS公式をそのまま適用した場合，状態変数とLagrangianの勾配の差分
$$\begin{align}
\boldsymbol{s}^{(n)} &=\boldsymbol{x}^{(n+1)}-\boldsymbol{x}^{(n)}\\
\boldsymbol{y}^{(n)} &= \nabla L_{\boldsymbol{x}}(\boldsymbol{x}^{(n+1)}, \boldsymbol{\lambda}^{(n+1)}, \boldsymbol{\mu}^{(n+1)}) - \nabla L_{\boldsymbol{x}}(\boldsymbol{x}^{(n)}, \boldsymbol{\lambda}^{(n+1)}, \boldsymbol{\mu}^{(n+1)})
\end{align}$$
およびそれを用いた係数：
$$\begin{align}
\beta^{(n)} &= \boldsymbol{s}^{(n)}, &
\gamma^{(n)} &= \boldsymbol{B}^{(n)} \boldsymbol{s}^{(n)}
\end{align}$$
を用いて正定値行列$\boldsymbol{B}^{(n)}$の改訂式は
$$
\boldsymbol{B}^{(n+1)} := \boldsymbol{B}^{(n)} + \frac{1}{\beta^{(n)}} \boldsymbol{y}^{(n)} - \frac{1}{\gamma^{(n)}}
\boldsymbol{B}^{(n)} \boldsymbol{s}^{(n)} \boldsymbol{B}^{(n)}
$$
と書ける．ただし，逐次二次計画法では$\beta^{(n)} > 0$とは限らないため，上記公式をそのまま用いるだけでは，$\boldsymbol{B}^{(n+1)}$が**正定値となることを保証できない**．

そこで，$\beta^{(n)}$が小さい時にだけ以下の修正を行なう**Powellの修正BFGS公式**が知られている：まず，パラメタ$0 < \xi < 1$を決めておき，$\beta^{(n)} < \xi \gamma^{(n)}$の時には
$$\begin{align}
\phi^{(n)} &= \frac{(1-\xi) \gamma^{(n)}}
{\gamma^{(n)}-\beta^{(n)}}, &
\hat{\boldsymbol{y}}^{(n)} &= \phi^{(n)} \boldsymbol{y}^{(n)} + (1-\phi^{(n)}) \boldsymbol{B}^{(n)} \boldsymbol{s}^{(n)}
\end{align}$$
とし，この$\hat{\boldsymbol{y}}^{(n)}$を用いて$\hat{\beta}^{(n)} = (\hat{\boldsymbol{y}}^{(n)})^{\top}\boldsymbol{s}^{(n)}$として以下の式で正定値行列$\boldsymbol{B}^{(n)}$を改訂する：
$$
\boldsymbol{B}^{(n+1)} := \boldsymbol{B}^{(n)} + \frac{1}{\hat{\beta}^{(n)}} \boldsymbol{y}^{(n)} - \frac{1}{\gamma^{(n)}}
\boldsymbol{B}^{(n)} \boldsymbol{s}^{(n)} \boldsymbol{B}^{(n)}
$$
この更新公式を用いると$\boldsymbol{B}^{(n+1)}$は正定値となることが保証できる．

降下方向$\boldsymbol{d}^{(n)}$へ進む量(ステップ・サイズ)$\alpha>0$を決定する際には，制約条件を満たさなければならないため，目的関数そのものではなく，制約条件に応じたペナルティ関数やバリア関数を組み込んだ**メリット関数** (後述)$\phi(\boldsymbol{x})$を最小化するように一次探索を行う：
$$
\alpha^{(n)} := \arg\min_{\alpha \geq 0} \phi(\boldsymbol{x}^{(n)} + \alpha \boldsymbol{d}^{(n)})
$$

逐次二次計画法のアルゴリズムは以下のように整理される：
- [Step 0] 初期実行可能解$\boldsymbol{x}^{(0)} \in \Omega$および初期正定値行列$\boldsymbol{B}^{(0)}$を選ぶ．パラメタ$0 < \xi < 1$を選ぶ．$n:=0$とする．
- [Step 1] $\boldsymbol{x}^{(n)}$が最適解と判断できるなら終了(収束判定)．
- [Step 2] 補助問題[QP-Aux]を解き，降下方向$\boldsymbol{d}^{(n)}$を求める．
- [Step 3] メリット関数$\phi$に対する直線探索問題$\displaystyle\min_{0\leq \alpha } \phi(\boldsymbol{x}^{(n)} + \alpha \boldsymbol{d}^{(n)})$ を問いてステップ・サイズ$\alpha^{(n)}$を決定する．
- [Step 4] 解を$\boldsymbol{x}^{(n+1)} := \boldsymbol{x}^{(n)} + \alpha^{(n)} \boldsymbol{d}^{(n)}$と改訂する．
- [Step 5] $\beta^{(n)} = \boldsymbol{s}^{(n)}$と$\gamma^{(n)} = \boldsymbol{B}^{(n)} \boldsymbol{s}^{(n)}$を求め，$\beta^{(n)} < \xi\gamma^{(n)}$か否かに応じて$\boldsymbol{B}^{(n+1)}$を決定(Powellの修正BFGS公式)．
- [Step 5] $n:=n+1$として [Step 1]へ．

### メリット関数
不等式制約つき非線形最適化問題：
$$
\min_{\boldsymbol{x}} \left\{f(\boldsymbol{x})\left|g_{j}(\boldsymbol{x}) \geq 0, \quad j = 1, \cdots, M\right.\right\}
$$
を制約なし最適化問題に変換して解くことを試みる．その代表的方法として，**ペナルティ関数法**および**バリア関数法**が知られている．

**ペナルティ関数**とは，状態変数$\boldsymbol{x}$が許容領域$\Omega=\left\{\boldsymbol{x}| g_{j}(\boldsymbol{x}) \geq 0, \quad j = 1, \cdots, M\right\}$内にあるときは0，そうでない時に正の値を取るような関数である．例えば，
$$
P(\boldsymbol{x}) = - \sum_{j=1}^{M} \min\left\{g_{j}(\boldsymbol{x}), 0\right\}
$$
などが挙げられる．

**バリア関数**とは，許容領域$\Omega$の内部で定義され，$\Omega$の境界に近づくと無限大に発散するような関数である．例えば，
$$
B(\boldsymbol{x}) = \sum_{j=1}^{M} 1/g_{j}(\boldsymbol{x})
$$
が挙げられる．なお，非負制約$x_{i}\geq0, \quad i =1, \cdots, N$に対しては，以下のバリア関数がよく利用される：
$$
B(\boldsymbol{x}) = - \sum_{i=1}^{N} \ln x_{i}
$$

**メリット関数**は目的関数にペナルティ関数やバリア関数を組み込んだもので，最もシンプルなものは，適当なパラメタを$\mu > 0$として
$$\begin{align}
\phi(\boldsymbol{x}) &= f(\boldsymbol{x}) + \mu P(\boldsymbol{x}),&
\phi(\boldsymbol{x}) &= f(\boldsymbol{x}) + \mu B(\boldsymbol{x})
\end{align}$$
などとしたものである．