2019-11-20

基本情報技術者試験に対する試験勉強

日記

はじめに

基本情報技術者試験に合格したので試験勉強として行ったことを書く．試験の概要などについては公式のページを参照されたい．

時系列

2019-8-14：受験申し込み(申し込みの期限日)
2019-9-10：試験勉強開始
2019-10-20：受験
2019-11-20：合格発表

試験勉強

基本的には，参考書を読むことと過去問を解いた．参考書は『キタミ式イラストIT塾基本情報技術者』(きたみりゅうじ)を使用した．参考書に関しては賛否両論あると思うが，書店で気に入ったものを購入すれば良いと思う．買わないという選択肢もあると思うが，体系的に纏まった書籍が一冊くらい手元にあると勉強が楽になると思う．過去問は IPA のページから入手可能．過去問の解説に関しては，基本情報技術者試験トットコムにほとんど載っている*1．午前の問題は 13 回分(平成24年春季〜平成31年春季)，午後の問題は 7 回分 (平成28年春季〜平成31年春季)をそれぞれ解いた．プログラミング言語は C 言語を選択した．言語選択に関しては色々な意見があると思うが，使用したことのある言語がある場合はそれらを選択すれば良いと思う．なければ，表計算を選択するのが無難だと思う．

おわりに

試験後の記事でも書いたように，午前の試験問題の半数程度は過去に出題された問題と全く同じものが出題される．そのことをもっと早い段階で理解していれば，参考書を隅々読むような無駄な時間を過ごさずに済んだかもしれない．また，午後の試験においては，時間を計測して解かないまま本番を迎えてしまったため，本番で時間が足りなくなるという惨事に見舞われてしまった．日々の勉強にもう少し緊張感を持って取り組むべきだった．また，午後の問題は選択式となっているため，早い段階でどの問題を選択するのかを決めて試験に臨むべきだった．

*1:ただし一部の午後の試験については解説がないので注意．

2019-11-10

リーマン多様体上の最適化の初歩と Pymanopt による数値実験

数学 Python 最適化

新しい記事を書きました．Pymanopt による実装は以下の記事を参照してください．

mirucacule.hatenablog.com

はじめに
予備知識
- リーマン多様体
- リーマン多様体の例
定式化
- 一般の最適化問題
- リーマン多様体上の最適化問題
リーマン多様体上の最適化問題の例
アルゴリズム
- $\R ^n$ における最適化アルゴリズム(直線探索)
- 多様体上の最適化アルゴリズム
実装
数値実験
- レイリー商最小化問題
  - 実装
  - 結果
- 低ランク半正定値行列の可解問題
  - 実装
  - 結果
おわりに
References

はじめに

この記事では，リーマン多様体上の最適化問題について扱う．まず，リーマン多様体に関する簡単な説明を行い，通常の最適化問題と比較する形でリーマン多様体上の最適化問題を定式化する．そして，Pymanoptと呼ばれるソルバーを用いて数値実験を行う．

予備知識

リーマン多様体

「リーマン多様体」とタイトルで使っているので，多様体の定義くらいは述べておきたいと思う．しかしながら，厳密に定義しようとすると，位相空間論の知識が必要になってしまうので，厳密性にこだわらない場合は，曲がった空間を抽象化したものくらいの認識でいいと思われる．さて，定義を述べる．

定義(リーマン多様体) $C ^{\infty}$ 級多様体 $\mathcal{M}$ 上の $2$ 次の対称テンソル場 $\omega$ が，$\mathcal{M}$ の各点 $p$ において正定値であるとき，$\omega$ を $\mathcal{M}$ 上のリーマン計量(Riemannian metric)という．また，リーマン計量 $g$ がひとつ与えられた多様体 $(\mathcal{M},g)$ のことをリーマン多様体(Riemannian manifold)という．

正直，この定義だけでリーマン多様体を理解することは不可能だと思う．詳しくは『多様体入門』(松島)や『多様体の基礎』(松本)などを参照されたい．

リーマン多様体の例

リーマン多様体の最も基本的な例として，$(n-1)$ 次元球面 $S ^{n-1} := \{ x\in\R ^n \mid x ^\top x = 1 \}$ がある．他の例としては次のような空間が挙げられる．

$n$ 次直交群 $\mathcal{O}(n) = \{ X \in \R ^{n \times n} \mid X ^{\top} X = I_n \} \subset \R ^{n \times n}$
シュティーフェル多様体 $\mathrm{St}(p,t) = \{ Y \in \R ^{n \times p} \mid Y ^{\top} Y = I_p \} \subset \R ^{n \times p}$
グラスマン多様体 $\mathrm{Grass}(p,n)$：$\R ^n$ の $p$ 次元部分空間全体

定式化

ここでは，一般の( $\R ^n$ 上の )最適化問題とリーマン多様体上の最適化問題の定式化を行う．

一般の最適化問題

