非線形計画問題の主双対内点法

はじめに

主双対内点法とは、実行可能領域の内部を最適解に向けて探索する手法である。本記事では、非線形問題に対する主双対内点法のアルゴリズムについて述べる。線形問題に対する主双対内点法については以下の記事を参照のこと。
線形計画問題の主双対内点法

対象とする非線形問題

以下の制約付き非線形計画問題を考える。

$$ \begin{array}{l} \mathrm{min} \ f(\boldsymbol{x}) \\ \mathrm{s.t.} \ \boldsymbol{g}(\boldsymbol{x})=\boldsymbol{0}, \boldsymbol{x} \ge \boldsymbol{0} \end{array} $$

ここで、$\boldsymbol{x} \in \mathbb{R}^n, \boldsymbol{g}(\boldsymbol{x}) \in \mathbb{R}^m$である。

また、問題に$\boldsymbol{x} \ge \boldsymbol{0}$以外の不等式制約が含まれる場合にも、上記の式の形に変換することが出来る。例えば、
$\boldsymbol{h}(\boldsymbol{x}) \le \boldsymbol{0}$
なる不等式制約$\boldsymbol{h}(\boldsymbol{x})$があるとき、新たに変数ベクトル$\boldsymbol{s} \ge \boldsymbol{0}$を導入して、
$ \boldsymbol{h}(\boldsymbol{x}) + \boldsymbol{s} = \boldsymbol{0}$
という等式制約に変換できる。なお、この$\boldsymbol{s}$をスラック変数と呼ぶ。

最適解が満たす条件

次に、上記の最適化問題の最適解が満たす条件（最適性条件）を考える。最適性条件は、ラグランジュの未定乗数法により得られる。
まず、最適化問題のラグランジュ関数を次式で定義する。

$$ L(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}) = f(\boldsymbol{x}) + \boldsymbol{y}^{\top} \boldsymbol{g}(\boldsymbol{x}) - \boldsymbol{z}^{\top} \boldsymbol{x} $$

ここで、$\boldsymbol{y}, \boldsymbol{z}$はラグランジュ乗数ベクトルである。

局所解は、次の5つの式を満たす解である。

$$ \begin{array}{l} \nabla_\boldsymbol{x} L(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}) = \boldsymbol{0} \\ \nabla_\boldsymbol{y} L(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}) = \boldsymbol{g}(\boldsymbol{x}) = \boldsymbol{0} \\ \nabla_\boldsymbol{z} L(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}) = - \boldsymbol{x} \le \boldsymbol{0} \\ \boldsymbol{z} \ge \boldsymbol{0} \\ z_i x_i = 0 \ (i=1, …, m) \end{array} $$

上記の5つの条件を Karush-Kuhn-Tucker条件（KKT条件） と呼ぶ。また、$z_i x_i = 0 \ (i=1, …, m)$を相補性条件 (complementarity condition)と呼ぶ。

さらに、
$ \boldsymbol{w}=[\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}]^{\top} $
とおき、KKT条件を以下のように変形する。

$$ \begin{array}{c} \boldsymbol{r_0} (\boldsymbol{w}) = \left[ \begin{array}{c} \nabla_\boldsymbol{x} L(\boldsymbol{w}) \\ \boldsymbol{g}(\boldsymbol{x}) \\ XZ \boldsymbol{e} \end{array} \right] = \left[ \begin{array}{l} \boldsymbol{0} \\ \boldsymbol{0} \\ \boldsymbol{0} \end{array} \right] \\ \boldsymbol{x} \ge \boldsymbol{0} \\ \boldsymbol{z} \ge \boldsymbol{0} \end{array} $$

ここで、
$$ \begin{array}{c} X = \mathrm{diag}(x_1, x_2, …, x_n) \in \mathbb{R}^{n \times n}, \\ Z = \mathrm{diag}(z_1, z_2, …, z_n) \in \mathbb{R}^{n \times n}, \\ \boldsymbol{e} = [1, 1, …, 1 ] ^{\top} \in \mathbb{R}^n \end{array} $$
である。なお、diagは対角行列である。

