非線形な場合の一般論は扱わず，目的関数が2次形式，制約はアフィンなものに絞って主内点法，主双対内点法のアルゴリズムをメモる．
一般的な定式化
以降基本的に $\boldsymbol{x} \in \mathbb{R}^{n}$ とする．等式制約と不等式制約の下で目的関数を最小化する問題は一般的に以下のように表される．
$$\begin{align*} \min \ & f(\boldsymbol{x}) \\ {\rm s.t.}\ \ & \boldsymbol{g}(\boldsymbol{x}) = 0 \in \mathbb{R}^{m} \\ & \boldsymbol{h}(\boldsymbol{x}) \geq 0 \in \mathbb{R}^{l} \end{align*}$$
ここで不等式制約のスラック変数 $\boldsymbol{s} \in \mathbb{R}^{l}$ を導入すると
$$\begin{align*} \min \ & f(\boldsymbol{x}) \\ {\rm s.t.}\ \ & \boldsymbol{g}(\boldsymbol{x}) = 0 \in \mathbb{R}^{m} \\ & \boldsymbol{h}(\boldsymbol{x}) - \boldsymbol{s} = 0 \in \mathbb{R}^{l} \\ & \boldsymbol{s} \geq 0 \in \mathbb{R}^{l} \end{align*}$$
のように不等式制約を等式制約にできる．
標準形への変換
さらに
- $\boldsymbol{x} = \boldsymbol{x}_{+} - \boldsymbol{x}_{-}, \boldsymbol{x}_{\pm} \geq 0$
- $\tilde{\boldsymbol{x}}$ を $(\boldsymbol{x}_{+}, \boldsymbol{x}_{-}, \boldsymbol{s})$ を縦に並べたベクトル
- $\tilde{\boldsymbol{g}}(\boldsymbol{\boldsymbol{x}})$ を $(\boldsymbol{g}(\boldsymbol{x}), \boldsymbol{h}(\boldsymbol{x}) - \boldsymbol{s})$ を縦に並べたベクトル
とすると
$$\begin{align*} \min \ & f(\tilde{\boldsymbol{x}}) \\ {\rm s.t.}\ \ & \boldsymbol{g}(\tilde{\boldsymbol{x}}) = 0 \in \mathbb{R}^{m+l} \\ & \tilde{\boldsymbol{x}} \geq 0 \in \mathbb{R}^{n+n+l} \end{align*}$$
のように変形できる．
主内点法
バリア関数を用いる手法．上の標準形において $\boldsymbol{x} \geq 0$ を $\log$ を使って以下のように置き換える．
$$\begin{align*} \min \ & f(\boldsymbol{x}) - \nu \sum \log(x_i) \\ {\rm s.t.}\ \ & \boldsymbol{g}(\boldsymbol{x}) = 0 \in \mathbb{R}^{m} \end{align*}$$
もし制約 $\boldsymbol{x} \geq 0$ が破られそうになると $-\log$ の項が無限大に近づくため目的関数が大きくなってしまう．そのため制約を破らないよう目的関数を最小化できる．
こちらの講義資料 によるとこのバリア関数法により得られた最適解 $\boldsymbol{x}^{\star}$ による最小値と真の最小値 $f^{\star}$ の差は以下のように評価される．
$$\begin{align*} f(\boldsymbol{x}^{\star}) - f^{\star} \leq \nu \mathrm{dim}(\boldsymbol{x}) \end{align*}$$
よってある $\nu^k$ の値を使ってニュートン法などで最適値を求めてみた後，もし $\nu^k \mathrm{dim}(\boldsymbol{x})$ が $\epsilon$ 以下であればその解を採用すればよい．もし満たさなければ $\nu^{n+1} \leftarrow c\nu^{n}$ として $\nu$ を小さくし，再度最適化問題を解き直す．
2次計画法の場合
以下の問題を数値的に解いてみよう．
$$\begin{align*} \min \ & \frac{1}{2}\boldsymbol{x}^{\text{T}}Q\boldsymbol{x} + \boldsymbol{c}^{\text{T}}\boldsymbol{x} - \nu \sum \log(x_i) \\ {\rm s.t.}\ \ & A\boldsymbol{x} = \boldsymbol{b} \in \mathbb{R}^{m} \end{align*}$$
等式制約に対するラグランジュ乗数を $\boldsymbol{\lambda} \in \mathbb{R}^m$ とするとラグランジアンは
$$\begin{align*} L(\boldsymbol{x}, \boldsymbol{\lambda}) = \frac{1}{2}\boldsymbol{x}^{\text{T}}Q\boldsymbol{x} + \boldsymbol{c}^{\text{T}}\boldsymbol{x} - \nu \sum \log(x_i) + \boldsymbol{\lambda}^{\text{T}}(A\boldsymbol{x} + \boldsymbol{b}) \end{align*}$$
で与えられる．未定乗数法により
$$\begin{align*} Q \boldsymbol{x} + \boldsymbol{c} - \nu [1/x_i] + A^{\text{T}} \boldsymbol{\lambda} = 0 \end{align*}$$
が最適解を与える必要条件である．勾配法によりこれを満たす $(\boldsymbol{x}^{\star}, \boldsymbol{\lambda}^{\star})$ を求めたいのであるが， $1 / x_i$ を微分するのは都合が悪い．そこで $z_i = \nu / x_i$ としてベクトル $\boldsymbol{z} \in \mathbb{R}^n$ を定義する．すると元の等式制約を含めて，以下を満たす $(\boldsymbol{x}, \boldsymbol{\lambda}, \boldsymbol{z})$ を求める問題に変形される．
ただし $X = {\rm diag}(x_1, \ldots, x_n)$ ， $Z = {\rm diag}(z_1, \ldots, z_n)$ である．数値解法により各変数の勾配を求めて
のようにして更新を行うことを考えると， $(\boldsymbol{x}_{n+1}, \boldsymbol{\lambda}_{n+1}, \boldsymbol{z}_{n+1})$ を上の方程式に代入して以下の連立方程式が得られる．
ここで $dX_n$ は $d\boldsymbol{x}_n$ を対角成分とする行列である． また $dX_n dZ_n \approx 0$ とするとこれは以下のような連立方程式に書き換えられる．
初期値 $(\boldsymbol{x}_{0}, \boldsymbol{\lambda}_{0}, \boldsymbol{z}_{0})$ から左辺の行列と右辺のベクトルは具体的な数値が計算できるので，各勾配も数値的に求めることができる．さらに $\Sigma_n = X_{n}^{-1}Z_n$ とすると
$$\begin{align*} Z_n d\boldsymbol{x}_{n} + X_n d\boldsymbol{z}_{n} = \nu \boldsymbol{e} - X_n Z_n e \end{align*}$$
の両辺に $X_{n}^{-1}$ をかけることで
$$\begin{align*} \Sigma_n d\boldsymbol{x}_{n} + d\boldsymbol{z}_{n} = \nu X_{n}^{-1} \boldsymbol{e}-\boldsymbol{z}_n \end{align*}$$
となるので， $\nu X_{n}^{-1} \boldsymbol{e}$ の部分を無視して連立方程式に再代入すると以下の連立方程式が得られる．
$$\begin{align*} \begin{bmatrix} A & O \\ Q + \Sigma_n& A^{\text{T}} \\ \end{bmatrix} \begin{bmatrix} d\boldsymbol{x}_{n} \\ d\boldsymbol{\lambda}_{n} \end{bmatrix} = -\begin{bmatrix} A\boldsymbol{x}_{n} - \boldsymbol{b} \\ Q\boldsymbol{x}_{n} + A^{\text{T}} \boldsymbol{\lambda}_{n} \end{bmatrix} \end{align*}$$
この連立方程式を解いて得られた解 $(d\boldsymbol{x}_{n}, d\boldsymbol{\lambda}_{n})$ を用いて
$$\begin{align*} d\boldsymbol{z}_{n} = \nu X_{n}^{-1} \boldsymbol{e} - \boldsymbol{z}_n - \Sigma_n d\boldsymbol{x}_{n} \end{align*}$$
を求める．こうすることで求まった各変数の勾配に刻み幅を乗じて解を更新し，$(\boldsymbol{x}^{\star}, \boldsymbol{\lambda}^{\star}, \boldsymbol{z}^{\star})$ が収束したら終了すればよい．
双対問題
上の方法では不等式制約 $\boldsymbol{x} \geq 0$ は陽には用いられていなかった．ここでは不等式制約に対してもラグランジュ乗数をかけてラグランジアンを作り，KKT条件を導く．
まずはじめに以下のようにラグランジアンを作る．
$$\begin{align*} L(\boldsymbol{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}) = f(\boldsymbol{x}) + \boldsymbol{\lambda}^{\text{T}} \boldsymbol{g}(\boldsymbol{x}) + \boldsymbol{\mu}^{\text{T}} \boldsymbol{h}(\boldsymbol{x}) \end{align*}$$
弱双対定理
このラグランジアンを $\boldsymbol{x}$ が等式・不等式制約を満たすような範囲で最小化した以下の関数
$$\begin{align*} q(\boldsymbol{\lambda}, \boldsymbol{\mu}) = \inf_{x} L(\boldsymbol{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}) \end{align*}$$
をラグランジュ双対関数という．$\boldsymbol{\mu} \geq 0$ の下で
$$\begin{align*} f(\boldsymbol{x}) \geq q(\boldsymbol{\lambda}, \boldsymbol{\mu}) \end{align*}$$
が成立するというのが弱双対性定理．またこの最大化問題は双対問題と呼ばれている．
強双対定理
もともとの最小化問題の最小値 $f(\boldsymbol{x}^{\star})$ と双対問題の最大値 $q(\boldsymbol{\lambda}^{\star}, \boldsymbol{\mu}^{\star})$ が一致するというのが強双対性定理．
KKT条件
以上の定理から以下のKKT条件(+元々の制約)が導かれる．
$$\begin{align*} & \nabla f(x) + \sum_{i=1}^m \lambda_i \nabla g_i(x) + \sum_{i=1}^l \mu_i\ \nabla h_i(x)=0 \\ & \mu_i \geq 0 \\ & \mu_i h_i(\boldsymbol{x}) = 0 \end{align*}$$
最後の式は相補性条件と呼ばれている．
主双対内点法(パス追従法)
先ほどと同じ2次計画問題
$$\begin{align*} \min \ & \frac{1}{2}\boldsymbol{x}^{\text{T}}Q\boldsymbol{x} + \boldsymbol{c}^{\text{T}} \boldsymbol{x} \\ {\rm s.t.}\ \ & A\boldsymbol{x} = \boldsymbol{b} \in \mathbb{R}^{m} \\ & x_i \geq 0 \end{align*}$$
にKKT条件を適用すると
$$\begin{align*} A \boldsymbol{x} &= b \\ A^{\text{T}} \boldsymbol{\lambda} + \boldsymbol{\mu} &= Q \boldsymbol{x} + \boldsymbol{c} \\ \boldsymbol{\mu} &\geq 0 \\ x_i \mu_i &= 0 \end{align*}$$
が得られる．ここで $X = \mathrm{diag} (x_1, \cdots, x_n)$ とすると最後の式は
$$\begin{align*} X \boldsymbol{\mu} = 0 \end{align*}$$
と表される．先ほどと同じように現在の解に差分を足して数値的に解を更新していくので，
を代入して差分に対する方程式を作っていく．
第3式において $dX_n d \boldsymbol{\mu}_n \sim 0$ とする．ここで現在の解の値を用いて
$$\begin{align*} \nu_n = \sum_i (x_i \mu_i)_n \end{align*}$$
とし，適当な係数 $\sigma \in [0, 1]$ を用いて第3式を
$$\begin{align*} X_n \boldsymbol{\mu}_n+ X_n d \boldsymbol{\mu}_n + (dX_n) \boldsymbol{\mu}_n = \sigma \nu_n \boldsymbol{e} \end{align*}$$
として以下の連立方程式により差分を求める．
これを用いて解を更新する際に，不等式制約 $x_i \geq 0$，$\mu_i \geq 0$ を破らないように
としてこれを差分にかけて解を更新すれば良い．
バリア関数法との比較
バリア関数法は以下のように2重ループである．
while converge:
while converge:
推定解を初期化
バリア関数法で解を更新
収束したらbreak
end
$\nu \mathrm{dim}(\boldsymbol{x}) < \epsilon$ ならbreak
$\nu \leftarrow c\nu$
end
一方パス追従法は1重ループである．
while converge:
推定解を初期化
パス追従法で解を更新
相補性ギャップ $\nu$ が $\epsilon$ 以下ならbreak
end