メロートラの予測子修正子法

メロートラの予測子修正子法（メロートラのよそくししゅうせいしほう、英: Mehrotra's predictor–corrector method）とは数理最適化において線形計画問題に対する内点法の一種である。1989年にサンジェイ・メロートラによって提案された^[1]。

予測子修正子法では探索方向を求めるためにコレスキー分解を用いて大規模な行列を必要に応じて計算し、反復する内点法である。行列の分解ステップでは計算量の大きいステップである。このことは行列した分解を複数回の反復を通じて使用することで実用上の計算量を抑えることができる。

予測子修正子法の各反復では予測方向・修正方向の二つの探索方向を求めるために同一のコレスキー分解した行列を使用する。

予測子修正子法の基本的なアイデアとして、始めに最適解への探索方向を線形の方程式系を解いて決定する（予測子）。続いて最適解への探索方向に対するステップサイズを求めて中心への探索方向も求める。このとき中心方向と2次の方程式系によって決定される（修正子）。

予測子修正子法の探索方向は予測方向・修正方向の和をとることで求まる。

メロートラの予測子修正子法は大変実践向きの内点法ではあるが、多項式オーダーの証明は今のところされていない^[2]。修正方向では予測方向で用いたコレスキー分解した行列をもう一度使用して計算することから効率よく反復を行うことができるため、他の内点法と比べても計算にかかる時間はわずかとされている。しかしながら、反復における追加の計算も最適解に到達するまでに必要な反復回数も十分減少することから、計算にかかる全体の時間も大きなものではないとされる。また最適解に十分に近い点では非常に早く収束することが知られている。

導出

Nocedal、Wrightによってまとめられた導出の流れについて説明する^[3]。

予測ステップ - アフィンスケーリング方向

線形計画問題を以下の標準形に書き直す。これは任意の線形計画問題に対して変換することができる。

${\displaystyle {\begin{aligned}&{\underset {x}{\min }}&q(x)&=c^{T}x,\\&{\text{s.t.}}&Ax&=b,\\&\;&x&\geq 0,\end{aligned}}}$

ただし、${\displaystyle c\in \mathbb {R} ^{n\times 1}}$、${\displaystyle \;A\in \mathbb {R} ^{m\times n}}$、${\displaystyle b\in \mathbb {R} ^{m\times 1}}$ によって ${\displaystyle m}$ 個の制約と ${\displaystyle n}$ 個の等式制約が定義され、${\displaystyle x\in \mathbb {R} ^{n\times 1}}$ は変数ベクトルを表す。

上記の問題に対するKKT条件は以下のように表される:

${\displaystyle {\begin{aligned}A^{T}\lambda +s&=c,\;\;\;{\text{(Lagrange gradient condition)}}\\Ax&=b,\;\;\;{\text{(Feasibility condition)}}\\XSe&=0,\;\;\;{\text{(Complementarity condition)}}\\(x,s)&\geq 0,\end{aligned}}}$

ただし、${\displaystyle X={\text{diag}}(x)}$、${\displaystyle S={\text{diag}}(s)}$ （${\displaystyle x,s}$ を対角成分に並べた行列）であり、 ${\displaystyle e=(1,1,\dots ,1)^{T}\in \mathbb {R} ^{n\times 1}}$ (成分がすべて 1 のベクトル)である。

KKT条件を整理して以下のように ${\displaystyle F:\mathbb {R} ^{2n+m}\rightarrow \mathbb {R} ^{2n+m}}$ と表すとする。

${\displaystyle {\begin{aligned}F(x,\lambda ,s)={\begin{bmatrix}A^{T}\lambda +s-c\\Ax-b\\XSe\end{bmatrix}}&=0\\(x,s)&\geq 0\end{aligned}}}$

予測子修正子法はニュートン方程式を解いてアフィンスケーリング方向を求める。ニュートン方程式は以下のような線形方程式系で表される:

${\displaystyle J(x,\lambda ,s){\begin{bmatrix}\Delta x^{\text{aff}}\\\Delta \lambda ^{\text{aff}}\\\Delta s^{\text{aff}}\end{bmatrix}}=-F(x,\lambda ,s)}$

ただし、${\displaystyle J}$ は

${\displaystyle J(x,\lambda ,s)={\begin{bmatrix}\nabla _{x}F&\nabla _{\lambda }F&\nabla _{s}F\end{bmatrix}},}$

と、F のヤコビ行列である。

すなわち、線形方程式系は以下のように表される:

${\displaystyle {\begin{bmatrix}0&A^{T}&I\\A&0&0\\S&0&X\end{bmatrix}}{\begin{bmatrix}\Delta x^{\text{aff}}\\\Delta \lambda ^{\text{aff}}\\\Delta s^{\text{aff}}\end{bmatrix}}={\begin{bmatrix}-r_{c}\\-r_{b}\\-XSe\end{bmatrix}},\;\;\;r_{c}=A^{T}\lambda +s-c,\;\;\;r_{b}=Ax-b}$

中心化ステップ

