Benchmarking Quantum Annealing Controls with Portfolio Optimization

Quantum annealing offers a novel approach to finding the optimal solutions for a variety of computational problems, where the quantum annealing controls influence the observed performance and error mechanisms by tuning the underlying quantum dynamics. However, the influence of the available controls is often poorly understood, and methods for evaluating the effects of these controls are necessary to tune quantum computational performance. Here we use portfolio optimization as a case study by which to benchmark quantum annealing controls and their relative effects on computational accuracy. We compare empirical results from the D-Wave 2000Q quantum annealer to the computational ground truth for a variety of portfolio optimization instances. We evaluate both forward and reverse annealing methods and we identify control variations that yield optimal performance in terms of probability of success and probability of chain breaks.


INTRODUCTION
Optimization is integral to many scientific and industrial applications of applied mathematics including verification and validation, operations research, data analytics, and logistics, among others [1,2]. In many cases, exact methods of solution, including stochastic optimization and quadratic programming, are computationally intractable and novel heuristics are used frequently to solve problems in practice [3]. Quantum annealing (QA) offers a novel meta-heuristic that uses quantum mechanics for unconstrained optimization by encoding the problem cost function in a Hamiltonian [4,5]. Recovery of the Hamiltonian ground state solves the original optimization problem and this approach has been mapped to a variety of application areas [6][7][8][9]. Several experimental efforts have realized quantum annealers [10][11][12], and application benchmarking of these systems has shown QA is capable of finding the correct result with varying probability of success [13][14][15][16][17][18][19].
QA performance depends implicitly on the complexity of the underlying problem instance as well as the controls that implement the heuristic [20,21]. Presently, there are multiple controls available to program quantum annealers that may each impact the observed probability of success. Notionally, the controls may be categorized as pre-processing, annealing, and post-processing methods. Whereas pre-processing controls define the encoded Hamiltonian and embedding onto the quantum annealer [22,23], the annealing controls drive the time-dependent physics of the device and the underlying quantum state [20,24] while post-processing controls influence the readout and decoding of the observed results [25,26]. Collectively, the choice for each type of control may either enhance or impede the probability of reaching the encoded ground state and, therefore, impact the resulting solution state.
Here we benchmark a selection of pre-processing and annealing controls available in a programmable quantum annealer [10] using a well-defined class of unconstrained optimization problems derived from the application of Markowitz portfolio theory [27]. As a variant of binary optimization, Markowitz portfolio optimization selects the subset of investment assets expected to yield the highest return value and minimal risk while staying within a total budget constraint [27,28]. We cast this problem which forms a complete graph as unconstrained optimization and benchmark the probability of success for QA to recover the global optimum. In particular, we benchmark the pre-processing and annealing controls available in the 2000Q, a programmable quantum annealer from D-Wave Systems [10]. This includes controls for mapping the logical problem onto hardware and scheduling the annealing process. We gather insight into the underlying dynamics using multiple measures of success tested across an ensemble of randomly generated instances of portfolio optimization.
Previous research has benchmarked QA in comparison to classical heuristics for solving various optimization problems [14,29,30]. In particular, several variations of portfolio optimization have been used to benchmark QA performance [20,31,32]. Rosenberg et al. demonstrated several encodings of a multi-period Markowitz portfolio optimization formulation to be solvable using a quantum annealer and found promising initial results in probability to find the optimal result [32]. Venturelli et al. benchmarked a similar mean-variance model of portfolio optimization using a hybrid solver that couples quantum annealing with a genetic algorithm [20]. This hybrid algorithm was found to be 100x faster than forward annealing alone. In this work, we present a formulation of portfolio optimization to benchmark the behaviour of QA controls. We present studies focused on the variability in success with respect to available quantum annealing controls in an attempt to establish a methodology for finding an optimal set of controls which yield the highest solution quality [33,34].
The presentation is organized as follows. In Sec. II, we review quantum annealing and the the available controls. In Sec. III, we provide an overview of the benchmarking methods and the use of Markowitz portfolio selection for problem specification. In Sec. IV, we present the results from experimental testing using different quantum annealing controls with the 2000Q. We offer conclusions in Sec. V.

