リバース量子アニーリングを用いたポートフォリオ最適化

リバース量子アニーリングを用いたポートフォリオ最適化#

はじめに#

この章ではポートフォリオ最適化問題を例として、OpenJijを用いてリバース量子アニーリング（Reverse Quantum Annealing）という、古典的な最適化計算と、量子アニーリングを組み合わせた手法を紹介します。このチュートリアルでは、古典的なアルゴリズム、普通の量子アニーリング、そしてリバース量子アニーリングそれぞれで同じポートフォリオ最適化問題を解き、その結果を比較します。またその過程を通して、OpenJijでの量子アニーリングの実装や実験を行う時の注意点についても触れていきます。

注意：このポートフォリオの内容を正しく実行するにはOpenJij0.4.9以上が必要です。必要な場合は実行前に次のコマンドを実行し、OpenJijのアップデートを行なってください。

pip install -U openjij

ポートフォリオ最適化問題#

資産運用のための投資活動を行うとき、誰しも「できるだけリスクを回避しながら大きな収益を実現させたい」と考えるはずです。よって収益が限定的ですが、リスクが小さい（ないし全くない）資産と、見込める収益が大きい分リスクも大きい資産を組み合わせてポートフォリオを作り、分散投資を行うのが一般的な戦略となります。

この時、与えられた一定のリスクでは最大利益を実現したい、または同じ利益を実現する場合にできるだけ小さいリスクを取るには、最適なポートフォリオを選ぶ必要があります。現在、よく利用される手法にMarkowitzによる現代ポートフォリオ理論(Modern Portfolio Theory)があります。ポートフォリオを構成する銘柄の間の相関を考慮し、その共分散を最小にすることを目指す手法です[1]。

\min \sum_{i=1}^{N}\sum_{j=1}^{N}w_{i}w_{j}\sigma_{ij}; \quad \sum_{i=1}^{N}W_{i}=1; \quad \sum_{i=1}^{N}W_{i}\mu_{i}=R

$w_{i}$ は各銘柄がポートフォリオ内に占める重み、 $\sigma_{ij}$ は銘柄間の共分散を表します。 $\mu_{i}$ は各銘柄の期待収益で、 $R$ はこのポートフォリオの総収益です。

この時、ポートフォリオのパフォーマンスの評価によく使われる指標としてSharpe Ratioというものがあります。

Sharpe Ratio#

Sharpe Ratioは投資の効率性を評価する指標で、銘柄単体とポートフォリオに対して計算されます。1966年にウィリアム・シャープによって提案されました[3]。その基本的な定義を以下に示します。

S = \frac{\bar{r}-r_0}{\sigma}

ここで $\bar{r}$ は一定時間内の平均的なリターン率で、月間や年間リターンの平均を利用することが多くあります。 $r_0$ はリスクなしの場合のリターン率を表し、普段は安定した国債の利子率などを用います。そして $\sigma$ は資産価値の変動の激しさを意味するボラティリティです。

実際にSharpe Ratioを計算するときには、場面によって異なる計算方法でこれらの値を算出します。例えばリターンは単純リターンとして $r=P(t+1)/P(t) - 1$ で資産価格 $P$ を使って計算される場合と、 $r=log[P(t+1)/P(t)]$ のように対数リターン率として計算される場合が存在します。そしてボラティリティも価格の変動から計算する手法とリターン率からの計算するやり方があります。具体的に興味がある方は金融工学や市場分析などに関連する本とサイトに参照すると良いでしょう。

ポートフォリオの評価と最適化#

Sharpe Ratioの定義から、その値が大きいということは収益が大きいまたはリスク（ボラティリティ）が小さいことを意味します。ある与えられた銘柄の集合から作られたポートフォリオのうち、最大のSharpe Ratioを実現できるような組合せは、単位リスクに対して最大のリターンを得る組み合わせと言えます。よってポートフォリオの最適化問題が一定の考え方の元では、Shape Ratioを最大化する問題に置き換えることができます。

現実では、この最適化は銘柄を選別した上で、それぞれの銘柄に投入する資金の量を決定しなければなりません。すなわち、これは選ぶ銘柄と各銘柄のウェイトを決定する組合せ最適化問題となります。このチュートリアルは、参考にした論文[2]に従い、簡単のために各銘柄のウェイトを等しくする戦略について考えます。 $M$ 個ある銘柄のうち $N$ 個銘柄を選び出し、Sharpe Ratioを最大にする組合せを求めます。すると、次のようなQUBO形式に問題を翻訳することができます。

\mathcal{O}(\boldsymbol{q})=\sum_{i=1}^{N}a_{i}q_{i}+\sum_{i=1}^{N}\sum_{j=i+1}^{N}b_{ij}q_{i}q_{j}

ここで $q_{i}$ は各銘柄に対する選択を表すバイナリ変数で、 $1$ であればポートフォリオに組み入れ、 $0$ であれば組み入れないことを意味します。 $a_i$ は銘柄自体のパフォーマンスによる魅力スコアであり、 $b_{ij}$ は銘柄間のペアワイズ相関で決められた罰金または賞金の度合を表します。

具体的に $a_i$ と $b_{ij}$ は次の表1に従って決定されます。

../../../_images/reverse_annealing_RQA_QUBO.png

この表の $a_i$ のグループは、ポートフォリオを構成する人によって決められた基準であり、銘柄を魅力的であるものから順に並べられたものです。それらを順位付けした後に、均等に分けた場合に形成されたグループを示しています。この魅力的な基準はリターン率の高さなどの要素も取り入れられますが、このチュートリアルでは単純に銘柄単体のSharpe Ratio順を用いています。次の $b_{ij}$ は対数リターンの時系列から求められた相関行列の成分 $\rho_{ij}$ の値によって決められたものです。

このQUBOを用いて量子アニーリングを実行すれば、最適なポートフォリオを得ることができるはずです。また、選択銘柄の数について指定したい場合は、銘柄数に対して制約条件を付け加えます。ただし、参考にした論文[2]によると、制約条件を罰金法でQUBOに付け加えた場合、導入した大きい相互作用（ペナルティ）はアニーリング性能や最終的に取得した結果にも影響することが示されています。現実の問題では、最適なポートフォリオに含まれる銘柄の数は分からないため、このチュートリアルでも同様に銘柄数制限なしての最適化を行います。

Reverse Quantum Annealingによるポートフォリオ最適化#

現実のポートフォリオ最適化の問題は基本、大きな銘柄数を持つ市場からいくつかの銘柄を選択することになります。候補銘柄数が増加すると、最適なポートフォリオを作るための計算量も大きくなります。よって一定の条件の下に選別を行い、その次に最適化アルゴリズを適応するのが一般的です。そして、銘柄の組合せの自由度から、最終的に実現したSharpe Ratioが最大Sharpe Ratioと極めて近い値を持つような局所解が多く存在することもあり得ます。伝統的なアルゴリズムがそのようなエネルギーランドスケープの谷に囚われた場合、そこから抜け出すのにも時間とリソースが必要となります。そのため、以上の既存の問題点を回避しつつポートフォリオ最適化問題を解くには、量子アニーリングのような大きい数の相関を取り扱うことができ、かつ局所解から脱却できる計算手法の方が上手に解くことができます。

ただし参考論文[2]に示されたように、単純なForward Annealingのみを行う場合では、銘柄数の増加により、最適解までにたどり着く時間も大きく増加します。結果としてこれまでの遺伝的アルゴリズムと近いパフォーマンスになることが判明しています。その理由としては先述の通り、最適解と近い値のSharpe Ratioを持つ解と最適解はエネルギーの違いが小さく、最適解へのエネルギー変化がほぼない準位間に高いポテンシャル障壁があるような遷移となるため、遷移が起こるまでにより長い時間が必要となるためでし。また量子アニーリングマシン実機では、通常のアニーリングを行う途中で熱ノイズなどにより系が励起されて時間発展する場合があります。これらの要因により、最終結果が最適解とならないこともあります。

Reverse Quantum Annealing#

これら量子アニーリングを用いた時の問題点を解決するために提案された手法の一つに、Reverse Quantum Annealingがあります。名前通り、通常のForward Annealingと逆方向のステップを導入するような量子アニーリングを実行することを意味します。通常のForward Annealingの時間依存のハミルトニアンは次のように書かれます。

\mathcal{H}_{\mathrm{QA}}(t)=A[t]\sum_{i=1}^N\sigma_{i}^{\mathrm{x}}+B[t]\mathcal{H}_{\chi-\mathrm{Ising}}

ここで $A[t]$ はアニーリングを行うとき系にかける横磁場の大きさ、 $B[t]$ は問題のハミルトニアン $\mathcal{H}_{\chi-\mathrm{Ising}}$ の振幅を表します。次の図でのa)で示されるように、Forward Annealingの場合はアニーリングの進行に従って、 $A[t]$ を徐々に小さくし、同時に $B[t]$ を支配的にさせていきます。アニーリング完了時では $A[t]=0$ となります。

../../../_images/reverse_annealing_QAandRQA.png

それに対し、Reverse Annealingはb)のように $A[t]=0,B[t]=1$ の状態から出発します。このとき、まず通常とは逆に $B[t]$ を小さくしていきます。これをReverse Phaseと呼びます。この段階では系のエネルギーが逆に高くなるため、c)に示されているようなAB間のポテンシャル障壁を越えやすくなります。そして $A[t], B[t]$ をある値に到達させた後、系をそのままの状態を維持してしばらく待ちます。これをPause Phaseと呼びます。この間も高いエネルギーの状態が維持されるため、熱的ホッピングが起こりやすい状況となります。するとBの右側のようなポテンシャル障壁も越えやすくなっています。これにより最初にいた局所解から脱し、最適解付近にたどり着く可能性が上がることが示唆されます。最後にまた $A[t]$ を小さくし、 $B[t]$ を大きくさせるForward Phaseで通常のアニーリングを実行することで、最適解を取得します。このプロトコルでは、Reverse phadeで $A[t], B[t]$ を停止させた値や、各段階での継続時間などにより、最終的に得られる解のパフォーマンスにも影響することが考えられます。

Reverse Quantum Annealingによるポートフォリオ最適化の手順#

