線形相補性問題とは

線形相補性問題（Linear Complementary Problem, LCP） とは，次の条件を満たす変数$\boldsymbol{z}$を求める問題です．

\boldsymbol{0}\leq\boldsymbol{z}\perp\boldsymbol{M}\boldsymbol{z}+\boldsymbol{p}\geq\boldsymbol{0}

ただし，$\boldsymbol{M}$は非負定値対称行列，$\boldsymbol{p}$は定数ベクトルです．
あるいは次のようにも書けます．

\boldsymbol{w}=\boldsymbol{M}\boldsymbol{z}+\boldsymbol{p},\quad
\boldsymbol{w}\geq\boldsymbol{0},\quad
\boldsymbol{z}\geq\boldsymbol{0},\quad
\boldsymbol{w}^{\mathrm{T}}\boldsymbol{z}=0

$\boldsymbol{w}=[w_{1}~\cdots~w_{n}]^{\mathrm{T}}$と$\boldsymbol{z}=[z_{1}~\cdots~z_{n}]^{\mathrm{T}}$（$n$は$\boldsymbol{z}$の次元数）は同数の成分を持ちますが，$\boldsymbol{w}\geq\boldsymbol{0}$かつ$\boldsymbol{z}\geq\boldsymbol{0}$なので，$\boldsymbol{w}^{\mathrm{T}}\boldsymbol{z}=0$となるためには$w_{i}$と$z_{i}$（$\forall i\in{1,\cdots,n}$）のどちらか一方は必ず$0$になる，ということが「相補性」の意味です．
ちょっと変わったクラスの問題に見えるかも知れませんが，後ほど記すように，ポピュラーな凸2次計画問題も双対性に基づいて線形相補性問題に等価変換できるので，わりと様々な事例がこのように表せます．

線形相補性問題の標準的解法としては，相補掃き出し法（Complementary Pivotting Method, Lemke法） が知られています．

C. E. Lemke, "A Method of Solution for Quadratic Programs", Management Science, 8(4):442-453 (1962).

が，その説明がウェブ上になかなか見当たらないので自分で書くことにしました．なお，本記事の作成にあたって次の書籍を参考にしました．

岩波講座応用数学[方法7] 最適化法，藤田宏，今野浩，田邉國士著，岩波書店，1994

相補掃き出し法のアルゴリズム

まず，次の集合$\mathcal{Q}$を定義します．

\mathcal{Q}=\left\{(\boldsymbol{w},\boldsymbol{z})\left|
\boldsymbol{w}-\boldsymbol{M}\boldsymbol{z}=\boldsymbol{p},\boldsymbol{w}\geq\boldsymbol{0}, \boldsymbol{z}\geq\boldsymbol{0}
\right.\right\}

$\mathcal{Q}$の元のうち条件$\boldsymbol{w}^{\mathrm{T}}\boldsymbol{z}=0$を満たすものが線形相補性問題の解です．もし$\boldsymbol{p}\geq\boldsymbol{0}$ならば，$\boldsymbol{w}=\boldsymbol{p}$，$\boldsymbol{z}=\boldsymbol{0}$は自明解となることに注意して下さい．以降では$\boldsymbol{p}\ngeq\boldsymbol{0}$を仮定します．

天下り的ですが，ここで人為変数$v$を新たに導入して，上記の$\mathcal{Q}$と似た次の集合$\tilde{\mathcal{Q}}$を定義します．

\tilde{\mathcal{Q}}=\left\{(\boldsymbol{w},\boldsymbol{z},v)\left|
\boldsymbol{w}-\boldsymbol{M}\boldsymbol{z}-v\boldsymbol{l}=\boldsymbol{p}, \boldsymbol{w}\geq\boldsymbol{0}, \boldsymbol{z}\geq\boldsymbol{0}, v\geq 0
\right.\right\}

ただし，$\boldsymbol{l}\overset{\mathrm{def}}{=}[1~\cdots~1]^{\mathrm{T}}$とおきました．
$(\boldsymbol{w},\boldsymbol{z},0)\in\tilde{\mathcal{Q}}$であれば$(\boldsymbol{w},\boldsymbol{z})\in\mathcal{Q}$であることは，すぐにご理解頂けるともの思います．

ここで，$\boldsymbol{p}=[p_{1}~\cdots~p_{n}]^{\mathrm{T}}$に対し$s=\mathop{\mathrm{arg~min}}_{i}\left\{p_{i}\right\}$とおけば，$(\boldsymbol{w},\boldsymbol{z},v)=(\boldsymbol{p}-p_{s}\boldsymbol{l},\boldsymbol{0},-p_{s})$は$\tilde{\mathcal{Q}}$の元となり（$\boldsymbol{p}\ngeq\boldsymbol{0}$を仮定しているので$-p_{s}>0$です），かつ$\boldsymbol{w}^{\mathrm{T}}\boldsymbol{z}=0$を満たします．この操作によって，今度は$w_{s}$が$0$になりますので，それと相補関係にある$z_{s}$を$0$から増やしても$\boldsymbol{w}^{\mathrm{T}}\boldsymbol{z}=0$は維持されることになります．