主双対内点法のアルゴリズム

中心化KKT条件

変形後のKKT条件をそのまま内点法で解く場合、計算途中の解が$\boldsymbol{x} = \boldsymbol{0}, \boldsymbol{z} = \boldsymbol{0}$の境界に近づくと、収束が遅くなる可能性がある。そこで、バリアパラメータと呼ばれる変数を導入し、計算途中の解を境界から離しながら最適解に近づける。

バリアパラメータ$\mu (> 0)$を導入したKKT条件は次式のようになる。

$$ \begin{array}{c} \boldsymbol{r} (\boldsymbol{w}; \mu) = \left[ \begin{array}{c} \nabla_\boldsymbol{x} L(\boldsymbol{w}) \\ \boldsymbol{g}(\boldsymbol{x}) \\ XZ \boldsymbol{e} - \mu \boldsymbol{e} \end{array} \right] = \left[ \begin{array}{l} \boldsymbol{0} \\ \boldsymbol{0} \\ \boldsymbol{0} \end{array} \right] \\ \boldsymbol{x} \ge \boldsymbol{0} \\ \boldsymbol{z} \ge \boldsymbol{0} \end{array} $$

これを中心化KKT条件またはバリアKKT条件と呼ぶ。$\mu =0$とすると、中心化KKT条件の解はもとのKKT条件の解に等しくなる。主双対内点法では、始めに適当な$\mu$の初期値を与え、解の更新の度に徐々に$\mu$の値を小さくして0に近づけていく。

解の更新

中心化KKT条件の方程式は非線形であるため、最適解を直接求めることは困難である。よって、ニュートン法により解候補を徐々に最適解の方向に更新することによって最適解を求める。

解の更新方向を
$ \boldsymbol{\Delta w} = [\boldsymbol{\Delta x}, \boldsymbol{\Delta y}, \boldsymbol{\Delta z}]^{\top}$
とすると、これを求めるには次の1次方程式を解けばよい。
$J(\boldsymbol{w}) \boldsymbol{\Delta w} = - \boldsymbol{r} (\boldsymbol{w}; \mu)$
ここで、$J(\boldsymbol{w})$は$\boldsymbol{r} (\boldsymbol{w}; \mu)$のヤコビ行列であり、

$$ J(\boldsymbol{w}) = \left[ \begin{array}{ccc} \nabla_x^2 L(\boldsymbol{w}) & \nabla \boldsymbol{g}(\boldsymbol{x}) & -I \\ \nabla \boldsymbol{g}(\boldsymbol{x})^{\top} & O & O \\ Z & O & X \end{array} \right] $$
で与えられる。ただし、$I$は単位行列である。

解の探索方向が得られれば、適当なステップサイズ$\alpha > 0$により、解を以下のように更新する。
$ \boldsymbol{w}_{k+1} = \boldsymbol{w}_k + \alpha \boldsymbol{\Delta w}_k $
ただし、添え字の$k$は更新回数を示す。

以上が、主双対内点法の基本的なアルゴリズムである。

補足

解の更新幅

解の更新幅$\alpha$の決め方として、Armijo（アルミホ）条件を用いた直線探索が参考文献の「最適化とその応用」で紹介されている。また、Ipoptという最適化ソルバでは、フィルター法による直線探索が実装されている。

ヘッセ行列の計算

中心化KKT条件のヤコビ行列$J(\boldsymbol{w})$には、ラグランジュ関数のヘッセ行列
$ \nabla^2_x L(\boldsymbol{w}) $
が現れている。これは、元の最適化問題の目的関数$f(\boldsymbol{x})$のヘッセ行列の計算が必要なことを意味する。しかしながら、ヘッセ行列が得られない場合は、準ニュートン法(L-BFGS法)による近似行列で代用することも可能である。