ここまでの議論から、Reverse annealingではReverse Phaseの初期状態が必要です。ランダムに生成された初期状態を用いることも可能ですが、手法が提案された背景からもわかるように、ある種の局所解の状態を用意したほうが最適解探索の効率が向上することが予想されます。参考論文[2]では初期状態となる局所解を得るために、伝統的なアルゴリズムで得られた出力を利用しています。このチュートリアルもそれに従って、Revearse annealingを実行してみましょう。

以上から、これから行う実装は基本論文[2]の流れを参考にして以下のようになります。

最適化を行う銘柄データの生成
古典なアルゴリズムの実装とその結果の確認
通常のForward Annealingによる最適解探索
Reverse Quantum Annealingによる最適解探索
Reverse Quantum Annealingのパラメータ探索

銘柄データの生成と古典アルゴリズムの実装#

最適化を行う銘柄データの生成#

参考論文[2]にある方法に従い、与えられた初期値を用いて、ブラウン運動による銘柄チャートを生成します。

import numpy as np
import matplotlib.pyplot as plt
import random

#Magic numbers to generate assets
rho = 0.1 # input uniform correlation
mu = 0.075 # expected value
sigma = 0.15 # volatility/standard error
r0 = 0.015 # no risk return

参考論文[2]の付録Aから、各時刻においてチャートの運動は前時刻の運動を用いて

S(t_{n+1})=S(t_n)\exp(\mu-\frac{1}{2}\sigma^2)\Delta t + \sigma Z_n\sqrt{\Delta t}

のように与えられます。ここで $Z_n$ はCholesky分解で作られた一様相関行列 $\rho$ を従う多変量正規分布です。

これを実行し、チャートをプロットして様子を確認しましょう。初期に同じ価値から始め12ヶ月の推移を見ると、全体が広がっていく様子がわかります。一部銘柄では大きく上昇または降下する振る舞いも確認できます。

#相関正規確率変数Znの生成
def createZvariables(N, rho):
  rho_mat = np.full((N,N), rho)
  rho_mat[range(N), range(N)] = 1.0
  rho_chole = np.linalg.cholesky(rho_mat)
  zNs_temp = np.random.normal(0, 1, (10000, N))
  zNs = zNs_temp @ rho_chole
  return zNs

#チャートのブラウン運動の生成
def GetNextSt(St, mu, sigma, zN):
  Deltat = 1
  scale = np.exp((mu-0.5*sigma*sigma)*Deltat + sigma*zN*np.sqrt(Deltat))
  NextSt = St * scale
  return NextSt

Nassets = 48 #生成する銘柄の数
chart = list()
ZList = list()
Zvariables = createZvariables(Nassets, rho)
Zlabels = random.sample([x for x in range(10000)], 12) #Z variable shuffle

for label in Zlabels:
  ZList.append(Zvariables[label])

for iasset in range(Nassets):
  chart_asset = [1.0] #初値 相対価格で1.0

  #12ステップ（か月）シミュレーションを行う
  for month in range(12):
    chart_asset.append(GetNextSt(chart_asset[month], mu, sigma, ZList[month][iasset]))

  chart.append(chart_asset)
  #print(chart_asset) #各銘柄の具体的な数値をプリント
  plt.plot(list(range(13)), chart_asset)

plt.xlabel("Month", loc="right")
plt.ylabel("Relative Price", loc="top")
plt.show()

../../../_images/9bb7c157c2a447411462762ab355bacd6998ee317a2d364cc0c28b90488eda0e.png

次に各銘柄のSharpe Ratioを計算します。ここでは先程の結果を評価するに実現Sharpe Ratioを用います。その定義に従い、無リスク利回り率を超過した超過収益率の平均 $\bar{r}$ と超過リターン率の標準偏差 $\sigma$ から求めましょう。チャートの確認から判明したように、銘柄数が少ない場合は偶然による偏りが大きい・生成条件に銘柄間の相関による影響があるため、銘柄を100回生成しその平均Shape Ratioを確認することでそれらの影響を減らします。結果として、銘柄の平均Sharpe Ratioが期待通り $0.4$ になるのが確認できます。またSharpe Ratioの計算定義を変更した場合には、当然Sharpeの値が変化することが確認できますが、おおむね $0.40\pm0.05$ の範囲内には収まることもわかります。

#銘柄生成関数
def CreateAssets(Nassets): 
  chart = list()
  ZList = list()
  Zvariables = createZvariables(Nassets, rho)
  Zlabels = random.sample([x for x in range(10000)], 12)
  for label in Zlabels:
    ZList.append(Zvariables[label])
  for iasset in range(Nassets):
    chart_asset = [1.0]
    for month in range(12):
      chart_asset.append(GetNextSt(chart_asset[month], mu, sigma, ZList[month][iasset]))
    chart.append(chart_asset)
  return chart

#オーソドックスなSharpe Ratio
def CalculateSharpeRatio(asset):
    monthly_return = list()
    for month in range(12):
      valueChange = asset[month+1]/asset[month] - 1.0 - r0
      monthly_return.append(valueChange)

    annualized_excess_return = np.mean(monthly_return)
    volatility = np.std(monthly_return, ddof=1)

    return(annualized_excess_return/volatility)

#log-returnに基づいたSharpe Ratio
def CalculateSharpeRatio_log(asset):
  monthly_log_return = list()
  for month in range(12):
    valueChange = np.log(asset[month+1]/asset[month])
    monthly_log_return.append(valueChange)
  mean_log_return = np.mean(monthly_log_return)
  volatility = np.std(monthly_log_return, ddof=1)

  return (mean_log_return-r0)/volatility

#少し変わった定義のSharpe Ratio
def CalculateSharpeRatio_2(asset):
  monthly_log_return = list()
  monthly_return = list()
  for month in range(11):
    log_valueChange = np.log(asset[month+1]/asset[month])
    monthly_log_return.append(log_valueChange)
    valueChange = asset[month+1]/asset[month] - 1.0
    monthly_return.append(valueChange)
  annualized_excess_return = np.mean(monthly_return) - r0
  volatility = np.std(monthly_log_return, ddof=1)

  return (annualized_excess_return/volatility)

allmean = 0
for ntry in range(100): #複数回試行で同じセットないの相関による影響を消す
  Chart = CreateAssets(48) #指定した数の銘柄を作る
  mean_SR = 0.0
  n = 0.
  for asset in Chart:
    assetSR = CalculateSharpeRatio(asset) #ここを変更すればShape Ratioの計算方法を変えられる
    #assetSR = CalculateSharpeRatio_log(asset)
    #assetSR = CalculateSharpeRatio_2(asset)
    mean_SR = ((mean_SR*n)+assetSR) / (n+1) #平均Shape Ratio
    n+=1
  print("SubSet "+ str(ntry)+ " average Sharpe Ratio: " + str(mean_SR))
  allmean = ((allmean*ntry)+mean_SR) / (ntry+1)
  print("Average Sharpe Ratio of all generated: "+ str(allmean))