一般の( $\R ^n$ 上の )最適化問題は次のように定式化される． $$ \begin{align} \mathrm{Minimize} \quad &f(x) \\ \mathrm{subject\ to} \quad &x \in S \end{align} $$ ここに，実数値関数 $f:\R ^n \to \R$ は目的関数を表す．また，集合 $S$ は $n$ 次元ユークリッド空間 $\R ^n$ の部分空間であり，特に $S=\R ^n$ の場合，最適化問題は無制約最適化問題と呼ばれる．

リーマン多様体上の最適化問題

リーマン多様体 $\mathcal{M}$ 上の最適化問題は次のように定式化される． $$ \begin{align} \mathrm{Minimize} \quad &f(x) \\ \mathrm{subject\ to} \quad &x \in \mathcal{M} \end{align} $$ ここに，実数値関数 $f:\mathcal{M} \to \R$ は目的関数を表す．一見，制約付き最適化問題にように見えるかもしれないが， $$ \begin{align} \underset{x\in \mathcal{M}}{\mathrm{Minimize}} \quad &f(x) \end{align} $$ と同じ意味であり，リーマン多様体上で考えれば制約なしと見なすことができる．

リーマン多様体上の最適化問題の例

ここでは，リーマン多様体上の最適化問題として定式化される問題の例を挙げる．

レイリー商最小化問題

レイリー商(Rayleigh quotient)最小化問題は次のように定式化される問題である． $$ \begin{align} \mathrm{Minimize} \quad &\tilde{f} (x) = \frac{x ^\top A x}{x ^\top x} \\ \mathrm{subject\ to} \quad &x \in \R ^n - \{ 0 \} \end{align} $$ ここに，$A$ は所与の $n$ 次元実対称行列を表す．制約条件で原点 $\{0\}$ を除いているため完全には無制約ではないが，原点から十分に離れたところで考えれば無制約と見なすことができる．関数 $\tilde{f}$ は任意の非零な実数 $\lambda$ に対して，次の性質を満たす． $$ \begin{align} \tilde{f}( \lambda x) &= \frac{(\lambda x) ^\top A (\lambda x)}{(\lambda x) ^\top (\lambda x)} \\ &= \frac{\lambda ^2 \cdot x ^\top A x}{\lambda ^2 \cdot x ^\top x} \\ &= \frac{x ^\top A x}{x ^\top x} \\ &= \tilde{f}(x) \end{align} $$ これより，関数 $\tilde{f}$ の値は引数 $x$ の方向のみに依存しており，ユークリッドノルムには依存していないことが分かる．したがって，$x$ のノルムを $1$ に固定したとしても，関数 $\tilde{f}$ を最小にする $x$ には影響がない．そこで，探索領域を $x ^\top x = 1$ 上に制限すると，上記最適化問題は次のように書き換えることができる． $$ \begin{align} \underset{x\in \R ^n}{\mathrm{Minimize}} \quad &f (x) = x ^\top A x \\ \mathrm{subject\ to} \quad &x ^\top x = 1 \end{align} $$ この問題は，通常の最適化の文脈では制約付き最適化問題となる．さて，$x ^\top x = 1$ を満たすようなベクトル $x$ の集合は，$(n-1)$ 次元球面 $S ^{n-1}$ で表されるので，上記の問題は $$ \begin{align} \mathrm{Minimize} \quad &f (x) = x ^\top A x \\ \mathrm{subject\ to} \quad &x \in S ^{n-1} \end{align} $$ と定式化することもできる．これは，$(n-1)$ 次元球面上の無制約最小化問題と見なすことができる．

所与の正方行列を近似する低ランク半正定値行列の可解問題

与えられた $n$ 次元実対称行列 $A$ に対して，$A$ を最も良く近似するようなランク $k<n$ の正定値行列 $S$ を求める問題を考える．ここで，行列 $A$ と半正定値行列 $S$ との近似度は次式で測られるとする． $$ \begin{align} L_{\delta} (S, A) := \sum_{i=1} ^{n} \sum_{j=1} ^{n} H_{\delta} (s_{i,j} - a_{i,j}) \end{align} $$ ここに，$\delta > 0$ は正の定数であり，$H_{\delta}(x) := \sqrt{x ^2 + \delta ^2} - \delta$ て定義される関数であり，"pseudo-Huber loss function" と呼ばれる．この問題を，$\R ^n$ 上で定式化すると次のようになる． $$ \begin{align} \underset{S\in \R ^{n\times n}}{\mathrm{Minimize}} \quad &L_{\delta} (S,A) \\ \mathrm{subject\ to} \quad &S \succcurlyeq O , \, \mathrm{rank}(S) = k \end{align} $$ さて，ランク $k$ の $n$ 次元半正定値行列全体の集合を $\mathcal{PSD}_{k} ^{n}$ で表す．すなわち，$\mathcal{PSD}_{k} ^{n} := \{ X \in \R ^{n \times n} \mid X \succcurlyeq O ,\, \mathrm{rank}(X) = k \}$ である．$\mathcal{PSD}_{k} ^{n}$ はリーマン多様体であるため，$\mathcal{PSD}_{k} ^{n}$ 上での最適化として定式化し直すと次のようになる． $$ \begin{align} \underset{S\in \mathcal{PSD}_{k} ^{n}}{\mathrm{Minimize}} \quad &L_{\delta} (S,A) \end{align} $$