II. QUANTUM ANNEALING
Under ideal conditions, forward annealing evolves a quantum state |Ψ(t) under the time-dependent Schrödinger equation where T is the total forward annealing time and the timedependent Hamiltonian is where s(t) ∈ [0, 1] is the control schedule and timedependent amplitudes A(s) and B(s) satisfy the conditions A(0) B(0) and A(1) B(1). We consider the initial Hamiltonian H 0 = − n i σ x i as a sum of Pauli-X operators σ x i over n spins. The final Hamiltonian H 1 represents the unconstrained optimization problem with a corresponding ground state that encodes the computational solution. We will consider below only problems represented using the Ising Hamiltonian where h i is the bias on the i th spin, J i,j is the coupling strength between the i th and j th spin, σ z i is the Pauli-Z operator for the i th spin, and β is a problem-specific constant. The Ising Hamiltonian is well known for representing a variety of unconstrained optimization problems [35].
Instantaneous eigenstates at time t are defined as where j ranges from 0 to N − 1 with N = 2 n the dimension of the Hilbert space. For an initial quantum state prepared in the lowest-energy eigenstate at time t = 0, i.e, |Ψ(0) = |Φ 0 (0) , adiabatic evolution under the Hamiltonian in Eq.
(2) to time T will prepare the final state |Ψ(T ) = |Φ 0 (T ) with high probability provided T is sufficiently large. In particular, the evolution must be much longer than the inverse square of the minimum energy gap between the ground and first excited states [4]. At time T , the prepared quantum state is measured in the computational basis to generate a candidate solution for the encoded problem. Another variation of quantum annealing reverses the time-evolution process by beginning in an eigenstate of H 1 . Known as reverse annealing, the initial quantum state evolves under Eq. (2) in the reverse direction. The Hamiltonian starts as H 1 at time t = 0 and evolves backward to a point s p in the control schedule that corresponds to time t 1 . The Hamiltonian then pauses for a time t p = t 2 − t 1 before evolving in the forward direction from the value s p at time t 2 back to the final Hamiltonian at time T , where the latter time represent the reverse annealing time. The control schedule for reverse annealing is then defined as [36,37] The differences in the control schedules of forward and reverse annealing are demonstrated in Fig. 1, where a linear reverse annealing schedule is compared to a linear forward annealing schedule using the amplitudes A(s) = (1 − s) and B(s) = s. Notably, forward annealing controls increase monotonically with time whereas reverse annealing controls include a change in the direction of the control schedule where the ramp time from s = 1 to s p is t r = t 1 , the time paused at s p is t p , and the quench time back from s p to s = 1 is t q = T − t 2 .

A. Quantum Annealing Controls
In practice, there are non-ideal behaviours that arise in practical implementations of quantum annealing. Equations (1)-(5) represent quantum annealing under ideal adiabatic conditions that are difficult to realize in actual quantum devices. Real-world quantum annealers have limits in the ability to control the Hamiltonian and quantum dynamics [38]. In addition, the presence of illcharacterized environmental couplings give rise to flux noise [39]. The imperfect setting of the Hamiltonian parameters (h, J i,j ) by the analog control circuits gives rise to a small intrinsic control error [34]. These errors undermine the accuracy of the physical hardware [22,38]. Finally, annealing too quickly may violate the essential adiabatic condition [4], while annealing too slowly may lead to undesired thermal excitations of the quantum state due non-zero temperature fluctuations [40]. This multitude of effects complicates both the description of quantum annealing as well as the assessment of its performance.
Given the implicit dependence on several competing factors, a variety of strategies have emerged for controlling quantum annealing to maximize probability of success in recovering the ground state and minimizing errors in the quantum computational solution. These control strategies include efficiently mapping the problem Hamiltonian onto the physical hardware Hamiltonian, tuning annealing schedule, applying variable transformations to mitigate control biases, and using reverse annealing to refine initial solutions [34,36].
We investigate a subset of controls available in the D-Wave 2000Q, a programmable quantum annealer composed from an array of superconducting flux qubits operated at cryogenic temperatures [41]. The 2000Q consists of up to 2048 physical qubits arranged in a sparsely connected array whose governing Hamiltonian is described by a time-dependent, transverse Ising Hamiltonian [42] for which with the Hamiltonian parameters in the device can be programmed individually. This enables a broad variety of computational problems, including portfolio optimization, to be realized. We briefly review some of the controls available to influence the success of solving these problems using quantum annealing.

