Portfolio Optimization using Reverse Quantum Annealing

Introduction

This chapter presents the utilization of Reverse Quantum Annealing for resolving portfolio optimization issues through the utilization of OpenJij. Reverse Quantum Annealing is an approach that amalgamates classical optimization methodologies with quantum annealing. The chapter is structured in the following manner: We will elaborate on the methods of resolving portfolio optimization problems through classical algorithms, standard quantum annealing, and reverse quantum annealing. Subsequently, we contrast the outcomes obtained from these techniques. Throughout the chapter, we will describe the implementation of quantum annealing via OpenJij and accentuate any potential limitations or caveats of the numerical experiments.

Note: Note: Execution of this content requires OpenJij version 0.4.9 or higher. If necessary, please execute the following command to update OpenJij prior to proceeding with this notebook.

pip install -U openjij

Portfolio Optimization

In the realm of asset management, investors strive to achieve substantial returns while minimizing risk as much as feasible. Consequently, a prevalent approach is to construct a diversified investment portfolio by integrating assets with moderate returns but minimal (or no) risk and assets with substantial expected returns but also substantial risk.

In this scenario, an optimal portfolio must be chosen in order to achieve the maximal profit for a given level of risk or to minimize the risk while obtaining the same profit. One of the most widely utilized methodologies currently is the Modern Portfolio Theory proposed by Markowitz. This method takes into account the interdependencies among the stocks comprising the portfolio and aims to minimize their covariance [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$, and $R$ denote the weight of each stock in the portfolio, the covariance among the stocks, the expected return of each stock, and the total return of the portfolio, respectively.

Here, the Sharpe ratio is a widely-utilized metric for assessing the performance of a portfolio.

Sharpe Ratio

The Sharpe ratio, proposed by William Sharpe in 1966 [3], is a metric for evaluating the efficiency of an investment, applied to both individual stocks and portfolios. It is defined as follows

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

Here, $\bar{r}$ represents the average rate of return over a specified time period, typically calculated using the mean of monthly or annual returns. $r_0$ is the rate of return in the risk-free context, usually represented by the interest rate of a stable government bond. $\sigma$ denotes volatility, which is the degree of fluctuation in asset values.

The calculation of the Sharpe ratio varies depending on the situation, and the methodologies for determining return and volatility may include using simple returns calculated from asset price $P$ as $r=P(t+1)/P(t) - 1$, or log returns calculated as $r=log[P(t+1)/P(t)]$. Also volatility can be determined from price changes or from the rate of return. For a more comprehensive examination of these methodologies, please refer to literature and resources about financial engineering and market analysis.

Portfolio Evaluation and Optimization

As per the definition of the Sharpe ratio, a high value connotes high returns or low risk (volatility). Among portfolios constructed from a given set of stocks, the combination that attains the highest Sharpe ratio can be considered as the combination that yields the highest returns per a unit of risk. Thus, the problem of portfolio optimization can be rephrased as the problem of maximizing the Sharpe ratio under a given constraint.

In practice, this optimization process involves identifying stocks and determining the proportion of investment allocated to each stock. In other words, it is a combinatorial optimization problem that involves determining the selection of stocks and the weight assigned to each stock. For simplicity, this tutorial will adopt a strategy in which each stock is assigned an equal weight, as per the referenced paper [2]. We select $N$ stocks out of a total of $M$ stocks and find the combination that maximizes the Sharpe ratio. Then, the problem can be formulated into the QUBO form as follows.

$$\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}$$

Here, $q_{i}$ is a binary variable that denotes the selection of each stock, where a value of $1$ indicates inclusion in the portfolio and a value of $0$ indicates non-inclusion. $a_i$ is the score of attractiveness, derived from the performance of the individual stock, and $b_{ij}$ represents the fines or bonuses as determined by the pairwise correlation between stocks.

The specific values of $A_i$ and $B_{ij}$ are established in accordance with the Table 1 below.

The $a_i$ groups in this table are criteria established by the portfolio’s creators, with stocks ordered by their attractiveness. The table displays the groups formed when they are ranked and divided equally. The attractiveness criterion can also incorporate factors such as a high rate of return, but for this tutorial, we simply use the order of individual stocks based on the Sharpe ratio. $B_{ij}$ is determined by the value of the component $\rho_{ij}$ of the correlation matrix obtained from the time series of log returns.

By performing quantum annealing using this QUBO, we should be able to obtain the optimal portfolio. If multiple stocks are selected, a constraint on the number of stocks can be added. However, as described in the referenced paper [2], when constraints are added to the QUBO using this method, the large interaction (penalty) introduced can also affect the annealing performance and the final results obtained. In a real-world scenario, the number of stocks in the optimal portfolio is unknown, thus, this tutorial performs the same optimization without any restriction on the number of stocks.

Portfolio Optimization with Reverse Quantum Annealing

The problem of portfolio optimization in the real world is to select a subset of stocks from a market with a vast array of options. Additionally, due to the huge number of possible stock combinations, it is likely that there exist numerous local minima, which correspond portfolios with Sharpe ratios that are quite close to the maximum. Conventional algorithms, when trapped in the trough of such an energy landscape, need a significant amount of time and resources to escape. To circumvent the above issue and solve portfolio optimization problems more efficiently, computational methods such as quantum annealing, which can handle a large number of correlations and are less prone to getting trapped in local minima, are more effective.