上記のことを，次のような表（タブロー）を作ってもう一度考えてみます．

\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
1      & 2      & \cdots & s-1    & s      & s+1    & \cdots & n      & n+1       & \cdots & n+s       & \cdots & 2n        & 2n+1   & \\
\hline
1      & 0      &        & 0      & 0      & 0      &        & 0      & -m_{11}    &        & -m_{1s}    & \cdots & -m_{1n}    & -1     & p_{1} \\
0      & 1      &        & \vdots & \vdots & \vdots &        &        & \vdots    &        & \vdots    &        & \vdots    & \vdots & \vdots \\
       & 0      &        & 0      &        &        &        &        &           &        &           &        &           &        &        \\
       &        &        & 1      & 0      &        &        &        & -m_{(s-1)1} &        & -m_{(s-1)s} &        & -m_{(s-1)n} & -1     & p_{s-1} \\
\vdots & \vdots & \cdots & 0      & 1      & 0      & \cdots & \vdots & -m_{s1}    & \cdots & -m_{ss}    & \cdots & -m_{sn}    & -1     & p_{s} \\
       &        &        &        & 0      & 1      &        &        & -m_{(s+1)1} &        & -m_{(s+1)s} &        & -m_{(s+1)n} & -1     & p_{s+1} \\
       &        &        & \vdots &        & 0      &        &        & \vdots    &        & \vdots    &        & \vdots    & \vdots & \vdots \\
       &        &        &        & \vdots & \vdots &        & 0      &           &        &           &        &           &        &        \\
0      & 0      &        & 0      & 0      & 0      &        & 1      & -m_{n1}    &        & -m_{ns}    &        & -m_{nn}    & -1     & p_{n} \\
\hline
\end{array}

ただし，$m_{ij}$は$\boldsymbol{M}$の$i$行目$j$列目成分です．$1\sim n$列目が$\boldsymbol{w}$に，$n+1\sim 2n$列目が$\boldsymbol{z}$に，$2n+1$列目が$v$にそれぞれ対応します．最右列は$\boldsymbol{p}$です．
$s=\mathop{\mathrm{arg~min}}_{i}\left\{p_{i}\right\}$を求め，$s$行$2n+1$列目を軸として掃き出しを行うと，次のようになります．

\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
1      & 2      & \cdots & s-1    & s      & s+1    & \cdots & n      & n+1            & \cdots & n+s            & \cdots & 2n             & 2n+1   & \\
\hline
1      & 0      &        & 0      &-1      & 0      &        & 0      & m_{s1}-m_{11}    &        & m_{ss}-m_{1s}    & \cdots & m_{sn}-m_{1n}    & 0      & p_{1}-p_{s} \\
0      & 1      &        & \vdots & \vdots & \vdots &        &        & \vdots         &        & \vdots         &        & \vdots         & \vdots & \vdots \\
       & 0      &        & 0      &        &        &        &        &                &        &                &        &                &        &        \\
       &        &        & 1      &-1      &        &        &        & m_{s1}-m_{(s-1)1} &        & m_{ss}-m_{(s-1)s} &        & m_{sn}-m_{(s-1)n} & 0      & p_{s-1}-p_{s} \\
\vdots & \vdots & \cdots & 0      &-1      & 0      & \cdots & \vdots & m_{s1}          & \cdots & m_{ss}          & \cdots & m_{sn}          & 1      &-p_{s} \\
       &        &        &        &-1      & 1      &        &        & m_{s1}-m_{(s+1)1} &        & m_{ss}-m_{(s+1)s} &        & m_{sn}-m_{(s+1)n} & 0      & p_{s+1}-p_{s} \\
       &        &        & \vdots &        & 0      &        &        & \vdots         &        & \vdots         &        & \vdots         & \vdots & \vdots \\
       &        &        &        & \vdots & \vdots &        & 0      &                &        &                &        &                &        &        \\
0      & 0      &        & 0      &-1      & 0      &        & 1      & m_{s1}-m_{n1}    &        & m_{ss}-m_{ns}    &        & m_{sn}-m_{nn}    & 0      & p_{n}-p_{s} \\
\hline
\end{array}

そして，$2n+1$列目を$s$列目に，$s$列目を$n+s$列目に，$n+s$列目を$2n+1$列目にそれぞれ移動します．