Problem Embedding
The Hamiltonian encoding the computational problem must be mapped into the physical hardware while satisfying the constraints of limited connectivity. The 2000Q hardware supports a sparse Chimera graph in which physical qubits are not fully connected but have average degree 6. A fully connected graph, like in Fig. 2, must be mapped onto the more sparse Chimera graph. A single spin from the input Hamiltonian may be realized in hardware using multiple physical qubits that form a strongly interacting representative chain of spins. By judiciously choosing these chains and their interactions, the original input Hamiltonian may be constructed. This process, known as embedding, depends on the input problem as well as the target hardware connectivity. In general, embedding is NP-hard for arbitrary input graphs [43], and there are upper limits on the maximum graph that can be embedded [44]. For example, the largest fully connected problem that can be embedded onto the 2000Q has ∼ 60 spins, while the limit in practice depends on the number of faulty/inactive physical qubits in the device.
Embedding algorithms that optimize chain length may greatly reduce the number of physical qubits required by considering problem symmetry as well as the location of faults in the hardware. We highlight two embedding algorithms widely used in programming the 2000Q. The first method by Cai, Macready, and Roy is based on randomized placement and search for the weighted shortest path between spin chains [45]. This method, which we denote as CMR, applies to arbitrary input graphs but typically creates a distribution of chain lengths. By contrast, a second method by Boothby, King, and Roy based on a clique embedding typically generates shorter and uniform chain lengths of size for n logical spins [46]. A representative example of the output from these different methods is shown in Fig. 2 using a fully connected problem with 20 logical spins. Both methods are available in the D-Wave Ocean software library [47].
Ensuring an embedded chain of qubits collectively represents a single logical variable requires an intra-chain coupling that is larger in magnitude than the the interchain couplings between chains. In other words, the chain of physical qubits must be strongly coupled to remain a single logical spin. However, it is possible for chains to become "broken" in so far as individual physical spins within the chain differ in their final state. In general, chain breaks arise from non-adiabatic dynamics that lead to local excitation out of the lowest energy state with longer chains more susceptible to these effects [34,48].
An additional control is required for decoding embedded chains to recover the computed logical spin state. In the absence of chain breaks, the logical value is inferred directly from the unanimous selection of a single spin state by every physical qubit. In the presence of chain breaks, several strategies may be employed to decide the logical value including majority vote [34], which selects the logical spin value as the value that occurs with the highest frequency in a chain.  Figure a) is complete K20 graph which is fully connected with 20 nodes and 190 edges where each node represents a logical spin and each edge is a coupling between spins. Figure b) is the CMR algorithm which requires the allocation of 23 unit cells and c) is the clique embedding algorithm which requires the allocation of 15 unit cells. The nodes represent physical qubits, lines are the couplings between physical qubits, and each color is a different physical spin chain corresponding to a logic spin.

Spin Reversal
Interactions between embedded chains arise from the required coupling between the logical spins. However, imperfections in the control of these spins lead to small biases that can become non-negligible for larger qubit chains and contribute to the complex dynamics describing the device. In turn, the probability for finding the expected ground state solution can decrease do to these bias errors. The influence of these errors on the compu-tational result may be mitigated by using spin reversal transforms to average out biases.
As a gauge transformation, spin reversal redefines the Hamiltonian by replacing the biases and couplings for a subset of spins with their negated value [33,34]. This transformation maintains the ground state of the logical problem. However, this transformation flips the sign of randomly selected qubits so that on average their bias is reduced. This strategy mitigates errors on individual spins by balancing the noise on the device prior to annealing [26]. The number of selected spins as well as the parameter g that defines the number of times to perform the transformation.

Annealing Schedules
Tailoring the annealing amplitudes A(s) and B(s) is perhaps the most direct method to control forward annealing. The annealing schedules control the rate of change of the H(t), which must be sufficiently slow to approximate the adiabatic condition [49]. An example of the amplitudes in a D-Wave 2000Q is shown in Fig. 3. While forward annealing on the D-Wave 2000Q, A(s(t)) >> B(s(t)) at t = 0, A(s(t)) decreases and B(s(t)) increases for 0 < t < T , and B(s(t)) >> A(s(t)) at t = T . The optimal annealing time is problem dependent and inversely proportional to the minimum energy gap [4], and, in general, the value and position of the minimum energy gap for a given H(t) is typically unknown and hard to identify. Extending the annealing time T arbitrarily long may not only be limited by hardware parameters but also be counter-productive due to competing thermal processes that depopulate the ground state [25,50]. There is an upper limit to the total job time (N s T ≤ 1 s) as well as total annealing time (T ≤ 2 s) on the D-Wave 2000Q.
When reverse annealing, the three primary parameters for control are the initial state e i , the pause point s p , and the pause duration t p . The times t r and t q can also be manipulated, but we keep these constant and symmetric for our experiments. Reverse annealing uses e i to set the initial quantum state, which may be based on the output of forward annealing, a heuristically computed candidate, a random state or other methods. Our experiments use a pre-computed initial state, e.g., using forward annealing. An iterative procedure is then used which replaces the e i of each subsequent reverse annealing sample with the output from previous reverse annealing iteration.