However, as demonstrated in the reference paper [2], when only simple forward annealing is used, the time required to reach the optimal solution increases dramatically as the number of issues increases. This is because the difference in energy between the optimal solution and local minima with a Sharpe ratio close to the optimal solution is small. Due to a high potential barrier between these levels with almost no energy change, transition from local minima to optimal require a longer time. Furthermore, in actual quantum annealing machines, the system may be excited by thermal noise or other factors during standard annealing, resulting in a deviation from the intended time evolution and a final solution that is not truly optimal obtained.

Reverse Quantum Annealing

One of the proposed methods for addressing the challenges associated with using quantum annealing in portfolio optimization is reverse quantum annealing. As the name implies, this approach involves conducting quantum annealing in a manner that incorporates a step in the opposite direction of conventional forward annealing. The time-dependent Hamiltonian of standard forward annealing is represented as

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

where $A[t]$ is the magnitude of the transverse magnetic field applied to the system during annealing, and $B[t]$ represents the amplitude of the problem Hamiltonian $\mathcal{H}_{\chi-\mathrm{Ising}}$, respectively. As depicted in the following figure a), in the case of forward annealing, $A[t]$ is gradually decreased as the annealing progresses, and concurrently, $B[t]$ is made dominant. At the completion of annealing, $A[t] = 0$.

In contrast, reverse annealing starts from $A[t]=0, B[t]=1$, as depicted in Figure b). In this method, $B[t]$ is initially decreased, in contrast to the standard annealing. This is referred to as the reverse phase. In this phase, the system’s energy is inversely increased, making it easier to tunnel the potential barrier between A and B, as illustrated in c). After $A[t]$ and $B[t]$ reach a specific value, the system is kept in the same state and allowed to wait for a while. This is referred to as the pause phase. During this phase, a high energy state is maintained, increasing the likelihood of thermal hopping. The potential barrier on the right side of B is also easily overcome. This suggests that the possibility of escaping from the initial local minimum and reaching a near-optimal solution is increased. Finally, the optimal solution is attained by conducting standard annealing in the forward phase, which once again makes $A[t]$ smaller and $B[t]$ larger. In this protocol, the performance of the final solution may also be impacted by the value at which $A[t], B[t]$ is stopped in the reverse phase, the duration of each phase, etc.

Portfolio Optimization Procedure with Reverse Quantum Annealing

From the previous discussion, reverse annealing requires a reverse phase initial state. While it can be used randomly generated initial states, it is predicted that the efficacy of the optimal solution search will be enhanced when a specific local minumu state is implemented. The reference paper [2] employs the output generated by a conventional algorithm. This tutorial will follow the method outlined in the reference paper [2], which entails the following steps.

• Generation of data for stocks to be optimizaed

• Implementation of a classical algorithm and examination of its results

• Search for an optimal solution through standard forward annealing

• Search for an optimal solution through reverse quantum annealing

• Exploration of parameters for reverse quantum annealing

Generation of stock data and implementation of classical algorithms

Generation of stock data to be optimized

We generate a chart of stocks by Brownian motion using the given initial values according to the method in reference paper [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

The motion of the chart at each time, as outlined in Appendix A of reference paper [2], is derived from the motion at the preceding time, as follows:

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

where $Z_n$ is the multivariate normal distribution following the uniform correlation matrix ρ, constructed via the Cholesky decomposition.

We shall now implement and run this equation, and plot the chart to observe its appearance. By beginning with the same initial value and examining the 12-month trend, we can discern the overall behavior to spread out. Additionally, we can observe certain stocks note any significant increases or decreases.

# generate correlated normal probability variables, 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

# generate Brownian motion of the chart
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 # the number of stocks
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] # initial value: 1.0 at relative price

  # simulate 12 steps (12 months)
  for month in range(12):
    chart_asset.append(GetNextSt(chart_asset[month], mu, sigma, ZList[month][iasset]))

  chart.append(chart_asset)
  #print(chart_asset) # display specific values for each stocks
  plt.plot(list(range(13)), chart_asset)

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

Next, the Sharpe Ratio is calculated for each stock. Utilizing the realized Sharpe Ratio, we evaluate the earlier outcome. Following its definition, we determine the average return exceeded the risk-free rate of return, $\overline{r}$ and the standard deviation of the excess return rate $\sigma$. As plot of the chart , when the number of stocks is limited, there is a significant bias due to chance and existence of correlation among the stocks in generation condition. Therefore, we mitigate their impact by generating 100 iterations and assessing the average Shape ratio. Consequently, we can see that the average Sharpe ratio of the stocks is $0.4$. Additionally, when the definition of the Sharpe ratio calculation is changed, the Sharpe ratio will naturally vary, but it will generally remain within the range of $0.40 \pm 0.05$.

# function for generating assets
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

# function for computing standard 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)

# function for computing Sharpe ratio based on log-return
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

