Solving the Optimal Trading Trajectory Problem Using Simulated Bifurcation
SSolving the Optimal Trading Trajectory ProblemUsing Simulated Bifurcation
Kyle Steinhauer
Scientist at AlpacaJapanAlpacaJapan Co., Ltd.Email: [email protected]
Takahisa Fukadai
Scientist at AlpacaJapanAlpacaJapan Co., Ltd.Email: [email protected]
Sho Yoshida
Chief Data Scientist at AlpacaJapanAlpacaJapan Co., Ltd.Email: [email protected]
Abstract
We use an optimization procedure based on simulated bi-furcation (SB) to solve the integer portfolio and trading tra-jectory problem with an unprecedented computational speed.The underlying algorithm is based on a classical descriptionof quantum adiabatic evolutions of a network of non-linearlyinteracting oscillators. This formulation has already provento beat state of the art computation times for other NP-hardproblems and is expected to show similar performance forcertain portfolio optimization problems. Inspired by such weapply the SB approach to the portfolio integer optimizationproblem with quantity constraints and trading activities. Weshow first numerical results for portfolios of up to 1000 as-sets, which already confirm the power of the SB algorithmfor its novel use-case as a portfolio and trading trajectoryoptimizer.
The trading trajectory problem can be described as theproblem to find the optimal set of portfolios and respectivetrading activities, which maximize the future expected re-turn over a given time period while taking into account trad-ing costs and expected risks. The optimal trajectory doesn’tnecessarily maximize the return at each time point, but max-imizes the future return over the entire time-period includingall trading activities. We will follow the mean variance port-folio description [1, 2, 3] in order to express portfolio valuesat each point in time and will add cost terms that account forthe necessary rebalancing from one time step to the other.There are many scenarios in asset management wheretrading and investment activities are constraint. Some com-mon constraints lead to the problem of finding the optimaltrading trajectory with only integer-valued solutions. Thisoccurs for example in ETF block trades, where only a certain integer amount of a standard lot-size can be traded. The inte-ger portfolio formulation (in its general form) belongs to thenon-convex mixed integer quadratic problems and thereforefalls into the class of NP-hard problems [4, 5]. This stronglycomplicates the search of the optimal trajectory and partic-ularly effects the computation time with increasing systemsize. A further complication comes from quantity and cardi-nality constraints that are required in almost any real worldapplication. A variety of studies [6,7,8,9,10] have been con-ducted to push the understanding and computational perfor-mance of the integer portfolio optimization with constraints.A branch and bound method [11] for example, allowed tofind the exact solution to a single time period portfolio opti-mization problem with up to 200 assets.In the last couple of years we saw substantial progressin building quantum annealer systems and actual quantumcomputers, and the interest is rising to harness this technol-ogy for real world applications like the discussed mixed inte-ger quadratic problems. An important step into this directionhas been conducted by G. Rosenberg et al. [12], which solvedthe integer trading trajectory problem on D-Wave’s quantumannealer system. Very recent studies [13] have shown con-tinued interest and progress in this direction. Despite theincreasing accessibility and power of those machines, the ac-tual business applicability in finance is however still missingcertain technological advancements.H. Goto et al. however, have been studying quantumadiabatic evolution in detail, and inspired by such presenteda recent formulation [14] that can potentially bring quan-tum computer like speed for certain descriptions of NP-hardproblems to classical computers. In their recent study [14]they introduced a formulation of the Ising model using Kerr-nonlinear parametric oscillators, which is highly effective insolving the Ising optimization problem via its quantum adi-abatic evolution through its bifurcation point. This simu- a r X i v : . [ q -f i n . C P ] S e p ated bifurcation (SB) algorithm is highly suitable for parallelcomputing and beats (around 10 × faster) the current state ofthe art custom built machine to solve a fully connected 2000spin problem [15]. Since many NP-hard problems can beexpressed in form of an Ising problem [16], the simulatedbifurcation has the potential to boost computational perfor-mance for a large range of those combinatorially difficultproblems. The discovery of this highly promising formu-lation has also triggered larger institutions, like Toshiba forexample, to construct dedicated simulated bifurcation ma-chines (SBM) and already advertise potential real world use-cases, e.g. the detection of triangular arbitrage.In this work we make use of this novel SB-formulationto solve the optimal portfolio and optimal trading trajec-tory problem. For this we will first map the integer portfo-lio problem to the respective Ising problem, recapitulate thementioned SB-algorithm, and then use our implementationto generate first results of optimal portfolios and trajectoriesfound by simulated bifurcation. The Ising problem [17] can be formulated as the prob-lem to find a ground-state (spin configuration) that minimizesthe Ising energy defined by E = − ∑ i , j J i j s i s j + ∑ i h i s i , (1)where s i ∈ {− , } are commonly interpreted as spins point-ing down or upwards. Note that both J i j and h i are ∈ R and the coupling among the spins is symmetric, such that J i j = J ji . The sum runs over all N spin-variables in the sys-tem. This problem has been studied in connection to vari-ous fields of physics, mainly due to the spontaneous sym-metry breaking that occurs in 2 or more dimensions. Thisis something observed in a real world ferromagnet but alsothe mechanism behind the mass generation in the StandardModel of particle physics [18, 19, 20]. For this work how-ever, the more important point of interest is the combinatorialand computational aspect of the problem. Without any topo-logical structure, the fully connected Ising problem falls intothe class of NP-hard problems [21]. The number of possiblestates is 2 N , where N is the number of spins in the system.For decently sized systems it is impossible to deterministi-cally find the ground-state in a reasonable amount of time ona classical computer by iterating through all combinations.Monte-Carlo algorithm like the Metropolis algorithms andmany others, are successfully used to sample the differentspin states and to approach the ground-state without gener-ating all combinatorial possibilities. Note however, that forfully connected Ising systems (non-local actions), most ofthese algorithms show rather weak performance, are not suit-able for parallel computing, or simply not usable.The already mentioned work [14] by H. Goto et al. de-scribes the Ising system with Kerr-nonlinear parametric os-cillators and proposes a new optimization algorithm that sim- ulates adiabatic evolutions of classical nonlinear Hamilto-nian systems. If the system is initialized in a proper vacuum-state, the so called simulated bifurcation (SB) algorithm sim-ulates the adiabatic evolution of the system through its bifur-cation point and ends up in a ground-state that minimizesthe Ising energy. Using this formulation has shown to beatstate of the art computation times for large fully connectedsystems. Note that for smaller systems, different heuristicapproaches, like the so called Digital Annealer [22], mightstill outperform the SB-approach, however do not allow forseparate updates of the system’s variables and are thereforeless suited for parallel computing.The SB-algorithm is, to our knowledge, currently thefastest way to solve a fully connected Ising problem andtherefore also an ideal candidate to solve the optimal inte-ger portfolio problem in the Ising representation. In this section the optimal trading trajectory problem isformulated and mapped to a binary-bit representation andthen mapped to the Ising representation.
The optimal trading trajectory problem can be describedas a temporal sequence of mean variance portfolios [1, 2, 3],which are composed of an expected return, an expected riskand an expected trading/rebalancing cost term. Mathemati-cally this can be expressed by w = arg max w ∑ t (cid:16) w Tt µ t − γ w Tt Σ t w t − ∆ w Tt Λ t ∆ w t (cid:17) , (2)where w t is the portfolio weight vector of size N at time t and the individual entries fulfil w it ∈ N . The vector µ t is theestimated future return vector for time t , Σ t the estimated co-variance and γ a risk aversion parameter. The matrix Λ t holdsthe costs to perform the rebalancing of investments from timepoint t to t + ∑ Ni = w it ≤ ˜ K at all time-points, where˜ K ≥ K is constant over time. We will intro-duce a further constraint that allows only for limited invest-ments in each asset. This constraint could also be time de-pendent and different for each asset, for simplicity we how-ever choose a time independent and uniform constraint of theform ∀ t : 0 ≤ w it ≤ ˜ K / N , that allows only a certain fractionof our total units (total capital) to be invested in one asset. Both, the Ising problem and the portfolio problem, atleast in their native form, have the same structure and ask usto minimize an expression with a linear and a quadratic term.The formulation with spins, where s i ∈ {− , } , is closely re-lated to the perhaps more intuitive and well known binary bitepresentation where b i ∈ { , } , which can be obtained by b i = ( s i + ) /
2. Therefore we first map the portfolio prob-lem with integer weights to the binary-bit representation andthen to the binary-spin representation. Mapping the optimalporfolio problem to a binary representation has already beenconducted multiple times using various binary representa-tions [12].Let us first move to the binary representation by rewrit-ing w it = ∑ log ( K )+ k = k b ikt where b ikt ∈ { , } in which K isthe largest integer we want to represent in w (note the dif-ference to ˜ K ). All individual bits b ikt can be packed into asingle vector b , with which eq. 2, dropping Λ t for now, canbe reformulated as b = arg max b ∑ t (cid:16) b Tt ˆ µ t − γ b Tt ˆ Σ t b t (cid:17) , (3)where the hat simply indicates that these are the respectivequantities in the new basis. The future return vector and thecovariance matrix will need to be extended and appropriateentries multiplied with 2 d , the exact order and value dependson how the elements of b are ordered. Further the lengthof b is larger than w by a factor of (cid:98) log ( K ) (cid:99) +
1. In orderto eliminate all ambiguity let us for clarity write out all thesums that are hidden in the vector multiplication above. b = arg max b ∑ t (cid:32) ∑ i ∑ k k b ikt µ it − γ ∑ i , j ∑ k , l k l b ikt b ilt Σ i jt (cid:33) , (4)where the sum over k and l starts at 0 and goes until (cid:98) log ( K ) (cid:99) + i and j are the indices labelling all as-sets. The term Σ i jt stands for the covariance between asset i and j at time t . We chose to order b according to the bitendianness, meaning bundling by significance of respectivebit. In this basis the transformed future returns ˆ µ and ˆ Σ canbe easily constructed by a trivial block-wise extension of theoriginal quantities found in eq. 2. Each block (and respectiveentries) ˆ Σ klt are simply obtained by multiplying the Σ t matrixwith the corresponding bit significance, i.e. 2 k l Σ t , which isan N × N matrix block. For the linear term ˆ µ t we simply ex-pand µ t and multiply the respective entries with 2 k . We havepicked one binary representation and also one distinct choiceof basis, note however that there are many other options withtheir advantages and disadvantages. An interesting overviewon different binary representations was part of the study [12]around the portfolio optimization on a D-Wave system.The portfolio problem in its binary-bit basis in eq. 3 cannow be transform to the spin representation by introducing b = ( s + ) /
2, where denotes the unit vector with samelength as s . This introduces scalar correction terms inde-pendent of s , for the optimization step however they can bedropped and the overall expression can be compressed to s = arg min s ∑ t (cid:104) γ s Tt ˆ Σ t s t + (cid:16) γ ˆ Σ t − ˆ µ (cid:17) s t (cid:105) , (5) where we also flipped the sign in order to minimize the ex-pression. From here we can directly read off the correspond-ing couplings of the respective problem in the Ising represen-tation. The quadratic blocks (in our basis) appearing in theblock diagonal matrix J are J t = − γ ˆ Σ t . The linear term in theIsing representation can also be broken down into temporalsequences that are obtained via h t = γ ˆ Σ t − ˆ µ t . Note that theblock diagonal nature of J is extended with further blockslinking spins of different time points with each other, as soonas the transaction costs in Λ t are introduced. These blocksin the Ising representation ˆ Λ t , linking spins from differenttime-points can be constructed such that they introduce a ten-dency for identical spins on neighbouring time points to bealigned, unless the contribution of the expected future returnis strong enough to trigger a spin flip. A way how this can beachieved is to add terms of the form s ik , t s ik , t + c i ( t , t + ) k to the sum, where i denotes the asset, k the bit-significanceof the respective spin and c i ( t , t + ) the transaction penalty.We will follow exactly this approach, will however choosea constant cost c i ( t , t + ) = c for all assets and time points.Like this the diagonal matrices that appear in off-diagonalblocks of the interaction matrix J, reduce to ˆ Λ t = c , wherewe omitted writing the multiplication of the spin dependentfactors of 2. The just discussed approach is one way to gen-erate energy-gaps in units of c depending only on ∆ w .To include constraints in the Ising representation can bemore elaborate than in the integer representation of the prob-lem. Various different constraints can for example demandthe introduction of ancillary spins [23], this is why we chosea set of constraints that are directly embedded as hard con-straints in the formulation of the problem. If for a given asset i we choose to have only (cid:98) log ( ˜ K / N ) (cid:99) + w it ≤ ˜ K / N and hencealso ∑ Ni = w it ≤ ˜ K . In this section we discuss the core components of thesimulated bifurcation algorithm, introduced in a recent study[14], from which we will borrow heavily in this section (con-tent and also nomenclature) .
Following the exact steps from [14] we formulate theIsing energy defined in eq. 1 with a network of Kerr non-linear parametric oscillators. In a quantum mechanical for-mulation the Hamiltonian is given byˆ H ( t ) = ¯ h ∑ i (cid:18) K a † i a i − p ( t ) ( a †2 i + a i ) + ∆ i a † i a i (cid:19) − ¯ h ξ ∑ i , j J i j a † i a j + ¯ h ξ A ( t ) ∑ i h i ( a † i + a i ) , (6)where a † i is the creation and a i the annihilation operator forthe i-th oscillator. The parameter ∆ i is the detuning frequencyhich plays an important role when defining the initial vac-uum state. The time dependent parameter p ( t ) is the pump-ing amplitude, ξ is a constant parameter (in units of fre-quency) and K is the Kerr coefficient. Due to the same argu-ments given in [14] we assume that K , ∆ and ξ are positive.For our considerations the numerical value of the reducedPlanck constant ¯ h is irrelevant, note however that it also car-ries the unit of time which gets cancelled with the unit offrequency (inverse time) in ξ , guaranteeing that the entireexpression is in units of energy only. A ( t ) is a positive di-mensionless parameter that increases with p over time, suchthat A ≈ p (cid:28) ∆ and A ≈ (cid:112) ( p − ∆ ) / K when p (cid:29) ∆ .A quantum adiabatic evolution of this Hamiltonian is de-sired which will end up in the ground-state that will mini-mize the Ising energy. In order to achieve this we initializeall oscillators in their vacuum states and gradually increasethe pumping amplitude p ( t ) from zero to a sufficiently largevalue compared to ∆ and ξ . In order to initialize the systemin the vacuum state a proper tuning of ∆ i can be necessary.The reader is referred to H. Goto’s work [24] for more detailson the initialization of the vacuum-state. Also finding a proofin the appendix that shows that if the variation of p ( t ) is suf-ficiently slow, the final state will become the ground-state ofthe final Hamiltonian by the quantum adiabatic theorem.In the following we will formulate the correspondingclassical Hamiltonian which can be derived by approximat-ing the expectation values of a i via a complex amplitude x i + iy i (note that introductory literature often uses p and q here instead). The real and imaginary part form a conjugatevariable pair that correspond to position and momentum ofthe i-th oscillator. With x i = ( a i + a † i ) / y i = ( a i − a † i ) / H ( x , y , t ) = ∑ i (cid:20) K ( x i + y i ) − p ( t ) ( x i − y i ) + ∆ i ( x i + y i ) (cid:21) − ξ ∑ i j J i j ( x i x j + y i y j ) + ξ A ( t ) ∑ i h i x i (7)Note that in this formulation, after adiabatic evolution, thespin value s i corresponds to the sign of the amplitude x i .Since eq. 7 describes a classical Hamiltonian we can derivethe equations of motions for variables by taking the deriva-tive with respect to time, denoted with a dot, these are givenby following the classical time evolution formulas of Hamil-tonian mechanics:˙ x i = ∂ H ∂ y i = (cid:2) K ( x i + y i ) + p ( t ) + ∆ i (cid:3) y i − ξ ∑ j J i j y j (8)˙ y i = − ∂ H ∂ x i = − (cid:2) K ( x i + y i ) − p ( t ) + ∆ i (cid:3) x i + ξ ∑ j J i j x j − ξ A ( t ) h i (9)Note that there is a minus sign in front of the derivative of the second conjugate variable. These formulas describe kineticsthat allow us to simulate classically the quantum adiabaticevolution. The equations of motion, derived in eq. 8 and 9, can befurther simplified in order to be more suitable for fast numer-ical simulation. Again following [14] the terms proportionalto the momenta y , which vary around zero, can be droppedand the equations of motion can be reformulated as˙ x i = ∆ i y i (10)˙ y i = − (cid:2) Kx i − p ( t ) + ∆ i (cid:3) x i + ξ ∑ j J i j x j − ξ A ( t ) h i (11)Note that this approximation allows us to use the symplecticEuler method [25] to simulate the hamiltonian dynamics ofthe system because the two variables are now separable. Wetherefore discretize time with t = n ∆ t , where ∆ t is our timeincrement. Like this we can write the algorithmic update stepfor the position variable as x i ( t n + ) = x i ( t n ) + ∆ i y i ( t n ) ∆ t (12)and the update of the momentum variable as y i ( t n + ) = y i ( t n ) − (cid:2) Kx i ( t n + ) + ( ∆ i − p ( t n + )) x i ( t n + ) − ξ ∑ j J i j x j ( t n + ) + ξ A ( t n + ) h i (cid:35) ∆ t . (13)From here a number of further different simplifications andre-formulations can be made in order to arrive at an evenfaster algorithm. The expression above however describesthe core update steps of the SB algorithm. The results section is divided into four smaller sub-sections. First we will discuss results obtained for portfoliosat a single time point, then present solutions for trading tra-jectory problems, then discuss performance and finally con-sider problems where only close-to-optimal solutions werefound.
In this section we consider results, obtained by the SB-algorithm, for the optimal portfolio problem formulated inthe Ising representation. For the results in this section wecreate an artificial market situation with N different assets bysampling random returns from a geometric Brownian motionover 1000 time increments in order to estimate a future return ig. 1. Shows random portfolios (including the important edgecases) and those found via adiabatic evolution in the simulated bi-furcation algorithm for various different risk aversion parameters γ .This is for an artificial portfolio with N = assets, a total capital of ˜ K = units, where can be spent maximally per asset. Theaverage expected return of the sampled returns is 0.5% vector µ and a covariance matrix Σ . Optionally we might in-clude a drift in the market that results in an average expectedreturn of κ .In fig. 1 we show the configurations found by the SB-algorithm for a small portfolio of N = K =
75 units available and a maximal amountof ˜ K / N =
15 units to be distributed per asset. Further 10 random portfolios are added, with the same constraints, toillustrate the universe of different investment options. Thisis done in a market with a positive drift resulting in an aver-age expected return of κ = . γ from zero to a sufficiently large value.We repeat this experiment and replace a random assetwith a risk free asset with 1% future return. If we push therisk aversion parameter γ to large enough values we shouldend up in the scenario where we invest only in the risk freeasset and avoid any investments in risky assets. Exactly thisresult is shown in fig. 2, where we see that the SB-algorithmcorrectly finds the optimal portfolios (as also highlighted bythe optimal frontier line calculated via a quadratic solver),and ends up in the just described risk free scenario when γ reaches large enough values. Since ˜ K / N =
15 and the riskfree asset’s return is 1%, the risk free portfolio correspondsto the point (0, 0.15), as shown in the figure.In fig. 3 we are considering a scenario where the num-ber of assets is N =
400 and we can either buy one unit of theasset or not. This allows for very fast computation times ofless than half a second per point on a standard desktop CPU.For the extreme case of N = K / N =
1, the SB-algorithm uses roughly 1 seconds to find the optimal solu-tion, which even beats many out of the box quadratic solvers
Fig. 2. Shows random portfolios (including some edge cases)and those found via adiabatic evolution in the simulated bifurcationalgorithm for various different risk aversion parameters γ . This is foran artificial portfolio with N = assets, a total capital of ˜ K = units, where can be spent maximally per asset. The averagevalue of the sampled returns is 0%, but we introduced a risk freeasset with 1% return, visible by the off-set on the y-axis. Further theoptimal portfolio boundary is highlighted in green.Fig. 3. Portfolio with N = assets, a total capital of ˜ K = units, where can be spent maximally per asset, either purchase ornot. The average expected return of the sampled returns is 0.5%.Computation time of the approximation via the quadratic solver of thecontinous problem and the SB-approach were identical, both around0.3 seconds per point on a single Intel Core i5-7200. to tackle the respective continuous problem with the identicalconstraints. We investigate the performance in more detailin section 5.3 and in the next section consider the solutionsfound by the SB-approach for the optimal trading trajectoryproblem. .2 Optimal Trading Trajectory with SB-algorithm We now construct T future random market conditionswith different values of average expected returns and use theSB-algorithm to find the optimal asset allocation trajectoryover time. We first consider a setup with only N = K / N =
3, but T = differ-ent trajectories from which we want to find the optimal one.The trading costs are organised such that for a fixed asset weget a cost penalty of c = .
01 for each unit changed betweentwo time steps. We will consider the optimization for threedifferent γ values and consider the trajectories found by theSB-algorithm.In fig. 4 at the top we show the trajectory value w Tt µ t − γ w Tt Σ t w t − ∆ w Tt Λ t ∆ w t at each time point t . The redline again indicates the trajectory found by the SB-algorithm,whereas the blue samples indicate the universe of possibletrajectories. In the second row of fig. 4 only the return term w Tt µ t is shown and in the bottom only the risk component w Tt Σ t w t . For γ =
0, and also for sufficiently small values ofthe risk aversion parameters, the optimization procedure canignore the risk, which is exactly what is observed when look-ing at the solution found by the SB-algorithm. As we see inthe bottom of fig. 4, the trajectory is exposed to a lot of riskand solely focuses on optimizing the return and trading costs.
Fig. 4. Optimal trajectory for N = , ˜ K / N = and T = forzero and close to zero risk aversion. In fig. 5, we are looking at the same scenario as in fig.4, this time however, with an increased risk aversion param-eter set to γ = .
02. This change makes us more risk averseand we observe how the SB-algorithm finds a trajectory forwhich risk is only taken if the magnitude of expected returnis large enough, else, if the risk term dominates, portfoliosare found that have a sufficiently small risk value. In the topof fig. 5, we see that the overall portfolio value, includingreturn and properly weighted risk still looks meaningfullymaximized.In fig. 6 the risk aversion parameter was increased fur-ther in order to check if we observe the extreme scenario in
Fig. 5. Optimal trajectory for N = , ˜ K / N = and T = for γ = . . which the risk term is dominating completely and forcing azero investment trajectory. This is exactly what we observewhen looking at the solution found by the SB-algorithm, riskis minimized completely by not suggesting any investmentsduring the entire time period. In the middle panel of fig. 6 weagain see the potential future return, from which we howevernot benefit due to the increased risk awareness. Fig. 6. Optimal trajectory for N = , ˜ K / N = and T = suf-ficiently large risk aversion parameters. The SB-algorithm finds trajectories under different con-ditions that are meaningful and that seem to coincide withthe expected results, we have however not delivered proofthat, at finite trading costs, the found solutions correspond tothe true optimal trajectory. Let us in the following establishnumerical evidence that the SB-algorithm finds the trajectoryof configurations that correspond to the global optimum. Forthis we consider a smaller system with only 3 time points. Inthis small system it is possible to compare all possible trajec- ig. 7. Shows sorted values of all possible trajectories in asetup where N = , ˜ K / N = and T = . The In addition the con-figuration obtained with the SB-algorithm is plotted, giving numericalevidence of finding the global optimum for the mentioned system.Note that the corresponding trading costs of each trajectory havebeen included in the total portfolio value. tories with the solution obtained by the SB-algorithm. In fig.7 we show the ordered accumulated portfolio values of all2 possible trajectories for a system with N =
3, ˜ K / N = T = c = .
01. We observe that the SB-algorithm findsthe true global optimum out of all configurations.If we increase the temporal dimension to T = possible combinations, where our compar-ison with a brute force method is infeasible and the com-parison to random trajectories becomes meaningless. Forthose larger systems we were only able to check the cor-rect behaviour of the approach in certain limits. In the limitof vanishing trading costs for example, i.e. c →
0, the sys-tem reduces to a time independent optimization problem of T separately solvable portfolios. The numerical checks inthat limit show that the global optimization coincides withthe trajectory that optimizes the portfolios at each point intime separately. This is illustrated in fig. 8, where we seehow the globally optimized portfolio converges to the localoptimization when decreasing the trading penalty. Fig. 8. The blue line shows the total trajectory value of the solutionsof the global optimization problem for different finite trading costs.The red line indicates the value obtained when optimizing every port-folio at each point in time separately, i.e. performing local optimiza-tion. This is for a system with N = , ˜ K / N = and T = . Numerical checks are only possible in certain limits,and hence we can not guarantee the correctness of solutionsfound for much larger random markets with non-zero tradingcosts. The confirmation obtained by comparing smaller sys- tems with brute force methods and the correct behaviour oflarger systems in given limits, give us confidence, howevernot proof, that we find optimal or close-to-optimal resultsfor larger systems with finite trading costs. Another exam-ple where a direct cross check with another method was notpossible is illustrated in fig. 9, where we optimize a portfoliowith N =
10, ˜ K / N =
15 and T = possible combinations. Fig. 9. Optimal trajectory for N = , ˜ K / N = and T = . Note that for some scenarios we observed certain excep-tions where the SB-algorithm only found close-to-optimalsolutions, those will be discussed in section 5.4.
Fig. 10. Trading trajectory investment suggestion from the SB-algorithm for N = , ˜ K / N = and T = with a cross-assetseasonal and a random component in the expected future market.On top of the random component a seasonal affect is included. Inthe top panel the trading cost term is set to c = . allowing forfrequent trading, whereas in the bottom to c = . , making anychange extremely costly. In order to gain insight into the optimal trading strategyof a selected trajectory we can display the number of sug-gested units to be held for each asset for each point in time. Asimple illustration of this with a portfolio of N =
4, ˜ K / N = c and, as expected, observe that the optimal trajectoryfound by the SB-algorithm avoids changing the positions asfrequently as before. The future market was generated witha seasonality and a random effect in order to enforce alsoseasonal effects in the trading pattern. The result in fig. 10establishes confidence that the penalty term introduced in theIsing representation allows to control the trading costs, e.g.the amount of trades. This however is only a heuristic san-ity check and numerical evidence which proves that the SB-algorithm found the global optimum including trading costswere only possible in the already discussed setups. Despitethe successfully controllable trading activities, the transla-tion from real world trading costs and concepts to properlycalibrated penalty terms in the Ising formulation still formsan open challenge that is subject of ongoing investigations. Goto et. al. [14] took advantage of the fact that the SB-algorithm can be set up in a highly parallel manner on a GPUcluster and thereby solved a fully connected 2000 spin prob-lem 10 × faster than the current state of the art custom builtlaser machine. We are following the same path for the port-folio optimization by porting the discussed SB optimizationprocedure to a GPU cluster. Performance numbers basedon the GPU implementation are currently not available, wewill however consider the already convincing performanceon a single Intel Core i5-7200 with 2.5GHz (fully vectorizedeq. 12 and eq. 13 in Python 3.7). For this we measure thetime of the adiabatic evolution in relation to the system size.For comparison we keep the parameters of the SB-algorithm( p max , ε p , ∆ t , ∆ i ) fixed, when increasing the number of assets N and the number of maximal units per asset ˜ K / N . In fig.11 we see the level of increase in computation time whenincreasing N . Despite the fact that the computation timeshows an exponential increase with increasing system size,it is astonishing to observe that the simulated bifurcation ofa system with 256 assets, with maximal investment of 512units per asset is performed under 4 seconds. It cannot bedirectly compared with other methods that ran in differentsetups (e.g. different constraints), it is however noteworthyto mention that it clearly beats the currently existing com-petitive numbers, like the 200 asset optimization by a branchand bound method that showed an average run-time of 4800seconds [11]. The move to a GPU cluster will not only al-low to drastically decrease the already fast computation timebut will also allow to include an extremely large amount ofassets.The measurements displayed in fig. 11 create a valuableinsight into the proper scaling of the algorithm, it however isnot guaranteed that a ground-state is found in the displayedtime for all systems of that size. For smaller systems the pa-rameters in the algorithm could be set much more aggressive,resulting in an even faster simulated bifurcation that wouldstill result in the proper ground-state. In fig. 12 we showthe amplitudes of the conjugate variables describing posi-tion and momentum of the approximated Kerr-oscillators in Fig. 11. Shows the computation time of the SB-algorithm for in-creasing number of assets N . The different lines correspond to dif-ferent values of maximal investments per asset. Each point is es-timated by averaging 10 identical experiments to account for even-tual fluctuations. This was done for a fixed set of SB-parameters ε p = . , ∆ t = . , ∆ i = . the top panel and the maximized portfolio value in the bot-tom panel during the time in which the pumping amplitude p ( t ) is increased. In this setup, where the number of assetsare only 20, the optimal value has been reached before thepumping amplitude has reached its maximal value, theoreti-cally allowing us to reduce computational efforts whilst stillfinding the same state. Note that this effect also goes intothe other direction and for larger systems the properly tunedparameters can lead to an increase in computation time. Fig. 12. Top: Shows the average of the position amplitude x andmomentum amplitude y over the time in which the pumping amplitude p ( t ) is increased. Bottom: Shows the portfolio value during the firstpart of t . This is for N = , ˜ K = , γ = . , ε p = . and ∆ t = . . or large enough systems, or scenarios with extremevalues of µ and Σ , a fixed set of parameters will eventuallylead to meaningless results. As expected we observed that inscenarios where the steps of the pumping amplitude ε p andthe finite time increments ∆ t are picked too large, the algo-rithm will not end up in the ground-state. In such a scenariothe ’resolution’ of the algorithm needs to be increased andhence the respective computation time to simulate the adi-abatic evolution will increase as well. The exact parameterchoice has proven to be a delicate fine-tuning problem which,if done properly however, can also help reduce computationtime.Note that we have presented the performance measure-ments from an end-user perspective and will skip a detaileddiscussion around the bifurcation phenomenon and the prob-ability of success of the algorithm here. A highly valuablediscussion around such can be found in the original intro-duction [14] of the algorithm. Certain systems can have multiple configurations withIsing energy values, or equivalently portfolio values, that arenumerically very close to each other or even identical. Thelevel of degenerate and almost degenerate states depends onthe number of assets N , the amount of units to distribute ˜ K and on the expected return and risk of the assets. The amountof degenerate and almost degenerate configurations increasesfor larger values of ˜ K , particularly if the estimated returnsand risks are very similar among the assets. What we observeis that for systems with larger values of ˜ K , the SB-algorithmneeds a much more carefully picked set of parameters tofind the optimum. If the parameters of the SB-algorithm arenot tuned to an optimum, the algorithm seems to fail to de-tect those seemingly in-existing energy differences amongthe different configurations and ends up in one of the many(pseudo-)degenerate ground-states, depending on the initialrandom configurations of the oscillators. In the following wewill discuss an example where, without delicate fine-tuningof the SB-parameters, the SB-algorithm finds only close-to-optimal solutions, and where the increased level of degener-acy is expected to be strongly connected to the issue. Notehowever, that even for systems with a high level of degener-acy, we can still find close-to-optimal solutions for very large N and ˜ K faster than approaches from previous studies, as wasillustrated in fig. 11.To illustrate this we perform the portfolio optimizationat a fixed time point for 100 assets. If we limit the maxi-mal number of units to be spent per asset to just ˜ K / N = ε p ∈ [ . , . ] and ∆ t ∈ [ . , . ] all generate thesame results. The amount of different portfolios with sim-ilar values is comparably low due to the fact that only zero orone unit can be spent per asset. If we however increase thenumber of units to be spent to ˜ K = K / N = Fig. 13. N = , ˜ K / N = for various values of γ . Additionally random portfolios are shown.Fig. 14. N = , ˜ K / N = for various values of γ . Additionally random portfolios are shown. Calculation time for each pointobtained by the SB-algorithm is around 0.4 seconds. identical value increases. In such scenarios we can detectdeviations between the optimal solutions and those foundby the SB-algorithm if no specific fine-tuning of the SB-parameters is conducted. This is illustrated in fig. 14 wherethe SB-algorithm only finds close-to-optimal solutions with ε p = .
05 and ∆ t = .
02. This behaviour has a particularaccent for γ values that allow to balance risk and return, forextreme values of the risk aversion parameter however, thiswas not observed. In those limits also the number of portfo-lios with almost identical value decreases. Note that in bothfig. 13 and fig. 14, a set of random portfolios is plotted aswell. This however is for illustrative purposes only, sincethe number of possible combinations is extremely large andrandom samples cannot help us detect the entire universe ofpossible portfolios. What is nevertheless remarkable in thosetwo figures is the fact that without any fine-tuning the SB-algorithm finds close-to-optimal solutions in 0.4 seconds on single CPU out of 2 configurations and the global opti-mum even faster out of 2 possibilities.For smaller values of ˜ K / N , the parameter selection hasproven to be much more forgiving. Even for large systemsup to 1000 assets the SB-algorithm successfully finds the op-timal solutions. This is illustrated in fig. 15, where optimalportfolios with N = ε p ∈ [ . , . ] and ∆ t = [ . , . ] generates the same dis-played solutions. This flexibility allows to chose aggressiveSB-parameters, such that the computation of the optimal so-lution of systems of up to N = K / N = Fig. 15. Portfolio with N = assets, a total capital of ˜ K = units, where can be spent maximally per asset, either pur-chase or not. Computation time of the SB-approach is around 1 sec-onds per point, which can be pushed into the sub-second regime withthe mentioned parameter tuning. In the previous section 5.2 we discussed the behaviourof systems in the limit of zero trading costs. This is a par-ticularly useful limit since it allows us to test if the optimaltrajectory coincides with the set of individually optimizedportfolios at each point in time. Some trajectories found bythe SB-agorithm however, have not shown the success ob-served in fig. 8 and without careful fine-tuning do not co-incide exactly with the global optimum at c =
0. For thosecases this strongly indicates that the corresponding trajecto-ries obtained for non-zero trading costs are also only close-to-optimal. First preliminary experiments suggest that this isnot directly connected with the number of combinations, butagain to the amount of trajectories with values very close tothe optimum.Note that this numerical cross-check with the trajec-tory at zero trading cost can also assist in tuning the SB-parameters. If the optimal solution is not known due tothe high amount of possible configurations, a set of SB-parameters is selected that reproduces the reference config-urations in the limit of vanishing trading costs. Note that the trading costs are elements of the interaction matrix inthe Ising formulation, and hence this approach is assumed towork only for small enough values of c .In this section we have identified that without properfine-tuning of the SB-parameters we can get solutions thatare only close to the optimum. This was however only ob-served for systems that show a large amount of portfoliosor trajectories with values very close to the optimum. Wehave yet not established a proper dynamic parameter selec-tion framework and resorted to manual adjustments in thisinitial phase of research. A mathematically rigorous inves-tigation is necessary to understand the connection betweenthe algorithm’s parameter and the distribution of the differ-ent portfolio-values and the magnitude of their numerical dif-ferences. This is subject of ongoing investigations and willserve as the main key to construct a dynamic SB-parameterselection framework that allows for optimal portfolio opti-mization. We have shown that the SB-algorithm can successfullybe used to find optimal solutions for the integer portfolioproblem as well as the integer trading trajectory problem.We have investigated portfolios with up to 1000 assetsand our first performance investigation on a single desktopCPU has already confirmed the power of the SB-algorithmfor its novel use-case as a portfolio and trading trajectoryoptimizer. For the investigated portfolios the computationtime does not exceed a couple of seconds and truly showsan unprecedented speed in finding optimal and close to opti-mal solutions. In a next step we will move the optimizationframework to our GPU cluster for parallel computing andhope to share results in the near future that capture a sub-stantial amount of the actively traded assets.The presented formulation in the Ising representation in-corporates the upper quantity limits as hard constraints andwould also allow for individual asset dependent upper quan-tity constraints. With a proper penalty term in the Ising de-scription we are able to control the costs of rebalancing theportfolio from one point in time to the next, which can alsostraightforwardly be made asset and time dependent.For a large number of assets, in combination with alarge number of units/capital to distribute per asset, theSB-algorithm has detected optimal but sometimes also onlyclose-to-optimal solutions. For larger systems with a largeamount of almost degenerate ground-states, the global op-timum was often only detected after careful fine-tuning ofthe SB-parameters. Note however, that the rudimentary stan-dard set of parameters that was successfully used across awide range of problems already allowed to obtain solutionsvery close to the optimum. We have not established a frame-work which maps the problem settings to a suitable set ofSB-parameters, that guarantee finding the ground-state. Inthis work we have resorted to manual adjustments of this im-portant and delicate step, and consider a detailed parameterdiscussion the main key for further improvements in this di-rection.nother future challenge is the introduction of furtherconstraints, such as cardinality or fixed quantity constraints.The inclusion of further constraints can however quickly leadto the introduction of many ancillary spins and also has thepotential to disrupt the evolution of the SB-algorithm into theground-state. The constraints beyond the hard-wired limita-tions presented in this work are therefore considered one ofthe biggest challenges to make the SB-algorithm applicablefor a variety of real world constraints in portfolio optimiza-tion problems.The presented formulation allows to control also timeand asset dependent trading activities, lacks however aproper translation from real word trading concepts to con-crete numerical penalty values in the Ising formulation.There is therefore a strong interest to develop a more elab-orate trading cost framework and to introduce rebalancingscenarios that are closer to real world applications.Despite the open challenges, this work has shown thefirst highly successful and incredibly fast portfolio optimiza-tion with the simulated bifurcation algorithm. An approachthat we believe will see a wide range of applications in manyother fields of finance as well.
References [1] H. M. Markowitz
Portfolio selection
J. Financ. 7(1952), no. 1, 7791.[2] J. Wiley, Sons Inc.
Portfolio selection: Efficient diversi-fication of investments
Cowles Foundation for Researchin Economics at Yale University, Monograph 16, NewYork, 1959.[3] J. B. Blackwell
Mean-variance analysis in portfoliochoice and capital markets
Oxford, 1987.[4] D. Bienstock
Computational study of a family of mixed-integer quadratic programming problems
Math. Pro-gram 74 (1996), 121-140[5] D.X. Shawa, S. Liub, and L. Kopmanb
Lagrangian re-laxation procedure for cardinality-constrained portfo-lio optimization
Optim. Method. Softw. 23 (2008),no.3, 411420.[6] H. Kellerer, R. Mansini, and M. G. Speranza
Select-ing portfolios with fixed costs and minimum transac-tion lots
Annals of Operations Research, vol. 99, no.1-4, pp. 287 304, 2000[7] R. Mansini and M. G. Speranza
Heuristic algorithmsfor the portfolio selection problem with minimum trans-action lots
European Journal of Operational Research,vol. 114, no. 2, pp. 219233, 1999.[8] N. J. Jobst, M. D. Horniman, C. A. Lucas, G. Mitra etal
Computational aspects of alternative portfolio selec-tion models in the presence of discrete asset choice con-straints
Quantitative finance, vol. 1, no. 5, pp. 489501,2001[9] J. P. Vielma, S. Ahmed, and G. L. Nemhauser
A liftedlinear programming branch-and-bound algorithm formixedinteger conic quadratic programs
INFORMSJournal on Computing, vol. 20, no. 3, pp. 438450, 2008.[10] M. Corazza and D. Favaretto
On the existence of so- lutions to the quadratic mixed-integer meanvarianceportfolio selection problem
European Journal of Oper-ational Research, vol. 176, no. 3, pp. 19471960, 2007[11] P. Bonami and M. A. Lejeune
An exact solutionapproach for portfolio optimization problems understochastic and integer constraints
Operations research,vol. 57, no. 3, pp. 650670, 2009.[12] G. Rosenberg P. Haghnegahdar P. Goddard, P. Carr,Kesheng Wu, M. L. de Padro
Solving the Optimal Trad-ing Trajectory Problem Using a Quantum Annealer
IEEE Journal of Selected Topics in Signal Processing(JSTSP), Volume 10, Issue 6, 2016.[13] J. Cohen, A. Khan, C. Alexander
Portfolio Optimiza-tion of 40 Stocks Using DWaves Quantum Annealer https://arxiv.org/pdf/2007.01430.pdf, (2020)[14] H. Goto, K. Tatsumura and A. R. Dixon
Combinatorialoptimization by simulating adiabatic bifurcations innonlinear Hamiltonian systems
Sci. Adv. 5, eaav2372(2019).[15] L. McMahon, A. Marandi, Y. Haribara, R. Hamerly,C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S.Utsunomiya, K. Aihara, R. L. Byer, M. M. Fejer, H.Mabuchi, Y. Yamamoto
A fully-programmable 100-spin coherent Ising machine with all-to-all connections
Science 354, 614617 (2016)[16] L. Andrew
Ising formulations of many NP problems
Frontiers in Physics, 2014 Operations research, vol. 2[17] E. Ising ”Beitrag zur Theorie des Ferromagnetismus”
Zeitschrift f¨ur Physik, vol. 31, issue 1, pp. 253-258(1925)[18] P.W. Anderson
Plasmons, gauge invariance, and mass
Physical Review. 130 (1): 43942[19] F. Englert, R. Brout
Broken symmetry and the mass ofgauge vector mesons
Physical Review Letters. 13 (9):32123.[20] P. W. Higgs
Broken symmetries and the masses ofgauge bosons
Physical Review Letters. 13 (16): 50809[21] F. Barahona
On the computational complexity of Isingspin glass models
J. Phys. A 15, 32413253 (1982)[22] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa,H. Tamura, H. G. Katzgraber
Physics-inspired opti-mization for quadratic unconstrained problems usinga digital annealer arXiv:1806.08815.[23] Z. Bian, F. Chudak, R. Israel, B. Lackey, W. G.Macready, A. R. Miyaza
Mapping constrained opti-mization problems to quantum annealing with applica-tion to fault diagnosis
Frontiers in ICT, 3, 14.[24] H. Goto
Bifurcation-based adiabatic quantum compu-tation with a nonlinear oscillator network
Sci. Rep. 6,21686 (2016)[25] B. Leimkuhler, S. Reich