B. Quantum Annealing Metrics
We characterize quantum annealing using the probability of success defined as the overlap of the final, potentially mixed quantum state ρ prepared by QA with the pure-state describing the expected computational outcome Φ 0 (T ). Empirically, the probability of success is estimated from the frequency with which the observed solution state matches the expected outcome. When the expected ground state solution is known, we define the statistic δ i = 1 if the i-th sample matches the known ground state and δ i = 0 if it does not. For the k-th problem Hamiltonian instance, the estimated probability of success is then defined asp where N s is the total number of samples. The average over an ensemble of N p problem instances is defined as A second metric for characterizing quantum annealing performance, and especially the non-adiabatic dynamics, is the number of chain breaks observed in the recovered solution samples. As noted above, a chain break is observed when the chain of physical qubits embedding a logical spin has more than one unique spin value. We estimate the probability of chain breaks for a problem instancep where the statistic i = 1 when the i-th sample solution contains at least one broken chain for any of the logical spins and i = 0 when no embedded chain is broken. The average probability of chain breaks over an ensemble of N p problem instances is then defined as It is important to note that the effects of chain breaks can be mitigated by post-processing methods, such as majority vote, which make hard decisions on the logical spin value. While the above metrics quantify the probability with which quantum annealing recovers the correct solution, additional information about computational performance comes from the distribution of all solution samples obtained. In particular, the distribution over sample energies provides a representation for the weight of errors in the solution samples. A distribution concentrated around the lowest energy indicates a small number of errors in the computed solutions, while a broad or shifted distribution hints at a larger number of errors. We denote the energy computed from the i-th solution sample as E(i) and we define the j-th energy bin as h j . The bin h j counts the number of samples with an energy in the range [j, j +1]∆ where ∆ controls the granularity of binning the energies. The resulting set {(j∆, h j )} approximates the energy distribution of the sampled solutions.

III. BENCHMARKING METHODS
We benchmark performance of the quantum annealing controls presented in Sec. II using a variant of constrained optimization derived from Markowitz portfolio theory. We recast this problem as unconstrained optimization before reducing to quadratic unconstrained binary optimization (QUBO) form. The latter form is easily translated to the classical Ising spin Hamiltonian and, subsequently, to the problem Hamiltonian defined by Eq. (3).

