数理計画法入門～Markowitzの最小分散ポートフォリオの計算まで
この記事の立ち位置
この記事は数理計画法についてより詳細な教科書や資料で学ぶ前の、簡単なまとめである。計算や証明の詳細には立ち入らないが、理解する上での要点を記述し、円滑な理解の一助となるべく執筆した。必ずしも全ての記述、用語が正確であるとの保証はしないが、テキストを読む上での「気持ち」を十分に捉えたものだと自負している。
数理計画法を学ぶテキストとしては、東京大学工学教程の『最適化と変分法』が良い。
但し、このテキストは補足資料として同シリーズの『線形代数Ⅱ』を参照し、整数計画問題は『離散数学』の方に記載されている為、可能であればそれらも購入するのが望ましい。『離散数学』を読む上で『代数学』も準備すべきである。
ウェブ上で学べる資料としては、東京工業大学の水野研究室で公開されているpdfが役立つ。
数理計画法の定義
数理計画法に限らず広く一般に、最適化問題とは、取り得る選択肢の中からどれが一番良いか、何らかの基準を決めた上で決定する問題のことである。その定式化には様々な形が考えられるが、一般的な形で書くと例えば次のようになる。
$$
\begin{rcases}
\text{Minimize} & f(x)\\
\text{Subject to} & x \in S
\end{rcases}
$$
ここで、$${S}$$は実行可能領域、$${f(x)}$$は目的関数であり、実行可能領域の中から目的関数を最小化する解$${x^*}$$を見付けるというのが目標である。
ここで$${x}$$や$${S}$$としてどのようなものが考えられるか、本来は数値に限らず自由である。
具体例として対戦に使うポケモンのパーティの構築を考えてみる。1匹のポケモンについて対戦考察に必要な情報が種族名、性格、特性、覚えている技4つ、6つのパラメータそれぞれに対する努力値と個体値、持たせる道具で表現出来るものとする。この時、1匹のポケモンは8個の文字列(種族名、性格、特性、技4つ、道具)と12個の整数(努力値と個体値が6個ずつ)で表される。そして6匹のポケモンで1つのパーティが構成されるのであるから、48個の文字列と72個の整数、合計120成分のベクトル$${\bm{x}}$$によって1つのパーティを表すことが出来る。
しかし、ベクトルの成分として任意の文字列や整数が許容される訳ではない。例えばポケモンの種族名の成分として認められるのは「ピカチュウ」や「ヒトカゲ」であって、「けつばん」なり「繧ｫ繧､繝ｪ繝･繝ｼ」は許されないし、各パラメータの努力値に1023等の通常取り得ない値を入れて滅茶苦茶強いポケモンを作っては改造になってしまう。このように、ベクトルの成分の値が適切か判断するというのが$${S}$$の役割であり、適切なベクトル$${\bm{x}}$$は$${S}$$に含まれ、そうでないものは含まれないとする。
そしてベクトル$${\bm{x}}$$が与えられた時、自分が関心を持っているものを目的関数$${f(\bm{x})}$$として設定する。例えば強いパーティが欲しいのであれば勝率、或いは単位時間当たりに稼げるレートである。但し、最適化問題は最小化を考える為、今回のように目標とする数値が高い程望ましい場合には、マイナス符号を付ける必要がある(最適化問題で基本的に最小化を考えるのは、主問題が最小化、双対問題が最大化という形式に統一することで、どちらがどちらか分かり易くする為であろう。双対問題が何なのかは後述)。
しかし、上に具体例として挙げた問題は解き難い。ベクトルの成分として文字列が入っており、また目的関数もどのように計算すれば妥当なのか判然としないからである。ゲームのプレイヤーは感覚的にベクトルの更新を行い、その場合の勝率を体感的には把握しているが、それをいざ定量的かつ他人に分かる形で言語化するのは不可能であろう。このような問題はAlpha GOの開発のように、コンピュータ内で大量に対戦を繰り返して機械学習等によって最強を探すということは可能でも、最適なベクトルを解析的に求めるというのには向いていない。
そこで、数理計画法においては、ベクトルの成分として数のみ許容し、目的関数$${f(\bm{x})}$$も実行可能領域$${S}$$も数式によって書ける問題を考え、解は経験的に発見するのではなく数学的に最良のものを導出するものとする。実行可能領域は多数の方程式と不等式で構成する。
数理計画法の分類
とはいえ、数式で全て表現された最適化問題だからといって、それだけの情報で最適解を上手く計算出来る訳ではない。それは中学や高校の数学の授業において「方程式を解け」と言われた時、その方程式の種類が何なのか――1次方程式なのか2次方程式なのか、或いは更に高次であったり、三角関数や指数、対数関数を含むのか――等によって解法を変えるのと同じである。広い範囲の問題を解ける汎用的な手法があるにしても、簡単な問題であると分かっているならば簡単な手法を用いた方が効率が良いのは自明である。数理計画法では手計算では扱えないような規模の問題をコンピュータで解くのが普通である(例えば、先に挙げたポケモンの例は120次元のベクトルを考えており、これは明らかに手計算の範疇を超えている！)ので、計算効率は無視出来ない重要な問題と言えよう。これ故に、「○○問題」「○○計画法」という語が大量に生じているのである。
目的関数も、制約条件の等式も不等式も1次式で与えられる問題を線形計画問題という。即ち
$$
\begin{rcases}
\text{Minimize} & \bm{c}^{\rm T} \bm{x} \\
\text{Subject to} & A \bm{x} = \bm{b} \\
&A' \bm{x} \geq \bm{b}'
\end{rcases}
$$
と表記された問題のことである。尚、ベクトルの不等号は一方のベクトルの各成分がもう一方の対応する成分より値が大きいことを意味する。但し、この表記のまま使われることは少ない。基本的に、非負の変数(スラック変数)を幾つか使って式を書き直し、不等号条件がベクトルの各成分が非負という単純なものとなるようにする。
$$
\begin{rcases}
\text{Minimize} & \bm{c}^{\rm T} \bm{x} \\
\text{Subject to} & A \bm{x} = \bm{b} \\
&\bm{x} \geq \bm{0}
\end{rcases}
$$
これを等式標準形という。等式標準形を作る利点は、ペナルティ関数やバリア関数を目的関数に足し合わせることで、不等号条件を実質的に消せることであり、逐次2次計画法、主双対内点法、拡張ラグランジュ関数(乗数)法等の実用的なアルゴリズムの一部として組み込まれている。不等号条件をスラック変数で簡略化するのは線形計画問題に限らない一般の数理計画問題で行われる。
また、今まで述べてきた最適化問題(主問題と呼ばれる)に対し、対応する双対問題を同時に考えることがある。双対問題は主問題からある手続きによって生成される最適化問題で、双対問題の双対問題は主問題となり、両者の最適解の値が一致するものである(常にこのような性質が成り立つ訳ではなく、最適解の値に差(双対ギャップと呼ばれる)が存在することもあるが、線形計画問題のような性質の良い問題を考えている場合には気にしなくて良い。双対問題を作る手続きも複数あり、手続きによって異なる双対問題を作ることもあり得るが、これも気にしなくて良い)。主双対内点法等のアルゴリズムは主問題だけでなく双対問題も同時に計算し最適解を導出するので、双対問題という概念は知っておく必要がある。
等式標準形の線形計画問題の場合、対応する双対問題は以下のようになる。
$$
\begin{rcases}
\text{Maximize} & \bm{b}^{\rm T} \bm{y} \\
\text{Subject to} & A^{\rm T} \bm{y} + \bm{s} = \bm{c} \\
&\bm{s} \geq \bm{0}
\end{rcases}
$$
線形計画問題で表現出来る問題は多いが、当然、それだけで全ての最適化問題を表現出来る筈もない。線形計画問題以外の問題設定としては、目的関数や実行可能領域は同じく線形関数で表現されるが解として整数値しか許さない整数計画問題、実行可能領域は線形関数のままだが目的関数に2次式を許した2次計画問題等がある。等式標準形の2次計画問題の主問題と双対問題はそれぞれ次のように書かれる。尚、$${\bm{z}}$$は双対問題のスラック変数である。
$$
\begin{rcases}
\text{Minimize} & \frac 1 2 \bm{x}^{\rm T}Q\bm{x} + \bm{c}^{\rm T} \bm{x} \\
\text{Subject to} & A \bm{x} = \bm{b} \\
&\bm{x} \geq \bm{0}
\end{rcases}
$$
$$
\begin{rcases}
\text{Maximize} & \bm{b}^{\rm T} \bm{y} - \frac 1 2 \bm{x}^{\rm T}Q\bm{x} \\
\text{Subject to} & A^{\rm T} \bm{y} + \bm{z} = Q \bm{x} + \bm{c} \\
&\bm{z} \geq \bm{0}
\end{rcases}
$$
その他、様々な分類が(特に目的関数が凸関数、実行可能領域が凸集合となる凸計画問題で主に)なされているが、それらについては省略する。
しかしながら、使われている数式によって問題を分類しても、まだ解けるようにならない。というより、目的関数は方程式ではなく関数であるから、具体的なベクトル$${\bm{x}}$$を代入してその値を計算することは出来ても、それを最小化する特定の$${\bm{x}^*}$$を返すという機能が無いのである。そこで、$${f(\bm{x})}$$を最小化する$${\bm{x}^*}$$を解に持つような方程式を、元の最適化問題を変形することによって作りだす必要がある。ラグランジュの未定乗数法やKKT条件とは、最適化問題を方程式に帰着させる方法、結果であると言えよう。
ここで、そうして求められた方程式は微分記号が入っているとはいえ、ただの方程式であって微分方程式ではないことに注意が必要である。というのも、目的関数も実行可能領域も既知である為、それらの微分も所詮は関数であり、未知関数の微分が方程式内に出てこないからである。また、方程式であるからこそ解が数となり最適化問題の解と一致するのであって、もしこれが微分方程式であったならば、解が関数となり、その点でも不整合となってしまう。何かを最小(最大)化させる問題で求める解が関数である問題は変分法と呼ばれ、こちらは微分方程式と対応している。
数理計画法の使用
数理計画法の使用例として、金融工学における最適ポートフォリオの構築を扱う。価格が時間変化するリスク資産(例えば株のようなもの)が多数の種類与えられた時、どれをどの割合で持つのがリターン(資産価格の上昇率)を高めつつボラティリティ(資産価格変動の標準偏差)を低く出来るかという問題である。リスク資産の種類を$${n}$$種類、リスク資産$${i}$$の平均リターンを$${r_i}$$、リスク資産$${i}$$と$${j}$$の共分散を$${\sigma _{ij}}$$とする。第$${i}$$成分が$${r_i}$$となるリターンのベクトルを$${\bm{r}}$$、第$${(i,j)}$$成分が$${\sigma _{ij}}$$となる行列、即ち分散共分散行列を$${\Sigma}$$と置く。ポートフォリオはベクトル$${\bm{x}}$$で表すとし、その第$${i}$$成分$${x_i}$$がリスク資産$${i}$$の比率で、その合計は1と設定する($${\sum x_i = 1}$$)。
最小化すべき目的関数はポートフォリオ全体での分散(の$${\frac 1 2}$$倍である$${\frac 1 2 \bm{x}^{\rm T} \Sigma \bm{x}}$$)、制約条件が$${\sum x_i = \bm{1}^{\rm T} \bm{x} = 1}$$であることはすぐに分かる。但し、これだけではリターンの要素が入っていないので、何らかの形で付け加える必要がある。
リターンをモデルに付け加える方法は、主に3つある。
目的関数からリターンの定数倍を引く
制約条件にリターンが一定値であること(等号条件)を追加する
制約条件にリターンが一定値以上であること(不等号条件)を追加する
以上の前提を基に、資産保有として負の数量を許容するかしないか場合分けし、ポートフォリオ最適化問題を考察する。
負の資産保有を許容する場合
$${x_i < 0}$$を許容する場合を考える(株で言えば、これは空売りに相当する)。実はこの場合、ラグランジュの未定乗数法で最適化問題をただの連立1次方程式に変換出来る為、数理計画法を使わずとも行列演算のみによって最適解を導出することが出来る。リターンを制約条件に等号条件で付け加えた場合について、以下記述する。
最適化問題としてこれを記述すると、以下のようになる。
$$
\begin{rcases}
\text{Minimize} & \frac 1 2 \sum \sigma_{ij} x_i x_j \\
\text{Subject to} & \sum r_i x_i = r \\
& \sum x_i = 1
\end{rcases}
$$
ここで、$${r}$$はポートフォリオの目標リターンである。ラグランジュの未定乗数法で解く場合、目的関数の式に制約条件の未知数倍を足し合わせた式を最初に作る。
$$
\frac 1 2 \sum \sigma_{ij} x_i x_j + \lambda (r - \sum r_i x_i) + \mu (1 - \sum x_i)
$$
その後、$${x_i (i = 1,2, \dotsb ,n)}$$と$${\lambda}$$と$${\mu}$$、合計$${n+2}$$個の変数で微分し、それが0に等しいという$${n+2}$$個の方程式を作る。この方程式の解(の$${x_i (i = 1,2, \dotsb ,n)}$$の部分)が求める最適ポートフォリオである。
$$
\sum_j \sigma_{ij} x_j -\lambda r_i -\mu = 0 \ \ \ (i=1,2,\dotsb,n)\\
\sum r_i x_i = r \\
\sum x_i = 1
$$
これは連立1次方程式であるので、行列とベクトルを用いて書ける。
$$
\begin{bmatrix}
\sigma_{11} & \sigma_{12} & \dotsb & \sigma_{1n} & -r_1 & -1 \\
\sigma_{21} & \sigma_{22} & \dotsb & \sigma_{2n} & -r_2 & -1 \\
\vdots & \vdots & & \vdots & \vdots &\vdots \\
\sigma_{n1} & \sigma_{n2} & \dotsb & \sigma_{nn} & -r_n & -1 \\
r_1 & r_2 & \dotsb & r_n & 0 & 0 \\
1 & 1 & \dotsb & 1 & 0 & 0 \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
\vdots \\
x_n \\
\lambda \\
\mu \\
\end{bmatrix}
=
\begin{bmatrix}
0 \\
0 \\
\vdots \\
0 \\
r \\
1 \\
\end{bmatrix}
$$
線形関数の計算ライブラリを用いれば、左辺の行列と右辺のベクトルを与えることで求める解は得られる。
負の資産保有を許容しない場合
制約条件に$${x_i \geq 0}$$、即ち$${\bm{x} \geq \bm{0}}$$が加わる。この場合、ラグランジュの未定乗数法を用いても行列演算には帰着出来ない為、2次計画問題として数理計画法で解くことになる。幸いにも資産保有の非負制約がそのまま変数の非負制約であるから、不等号条件を追加しなければスラック変数による式変形をせずとも等式標準形で問題を記述することが出来る。そこで、今回はリターンについて、目的関数に組み込んで最適化問題を作るとすると、次のようになる。
$$
\begin{rcases}
\text{Minimize} & \frac 1 2 \bm{x}^{\rm T} \Sigma \bm{x} - c \bm{r}^{\rm T} \bm{x} \\
\text{Subject to} & \bm{1}^{\rm T} \bm{x} = 1 \\
&\bm{x} \geq \bm{0}
\end{rcases}
$$
ここで、$${c>0}$$はリターンの目的関数における重みであり、負の符号が付いているのはリターンが大きい方がより良い解であると設定する為である。等式標準形の式との対応関係は$${Q=\Sigma}$$、$${\bm{c} = -c \bm{r}}$$、$${A = \bm{1}^{\rm T}}$$、$${\bm{b} = 1}$$である。
2次計画問題の解き方であるが、今回は主双対内点法を使うものとする。主双対内点法では解$${(\bm{x},\bm{y},\bm{z})}$$を逐次的に更新する。$${k}$$番目の解を$${(\bm{x}_k,\bm{y}_k,\bm{z}_k)}$$とした時、そこから$${k+1}$$番目の解を得る手順は以下のようになる。但し、$${X_k = \text{diag}(\bm{x}_k)}$$、$${Z_k = \text{diag}(\bm{z}_k)}$$とし、$${\gamma, \lambda \in (0,1)}$$は定数である。
$$
\begin{gathered}
\mu_k = \frac {\bm{x}_k^{\rm T} \bm{z}_k}{n} \\
\nu_k = \gamma \mu_k \\
\Delta \bm{y}_k = (A(Z_k + X_k Q)^{-1} X_k A^{\rm T})^{-1} A (Z_k +X_k Q)^{-1} (X_k \bm{z}_k - \nu_k \bm{1}) \\
\Delta \bm{x}_k = (Z_k + X_k Q)^{-1}(X_k A^{\rm T} \Delta \bm{y}_k - X_k \bm{z}_k + \nu_k \bm{1}) \\
\Delta \bm{z}_k = Q \Delta \bm{x}_k - A^{\rm T} \Delta \bm{y}_k \\
\hat{\alpha} = \max \lbrace \alpha | \bm{x}_k + \alpha \Delta \bm{x}_k \geq \bm{0}, \bm{z}_k + \alpha \Delta \bm{z}_k \geq \bm{0} \rbrace \\
\alpha = \lambda \hat{\alpha} \\
\bm{x}_{k+1} = \bm{x}_k + \alpha \Delta \bm{x}_k \\
\bm{y}_{k+1} = \bm{y}_k + \alpha \Delta \bm{y}_k \\
\bm{z}_{k+1} = \bm{z}_k + \alpha \Delta \bm{z}_k \\
\end{gathered}
$$
この繰り返しの終了は、十分に小さい数$${\varepsilon > 0}$$を用意して以下の3つの不等式が全て成り立つか、反復回数$${k}$$が一定の上限値を超えた時とする。
$$
\begin{gathered}
||A \bm{x} - \bm{b}|| < \varepsilon \\
||A^{\rm T} \bm{y} + \bm{z} - Q \bm{x} - \bm{c}|| < \varepsilon \\
\mu_k < \varepsilon \\
\end{gathered}
$$
初期の実行可能解を決めてやれば、以上の手順により非負制約の最適ポートフォリオを導出することが出来る。
この記事が気に入ったらサポートをしてみませんか？