\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
1      & 2      & \cdots & s-1    & 2n+1   & s+1    & \cdots & n      & n+1            & \cdots & s      & \cdots & 2n             & n+s            & \\
\hline
1      & 0      &        & 0      & 0      & 0      &        & 0      & m_{s1}-m_{11}    &        & -1     & \cdots & m_{sn}-m_{1n}    & m_{ss}-m_{1s}    & p_{1}-p_{s} \\
0      & 1      &        & \vdots & \vdots & \vdots &        &        & \vdots         &        & \vdots &        & \vdots         & \vdots         & \vdots \\
       & 0      &        & 0      &        &        &        &        &                &        &        &        &                &                &        \\
       &        &        & 1      & 0      &        &        &        & m_{s1}-m_{(s-1)1} &        & -1     &        & m_{sn}-m_{(s-1)n} & m_{ss}-m_{(s-1)s} & p_{s-1}-p_{s} \\
\vdots & \vdots & \cdots & 0      & 1      & 0      & \cdots & \vdots & m_{s1}          & \cdots & -1     & \cdots & m_{sn}          & m_{ss}          &-p_{s} \\
       &        &        &        & 0      & 1      &        &        & m_{s1}-m_{(s+1)1} &        & -1     &        & m_{sn}-m_{(s+1)n} & m_{ss}-m_{(s+1)s} & p_{s+1}-p_{s} \\
       &        &        & \vdots &        & 0      &        &        & \vdots         &        & \vdots &        & \vdots         & \vdots         & \vdots \\
       &        &        &        & \vdots & \vdots &        & 0      &                &        &        &        &                &                &        \\
0      & 0      &        & 0      & 0      & 0      &        & 1      & m_{s1}-m_{n1}    &        & -1     &        & m_{sn}-m_{nn}    & m_{ss}-m_{ns}    & p_{n}-p_{s} \\
\hline
\end{array}

これで，列の入れ替えは起こったものの形式的には最初と同じ形になりました．
左から$n$列目までの部分は単位行列となっており，最右列の成分が全て$0$以上なので，

$w_{i}=p_{i}-p_{s}$（$i=1,\cdots,s-1,s+1,\cdots,n$），$v=-p_{s}$，$w_{s}=z_{i}=0$（$i=1,\cdots,n$，$i\neq s$）が$\tilde{\mathcal{Q}}$の元
$w_{s}=0$なので，$z_{s}$は非負となる範囲で自由に値を変えられる

というのがこの表の読み方です．

このことを一般化して考えてみましょう．$\boldsymbol{w}$，$\boldsymbol{z}$，$v$をまとめて$\boldsymbol{q}=[\boldsymbol{w}^{\mathrm{T}}~\boldsymbol{z}^{\mathrm{T}}~v]^{\mathrm{T}}=[q_{1}~\cdots~q_{2n+1}]^{\mathrm{T}}$とおきます．これらが，次の式を満たす変数群$\boldsymbol{q}_{\mathrm{B}}$，$\boldsymbol{q}_{\mathrm{N}}$，$q_{\mathrm{D}}$に分けられる状況を考えます．

\boldsymbol{q}_{\mathrm{B}}+\boldsymbol{H}\boldsymbol{q}_{\mathrm{N}}+q_{\mathrm{D}}\boldsymbol{d}=\boldsymbol{r}

ただし，$\boldsymbol{H}$は定数行列，$\boldsymbol{d}$は定数ベクトル，$\boldsymbol{r}$は全ての成分が非負の定数ベクトルで，$\boldsymbol{q}_{\mathrm{B}}$と$\boldsymbol{q}_{\mathrm{N}}$はともに$n$次ベクトルであるものとします．$\boldsymbol{q}_{\mathrm{B}}$を基底変数，$\boldsymbol{q}_{\mathrm{N}}$を非基底変数，$q_{\mathrm{D}}$を駆動変数とそれぞれ呼びます．
$\boldsymbol{q}_{\mathrm{B}}$と$\boldsymbol{q}_{\mathrm{N}}$，$q_{\mathrm{D}}$のインデックスの集合を$\mathcal{I}_{\mathrm{B}}=\left\{i_{\mathrm{B}1},\cdots,i_{\mathrm{B}n}\right\}$，$\mathcal{I}_{\mathrm{N}}=\left\{i_{\mathrm{N}1},\cdots,i_{\mathrm{N}n}\right\}$，$\mathcal{I}_{\mathrm{D}}=\left\{i_{\mathrm{D}}\right\}$とそれぞれおき，

\boldsymbol{q}_{\mathrm{B}}=\begin{bmatrix}
q_{i_{\mathrm{B}1}} \\
 \vdots       \\
q_{i_{\mathrm{B}n}} \\
\end{bmatrix}
,\quad
\boldsymbol{q}_{\mathrm{N}}=\begin{bmatrix}
q_{i_{\mathrm{N}1}} \\
 \vdots       \\
q_{i_{\mathrm{N}n}} \\
\end{bmatrix}
,\quad
\boldsymbol{H}=\begin{bmatrix}
h_{1} & \cdots & h_{1n} \\
 \vdots        & \ddots & \vdots \\
h_{n1} & \cdots & h_{nn} \\
\end{bmatrix}
,\quad
\boldsymbol{d}=\begin{bmatrix}
d_{1} \\
 \vdots       \\
