逐次$2$次計画(SOP; Successive Quadratic Programming)法は制約付き最適化問題の近似解を計算するにあたって用いられる手法です。当記事では逐次$2$次計画法の概要や数式、Pythonを用いた実際の計算例について取り扱いました。
「新版 数理計画入門(朝倉書店)」の$4.9$節の「逐次$2$次計画法」などの内容を参考に作成を行いました。
・用語/公式解説
https://www.hello-statisticians.com/explain-terms
前提の確認
問題設定
当記事では下記のような目的関数$f(\mathbf{x})$の制約付き最適化問題を定義し、以下取り扱う。
$$
\large
\begin{align}
\mathrm{Objective} : \quad & f(\mathbf{x}) \, \rightarrow \, \mathrm{Minimize} \\
\mathrm{Constraint} : \quad & c_i(\mathbf{x}) = 0 \quad (i=1, 2, \cdots ,l)
\end{align}
$$
上記に対し、KKT条件は下記のように表される。
$$
\large
\begin{align}
& \nabla f(\mathbf{x}) + \sum_{i=1}^{n} \lambda_i \nabla c_i(\mathbf{x}) = \mathbf{0} \quad (1) \\
& c_i(\mathbf{x}) = 0 \qquad (i=1, \cdots, l)
\end{align}
$$
また、ラグランジュ関数$L(\mathbf{x},\boldsymbol{\lambda})$を下記のように定義する。
$$
\large
\begin{align}
L(\mathbf{x},\boldsymbol{\lambda}) = f(\mathbf{x}) + \sum_{i=1}^{n} \lambda_i c_i(\mathbf{x})
\end{align}
$$
このときラグランジュ関数$L(\mathbf{x},\boldsymbol{\lambda})$の勾配ベクトル$\nabla L(\mathbf{x},\boldsymbol{\lambda})$は下記のように表すことができる。
$$
\large
\begin{align}
\nabla L(\mathbf{x},\boldsymbol{\lambda}) &= \left(\begin{array}{c} \nabla_{x} L(\mathbf{x},\boldsymbol{\lambda}) \\ \nabla_{\lambda} L(\mathbf{x},\boldsymbol{\lambda}) \end{array} \right) \\
&= \left(\begin{array}{c} \displaystyle \nabla f(\mathbf{x}) + \sum_{i=1}^{n} \lambda_i \nabla c_i(\mathbf{x}) \\ \mathbf{c}(\mathbf{x}) \end{array} \right) \quad (1) \\
\mathbf{c}(\mathbf{x}) &= \left(\begin{array}{c} c_1(\mathbf{x}) \\ \vdots \\ c_l(\mathbf{x}) \end{array} \right)
\end{align}
$$
ここでKKT条件は$(1)$式を元に下記のように表すことができる。
$$
\large
\begin{align}
\nabla L(\mathbf{x},\boldsymbol{\lambda}) = \left(\begin{array}{c} \displaystyle \nabla f(\mathbf{x}) + \sum_{i=1}^{n} \lambda_i \nabla c_i(\mathbf{x}) \\ \mathbf{c}(\mathbf{x}) \end{array} \right) = \left(\begin{array}{c} 0 \\ \vdots \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right) = \mathbf{0}
\end{align}
$$
多次元ベクトル値関数に基づく連立方程式とニュートン法
前項ではKKT条件が$\nabla L(\mathbf{x},\boldsymbol{\lambda}) = \mathbf{0}$のように表せることの確認を行なった。このKKT条件を$\mathbf{x},\boldsymbol{\lambda}$について解くにあたっては、$1$次のテイラー展開に基づくニュートン法を活用することで近似解が得られる。
当項では以下、$p$次元ベクトル$\mathbf{z}$を元に定義する$p$次元ベクトル値関数$\mathbf{F}(\mathbf{z})$に関する方程式$\mathbf{F}(\mathbf{z})=\mathbf{0}$について、ニュートン法に基づいて近似解を得る手順について確認を行う。
まず$p$次元ベクトル$\mathbf{z}$、$p$次元ベクトル値関数$\mathbf{F}(\mathbf{z})$を下記のように定義する。
$$
\large
\begin{align}
\mathbf{z} &= \left(\begin{array}{c} z_1 \\ \vdots \\ z_p \end{array} \right) \\
\mathbf{F}(\mathbf{z}) &= \left(\begin{array}{c} F_1(\mathbf{z}) \\ \vdots \\ F_p(\mathbf{z}) \end{array} \right)
\end{align}
$$
上記で定めた$p$次元ベクトル値関数$\mathbf{F}(\mathbf{z})$に対して、勾配ベクトル$\nabla F_i(\mathbf{z})$を第$i$列とする$p \times p$正方行列を$\nabla \mathbf{F}(\mathbf{z})$とおくと、$\nabla \mathbf{F}(\mathbf{z})$は下記のように表せる。
$$
\large
\begin{align}
\nabla \mathbf{F}(\mathbf{z}) &= \left(\begin{array}{ccc} \nabla F_1(\mathbf{z}) & \cdots & \nabla F_p(\mathbf{z}) \end{array} \right) \\
&= \left(\begin{array}{ccc} \displaystyle \frac{\partial F_1(\mathbf{z})}{\partial z_1} & \cdots & \displaystyle \frac{\partial F_p(\mathbf{z})}{\partial z_1} \\ \vdots & \ddots & \vdots \\ \displaystyle \frac{\partial F_1(\mathbf{z})}{\partial z_p} & \cdots & \displaystyle \frac{\partial F_p(\mathbf{z})}{\partial z_p} \end{array} \right)
\end{align}
$$
上記に基づいてニュートン法の漸化式が下記のように得られる。
$$
\large
\begin{align}
\mathbf{F}(\mathbf{z}^{(k)}) + \nabla \mathbf{F}(\mathbf{z}^{(k)})^{\mathrm{T}} (\mathbf{z}^{(k+1)} – \mathbf{z}^{(k)}) &= \mathbf{0} \\
\mathbf{z}^{(k+1)} – \mathbf{z}^{(k)} &= -\left[ \nabla \mathbf{F}(\mathbf{z}^{(k)})^{\mathrm{T}} \right]^{-1} \mathbf{F}(\mathbf{z}^{(k)}) \\
\mathbf{z}^{(k+1)} &= \mathbf{z}^{(k)} – \left[ \nabla \mathbf{F}(\mathbf{z}^{(k)})^{\mathrm{T}} \right]^{-1} \mathbf{F}(\mathbf{z}^{(k)})
\end{align}
$$
逐次2次計画法の数式
「多次元ベクトル値関数に基づく連立方程式とニュートン法」では$\mathbf{F}(\mathbf{z})=\mathbf{0}$を取り扱ったが、逐次$2$次計画法ではKKT条件の式$\nabla L(\mathbf{x},\boldsymbol{\lambda}) = \mathbf{0}$に対してニュートン法を適用する。
KKT条件の式$\nabla L(\mathbf{x},\boldsymbol{\lambda}) = \mathbf{0}$に関するニュートン法の漸化式は下記のように表すことができる。
$$
\large
\begin{align}
\nabla L(\mathbf{x},\boldsymbol{\lambda}) + \nabla^{2} L(\mathbf{x},\boldsymbol{\lambda}) \left(\begin{array}{c} \mathbf{x}^{(k+1)} – \mathbf{x}^{(k)} \\ \boldsymbol{\lambda}^{(k+1)} – \boldsymbol{\lambda}^{(k)} \end{array} \right) &= \mathbf{0} \\
\left(\begin{array}{c} \mathbf{x}^{(k+1)} \\ \boldsymbol{\lambda}^{(k+1)} \end{array} \right) &= \left(\begin{array}{c} \mathbf{x}^{(k)} \\ \boldsymbol{\lambda}^{(k)} \end{array} \right) – \left[ \nabla^{2} L(\mathbf{x},\boldsymbol{\lambda}) \right]^{-1} \nabla L(\mathbf{x},\boldsymbol{\lambda})
\end{align}
$$
ここで上記の$\nabla L(\mathbf{x},\boldsymbol{\lambda}), \nabla^{2} L(\mathbf{x},\boldsymbol{\lambda})$は下記のように表される。
$$
\large
\begin{align}
\nabla L(\mathbf{x},\boldsymbol{\lambda}) &= \left(\begin{array}{c} \displaystyle \nabla f(\mathbf{x}) + \sum_{i=1}^{n} \lambda_i \nabla c_i(\mathbf{x}) \\ \mathbf{c}(\mathbf{x}) \end{array} \right) \\
\nabla^{2} L(\mathbf{x},\boldsymbol{\lambda}) &= \left(\begin{array}{cc} \nabla_{x}^{2} L(\mathbf{x},\boldsymbol{\lambda}) & \nabla \mathbf{c}(\mathbf{x}) \\ \nabla \mathbf{c}(\mathbf{x})^{\mathrm{T}} & \mathbf{0} \end{array} \right) \\
\nabla \mathbf{c}(\mathbf{x}) &= \left(\begin{array}{ccc} \nabla c_1(\mathbf{x}) & \cdots & \nabla c_l(\mathbf{x}) \end{array} \right)
\end{align}
$$