金融のための量子アニーリングことはじめ3 | ポートフォリオ最適化 2

#はじめに

某社でクオンツアナリストとして働いておりますhalと申します。
本記事は自分の勉強も兼ねて、主に金融実務者の方向けに典型的なポートフォリオ最適化を通じて、
"量子アニーリング"の紹介&デモンストレーションをしていきます。
pythonの導入されたPCがあればどなたでも無料で試せるので、興味がある方は実際に動かしてみてください。

本記事ではD-Wave社のLeapという無料の量子アニーリングサービスを使用してきます。
ことはじめ1で導入方法を紹介しているのでこちらもご覧ください。

#おさらい
前回では最適化問題をleapで解かせるには問題をQUBO形式で与える必要があることを紹介しました。
またQUBO形式では0,1の解しか得られないため、任意のウェイトでの最適化は工夫が必要と書きましたが、今回はこの工夫について考えていきます。

#最適ポートフォリオ

N種類の資産を投資対象銘柄として、最適ポートフォリオを求める。
目的関数は以下のように定義する。

\begin{align}
\mathrm{minimize}~ ~ &\mathcal{H}=-\sum_i r_i w_i +\gamma \sum_{i,j} \sigma_{i,j} w_i w_j=-\boldsymbol{ r}^\mathrm{T}\cdot \boldsymbol{w}+\gamma \boldsymbol{w}^\mathrm{T}\boldsymbol{C}\boldsymbol{w}\\
\mathrm{subject ~ to}~~&||\boldsymbol{w}||=1
\end{align}

#解の定義

前回とは違い、１資産に複数の２値変数（量子ビット)を割り当てることで資産のウェイトの表現をしていきます。
今回は対応する1の個数を数えるような定義にします。
即ち、資産iにD個の変数が割り当てられてるとして、

w_i=\sum_{d=1}^D  q_d

とウェイトを表現し、制約条件として以下の式を与えます

\sum_{i=1}^N w_i=D

最後にDでウェイトを規格化すれば1/Dの精度で最適ポートフォリオが求まりますね。

QUBO形式に変換していくためにまずは制約条件を目的関数に取り込みましょう。
ラグランジュ乗数$\lambda$を用いて目的関数を以下のように書き換えます

\begin{align}
\mathrm{minimize}~ ~ &\mathcal{H}=-\boldsymbol{ r}^\mathrm{T}\cdot \boldsymbol{w}+\gamma \boldsymbol{w}^\mathrm{T}\boldsymbol{C}\boldsymbol{w}+\lambda\left(D-\sum_{i=1}^N w_i\right)^2
\end{align}

#QUBO形式に変換

まずはウェイトベクトルと変数（量子ビット)ベクトルの変換を考えて行きましょう
わかりやすさのために、資産数N=3、１資産あたりのビット数D=2の場合で考えていきます。

まずは変数ベクトル$q$を以下のように対応させます

    \boldsymbol{q}=(\underbrace{q_1,q_2}_{銘柄1},\underbrace{q_3,q_4}_{銘柄2},\underbrace{q_5,q_6}_{銘柄3})^T

このときウェイトベクトルは以下のように変換させることができます。

\begin{align}
   \left( \begin{array}{c}{w_1\\w_2\\w_3}\end{array} \right)&=
   \underbrace{\left(\begin{array}{cccccc}
       1 & 1 & 0 & 0 & 0 & 0 \\
       0 & 0 & 1 & 1 & 0 & 0 \\
       0 & 0 & 0 & 0 & 1 & 1  
    \end{array}\right) }_\boldsymbol{B}
    \left(\begin{array}{c}{q_1\\q_2\\q_3\\q_4\\q_5\\q_6}\end{array}\right)\\

    \boldsymbol{w}&=\boldsymbol{Bq}
\end{align}

この変換行列$B$を用いてQUBO形式に変換していきます

###第１項

無理矢理２次式に変換するため第１項は少しテクニカルです
まずそのまま行列$B$を用いて変換すると

\boldsymbol{r^Tw}=\boldsymbol{r^TBq}

２次形式にするためには左から$\boldsymbol{q^T}$を掛ける必要がありますが$q_i^2=q_i$であることを利用して６次元単位行列$\boldsymbol{E}$を用いて

