title: Hybrid quantum-classical optimization for financial index tracking
authors: Samuel Fernández-Lorenzo Diego Porras Juan José García-Ripoll
text: 
  
  
  
  
  BBVA Client Solutions Research & Patents, Calle Sauceda 28, 28050 Madrid, Spain.
  
  samuel.fernandez.lorenzo.contractor@bbva.com
  
  
  Instituto de Física Fundamental, IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain.
  
  
  Instituto de Física Fundamental, IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain.
  
  
  
  
  In this work we develop tools to address combinatorial optimization problems with a cardinality constraint, in which only a subset of variables end up having nonzero values. 
  Firstly, we introduce a new heuristic pruning method that iteratively discards variables through a hybrid quantum-classical optimization step. Secondly, we analyse the use of soft constraints in the form of “chemical potentials" to control the number of non-zero variables. 
  We illustrate the power of both techniques using the problem of index tracking, which aims to mimicking the performance of a financial index with a balanced subset of assets. 
  We also compare the performance of different state-of-the-art quantum variational optimization algorithms in our pruning method. 
  
  
  
  
  Hybrid quantum-classical optimization with cardinality constraints and applications to finance
      Juan José García-Ripoll
      December 30, 2023
  ==============================================================================================
  
  
  
  
  § INTRODUCTION
  
  Many relevant problems in quantitative finance translate into daunting computational tasks, such as combinatorial optimization problems and Monte Carlo simulations <cit.>, suffering from lack of parallelization or slow convergence. The emergent field of quantum finance, rooted in quantum physics and quantum computing, develops new algorithms and formulations of financial problems, to address these problems and limitations <cit.>. Recent works in quantum finance have explored the design of optimal trading trajectories <cit.>, credit scoring <cit.>, portfolio optimization <cit.>, Monte Carlo pricing of derivatives <cit.>, risk analysis <cit.> or financial crisis forecast <cit.>, among others.
  
  Among these challenges, in this work we focus on the family of quadratic optimization with cardinality and linear constraints, which we write as
  
  
      _𝐰 T_err :=_𝐰[𝐰^TΣ𝐰-2𝐰^T𝐠+ϵ_0] ,
         ||𝐰||_0=∑_j=1^N |w_j|^0= d, 
         A𝐰≤𝐛.
  
  
  Here, the w_j∈(-∞,∞) are N weights of which only d<N are nonzero, ‖w‖_0=d, using 0-norm notation, and Σ being a positive semi-definite matrix. Equality and inequality linear constraints are condensed in the term A𝐰≤𝐛. 
  
  Many canonical financial problems can be cast into the form (<ref>), like portfolio selection in a mean-variance framework <cit.> or  index tracking <cit.>, where a benchmark index with N possible securities is approximated by a basket of d≪ N assets properly selected and weighted. Beyond finance, we may encounter problem (<ref>) in other disciplines such as machine learning—like in future selection of a multi-variate linear regression—or supply-chain optimization. In this work, we shall rely on a simplified version of the index tracking problem as a workbench.
  Problem (<ref>) poses an implicit mixed-integer optimization as the selection of d variables is already a combinatorial optimization whose state space grows very fast as N d. 
  The limit N, d ≫ 1, and N ≫ d, is specially relevant for many practical optimization tasks, including the index tracking problem. 
  In this limit, the number of possible solutions actually grows exponentially as a function of the number of total variables, showing a dependence ∼ N^d. Exploring all configurations by brute force quickly becomes unfeasible, like in the case of financial problems composed by hundreds or thousands of securities, thereby approximation methods <cit.> and heuristic approaches come into play <cit.>. In particular, when addressing this problem in Noisy Intermediate Scale Quantum (NISQ) devices <cit.>, a typical route is to discretize the weights and transform it into a Quadratic Unconstrained Binary Optimization (QUBO) problem with constraints. Unfortunately, this leads to an algebraic growth in the number of qubits which severely limits the size and interest of real instances that can be addressed.
  
  Our solution to problem (<ref>) is a hybrid quantum-classical algorithm that we call heuristic k-step pruning (k-PA). This method nests two interdependent optimization problems: a classical convex optimization stage to determine the relative weights of the assets, combined with a quadratic binary optimization that determines the relevant subset of d variables that optimizes (<ref>). This approach makes best use of small NISQ devices, maximizing the number of weights that are described by a single quantum register—namely, N can include a universe of assets with up to N securities. We demonstrate that the pruning method can be executed in NISQ quantum computers (e.g. IBM-Q), using variational methods with a very small number of gates but very good performance. However, the same method can be extended to work with quantum annealers (D-Wave's) or quantum-inspired optimizers. 
  
  We stress that our pruning algorithm uses a classical computer for the convex optimization step, 
  which is not computationally hard, whereas it uses a quantum computer to address the hard combinatorial problem of selecting the d variables that take non-zero values. 
  Thus, our algorithm “concentrates” the hardness of combinatorial optimization into the step that is solved by the quantum computer. 
  Our algorithm has a further advantage since, even though the original problem is stated in terms of continuous variables, the quantum computer is used precisely in the step where that problem is reduced to binary variables, which can be more efficiently stored and processed by qubits. 
  In principle, the variable selection step could also be solved by means of classical QUBO solvers. However, in applications beyond hundred or thousand variables, it is computationally challenging to find the global optimum, and here quantum computers may offer an advantage over classical methods.
  Since we expect that near or mid-term quantum devices will have severe limitations in size, the advantages offered by our scheme can be decisive for finding useful practical applications of those devices. 
  
  In Section <ref> we describe and analyze our first contribution, the pruning algorithm. In Section <ref> we discuss the introduction of constraints in optimization problems. We compare the use of hard constraints with an approach based on Lagrange multipliers or “chemical potentials”. Section <ref> introduces the variational ansätze and classical optimizers employed in our study, as well as the index tracking problem that we will use to benchmark our solutions. Our algorithms and existing techniques are then numerically compared in Section <ref> by using a real dataset of intra-day data of the PHLX Oil Service Sector (OSX) index, composed of 15 companies involved in the oil services sector (cf. <ref>). 
  Finally, Section  <ref> discusses the conclusions of our study.
  
  
  
  
  
  § ITERATIVE ELIMINATION OF VARIABLES
   
  
  
  
  
   §.§ One-step pruning
  
  
  The optimization problem from Eq. (<ref>) is a mixed-integer optimization that implicitly joins a computationally hard task—the selection of d variables within an exponentially large family of choices—, with a classically tractable, convex optimization problem—assigning the weights of the d selected variables {w_j_n}_n=1^d—. This suggest that we split the whole problem into two problems, which we later address by a quantum and a classical computer, in line with the current philosophy of hybrid classical-quantum computations.
  
  The first iteration of our idea is the single-step selection algorithm (1-SA). Let us define a vector 𝐱 of decision variables, which are nonzero x_i=1 only when the corresponding weight w_i is included in the final basket of variables. We write down the combinatorial optimization problem of selecting the d variables that will have non-zero weights
  
  
      _𝐱 S_err := _𝐱[𝐱^TΣ𝐱-2𝐱^T𝐠] 
         ||𝐱||_0= d.
  
  
  The solution of this is used to build a cost function with smaller matrices, Σ_r(𝐱) and g_r(𝐱), that only contain the rows and columns where x_i≠ 0. The reduced problem gives us a smaller vector of weights 𝐰_r∈ℝ^d by solving the following convex optimization problem
  
  
      _𝐰_r T_err^r(𝐱) :=_𝐰_r[𝐰_r^TΣ_r(𝐱)𝐰_r-2𝐰_r^T𝐠_r(𝐱)+ϵ_0 ] 
         A𝐰≤𝐛.
  
  
  This 1-SA formulation is problematic because the selection of variables does not consider their relative importance in the final weights. A better strategy is the single-step pruning algorithm (1-PA). In this reversed strategy, the cost function of the selection problem incorporates the weights, through the transformation x_i→ w_ix_i. More precisely, we first solve the convex optimization problem
  
  
      _𝐰 T_err^full :=_𝐰[𝐰^TΣ𝐰-2𝐰^T𝐠+ϵ_0]  
         A𝐰≤𝐛.
  
  
  After solving Eq. (<ref>), we obtain a vector of optimal weights, w_opt. This information is subsequently incorporated into a diagonal matrix D_ij( w_opt)= 
  w_opt,iδ_ij, with which we select the variables in a modified purely combinatorial optimization problem
  
  
      _𝐱 P_err(𝐰_opt) :=_𝐱[𝐱^T𝐃(𝐰_opt)Σ 𝐃(𝐰_opt)𝐱
      -2𝐱^T𝐃(𝐰_opt)𝐠] 
         ||𝐱||_0= d.
  
  
  Finally, we update the optimal weights for this combination through another convex optimization problem on a reduced subspace like in Eq. (<ref>).
  
  
  
   §.§ k-step pruning
  
  
  We can use 1-PA as the basis for an iterative pruning algorithm that constructs solutions for decreasing basket sizes. More precisely one would design a discrete schedule with k steps and decreasing universe sizes[One possibility is to decrease the size of the search spaces linearly, N_i=N-s× i, but other schedules are also feasible.] N_1=N > d_1=N_2 > …≥ N_k > d_k=d. For the k-th step, the 1-PA algorithm assigns a series of weights to all variables in the universe of size N_k, constructs the reduced matrices Σ_r and g_r and selects a basket of d_k nonzero variables that will move on to the next step.
  
  We can combine the iterative pruning with stochastic optimization methods, such as a quantum optimization algorithm, to solve each of the 1-PA variable selection steps. In this scenario, we have to fine tune not only the universe size, but also the number of repetitions of the stochastic method. We suggest to gradually increase the number of repetitions as the universe size shrinks, according to the schedule r→ r_0 +α r. Heuristically, in early stages where d∼ N, we are more likely to retain variables in the optimal portfolio, so it makes sense to spend less repetitions. In later stages, each repetition will be exponentially cheaper and more accurate, due to the decrease in universe size, and it pays off to use the same number of resources to get more accurate intermediate solutions. We refer to the combination of universe size and repetition schedule as our heuristic k-step pruning algorithm (k-PA). It is described in Algorithm <ref> as a pseudocode function to calculate a target basket of size d_target and weights w_j from an universe of size N in steps of size s.
  
  
  
  In Algorithm <ref>, r_0, α and s are hyperparameters that have to be chosen depending on the specific problem at hand. 
  Parameter s is particularly relevant, since it determines the number of pruning steps, as well as the dimension of the configuration space at each step. 
  The value of s should be large enough so that at each step of the pruning process,
  a subset of weights are selected that takes into account significant correlations with the rest of the weights. 
  However, a large value of s can also induce errors related to the reduction of the original weights, w_j, to binary values, x_j. 
  The optimal value of s will thus inevitably have to be calculated for particular problems by directly applying our method and comparing the quality of the solutions.
  
  
  
  
  § OPTIMIZATION WITH CONSTRAINTS
  
  
  
  The pruning algorithm consists of two stages: a classically simple step in which the N variables are assigned prospective weights by convex optimization, and a second one involving a quadratic binary optimization to select d nonzero variables (<ref>). We will now discuss two methods to treat this second stage.
  
  
  
   §.§ Hard-constraint formulation
   
  
  Quantum variational algorithms where originally devised for tackling unconstrained optimization problems. One alternative to circumvent this limitation is to include the cardinality constraints in Eq. (<ref>) into a Quadratic Unconstrained Binary Optimization (QUBO) cost function, as quadratic penalties.
  
  
      _𝐱 P_err^QUBO = _𝐱𝐱^T𝐃Σ 𝐃𝐱-2𝐱^T𝐃𝐠+P(1^T𝐱 - d)^2.
  
  
  Here, 1 is a vector of ones and P is a penalty weight—large enough to prevent unfeasible configurations, but small enough to allow tunneling between different baskets that satisfy the constraints. The optimization problem (<ref>) can be brought to a quadratic form P_err^QUBO∼𝐱^T Q 𝐱, including the linear terms[Note that for binary variables x_i^2 = x_i x_i = x_i and
    ∼𝐛^t 𝐱=∑_ijx_ix_j b_iδ_ij.]. The QUBO formulation with hard constraints allows us to control very well the basket size, but it creates very rough energy landscapes where classical optimizers such as COBYLA get easily trapped. For that reason, we introduce now a technique inspired in the notion of “chemical potential” from statistical and condensed matter physics.
  
  
  
   §.§ Soft-constraint formulation
   
  
  We can improve the convergence of the optimizer by relaxing our control over the basket size, changing the cost function to include Lagrange multipliers instead of quadratic penalties. The new regularization term acts as a chemical potential that favors different values of the constraint, depending on the regularization parameter λ,
  
      P_err^QUBO_soft∼𝐱^T Q 𝐱+λ||𝐱||_0.
  
  Although regularization is widely used in machine learning (such as lasso and ridge regularization), this particular form is not so common, because it leads to non-differentiable functions and NP-hard formulations. In our case the regularization amounts to introducing a multiple of the identity matrix 1
  
  
      _𝐱 P_err^QUBO_soft =_𝐱[𝐱^T𝐃(Σ+λ1) 𝐃𝐱-2𝐱^T𝐃𝐠],
  
  
  updating the diagonal elements of Q. Note that in the case of financial applications, λ can also be interpreted as a constant transaction cost for each security. The cardinality constraint is thus not explicitly imposed on Eq. 
  (<ref>), but the term λ is used to control the number of non-zero variables.
  
  A drawback of this method is that we cannot predict a priori the number of variables d(λ) selected by the regularization parameter λ. In practical applications, such as index tracking, this is not a problem, because the same optimization can be repeated multiple times, with minor variations in the data, but not the structure of the problem. Under such circumstances, one can find a monotonic relationship between the regularization parameter and the average number of variables selected (as exemplified in figure <ref>).
  
  
  
  
  § BENCHMARKS
   
  
  
  In the remaining of this work we will study the application of the pruning, with hard and soft constraints, to a model problem of index tracking that uses data from a real-world index. For the variable selection stages we will explore the use of both quantum and classical algorithms. In the following sections we describe the quantum variational algorithm, with the ansätze that we will apply, the classical algorithms that we compare with, and the specific instance of the index tracking problem which we use to define our numerical simulations.
  
  Applications of quantum algorithms to similar problems, like portfolio optimization have been published in the last years. 
  Quantum algorithms which offer a speed up in portfolio optimization can be found in <cit.>, although they are not specifically designed for the index tracking problem discussed here. Other algorithms for portfolio optimization using hybrid quantum-classical methods like <cit.> are closer to our own approach. The distinctive feature of our approach is that our scheme requires fewer qubits and thus it is more feasible with near term quantum devices.
  
  
  
   §.§ Index-tracking problem
  
  
  
  
  The goal of the index tracking problem is to replicate, over a given look-back window, a benchmark index by selecting a portfolio of d components out of a large universe N≫ d of possible assets. We expect that the smaller tracking portfolio will closely match the returns of the benchmark index, while at the same time reducing the transaction costs and simplifying the management process required to keep both the portfolio in-sync with the index (yet, excessive low values of d would result in poor tracking). We will now show that this goal can be mapped to a particular case of the generic optimization task (<ref>).
  
  Let us assume that we know the prices P_j(t) of the assets in our universe over an interval of time t∈[0,T], which we divide into periods with fixed size t_n=n×Δt. At any time, the composition of our portfolio is determined by a set c_j(t) of product units that are adapted during the tracking period. The value of the portfolio during this interval is
  
  
      P(t) = ∑_j=1^d c_j P_j(t),    t ∈ [0 , T].
  
  
  Within each sub-interval, we can define the overall return of the portfolio r_p(t_n) as the relative increment in portfolio value. It can be related to the returns of individual assets r_i(t_n)=P_i(t_n+Δt)/P(t_n)-1
  
  
      r_p(t_n)=P(t_n + Δ t)-P(t_n)/P(t_n) = ∑_j^N w_j  r_j(t_n) = 𝐰^T𝐫(t_n),
  
  
  through a vector of weights 𝐰∈ℝ^N characterizing the portfolio
  
  
      w_j=c_j(t) P_j(t)/∑_j^n c_j(t) P_j(t)∈[0,1].
  
  
  The benchmark index I(t) is known through a time series of returns, r_I(t_n)=I(t_n+Δt)/I(t_n)-1, taken by consecutive snapshots at periodic intervals. We will use this information, to construct the portfolio with d assets (and fixed weights) that best replicates the index' trajectory <cit.>. A reasonable assumption is that this optimal portfolio will perform similarly to the index over a short time horizon, t∈ [T, T + T_H]. A more sophisticated approach would involve a forecast of future returns of each index component, and to optimize our tracking portfolio over such forecast instead of optimizing over historical returns. 
  
  To implement the search of the optimal portfolio, we introduce a tracking error T_err that measures the distance between the sequence of returns of the index and our portfolio. We use the mean squared error between both time series, a quadratic form
  
  
      T_err=∑_n^T(r_p(t_n)-r_I(t_n))^2 =
      ∑_i,j = 1^N w_i Σ_i j w_j - 2 ∑_j=1^N w_j g_j + ϵ_0,
  
  
  with a matrix Σ_ij, a vector g_i and an offset ϵ_0, defined as
  
      Σ_i j = ∑_n=0^n_T r_i(t_n)r_j(t_n),   g_j = ∑_n=0^n_T r_j(t_n) r_I(t_n),  ϵ_0 = ∑_n=0^n_T r_I(t_n)^2.
  
  The index can also be expressed in terms of a vector of weights 𝐖∈ℝ^N valid for the specified time interval[Even if those weights are not explicitly provided or available in financial data, such vector can be computed by convex optimization of the universe of N index assets with no other restriction that demanding positive weights for all of them and proper normalization (budget constraint).], giving us explicit formulas for the components as
  
      g_i = ∑_j=1^N Σ_ijW_j  ϵ_0 = ∑_i,j=1^N W_iΣ_ijW_j.
  
  The search for the optimal portfolio results in the following mixed-integer optimization problem 
  
  
      _𝐰 T_err :=_𝐰[𝐰^TΣ𝐰-2𝐰^T𝐠+ϵ_0] 
         ||𝐰||_0=∑_j=1^N |w_j|^0= d,  
         ||𝐰||_1=∑_j=1^N |w_j|= 1, 
         w_j≥ 0   ∀ j .
  
  
  Here, the w_j∈[0,1] are N normalized weights with ‖w‖_1=1, of which only d are nonzero, ‖w‖_0=d, using the 1-norm and 0-norm notations. Eq. (<ref>)  is introduced to avoid short selling, and the integer nature of the problem is manifested when selecting a portfolio of d assets so that Eq. (<ref>) is respected. Notice that the index tracking problem here defined is an instance of the generic optimization problem (<ref>).
  
  Our formulation assumes that 𝐰 is a vector of continuous values, which is a good approximation commonly employed in practice. However, one could also tackle the problem as fully combinatorial recalling that equities are bought and sold in fixed lot sizes, thereby requiring a discretization in the values of w_j as discussed in <cit.>.
  
  Notice that the symmetric matrix Σ_ij is similar to a covariance matrix among assets,
  whereas the vector g_j accounts for covariance-like components with respect to the benchmark.
  Hence, the cost function (<ref>) can be interpreted similarly to a Markowitz's portfolio optimization: the mission of the quadratic term is to reduce the volatility associated to the uncertain evolution of the assets while the linear term favors those assets more correlated with the benchmark <cit.>.
  
  The problem, as stated, tends to favor assets with larger weights in the index. In realistic scenarios there may be further constraints to ensure a well-diversified portfolio, such as a capital concentration constraint of the form,
  
  
      𝐥≤ C𝐰≤𝐮,
  
  
  where 𝐥 and 𝐮 are vectors indicating the lower and upper bounds of certain linear combinations of investments gathered in the matrix C. Although this type of constraint is also linear and can be tackled by our heuristic pruning algorithm, we shall not consider them in this work. Instead, we focus on drawing the quantities Σ_ij and g_i from a data set that leads to a computationally demanding problem. More precisely, we will look for indices with assets that are strongly correlated, in a way that makes the quadratic terms are comparable to the linear terms. Otherwise, the problem can be solved in good approximation retaining just the linear contribution, which comes down to sorting assets by their contribution to the index and picking the top d from the list. For our numerical studies, we picked an index consisting of a cluster of 15 oil companies: the PHLX Oil Service Sector (OSX) index (N=15 consequently). The dataset gathers prices of all the instruments with interval Δ t of 20 minutes over several days. We take look-back windows corresponding to a single day, and average over several days. An example of the correlation heat map of these equities and their corresponding (normalized) component g_i for a random day (2019-06-17) is shown in figure <ref>.
  
  Let's clarify that such a small index tracking problem can be efficiently solved even by brute force. We use this problem simply as a workbench that allows us to classically simulate small quantum variational circuits and to demonstrate our algorithms. In a real scenario, a financial index composed of hundreds or thousands of securities along with capital concentration constraints would pose a real computationally hard task. Since we expect NISQ devices with hundreds or thousands over the next years, our work is intended to illustrate potential use cases of this technology. 
  
  
  
   §.§ Quantum variational algorithms
  
  
  We shall now address how to solve the integer optimization problem from Eq. (<ref>) using variational quantum circuits for combinatorial optimization <cit.>. Generally speaking, these methods represent the solution of the optimization problem as a quantum superposition over all possible configurations. The wavefunction is parameterized by the angles and phases of the quantum gates used to build it. These real parameters are processed by some classical optimizer, in a hybrid classical-quantum model that starts from rather uniform superpositions and aims at concentrating probability on the best arguments to our cost function. The use of quantum variational methods opens new questions, such as the shape and structure of the ansatz, the type of classical optimizers to use, or the role of quantum fluctuations and entanglement in exploring the complex and noisy space of solutions.
  
  We can write the QUBO cost function as an Ising model,
  
      H_err = ∑_ij J_ijσ^z_iσ^z_j + ∑_i σ^z_i h_i + E_0,
  
  by means of the mapping x_i →1/2(1-σ^z_i), from bits to Pauli Z operator σ^z_i. Our optimization problem translates to finding one of the ground states |s_1,s_2,…,s_N⟩ of H_err. We will approximate this as a variational method
  
      _θ⃗⟨Ψ(θ⃗)|H_err|Ψ(θ⃗)|⟩
  
  finding the best choice θ⃗ over a family of wavefunctions |Ψ(θ)⟩, constructed with gates that are parameterized by the angles θ_i. The optimization is done classically, constructing the wavefunction in the quantum computer, estimating the energy and estimating an update for the parameters. 
  
  We expect that, if the variational family is dense enough, |Ψ (θ⃗⟩) 
  will approximate the ground state with arbitrary accuracy, or at least the spin configuration that 
  minimizes H_ err can be detected with a
  significant probability after sampling  |Ψ (θ⃗⟩).
  In all the numerical experiments below we assume that variational wavefunctions are sampled by repeating experiments a number of times N_ meas = 100.
  
  
  
   §.§ Variational states
  
  
  Given the limitations of NISQ devices concerning circuit depth, a first approach is to restrict ourselves to hardware-efficient trial states <cit.>, variational wavefunctions created within the limitations of existing devices. More precisely, we will use the VQE Ry ansatz (henceforth just VQE) introduced in Ref. <cit.>. Formally, the wavefunction reads
  
  
      |Ψ_VQE (θ⃗) ⟩  =
       (∏_j=2^p  exp(-i∑_l θ_jlσ^y_i) U_ent) exp(-i∑_l θ_1lσ^y_i) |Ψ_0⟩.
  
  
  The ansatz starts from the zero state of the quantum register
  |Ψ_0⟩= |0⟩_1 ⋯|0⟩_N, and then subsequent layers of local rotations and entangling operations are applied. The variational parameters θ_jl determine the local rotation on the l-th qubit in the j-th layer. As entangling operations we use a sequence of control-Z gates between nearest-neighbor qubits, assuming the qubits form a 1D register.
  
  The Quantum Approximate Optimization Algorithm <cit.> (QAOA) is more hardware-demanding ansatz, that implements Trotterized version of quantum annealing. The variational parameters are the angles β_n of global x-rotations and of γ_n entangling operations implemented by the problem Hamiltonian,
  
  
      |Ψ_QAOA (γ⃗,β⃗)⟩  =
      ∏_n=1^p exp(- i β_n ∑_i σ^x_i ) exp(- i γ_n H_err)
      |Ψ_+⟩.
  
  
  The initial state is the uniform superposition of all register states | Ψ_+ ⟩ = ⊗_j=1^N |+⟩_j, with |+⟩_j = (1/√(2)) (|0 ⟩_j + |1⟩_j ). In the limit of very large number of layers, QAOA may converge to the dynamics of a quantum annealing process and therefore prepare the ground state with high fidelity. In practice, the number of layers p is smaller, but the ansatz captures the structure of the problem, which is why it is expected to perform better than the generic VQE method, but with an added cost of more gates to implement H_err.
  
  Another strategy is to force the quantum circuit to search over subspaces of the Hilbert spaces that respect the cardinality constraint  without the need of imposing an energy penalty.
  The idea is to start with an initial wavefunction with a well defined number of nonzero qubits and generate the variational ansatz with a sequence of quantum gates which do conserve spin-excitation number <cit.>.
  
  
      |Ψ_SWAP,d(θ⃗)⟩  =
      
            ∏_j=2^p [exp(-i ∑_i θ_jiσ^z_i ) U_√(SWAP)] exp(-i ∑_i θ_1iσ^z_i )
      |Ψ(d)⟩.
  
  
  The initial state |Ψ(d)⟩ has d qubits with value 1, equispaced over the quantum register. Local rotations around the z-axis are interleaved with U_√(SWAP) gates that implements a 50-50 beam-splitter, partially swapping the qubit excitations
  
  
      U_√(SWAP) = ∏_j=1^N exp( i π/4( σ^+_j σ_j+1^- + σ^+_j+1σ_j^-  )  ).
  
  
  
  
   §.§ Classical Optimization
  
  
  As for the classical optimization part, derivative-free optimizers like Constrained Optimization BY Linear Approximation (COBYLA) have proven to be more effective and robust in this task than gradient-based optimizers like gradient descent <cit.>. Nevertheless, all these methods tend to easily get stuck in local minima when the parameter landscape is highly nonconvex. 
  We checked that this problem is specially acute for the case of low-depth QAOA, in which all the energy structure plus the energy barriers penalizing unfeasible combinations get "compressed" into a few parameters. 
  An example of the poor performance of COBYLA when working with QAOA can be found in our Fig. <ref>a below
  For that reason, we chose the implementation of a global derivative-free optimizer based on dual annealing from scipy in Python, which combines a generalization of simulated annealing coupled to a local search. Unlike COBYLA, this method finalizes when the maximum number of iterations is reached; this is limited by the parameter maxiter, which regulates the maximum number of global search iterations.
  
  When using COBYLA, the parameter tol marks a convergence criteria; we found that tol=0.01 was enough to ensure convergence and avoid unnecessary iterations at the same time. Yet, the method has to be limited by a total number of iterations, which was set to 2000 because our numerical experiments normally converge within that limit. To attempt a fair comparison with dual annealing, we should specify a maxiter parameter provided by this method from which the quality of results plateaus. Finally, we have to establish some bounds where the global optimizer is going to be searching: we set the intervals [0,2π] and [0,π] for γ⃗ and β⃗ respectively. No further fine tuning was considered apart from the already mentioned.
  
  
  
  
  § NUMERICAL RESULTS
   
  
  
  
  
   §.§ Utility of heuristic pruning algorithm
  
  
  
  
  To test our heuristic pruning algorithm <ref>, we have first compared the global minima provided by the 1-SA and 1-PA algorithms against the brute search solution T_err^opt of the original index tracking problem defined in <ref>. Figures <ref>a-b compare the value of cost functions of each algorithm with the value of the true optimal solution, for all ratios of tracker size d vs universe size N and multiple dates of our dataset. Points located on the graph diagonal represent scenarios where the heuristic algorithm is able to capture the true optimal solution of the corresponding index tracking problem. The concentration of points along the diagonal is an evidence of the quality of the 1-PA method, which offers a much better heuristic than the 1-SA. We quantify the difference using the Pearson's correlation coefficient, which is 0.92 and 0.68 for 1-PA and 1-SA, respectively. Another way to compare both methods is to study the probability that an algorithm finds a portfolio with relative error with respect to the true optimal one Δ=(T_err-T_err^opt)/T_err^opt. As shown in figure <ref>c, 1-PA achieves an error Δ≤ 20% in 62.5% of the runs, while 1-SA only does it with 17.7% probability. Our analysis indicates that the heuristic pruning is able to incorporate the optimal or close to optimal portfolio for any given size d. 
  
  
  
  
  
  We have studied also the performance of the k-PA method, in combination with classical simulated annealing, a well-known classical stochastic algorithm. As publicly available reference, we used the dimod library <cit.> for QUBO problems, without fine tuning, but solving the problem with the hard-constraint penalty described in section <ref>. For a fair comparison, the experiment is designed to arrive to the same portfolio size (d=5) with the same number of total overall repetitions (120) through 1-PA (r_0=120,α=0), 2-PA  (r_0=120,α=4), and 3-PA (r_0=20,α=1) respectively.
  
  The series of figures <ref> illustrates an example of the 3-PA sequence for a random day (2016-06-17), with a final relative error of Δ=1.3%. 
  The bar plot in figure <ref>a displays a comparison of the median of the relative error Δ for the multi-step pruning algorithm.
  We observe that the error Δ decreases as we increase the number of steps, k. 
  This effect can be explained by the fact that a larger k leads to a smaller value s of neglected variables at each step, thus diminishing the error caused by the conversion from continuously valued weights, w_j, to binary variables, x_j.
  The boxplot <ref>b compares the time spent in each algorithm. 
  We observe the somehow counter-intuitive effect that a pruning algorithm with more steps takes a longer time to be executed. 
  However, note that a larger value of k in k-PA
  algorithms involves that a lower number of variables, s, has to be eliminated at each step. Thus each individual step is defined over a smaller configuration space, which explains the overall speedup effect.
  To summarize, our numerical experiment indicates that the k-PA is not only faster as we increase the number of steps, but also more accurate. 
  
  The resulting algorithm offers an alternative appealing heuristic with respect to others local-search algorithms studied in the literature, like the threshold accepting method <cit.>, because the combinatorial optimization stage can be tackled by means of several global-search strategies, including quantum optimization heuristics like the ones described in section <ref>.
  
  
  
  
  
  
   §.§ Optimization with hard constraints
  
  
  
  
  In solving the combinatorial optimization (<ref>), we have compared numerically the performance of  VQE (COBYLA), QAOA (COBYLA), SWAP (COBYLA) and QAOA (dual annealing) using the following procedure. We selected a set of random dates from the whole dataset. For each date we ran each algorithm with a random initial point for every portfolio size (in the case of QAOA the initial point lies within the region previously indicated). 
  Once the classical optimization converged, we sampled the variational wavefunction, simulating the outcomes of 100 measurements. 
  Out of the resulting binary numbers, we kept the configuration 𝐱 with the lowest energy that satisfies the cardinality constraint ‖𝐱‖_0=d. 
  As a classical reference, we also calculated the best result among 100 repetitions of simulated annealing using the dimod library.
  
  We quantify the performance of the stochastic methods, quantum and classical, using two metrics. Figure <ref>a compares the relative error obtained with different circuits and optimization methods, taking as basis the error of the optimal solution obtained through diagonalization,
  Δ=P_err-P^opt_err/P^opt_err. Figure <ref>b focuses on the number of evaluations for each quantum ansatz and optimization method. In all cases, we explore portfolio sizes in the interval d = [5,10], representing between 1/3 and 2/3 of the universe of equities with N = 15 (reasonable for the multi-step pruning algorithm). 
  Finally, to account for statistical fluctuations, we apply the non-parametric 
  Wilcoxon rank-sum test, which tests the null hypothesis that two sets of measurements are drawn from the same distribution—the alternative hypothesis is that values in one sample are more likely to be larger than the values in the other one. Let us point out that the dual annealing method depends on a  maxiter
  (maximum number of iterations) parameter. 
  To attempt a fair comparison with COBYLA, we numerically found a value from which the quality of results plateaus. 
  In view of figure <ref>, that shows the evolution of the relative error for QAOA with one layer and d = 5, we conclude
  that the method rapidly plateaus, being maxiter=10 a reasonable choice.
  
  
  
  
  
  We can offer some conclusions from the results 
  shown in figures <ref> and <ref>.
  First, the SWAP ansatz offers the lowest average relative error (⟨Δ|=⟩4%), followed closely by QAOA with dual annealing (6%). This is confirmed by the rank-sum test with a p-value p=0.0072 for a significance level α=0.05. VQE, in comparison, offers poor performance (31%). An additional VQE layer improves the quality, but requires twice the number of function evaluations. We can see that COBYLA quickly gets trapped in local minima for QAOA, leading to worse relative error (16%). 
  Finally, an advantage of QAOA over SWAP stems from the fact that it requires less function evaluations to achieve similar relative errors. To summarize, our main conclusion is that, on average, the SWAP ansatz leads to the smallest errors, however, QAOA with dual annealing requires much fewer function evaluations, while still showing reasonably small errors.
  
  We can also look into how the quality of results vary as we increase the ratio d/N. Figure <ref> displays swarm plots for each method whereby we can graphically compare the dispersion of results, their average and the probability of finding the global optimum.
  In the swarm plots, the accuracy of the optimization methods 
  can be visualized as the concentration of solutions exactly on Δ = 0, which indicates optimal solutions.
  All methods are consistent in showing a betterment of all these metrics as the portfolio grows. Interestingly, the probability of finding the optimal portfolio for SWAP and QAOA is better than simulated annealing for the hardest instance (d=5); the rank-sum confirms this shift in the distribution with p-values p=8.3× 10^-5 and p=2.7× 10^-3 respectively. 
  
  In comparison, the test does not find statistical difference for portfolios with a number of assets d > 7.
  
  
  
  
   §.§ Optimization with soft constraints
  
  
  In Sect. <ref> we introduced the idea of a Lagrange multiplier or chemical potential λ to enforce a desired value of the soft constraint d(λ). This assumed the monoticity of the function ⟨ d(λ)⟩ and the existence of an optimal λ^* satisfying ⟨ d(λ^*)⟩≈ d_target for each constraint—in our case, the portfolio size. Therefore, to generate portfolios of a target size d, we must tune[Remarkably, the values seem to be independent of the ansatz, perhaps suggesting that we could even train a supervised machine learning model that outputs a suitable λ given a certain matrix Q so as to automatize the parameter tuning.] λ so that the average portfolio produced by our quantum state coincides with d as much as possible. Fig. <ref> indicates that such tuning only needs to happen at the beginning of a series of optimizations over different time windows, because we have found that indices are stable enough to develop a monotonic relation between the regularization parameter and the average portfolio size. For a given λ, the average size remains fairly stable for different random initial points and different dates. 
  
  Small deviations from the optimal λ^* will not cause a failure to recover the right portfolios—they may affect the quality of the results and the likelihood of finding low energy configurations with the right size d—, and one might also design search strategies by which the value of λ may be dynamically adjusted.
  
  
  
  