その他の例

上記で挙げた例の他にも数多くの応用先があることが知られている．特に近年では，機械学習の分野における諸問題がリーマン多様体上の最適化問題として定式化することによって効率的に最適解を求解できることもあり，理論・応用の両面から盛んに研究されているようである．

特異値分解への応用
正準相関分析への応用
テンソル補完問題への応用

アルゴリズム

$\R ^n$ における最適化アルゴリズム(直線探索)

$\R ^n$ 上の無制約最適化問題に対するアルゴリズムを述べる．ここでは，最も基本的な反復法について述べる．

Step0 : 初期点 $x_0$ を設定し，$k=0$ とする
Step1 : 終了判定を行い，条件を満たすならば $x_k$ を最適解として出力する
Step2 : 探索方向 $\xi_k \in \R ^n$ とステップ幅 $\alpha_k > 0$ を決定する
Step3 : 暫定解を $x_{k+1} = x_k + \alpha_k \xi_k$ で更新し，Step1 へ戻る

多様体上の最適化アルゴリズム

$\R ^n$ 上の無制約最適化問題に対するアルゴリズムを拡張する形でリーマン多様体上の最適化アルゴリズムを構築する．ここで注意しなければならない点として，多様体上では Step3 における "$+$" が定義されていないということである．つまり，通常の意味での加算の結果，考えている多様体からはみ出してしまう場合があるということである．そこで，はみ出したところから多様体に引き戻す操作を行う必要があるわけだが，そのような概念を一般化したレトラクションと呼ばれる写像 $R$ を定義することができる．なお，本記事では具体的な探索方向ベクトルの求め方やレクトラクションの構成方法については触れない．

Step0 : 初期点 $x_0$ を設定し，$k=0$ とする
Step1: 終了判定を行い，条件を満たすならば $x_k$ を最適解として出力する
Step2 : 探索方向 $\xi_k \in \R ^n$ とステップ幅 $\alpha_k > 0$ を決定する
Step3 : 暫定解を $x_{k+1} = R_{x_k} (\alpha_k \xi_k)$ で更新し，Step1 へ戻る

実装

リーマン多様体上の最適化問題を Python で解くことを考える．本記事では自前で実装することはせず，Pymanoptという Python の toolbox を使って求解することを考える．

Pymanopt の概要

Pymanoptは多様体上での最適化を行うための Python toolbox であり，勾配やヘッセ行列を自動で計算することができる．ユーザは以下の3点を定義するだけで最適化を行うことができる．

最適化を行う多様体 $\mathcal{M}$ を定義する
目的関数 $f:\mathcal{M} \to \R$ を定義する
Pymanoptをインスタンス化する

より詳しい解説は公式サイトや文献を参照されたい．

Pymanopt の導入

Pymanoptはpipを使うことによって簡単にインストールできる．

pip install --user pymanopt

また，自動微分を行うために，目的関数をAutograd Theano TensorFlowのいずれかで定義する必要がある．これらのいずれにも馴染みのない場合，Autogradを使用することが推奨されている．AutogradはNumPyをラップしているため，NumPyに馴染みのある場合，簡単に実装できると思われる．なお，Autogradはpipによってインストールできる．

pip install autograd

実装

多様体の定義

Pymanoptでサポートされている多様体は以下の通りである．

Euclidean manifold：ユークリッド空間
Symmetric matrices：対称行列全体の集合
Sphere manifold：球面
Complex circle：複素平面上の円 $\{x \in \mathbb{C}^n \mid |x_1| = |x_2| =\dots = |x_n| = 1\}$
Rotations manifold：$\mathrm{SO}(n) := \{X\in\mathbb{R}^{n\times n} \mid X ^\top X=I_n, \mathrm{det}(X)=1\}$
Stiefel manifold：シュティーフェル多様体
Grassmann manifold：グラスマン多様体
Oblique manifold：Manifold of $m \times n$ matrices with unit-norm columns
Positive-definite matrices：正定値行列全体の集合
Elliptope manifold：Manifold of $n \times n$ positive-semidefinite matrices with rank $k$ and unit diagonal elements.
Fixed rank positive-semidefinite matrices：ランク $k$ の半正定値行列全体の集合
Fixed rank matrices：ランク $k$ の実対称行列全体の集合
Product manifold：積多様体

例えば，Sphere manifold を使用したい場合，以下のようにモジュールをインポートする．

from pymanopt.manifolds import Sphere

その上で，多様体を次のように定義する．例えば，$3$ 次元球面を定義したい場合には次のようにする．