d_{n} \\
\end{bmatrix}
,\quad
\boldsymbol{r}=\begin{bmatrix}
r_{1} \\
 \vdots       \\
r_{n} \\
\end{bmatrix}

と成分表示すると，タブローは次のように書けます．

\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
i_{\mathrm{B}1} & i_{\mathrm{B}2} & \cdots & i_{\mathrm{B}(s-1)} & i_{\mathrm{B}s} & i_{\mathrm{B}(s+1)} & \cdots & i_{\mathrm{B}n} & i_{\mathrm{N}1} & \cdots & i_{\mathrm{N}s} & \cdots & i_{\mathrm{N}n} & i_{\mathrm{D}} & \\
\hline
1      & 0      &        & 0      & 0      & 0      &        & 0      & h_{11}    &        & h_{1s}    & \cdots & h_{1n}    & d_{1}   & r_{1} \\
0      & 1      &        & \vdots & \vdots & \vdots &        &        & \vdots   &        & \vdots   &        & \vdots   & \vdots & \vdots \\
       & 0      &        & 0      &        &        &        &        &          &        &          &        &          &        &        \\
       &        &        & 1      & 0      &        &        &        & h_{(s-1)1} &        & h_{(s-1)s} &        & h_{(s-1)n} & d_{s-1}  & r_{s-1} \\
\vdots & \vdots & \cdots & 0      & 1      & 0      & \cdots & \vdots & h_{s1}    & \cdots & h_{ss}    & \cdots & h_{sn}    & d_{s}   & r_{s} \\
       &        &        &        & 0      & 1      &        &        & h_{(s+1)1} &        & h_{(s+1)s} &        & h_{(s+1)n} & d_{s+1}  & r_{s+1} \\
       &        &        & \vdots &        & 0      &        &        & \vdots   &        & \vdots   &        & \vdots   & \vdots & \vdots \\
       &        &        &        & \vdots & \vdots &        & 0      &          &        &          &        &          &        &        \\
0      & 0      &        & 0      & 0      & 0      &        & 1      & h_{n1}    &        & h_{ns}    &        & h_{nn}    & d_{n}   & r_{n} \\
\hline
\end{array}

軸をどのように選ぶかはいったん置いておくとして，$s$行目$i_{\mathrm{D}}$列目に選んで掃き出しを行うと，次のようになります．

\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
i_{\mathrm{B}1} & i_{\mathrm{B}2} & \cdots & i_{\mathrm{B}(s-1)} & i_{\mathrm{B}s} & i_{\mathrm{B}(s+1)} & \cdots & i_{\mathrm{B}n} & i_{\mathrm{N}1} & \cdots & i_{\mathrm{N}s} & \cdots & i_{\mathrm{N}n} & i_{\mathrm{D}} & \\
\hline
1      & 0      &        & 0      &-d_{1}/d_{s}  & 0      &        & 0      & h_{11}^{\prime}    &        & h_{1s}^{\prime}    & \cdots & h_{1n}^{\prime}    & 0      & r_{1}^{\prime} \\
0      & 1      &        & \vdots & \vdots      & \vdots &        &        & \vdots          &        & \vdots          &        & \vdots          & \vdots & \vdots \\
       & 0      &        & 0      &             &        &        &        &                 &        &                 &        &                 &        &        \\
       &        &        & 1      &-d_{s-1}/d_{s} &        &        &        & h_{(s-1)1}^{\prime} &        & h_{(s-1)s}^{\prime} &        & h_{(s-1)n}^{\prime} & 0      & r_{s-1}^{\prime} \\
\vdots & \vdots & \cdots & 0      & 1/d_{s}      & 0      & \cdots & \vdots & h_{s1}^{\prime}    & \cdots & h_{ss}^{\prime}    & \cdots & h_{sn}^{\prime}    & 1      & r_{s}^{\prime} \\
       &        &        &        &-d_{s+1}/d_{s} & 1      &        &        & h_{(s+1)1}^{\prime} &        & h_{(s+1)s}^{\prime} &        & h_{(s+1)n}^{\prime} & 0      & r_{s+1}^{\prime} \\
       &        &        & \vdots &             & 0      &        &        & \vdots          &        & \vdots          &        & \vdots          & \vdots & \vdots \\
       &        &        &        & \vdots      & \vdots &        & 0      &                 &        &                 &        &                 &        &        \\
0      & 0      &        & 0      &-d_{n}/d_{s}  & 0      &        & 1      & h_{n1}^{\prime}    &        & h_{ns}^{\prime}    &        & h_{nn}^{\prime}    & 0      & r_{n}^{\prime} \\
\hline
\end{array}

ただし，