# function for computing Sharpe ratio with a slightly different definition
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): # multiple execution to eliminate the effect of correlations that are not in the same set
  Chart = CreateAssets(48) # create a specified number of assets
  mean_SR = 0.0
  n = 0.
  for asset in Chart:
    assetSR = CalculateSharpeRatio(asset) # choose the way of computation of Sharpe ratio
    #assetSR = CalculateSharpeRatio_log(asset)
    #assetSR = CalculateSharpeRatio_2(asset)
    mean_SR = ((mean_SR*n)+assetSR) / (n+1) # average of 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.41952730798198146
Average Sharpe Ratio of all generated: 0.41952730798198146
SubSet 1 average Sharpe Ratio: 0.5188854594500679
Average Sharpe Ratio of all generated: 0.4692063837160247
SubSet 2 average Sharpe Ratio: 0.3892844885364774
Average Sharpe Ratio of all generated: 0.4425657519895089
SubSet 3 average Sharpe Ratio: 0.26737662609469687
Average Sharpe Ratio of all generated: 0.3987684705158059
SubSet 4 average Sharpe Ratio: 0.5171695546459764
Average Sharpe Ratio of all generated: 0.42244868734183993
SubSet 5 average Sharpe Ratio: 0.306273107623704
Average Sharpe Ratio of all generated: 0.4030860907221507
SubSet 6 average Sharpe Ratio: 0.4160148427860553
Average Sharpe Ratio of all generated: 0.40493305530270846
SubSet 7 average Sharpe Ratio: 0.49292110118053906
Average Sharpe Ratio of all generated: 0.4159315610374373
SubSet 8 average Sharpe Ratio: 0.5074406259584451
Average Sharpe Ratio of all generated: 0.42609923491754925
SubSet 9 average Sharpe Ratio: 0.46631157693548514
Average Sharpe Ratio of all generated: 0.4301204691193428
SubSet 10 average Sharpe Ratio: 0.3743293171514758
Average Sharpe Ratio of all generated: 0.4250485462131731
SubSet 11 average Sharpe Ratio: 0.5338384819395325
Average Sharpe Ratio of all generated: 0.43411437419036975
SubSet 12 average Sharpe Ratio: 0.24616701533484453
Average Sharpe Ratio of all generated: 0.41965688504763704
SubSet 13 average Sharpe Ratio: 0.3382598699167572
Average Sharpe Ratio of all generated: 0.4138428125382885
SubSet 14 average Sharpe Ratio: 0.36530136964923837
Average Sharpe Ratio of all generated: 0.41060671634568513
SubSet 15 average Sharpe Ratio: 0.36662023256742665
Average Sharpe Ratio of all generated: 0.40785756110954396
SubSet 16 average Sharpe Ratio: 0.39786060950415325
Average Sharpe Ratio of all generated: 0.4072695051327563
SubSet 17 average Sharpe Ratio: 0.38338657251276187
Average Sharpe Ratio of all generated: 0.40594267554275665
SubSet 18 average Sharpe Ratio: 0.2756191797386283
Average Sharpe Ratio of all generated: 0.3990835441846446
SubSet 19 average Sharpe Ratio: 0.3868552012811342
Average Sharpe Ratio of all generated: 0.3984721270394691
SubSet 20 average Sharpe Ratio: 0.46993312200189935
Average Sharpe Ratio of all generated: 0.40187503156148957
SubSet 21 average Sharpe Ratio: 0.3032447143907356
Average Sharpe Ratio of all generated: 0.39739183532645533
SubSet 22 average Sharpe Ratio: 0.33469358355099516
Average Sharpe Ratio of all generated: 0.3946658243796962
SubSet 23 average Sharpe Ratio: 0.4694551685379134
Average Sharpe Ratio of all generated: 0.39778204705295533
SubSet 24 average Sharpe Ratio: 0.37218504664863317
Average Sharpe Ratio of all generated: 0.39675816703678246
SubSet 25 average Sharpe Ratio: 0.3818505148173172
Average Sharpe Ratio of all generated: 0.39618479579757226
SubSet 26 average Sharpe Ratio: 0.37082286358565036
Average Sharpe Ratio of all generated: 0.3952454649749085
SubSet 27 average Sharpe Ratio: 0.23708338156208067
Average Sharpe Ratio of all generated: 0.389596819138736
SubSet 28 average Sharpe Ratio: 0.4011810965944176
Average Sharpe Ratio of all generated: 0.3899962769820354
SubSet 29 average Sharpe Ratio: 0.31549901042040657
Average Sharpe Ratio of all generated: 0.3875130347633144
SubSet 30 average Sharpe Ratio: 0.3247377570455879
Average Sharpe Ratio of all generated: 0.38548802580467806
SubSet 31 average Sharpe Ratio: 0.47290754065242563
Average Sharpe Ratio of all generated: 0.38821988564367016
SubSet 32 average Sharpe Ratio: 0.5126995854425292
Average Sharpe Ratio of all generated: 0.3919919977587871
SubSet 33 average Sharpe Ratio: 0.3349977979677347
Average Sharpe Ratio of all generated: 0.39031569776493263
SubSet 34 average Sharpe Ratio: 0.24991519132194343
Average Sharpe Ratio of all generated: 0.38630425472370433
SubSet 35 average Sharpe Ratio: 0.5021355665729327
Average Sharpe Ratio of all generated: 0.3895217911639607
SubSet 36 average Sharpe Ratio: 0.43159588100836116
Average Sharpe Ratio of all generated: 0.3906589287273229
SubSet 37 average Sharpe Ratio: 0.4575539909386306
Average Sharpe Ratio of all generated: 0.3924193251013047
SubSet 38 average Sharpe Ratio: 0.2373653509680486
Average Sharpe Ratio of all generated: 0.3884435821748109
SubSet 39 average Sharpe Ratio: 0.28390137996064174
Average Sharpe Ratio of all generated: 0.3858300271194567
SubSet 40 average Sharpe Ratio: 0.2972174936099685
Average Sharpe Ratio of all generated: 0.38366874581434723
SubSet 41 average Sharpe Ratio: 0.41958902649055996
Average Sharpe Ratio of all generated: 0.3845239905923523
SubSet 42 average Sharpe Ratio: 0.4470866640087215
Average Sharpe Ratio of all generated: 0.38597893648575626
SubSet 43 average Sharpe Ratio: 0.29623936707558857
Average Sharpe Ratio of all generated: 0.3839394008173434
SubSet 44 average Sharpe Ratio: 0.2784424687757346
Average Sharpe Ratio of all generated: 0.3815950245497521
SubSet 45 average Sharpe Ratio: 0.5595418778979847
Average Sharpe Ratio of all generated: 0.3854634344051484
SubSet 46 average Sharpe Ratio: 0.29962483833237585
Average Sharpe Ratio of all generated: 0.38363708129721713
SubSet 47 average Sharpe Ratio: 0.34513372331701014
Average Sharpe Ratio of all generated: 0.3828349280059628
SubSet 48 average Sharpe Ratio: 0.3489474431121186
Average Sharpe Ratio of all generated: 0.38214334668159866
SubSet 49 average Sharpe Ratio: 0.468581329634542
Average Sharpe Ratio of all generated: 0.38387210634065755
SubSet 50 average Sharpe Ratio: 0.4892250105607534
Average Sharpe Ratio of all generated: 0.38593784956065935
SubSet 51 average Sharpe Ratio: 0.35818600932399797
Average Sharpe Ratio of all generated: 0.385404160325339
SubSet 52 average Sharpe Ratio: 0.32695883491146954
Average Sharpe Ratio of all generated: 0.384301418336398
SubSet 53 average Sharpe Ratio: 0.30571189507347873
Average Sharpe Ratio of all generated: 0.3828460567944921
SubSet 54 average Sharpe Ratio: 0.43869940090759973
Average Sharpe Ratio of all generated: 0.3838615721420031
SubSet 55 average Sharpe Ratio: 0.47122449077188405
Average Sharpe Ratio of all generated: 0.3854216242603939
SubSet 56 average Sharpe Ratio: 0.31792508635103
Average Sharpe Ratio of all generated: 0.3842374744725103
SubSet 57 average Sharpe Ratio: 0.26755589243437267
Average Sharpe Ratio of all generated: 0.38222572305805963
SubSet 58 average Sharpe Ratio: 0.436035953711601
Average Sharpe Ratio of all generated: 0.3831377608657468
SubSet 59 average Sharpe Ratio: 0.39784574563113595
Average Sharpe Ratio of all generated: 0.38338289394516994
SubSet 60 average Sharpe Ratio: 0.27758718173766755
Average Sharpe Ratio of all generated: 0.3816485380073421
SubSet 61 average Sharpe Ratio: 0.5005567232270803
Average Sharpe Ratio of all generated: 0.3835664119624992
SubSet 62 average Sharpe Ratio: 0.37534763227404544
Average Sharpe Ratio of all generated: 0.38343595514204754
SubSet 63 average Sharpe Ratio: 0.39595572153593483
Average Sharpe Ratio of all generated: 0.383631576491952
SubSet 64 average Sharpe Ratio: 0.2769412389044738
Average Sharpe Ratio of all generated: 0.3819901866829139
SubSet 65 average Sharpe Ratio: 0.32518582362698684
Average Sharpe Ratio of all generated: 0.38112951451539984
SubSet 66 average Sharpe Ratio: 0.4783506573888539
Average Sharpe Ratio of all generated: 0.382580576349332
SubSet 67 average Sharpe Ratio: 0.35857981628943714
Average Sharpe Ratio of all generated: 0.3822276239955101
SubSet 68 average Sharpe Ratio: 0.3218432278350593
Average Sharpe Ratio of all generated: 0.38135248781927167
SubSet 69 average Sharpe Ratio: 0.4026345921987826
Average Sharpe Ratio of all generated: 0.38165651788183613
SubSet 70 average Sharpe Ratio: 0.3103267773595729
Average Sharpe Ratio of all generated: 0.38065187364912817
SubSet 71 average Sharpe Ratio: 0.5710399960548177
Average Sharpe Ratio of all generated: 0.38329615312698495
SubSet 72 average Sharpe Ratio: 0.3756727977432663
Average Sharpe Ratio of all generated: 0.3831917236011806
SubSet 73 average Sharpe Ratio: 0.39117837718983095
Average Sharpe Ratio of all generated: 0.3832996513523786
SubSet 74 average Sharpe Ratio: 0.3432956830991021
Average Sharpe Ratio of all generated: 0.3827662651090015
SubSet 75 average Sharpe Ratio: 0.5402088003817597
Average Sharpe Ratio of all generated: 0.384837877415222
SubSet 76 average Sharpe Ratio: 0.5325663015423177
Average Sharpe Ratio of all generated: 0.3867564283779116
SubSet 77 average Sharpe Ratio: 0.24399214739846198
Average Sharpe Ratio of all generated: 0.38492611708330327
SubSet 78 average Sharpe Ratio: 0.36508745971334616
Average Sharpe Ratio of all generated: 0.3846749948381139
SubSet 79 average Sharpe Ratio: 0.48577088410347297
Average Sharpe Ratio of all generated: 0.3859386934539309
SubSet 80 average Sharpe Ratio: 0.3992112559568508
Average Sharpe Ratio of all generated: 0.3861025522502632
SubSet 81 average Sharpe Ratio: 0.48239592473034576
Average Sharpe Ratio of all generated: 0.387276861670752
SubSet 82 average Sharpe Ratio: 0.43089396930324037
Average Sharpe Ratio of all generated: 0.38780236899162535
SubSet 83 average Sharpe Ratio: 0.3435007363220112
Average Sharpe Ratio of all generated: 0.3872749686027014
SubSet 84 average Sharpe Ratio: 0.5110329701523618
Average Sharpe Ratio of all generated: 0.38873094509152095
SubSet 85 average Sharpe Ratio: 0.32874385087389524
Average Sharpe Ratio of all generated: 0.38803342074015323
SubSet 86 average Sharpe Ratio: 0.36240281381239775
Average Sharpe Ratio of all generated: 0.3877388160628227
SubSet 87 average Sharpe Ratio: 0.4122849607539785
Average Sharpe Ratio of all generated: 0.38801774952522217
SubSet 88 average Sharpe Ratio: 0.37921503154161357
Average Sharpe Ratio of all generated: 0.3879188425815861
SubSet 89 average Sharpe Ratio: 0.26905075564182396
Average Sharpe Ratio of all generated: 0.38659808606003315
SubSet 90 average Sharpe Ratio: 0.419817874621962
Average Sharpe Ratio of all generated: 0.3869631386815928
SubSet 91 average Sharpe Ratio: 0.548206430226296
Average Sharpe Ratio of all generated: 0.3887157831549048
SubSet 92 average Sharpe Ratio: 0.29546389974672704
Average Sharpe Ratio of all generated: 0.38771307473116096
SubSet 93 average Sharpe Ratio: 0.3981392177661411
Average Sharpe Ratio of all generated: 0.38782399114642674
SubSet 94 average Sharpe Ratio: 0.31298738000994497
Average Sharpe Ratio of all generated: 0.3870362373449901
SubSet 95 average Sharpe Ratio: 0.320383077107442
Average Sharpe Ratio of all generated: 0.3863419335925156
SubSet 96 average Sharpe Ratio: 0.4812804743445556
Average Sharpe Ratio of all generated: 0.3873206814353202
SubSet 97 average Sharpe Ratio: 0.38017584381363273
Average Sharpe Ratio of all generated: 0.38724777492897644
SubSet 98 average Sharpe Ratio: 0.33724981086237676
Average Sharpe Ratio of all generated: 0.38674274498890976
SubSet 99 average Sharpe Ratio: 0.5139184002313819
Average Sharpe Ratio of all generated: 0.3880145015413345