SubSet 0 average Sharpe Ratio: 0.3405803595058439
Average Sharpe Ratio of all generated: 0.3405803595058439
SubSet 1 average Sharpe Ratio: 0.3892361471066652
Average Sharpe Ratio of all generated: 0.36490825330625454
SubSet 2 average Sharpe Ratio: 0.39545623751586306
Average Sharpe Ratio of all generated: 0.3750909147094574
SubSet 3 average Sharpe Ratio: 0.4798707033364349
Average Sharpe Ratio of all generated: 0.40128586186620174
SubSet 4 average Sharpe Ratio: 0.5200192047618962
Average Sharpe Ratio of all generated: 0.4250325304453407
SubSet 5 average Sharpe Ratio: 0.29779873253384015
Average Sharpe Ratio of all generated: 0.40382689746009054
SubSet 6 average Sharpe Ratio: 0.6332118455733746
Average Sharpe Ratio of all generated: 0.4365961757619883
SubSet 7 average Sharpe Ratio: 0.3567160464145213
Average Sharpe Ratio of all generated: 0.4266111595935549
SubSet 8 average Sharpe Ratio: 0.43047549262701645
Average Sharpe Ratio of all generated: 0.42704052993060615
SubSet 9 average Sharpe Ratio: 0.4449095866367027
Average Sharpe Ratio of all generated: 0.4288274356012158
SubSet 10 average Sharpe Ratio: 0.49487010133589027
Average Sharpe Ratio of all generated: 0.43483131430436806
SubSet 11 average Sharpe Ratio: 0.4655791900224186
Average Sharpe Ratio of all generated: 0.4373936372808722
SubSet 12 average Sharpe Ratio: 0.49768825617610396
Average Sharpe Ratio of all generated: 0.4420316848881977
SubSet 13 average Sharpe Ratio: 0.4111569709928875
Average Sharpe Ratio of all generated: 0.4398263481813898
SubSet 14 average Sharpe Ratio: 0.31836607129949707
Average Sharpe Ratio of all generated: 0.4317289963892636
SubSet 15 average Sharpe Ratio: 0.24272766235135687
Average Sharpe Ratio of all generated: 0.4199164130118945
SubSet 16 average Sharpe Ratio: 0.48880125742023917
Average Sharpe Ratio of all generated: 0.4239684626829736
SubSet 17 average Sharpe Ratio: 0.36776075046979634
Average Sharpe Ratio of all generated: 0.42084581200446375
SubSet 18 average Sharpe Ratio: 0.39158143864654593
Average Sharpe Ratio of all generated: 0.4193055818277312
SubSet 19 average Sharpe Ratio: 0.36792052818624854
Average Sharpe Ratio of all generated: 0.41673632914565706
SubSet 20 average Sharpe Ratio: 0.5142335911180681
Average Sharpe Ratio of all generated: 0.421379055906248
SubSet 21 average Sharpe Ratio: 0.4168616287280342
Average Sharpe Ratio of all generated: 0.4211737183072383
SubSet 22 average Sharpe Ratio: 0.7586123339463494
Average Sharpe Ratio of all generated: 0.4358449624654605
SubSet 23 average Sharpe Ratio: 0.31548552953471654
Average Sharpe Ratio of all generated: 0.4308299860933462
SubSet 24 average Sharpe Ratio: 0.34710371637581416
Average Sharpe Ratio of all generated: 0.4274809353046449
SubSet 25 average Sharpe Ratio: 0.38280042932643976
Average Sharpe Ratio of all generated: 0.42576245430548315
SubSet 26 average Sharpe Ratio: 0.599196175106843
Average Sharpe Ratio of all generated: 0.4321859254462743
SubSet 27 average Sharpe Ratio: 0.4134693918661277
Average Sharpe Ratio of all generated: 0.43151747781841193
SubSet 28 average Sharpe Ratio: 0.4779120342518383
Average Sharpe Ratio of all generated: 0.43311729010921973
SubSet 29 average Sharpe Ratio: 0.6207867616987124
Average Sharpe Ratio of all generated: 0.4393729391622028
SubSet 30 average Sharpe Ratio: 0.4092133768439145
Average Sharpe Ratio of all generated: 0.43840005005516125
SubSet 31 average Sharpe Ratio: 0.5463353502043146
Average Sharpe Ratio of all generated: 0.4417730281848223
SubSet 32 average Sharpe Ratio: 0.314524346679724
Average Sharpe Ratio of all generated: 0.4379170075331526
SubSet 33 average Sharpe Ratio: 0.4348002318379652
Average Sharpe Ratio of all generated: 0.4378253376597647
SubSet 34 average Sharpe Ratio: 0.4185786684737743
Average Sharpe Ratio of all generated: 0.43727543282587933
SubSet 35 average Sharpe Ratio: 0.3985495346596739
Average Sharpe Ratio of all generated: 0.4361997134323736
SubSet 36 average Sharpe Ratio: 0.44006339546364037
Average Sharpe Ratio of all generated: 0.4363041372710565
SubSet 37 average Sharpe Ratio: 0.34392044839655816
Average Sharpe Ratio of all generated: 0.4338729875638329
SubSet 38 average Sharpe Ratio: 0.5572400389560374
Average Sharpe Ratio of all generated: 0.43703624529183815
SubSet 39 average Sharpe Ratio: 0.4165844170463176
Average Sharpe Ratio of all generated: 0.4365249495857001
SubSet 40 average Sharpe Ratio: 0.4610683975869365
Average Sharpe Ratio of all generated: 0.43712357026865706
SubSet 41 average Sharpe Ratio: 0.4739321768875668
Average Sharpe Ratio of all generated: 0.43799996566434535
SubSet 42 average Sharpe Ratio: 0.3337306742853509
Average Sharpe Ratio of all generated: 0.43557509842297343
SubSet 43 average Sharpe Ratio: 0.39090459268452404
Average Sharpe Ratio of all generated: 0.4345598596561905
SubSet 44 average Sharpe Ratio: 0.4900628824390543
Average Sharpe Ratio of all generated: 0.4357932601624764
SubSet 45 average Sharpe Ratio: 0.4250312557461802
Average Sharpe Ratio of all generated: 0.4355593035447308
SubSet 46 average Sharpe Ratio: 0.25225098889645875
Average Sharpe Ratio of all generated: 0.4316591266373208
SubSet 47 average Sharpe Ratio: 0.37656606084815075
Average Sharpe Ratio of all generated: 0.4305113544333798
SubSet 48 average Sharpe Ratio: 0.5245550873009045
Average Sharpe Ratio of all generated: 0.43243061428781904
SubSet 49 average Sharpe Ratio: 0.2933288773009551
Average Sharpe Ratio of all generated: 0.4296485795480817
SubSet 50 average Sharpe Ratio: 0.4506599747848008
Average Sharpe Ratio of all generated: 0.43006056768997825
SubSet 51 average Sharpe Ratio: 0.3248013768074119
Average Sharpe Ratio of all generated: 0.4280363524806981
SubSet 52 average Sharpe Ratio: 0.47764556581182505
Average Sharpe Ratio of all generated: 0.4289723753737382
SubSet 53 average Sharpe Ratio: 0.43525180045212314
Average Sharpe Ratio of all generated: 0.4290886610233379
SubSet 54 average Sharpe Ratio: 0.4017105105087737
Average Sharpe Ratio of all generated: 0.42859087646852767
SubSet 55 average Sharpe Ratio: 0.36164210045453293
Average Sharpe Ratio of all generated: 0.4273953626111349
SubSet 56 average Sharpe Ratio: 0.2434897313837182
Average Sharpe Ratio of all generated: 0.42416894802819777
SubSet 57 average Sharpe Ratio: 0.415846466965053
Average Sharpe Ratio of all generated: 0.4240254569753849
SubSet 58 average Sharpe Ratio: 0.4082445688772654
Average Sharpe Ratio of all generated: 0.4237579842957558
SubSet 59 average Sharpe Ratio: 0.3839215999712577
Average Sharpe Ratio of all generated: 0.42309404455701416
SubSet 60 average Sharpe Ratio: 0.5040699500602644
Average Sharpe Ratio of all generated: 0.4244215184177232
SubSet 61 average Sharpe Ratio: 0.3039919439827901
Average Sharpe Ratio of all generated: 0.4224791059268372
SubSet 62 average Sharpe Ratio: 0.36259737906779366
Average Sharpe Ratio of all generated: 0.42152860232589995
SubSet 63 average Sharpe Ratio: 0.31408304541988047
Average Sharpe Ratio of all generated: 0.4198497654992434
SubSet 64 average Sharpe Ratio: 0.35246197722780886
Average Sharpe Ratio of all generated: 0.4188130302950675
SubSet 65 average Sharpe Ratio: 0.2848564854560174
Average Sharpe Ratio of all generated: 0.416783385676294
SubSet 66 average Sharpe Ratio: 0.43903855864761104
Average Sharpe Ratio of all generated: 0.4171155524370599
SubSet 67 average Sharpe Ratio: 0.5605391929856814
Average Sharpe Ratio of all generated: 0.41922472362159846
SubSet 68 average Sharpe Ratio: 0.493239760652098
Average Sharpe Ratio of all generated: 0.42029740531769266
SubSet 69 average Sharpe Ratio: 0.5184072461552295
Average Sharpe Ratio of all generated: 0.42169897447251464
SubSet 70 average Sharpe Ratio: 0.35468130482298393
Average Sharpe Ratio of all generated: 0.42075506363238036
SubSet 71 average Sharpe Ratio: 0.4091833415929928
Average Sharpe Ratio of all generated: 0.4205943452707222
SubSet 72 average Sharpe Ratio: 0.56342811348888
Average Sharpe Ratio of all generated: 0.42255097223261473
SubSet 73 average Sharpe Ratio: 0.3876100574708974
Average Sharpe Ratio of all generated: 0.42207879770880774
SubSet 74 average Sharpe Ratio: 0.3618673627011537
Average Sharpe Ratio of all generated: 0.42127597857537236
SubSet 75 average Sharpe Ratio: 0.5777980510678504
Average Sharpe Ratio of all generated: 0.42333547952922074
SubSet 76 average Sharpe Ratio: 0.4328994798954758
Average Sharpe Ratio of all generated: 0.42345968732618505
SubSet 77 average Sharpe Ratio: 0.47812302634572007
Average Sharpe Ratio of all generated: 0.424160499364897
SubSet 78 average Sharpe Ratio: 0.4764553233318796
Average Sharpe Ratio of all generated: 0.42482245916194744
SubSet 79 average Sharpe Ratio: 0.32957575041323806
Average Sharpe Ratio of all generated: 0.4236318753025886
SubSet 80 average Sharpe Ratio: 0.44673295423458775
Average Sharpe Ratio of all generated: 0.42391707380792193
SubSet 81 average Sharpe Ratio: 0.2903969020582399
Average Sharpe Ratio of all generated: 0.42228877903048684
SubSet 82 average Sharpe Ratio: 0.22637733859733908
Average Sharpe Ratio of all generated: 0.41992840023008743
SubSet 83 average Sharpe Ratio: 0.37291430272985204
Average Sharpe Ratio of all generated: 0.4193687085931799
SubSet 84 average Sharpe Ratio: 0.4758754938895337
Average Sharpe Ratio of all generated: 0.4200334943025487
SubSet 85 average Sharpe Ratio: 0.3695486159580374
Average Sharpe Ratio of all generated: 0.4194464608334265
SubSet 86 average Sharpe Ratio: 0.35843180041903894
Average Sharpe Ratio of all generated: 0.418745142897629
SubSet 87 average Sharpe Ratio: 0.3680448562580317
Average Sharpe Ratio of all generated: 0.41816900327672446
SubSet 88 average Sharpe Ratio: 0.3652359460862152
Average Sharpe Ratio of all generated: 0.41757424982514574
SubSet 89 average Sharpe Ratio: 0.32706073577912537
Average Sharpe Ratio of all generated: 0.4165685441135233
SubSet 90 average Sharpe Ratio: 0.48103209487848725
Average Sharpe Ratio of all generated: 0.41727693478127015
SubSet 91 average Sharpe Ratio: 0.39665342617983185
Average Sharpe Ratio of all generated: 0.41705276620951537
SubSet 92 average Sharpe Ratio: 0.12044458123320058
Average Sharpe Ratio of all generated: 0.4138634308871894
SubSet 93 average Sharpe Ratio: 0.3900591522216232
Average Sharpe Ratio of all generated: 0.41361019388010895
SubSet 94 average Sharpe Ratio: 0.4309498445910629
Average Sharpe Ratio of all generated: 0.4137927165191716
SubSet 95 average Sharpe Ratio: 0.39312633042930073
Average Sharpe Ratio of all generated: 0.4135774416640688
SubSet 96 average Sharpe Ratio: 0.5220227422686242
Average Sharpe Ratio of all generated: 0.4146954344538065
SubSet 97 average Sharpe Ratio: 0.47234790254148634
Average Sharpe Ratio of all generated: 0.4152837249444971
SubSet 98 average Sharpe Ratio: 0.4339074541663541
Average Sharpe Ratio of all generated: 0.4154718434214855
SubSet 99 average Sharpe Ratio: 0.4912775785220973
Average Sharpe Ratio of all generated: 0.4162299007724916

古典アルゴリズムによる探索#

次にQUBO行列の準備及び古典なアルゴリズムに必要な関数を実装します。ここで、ペアワイズ相関行列は論文に従い対数リターンで計算を行います。

import heapq
import pandas as pd
import copy