${\displaystyle x_{i}s_{i},\;i=1,2,\dots ,n}$ の積の平均値は現在の点 ${\displaystyle (x^{k},s^{k})}$ が（ここで ${\displaystyle k}$ は現在の反復の回数を表す。）最適解にどれ程近づいたかを表す重要な指標となる。この指標は双対ギャップを表しており、以下のように定義される:

${\displaystyle \mu ={\frac {1}{n}}\sum _{i=1}^{n}x_{i}s_{i}={\frac {x^{T}s}{n}}.}$

中心化パラメータ ${\displaystyle \sigma \in [0,1],}$ を用いて以下の方程式系の解を求める:

${\displaystyle {\begin{bmatrix}0&A^{T}&I\\A&0&0\\S&0&X\end{bmatrix}}{\begin{bmatrix}\Delta x^{\text{cen}}\\\Delta \lambda ^{\text{cen}}\\\Delta s^{\text{cen}}\end{bmatrix}}={\begin{bmatrix}-r_{c}\\-r_{b}\\-XSe+\sigma \mu e\end{bmatrix}}}$

修正ステップ

は上記の方程式系によって求めたアフィンスケーリング方向にそのまま進もうとすると相補性条件が満たされなくなることが分かる。

${\displaystyle \left(x_{i}+\Delta x_{i}^{\text{aff}}\right)\left(s_{i}+\Delta s_{i}^{\text{aff}}\right)=x_{i}s_{i}+x_{i}\Delta s_{i}^{\text{aff}}+s_{i}\Delta x_{i}^{\text{aff}}+\Delta x_{i}^{\text{aff}}\Delta s_{i}^{\text{aff}}=\Delta x_{i}^{\text{aff}}\Delta s_{i}^{\text{aff}}\neq 0.}$

そのため、相補性条件の誤差を修正するための方程式系は以下ののように定義される。ただし、この方程式系の係数行列はアフィンスケーリング方向で用いた方程式系の係数行列と等価であることから、ここでの計算量は以前求めた行列に依存する:

${\displaystyle {\begin{bmatrix}0&A^{T}&I\\A&0&0\\S&0&X\end{bmatrix}}{\begin{bmatrix}\Delta x^{\text{cor}}\\\Delta \lambda ^{\text{cor}}\\\Delta s^{\text{cor}}\end{bmatrix}}={\begin{bmatrix}0\\0\\-\Delta X^{\text{aff}}\Delta S^{\text{aff}}e\end{bmatrix}}}$

中心化・修正方向を集約した方程式系

修正方向では予測・中心化方向の方程式系の右辺を一つの方程式系に集約する。この方程式系の計算量についてもアフィンスケーリング方向で求めた係数行列の計算によって決定される。したがってこの方程式系の係数行列は予測方向で使用した分解した行列を再度用いる。

集約した方程式系は:

${\displaystyle {\begin{bmatrix}0&A^{T}&I\\A&0&0\\S&0&X\end{bmatrix}}{\begin{bmatrix}\Delta x\\\Delta \lambda \\\Delta s\end{bmatrix}}={\begin{bmatrix}-r_{c}\\-r_{b}\\-XSe-\Delta X^{\text{aff}}\Delta S^{\text{aff}}e+\sigma \mu e\end{bmatrix}}}$

予測子修正子法は始めにアフィンスケーリング方向を求める。続いて現在の反復点からの探索方向を求める。

中心パラメータの決定方法

アフィンスケーリング方向において中心化パラメータ ${\displaystyle \sigma }$ はヒューリスティックに決定される:

${\displaystyle \sigma =\left({\frac {\mu _{\text{aff}}}{\mu }}\right)^{3}}$

ただし、

${\displaystyle {\begin{aligned}\mu _{\text{aff}}&=(x+\alpha _{\text{aff}}^{\text{pri}}\Delta x^{\text{aff}})^{T}(s+\alpha _{\text{aff}}^{\text{dual}}\Delta s^{\text{aff}})/n,\\\alpha _{\text{aff}}^{\text{pri}}&=\min \left(1,{\underset {i:\Delta x_{i}^{\text{aff}}<0}{\min }}-{\frac {x_{i}}{\Delta x_{i}^{\text{aff}}}}\right),\\\alpha _{\text{aff}}^{\text{dual}}&=\min \left(1,{\underset {i:\Delta s_{i}^{\text{aff}}<0}{\min }}-{\frac {s_{i}}{\Delta s_{i}^{\text{aff}}}}\right),\end{aligned}}}$

であり、${\displaystyle \mu _{\text{aff}}}$ はアフィン関された空間における双対ギャップを表す指標であり、${\displaystyle \mu }$ は更新前の点における双対ギャップを表す指標である^[3]。

ステップ長

実用上で使用されるステップ長は非負条件 ${\displaystyle (x,s)\geq 0}$ を満たす最大ステップ長を採用するため、直線探索によって求められる^[3]。

二次計画問題への適用

メロートラ内点法はの線形計画問題に対する解法ではあるが、その修正版が二次計画問題に対しても拡張することができる^[3]。