Search with Classical Algorithms

The subsequent step is to prepare the QUBO matrix and implement the functions required for the classical algorithm. Here, the pairwise correlation matrix is computed with logarithmic returns, as following the paper.

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

# compute Sharpe ratio from obtained portfolio
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

# introduce mutations to make next generation in genetic algorithm
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

# regenerate set of assets
# we use this set of assets from here
Nassets = 48

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

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

# convert from Bucket
SRBucket(SR_list)
CorrelationBucket(PairwiseCorrMat)

#print(PairwiseCorrMat)　# display pairwise correlation matrix

Greedy Method

First, we conduct a greedy search algorithm. Commencing from a randomly generated initial state, we select stocks that are likely to increase the Sharpe ratio of the portfolio, by examining their impact on overall performance. By iterating, we can observe that the algorithm converges from a random initial state to a particular state.

After execution, we compare the initial state with the state obtained by the greedy method. The reference of comparison is the Sharpe Ratio of the portfolio. In this tutorial, the calculation is based on the scenario where all stocks have the same weight, thus the average chart of the selected stocks is used to calculate the Sharpe ratio. We can also observe that if we reset the set of stocks by running the previous code block multiple times, the results of the greedy method may be inferior to the initial state or remain almost unchanged.

There also appear final states where the final returns are inferior to the initial state, but the Sharpe Ratio is larger. Comparing these final states, we can see that the portfolio with a larger Sharpe ratio has a smoother chart and smaller fluctuations (less risky).