#Attractiveness Coversion
def SRBucket(SR_list):
    Buckets=sorted(SR_list)
    Buckets.reverse()
    GroupedList = list(np.array_split(Buckets,11))
    for i in range(len(SR_list)):
        if   SR_list[i] in GroupedList[0]: SR_list[i]=15
        elif SR_list[i] in GroupedList[1]: SR_list[i]=12
        elif SR_list[i] in GroupedList[2]: SR_list[i]=9
        elif SR_list[i] in GroupedList[3]: SR_list[i]=6
        elif SR_list[i] in GroupedList[4]: SR_list[i]=3
        elif SR_list[i] in GroupedList[5]: SR_list[i]=0
        elif SR_list[i] in GroupedList[6]: SR_list[i]=-3
        elif SR_list[i] in GroupedList[7]: SR_list[i]=-6
        elif SR_list[i] in GroupedList[8]: SR_list[i]=-9
        elif SR_list[i] in GroupedList[9]: SR_list[i]=-12
        elif SR_list[i] in GroupedList[10]: SR_list[i]=-15

#Penalty/Reward Conversion 
def CorrelationBucket(Corr):
    for i in range(len(Corr)):
        for j in range(len(Corr)):
            if Corr[i][j] >= -1.00 and  Corr[i][j] < -0.25: Corr[i][j] = -5
            elif Corr[i][j] >= -0.25 and  Corr[i][j] < -0.15: Corr[i][j] = -3
            elif Corr[i][j] >= -0.15 and  Corr[i][j] < -0.05: Corr[i][j] = -1
            elif Corr[i][j] >= -0.05 and  Corr[i][j] < 0.05: Corr[i][j] = 0
            elif Corr[i][j] >= 0.05 and  Corr[i][j] < 0.15: Corr[i][j] = 1
            elif Corr[i][j] >= 0.15 and  Corr[i][j] < 0.25: Corr[i][j] = 3
            elif Corr[i][j] >= 0.25 and  Corr[i][j] < 1.00: Corr[i][j] = 5

#Ising component for classical algorithms
def hi(SR_list, Corr, i):
    h = 0.5*SR_list[i] + np.sum(Corr[i])
    return h

def jij(Corr, i, j):
    return 1./4.*Corr[i][j]

#Create Pairwise Correlation Matrix
def CreateCorrMat(Chart):
    assets = list()
    for iasset in range(len(Chart)):
        returns = list()
        for month in range(12):
            log_return = np.log(Chart[iasset][month+1]/Chart[iasset][month])
            returns.append(log_return)
        assets.append(returns)
    Chart_pd = pd.DataFrame(assets).T
    pairwise_corr = Chart_pd.corr(method='pearson')
    return pairwise_corr

#Initial state for Greedy Search and Genetic Algorithm
def GenerateRandomSolution(Nassets):
    Solution = np.random.randint(2, size=Nassets)
    for i in range(Nassets):
        Solution[i] = 2*Solution[i] - 1
    return Solution

#得られた解のポートフォリオのSharpe Ratioを計算する
def EvaluateSolution(Solution,Chart):
    selected_assets = list()
    for i in range(len(Solution)):
        if Solution[i] == 1:
            selected_assets.append(Chart[i])
    portfolioChart = np.mean(selected_assets, axis=0)
    portfolioSR = CalculateSharpeRatio(portfolioChart)
    return portfolioSR

#遺伝アルゴリズム、突然変異を導入して子孫を作る
def CreateDescendant(Ancestor, Ndescendants, MaxMutation):
    n = 0
    index_list = range(len(Ancestor))
    Descendants = list()
    while n < Ndescendants:
        Nmutaion = np.random.randint(MaxMutation)
        Place_to_change = np.random.choice(index_list, size=Nmutaion, replace=False)
        Descendant = list()
        for place in Place_to_change:
            Descendant = copy.deepcopy(Ancestor)
            Descendant[place] = Ancestor[place] * -1
            Descendants.append(Descendant)
        n += 1
    return Descendants

#銘柄セットを再生成する
#ここからさっきのステップはこの銘柄セットをずっと使う
Nassets = 48

Chart = CreateAssets(Nassets)
PairwiseCorrMat = CreateCorrMat(Chart)
SR_list = list()

for asset in Chart:
    SR_list.append(CalculateSharpeRatio(asset))

#Bucketの翻訳
SRBucket(SR_list)
CorrelationBucket(PairwiseCorrMat)

#print(PairwiseCorrMat)　#ペアワイズ相関行列の確認

貪欲法による探索#

まずは貪欲法アルゴリズムによる探索を行います。ランダムに生成した初期状態から出発し、全体のパフォーマンスに対する影響を見ながら、ポートフォリオのSharpe Ratioを高くできそうな銘柄を選んでいきます。反復を行うことで、ランダムな初期状態からある一定の状態に落ち着くことが観測できます。

実行後、初期状態と貪欲法で得られた状態の比較を行います。このときの比較基準はポートフォリオのSharpe Ratioとします。このチュートリアルでは全銘柄のウエイトが同じというシナリオで計算を行なっているため、選ばれた銘柄の平均的なチャートを得てSharpe Ratioの計算を行っています。

また、一個前のコードブロックを複数回実行して銘柄セットをリセットすると、貪欲法の結果が初期状態よりも悪くなったり、ほぼ変わらなかったりする場合もあることがわかります。また最終的なリターンが初期状態に負けているが、Sharpe Ratioが大きいという最終状態も出現します。それらを比較すると、Sharpe Ratioが大きい場合のポートフォリオはチャートが滑らかで、変動が小さい(リスクが小さい)ことが分かります。

#Greedy Search

#初期状態生成
Solution = GenerateRandomSolution(Nassets)
print("Initial random selection:")
print(Solution, EvaluateSolution(Solution, Chart))
print("")

selected_charts = list()
for i in range(Nassets):
    if Solution[i] == 1:
        selected_charts.append(Chart[i])
portfolioChart = np.mean(selected_charts, axis=0)

#初期状態に基づいてエネルギーの初期化を行う
Energies = list()
for iasset in range(Nassets):
    h = hi(SR_list, PairwiseCorrMat, iasset)
    energyTuple = [-1*abs(h), h , iasset]
    Energies.append(energyTuple)

#貪欲サーチの実行
NGreedyLoop = 10
for i in range(NGreedyLoop):
    heapq.heapify(Energies)

    ntry = 0
    #print(Energies) #エネルギー変化を確認
    while(ntry < len(Energies)):
        x, e, i = heapq.heappop(Energies)
        if e > 0:
            Solution[i] = -1.
        else:
            Solution[i] = 1.
        for ie in Energies:
            n = ie[2]
            ie[1] = ie[1] + Solution[i]*(jij(PairwiseCorrMat, i, n) + jij(PairwiseCorrMat, n, i))
            ie[0] = -ie[1]
        ntry+=1


#貪欲サーチの出力
print("Selection After Greedy Search:")
print(Solution, EvaluateSolution(Solution, Chart))
print("")

selected_charts2 = list()
for i in range(Nassets):
    if Solution[i] == 1:
        selected_charts2.append(Chart[i])
portfolioChart2 = np.mean(selected_charts2, axis=0)

plt.plot(list(range(13)), portfolioChart,color="r",label="Before Greedy Search")
plt.plot(list(range(13)), portfolioChart2,color="b",label="After Greedy Search")
plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()

Initial random selection:
[-1 -1  1 -1  1  1  1 -1 -1 -1 -1 -1  1 -1 -1  1  1 -1  1  1 -1 -1  1 -1
  1 -1 -1 -1 -1  1 -1 -1 -1  1 -1  1 -1  1  1 -1  1  1 -1  1  1  1  1  1] 1.761878543085147

Selection After Greedy Search:
[-1 -1 -1  1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1
  1  1 -1 -1  1  1 -1 -1  1 -1  1  1 -1 -1 -1  1 -1  1 -1  1 -1 -1 -1  1] 1.294291291375176

../../../_images/2d7c96b9c27673d4eadaea84b75f19877819c9c0f956656fdeb02e2821e8347a.png

遺伝的アルゴリズム#

貪欲法のほかに使われた古典アルゴリズムとして、遺伝的アルゴリズムがあります。ランダムに生成した複数のポートフォリオからパフォーマンスが良いものを選択し、その良いポートフォリオをもとにして次の世代のポートフォリオを生成します。次の世代のポートフォリオは、前の世代のポートフォリオにランダムな変更を加える突然変異と、ポートフォリオ間の遺伝子を交換させる操作などから生成されます。ここでは銘柄数が比較的小さいため、突然変異のみを考えて、遺伝子間の交換などは行わないものとします。

結果を比較することで、最後に計算された世代のパフォーマンスは第一世代よりも良くなっていることがわかります。ただし、ここでは試行した世代数が限られているため、パフォーマンスが改善されないまたは改善が小さい場合もあります。しかし、総じて貪欲法よりも安定な出力が得られることも観察できます。

注意：遺伝子の組数や子孫数そして変異の数を増やした場合、計算量の増加により計算が終わらない場合もあります。実行環境に合わせて設定を変更してください。

#Genetic Algorithm
Ngenes = 10 #遺伝子の組数
NGenerations = 5 #世代数
#ランダム初期化
Genes = list()
for i in range(Ngenes):
    gene = GenerateRandomSolution(Nassets)
    geneSR = EvaluateSolution(gene, Chart)
    Genes.append([geneSR, gene])
Genes = sorted(Genes,reverse=True)

#第一世代最優良ポートフォリオ
print("Best Gene 1st generation:")
print(Genes[0][1], EvaluateSolution(Genes[0][1], Chart))
print("")

selected_charts3 = list()
for i in range(Nassets):
    if Genes[0][1][i] == 1:
        selected_charts3.append(Chart[i])
portfolioChart3 = np.mean(selected_charts3, axis=0)


#遺伝的アルゴリズム実行
K_best = 5 #上位抽出遺伝子組数
Ndescendants = 2 #子孫数
MaxMutation = 2 #最大変異数
for iIter in range(NGenerations):
    selected_Genes = Genes[:K_best]
    for gene in selected_Genes:
        Descendants = CreateDescendant(gene[1], Ndescendants, MaxMutation)
        for Descendant in Descendants:
            DescendantSR = EvaluateSolution(Descendant, Chart)
            selected_Genes.append([DescendantSR, Descendant])
    selected_Genes = sorted(selected_Genes, key=lambda x: x[0], reverse=True)
    selected_Genes = copy.deepcopy(selected_Genes[:K_best])

#ラスト世代最優良ポートフォリオ
print("Best Gene last generation:")
print(selected_Genes[0][1], EvaluateSolution(selected_Genes[0][1], Chart))
print("")
selected_charts4 = list()
for i in range(Nassets):
    if selected_Genes[0][1][i] == 1:
        selected_charts4.append(Chart[i])