manifold = Sphere(3)

他の多様体を用いる際の具体的な定義の仕方については公式サイトを参照されたい．

目的関数の定義

Autogradを使う場合，目的関数は次のように定義すればよい．

def cost(X):
    return np.exp(-np.sum(X**2))

多様体と目的関数の統合

多様体と目的関数は Pymanopt Problem オブジェクトによって結合し，その上でPymanoptソルバーに渡す必要がある．具体的には次のようにすればよい．

from pymanopt import Problem
problem = Problem(manifold=manifold, cost=cost)

ソルバー

以上で多様体の定義と目的関数の定義を終えた．後はソルバーに渡して解かせるだけである．ソルバーのパラメータは，ソルバーをインスタンス化したときするときに指定される．以下にデフォルト値を載せておく．

maxtime=1000：ソルバーが実行される最大時間(秒)
maxiter=1000：反復回数の最大回数
mingradnorm=1e-6：終了判定に用いられる勾配のノルム
minstepsize=1e-10：終了判定に用いられるステップサイズ
maxcostvals=5000：目的関数を評価する最大回数
logverbosity=0：ソルバーの動作中にログに記録される情報のレベル．0を指定すると最終的な最適解の情報のみを返す．2を指定すると最終的な最適解の情報とログ情報を保持する辞書を返す．

具体的な実装方法は次のように記述する．

xoptimal = solver.solve(problem, logverbosity=0)
xoptimal, optlog = solver.solve(problem, logverbosity=2)

Pymanoptに実装されているソルバーを以下に示す．これらのソルバーは Manopt をベースにしている．なお，ソルバーを自前で実装する場合，Solver abstract base class から継承する必要がある．

Trust regions()
SteepestDescent()
ConjugateGradient()
NelderMead()
ParticleSwarm()

例えば，Trust regions を用いる場合，次のように記述する．

from pymanopt.solvers import TrustRegions
solver = TrustRegions()
Xopt = solver.solve(problem)

数値実験

Pymanoptの説明が終わったので，Pymanoptを使って数値実験を行う．ここでは，上で例として挙げた問題を実験に用いる．以下では，必要なモジュールは次のようにインポート済みであるとする．

import autograd.numpy as np
from autograd import grad
from pymanopt import Problem
from pymanopt.solvers import TrustRegions

レイリー商最小化問題

$2$ 次対称行列 $A = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}$ におけるレイリー商最小化問題を扱う．実は，レイリー商最小化問題の最適化 $x ^{\star}$ は行列 $A$ の最小固有値に属する固有ベクトルであることが知られている．そのため，この問題に対する最適解は固有値および固有ベクトルを計算することによって解析的に計算することができる．そのため，アルゴリズムによって求まった最適化と比較することができる．なお，行列 $A = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}$ に対する固有値は，$\lambda_1 = 4 ,\, \lambda_2 = 2$ であり，$\lambda_2 = 2$ に属する固有ベクトルは，$t \left( -\frac{1}{\sqrt{2}} , \, \frac{1}{\sqrt{2}} \right) ^{\top} \; (\forall t \in \R)$ となる．

実装

from pymanopt.manifolds import Sphere
A = np.array([[3,1],[1,2]])
manifold = Sphere(2)
def cost(x):
    return np.dot(x, np.dot(A, x))
problem = Problem(manifold=manifold, cost=cost)
solver = TrustRegions(logverbosity=2)

結果

実行結果は以下のようになる．

xopt, optlog = solver.solve(problem)
Compiling cost function...
Computing gradient of cost function...
Computing Hessian of cost function...
Optimizing...
                                            f: +2.637939e+00   |grad|: 1.864309e+00
acc TR+   k:     1     num_inner:     0     f: +2.100392e+00   |grad|: 8.733960e-01   exceeded trust region
acc       k:     2     num_inner:     0     f: +2.000295e+00   |grad|: 4.859006e-02   maximum inner iterations
acc       k:     3     num_inner:     0     f: +2.000000e+00   |grad|: 7.173227e-06   maximum inner iterations
acc       k:     4     num_inner:     0     f: +2.000000e+00   |grad|: 4.440892e-16   maximum inner iterations
Terminated - min grad norm reached after 4 iterations, 0.02 seconds.

optlog の中身は次のように確認できる．

print(optlog)
{'solver': 'TrustRegions',
 'stoppingcriteria': {'maxtime': 1000,
  'maxiter': 1000,
  'mingradnorm': 1e-06,
  'minstepsize': 1e-10,
  'maxcostevals': 5000},
 'solverparams': None,
 'stoppingreason': 'Terminated - min grad norm reached after 4 iterations, 0.02 seconds.',
 'final_values': {'x': array([-0.70710678,  0.70710678]),
  'f(x)': 2.0,
  'time': 0.024547100067138672,
  'stepsize': inf,
  'gradnorm': 4.440892098500626e-16,
  'iterations': 4}}