#Greedy Search

# generate initial state
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)

# initialize energy based on initial state
Energies = list()
for iasset in range(Nassets):
    h = hi(SR_list, PairwiseCorrMat, iasset)
    energyTuple = [-1*abs(h), h , iasset]
    Energies.append(energyTuple)

# execute greedy search
NGreedyLoop = 10
for i in range(NGreedyLoop):
    heapq.heapify(Energies)

    ntry = 0
    #print(Energies) # display energy change
    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


# output of greedy search 
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.025433728324847

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] 0.8195667471626309

Genetic Algorithm

An alternate classical algorithm to the greedy method is the genetic algorithm. It selects the most optimal solution from a number of randomly generated portfolios and generates the next iteration of portfolios based on that optimal solution. The subsequent step of portfolios is generated through mutations that introduce random modifications to the previous step portfolios, and through crossover which exchanges genetic traits between portfolios. Given the relatively small number of stocks in this scenario, we only consider mutations and abstain from genetic crossover. By comparing the results, it is evident that the performance of the final computed generation surpasses that of the initial one. As the number of iterations attempted here is limited, performance may not improve or the improvement may be trivial. However, it is also apparent that, in general, the output is more stable than that of the greedy method.

Note: Increasing the number of genetic pairs, the number of offspring, and the number of mutations may result in an increase in computational complexity, causing the calculation to not be completed. Therefore, please adjust the settings according to your execution environment.