portfolioChart4 = np.mean(selected_charts4, axis=0)

#比較プロットを作る
plt.plot(list(range(13)), portfolioChart,color="r",label="Before Greedy Search")
plt.plot(list(range(13)), portfolioChart2,color="b",label="After Greedy Search")
plt.plot(list(range(13)), portfolioChart3,color="g",label="1st Gen Genetic")
plt.plot(list(range(13)), portfolioChart4,color="y",label="Last Gen Genetic")
plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()

Best Gene 1st generation:
[ 1  1  1 -1  1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1  1 -1 -1 -1  1  1  1  1
  1  1 -1 -1 -1 -1 -1 -1  1 -1 -1  1 -1 -1 -1 -1  1 -1 -1  1 -1 -1  1  1] 2.2242795097842576

Best Gene last generation:
[ 1 -1  1 -1  1 -1  1 -1 -1 -1 -1 -1 -1  1 -1 -1  1 -1 -1 -1  1  1  1  1
  1  1 -1 -1 -1 -1  1 -1  1 -1 -1  1 -1 -1 -1 -1 -1 -1 -1  1 -1 -1  1 -1] 2.7340523817811446

../../../_images/9c9707869bcccbfa671f828d69b851f579762d8763b969f5d0f6c1a55ce283c1.png

OpenJijを用いた量子アニーリング手法の実装#

Quantum Annealingによる解法#

Reverse Quantum Annealingでこの最適問題解くために、古典なアルゴリズムで得られたスピン形式の結果をアニーリングが使う形式に変換しましょう。

def ConvertSolutionToQuboState(solution):
    output = list()
    for i in range(len(solution)):
        if solution[i] == 1:
            output.append(1)
        else: output.append(0)
    return output

QA_init_state = ConvertSolutionToQuboState(Solution) #貪欲サーチの結果を初期状態として使う
#QA_init_state = ConvertSolutionToQuboState(selected_Genes[0][1]) #遺伝的アルゴリズムラスト世代の最優良ポートフォリオを初期状態として使う
#QA_init_state = GenerateRandomSolution(Nassets) #ランダムな初期状態を使う
print(QA_init_state)
    

[0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1]

Forward Annealingの場合#

まずはForward Annealingを実行してみましょう。銘柄の魅力度とペアワイズ相関による罰金と賞金の度合を使ってQUBO行列を作ります。以下では量子アニーリングを古典計算機上でシミュレートするSimulated Quantum Annealing (SQA)を用いて計算を行います。最終リターン率の改善は場合によりますが、貪欲法で得られた結果よりも安定したチャートを得られていることがわかります。

from openjij import SQASampler
sampler = SQASampler()

QUBO = np.random.rand(Nassets**2).reshape(Nassets, Nassets)
for i in range(Nassets):
    for j in range(Nassets):
        QUBO[i][j] = PairwiseCorrMat[i][j]
for i in range(Nassets):
    QUBO[i][i] = QUBO[i][i] + SR_list[i]

print(type(QUBO))
import matplotlib.pyplot as plt
plt.imshow(QUBO)
plt.colorbar()
plt.show()

sampleset_FQA = sampler.sample_qubo(QUBO,num_reads=10)
print(sampleset_FQA.record)
print(sampleset_FQA.record[0][0], EvaluateSolution(sampleset_FQA.record[0][0], Chart))
selected_charts = list()
for i in range(Nassets):
    if sampleset_FQA.record[0][0][i]:
        selected_charts.append(Chart[i])
portfolioChart_FQA = np.mean(selected_charts, axis=0)
plt.plot(list(range(13)), portfolioChart_FQA, color="c", label="Forward Annealing")
plt.plot(list(range(13)), portfolioChart,color="r",label="Before Greedy Search")
plt.plot(list(range(13)), portfolioChart2,color="b",label="After Greedy Search")
plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()

<class 'numpy.ndarray'>
[([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)]
[1 0 0 1 1 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 1 0 1 1 1 1 1 0 1 1 0 0 1 0 1 1 0
 0 0 1 0 1 1 1 0 0 0 1] 3.3954993443502355

../../../_images/39f1da2c41a9f12b7278cb46bd94a277938b318f5dd88d8e3f2c3667ba795895.png

../../../_images/eed689419eb2527cf1136b9b2e2bc08dc4fd056c91f2f25aa7ec14cff296143a.png

Reverse Quantum Annealing(RQA)の場合#

Reverse Annealingの場合、まずは古典的アルゴリズムから得られた状態を初期状態として入力してからアニーリングを行います。まずはアニーリングのスケジュールの設定です。ここではOpenJijを用いてSQAのハミルトニアンは次のように定めます。

\mathcal{H}_{\mathrm{SQA}}(s) = s\mathcal{H}_{p}-(1-s)\sum_{i}\sigma_{i}^{x}, 0\leq x\leq 1

ここで $\mathcal{H}_{p}$ は最適解を求めたい問題ハミルトニアンで、 $s$ と $1-s$ はそれぞれ参考論文[2]の $B[t]$ と $A[t]$ に対応します。この $s$ は

sample_qubo(QUBO, num_reads=10)

に追加の引数を指定することによって設定できます。初期状態の設定や更新を含めて、先程のsample_qubo部分を

sample_qubo(QUBO, schedule = user_schedule, initial_state = user_initial_state, num_reads=10, reinitialize_state = False)

のように変更すれば良いでしょう。ここのscheduleはアニーリングスケジュールを指しており

user_schedule = [
    [0, 0.1, 4], #s=0, beta=0.1, 4 steps
    [0.5, 1, 3],  #s=0, beta=1, 3 steps
    [1, 10, 3]  #s=1, beta=10, 3 steps
]

のような構造を持ちます。格納されたリストの[0]成分は $s$ の値で、アニーリング時の横磁場の強さを反映しています。[1]成分は逆温度 $\beta$ を表しています。大きければ低い温度を意味し、より基底状態に到達しやすくなります（デフォルトでは $\beta = 5$ で設定されています。）そして最後の[2]成分は量子モンテカルロ法のシミュレーションステップ数を表しており、各スケジュールの長さを調整できます。

initial_stateの引数に代入する初期状態を指定します。reinitialize_stateのオプションをFalseにすると、反復実行において前ステップでのアニーリング出力が次のステップのアニーリング初期状態として代入されます。上述のスクリプト例では、古典アルゴリズムで準備された初期状態から出発し、RQAを10回実行することになります。このオプション引数のデフォルトであるTrueでは、実行するたびに初期状態の再設定を行います。RQAを一回実行し、複数回の結果を比較したい場合はreinitialize_state = Trueにすると良いでしょう。

ではReverse Annealingのスケジュールを作成しましょう。OpenJijのsamplerは指定された $s$ と $\beta$ を指定したMC step数で量子モンテカルロを実行するため、短い定数をつなぎ合わせて $s$ や $\beta$ が少しずつ変化するようなスケジュールを作り、RQAを動作させる必要があります。ここでは $s$ の変化に注目してそれをプロットすると逆さまの台形のスケジュールを確認できます。横軸は量子モンテカルロ法のモンテカルロステップ数を表します。スケジュール作成関数の中身を書き換えれば、そのほかの形のスケジュールを自在に作成することができます。

def ScheduleFunction(x, RQAschedule):
    for step in RQAschedule:
        length = step[2]
        if x - length > 0:
            x = x-length
            continue
        else:
            s = step[0]
            return s


def ConvertScheduleToPlot(RQAschedule):
    TotalLength = np.sum(RQAschedule, axis=0)[2]
    x = np.arange(0, TotalLength+1)
    s = [ScheduleFunction(n, RQAschedule) for n in x]
    plt.plot(x, s ,label="s")
    plt.xlabel("MC step", loc="right")
    plt.ylabel("$s$, progress of quantum annealing", loc="top")
    plt.legend()
    plt.show()

#Create RQA schedule
RQAschedule = []
NReverseStep = 50
NPauseStep = 50
NForwardStep = 50
NMCStep = 20
TargetS = 0.38
ReverseStep = (1.0 - TargetS) / NReverseStep
ForwardStep = (1.0 - TargetS) / NForwardStep
beta = 5.
#Reverse Step
for i in range(NReverseStep):
    step_sche = [1.0-i*ReverseStep, beta, NMCStep]
    RQAschedule.append(step_sche)   

#Pause Step
RQAschedule.append([TargetS, beta, NPauseStep*NMCStep])

#Forward Step
for i in range(NForwardStep):
    step_sche = [TargetS+(i+1)*ForwardStep, beta, NMCStep]
    RQAschedule.append(step_sche)

#Plot Annealing Schedule
ConvertScheduleToPlot(RQAschedule)

../../../_images/859d0309ff9d42cdce88a448869326abd20fa34f2ccfde170b3727456dbaffe4.png

このようにして作成されたスケジュールと初期状態を指定したsample_qubo関数を用いて、RQAを実行しましょう。同時に貪欲法の結果やForward Annealingの結果の比較を行います。結果を見ると、このチュートリアルが行う設定では銘柄セットのサイズ制限によりForward Annealingと同じ最適解にたどり着く場合が多いことがわかります。

#初期状態を準備
init_state = QA_init_state

#Reverse Annealing
sampleset_RQA = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=False)
print(sampleset_RQA.record)

#チャートの比較
selected_charts = list()
for i in range(Nassets):
    if sampleset_RQA.record[0][0][i]:
        selected_charts.append(Chart[i])
portfolioChart_RQA = np.mean(selected_charts, axis=0)
plt.plot(list(range(13)), portfolioChart_FQA, color="c", label="Forward Annealing")
plt.plot(list(range(13)), portfolioChart_RQA, color="m", label="Reverse Annealing")
plt.plot(list(range(13)), portfolioChart,color="r",label="Before Greedy Search")
plt.plot(list(range(13)), portfolioChart2,color="b",label="After Greedy Search")
plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()
#結果の比較
print("Forward Quantum Annealing Result:")
print(sampleset_FQA.record[0][0], EvaluateSolution(sampleset_FQA.record[0][0], Chart))
print("Reverse Quantum Annealing Result:")
print(sampleset_RQA.record[0][0], EvaluateSolution(sampleset_RQA.record[0][0], Chart))