A. Markowitz Portfolio Selection
Portfolio optimization selects the best allocation of assets to maximize expected returns while staying within the budget and minimizing financial risk. The Markowitz theory for portfolio selection focuses on diversification of the portfolio for risk mitigation [27]. Instead of allocating high percentages of a budget toward assets with the highest projected returns, the budget is distributed over assets that minimize correlation between the asset's historical prices. In this model, the covariance between purchasing prices serves as a proxy for risk in which positively correlated assets are considered to be more risky. We review the methods by which the benchmark problems are generated and solved in this section.
We consider Markowitz portfolio optimization as a quadratic programming problem that determines the fraction of available budget b to allocate toward purchasing assets with the goal of maximizing returns while minimizing risk. By selecting a partition number w, the fraction p w = 1 2 (w−1) represents the granularity of the partition. The portfolio optimization problem selects how many of those partitions to allocate toward each asset with an integer z u . Thus, the fraction of b to invest in each u th asset is given by p w bz u , and portfolio optimization identifies how much of the m assets to select given the budget b and a risk threshold c. Thus, portfolio selection is cast as where for the u th asset r u is the expected return and c u,v is the historical price correlation between assets u, v. In Eq. (12), the first term represents maximization of the expected returns over the available assets. There are many methods for forecasting expected returns, e.g., based on market price, expert judgement, and historical price data [51,52]. For simplicity, we model expected returns as whereā u is the average of a u , the history of price data for the u th asset. The first constraint in Eq. (12) places a hard constraint on the total allocation of assets to sum to b. We emphasize that this constraint penalizes portfolios that do not allocate the entire budget as well as those that over commit. Finally, the second constraint accounts for diversification by asserting that the sum of covariance between asset prices c u,v be less than or equal to the risk threshold c. The historical price covariance is calculated as the correlation between pairs of assets by comparing the p w fraction of each asset's historical price data. Here covariance is defined as where a u,l is the l th historical price value for asset u and N f is the number of price points in the historical data. We solve this variation of Markowitz portfolio selection using quantum annealing by casting the formulation in Eq. (12) into quadratic unconstrained binary optimization (QUBO). We express the integer variable z u as a w-bit binary expansion with x i ∈ {0, 1} and the composite index i(u, k) = (u − 1)w + k. The expected returns are then expressed as while the allocation constraint becomes the penalty term We consider a correlation threshold c = 0 such that the correlation constraint becomes Our formulation of Markowitz portfolio selection as an unconstrained optimization problem then becomes where the problem size n = mw, r i = 2 k−1 r u , c i,j = 2 k−1 2 k −1 c u,v , and θ 1 , θ 2 and θ 3 are Lagrange multipliers used to weight each term for maximization or penalization.
For purposes of benchmarking, we generate an ensemble of problem instances by sampling from uniform random price data with a seed of b/5 . A random number is drawn as the initial price a u,1 and every subsequent historical price point up to the purchasing price is −25% to +25% of the previous price a u,l . The price range was set to be between b/10 and b with N f = 100 historical price points per asset. In addition, we normalize all a u,l by a u,N f to keep all asset prices to a similar range.
We set θ 1 = 0.3, θ 2 = 0.5, θ 3 = 0.2 in the problem instances where θ 2 is set higher to enforce the budget constraint. These weights were chosen after testing which combination stayed on budget and gave some diversity. By keeping θ 2 constant and increasing θ 3 while decreasing θ 1 , an investor could increase the diversity relative to the potential returns and vice versa when decreasing θ 3 relative to θ 1 . We generate 1000 problems for each problem size with m = 2, 3, 4, 5 assets and w = 4 slices.

B. QUBO to Ising Hamiltonian
We formalize the unconstrained portfolio optimization problem in Eq. (19) to quadratic unconstrained binary optimization (QUBO) as where q i is the linear weight for the i th spin, Q i,j is the quadratic weight for interactions between the i th and j th bits, and γ is a constant. Note that our definition of QUBO expresses optimization as minimization by switching the sign of Eq. (19) to be consistent with the use of quantum annealing to recover the lowest-energy state. The corresponding relationships with the original problem instance are given as Similarly, the quadratic binary form may be reduced to a classical Ising Hamiltonian where spin s i ∈ {−1, 1} is defined by s i = 2x 1 − 1 with s = (s 1 , s 2 , . . . , s n ) while h i is the spin weight, J ij is the coupling strength, and β is a problem-specific constant. The parameters for the Ising Hamiltonian are given as The classical Ising formulation is then converted into a corresponding quantum Ising Hamiltonian given by Eq. (3) using the correspondence s i → σ z i .

C. Computational Methods
We used a D-Wave 2000Q quantum annealer for our experiments. We calculate the probability of success, the probability of chain breaks, and the energy distribution across each problem instance. For each instance, we estimated these metrics by collecting N s = 1000 samples of the computed solution. We used D-Wave's solver API (SAPI) with Python 2.7 to solve each instance of Markowitz portfolio selection using the hardware controls outlined in Sec. II A. We ran 1, 000 samples per problem over a set of 1, 000 problems for forward annealing examples an 100 problems for revere annealing examples. We implement the majority vote post-processing technique for any broken chains to interpret raw solutions returned by the 2000Q. The program implementation and data sets collected from these experiments are available online [53].
For benchmarking purposes, we also solved each problem instance using brute force search for the minimal energy solutions of the QUBO formulation. We computed the complete energy spectrum for each portfolio instance. These energy spectrum and the corresponding states were then used as ground truth for testing the accuracy of results obtained from quantum annealing. By sorting the spectrum, we benchmarked the success of reverse annealing using initial states e i sampled from these different parts of the spectrum.