#Genetic Algorithm
Ngenes = 10 # the number of pairs of genes
NGenerations = 5 # the number of generations
# initialize randomly
Genes = list()
for i in range(Ngenes):
    gene = GenerateRandomSolution(Nassets)
    geneSR = EvaluateSolution(gene, Chart)
    Genes.append([geneSR, gene])
Genes = sorted(Genes,reverse=True)

# get the best portfolio in first generation
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)


# execute genetic algorithm
K_best = 5 # the number of top genes which are extracted
Ndescendants = 2 # the number of offspring
MaxMutation = 2 # the maximum number of mutations
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])

# get the best portfolio in the final generation
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)

# make a comparison plot
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] 1.2756964190448523

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] 1.6587225357221986

Implementation of quantum annealing method using OpenJij

Method for quantum annealing

To solve the optimization problem in reverse quantum annealing, we will translate the classical algorithm’s spin-form results into the formulation used in 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) # use the result from greedy search as an initial state
#QA_init_state = ConvertSolutionToQuboState(selected_Genes[0][1]) # use the best portfolio in the final generation of genetic algorithm as an initial state
#QA_init_state = GenerateRandomSolution(Nassets) # use random initial state
print(QA_init_state)

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

For forward annealing

Initially, we will execute Forward Annealing. We construct a QUBO matrix incorporating the stocks’ attractiveness and the pairwise correlations of penalties and rewards. In this following, we perform our computations through Simulated Quantum Annealing (SQA), which emulates quantum annealing on a classical computer. The enhancement in the final rate of return varies based on the case, yet we can observe that the charts demonstrate greater stability than the results procured by the greedy method.

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'>
[([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1], -261., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)]
[0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 0 1 0 0 1 0
 0 0 0 1 0 0 1 0 0 1 1] 1.723239194753436

For reverse quantum annealing

In reverse annealing, we input the state generated from the classical algorithm as the initial state before annealing. The first step is to set up the annealing schedule. In this tutorial, using OpenJij, the SQA Hamiltonian is defined as follows:

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

where $\mathcal{H}_{p}$ is the problem Hamiltonian that we seek to obtain the optimal solution for, and $s$ and $1-s$ correspond to $B[t]$ and $A[t]$ in reference paper [2] respectively. The value of $s$ can be set by adding supplementary arguments to sample_qubo.

sample_qubo(QUBO, num_reads=10)

We can modify the previous section’s sample_qubo to input the initial state settings and updates as follows:

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

where the schedule refers to the annealing schedule. It has the following structure: [0] component represents the value of $s$, reflecting the transverse magnetic field’s intensity during annealing; [1] component is the inverse temperature $\beta$, with a larger value denoting a lower temperature, facilitating the attainment of the ground state (default value of $\beta = 5$); and [2] component denotes the number of quantum Monte Carlo simulation steps, enabling us to alter the length of each schedule.

The initial state should be assigned to the initial_state argument. If the reinitialize_state option is set to False, the previous step’s annealing output is designated as the next step’s annealing initial state in the iteration. In the sample script above, commencing from the classical algorithm’s prepared initial state, RQA would execute ten times. The default of the optional argument, True, reset the initial state every time it runs; if we want to execute RQA once and contrast the results multiple times, we set reinitialize_state = True.