[([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)]
Forward Quantum Annealing Result:
[1 0 0 1 1 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 1 0 1 1 1 1 1 0 1 1 0 0 1 0 1 1 0
 0 0 1 0 1 1 1 0 0 0 1] 3.3954993443502355
Reverse Quantum Annealing Result:
[1 0 0 1 1 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 1 0 1 1 1 1 1 0 1 1 0 0 1 0 1 1 0
 0 0 1 0 1 1 1 0 0 0 1] 3.3954993443502355

../../../_images/c7f3a4623811288f9dfc265e6688777e736a3e1e5b4a4002579dc901fe7b4df7.png

Reverse Quantum Annealingの確認#

最終出力だけではRQAの振る舞いが分からないため、Phaseを分解してRQAを実行し、その挙動を見てみましょう。まずはReverse Phaseの確認です。同様にアニーリングスケジュールを設定し、スケジュールと初期状態を指定したアニーリングを行います。ここでは分かりやすくするためにReverse Phaseが止まる $s$ の値や逆温度の設定を変更しています。またアニーリングも反復実行するのではなく、毎回同じ初期状態から始まるように設定しました。結果を見ると、Reverse phaseの過程により、既に得られた解からよりエネルギーの高い状態に登っていることがわかります。また実行ごとに異なる状態で留まるため、高いエネルギー準位で解が定まっていない様子も伺えます。

#Create RQA schedule
RQAschedule = []
NReverseStep = 50
TargetS = 0.15 #よりランダムな状態に戻す
ReverseStep = (1.0 - TargetS) / NReverseStep
beta = 5.0
MC_step = 20
#Reverse Step
for i in range(NReverseStep):
    step_sche = [1.0-i*ReverseStep, beta, MC_step]
    RQAschedule.append(step_sche)   

ConvertScheduleToPlot(RQAschedule)
init_state = QA_init_state 
    

sampleset_RQA_Reverse = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
for state in sampleset_RQA_Reverse.record:
    selected_charts = list()
    for i in range(Nassets):
        if state[0][i]:
            selected_charts.append(Chart[i])
    portfolioChart = np.mean(selected_charts, axis=0)
    plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))

plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()
print(sampleset_RQA_Reverse.record)

../../../_images/1e8a9e72a368f5fd18635e2ff84b74b07749b9cdad30283589846a89432d0fd0.png

../../../_images/fda1ef818899436fe085e1cef94453edc2a0c85c3065ceda689421ccc7951a0b.png

[([0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -377., 1)
 ([1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -377., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -375., 1)
 ([0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -377., 1)
 ([0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -373., 1)
 ([1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -380., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -379., 1)]

次はPause Phaseの追加を行いましょう。追加後は同様にアニーリングを実行し、その結果を確認します。このPhaseの追加は実装上において上のReverse Phaseの一番最後のステップのMCステップ数を伸ばしたものに等しいため、上のReverseのみのケースと同じようにバラバラの結果となります。

#Create Pause Step
NPauseStep = 50
TargetS = 0.15 #よりランダムな状態に戻す
beta = 1.0 #より高い温度で熱的動きがしやすくようにする
MC_step = 20
step_sche = [TargetS, beta, MC_step*NPauseStep]
RQAschedule.append(step_sche)

#print(RQAschedule)
ConvertScheduleToPlot(RQAschedule)
init_state = QA_init_state 

sampleset_RQA_Reverse_Pause = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
for state in sampleset_RQA_Reverse_Pause.record:
    selected_charts = list()
    for i in range(Nassets):
        if state[0][i]:
            selected_charts.append(Chart[i])
    portfolioChart = np.mean(selected_charts, axis=0)
    plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))

plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()
print(sampleset_RQA_Reverse_Pause.record)

../../../_images/a27511f4bff0a1a5025571d03d9784c02a1e63fa4de2be62b7b1f281e1bc381f.png

../../../_images/cf08be5c2d700e0211d50f6cb18441f8b276a9369968b89de4f4931ad0eff9f5.png

[([1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1], -316., 1)
 ([0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1], -314., 1)
 ([0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1], -329., 1)
 ([0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1], -322., 1)
 ([0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -330., 1)
 ([1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -324., 1)
 ([0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1], -306., 1)
 ([1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1], -323., 1)
 ([1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1], -307., 1)
 ([0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1], -294., 1)]

最後にForward Phaseを追加し、アニーリングを実行します。この設定は先程のサンプルよりも少し緩い温度や磁場の条件で行いました。結果として、この銘柄セットで既知の量子アニーリングによる最適解が得られたことがわかります。RQAの節の最初の例と同様に、RQAによって最適解を発見できていることがわかります。

条件が少し緩いため解が収束せず、エネルギーが近い複数のsub-optimalの解の存在を確認することができます。それらの解の中には、場合によりForwardで得られた最適解よりも最終的なリターンが良く、形も比較的に滑らかなポートフォリオの存在も確認できます。そのポートフォリオの詳細を見ると、その大きなリターンは上昇率が大きい特定の銘柄によるもので、逆相関を持つ変動率が比較的に小さい銘柄と組み合わせて実現されたものだと分かります。実際の投資においてはそのようなポートフォリオが好ましいかもしれませんが、このチュートリアルが設定した「Shape Ratioがとにかく最大になる」条件では、相対的に小さいポートフォリSharpe Ratioが実現されたことから、量子アニーリングによって最適ではない解と判断されています。

#Create Forward Step
NForwardStep = 50
TargetS = 0.15 #よりランダムな状態に戻す
ForwardStep = (1.0 - TargetS) / NForwardStep
beta = 5.0 #より高い温度で熱的動きがしやすくようにする
MC_step = 20

#Forward Step
for i in range(NForwardStep):
    step_sche = [TargetS+(i+1)*ForwardStep, beta, MC_step]
    RQAschedule.append(step_sche)

ConvertScheduleToPlot(RQAschedule)
init_state = QA_init_state 

sampleset_RQA_Reverse_Pause_Forward = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
for state in sampleset_RQA_Reverse_Pause_Forward.record:
    selected_charts = list()
    for i in range(Nassets):
        if state[0][i]:
            selected_charts.append(Chart[i])
    portfolioChart = np.mean(selected_charts, axis=0)
    plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))

plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()
print(sampleset_RQA_Reverse_Pause_Forward.record)

../../../_images/c4cc120caf66c852f1e64908caa5aca5ddee0c983dfc47ed839beb414ca9acc4.png

../../../_images/e917d2c4b5a658b4920e877a13403946f9e0a8be0c5b4a5841b2967fda6aac15.png

[([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)
 ([1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1], -383., 1)]

RQAスケジュールを指定するときのパラメータ探索#

Phaseを分けてReverse Quantum Annealingの様子を確認するところで、スケジュールの設定に指定した $s$ , $\beta$ 、そして量子モンテカルロ法のステップ数によって結果が変化することが分かりました。このことから、Reverse Annealingでパフォーマンスを良くするには適切な値を設定する必要があるとわかります。

Reverse Phaseの場合#

まずはそれらがReverse Phaseに対する影響を見てみましょう。まずは量子アニーリングの度合を表す $s$ の影響から解析します。これは横磁場をどの程度強く戻すのかを表す量であり、最初のアニーリングで辿り着いた局所解の周りにあるポテンシャル障壁を越えて状態がばらけるようにするために、適切な強さ（ $s$ の適切な小ささ）を設定しなければなりません。このノートの結果では $s=0.20$ まではReverse Phaseの終状態のばらつきがあっても不定にはならず、 $s \leq 0.20$ では熱的揺らぎが支配的になり、ばらつきが大きくなる様子が分かります。よって $s \sim 0.2$ の近傍をさらに細かくスキャンすれば、最適の $s$ 値が得られると予想されます。

for TargetS in reversed(np.arange(0.05, 1.0, 0.05)):
    #Create RQA schedule
    RQAschedule = []
    NReverseStep = 50
    ReverseStep = (1.0 - TargetS) / NReverseStep
    beta = 5.0
    MC_step = 20
    #Reverse Step
    #for i in range(NReverseStep):
    for i in range(NReverseStep):
        step_sche = [1.0-i*ReverseStep, beta, MC_step]
        RQAschedule.append(step_sche)   

    init_state = QA_init_state 


    sampleset_RQA_Reverse = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
    for state in sampleset_RQA_Reverse.record:
        selected_charts = list()
        for i in range(Nassets):
            if state[0][i]:
                selected_charts.append(Chart[i])
        portfolioChart = np.mean(selected_charts, axis=0)
        plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))
    
    plt.xlabel("Month", loc="right")
    plt.ylabel("Portofolio Asset Amount", loc="top")
    plt.text(0.2,1.4, "Target s="+'{:.2f}'.format(TargetS))
    plt.show()

../../../_images/e40e448ba492779c09f459f1d4039bb6a1a4fd69298b97b24f8a5aa78c98686d.png

../../../_images/abd6fe5a09dfe8fe1073f411e7c7e8ca9e772669d0ee1b9968cf9eb406d173f6.png

../../../_images/2118642ea1fab2d8c7f6c877c69f6b007cff869b9102ad7a48aea063e6758419.png

../../../_images/4f96ee1039d0758d26e7edff2e05ce2c961628ba6d29d01fbe9e3b086530f179.png

../../../_images/91faf28821ae174de33925bf36cce76b7fd24f3753f6a0f53a49a336fe186709.png

../../../_images/08d4f8fecfee72194cc591ae20b06c2aff6931ec61e3f9c45bc9fc0b99ab1a61.png

../../../_images/008f0fd262e464818b33dca9147f588154979fbdb1372efecc5d04a4423f92b3.png

../../../_images/926fe2ba6b94d30de8caf000a53eb2ed1731d59091df6fefda863252a2b7de75.png

../../../_images/3de18a208bf49b1d2f0cced830a160a8a389fb311ef6addf47de71fef469f15a.png

../../../_images/dd82ed916bea673c26393de79bbf9b678aa92397a0b2025c50bdd62cbb84778e.png

../../../_images/c1bb2c1dab08cdaaeb2c6324ddad33093ccf008254d95c68e1c9fc4fae44ec3a.png

../../../_images/8139e50bd89d20dd156f970e3129984c2478113e430c5f203f772e492817c783.png

../../../_images/64ec09c2db76fe1d943e848e4941f9f767df2ac87941f072338b4ff4e829efae.png

../../../_images/798ab4779144180eff6ff5b0ba1069842cc3b6d25000367a9cee84c8bb3535be.png

../../../_images/6d03fce109ea23db131695e6d73f5d752d2fa59b557b121886c62ddd03d8fa86.png

../../../_images/5e23382a1b6f4ed22bb500c95ee32867dd0cac1cd4eabcd63bd025fe1202d312.png

../../../_images/06d1eb40c2650005d120007af230e4c3ff7079e68d536aea9eb18bda04f71cfd.png

../../../_images/e6107541dfb652a01884e03bea0c9f6b08f41d5480d2d20db0c193a91c4e0d00.png

../../../_images/f2a5949146617db170093ce02421e91c439558eda0042184a1a338962bac343a.png

次に $\beta$ による影響を調べてみましょう。 $\beta = \frac{1}{k_{\mathrm{B}}T}$ は逆温度を表しており、 $\beta$ が大きければ大きいほど系がより低温であることを意味します。 $\beta$ が大きいほど全体としては状態がより基底状態に落ちやすく、より「Annealed」な状態に近い振る舞いをします。言い換えると横磁場が一定の場合、 $\beta$ が大きくなると横磁場が効きにくくなり、実効的には大きい $s$ に近い振る舞いとなります。OpenJijにおいてその影響の詳細を知りたい方はこのQiita記事を参照すると良いでしょう。

また、その影響はQUBO行列に存在するエネルギーの典型的なスケール付近で変わることがわかっています。このチュートリアルの場合、QUBOに入っている典型な値の絶対値は大体 $5$ 程度で収まるため、 $\beta=5.0$ 程度で変化が起こると予想されます。

plt.hist((QUBO).flatten())
plt.xlabel("Values in QUBO", loc="right")
plt.ylabel("Number of entries", loc="top")

Text(0, 1, 'Number of entries')

../../../_images/36ce6d112eab4d3b0f0d52901b14e9c2c0ac85f3072cd7a812a04583b2ed53c4.png

対数的な $\beta$ 配列を用意し、 $s$ 以外の量は固定してReverse Phaseを行います。予想通り $\beta=5$ あたりで終状態のばらつきが小さくなり、それよりも大きい $\beta$ だと数少ない状態に収束してしまう様子がわかります。

for beta in [0.1,0.5,1.0,2.5,5.0,10.0,20.0,50.0]:
    #Create RQA schedule
    RQAschedule = []
    NReverseStep = 50
    TargetS = 0.18
    ReverseStep = (1.0 - TargetS) / NReverseStep
    MC_step = 20
    #Reverse Step
    #for i in range(NReverseStep):
    for i in range(NReverseStep):
        step_sche = [1.0-i*ReverseStep, beta, MC_step]
        RQAschedule.append(step_sche)   

    init_state = QA_init_state 


    sampleset_RQA_Reverse = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
    for state in sampleset_RQA_Reverse.record:
        selected_charts = list()
        for i in range(Nassets):
            if state[0][i]:
                selected_charts.append(Chart[i])
        portfolioChart = np.mean(selected_charts, axis=0)
        plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))
    
    plt.xlabel("Month", loc="right")
    plt.ylabel("Portofolio Asset Amount", loc="top")
    plt.text(0.2,1.8, "beta ="+'{:.2f}'.format(beta))
    plt.show()