\displaylines{
h_{ij}^{\prime}=\begin{cases}
h_{ij}-(d_{i}/d_{s})h_{sj} & (i=1,\cdots,n, i\neq s) \\
h_{sj}^{\prime}=h_{sj}/d_{s} & (i=s) \\
\end{cases}
\\
r_{i}^{\prime}=\begin{cases}
r_{i}-(d_{i}/d_{s})r_{s} & (i=1,\cdots,n, i\neq s) \\
r_{s}^{\prime}=r_{s}/d_{s} & (i=s) \\
\end{cases}
}

とおきました．
このように変形しても$r_{i}^{\prime}\geq 0$（$i=1,\cdots,n$）が満たされるためには，$s=\mathop{\mathrm{arg~min}}_{i}\left\{\left.r_{i}/d_{i}\right|d_{i}>0 \right\}$でなければなりません．もし全ての$i=1,\cdots,n$に対して$d_{i}\leq 0$ならば，元の相補性問題に解は無いということになるので，計算を停止します．

軸となる$s$が見つかり上記の掃き出しを行ったら，$i_{\mathrm{D}}$列目を$i_{\mathrm{B}s}$列目に，$i_{\mathrm{B}s}$列目を$i_{\mathrm{N}s}$列目に，$i_{\mathrm{N}s}$列目を$i_{\mathrm{D}}$列目にそれぞれ移動します．

\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
i_{\mathrm{B}1} & i_{\mathrm{B}2} & \cdots & i_{\mathrm{D}} & i_{\mathrm{B}s} & i_{\mathrm{B}(s+1)} & \cdots & i_{\mathrm{B}n} & i_{\mathrm{N}1} & \cdots & i_{\mathrm{B}s} & \cdots & i_{\mathrm{N}n} & i_{\mathrm{N}s} & \\
\hline
1      & 0      &        & 0      & 0           & 0      &        & 0      & h_{11}^{\prime}    &        &-d_{1}/d_{s}    & \cdots & h_{1n}^{\prime} & h_{1s}^{\prime}  & r_{1}^{\prime} \\
0      & 1      &        & \vdots & \vdots      & \vdots &        &        & \vdots          &        & \vdots          &        & \vdots          & \vdots & \vdots \\
       & 0      &        & 0      &             &        &        &        &                 &        &                 &        &                 &        &        \\
       &        &        & 1      & 0           &        &        &        & h_{(s-1)1}^{\prime} &        &-d_{s-1}/d_{s}     &        & h_{(s-1)n}^{\prime} & h_{(s-1)s}^{\prime} & r_{s-1}^{\prime} \\
\vdots & \vdots & \cdots & 0      & 1      & 0      & \cdots & \vdots & h_{s1}^{\prime}    & \cdots & 1/d_{s}   & \cdots & h_{sn}^{\prime}    & h_{ss}^{\prime}     & r_{s}^{\prime} \\
       &        &        &        & 0 & 1      &        &        & h_{(s+1)1}^{\prime} &        &-d_{s+1}/d_{s} &        & h_{(s+1)n}^{\prime} & h_{(s+1)s}^{\prime}     & r_{s+1}^{\prime} \\
       &        &        & \vdots &             & 0      &        &        & \vdots          &        & \vdots          &        & \vdots          & \vdots & \vdots \\
       &        &        &        & \vdots      & \vdots &        & 0      &                 &        &                 &        &                 &        &        \\
0      & 0      &        & 0      & 0  & 0      &        & 1      & h_{n1}^{\prime}    &        &-d_{n}/d_{s}   &        & h_{nn}^{\prime}    & h_{ns}^{\prime}   & r_{n}^{\prime} \\
\hline
\end{array}

その上で，

\displaylines{
\mathcal{I}_{\mathrm{B}}\leftarrow \mathcal{I}_{\mathrm{B}}/\left\{i_{\mathrm{B}s}\right\}\cup\left\{i_{\mathrm{D}}\right\} \\
\mathcal{I}_{\mathrm{N}}\leftarrow \mathcal{I}_{\mathrm{N}}/\left\{i_{\mathrm{N}s}\right\}\cup\left\{i_{\mathrm{B}s}\right\} \\
\mathcal{I}_{\mathrm{D}}\leftarrow \mathcal{I}_{\mathrm{D}}/\left\{i_{\mathrm{D}}\right\}\cup\left\{i_{\mathrm{N}s}\right\}
}

とすれば，形式的に最初と同じ形になります．これにより新たに非基底となった変数インデックス$i_{\mathrm{N}s}$が$2n+1$ならば，$q_{i_{\mathrm{B}s}}=v=0$ということになるので，このときの$\left(\boldsymbol{w},\boldsymbol{z}\right)=\left(\left[q_{1}~\cdots~q_{n}\right]^{\mathrm{T}},\left[q_{n+1}~\cdots~q_{2n}\right]^{\mathrm{T}}\right)$は$\mathcal{Q}$の元となり，かつ$\boldsymbol{w}^{\mathrm{T}}\boldsymbol{z}=0$を満たす，すなわちこれが元の相補性問題の解であるということになります．