Now, let’s fabricate a schedule for reverse annealing. Since OpenJij sampler performs quantum Monte Carlo with a given $s$ and $\beta$ over a specified number of MC steps, we need to create a schedule in which $s$ and $\beta$ gradually change by linking short constants to facilitate RQA’s functionality. Here, we focus on the change in $s$ and plot it to observe an inverted trapezoidal schedule. The horizontal axis represents the quantum Monte Carlo method’s number of Monte Carlo steps. By rewriting the schedule creation function’s contents, we can generate any schedule shapes.

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)

Let’s execute the RQA utilizing the sample_qubo function in conjunction with the schedule and initial conditions set up thus far. Simultaneously, we compare the outcomes of the greedy method, forward annealing, and RQA. The results indicate that the setup in this tutorial frequently yields the same optimal solution as forward annealing, owing to the constraints imposed by the size of the set of assets.

# prepare initial state 
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)

# compare charts
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()
# compare the results
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))

[([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)]
Forward Quantum Annealing Result:
[0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 0 1 0 0 1 0
 0 0 0 1 0 0 1 0 0 1 1] 1.723239194753436
Reverse Quantum Annealing Result:
[0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 0 1 0 0 1 0
 0 0 0 1 0 0 1 0 0 1 1] 1.723239194753436

Checking reverse quantum annealing

Since the final output alone fails to provide insight into the behavior of RQA, we dissect the phase and observe its behavior. Initially, we investigate the reverse phase. Similarly, an annealing schedule is set up and annealing is carried out using the designated schedule and initial state. Here, for clarity, we have altered the $s$ value and inverse temperature at which the reverse phase terminates. Additionally, annealing is arranged to commence from the same initial state on each occasion, as opposed to undergoing iteration. The results show that the reverse phase process ascends from the previously attained solution to a higher energy state. Furthermore, as the solution stays in a different state during each execution, we can deduce that the solution differs at a higher energy level.

#Create RQA schedule
RQAschedule = []
NReverseStep = 50
TargetS = 0.15 # return to more random state よりランダムな状態に戻す
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)

[([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -262., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1], -253., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -262., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1], -261., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -262., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1], -261., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -261., 1)]

Next, we add a pause phase. After the addition of the pause phase, annealing is performed similarly, and the results are evaluated. As the addition of this phase is tantamount to extending the number of MC steps of the final stage of the reverse phase described above in terms of implementation, the results are as dispersed as those obtained in the reverse phase-only scenario depicted previously.

#Create Pause Step
NPauseStep = 50
TargetS = 0.15 # return to more random state 
beta = 1.0 # allow for easier thermal hopping at higher temperatures
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) # execute annealing from the same initial state every time 
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)

[([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1], -237., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1], -240., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1], -194., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1], -240., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1], -227., 1)
 ([0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1], -192., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -249., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1], -211., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1], -227., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -258., 1)]

Finally, we incorporate a forward phase and execute annealing. This configuration is set up under conditions of slightly relaxed temperature and magnetic field relative to the previous sample. As a result, we can observe that the known optimal solution through quantum annealing is obtained for this particular set of stocks. Similar to the first example in the RQA section, we can observe that the optimal solution is obtained by RQA.

Since the conditions are slightly lenient, the solutions do not converge, and we can see the presence of multiple sub-optimal solutions with comparable energies. Among these solutions, we can also detect the presence of a portfolio whose final return outperforms the optimal solution obtained through forward annealing in certain cases and whose profile is relatively smooth. A closer examination of that portfolio details reveals that the substantial returns were achieved through a combination of specific stocks exhibiting substantial appreciation rates and stocks with relatively low volatility that possess inverse correlations. Although such a portfolio may be desirable in actual investment scenarios, the “Shape ratio is maximized” condition set up in this tutorial resulted in a relatively modest portfolio Sharpe ratio, which was deemed suboptimal through quantum annealing.

#Create Forward Step
NForwardStep = 50
TargetS = 0.15 # return to more random state 
ForwardStep = (1.0 - TargetS) / NForwardStep
beta = 5.0 # allow for easier thermal hopping at higher temperatures
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) # execute annealing from the same initial state every time 
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)

[([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -262., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1], -263., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1], -261., 1)
 ([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1], -263., 1)]

Parameter search for specifying RQA schedule

Examining the effects of the parameters in separate phases reveals that results vary based on the values of $s$ and $\beta$, as well as the number of steps in the Quantum Monte Carlo method. From this perception, to achieve optimal performance with reverse annealing, it is necessary to determine appropriate values for these parameters.

In reverse phase

Initially, we consider their impact on the reverse phase. We start by analyzing the effect of $s$, which represents the degree of quantum annealing. The transverse magnetic field intensity represented by $s$ must be set to an optimal value to allow for dispersion beyond the potential barrier around the local solution reached through initial annealing. Results show that for $s \leq 0.20$, thermal fluctuations become dominant and variations increase, so a finer scan around $s \sim 0.2$ is expected to yield the optimal value.

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) # execute annealing from the same initial state every time 
    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()