../../../_images/aaa96b312ad95e6e5b29d6bb1307b8f8b1d778fe7b763152b911c4c362780885.png

../../../_images/a0a7d7995646a5786aaa9e7deead8f81349a59e898857f61e3b254f28a8ce5c7.png

../../../_images/73e6cd9c3cf728e71b20563f54c66776cbe4198ee71f6b09f12b236b88dfb751.png

../../../_images/0ea00c25dc5f30de9e984c6ed1914df4b2358920cdc3546c00bf7a01fc8801f8.png

../../../_images/d65c3fea322e3ed351ea5b3c9db70d559358c866eed7a3df960fe9174920806e.png

../../../_images/abbc40f410c464362c05ee0f661402e6c5d5cff88d0095cd766d253d70558e2c.png

../../../_images/7dcb1d37d167fe5cfeadd8973f467d41419ebfe2c96cabb1c8f2924368a97d6b.png

../../../_images/78b6ec413f57d8e95de7bd5c1e6c16f17620b743d0758469333f185c01689f97.png

そしてスケジュール分割の影響を調べましょう。ここでは階段関数のような極端的な場合から、細かく分割して緩やかな変化にする場合まで調査を行いました。結果から分かるように、分割数が少ない場合には状態があまり変化せず、逆に分割数を増やすとばらつきが見えてきます。また、大きく増やしても計算時間が長くなるだけで、結果の違いが小さいこともも分かります。このチュートリアルが取り扱うポートフォリオ最適化問題の場合は $50$ 程度の設定で問題ないこともわかります。

for NReverseStep in [1,2,5,10,15,20,50,100,500,1000]:
    #Create RQA schedule
    RQAschedule = []
    TargetS = 0.18
    ReverseStep = (1.0 - TargetS) / NReverseStep
    beta = 5.0
    MC_step = 20
    #Reverse Step
    #for i in range(NReverseStep):
    for i in range(NReverseStep):
        step_sche = [1.0-i*ReverseStep, beta, MC_step]
        RQAschedule.append(step_sche)   

    init_state = QA_init_state 


    sampleset_RQA_Reverse = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
    for state in sampleset_RQA_Reverse.record:
        selected_charts = list()
        for i in range(Nassets):
            if state[0][i]:
                selected_charts.append(Chart[i])
        portfolioChart = np.mean(selected_charts, axis=0)
        plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))
    plt.xlabel("Month", loc="right")
    plt.ylabel("Portofolio Asset Amount", loc="top")
    plt.text(0.2,1.8, "NReverseStep ="+'{:.2f}'.format(NReverseStep))
    plt.show()

../../../_images/7a528f35601e1ec244979111de339508d42a6216e25e418cab101dadb082e414.png

../../../_images/53329fac74d28003ae635d55acbefaedeeb9a2426a0d3ad9a9970f02d5017c39.png

../../../_images/20671ec23f171d3d2af12dc57b1c2179718f019c6e607ac729f2252871e7be28.png

../../../_images/ebb0e9032f0f189148f87e7b3d6a16dcb097be4b51e3a61a01d1cf2c69e07679.png

../../../_images/a0b180cd2827a932dfd3c340d583613c3d9fd1fb32abfc28f0ee0084a52ee00e.png

../../../_images/792fe8d1b3335ec03055b73a80249ee7bdac42cd1b90cc7b693abb60a8bb693d.png

../../../_images/80980d0ea9cad11072cd8a68cb0a470da125cfba882ad0d9e40531e442ee64fc.png

../../../_images/f5bd3275858335eacfda9c3e341ad1face7cdba7dd9ac4cdbb9c6aea31471dbe.png

../../../_images/daa9c1be04a56796bceaa6443bf44a4b4fcff2283328ec2d2fe88eb9c6577c52.png

../../../_images/a7d7d628323e3fc8631eef22f53bc9be78f950792590e0d5b7d1162d7d574964.png

またMC stepsに関しては、(分割数)×(MC steps)が極端に少ない場合でなければ、最終の結果にはほぼ影響がないことが分かりました。SQASampler()でデフォルト設定の場合では、目安として $1000$ 程度で問題なく動作するはずです。

Pause Phaseの場合#

既に説明したように、Pause PhaseはReverse Phaseの最後のステップを延長させたものであるため、アニーリング全体の $s$ と $\beta$ を変化させたところでReverse Phaseと同じようなふるまいをします。また、Pause phaseのステップ数はどのぐらい最後のステップを伸ばすのかと言う意味にに等しいため、変更してもほぼReverse Phaseの状態を保持したままであり、影響は少ないと予想されます。

ただし、Pause Phaseだけ $\beta$ を変化させることは量子アニーリングが不完全な状態で古典的なアニーリングを行うことに近い行動です。ここでは $\beta=5.0$ から緩やかに目標とした $\beta$ に向かって系の温度を変化させた場合の影響を調べましょう。Reverse Phaseで異なる $\beta$ でアニーリングを行った場合の終状態に近い結果を得ることができます。

for beta in [0.1,0.5,1.0,2.5,5.0,10.0,20.0,50.0]:
    #Create RQA schedule
    RQAschedule = []
    NReverseStep = 50
    TargetS = 0.18
    ReverseStep = (1.0 - TargetS) / NReverseStep
    MC_step = 20
    #Reverse Step
    #for i in range(NReverseStep):
    for i in range(NReverseStep):
        step_sche = [1.0-i*ReverseStep, 5.0, MC_step]
        RQAschedule.append(step_sche) 
    
    #Pause Phase
    NPauseStep = 50
    betaStep = (5.0-beta) / NPauseStep
    for i in range(NReverseStep):
        step_sche = [TargetS, 5.0-i*betaStep, MC_step]
        RQAschedule.append(step_sche) 

    init_state = QA_init_state 


    sampleset_RQA_Reverse_Pause = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
    for state in sampleset_RQA_Reverse_Pause.record:
        selected_charts = list()
        for i in range(Nassets):
            if state[0][i]:
                selected_charts.append(Chart[i])
        portfolioChart = np.mean(selected_charts, axis=0)
        plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))
    
    plt.xlabel("Month", loc="right")
    plt.ylabel("Portofolio Asset Amount", loc="top")
    plt.text(0.2,1.8, "Target beta ="+'{:.2f}'.format(beta))
    plt.show()

../../../_images/9c99dd11049d4197a5df05c61d9fbee1b587e344cc2944d3b5376d62224aa169.png

../../../_images/d660a053b0e77bf74b9b1c526401c02422f977fea3c1008724642413985fece7.png

../../../_images/f9799a11b8b54e1e796c9425faf8bb434aef9d89b86d5476d2e98b1685973068.png

../../../_images/744adedb0eb742bb7f9dcfe69a306807183a5121aabbb5083c07a9ca507b0ca0.png

../../../_images/782107f2a6e421014e4356c9a288f7dd13247ccf642c8ef69c3efdb86f66a4ac.png

../../../_images/c64c689854477ecfec94bc34923fbe5266bbbe2e4e8a79914fe994829d62fb8c.png