以上をまとめれば，相補掃き出し法のアルゴリズムは次のようになります．

InitTableau$(\boldsymbol{M},\boldsymbol{p})$はタブローを初期化する処理
Sweep$(s;\mathcal{T},\mathcal{I}_{\mathrm{B}},\mathcal{I}_{\mathrm{N}},\mathcal{I}_{\mathrm{D}})$はタブローの$s$行目$i_{\mathrm{D}}$列目を軸として掃き出しを行う処理
SwapIndex$(s;\mathcal{I}_{\mathrm{B}},\mathcal{I}_{\mathrm{N}},\mathcal{I}_{\mathrm{D}})$は基底・非基底・駆動変数インデックスを入れ替える処理
FindPivot$(\mathcal{T},\mathcal{I}_{\mathrm{D}})$は軸を見つける処理
BuildAnswer$(\boldsymbol{p},\mathcal{I}_{\mathrm{B}})$は$\boldsymbol{p}$と$\mathcal{I}_{\mathrm{B}}$から解$\boldsymbol{z}$を構成する処理

です．

2次計画問題から線形相補性問題への等価変換

冒頭に述べた通り，凸2次計画問題も双対性に基づいて線形相補性問題に等価変換できます．これを示しましょう．

次の凸2次計画問題を考えます．

\boldsymbol{x}^{*}=\mathop{\mathrm{arg~min}}_{\boldsymbol{x}}\left\{
\left.
z=\frac{1}{2}\boldsymbol{x}^{\mathrm{T}}\boldsymbol{Q}\boldsymbol{x}+\boldsymbol{c}^{\mathrm{T}}\boldsymbol{x}
\right|
\boldsymbol{A}\boldsymbol{x}\geq\boldsymbol{b},
\boldsymbol{x}\geq\boldsymbol{0}
\right\}

ただし，$\boldsymbol{Q}$は実対称非負定値行列，$\boldsymbol{A}$は定数行列，$\boldsymbol{c}$と$\boldsymbol{b}$はどちらも定数ベクトルです．

まず，$\boldsymbol{x}^{*}$が存在するならば，次を満たす$\boldsymbol{\lambda}$も存在します．

\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{x}^{*}=0, \quad
\boldsymbol{y}^{\mathrm{T}}\boldsymbol{\lambda}=0, \quad
\boldsymbol{\mu}\geq\boldsymbol{0}, \quad
\boldsymbol{y}\geq\boldsymbol{0}, \quad
\boldsymbol{x}^{*}\geq\boldsymbol{0}, \quad
\boldsymbol{\lambda}\geq\boldsymbol{0}

ただし，

\displaylines{
\boldsymbol{\mu}\overset{\mathrm{def}}{=}\boldsymbol{Q}\boldsymbol{x}^{*}-\boldsymbol{A}^{\mathrm{T}}\boldsymbol{\lambda}+\boldsymbol{c} \\
\boldsymbol{y}\overset{\mathrm{def}}{=}\boldsymbol{A}\boldsymbol{x}^{*}-\boldsymbol{b} \\
}

とおきました．証明は次の通りです．

$\boldsymbol{A}(\boldsymbol{x}^{*}+\boldsymbol{\varepsilon})\geq\boldsymbol{b}$および$\boldsymbol{x}^{*}+\boldsymbol{\varepsilon}\geq\boldsymbol{0}$を満たす任意の微小なベクトル$\boldsymbol{\varepsilon}$を考えると，$z$は$\boldsymbol{x}=\boldsymbol{x}^{*}$において最小となるので，

\left.z\right|_{\boldsymbol{x}=\boldsymbol{x}^{*}+\boldsymbol{\varepsilon}}\geq\left.z\right|_{\boldsymbol{x}=\boldsymbol{x}^{*}}

すなわち

\frac{1}{2}\boldsymbol{\varepsilon}^{\mathrm{T}}\boldsymbol{Q}\boldsymbol{\varepsilon}+
(\boldsymbol{c}+\boldsymbol{Q}\boldsymbol{x}^{*})^{\mathrm{T}}\boldsymbol{\varepsilon}\geq 0

となるはずです．これが満たされるためには

(\boldsymbol{c}+\boldsymbol{Q}\boldsymbol{x}^{*})^{\mathrm{T}}\boldsymbol{\varepsilon}\geq 0

でなければなりません．
ここで，$\boldsymbol{\varepsilon}$についての線形計画問題

\boldsymbol{\varepsilon}^{*}=\mathop{\mathrm{arg~min}}_{\boldsymbol{\varepsilon}}\left\{
w=(\boldsymbol{c}+\boldsymbol{Q}\boldsymbol{x}^{*})^{\mathrm{T}}\boldsymbol{\varepsilon}
\left|
\boldsymbol{A}\boldsymbol{\varepsilon}\geq-\boldsymbol{y},
\boldsymbol{\varepsilon}\geq-\boldsymbol{x}^{*}
\right.
\right\}

