はじめに
本記事の目的は，数理計画の分野において重要な役割を果たす相補性問題の基本について解説し，最も基本となる線形相補性問題に対するアルゴリズムを提示*1し，Python によって数値実験を行なうことである．数理最適化の分野の発展に伴い，最適化問題に馴染みの深い読者は多いと思われるが，相補性問題はその重要性の割にあまり知れ渡っていないという印象である．以下，本記事における注意点を簡単に述べておく：
注意1： 本記事では，ベクトルを表す文字であっても太字で表記していない．すなわち，$x \in \mathbb{R}^n$ のような書き方をしており，$\mathbf{x} \in \mathbb{R}^n$ のような書き方をしていない．これは，行列やベクトル関数に対しても同様であり，すべて通常の字体で表記する．
注意2： 本記事の内容の多くは "A semismooth equation approach to the solution of nonlinear complementarity problems"(T. De Luca, F. Facchinei and C. Kanzow) に基づいている．したがって，細かい定理の証明などに関してはそちらを参照されたい．
注意3： 相補性問題を解くためのコードを提示しているが，あくまで線形相補性問題に対するコードであり，一般の相補性問題には対応していない．また，アルゴリズムの収束性に関する証明は扱っていないので，そのような内容が気になる場合は原著を参照されたい．
相補性問題の定式化
はじめに，相補性問題と呼ばれる問題がどのような問題なのかについて簡単に述べる．
与えられたベクトル関数 $F:\mathbb{R}^n\to\mathbb{R}^n$ *2に対して，相補性条件*3を満たすような $n$ 次元ベクトル $x\in\mathbb{R}^{n}$ を求める問題を相補性問題(Complementarity Problem)と呼び，次のように定式化される：
$$ \begin{align} \mathrm{Find}\quad &x \in \mathbb{R}^n \\ \mathrm{such\ that}\quad & x \geqq 0 ,\, F(x) \geqq 0 ,\, \left< x, F(x) \right> = 0 \end{align} $$
ここに，記号 $\left< \cdot, \cdot \right>$ は内積を表す． 特に，ベクトル関数 $F:\mathbb{R}^n\to\mathbb{R}^n$ が $F(x) = Mx+q$ と書ける，すなわち $F$ がアフィン関数のとき，相補性問題は線形相補性問題(Linear Complementarity Problem; LCP)と呼ばれ，次のように定式化される：
$$ \begin{align} \mathrm{Find}\quad &x \in \mathbb{R}^n \\ \mathrm{such\ that}\quad & x \geqq 0 ,\, Mx+q \geqq 0 ,\, \left< x, Mx+q \right> = 0 \end{align} $$ ここに，$M \in \mathbb{R}^{n\times n}$ および $q\in\mathbb{R}^n$ は所与の正方行列およびベクトルをそれぞれ表す．
さて，相補性問題を初めて見ると，なぜこのような問題を考える必要があるのかと不思議に思うかもしれないが，実は相補性問題の定式化能力は非常に高く，最適化問題のみならず均衡問題と呼ばれる問題*4などに対しても定式化することができ，非常に有用なのである．より広いクラスの問題として 変分不等式問題(Variational Inequality Problem; VIP)と呼ばれる問題があり，理論および応用どちらの面でも非常に有用なのだが，本稿では扱わないこととする．変分不等式問題が気になる方は "variational inequality problem" で検索してみると様々な情報が得られるので，それらを参照されたい*5．
相補性問題の再定式化
さて，相補性問題を定式化することができたが，どのように解いたらよいのだろうか．
相補性問題に対する解法アルゴリズムに関しては様々な研究があり，それらの多くは，等価なベクトル方程式に帰着させ，そのベクトル方程式を数理最適化の枠組みで解くというものである*6．例えば，Fischer によって初めて提唱された関数*7 $\varphi(a,b):= \sqrt{a^2 + b^2} - a - b$ を用いることによって相補性問題を再定式化することがよくある．この関数の特徴は次のような関係が成り立つことである： $$ \begin{align} \varphi(a,b) = 0 &\Leftrightarrow \sqrt{a^2+b^2} - a-b = 0 \\ &\Leftrightarrow a + b = \sqrt{a^2+b^2} \\ &\Leftrightarrow a+b \geqq 0, \, (a+b)^2 = a^2+b^2 \\ &\Leftrightarrow a\geqq 0 ,\, b \geqq 0,\, ab=0 \end{align} $$ この関数 $\varphi$ を用いて，関数 $\Phi_i(x):\mathbb{R}^n \to \mathbb{R} : (i=1,2,\dots,n)$ を次のように定義する： $$ \begin{align} \Phi_i(x) = \varphi(x_i, F_i(x)) \end{align} $$ すると，相補性問題は次のベクトル方程式と等価になる： $$ \Phi(x) = \left( \Phi_1(x) , \Phi_2(x) , \ldots, \Phi_n(x) \right)^\top = 0 $$ このベクトル方程式の解は次の最適化問題を解くことによって得られる*8： $$ \begin{align} \mathrm{min} \quad &\Psi(x) := \frac{1}{2} \| \Phi(x) \|_2 = \frac{1}{2} \sum_{i=1}^{n} \varphi(x_i,F_i(x))^2 \\ \mathrm{s.t.} \quad &x \in \mathrm{R}^n \end{align} $$ 次に，この無制約最適化問題に対するアルゴリズムについて述べる．アルゴリズムは Newton 法をベースにしたものであるが，関数 $\Psi(x)$ がある $i$ に対して，$(x_i,F_i(x))=(0,0)$ であるような点において微分不可能であることに注意が必要である．
相補性問題に対するアルゴリズム
早速だが，アルゴリズムを提示する．このアルゴリズムは Luca et.al. (1996) によって提唱されたものである．
Step0: パラメータとして，$\rho > 0, \, p > 2 ,\, \varsigma \in (0,\frac{1}{2}) ,\, \varepsilon > 0$ を設定する．また，初期点 $x_0 \in \mathbb{R}^n$ を任意に設定し，$k \gets 0$ とする．
Step1: 終了条件: $ \| \nabla \Psi (x^k) \|_2 \leqq \varepsilon $ を満たすならば反復を終了し，暫定解 $x^{k}$ を相補性問題の解として出力する．
Step2: ベクトル値関数 $\Phi$ の点 $x^k$ における B 劣微分 $V_k \in \partial_B \Phi(x^{k})$ を計算し，方程式： $$V_k d = - \Phi (x^{k})$$ の解 $d^{k} \in \mathbb{R}^n$ を求める．このとき，方程式の解が存在しない，もしくは，(解が存在したとしても)以下の関係式： $$\nabla \Psi (x^{k})^\top d^{k} \leqq -\rho \|d^{k}\|^{p}$$ が成り立たない場合*9，$d^{k} \in \mathbb{R}^n$ を次のように定める： $$d^{k} := - \nabla \Psi (x^{k}) $$
Step3: 次の関係式を満たすような最小の添え字 $i^k \in \{0,1,2,\ldots \}$ を求める： $$ \begin{align} \Psi \left( x^{k} + 2^{-i^k} d^{k} \right) \leqq \Psi ( x^{k}) + \frac{ \varsigma }{ 2^{i^k} } \nabla \Psi( x^{k} )^{\top} d^{k} \end{align} $$ 暫定解を次式で更新する： $$ \begin{align} x^{k+1} \gets x^{k} + 2^{-i^k} d^{k} \end{align} $$ $k \gets k+1$ として，Step1 へ戻る．
アルゴリズムに関するコメント：
- Step0におけるパラメータの推奨値に関する記述は少ないようだが，Luca et.al. (1996) では $\rho = 10^{-8},\, p=2.1$ としている．ただし，$\varsigma$ に関する記述は見つけることができなかった．また，終了条件に用いる $\varepsilon>0$ だが，十分小さな実数として $10^{-6}$ などを用いればよいだろう．
- Step1の終了条件に関して，Luca et.al. (1996) における数値実験では，$\| \min \{ x^k ,\, F(x^k) \} \|_2 \leqq 10^{-6}$ としているようである．
- Step2の $V_k$ はベクトル関数 $\Phi$ の点 $x^k$ における B-subdifferential (B 劣微分)と呼ばれ，$n\times n$ の正方行列となる．具体的な定義に関しては，[2] や Qi (1993) を参照されたいが，要するに Clarke 劣微分(subdifferential)*10をベクトル関数に拡張したようなものである．なお，B 劣微分 $\partial_B F(x)$ の凸包(convex hull)*11 $\mathrm{co} \, \partial_B F(x)$ の要素は一般化ヤコビ行列 (generalized Jacobian matrix) と呼ばれる．もし，ベクトル値関数 $F$ が微分可能であれば，B 劣微分 $\partial_B F(x)$ および 一般化ヤコビ行列はヤコビ行列 $\nabla F(x)$ に一致する．
アルゴリズム実装における注意点
アルゴリズムの実装に際して求める必要がある
- ベクトル値関数 $\Phi:\mathbb{R}^n\to\mathbb{R}^n$ の点 $x\in\mathbb{R}^n$ における B 劣微分 $\partial_B \Phi(x)$
- 実数値関数 $\Psi:\mathbb{R}^n\to\mathbb{R}$ の点 $x\in\mathbb{R}^n$ における勾配ベクトル $\nabla\Phi(x)$
の値を具体的に与える．
B 劣微分 $V \in \partial_B \Phi(x)$ の導出
結論から述べると，ベクトル値関数 $\Phi:\mathbb{R}^n\to\mathbb{R}^n$ の点 $x\in\mathbb{R}^n$ における B 劣微分 $\partial_B \Phi(x)$ は次のように計算することができる．
また，同じことだが，$[\partial \Phi_1(x),\ldots,\partial\Phi_n(x)]$ は対角行列を用いて次のように記述することもできる． $$ \begin{align} [\partial \Phi_1(x),\ldots,\partial\Phi_n(x)] = \mathrm{diag}[a_i(x)] + \nabla F(x) \mathrm{diag}[b_i(x)] \end{align} $$ ここに，$(a_i(x),b_i(x))$ は次のように与えられる． $$ \begin{align} (a_i(x), b_i(x)) = \left\{ \begin{array}{ll} \left( \xi_{i} - 1 ,\, \mu_{i} - 1 \right) & \mathrm{if} \quad (x_{i},F_{i}(x)) = (0,0) \\ \left( \frac{x_i}{\sqrt{x_i^{2} + F_i(x)^{2}}}-1 ,\, \frac{F_i(x)}{\sqrt{x_i^{2} + F_i(x)^{2}}}-1 \right) & \mathrm{if}\quad(x_i,F_i(x)) \neq (0,0) \end{array} \right. \end{align} $$ 上式における，$\xi_i$ および $\mu_i$ は $\xi_{i}^{2} + \mu_{i}^{2} \leqq 1$ を満たすような実数である．
さて，以下ではベクトル値関数 $\Phi:\mathbb{R}^n\to\mathbb{R}^n$ の点 $x\in\mathbb{R}^n$ における B 劣微分 $\partial_B \Phi(x) \in \mathbb{R}^{n\times n}$ が上のように記述される説明を行なう．アルゴリズムに関するコメントでも述べたように， B 劣微分 $\partial_B \Phi(x)$ に関する具体的な定義については文献を参照されたいが，以下の命題は以後で必要なので述べておく．
証明．Clarke (1990) [3] を参照のこと*14．
次に，局所 Lipschits 連続なベクトル値関数に対する Clarke 劣微分について述べた定理を挙げる．
さて，アルゴリズムの実装において必要なのは，ベクトル値関数 $\Phi : \mathbb{R}^n \to \mathbb{R}^n$ の点 $x \in \mathbb{R}^n$ における B 劣微分 $\partial_B \Phi(x)$ である．しかし，以上の命題および定理を用いることによって $$ \begin{align} \partial_B \Phi(x) \subseteq \mathrm{co}\, \partial_B &\Phi(x) \\ &\parallel \\ \partial &\Phi(x) \subseteq [\partial \Phi_1(x),\ldots,\partial\Phi_n(x)] \end{align} $$ が成立するので，結局のところ，$[\partial \Phi_1(x),\ldots,\partial\Phi_n(x)]$ を計算すればよいことになる． よって，各 $i$ に対する $\partial \Phi_i(x)$ の値を計算すればよいことになる．
なお，定理の仮定である局所 Lipschits 連続性だが，Fischer-Burmeister 関数 $\varphi$ は局所 Lipschits 連続であるので，関数 $\Phi$ もまた局所 Lipschits 連続となる．
勾配ベクトル $\nabla\Psi(x)$ の導出
次に，$\nabla \Psi (x)$ を導出する．こちらに関しては，具体的な計算によって以下のように表現される． $$ \begin{align} \nabla \Psi (x) &= \mathrm{diag}[a_i(x)] \Phi(x) + \nabla F(x) \mathrm{diag}[b_i (x)] \Phi(x) \\ &= \left( \mathrm{diag}[a_i(x)] + \nabla F(x) \mathrm{diag}[b_i (x)] \right) \Phi(x) \\ &= \left\{ \partial_B \Phi (x) \right\} \Phi(x) \end{align} $$ なお，$(a_i(x), b_i(x))$ は先に定義したものと同じである． 上式の第一式目は自明ではないのだが，簡単な計算によって成立が確認できる*15．また，$\nabla \Psi (x)$ を定義通りに書き下すと次のようになる： $$ \begin{align} \nabla \Psi (x) = \left( \frac{\partial}{\partial x_1} \Psi (x), \frac{\partial}{\partial x_2} \Psi (x), \ldots, \frac{\partial}{\partial x_n} \Psi (x) \right)^{\top} \end{align} $$ 線形相補性問題の場合，すなわち $F(x) = Mx+q$ の場合，各成分 $\frac{\partial}{\partial x_i} \Psi (x)$ は次のように与えられる： $$ \begin{align} \frac{\partial}{\partial x_i} \Psi(x) = a_i(x) \Phi_i(x) + \sum_{j=1}^{n} M_{ij} b_j(x) \Phi_j(x) \end{align} $$
Python による実装
理論編が非常に長くなってしまったが，以下に Python による実装例を示す．動作環境は Jupyter notebook で行なっている．なお，コードは線形相補性問題を解く内容となっているので注意されたい．少々長いので，関数の定義の部分と LCP を解く部分とで分けておく．また，numpyは以下のように定義済みであるとする．
import numpy as np
実装: 関数の定義
def FB(a, b):
return np.sqrt(a**2 + b**2) - a - b
def nablaF(x, M, q):
n = len(x)
A = np.zeros(n*n).reshape(n, n)
B = np.zeros(n*n).reshape(n, n)
for i in range(n):
if x[i]<=10**(-10) and np.dot(M[i],x)+q[i]<=10**(-10):
A[i][i] = 1/np.sqrt(2) - 1
B[i][i] = 1/np.sqrt(2) - 1
else:
A[i][i] = x[i] / np.sqrt(x[i]**2 + (np.dot(M[i], x) + q[i])**2) - 1
B[i][i] = (np.dot(M[i],x) + q[i]) / np.sqrt(x[i]**2 + (np.dot(M[i], x) + q[i])**2) - 1
ret = A + np.dot(M, B)
return ret
def Phi(x, M, q):
n = len(x)
ret = np.zeros(n)
for i in range(n):
ret[i] = FB(x[i], np.dot(M[i], x)+q[i])
return ret
def Psi(x, M, q):
ret = 0
n = len(x)
for i in range(n):
ret += FB(x[i], np.dot(M[i], x) + q[i])**2
return ret * 0.5
def nablaofPsi(x, M, q):
ret = np.dot(nablaF(x, M, q), Phi(x, M, q))
return ret
実装: 相補性問題を解くためのメインコード
def solve(M, q):
#step0
rho = 10**(-8)
p = 2.1
sigma = 10**(-3)
eps = 10**(-9)
np.random.seed(1)
n = len(q)
x = np.random.rand(n) #初期ベクトル
vecnabPhi = [] #グラフ用
vecnabPhi.append(np.linalg.norm(nablaofPsi(x, M, q)))
itermax = 100 #繰り返し回数の上限
for k in range(itermax):
#---step1---
nabPsi = nablaofPsi(x, M, q)
#終了判定
if np.linalg.norm(nabPsi) <= eps:
print('the optimum x is found!')
print('nablaPsi(x) = ', np.linalg.norm(nablaofPsi(x, M, q)))
return x, vecnabPhi
#---step2---
V = nablaF(x, M, q)
#Vが正則行列かチェック
det_V = np.linalg.det(V)
if det_V == 0:
d = -nablaofPsi(x, M, q)
else:
d = -np.dot(np.linalg.inv(V), Phi(x, M, q)) #本当はガウスの消去法を使うべき
if np.dot(nabPsi, d) > -rho*np.linalg.norm(d)**p:
d = -nablaofPsi(x, M, q)
#---step3---
ik = 0
while True:
term1 = Psi(x + d/(2**ik), M, q)
term2 = Psi(x, M, q)
term3 = np.dot(nablaofPsi(x, M, q), d) * sigma / (2**ik)
if term1 <= term2 + term3:
x = x + 2**(-ik) * d #暫定解の更新
vecnabPhi.append(np.linalg.norm(nablaofPsi(x, M, q))) #グラフ用
break
else:
ik += 1
if k%5 == 0: #5回ごとに出力
print('k = ',k, ' : ', Psi(x, M, q), ' : ', np.linalg.norm(nablaofPsi(x, M, q)))
return x,vecnabPhi #itermax以内に解が得られなかった場合
実装上のコメント：
- パラメータは，$\rho = 10^{-8}, \, p = 2.1 ,\, \varsigma = 10^{-3},\, \varepsilon = 10^{-9} $ とした．
- 収束判定は，$\| \nabla \Psi(x^k) \|_2 \leqq 10^{-9}$ とした．
- 初期点 $x_0 \in \mathbb{R}^n$ は，各要素が $[0,1]$ 上の一様乱数で生成した $n$ 次元ベクトルを用いた．
- B 劣微分の計算において，厳密に $(x_i,F_i(x)) = (0,0)$ が成り立つ場合は稀であると思われるので，$x_i \leqq 10^{-10}$ かつ $F_i(x) \leqq 10^{-10}$ が成り立つ場合には微分不可能であるとし，$\xi_i = \mu_i = \frac{1}{\sqrt{2}}$ とした．
数値実験
検証として，テスト問題に対して数値実験を行う．テスト問題には，次のような線形相補性問題を扱う． $$ \begin{align} \mathrm{Find}\quad &x \in \mathbb{R}^n \\ \mathrm{such\ that}\quad & x \geqq 0 ,\, Mx+q \geqq 0 ,\, \left< x, Mx+q \right> = 0 \end{align} $$ where $$ \begin{align} M = \begin{pmatrix} 4 & -1 & 0 & \cdots & 0 \\ -1 & 4 & -1 & \cdots & 0 \\ 0 & -1 & 4 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & -1 \\ 0 & \cdots & 0 & -1 & 4 \end{pmatrix} \end{align},\quad q = \begin{pmatrix} \sin \left( \frac{\pi}{2} \times 1 \right) \\ \sin \left( \frac{\pi}{2} \times 2 \right) \\ \vdots \\ \sin \left( \frac{\pi}{2} \times n \right) \end{pmatrix} $$
決定変数 $x\in\mathbb{R}^n$ の次元については任意に設定できるので，今回は $n=10$ で実験を行なった．
$10$ 回ごとの反復結果および最終的な出力結果を以下に示す．
x_opt, vecnabPhi = solve(M, q)
k = 0 : 2.4281593124919745 : 11.862373172208494
k = 10 : 7.198465316220975e-05 : 0.028843835531730776
k = 20 : 6.282334016218616e-07 : 0.004470012536269575
k = 30 : 2.796871979253169e-09 : 0.0002894708910221584
k = 40 : 5.2063123511471965e-12 : 1.1013249891757177e-05
k = 50 : 4.5779853670024516e-14 : 1.2492790866423237e-06
k = 60 : 1.8193030804420913e-16 : 7.250400110163342e-08
the optimum x is found!
nablaPsi(x) = 9.220883541289054e-10
さて，求まったベクトル $x$ および $Mx+q$ を見てみましょう．
print(x_opt)
array([6.29161933e-11, 7.14285715e-02, 2.85714286e-01, 7.14285717e-02,
9.53368225e-10, 7.14285717e-02, 2.85714286e-01, 7.14285715e-02,
9.80580150e-11, 9.43874375e-11])
print(np.dot(M, x_opt) + q)
array([ 9.28571429e-01, 3.84882160e-12, 0.00000000e+00, 5.70934662e-11,
8.57142860e-01, 5.70934199e-11, 2.22044605e-16, -2.89893656e-11,
9.28571429e-01, 2.79492347e-10])
非負制約 $x \geqq 0$ および $Mx+q \geqq 0$ が満たされていることが確認できる． 次に，条件: $\left< x, Mx+q \right> = 0$ が満たされているか見てみよう．
np.dot(x_opt, np.dot(M,x_opt) + q)
9.730093346493037e-10
...イマイチですね．やはり収束判定として $\| \nabla \Psi(x^k) \|_2$ を用いるのはよろしくないのだろうか．その辺りの考察については今回は行わないことにします．
いい感じでした．$\partial_B \Phi(x)$ を求めるプログラムが間違っていたことが原因だったようです．
収束過程の様子
折角なので，アルゴリズムの収束の様子を観察してみよう．具体的には，反復回数 $k$ を横軸に，収束条件として利用した $\| \nabla \Psi(x^k)\|_2$ の値を縦軸にとったグラフを作成する．グラフ作成用のコードは以下の通り*16：
import matplotlib.pyplot as plt
num = len(vecnabPhi) #データ数
ax_x = np.arange(num) #横軸(反復回数)
ax_y = vecnabPhi #縦軸(目的関数の勾配ベクトルの2ノルム)
plt.figure(figsize=(12, 10), dpi=50) #画像の大きさと解像度を指定
plt.yscale('log') #縦軸をlogスケールに
plt.grid(which='major',color='black',linestyle='-') #grid線を追加
plt.xlabel('$k$', fontsize=18)
plt.ylabel(r'$\|\|\nabla\Psi(x^k)\|\|_2$', fontsize=18)
plt.tick_params(labelsize=18)
plt.plot(ax_x,ax_y,color=(0,0,1),linewidth=1.5) #青でプロット
plt.savefig('./Image/exp1_n10.png') #画像の保存
出力されたグラフを以下に示す：
おわりに
少々専門的な内容になってしまったが如何だっただろうか．事前に Web 検索したところ，相補性問題に関する和文の記事はほとんど存在せず，あったとしても定式化や簡単な応用例の紹介くらいで，解法に関する内容を扱ったものは非常に少なかった．そういう意味では，このような内容を扱うことができたのは良かったと思っている．また，私は数理最適化に関してはそれなりに勉強してきたつもりだが，相補性問題に関してはほとんど無知な状態から初めて今回の記事を書いたので，起案してから投稿するまでにかなりの時間を要してしまった*17．今後ともより理解を深めていきたい．
Reference
[1] De Luca, Tecla, Francisco Facchinei, and Christian Kanzow. "A semismooth equation approach to the solution of nonlinear complementarity problems." Mathematical programming 75.3 (1996): 407-439.
[2] 福島雅夫. "非線形最適化の基礎." (2001).
[3] Clarke, Frank H. Optimization and nonsmooth analysis. Vol. 5. Siam, 1990.
[4] Cottle, Richard W., and George B. Dantzig. Complementary pivot theory of mathematical programming. No. TR-67-2. STANFORD UNIV CA OPERATIONS RESEARCH HOUSE, 1967.
[5] Kojima, Masakazu, et al. A unified approach to interior point algorithms for linear complementarity problems. Vol. 538. Springer Science & Business Media, 1991.
追記・修正情報
- Fischer-Burmeister 関数に関する部分を加筆．(191031)
- Reference の項目を追加．(191031)
- Reference の項目に，[1][2][3] を追加．(191031)
- 追記情報の項目を追加．(191031)
- 「線形相補性問題の場合におけるアルゴリズムに関するコメント」の項目を削除し，「B 劣微分 $\partial_B \Phi(x)$ と勾配ベクトル $\nabla\Phi(x)$ の導出」を追加．(191101)
- Reference に[4][5] を追加．(191101)
- 変分不等式問題，劣微分，局所 Lipschits 連続に関する脚注を追加．(191101)
- 目次(Contents)を追加．(191103)
- コードのデザインを修正．(191103)
- 命題，定理に対して囲みブロックを適用．囲みブロック内ではバックスラッシュ
\の個数が2つでよかったり，指数を表記するときに\でのエスケープが適用されなかったりするみたいなので注意．(191104) - $\nabla \Psi (x) =\left\{ \partial_B \Phi (x) \right\} \Phi(x)$ であることを追記．(191105)
- プログラムの間違いおよび出力結果を修正．takala 氏に感謝します． ．(191105)
*1:提示するアルゴリズム自体は一般の相補性問題に対しても適用できる．
*2:ベクトル関数 $F:\mathbb{R}^n\to\mathbb{R}^n$ に対する仮定として，連続的微分可能である条件を課しておく．
*3:2つの $n$ 次元ベクトル $x,\, y \in \mathbb{R}^n$ に対して，条件: $x\geqq 0,\, y\geqq 0,\, \left<x,\, y\right> = x^\top y = 0$ を満たすとき，$x,\, y$ は相補性条件を満たすという．相補性条件はしばしば $0 \leqq x \perp y \geqq 0$ のようにまとめて記述される．
*4:均衡問題として有名なものとしては，ナッシュ均衡問題("A procedure for finding Nash equilibria in bi-matrix games")や交通ネットワークの均衡問題(e.g., "NE/SQP: A robust algorithm for the nonlinear complementarity problem")，Spatial Price 均衡問題(e.g., "A solution condition for complementarity problems: with an application to spatial price equilibrium")などが挙げられる．
*5:折角素晴らしいスライドが存在するのでここに紹介しておきます．均衡問題の数理解析: 解の存在，一意性，安定性，分岐．
*6:その他の解法としては，Lemke 法などの不動点アルゴリズムをベースにしたもの[4]や二次計画問題に帰着させる手法[5]などが存在する．
*7:現在では，"Fischer-Burmeister" 関数と呼ばれることが多いようである．原著は"A special newton-type optimization method"．
*8:自明ではないとい思われる方もいらっしゃるかもしれないので述べておくと，任意の $x\in\mathbb{R}^n$ に対して $\Psi(x) \geqq 0$ であり，$x^{\star}$ が相補性問題の解であることと $\Psi(x^{\star})=0$ であることが等価であることより従う．
*9:すなわち，関係式 $\nabla \Psi (x^{k})^\top d^{k} > -\rho \|d^{k}\|^{p}$ が成り立つ場合である．
*10:そもそも劣微分を耳にしたことがない方も多いかもしれないが，劣微分とは微分不可能な凸関数に対して微分の概念を拡張したものである．劣微分の特徴として，一般の微分可能な関数に対する微分が関数で与えられるのに対して，劣微分は集合となる．数理計画用語集に，点 $x=1$ において微分不可能な関数 $f(x) = |x-1|$ に対する劣微分が例として挙げられている．また，凸関数に対する劣微分と区別するために，非凸関数に対する劣微分を Clarke 劣微分と呼ぶことがある．
*11:凸包とは，与えられた集合に対して，その集合を含む最小の凸集合と定義される．wikipedia にそれなりに詳しく書かれているので参照されたい．
*12:実数値関数 $f:\mathbb{R} ^n \to \mathbb{R}$ が局所 Lipschitz 関数であるとは，任意の有界な集合 $\Omega \subset \mathbb{R} ^n$ に対して，$|f(x) - f(y)| \leqq \kappa \|x-y\|_2 \; (x, y \in \Omega)$ を満たすような正の定数 $\kappa > 0$ が存在することをいう．
*13:相補性問題の $F$ との混同を避けるために，$\tilde{F}$ を用いている．
*14:文献 [2] でも証明は Clarke を参照している．
*15:簡単に確認できるのでぜひ手を動かしてみて欲しい．
*16:もっと上手い方法など御座いましたら教えてください．
*17:1週間くらいかかった...．