IV. RESULTS
We benchmark quantum annealing controls by evaluating their influence on the probability of success and probability of chain breaks across problem instances. We first characterize how problem parameters influence the baseline performance by estimating the probability of success for forward annealing using T = 15 µs, g = 0, and a randomized embedding strategy. As shown in Fig. 4, we comparep s for two cases of w = 1 and w = 4 across increasing n. The estimated probability of success for problems with w = 4 is consistently higher for problems with no slicing. These results are explained by the energy spectra for the different problem parameters, which indicate sharp differences in the density of states. As shown in Fig. 5, a typical problem instance with w = 4 has a much higher density of states than those with no slicing (w = 1). Intuitively, the single-slice behavior results from the specification that the price for each asset is proportional to budget, and, therefore, only a single asset may be selected Consequently, the probability to recover the lowestenergy state competes with these closely spaced, higher energy solutions, which leads to a corresponding decrease in the probability of success. For the remaining benchmark tests below, we chose w = 4 as it represents a more challenging test for the quantum annealer as well as a greater interest to real-world financial applications.
A. Benchmarking Forward Annealing Controls

Embedding
Embedding generates and places the physical spin chains for each logical spin on the quantum annealing hardware. We evaluated the CMR and clique embedding algorithms described in Sec. IV A 1 by estimating the probability of success across problem sizes of m = 8, 12, 16, and 20 logical spins. For all problem instances of a same problem size, we use the same embedding because they require the same number of fully connected logical spins. We set the parameters of the embedded Ising Hamiltonian by scaling the inter-chain couplings J i,j to lie in the range [−1, +1]. We scale all J i,j using a rescale factor of 1 jmax where j max is the largest J i,j so all embedded J i,j = 1 jmax J i,j . This scales all J i,j to be between + − 1. The intra-chain coupling strength is set to −1 to have a negative bias stronger than the J i,j values which range −10 −1 ≤ J i,j ≤ 10 −1 due to our data generation and normalization techniques.
The average chain length l c from CMR and clique embedding methods grows with the number of logical spins n. The average is computed with respect to all chains in an embedding and plotted with respect to n in Fig. 6. As expected by Eq. (6), the clique embedding method has a uniform chain length for each n. By contrast, the CMR method generates chains of variable length as indicated by the the average chain length and variance shown in the plot. From each of the embedding methods, we estimate the probability of success and probability of broken chains. As shown in Fig. 7, we observe very small differences in both metrics with increasing problem size. From fitting the resulting point to an exponential, we findp s decays sub-exponentially with respect to n with rate −0.523 for the CMR embedding and rate −0.528 for the clique embedding. We find thatp b grows at a sub-exponential rate of 0.1824 for CMR embedding and 0.1656 for clique embedding as n increases. There is not a significant difference in thep s performance between CMR and clique embedding, but clique embedding requires a fewer number of spins as n increases and shows a slight improvement inp b . Therefore, we chose to use clique embedding for subsequent benchmarks.

Forward Annealing Time
According to the adiabatic theorem, forward annealing more slowly should increase the probability of the system remaining in the ground state and thus increase the probability of success. We varied the forward annealing time T from 1 µs to 999µs, which is the broadest range accessible on the D-Wave 2000Q. As shown in the upper panel of Fig. 8, we observed statistically insignificant changes in the probability of success as annealing time increased at each problem size. Fitting the average probability of success with respect to problem size for the annealing time T = 100 µs, yields a sub-exponential decay rate for p s given by −0.528 and a sub-exponential growth rate for p b given by 0.1628 as n increases. We do observe a statistically significant difference in the estimated probability of chain breaksp b with respect to forward annealing time as shown in the lower panel of Fig. 8. For T = 100 µs, we recover a growth rate of 0.1656 for the probability of chain breaks with respect to problem size.

Spin Reversal
As discussed in Sec. IV A 1, embedding maps a logical spin to many physical spins by creating strongly coupled chains. Coupling of these embedded spins via J i,j in Eq. (3) can lead small biases that may be amplified by imperfections in the hardware. A spin reversal transform mitigates against bias errors by reversing the sign of a spin in the Ising Hamiltonian. This transform preserves the logical problem but reverses the bias error on the embedded spin chain. By randomly selecting a subset of spins to revise, we evaluate the influence of spin-reversal transform on the probability of success. We use g transforms when estimating the probability of success for a given problem instance, such that there are N s /g samples per transform. We observed nominal improvements in Fig. 9 by using at least g = 2 with no advantage to using g > 2. For g = 2, we observe an sub-exponential decay rate of −0.505 forp s and a sub-exponential growth