\boldsymbol{r^TBq}=\boldsymbol{q^T \mathrm{diag}(r^TB )q}

###第２項

$\boldsymbol{w^T}=\boldsymbol{q^TB^T}$であることを利用すれば簡単に変換できます。

\gamma\boldsymbol{w^TCw}=\gamma\boldsymbol{q^T(B^TCB)q}

###第３項

少し複雑ですが第1項と第2項でやった事を応用すれば変換できます。

\begin{align}
\left(D-\sum_i^N w_i\right)^2 &=w_1^2+w_2^2+w_3^2
                               +2w_1w_2+ 2w_2w_3+2w_3w_1
                               -2Dw_1-2Dw_2-2Dw_3
                               +D^2\\


                             &=\boldsymbol{w}^T
                             \underbrace{ \left( \begin{array}{ccc}
                             1 & 1 & 1\\
                             1 & 1 & 1 \\
                             1 & 1 & 1
                             \end{array} \right) }_\boldsymbol{A}
                             \boldsymbol{w}-2D\underbrace{(1,1,1)}_{\mathbb{1}^T}\boldsymbol{w}\\

                             &=\boldsymbol{q}^T(\boldsymbol{B}^T \boldsymbol{AB)q}-2D\boldsymbol{q^T}\mathrm{diag}(\mathbb{1}^T\boldsymbol{B})\boldsymbol{q}\\

                            &=\boldsymbol{q}^T(\boldsymbol{B^TAB}-2D~\mathrm{diag}(\boldsymbol{\mathbb{1}^TB}) )\boldsymbol{q}\\ \\



\end{align}

###QUBO化
さあこれで全部QUBO形式に変換できました。以上の式をまとめると最適化問題は

\mathcal{H}=\boldsymbol{q}^\mathrm{T}  \boldsymbol{Q} \boldsymbol{q}\\
\boldsymbol{Q}=-\mathrm{diag}(\boldsymbol{r^TB})+\gamma\boldsymbol{B^TCB}+\lambda (\boldsymbol{B^TAB}-2D~\mathrm{diag}(\boldsymbol{\mathbb{1}^TB}) )

となります。Bが増えてますが、基本前回と同じですね。あとはこれをLeapで計算してもらうだけです。

#実際に解いてみる
ではことはじめ2と同じ設定の問題を解いてみましょう.

$N=5,D=10$すなわち5資産での最適ポートフォリオで、ウェイトは10%刻みであるとします。
分析データは以下のように定義します。（適当です)
次回の記事と整合性を取るために共分散行列を通常の対角行列として表示します。問題としては同値です。

\begin{align}
    \boldsymbol{r}&=(1\%,2\%,3\%,4\%,5\%)\\　\\
    \boldsymbol{C}&=\left( \begin{array}{ccccc}
                           0.01   & -0.0062  & 0.0238  & 0.0209 &  0.0096 \\
                          -0.0062 &  0.0225  & 0.01455 & 0.0264 &  0.0147 \\
                           0.0238 &  0.01455 & 0.04    & 0.0114 &  0.0049 \\
                           0.0209 &  0.0264  & 0.0114  & 0.0625 &  0.0057 \\
                           0.0096 &  0.0147  & 0.0049  & 0.0057 &  0.09   
                             \end{array} \right)
\end{align}

また$\gamma=1.2 ,\lambda =0.2$とします。
この問題をLeapで最適化を行うコードは以下のようになります。

import numpy as np
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite


#リターンベクトル(縦ベクトルとして定義)
r=np.array([[0.01,0.02,0.03,0.04,0.05]]).T
#共分散行列
C=np.array([[ 0.01   , -0.0062 ,  0.0238 ,  0.0209 ,  0.0096 ],\
            [-0.0062 ,  0.0225 ,  0.01455,  0.0264 ,  0.0147 ],\
            [ 0.0238 ,  0.01455,  0.04   ,  0.0114 ,  0.0049 ],\
            [ 0.0209 ,  0.0264 ,  0.0114 ,  0.0625 ,  0.0057 ],\
            [ 0.0096 ,  0.0147 ,  0.0049 ,  0.0057 ,  0.09   ]])

#N=5資産D=10量子ビットを各資産に対応
N=5
D=10