final_values が最適解に対応しており，$x ^{\star} = \begin{pmatrix} -0.70710678 \\ 0.70710678 \end{pmatrix}$ であることが分かる．また，f(x) は最適値に対応しており，$f(x ^{\star}) = 2.0$ であることが分かる．これらは，確かに行列 $A$ の最小固有値とそれに属する固有ベクトルであることが確認できる．

低ランク半正定値行列の可解問題

例として，$3$ 次対称行列 $A = \begin{pmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \\ 3 & 5 & 6 \end{pmatrix}$ を最も良く近似するような正定値行列 $S$ を求める問題を扱う．なお，$S$ のランクは $2$ となるようにする．

実装

from pymanopt.manifolds import PSDFixedRank
A = np.array([[1,2,3],[2,4,5],[3,5,6]])
manifold = PSDFixedRank(A.shape[0], A.shape[0]-1)
def cost(Y):
    S = np.dot(Y, Y.T)
    delta = 0.5
    return np.sum(np.sqrt((S-A)**2 + delta**2) - delta)
problem = Problem(manifold=manifold, cost=cost)
solver = TrustRegions(logverbosity=2)

結果

実行結果は以下のようになる．

Yopt, optlog = solver.solve(problem)
Compiling cost function...
Computing gradient of cost function...
Computing Hessian of cost function...
Optimizing...
                                            f: +2.401402e+01   |grad|: 5.877665e+00
acc TR+   k:     1     num_inner:     0     f: +9.994776e+00   |grad|: 9.898422e+00   negative curvature
REJ TR-   k:     2     num_inner:     0     f: +9.994776e+00   |grad|: 9.898422e+00   exceeded trust region
acc TR+   k:     3     num_inner:     0     f: +3.883474e+00   |grad|: 1.009627e+01   exceeded trust region
REJ TR-   k:     4     num_inner:     1     f: +3.883474e+00   |grad|: 1.009627e+01   exceeded trust region
acc TR+   k:     5     num_inner:     0     f: +1.049904e+00   |grad|: 7.140704e+00   exceeded trust region
REJ TR-   k:     6     num_inner:     3     f: +1.049904e+00   |grad|: 7.140704e+00   reached target residual-kappa (linear)
acc       k:     7     num_inner:     0     f: +4.655802e-01   |grad|: 4.194963e+00   exceeded trust region
acc TR+   k:     8     num_inner:     3     f: +2.900976e-01   |grad|: 1.679291e+00   exceeded trust region
acc       k:     9     num_inner:     3     f: +2.542600e-01   |grad|: 6.877936e-02   reached target residual-kappa (linear)
acc       k:    10     num_inner:     1     f: +2.542176e-01   |grad|: 4.442510e-03   reached target residual-theta (superlinear)
acc       k:    11     num_inner:     4     f: +2.542130e-01   |grad|: 1.874113e-05   reached target residual-theta (superlinear)
acc       k:    12     num_inner:     4     f: +2.542130e-01   |grad|: 6.544584e-10   reached target residual-theta (superlinear)
Terminated - min grad norm reached after 12 iterations, 0.43 seconds.

optlog の中身は次のように確認できる．

print(optlog)
{'solver': 'TrustRegions',
 'stoppingcriteria': {'maxtime': 1000,
  'maxiter': 1000,
  'mingradnorm': 1e-06,
  'minstepsize': 1e-10,
  'maxcostevals': 5000},
 'solverparams': None,
 'stoppingreason': 'Terminated - min grad norm reached after 12 iterations, 0.13 seconds.',
 'final_values': {'x': array([[0.77928005, 0.82330077],
         [0.74491914, 1.87118565],
         [1.37193605, 2.07299432]]),
  'f(x)': 0.25421300993652163,
  'time': 0.13328790664672852,
  'stepsize': inf,
  'gradnorm': 2.2973372294669235e-11,
  'iterations': 12}}

最適解はfinal_valuesを参照すると， $$ \begin{align} Y_{\star} = \begin{pmatrix} 0.77928005 & 0.82330077 \\ 0.74491914 & 1.87118565 \\ 1.37193605 & 2.07299432 \end{pmatrix} \end{align} $$ であることが確認できる．それでは，最適解における $Y_{\star} Y_{\star} ^{\top}$ を見てみよう．

print(np.dot(Yopt, Yopt.T))
array([[1.28510159, 2.12104918, 2.77582023],
       [2.12104918, 4.05624029, 4.90093862],
       [2.77582023, 4.90093862, 6.17951397]])

いかがだろうか．与えられた行列 $A = \begin{pmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \\ 3 & 5 & 6 \end{pmatrix}$ と見比べると確かに似通っていることが見て取れる．では，実際にこの行列のランクが $2$ になっているのか確認してみよう．

S = np.dot(Yopt, Yopt.T)
print(np.linalg.matrix_rank(S))
2

numpy.linalg.matrix_rankは特異値分解をベースに与えられた行列のランクを計算してくれる関数である．この結果より，確かに得られた行列のランクは $2$ になっていることが分かる．

おわりに

今回はリーマン多様体上の最適化を簡単に説明した上で，Pymanoptを導入し，ソルバーを用いて数値実験を行なった．ほとんどソルバーに丸投げしたため，多様体上での最適化アルゴリズムを構築する理論的な解説を行うことができなかったが，自分自身も理解が不十分であるためご容赦いただければ幸いである．自前で実装する例に関しては，ksknw 様のリーマン多様体の最適化(シュティーフェル多様体上の最適化まで)を参照されたい．本記事では詳しく触れなかったシュティーフェル多様体上の最適化やリーマン多様体上での勾配について説明されている．また，本格的にこの分野を勉強したい場合は，[2] を参照されたい．リーマン多様体上の最適化に関して体系的にまとめられている数少ない書籍である．

References

[1] Townsend, James, Niklas Koep, and Sebastian Weichwald. "Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation." The Journal of Machine Learning Research 17.1 (2016): 4755-4759.
[2] Absil, P-A., Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2009.

2019-11-06

目次の生成・コードの書き方・囲みブロック【はてなブログ】

ブログ整備

はじめに
目次の生成
コードの挿入
囲みブロック
おわりに
参考

はじめに

今回はブログ整備のお話である．はてなブログで記事を書くにあたって，デフォルトで設定されているデザインだけでは見やすさの観点からイマイチである場合がある．そこで今回は，コードを(比較的)キレイに載せる方法と囲みブロックを適用する方法について述べる．ついでに，目次を記事中に載せる方法についても述べる．

目次の生成

目次の生成に関しては，[:contents]を目次を載せたい部分に書けば O.K. である．これほど簡単なことはない．目次に採用されるのは，Markdown では，ハッシュタグ#が3つからである．決して，#を2つ以下で記事を作らないように．目次のデザインを変更することや，目次を番号付けで表示させることもできるがここでは割愛する．参考で挙げたページを参照されたい．

コードの挿入

コードに関しては，バッククォートを3つ用いて記述することが一般的だと思われる*1．

#!/bin/python3
print("Hello, World!")

~~簡単ではあるが，あまり美しくない．そこで代替策を考える．~~（デザインCSSを修正することによって，“いい感じ”に表示させるようにしました．）

導入方法

導入の仕方は以下の通り：

「管理画面」→「デザイン設定」→「フッタ」の順に推移し，HTML を記述する欄に以下のコードを挿入する．
- 以下のコードはvsというテーマを使用した場合の記述例である．
- 他のテーマとしては，githubやatom-one-darkなどがあり，highlight.jsで確認できる．

<link rel="stylesheet"
     href="//cdnjs.cloudflare.com/ajax/libs/highlight.js/9.13.1/styles/vs.min.css">
<script src="//cdnjs.cloudflare.com/ajax/libs/highlight.js/9.13.1/highlight.min.js"></script>
<script>hljs.initHighlightingOnLoad();</script>

使い方

preタグとcodeタグで記述したいコードを次のように囲めばよい．

<pre style="background-color: rgb(252, 252, 252) !important;">
<code class="html">《記述したいコード》
</code></pre>

例えば，先ほどのPythonのコードは次のようになる．

#!/bin/python3
print("Hello, World!")

コードの中に< > &などがある場合，HTML Escape / Unescapeなどでエスケープする
<code>タグで改行すると，一行目が空白になってしまうため，<code>タグの直後からコードを書く必要がある．

囲みブロック

例えば，次のように定理や命題をブロックで囲みたい場合がある．

定理(留数定理) 複素関数 $f(z)$ の $n$ 個の極 $a_k \; (k=1,2,\dots,n)$ を囲む閉曲線 $C$ について， $$ \begin{align} \oint_{C} f(z) \mathrm{d}z = 2 \pi i \sum_{k} \mathrm{Res} (a_k, f) \end{align} $$が成り立つ．ここに，$\mathrm{Res}(a,f)$ は複素関数 $f$ の極 $z=a$ における留数を表す．

導入方法

「管理画面」→「デザイン設定」→「デザインCSS」の順に推移し，CSS を記述する欄に以下のコードを挿入する．

.theorem {
 position: relative;
 margin: 2em auto;
 padding: 1.2em;
 color: #555555; /* 文章色 */
 background-color: #f7f7f7; /* 背景色 */
 border: 1px solid #555555; /* 枠線の太さ・色 */
 width: 90%;
}
.title-theorem {
 position: absolute;
 padding: 0 .5em;
 left: 20px;
 top: -15px;
 font-weight: bold;
 background-color: #fff; /* タイトル背景色 */
 color: #555555; /* タイトル文字色 */
}

使い方

ブロックで囲みたい定理等を次のようにdivタグで囲めばよい．例えば，先ほどの留数定理は次のように記述する．

<div class="theorem"><span class="title-theorem">定理(留数定理)</span>
複素関数 $f(z)$ の $n$ 個の極 $a_k \; (k=1,2,\dots,n)$ を囲む閉曲線 $C$ について，
$$
\begin{align}
\oint_{C} f(z) \mathrm{d}z = 2 \pi i \sum_{k} \mathrm{Res} (a_k, f)
\end{align}
$$
が成り立つ．ここに，$\mathrm{Res}(a,f)$ は複素関数 $f$ の極 $z=a$ における留数を表す．
</div>

ブロックの中では，数式の改行が\2つでよかったりするなど，特有の注意が必要なようである*2．

おわりに

はてなブログにおける

目次の生成の仕方
コードの(比較的)キレイな載せ方
囲みブロックの使い方

について扱った．先行記事が多く存在するため，この記事を書くべきかどうか迷ったのだが，自分用の備忘録して書いておいた．

参考

はてなブログ「目次」の見た目をCSSで変える（コピペOK）
- 目次を番号付けしたい方はこちらを参考にするとよい．他にも同様の記事が多くあるが，一番よかった．
はてなブログのシンタックスハイライトを綺麗にする
- コードの書き方に関しては，こちらのブログを参考にした．コード内の行間の調整についても扱っているので，行間を調整したいという方は参照されたい．
【CSS】おしゃれなボックスデザイン（囲み枠）のサンプル30
- Web 検索で「CSS 囲みブロック」などと検索すれば，同様の記事が多く見つかると思う．今回はトップに出てきたこちらのサイトを参考にした．

*1:コード例を載せたかったのだが，バッククォートをエスケープする方法がよく分からなかったので割愛．

*2:その他の注意点については今後追記するかもしれない．

2019-10-28

相補性問題の初歩と Python による数値解法

Python 最適化数学

はじめに
相補性問題の定式化
相補性問題の再定式化
相補性問題に対するアルゴリズム
アルゴリズム実装における注意点
- B 劣微分 $V \in \partial_B \Phi(x)$ の導出
- 勾配ベクトル $\nabla\Psi(x)$ の導出
Python による実装
- 実装: 関数の定義
- 実装: 相補性問題を解くためのメインコード
数値実験
- 収束過程の様子
おわりに
Reference
追記・修正情報

はじめに

本記事の目的は，数理計画の分野において重要な役割を果たす相補性問題の基本について解説し，最も基本となる線形相補性問題に対するアルゴリズムを提示*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)$ は次のように計算することができる．

定理ベクトル値関数 $\Phi:\mathbb{R} ^n\to\mathbb{R} ^n$ の点 $x\in\mathbb{R} ^n$ における B 劣微分 $\partial_B \Phi(x) \in \mathbb{R} ^{n\times n}$ は以下のように与えられる． $$ \begin{align} \partial_B \Phi(x) \subseteqq [\partial \Phi_1(x),\ldots,\partial\Phi_n(x)] \end{align} $$ ここで，上式の右辺は，$\partial \Phi_i(x) \in \mathbb{R} ^n$ に属するベクトルを第 $i$ 列とするような $n\times n$ 行列全体を表し，各 $i$ に対する $\partial \Phi_i(x)$ は次のように与えられる． $$ \begin{eqnarray} \partial \Phi_i (x) = \begin{cases} \left\{ (\xi_{i} - 1)e ^{i} + (\mu_{i} - 1)\nabla F_{i}(x) \mid \xi_{i} ^{2} + \mu_{i} ^{2} \leqq 1 \right\} & \mathrm{if} \quad (x_{i},F_{i}(x)) = (0,0) \\ \left\{ \left(\frac{x_i}{\sqrt{x_i ^{2} + F_i(x) ^{2}}}-1\right)e ^{i} + \left(\frac{F_i(x)}{\sqrt{x_i ^{2} + F_i(x) ^{2}}}-1\right)\nabla F_i(x) \right\} & \mathrm{if}\quad(x_i,F_i(x)) \neq (0,0) \end{cases} \end{eqnarray} $$ ここに，$e ^i$ は第 $i$ 成分のみ $1$ でそれ以外の成分は $0$ であるような単位ベクトルを表す．

また，同じことだが，$[\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)$ に関する具体的な定義については文献を参照されたいが，以下の命題は以後で必要なので述べておく．

命題ベクトル値関数 $\tilde{F}:\mathbb{R} ^{n}\to\mathbb{R} ^{m}$ を局所 Lipschitz 関数*12とする．このとき，$\tilde{F}$ の点 $x$ における Clarke 劣微分 $\partial\tilde{F}(x)$ は，$\tilde{F}$ の B 劣微分 $\partial_B \tilde{F}(x)$ を用いて $$ \partial\tilde{F}(x) := \mathrm{co}\,\partial_B \tilde{F}(x) $$ のように定義される*13．

証明．Clarke (1990) [3] を参照のこと*14．

次に，局所 Lipschits 連続なベクトル値関数に対する Clarke 劣微分について述べた定理を挙げる．

定理(定理2.64[2]) 局所 Lipschits 連続なベクトル値関数 $\tilde{F}:\mathbb{R} ^{n}\to\mathbb{R} ^{m}$ とその成分である実数値関数 $\tilde{F}_{i} : \mathbb{R} ^{n}\to\mathbb{R} \; (i=1,2,\dots,m)$ の Clarke 劣微分のあいだに次の関係が成立する． $$ \begin{align} \partial \tilde{F}(x) \subseteq [\partial \tilde{F}_{1}(x),\partial \tilde{F}_{2}(x),\dots,\partial \tilde{F}_{m}(x)] \end{align} $$ ここで，上式における右辺は，$\partial \tilde{F}_{i}(x)$ に属するベクトルを第 $i$ 列とする $n \times {m}$ 行列全体の集合を表す．

さて，アルゴリズムの実装において必要なのは，ベクトル値関数 $\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') #画像の保存

出力されたグラフを以下に示す：

f:id:mirucacule:20191105200304p:plain — 収束過程

おわりに

少々専門的な内容になってしまったが如何だっただろうか．事前に 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週間くらいかかった...．

2019-10-25

【証明】半正定値行列のトレースの非負性と性質

数学

はじめに
半正定値行列のトレースの非負性
- 証明
半正定値行列のトレースに関する性質
- 証明
おわりに

はじめに

本記事では，以下の二つの命題に対する証明について述べる：

半正定値行列のトレースは非負である．
半正定値行列のトレースが $0$ であるならば，その行列は零行列である．

予備知識としては簡単な線型代数の基礎程度である．

半正定値行列のトレースの非負性

命題正方行列 $A\in\R^{n\times n}$ を半正定値行列，すなわち，任意のベクトル $d\in \mathbb{R}^n$ に対して： $$ \begin{align} \left< d,\, Ad \right> = d^\top A d \geqq 0 \end{align} $$ が成り立つとする*1．このとき，行列 $A$ のトレースは非負である，すなわち， $$ \begin{align} \mathrm{tr} (A) = \sum_{i=1} ^{n} A_{ii} \geqq 0 \end{align} $$ が成り立つ．

証明

証明行列のトレースの性質より，$A$ の固有値を $\lambda_1,\,\lambda_2,\dots,\lambda_n$ とすると， $$ \begin{align} \mathrm{tr} (A) = \sum_{i=1} ^{n} \lambda_i \end{align} $$ が成り立つ*2。また，行列 $A$ は半正定値であるので，すべての固有値は非負である*3。したがって，$\mathrm{tr}(A) \geqq 0$である．

半正定値行列のトレースに関する性質

命題半正定値行列 $A\in\mathbb{R}^{n\times n}$ のトレースが $0$ とする．このとき，行列 $A$ は零行列である．

証明

証明行列 $A$ は半正定値行列なので，$A = B^\top B$ を満たすような行列 $B \in \mathbb{R}^{n\times m}$ が存在する*4．よって， $$ \begin{align} \mathrm{tr} (A) = \mathrm{tr} (B^\top B) \end{align} $$ が自然に成り立つ．命題の条件より，$\mathrm{tr}(A)=0$ なので $\mathrm{tr} (B^\top B) = 0$ である．さて，$n$ 次正方行列 $B^\top B$ の対角成分(第 $(i,i)$ 成分)は次のように与えられる： $$ \begin{align} (B^\top B)_{ii} = B_{i1} ^2 + B_{i2} ^2 + \dots + B_{in} ^2 \end{align} $$ したがって，$\mathrm{tr} (B^\top B)$ は次のように表すことができる： $$ \begin{align} \mathrm{tr} (B^\top B) &= \sum_{i=1} ^{n} (B^\top B)_{ii} \\\\ &= \sum_{i=1} ^{n} \sum_{j=1} ^{n} B_{ij}^2 \end{align} $$ 上式と $\mathrm{tr} (B^\top B) = 0$ を合わせて，$B_{ij} = 0$ を得る．したがって，$B$ は零行列である．いま，$A=B^\top B$ であるので $A$ も零行列である．

おわりに

半正定値行列に関する性質について扱った．半正定値行列は数理最適化の分野で頻出する*5ので，その性質について熟知しておくと何かと便利である．

*1:このような二次の項のみからなる多項式のことを「二次形式」という．

*2:証明は，高校数学の美しい物語を参照されたい．

*3:半正定値行列 $A$ の二次形式 $\left< d,\, Ad \right>$ が非負であることと，半正定値行列 $A$ の固有値が非負であることは同値な条件である．これに関しても、高校数学の美しい物語に証明が載っている．

*4:これは決して自明ではない．理系アラカルトに証明が載っているので参照されたい．

*5:例えば，制御の分野などの諸問題を定式化する上で半正定値行列を変数とした最適化問題は重要であり，ソルバーで効率的に解くことができる．