B. Benchmarking Reverse Annealing Controls
From the reverse annealing controls listed in Sec. II A 3, we designed three experiments based on the e i for the reverse annealing heuristic that include (i ) starting in the known ground state e 0 , (ii ) starting in the known first excited state e 1 , and (iii ) starting in the lowest-energy state obtained from 1000 forward annealing samples e f . We then sweep over various schedules to find the optimal s p with a range of [0.1, 0.9] and t p with a range of [15 → 800]µs. The t r and t q parameters were set to be constant and symmetric at 5µs each. Thus, the total anneal time is T = t r + t p + t q where t p is the time parameter that we chose to analyze. For all experiments, we ran the reverse annealing iterative heuristic with 1000 samples for 100 random problems were also used in the forward annealing experiments. We estimated the probability of success for reverse annealing with respect to different choices for e i , s p , and t p . We compared the combined heuristic of reverse annealing with forward annealing to forward annealing alone withp s ,p b , as well as the frequency of finding energies in excited states to forward annealing alone [54].
By setting e i to the ground state, we tested for parameters s p and t p that decreasep s when the quantum annealer is fed the correct solution. For this experiment, p s can be thought of as the probability of staying in e 0 where p f is the probability that forward annealing found the ground state, α i ∈ {0, 1} indicates whether forward annealing found the ground state for the i th problem prior to reverse annealing, and δ j ∈ {0, 1} is a variable indicating whether the j th sample of the i th problem was measured to be the ground state with reverse annealing. By setting e i = e 1 , we tested whether reverse annealing enhances the probability to populate the ground state. For these tests,p s estimates the probability of moving from an excited state to the ground statep s (e e → e 0 ) = (1 − p f ) * p s (27) ( In addition to testing reverse annealing at e i = e 0 and e 1 , We tested reverse annealing in combination with forward annealing for whichp s estimates the cumulative probability of finding the correct solution state. For these experiments, we found it useful to primarily analyzep s (R) −p(e 0 → e 0 ) =p(e e → e 0 ) to determine if reverse annealing improved upon thep s of forward annealing. The results from setting e i = e 0 for each problem with a problem size of n = 20 where m = 5 and w = 4 is shown in Fig. 10. Because the computation begins in the correct solution state, this test measures the probability by which reverse annealing introduces errors into the correct solution. Ideally,p s will remain near unity for all s p and t p . We observe that reverse annealing causes the system to leave the ground state withp s reducing to on the order of 10 −5 by annealing back to at least s = .6 and increasing t p ≥ 200µs The results from setting e i = e 1 with a problem size of n = 20 where m = 5 and w = 4 for each problem is shown in Fig. 11. A maximal value of 4.8 × 10 −4 for p s is found with parameters s = 0.7 and t p = 800 µs. This is is ap s one order of magnitude higher than what is observed with forward annealing. This suggests that if e i is very close to e 0 , there may be some benefit to choosing reverse annealing over forward annealing. When solving optimization problems for applications in practice, the ground state and excited state will be unknown. However, one approach is to use reverse annealing in addition to forward annealing by using the lowest energy state found with 1, 000 forward annealing samples e f as e i for another 1, 000 samples of reverse annealing. The next experiment tests whether reverse annealing used in combination with forward annealing increasesp s with a problem size of n = 20 where m = 5 and w = 4 . The experimental results from setting e i = e f is shown in Fig. 12. These tests were constructed to determine when combining reverse annealing with forward annealing can improve upon forward annealing. Therefore, we removed the 6 problems forward annealing provided an e i = e 0 and thusp s for this experiment is given bỹ p s (R) − p(e 0 → e 0 ) in this analysis. Similar to the previous experiment in Fig. 11, thep s is at best on the order of 10 −4 at parameters s = 0.7 and t p = 400µs which is one order of magnitude greater than the forward annealing experiments.   12 shows a potential for reverse annealing to improve upon results found with forward annealing inp s . Therefore, we take a set of 100 problems solved with reverse annealing and forward annealing and compare thẽ p s of forward annealing (orange) alone to thep s of reverse annealing alone (blue) to thep s with a selection of either forward annealing or reverse annealing (green). If for a problem forward annealing found at least one ground state the forward annealingp s was plotted for that problem (6 problems) and otherwise the reverse annealingp s was plotted (94 problems). Thep s is measured over n ranging from [8,20]. The reverse annealing parameters are set to have an e i = e f , s = .7, and t p = 400µs. As shown in Fig. 13, we observe that when taking the combination of best results from forward annealing and reverse annealing with e i = e f , we get ap s that improves by an order of magnitude over forward annealing alone for n = [16,20] with a sub-exponential decay at a rate of −0.309. Note that although the blue reverse annealing trend looks to perform the best, however this trend is artificially inflated because 6 of the problems have e i = e 0 which has been demonstrated in Fig. 10 to yield ap s on the order of 10 −2 at s = .7 and t p = 400µs. We next visualize a histogram, as seen in Fig. 14, of all energies recorded from 1000 samples returned for a set of 94 FIG. 13. Theps as a function of n over a set of 100 problems each with 1000 samples. We compare reverse annealing (blue) with ei = e f , s = .7, and tp = 400µs to forward annealing (orange) with clique embedding, g = 0, and annealing time = 100 µs. We also compare the combination of forward annealing and reverse annealing where theps is chosen by problem (green). In this green trend, theps is calculated using the forward annealingp  problems where forward annealing did not find e 0 with n = 20. We compare forward annealing to reverse annealing where e i = e f . We observe even for problems where neither reverse annealing or forward annealing found e 0 , reverse annealing still on average finds a lower energy solution than forward annealing.

V. CONCLUSIONS
We have benchmarked quantum annealing using Markowitz portfolio selection to evaluate the effects of various controls on probability of success and chain breaks. We have explored a variety of quantum annealing controls including the embedding algorithm, the forward annealing time T , and the number spin reversal transforms g. When comparing clique embedding against CMR embedding, we found little difference in the estimated probability of successp s as both techniques yielded a sub-exponential decay forp s with exponents of −0.528 and −0.523, respectively. We did observe that CMR demonstrated a slightly higher probability of chain breaksp b , and we considered this a sufficient justification to use the clique embedding for studying the fully connected problems Markowitz portfolio selection problem.
When varying the forward annealing time T ∈ [1µs, 999µs], we found thatp b was slightly higher in the range T = [1µs, 5µs] while increasing the annealing time further yielded little to no improvement. For this reason, we chose to continue all future forward annealing experiments using T = 100µs where the exponential decay rate in p s was −0.528. When varying g = [0, 10], we found small improvements inp s between g = 0 and g = 2 where the exponential decay rate became −0.505 without much change from increasing the value of g further, and there was no consistent difference inp b .
We benchmarked reverse annealing controls with respect to the parameters e i , s, and t p . We began by observing the results inp s andp b at n = 20. We consistently observed thatp b was the same order of magnitude as with the forward annealing experiments andp b was consistently highest for s = 0.8. By setting e i = e 0 , we observed that thep s decreases exponential as s increased. By setting e i = e 1 , we observed that reverse annealing had ap s an order of magnitude higher than forward annealing. From these results, we conclude that when e i is close to the ground state, reverse annealing provided some advantage over forward annealing. Because in general the ground state won't be known for a problem, we developed a heuristic which sets e i = e f where we again observedp s to be an order of magnitude higher than using forward annealing alone.
We further evaluatedp s as a function of n to compare reverse annealing with e i = e f , s = 0.7, and t p = 400µs to forward annealing with clique embedding, T = 100µs, and g = 0 alone. In particular, we used thep (k) s of forward annealing for the 6 problem instances in which e i = e 0 and thep (k) s of reverse annealing for the 94 problems where e i = e 0 . We continued to observe reverse annealing demonstrate an order of magnitude increase in p s over forward annealing alone. Lastly, by creating a histogram which plots the lowest energies found across 1000 samples for the 94 problems where e i = e 0 , we found that reverse annealing(e i = e f ) on average finds lower energy solutions as compared to forward annealing.
In summary, the benchmarks presented here evaluate a variety of quantum annealing controls with respect to the baseline ground truth for portfolio selection. By comparing the observed influence of these controls on the performance of solution accuracy, we have developed insights into the best selections of controls for solving these problems with the highest accuracy which may help guide the future use of quantum annealing as a meta-heuristics for optimization.