Next, we consider the effect of $\beta$. $\beta = \frac{1}{k_\mathrm{B} T}$ represents the inverse temperature, where a larger $\beta$ means that the system is cooler. In other words, the larger $\beta$ is, the more easily the overall state reaches to the ground state and behaves more likely an “annealed” state. For a constant transverse magnetic field, as $\beta$ increases, the effectiveness of the transverse magnetic field decreases and begins to resemble that of a large $s$. For more details on this effect in OpenJij, please refer to this Qiita article.

It is acknowledged that the effect change by the typical scale of energy present within the QUBO matrix. In this tutorial, the typical value of the QUBO matrix roughly equates to $5$, thus we expect the variation to occur around $\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')

We generate a logarithmic array of $\beta$ and implement a Reverse Phase with all parameters fixed, excluding $s$. As expected, the dispersion of the final state becomes minimal around $\beta=5$, and for higher values of $\beta$, we observe that the states converge to a few states.

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) # execute annealing from the same initial state every time 
    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()

Next, we analyze the impact of schedule partitioning. Here we examine extreme cases, such as the step function, as well as the scenario of a gradual alteration achieved by dividing the schedule into smaller segments. As the results indicate, the state does not significantly change when the number of divisions is limited, whereas an increase in the number of divisions leads to noticeable variation. Additionally, a substantial increase in the number of divisions only increases the computational time and trivial differences in the outcomes. For the portfolio optimization problem addressed in this tutorial, approximately $50$ is found to not pose a problem.

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) # execute annealing from the same initial state every time 
    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()

Also, with regards to MC steps, it is determined that its impact on the final result is negligible unless the product of (number of divisions)x(MC steps) is excessively small. In the case of the default setting with SQASampler(), about $1000$ should work.

For pause phase

As already explained, the pause phase is an extension of the final step of the reverse phase, behaves similarly to the reverse phase when altering $s$ and $\beta$ of the overall annealing. Also, since the steps of the pause phase is equivalent to the extension of the final step, its impact on the reverse phase is expected to be limited, even if it is changed.

However, modifying only $\beta$ during the pause phase resembles classical annealing with incomplete quantum annealing. Here we examine the effect of changing the system temperature from $\beta=5.0$ to a gradual target $\beta$. We obtain the results that are similar to the final state when annealing is performed at a different $\beta$ in reverse phase.

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) # execute annealing from the same initial state every time 
    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()

For forward phase

The forward phase is the stage in which standard quantum annealing is executed, starting from the states acquired through the reverse and pause Phases. The conditions for the parameters in this phase we have to consider are essentially equivalent to those required for standard annealing to work correctly. We examine $\beta$ and the number of division steps in the same way as reverse and pause phases , we find that even with a modest number of divisions and $\beta$, the results converged to the optimal solution.

An example script is shown below. In the reverse phase and pause phase, we employ the settings obtained from previous investigations, with only the forward phase undergoing annealing with a coarse resolution and small $\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) # execute annealing from the same initial state every time 
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()

Concluding remarks

We implemented reverse quantum annealing utilizing OpenJij to solve a portfolio optimization problem. In this tutorial, we elucidated the RQA problem formulation, the procedure for specifying the schedule, and the parameters to take into account when performing RQA. We showed that the RQA method can be used to obtain the correct optimal solution.

Similar methods can be employed to discover a better solution for other difficult problem with numerous suboptimal solutions close to the optimal one. If we only obtain a predicament with local minimum through standard quantum annealing, attempting RQA in this manner may assist in resolving the predicament.

Additionally, there were some aspects of RQA that we did not discuss in this tutorial, such as the effect of the schedule from a linear function to a nonlinear form, yet we can be easily implemented it by modifying the code in this treatise. If you are interested, please give it a try.

Appendix A: if the simulation results are not good…

The script in this tutorial to generate a stock set is not guaranteed to yield an “admirable” set of stocks. Therefore, there exists the possibility of experimenting with sets that possess a biased distribution in certain instances. In such situations, the genetic algorithm calculation might experience an extended execution time or RQA might not converge with parameters such as $s$ or $\beta$ that were utilized by default in this tutorial. In such cases, we might consider regenerating the stock set or manually adjusting the RQA parameters to reach the optimal resolution.

Appendix B: Executing RQA on the D-Wave machine

RQA can also be performed on D-Wave machines by using the method described in this tutorial. We show a sample script below. The schedule settings and other arguments differ based on the implementation and the discrepancy between the simulation and the actual machine, however, the implementation is fundamentally the same. For more details, please refer to the documentation on the official D-Wave website. Moreover, the difference in outcomes between standard annealing and RQA is more likely to become apparent on an actual device, given that the effects of noise and the like are more likely to be observed.

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) # start and end time of each phase 
Ratio = (1.0,0.38,0.38,1.0) # value of s when each phase changes
schedule = list(zip(timing,Ratio)) # put them together in a schedule 
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

References

1. Harry Markowitz, “Portfolio selection”, The journal of finance, 7(1):77–91 (1952)

2. Davide Venturelli, Alexei Kondratyev, “Reverse Quantum Annealing Approach to Portfolio Optimization Problems”, Quantum Machine Intelligence volume 1, pages17–30 (2019)

3. Sharpe, William F., “Mutual fund performance”, The Journal of Business 39 (1), 119-138 (1966)