../../../_images/5c5bb5fa609ec34fdd6251c5b4b7f01809e6b948c72d4378d523d30a5e9fa749.png

../../../_images/5c0bea506caadf5262de584e55e72c3985276633fce1e8aacfed5619e25e8e6b.png

Forward Phaseの場合#

Forwar PhaseはReverse PhaseとPause Phaseを通して得られた状態から普通の量子アニーリングを行うフェーズです。このとき考えるぺきパラメータの条件は、通常のアニーリングが正しく動作するために必要な条件と基本的に同じです。Reverse PhaseとPause Phaseと同じ方法で $\beta$ と分割したステップ数のスキャンを行った結果、比較的小さい分割数と $\beta$ でも結果がちゃんと最適解に収束することが分かりました。

以下にコード例を示します。Reverse PhaseとPause Phaseにおいては、これまでの調査で得られた設定を用いており、Forward Phaseだけ粗い分解と小さい $\beta$ でアニーリングを行っています。

#Create RQA schedule
RQAschedule = []
NReverseStep = 50
TargetS = 0.18
ReverseStep = (1.0 - TargetS) / NReverseStep
MC_step = 20
#Reverse Step
#for i in range(NReverseStep):
for i in range(NReverseStep):
    step_sche = [1.0-i*ReverseStep, 5.0, MC_step]
    RQAschedule.append(step_sche) 
    
#Pause Phase
NPauseStep = 50
step_sche = [TargetS, 5.0, MC_step*NPauseStep]
RQAschedule.append(step_sche) 

#Forward Phase
NForwardStep = 5
ForwardStep = (1.0 - TargetS) / NForwardStep
for i in range(NForwardStep):
    step_sche = [TargetS+(i+1)*ForwardStep, 2.5, MC_step]
    RQAschedule.append(step_sche)

init_state = QA_init_state 

sampleset_RQA_Reverse_Pause_Forward = sampler.sample_qubo(QUBO, schedule=RQAschedule, initial_state = init_state, num_reads=10, reinitialize_state=True) #毎回同じ初期状態からアニーリング
for state in sampleset_RQA_Reverse_Pause_Forward.record:
    selected_charts = list()
    for i in range(Nassets):
        if state[0][i]:
            selected_charts.append(Chart[i])
    portfolioChart = np.mean(selected_charts, axis=0)
    plt.plot(list(range(13)), portfolioChart, label=("Energy="+str(state[1])))

plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.text(0.2,1.8, "NForwardStep ="+'{:.2f}'.format(NForwardStep))
plt.show()

../../../_images/7da1a3f822fba1642a933d28c34e6a786833507954d30ec2e3b865d2fccde6c4.png

終わりに#

OpenJijを用いてReverse Quantum Annealingを実装し、ポートフォリオ最適化問題を解いてみました。このチュートリアルでは、RQAの問題定式化やスケジュール指定のやり方、そしてRQAを行う際に考慮すべきパラメータなどの説明を行いました。RQAの手法を用いることで、正しく最適解を得られることを示しました。

最適解付近にsub-optimalな解が多く存在する他の問題に関しても、同様の手法でより良い解を求めることができます。通常の量子アニーリングで局所解に囚われた場合には、このようにReverse Quantum Annealingを試して見ると解決できるかもしれません。

また、例えばスケジュールの変化が線形関数から非線形な形にした場合の影響など、RQAについてこのチュートリアルで議論しなかった部分もありますが、このチュートリアルにあるコードを変更することそれらを簡単に実行できます。興味がある方は是非試してみてください。

付録A シミュレーションの結果が上手くいかない場合#

このチュートリアルで使用された銘柄セットの生成コードは、必ずしも構成が「良い」銘柄を生成する保証はありません。そのため、場合によって偏った分布を持ったセットに対して実験する可能性もあります。その時には、遺伝的アルゴリズム計算が非常に長い実行時間になることや、チュートリアルがデフォルトで使用した $s$ や $\beta$ などのパラメータでRQAが収束しないこともあり得ます。そうした場合は銘柄セットを再生成するか、手動でRQAのパラメータ調整を行って最適解に収束するようにすると良いでしょう。

付録B DwaveマシンでRQA実験をする場合のコード#

Dwaveマシンにおいてもこのチュートリアルが実行した手法でRQAを実行することができます。ここでは参考としてそのサンプルコードを掲載します。実装及びシミュレーションと実機の違いでスケジュールの設定方法や引数が異なりますが、基本的には同じような過程となります。詳細はDwave公式サイトなどにあるドキュメントを参考にすると良いでしょう。また実機では現実のノイズなどによる影響が見られやすいため、通常のアニーリングとRQAによる結果の違いは見られやすくなります。

from dwave.system import DWaveSampler, EmbeddingComposite

token = '*** your user token ***'
endpoint = '*** your dwave endpoint***' #'https://cloud.dwavesys.com/sapi/' as defaults

dw_sampler = DWaveSampler(solver='Advantage_system4.1', token=token) #Choose your dwave machine/simulator
sampler = EmbeddingComposite(dw_sampler)

Nassets = 48


# RQA schedule
timing = (0,1,2,3)#各フェーズの開始と終了時刻、フェーズの数を増やしても問題ない
Ratio = (1.0,0.38,0.38,1.0)#各フェーズが変わる時のs
schedule = list(zip(timing,Ratio))#スケジュールにまとめる
plt.plot(timing,Ratio)
plt.show()

#Create QUBO
QUBO = np.random.rand(Nassets**2).reshape(Nassets, Nassets)
for i in  range(Nassets):
    for j in range(Nassets):
        QUBO[i][j] = PairwiseCorrMat[i][j]
for i in range(Nassets):
    QUBO[i][i] = QUBO[i][i] + SR_list[i]

import matplotlib.pyplot as plt
plt.imshow(QUBO)
plt.colorbar()
plt.show()

#FQA
sampleset = sampler.sample_qubo(QUBO,num_reads=10)
print(sampleset.record)

min_forword = 0
for result in sampleset.record:
    if result[1] < min_forword:
        min_forword = result[1]
        best_forword = result[0]
selected_charts = list()
for i in range(Nassets):
    if best_forword[i]:
        selected_charts.append(Chart[i])
portfolioChart_FQA = np.mean(selected_charts, axis=0)
plt.plot(list(range(13)), portfolioChart_FQA,color="r",label="Forward Annealing")

#RQA
sampleset_RQA = sampler.sample_qubo(QUBO,num_reads=10,  anneal_schedule=schedule, initial_state = QA_init_state)
print(sampleset_RQA.record)

min_RQA = 0
for result in sampleset_RQA.record:
    if result[1] < min_RQA:
        min_RQA = result[1]
        best_RQA = result[0]
selected_charts = list()
for i in range(Nassets):
    if best_RQA[i]:
        selected_charts.append(Chart[i])
portfolioChart_RQA = np.mean(selected_charts, axis=0)
plt.plot(list(range(13)), portfolioChart_RQA,color="b",label="Reverse Annealing")
plt.xlabel("Month", loc="right")
plt.ylabel("Portofolio Asset Amount", loc="top")
plt.legend()
plt.show()

---------------------------------------------------------------------------
SolverAuthenticationError                 Traceback (most recent call last)
Cell In[21], line 6
      3 token = '*** your user token ***'
      4 endpoint = '*** your dwave endpoint***' #'https://cloud.dwavesys.com/sapi/' as defaults
----> 6 dw_sampler = DWaveSampler(solver='Advantage_system4.1', token=token) #Choose your dwave machine/simulator
      7 sampler = EmbeddingComposite(dw_sampler)
      9 Nassets = 48

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/dwave/system/samplers/dwave_sampler.py:185, in DWaveSampler.__init__(self, failover, retry_interval, **config)
    182 self._solver_penalty = defaultdict(int)
    184 self.client = Client.from_config(**config)
--> 185 self.solver = self._get_solver(penalty=self._solver_penalty)

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/dwave/system/samplers/dwave_sampler.py:200, in DWaveSampler._get_solver(self, refresh, penalty)
    198 filters = copy.deepcopy(self.client.default_solver)
    199 order_by = filters.pop('order_by', 'avg_load')
--> 200 solvers = self.client.get_solvers(refresh=refresh, order_by=order_by, **filters)
    202 # we now just need to de-prioritize penalized solvers
    203 solvers.sort(key=lambda solver: penalty.get(solver.id, 0))

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/dwave/cloud/events.py:105, in dispatches_events.__call__.<locals>.wrapped(*pargs, **kwargs)
    103 dispatch_event(self.before_eventname, obj=obj, args=args)
    104 try:
--> 105     rval = fn(*pargs, **kwargs)
    106 except Exception as exc:
    107     dispatch_event(self.after_eventname, obj=obj, args=args, exception=exc)

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/dwave/cloud/client/base.py:1173, in Client.get_solvers(self, refresh, order_by, **filters)
   1170     query['name'] = filters['name__eq']
   1172 # filter
-> 1173 solvers = self._fetch_solvers(**query)
   1174 solvers = [s for s in solvers if all(p(s) for p in predicates)]
   1176 # sort: undefined (None) key values go last

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/dwave/cloud/utils.py:489, in cached.__call__.<locals>.wrapper(*args, **kwargs)
    487     val = data.get('val')
    488 else:
--> 489     val = fn(*args, **kwargs)
    490     self.cache[key] = dict(expires=now+self.maxage, val=val)
    492 return val

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/dwave/cloud/client/base.py:807, in Client._fetch_solvers(self, name)
    804     url = 'solvers/remote/'
    806 try:
--> 807     data = Client._sapi_request(self.session.get, url)
    808 except SAPIError as exc:
    809     if name is not None and exc.error_code == 404:

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/dwave/cloud/client/base.py:1798, in Client._sapi_request(meth, *args, **kwargs)
   1796 else:
   1797     if response.status_code == 401:
-> 1798         raise SolverAuthenticationError(error_code=401)
   1800     try:
   1801         msg = response.json()

SolverAuthenticationError: Invalid token or access denied

参考文献#

Harry Markowitz, “Portfolio selection”, The journal of finance, 7(1):77–91 (1952)
Davide Venturelli, Alexei Kondratyev, “Reverse Quantum Annealing Approach to Portfolio Optimization Problems”, Quantum Machine Intelligence volume 1, pages17–30 (2019)
Sharpe, William F., “Mutual fund performance”, The Journal of Business 39 (1), 119-138 (1966)