#ガンマとラムダを定義
gamma=1.2
lam=0.2

#QUBO matrix を計算

#変換の為の配列を作成

#変換行列B(表記はわかりにくですね..)
B=np.vstack([np.roll([(1 if i<D else 0) for i in range(N*D)],j*D) for j in range(N)])

#要素が全て1のベクトル
l=np.ones(N)

A=np.ones((N,N))

Q=-np.diag(np.dot(r.T,B))+gamma*np.dot(B.T,np.dot(C,B))+lam*(np.dot(B.T,np.dot(A,B))-2*K*np.diag(np.dot(l,B)))

Q=np.triu(Q)+np.triu(Q,k=1) #非対角成分を二倍し上三角行列に
Q=Q/np.max(np.abs(Q))  #規格化しないと桁落ちにより正しい解が得られないことがあります

#問題を与える時は辞書にして渡す必要があるので変換

Qdict=dict()

for i in range(N*D):
    for j in range(N*D):
        if i<=j:
            Qdict.update({(i,j):Q[i,j]})

#問題をapiに投げ、結果を取得 (*)
sampler = EmbeddingComposite(DWaveSampler()) 
response = sampler.sample_qubo(Qdict,num_reads=10000)

#目的関数が最小な解を取得し、ウェイトに変換
answer=response.lowest().record[0][0]
weight=np.dot(B,answer)/D
print(weight)

また、問題をapiに投げる部分を以下のように変更すると、自分のPC上で量子コンピューターのシミュレーションを行い計算ができます。
シミュレーションだと少々時間はかかるものの、マシンタイムは消費されません。

#問題をapiに投げ、結果を取得 (*)
import neal
sampler = neal.SimulatedAnnealingSampler() 
response = sampler.sample_qubo(Qdict,num_reads=10000)

さて結果をみてみると以下のようになります(leapでの結果は異なっていることがあります)

	結果	（元々の)目的関数の値
leapでの結果	[0.6 ,0.3 ,0 ,0.1 ,0 ]	-0.006268
シミュレーションでの結果	[0.6 ,0.4 ,0 ,0 ,0 ]	-0.008931

結果が違っていますね、これはどういう事でしょう

##結果が違う！

なぜleapでの結果とシミュレーションでの結果が違うのでしょうか?
目的関数の値をみてみるとシミュレーションでの結果の方が本当の最適解であるとわかります。
なぜleapでは最適解が得られなかったのでしょう?

Qの精度

実はQの成分は-1から+1までで、精度には小数点第2位までの制限があります。
制約条件$\lambda$を大きくするとウェイト制約を満たした解が得られやすくなりますが、本来の目的関数の値が小さくなってしまい、値が切り捨てられてしまい本当の最適解に到達できなくなります。
逆に$\lambda$を小さくすると、目的関数の精度は上がりますが制約を満たした解が得られなくなってしまいます。

即ち精度と制約条件はトレードオフの関係であり、両方を満足するような$\lambda$を定めてあげる必要がありますが、今回の問題ではどうやっても両立する$\lambda$を見つけることはできませんでした。
そのため、leapでは最適な結果を得ることに失敗してしまいました。対してシミュレーションであれば精度の限界はないために、最適な結果が得られたということです。

D-Wave leapは量子アニーリングマシーンであり、組み合わせ最適化に特化したコンピューターです。
対してポートフォリオ最適化は連続最適化であり、今回は半ば強引に組み合わせ最適化に変換して解くことを試みました。
今回の試みでより高精度な量子アニーリングマシーンが出てくれば解いていくことが可能でしょうが、現状では難しいみたいですね。

#おわりに
いかがでしたでしょうか？当初ではちゃんと最適な結果を得られる予定だったのですが、不本意ながら精度の問題でうまくいきませんでしたね。
より拡張したパターンも用意していたのですが...
何か改善方法をご存知でしたら、コメント等で教えていただければ幸いです。

今回の記事では出来るだけ丁寧に書いてみたつもりですので、ご自分の環境でもぜひ試してみてください。

#リンク
金融のための量子コンピューターことはじめ1 | D-wave Leap導入編
金融のための量子コンピューターことはじめ2 | ポートフォリオ最適化 1
金融のための量子コンピューターことはじめ3 | ポートフォリオ最適化 2(本記事)