を考えると，$\boldsymbol{\varepsilon}=\boldsymbol{0}$は明らかにこの問題の実行可能解です．
このとき$w=0$ですが，上記より$w\geq 0$ですので，$\boldsymbol{\varepsilon}=\boldsymbol{0}$はこの問題の最適解となっている，つまり$\boldsymbol{\varepsilon}^{*}=\boldsymbol{0}$であると分かります．したがって同時に，上記の線形計画問題の双対問題

(\boldsymbol{\lambda},\boldsymbol{\mu})=\mathop{\mathrm{arg~min}}\_{(\boldsymbol{\lambda},\boldsymbol{\mu})}\left\{
\boldsymbol{y}^{\mathrm{T}}\boldsymbol{\lambda}+\boldsymbol{x}^{*\mathrm{T}}\boldsymbol{\mu}
\left|
\boldsymbol{A}^{\mathrm{T}}\boldsymbol{\lambda}+\boldsymbol{\mu}=\boldsymbol{c}+\boldsymbol{Q}\boldsymbol{x}^{*},
\boldsymbol{\lambda}\geq\boldsymbol{0},
\boldsymbol{\mu}\geq\boldsymbol{0}
\right.
\right\}

も最適解を持ち，しかもそれは$\boldsymbol{y}^{\mathrm{T}}\boldsymbol{\lambda}+\boldsymbol{x}^{*\mathrm{T}}\boldsymbol{\mu}=0$を満たすと分かります．$\boldsymbol{y}\geq\boldsymbol{0}$，$\boldsymbol{\lambda}\geq\boldsymbol{0}$，$\boldsymbol{x}^{*}\geq\boldsymbol{0}$，$\boldsymbol{\mu}\geq\boldsymbol{0}$ですので，結局$\boldsymbol{y}^{\mathrm{T}}\boldsymbol{\lambda}=0$かつ$\boldsymbol{x}^{*\mathrm{T}}\boldsymbol{\mu}=0$が成り立つことになります（証明終わり）．

逆に，上記を満たす$\boldsymbol{\lambda}$が存在するならば$\boldsymbol{x}^{*}$は元の凸2次計画問題の最適解となります．証明は次の通りです．

z-\left.z\right|_{\boldsymbol{x}=\boldsymbol{x}^{*}}
=
\frac{1}{2}(\boldsymbol{x}-\boldsymbol{x}^{*})^{\mathrm{T}}\boldsymbol{Q}(\boldsymbol{x}-\boldsymbol{x}^{*})
+(\boldsymbol{c}+\boldsymbol{Q}\boldsymbol{x}^{*})^{\mathrm{T}}(\boldsymbol{x}-\boldsymbol{x}^{*})

であるので，常に

z-\left.z\right|_{\boldsymbol{x}=\boldsymbol{x}^{*}}
\geq
(\boldsymbol{c}+\boldsymbol{Q}\boldsymbol{x}^{*})^{\mathrm{T}}(\boldsymbol{x}-\boldsymbol{x}^{*})

が言えます．$\boldsymbol{x}$が元の問題の実行可能解ならば，$\boldsymbol{A}\boldsymbol{x}\geq\boldsymbol{b}$，$\boldsymbol{x}\geq\boldsymbol{0}$，$\boldsymbol{\mu}\geq\boldsymbol{0}$，$\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{x}^{*}=0$かつ

\boldsymbol{c}+\boldsymbol{Q}\boldsymbol{x}^{*}=\boldsymbol{A}^{\mathrm{T}}\boldsymbol{\lambda}+\boldsymbol{\mu}

ですので，

\displaylines{
z-\left.z\right|_{\boldsymbol{x}=\boldsymbol{x}^{*}}
\geq
(\boldsymbol{A}^{\mathrm{T}}\boldsymbol{\lambda}+\boldsymbol{\mu})^{\mathrm{T}}(\boldsymbol{x}-\boldsymbol{x}^{*})
=\boldsymbol{\lambda}^{\mathrm{T}}\boldsymbol{A}(\boldsymbol{x}-\boldsymbol{x}^{*})
+\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{x}-\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{x}^{*}
\\
\geq-\boldsymbol{\lambda}^{\mathrm{T}}(\boldsymbol{A}\boldsymbol{x}^{*}-\boldsymbol{b})=-\boldsymbol{\lambda}^{\mathrm{T}}\boldsymbol{y}=0
}

すなわち

z\geq\left.z\right|_{\boldsymbol{x}=\boldsymbol{x}^{*}}

となります（証明終わり）．

以上より，元の凸2次計画問題を解くことは，

\displaylines{
\boldsymbol{\mu}=\boldsymbol{Q}\boldsymbol{x}^{*}-\boldsymbol{A}^{\mathrm{T}}\boldsymbol{\lambda}+\boldsymbol{c} \\
\boldsymbol{y}=\boldsymbol{A}\boldsymbol{x}^{*}-\boldsymbol{b} \\
\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{x}^{*}=0, \quad
\boldsymbol{y}^{\mathrm{T}}\boldsymbol{\lambda}=0 \\
\boldsymbol{\mu}\geq\boldsymbol{0}, \quad
\boldsymbol{y}\geq\boldsymbol{0}, \quad
\boldsymbol{x}^{*}\geq\boldsymbol{0}, \quad
\boldsymbol{\lambda}\geq\boldsymbol{0}
}

を満たす$\boldsymbol{x}^{*}$と$\boldsymbol{\lambda}$の組を求めることと等価であると分かります．これは実は

\boldsymbol{w}\overset{\mathrm{def}}{=}\begin{bmatrix}
\boldsymbol{\mu} \\ \boldsymbol{y}
\end{bmatrix}, \quad
\boldsymbol{z}\overset{\mathrm{def}}{=}\begin{bmatrix}
\boldsymbol{x}^{*} \\ \boldsymbol{\lambda}
\end{bmatrix}, \quad
\boldsymbol{M}\overset{\mathrm{def}}{=}\begin{bmatrix}
\boldsymbol{Q} & -\boldsymbol{A}^{\mathrm{T}} \\
\boldsymbol{A} & \boldsymbol{O}
\end{bmatrix}, \quad
\boldsymbol{p}\overset{\mathrm{def}}{=}\begin{bmatrix}
\boldsymbol{c} \\ -\boldsymbol{b}
\end{bmatrix}

とおいたときに，本記事冒頭に示した線形相補性問題になっています．

変数に非負条件$\boldsymbol{x}\geq\boldsymbol{0}$が課されない場合は，二つの非負ベクトル$\boldsymbol{x}_{+}\geq\boldsymbol{0}$，$\boldsymbol{x}_{-}\geq\boldsymbol{0}$を用いて$\boldsymbol{x}=\boldsymbol{x}_{+}-\boldsymbol{x}_{-}$とおけば良いです．すなわち

\displaylines{
\boldsymbol{x}^{*}=\mathop{\mathrm{arg~min}}_{\boldsymbol{x}}
\left\{
\left.
z=\frac{1}{2}(\boldsymbol{x}_{+}-\boldsymbol{x}_{-})^{\mathrm{T}}\boldsymbol{Q}(\boldsymbol{x}_{+}-\boldsymbol{x}_{-})+\boldsymbol{c}^{\mathrm{T}}(\boldsymbol{x}_{+}-\boldsymbol{x}_{-})
\right|
\boldsymbol{A}(\boldsymbol{x}_{+}-\boldsymbol{x}_{-})\geq\boldsymbol{b},
\boldsymbol{x}_{+}\geq\boldsymbol{0},
\boldsymbol{x}_{-}\geq\boldsymbol{0}
\right\}
\\
=\mathop{\mathrm{arg~min}}_{\boldsymbol{x}}\left\{
\left.
z=\frac{1}{2}\tilde{\boldsymbol{x}}^{\mathrm{T}}\tilde{\boldsymbol{Q}}\tilde{\boldsymbol{x}}+\tilde{\boldsymbol{c}}^{\mathrm{T}}\tilde{\boldsymbol{x}}
\right|
\tilde{\boldsymbol{A}}\tilde{\boldsymbol{x}}\geq\boldsymbol{b},
\tilde{\boldsymbol{x}}\geq\boldsymbol{0}
\right\}
\\
\tilde{\boldsymbol{x}}\overset{\mathrm{def}}{=}\begin{bmatrix} \boldsymbol{x}_{+}^{\mathrm{T}} & \boldsymbol{x}_{-}^{\mathrm{T}}\end{bmatrix}^{\mathrm{T}}
\\
\tilde{\boldsymbol{Q}}\overset{\mathrm{def}}{=}\begin{bmatrix} \boldsymbol{Q} & -\boldsymbol{Q} \\ -\boldsymbol{Q} & \boldsymbol{Q} \end{bmatrix}
\\
\tilde{\boldsymbol{c}}\overset{\mathrm{def}}{=}\begin{bmatrix} \boldsymbol{c}^{\mathrm{T}} & -\boldsymbol{c}^{\mathrm{T}} \end{bmatrix}^{\mathrm{T}}
\\
\tilde{\boldsymbol{A}}\overset{\mathrm{def}}{=}\begin{bmatrix} \boldsymbol{A} & -\boldsymbol{A} \end{bmatrix}
}

と変形することで，相補掃き出し法を適用できます．

ZMでの実装

ZMでは，相補掃き出し法による線形相補性問題求解関数zLCPSolveLemke()をzm_opt_lcp_lemke.cに，またこれを用いた凸2次計画問題求解関数zQPSolveLemke()をzm_opt_qp.cに，それぞれ実装しています．