Dataset
The data are collected from Yahoo!@finance29 using yfinance30, an open-source tool that uses Yahoos publicly available APIs. This tool, according to its creator, is intended for research and educational purposes.
To explore the efficiency of the proposed approach, small-sized examples are considered by extracting at most (N=4) different assets: (text {Apple}), (text {IBM}), (text {Netflix}) and (text {Tesla}). These are representative global assets with interesting dynamics influenced by financial and social events. For each asset i, with (1le i le N), the temporal range between 2011/12/23 and 2022/10/21 is considered. For each day t in this range ((0le t le T)), the performance of an asset is well represented by its closing price (p^t_i). A sub-interval of dates considered is shown in Table1. Additional experiments, performed on different dataset and falling within the same time interval considered here, are available in the supplementary information.
The first information extracted from this data set consists in the list P of current prices (P_i) of the considered assets.
$$begin{aligned} P_{i}=p^T_i. end{aligned}$$
(1)
Moreover, for each asset, the return (r^t_i) between the days (t-1) and t can be calculated:
$$begin{aligned} r^t_i=frac{p^t_i-p^{t-1}_i}{p^{t-1}_i} end{aligned}$$
(2)
These returns, calculated for days when the initial and the end prices are known, cannot be used for inference. Instead, it is convenient to define the expected return of an asset as an educated guess of its future performance. Assuming a normal distribution of the returns, the average of their values at each time t on the set of historical observations is a good estimator of the expected return. Therefore, given the entire historical data set, the expected return of each asset (mu _i) is calculated by:
$$begin{aligned} mu _i=E[r_i]=frac{1}{T}sum _{t=1}^{T}r^t_i. end{aligned}$$
(3)
Following the same principle, the variance of each asset return and the covariance between returns of different assets over the historical series can be calculated as follows:
$$begin{aligned}&sigma ^{2}_{i}=E[(r_{i}-mu _{i})^2]=frac{1}{T-1}sum _{t=1}^{T}(r^{t}_{i}-mu _{i})^{2}, \&sigma _{ij}=E[(r_{i}-mu _{i})(r_{j}-mu _{j})]=frac{1}{T-1}sum _{t=1}^{T}((r^{t}_{i}-mu _{i})(r^{t}_{j}-mu _{j})) nonumber . end{aligned}$$
(4)
The traditional theory of PO was initially formulated by Markowitz1. There are multiple possible formulations of PO, all embodying different degrees of approximation of the real-life problem. This work deals with Multi-Objective Portfolio optimization: this approach tries to simultaneously maximize the return and minimize the risk while investing the available budget. Even if other formulations include more objectives, the aim is still the solution of a constrained quadratic optimization problem; therefore, the formulation considered here is general enough to test the performances of the proposed approach.
A portfolio is defined as the set of investments (x_{i}) (measured as a fraction of the budget or number of asset units) allocated for each ith asset of the market. Therefore, the portfolio consists of a vector of real or integer numbers with dimensions equal to the number of assets considered. An optimal strategy for portfolio allocations aims to achieve the maximum portfolio return (mu ^{text {T}} x) while minimizing risk, defined as the portfolio variance (x^{text {T}}Sigma x) (whose square root is the portfolio volatility), where (mu ) is the vector of mean asset returns for each asset i calculated by (3), (Sigma ) is the covariance matrix calculated by (4), and x is the vector of investments measured as fractions of budget. Hence, the task of finding the optimal portfolio aims at finding the x vector that maximizes the following objective function:
$$begin{aligned} {{mathscr {L}}}(x): mu ^{text {T}} x - qx^{text {T}}Sigma x, end{aligned}$$
(5)
where the risk aversion parameter q expresses the propensity to risk of the investor (a trade-off weight between the risk and the return).
In a realistic scenario, the available budget B is fixed. Therefore, the constraint that the sum of (x_i) equals 1 must hold. Moreover, if only buying is allowed, each (x_ige 0), this constraint does not hold if either buying or selling is possible. As a consequence, in the general case, the problem can be stated as follows:
$$begin{aligned}&underset{x}{max }{mathscr {L}}(x): underset{x}{max }(mu ^{text {T}} x - qx^{text {T}}Sigma x),\&text {s.t.} quad sum ^{N}_{i=1}x_i=1 nonumber end{aligned}$$
(6)
However, if x is a possible solution to the problem with continuous variables, each product (x_iB) must be an integer multiple of the corresponding price (P_i) calculated by (1) since an integer number of units of each asset can be exchanged. Therefore, only a subset of the possible solutions corresponding to integer units is acceptable, and the problem is better stated as follows:
$$ begin{gathered} mathop {max }limits_{n} {mathcal{L}}(n):mathop {max }limits_{n} (mu ^{{prime {text{T}}}} n - qn^{{text{T}}} Sigma ^{prime } n), hfill \ {text{s}}{text{.t}}{text{.}}quad P^{{prime {text{T}}}} n = 1 hfill \ end{gathered} $$
(7)
where n is the vector of (n_i) integer units of each asset, while (P'=P/B), (mu '=P'circ mu ) and (Sigma '=(P'circ Sigma )^{text {T}}circ P') are appropriate transformations of (mu ) and (Sigma ). The latter formulation (7) is an integer constrained quadratic optimization problem.
Possible solutions to the problem (6) are those satisfying the constraint. Among them, some correspond to possible solutions to problem (7). The collection of possible solutions corresponding to portfolios with maximum return for any risk is called Markowitz efficient frontier. The solution of the constrained quadratic optimization problem lies on the efficient frontier, and the distance from minimum risk depends on q.
The general problem, if regarded in terms of continuous variables, can be solved exactly by Lagrange multipliers in case of equality constraints, or by KarushKuhnTucker conditions, which generalize the method of Lagrange multipliers to include inequality constraint31, as the covariance matrix is positive semi-definite32. Optimizing a quadratic function subject to linear constraints leads to a linear system of equations, solvable by Cholesky decomposition33 of the symmetrical covariance matrix. The exact solution involves the computation of the inverse of an (N times N) matrix, where N is the number of assets, thus requiring about (O(N^3)) floating-point operations34.
As long as integer or binary variables are considered, the problem turns into combinatorial optimization. The computational complexity is known to be high since the optimization problem is NP-hard35,36, while the decision version is NP-complete37. Indeed, a search approach should find the optimal one among possible solutions whose number increases exponentially with the number of assets (e.g., for b binary variables, (2^b) possible solutions, while for N integer variables ranging from 0 to (n_{max}), ({left( n_{max}+1right) }^N) possible solutions).
In practice, various methods are currently employed, either based on geometric assumptions, such as the branch-and-bound method2,3, or rather heuristic algorithms4,5,6, such as Particle Swarms, Genetic Algorithms, and Simulated Annealing. These have some limitations but allow to obtain approximate solutions. However, in all cases, the exact or approximate solution is feasible only for a few hundreds of assets on current classical computers.
Using quantum mechanical effects, like interference and entanglement, quantum computers can perform computational operations within the Bounded-error Quantum Polynomial (BQP) class of complexity, which is the quantum analogue of the Bounded-error polynomial probabilistic (BPP) class. Even if there is no NP problem for which there is a provable quantum/classical separation, it is widely believed that BQP (not subset ) BPP, hence when considering time complexity, quantum computers are more powerful than classical computers. More generally, it is conjectured that P is a subset of BQP. Therefore, while all problems that can be efficiently solved classically, are efficiently solvable by quantum computers as well, some problems exist that are considered intractable, until nowadays, by classical computers in polynomial space, and can they be solved with quantum machines. These facts are still matter of investigation but there are good reasons to believe that there are problems solvable by QC more efficiently than classical computers, thus QC will have a disruptive potential over some hard problems38, among which constrained quadratic optimization problems, including PO.
The branch-and-bound method3,7 is used in this work as a classical benchmark to compare the results of the proposed approach. It is based on the Lagrangian dual relaxation and continuous relaxation for discrete multi-factor portfolio selection model, which leads to an integer quadratic programming problem. The separable structure of the model is investigated by using Lagrangian relaxation and dual search. This algorithm is capable of solving portfolio problems with up to 120 assets.
Specifically, the library CPLEX freely available on Python provides a robust implementation of the aforementioned classical solving scheme.
As formulated in Eq. (7), the PO problem lies within the class of quadratic optimization problems. To be quantum-native, it has to be converted into a Quadratic Unconstrained Binary Optimization (QUBO) problem, i.e., the target vector to be found has to be expressed as a vector of zeros and ones, and constraints have to be avoided.
Therefore, the binary conversion matrix C is constructed with a number of binarizing elements (d_i) for each asset i depending on the price (P_i). Hence
$$begin{aligned} n^{max}_{i}=Intleft( frac{B}{P_i}right) , end{aligned}$$
(8)
where the operation Int stands for the integer part, and
$$begin{aligned} d_i=Intleft( log _{2}{n^{max}_i}right) , end{aligned}$$
(9)
such that
$$begin{aligned} n_{i}=sum _{j=0}^{d_i}2^{j}b_{i,j}. end{aligned}$$
(10)
In this way, the overall dimension of the binarized target vector, (b=left[ b_{1,0},dots ,b_{1,d_1},dots ,b_{N,0},dots ,b_{N,d_N}right] ), is (text {dim}(b) =sum _{i=1}^{N}left( d_i+1right) ), which is lower than that used in implementation available in Qiskit15. Conveniently, the encoding matrix C is defined as follows:
$$begin{aligned} C= begin{pmatrix} 2^{0} &{} dots &{} 2^{d_1} &{} 0 &{} dots &{} 0 &{} dots &{} 0 &{} dots &{} 0 \ 0 &{} dots &{} 0 &{} 2^{0} &{} dots &{} 2^{d_2} &{} dots &{} 0 &{} dots &{} 0 \ vdots &{} ddots &{} vdots &{}vdots &{} ddots &{}vdots &{} ddots &{}vdots &{}ddots &{}vdots \ 0 &{} dots &{} 0 &{} 0 &{} dots &{} 0 &{} dots &{} 2^{0} &{} dots &{} 2^{d_N} end{pmatrix}, end{aligned}$$
(11)
and thus, the conversion can be written in short notation as (n = Cb). It is possible to redefine the problem (7), in terms of the binary vector b, applying the encoding matrix by (mu ''=C^{text {T}}mu '), (Sigma ''=C^{text {T}}Sigma 'C) and (P''=C^{text {T}}P'):
$$ begin{gathered} mathop {max }limits_{b} {mathcal{L}}(b):mathop {max }limits_{b} left( {mu ^{{prime prime {text{T}}}} b - qb^{{text{T}}} Sigma ^{{prime prime }} b} right), hfill \ {text{s}}{text{.t}}{text{.}}quad P^{{prime prime {text{T}}}} b = 1 hfill \ quad quad b_{i} in { 0,1} quad forall i in left[ {1, ldots ,dim(b)} right]. hfill \ end{gathered} $$
(12)
The problem (12) falls into the wide set of binary quadratic optimization problems, with a constraint, given by the total budget. In this form, the problem cannot be cast directly into a suitable set of quantum operators that run on quantum hardware: the constraint, in particular, is troublesome, as it poses a hard limitation on the sector of Hilbert space that needs to be explored by the algorithm, to find a solution. It is thus necessary to convert the problem into a QUBO (Quadratic Unconstrained Binary Optimization) by transforming the constraint into a penalty term in the objective function. Each kind of constraint can be converted into a specific penalty term39, and the one considered in (12), which is equality, linear in the target variable, maps into (lambda (P''^{text {T}} b-1)^{2}), such that (12) can be written in terms of the following QUBO problem:
$$ mathop {max }limits_{b} {mathcal{L}}(b):mathop {max }limits_{b} left( {mu ^{{prime prime {text{T}}}} b - qb^{{text{T}}} Sigma ^{{prime prime }} b - lambda (P^{{prime prime {text{T}}}} b - 1)^{2} } right). $$
(13)
The penalty coefficient (lambda ) is a key hyperparameter to state the problem as the QUBO of the objective function (13).
There is a strong connection, technically an isomorphism, between the QUBO and the Ising Hamiltonian40: Ising Hamiltonian was originally constructed to understand the microscopic behavior of magnetic materials, particularly to grasp the condition that leads to a phase transition. However, its relative simplicity and natural mapping into QUBO have made the Ising model a fundamental benchmark well beyond the field of quantum physics. To convert (13) into an Ising, it is convenient to expand it in its components:
$$ {mathcal{L}}(b):sumlimits_{i} {mu _{i}^{prime } b_{i} } - qsumlimits_{{i,j}} {Sigma _{{i,j}}^{prime } } b_{i} b_{j} - lambda left( {sumlimits_{i} {P_{i}^{prime } b_{i} - 1} } right)^{2} , $$
(14)
where (mu ''_{i}, Sigma ''_{i,j}, P''_{i} ), are the components of the transformed return, covariance, and price, respectively, and (i,jin left[ 1,dim(b)right] ). Since the Ising represents spin variables (s_{i}), which have values ({-1,1}), the transformation (b_{i}rightarrow frac{1+s_{i}}{2}) is applied and coefficients are re-arranged, to obtain the Ising objective function to minimize:
$$begin{aligned}&underset{s}{min }{mathscr {L}}(s): underset{s}{min }left( sum _{i}h_{i}s_{i}+ sum _{i,j} J_{i,j}s_{i}s_{j}+lambda (sum _{i}pi _{i}s_{i}-beta )^{2}right) ,\&text {s.t.} quad s_{i,j}in {-1,1} quad forall i nonumber , end{aligned}$$
(15)
with (J_{i,j}) being the coupling term between two spin variables. It is now straightforward to obtain the corresponding quantum Hamiltonian, whose eigenvector corresponding to the minimum eigenvalue corresponds to the solution: in fact, the eigenvalues of the Pauli operators Z are (pm 1). Thus they are suitable for describing the classical spin variables (s_{i}). Furthermore, the two-body interaction term can be modeled with the tensor product between two Pauli operators, i.e., (Z_{i}otimes Z_{j}). The quantum Ising Hamiltonian reads:
$$begin{aligned} H= sum _{i}h_{i}Z_{i} + sum _{i,j} J_{i,j}Z_{i}otimes Z_{j}+lambda (sum _{i}pi _{i}Z_{i}-beta )^{2}. end{aligned}$$
(16)
With the procedure described above, the integer quadratic optimization problem of a portfolio allocation with budget constraints is expressed first as a binary problem via the binary encoding, then it is translated into a QUBO, transforming the constraints into a penalty term by the chosen penalty coefficient, and finally into a quantum Hamiltonian written in term of Pauli gates. Hence, the PO problem (7) is now formulated as the search of the ground state, i.e., the minimum energy eigenstate, of the Hamiltonian (16). Therefore, it is possible to use the VQE, employing real quantum hardware, and iteratively approximate such a state, as described in the following section, which corresponds to the optimal portfolio.
The VQE is a hybrid quantum-classical algorithm41, which is based on the variational principle: it consists in the estimation of the upper bound of the lowest possible eigenvalue of a given observable with respect to a parameterized wave-function (ansatz). Specifically, given a Hamiltonian H representing the observable, and a parameterized wave-function ({|{psi (theta )}rangle }), the ground state (E_{0}) is the minimum energy eigenstate associated s
$$begin{aligned} E_{0}le frac{{langle {psi (theta }|}H{|{psi (theta )}rangle }}{{langle {psi (theta )|psi (theta )}rangle }}, quad forall quad theta . end{aligned}$$
(17)
Hence, the task of the VQE is finding the optimal set of parameters, such that the energy associated with the state is nearly indistinguishable from its ground state, i.e., finding the set of parameters (theta ), corresponding to energy (E_{min}), for which (|E_{min}-E_{0}| $$begin{aligned} E_{text {min}}=underset{theta }{min }{langle {textbf{0}}|}U^{dagger }(theta )HU(theta ){|{textbf{0}}rangle }. end{aligned}$$ (18) where (U(theta )) is the parametrized unitary operator that gives the ansatz wave-function when applied on the initial state, (E_{min}) is the energy associated with the parametrized ansatz. The Hamiltonian H, defined for the specific problem, and in this case corresponding to (16), can be written in a specific operator basis that makes it naturally measurable on a quantum computer: this choice depends on the architecture considered. In this work, given the extensive use of the IBM quantum experience42, it is convenient to map the Hamiltonian into spin operators base. This base is formed by the tensor product of Pauli strings: (P_{l}in {I,X,Y,Z}^{otimes N}). In this base the Hamiltonian can always be written in the general form, (H=sum _{l}^{D}c_{l}P_{l}), where D is the number of Pauli strings that define the Hamiltonian and (c_{l}) is a suitable set of weights. It follows that the VQE in Eq. (18) can be written as: $$begin{aligned} E_{text {min}}=underset{theta }{min }sum _{l}^{D}c_{l}{langle {textbf{0}}|}U^{dagger }(theta )P_{l} U(theta ){|{textbf{0}}rangle }. end{aligned}$$ (19) Each term in Eq. (19) corresponds to the expectation value of the string (P_{l}) and is computed on quantum hardware (or a simulator). The summation and the optimization of the parameters are computed on a classical computer, choosing an ad-hoc optimizer. The eigenvector corresponding to the ground state corresponds to the solution of the problem(13), thus to the optimal portfolio. Schematic of the VQE algorithm. The ansatz wave-function ({|{psi (theta )}rangle })) is initialized with random parameters and encoded in a given set of quantum gates. The PO problem is translated into an Ising Hamiltonian and encoded into a set of Pauli gates. The collection of output measurement allows the reconstruction of the expectation value of the Hamiltonian H, which is the energy that needs to be minimized. A classical optimization algorithm provides an update rule for the parameters of the wave-function, which ideally moves iteratively towards the ground state of the problem, thus providing an estimation of the corresponding eigenstate. This corresponds to the solution of the original PO problem. In light of what is stated above, the complete VQE estimation process can be decomposed in a series of steps, as depicted in Fig. 1. First, it is necessary to prepare a trial wave-function (ansatz) on which the expectation value needs to be evaluated and realized via a parameterized quantum circuit. Then, it is necessary to define the Hamiltonian (16), whose ground state is the solution to the problem to be addressed, and convert it into the Pauli basis so that the observable can be measured on the quantum computer. Finally, the parameters are trained using a classical optimizer. This hybrid system ideally converges to a form that produces a state compatible with the ground state of the Hamiltonian. This procedure includes two hyperparameters that have to be settled, i.e., the type of ansatz and the optimizer. When defining the ansatz, two main features have to be taken into account: its expressivity, i.e., the set of states that can be spanned by the ansatz itself, and the trainability, i.e., the ability of the ansatz to be optimized efficiently with available techniques. It is worth pointing out the problem of the barren plateau43, related to the possibility of vanishing gradients when the cost function gradients converge to zero exponentially, as a function of the specific characteristic of the problem to be solved. The barren plateau depends on the number of qubits, the high expressivity of the ansatz wave-function, the degree of entanglement, and the quantum noise44. There are several methods to avoid or mitigate the effect of the barren plateau, especially in the context of VQE, most of which consist in finding a trade-off between the expressivity of the ansatz and its trainability and reducing the effective size of the Hilbert space of the problem formulation45. The following ansatzes are available in Qiskit and are analyzed in this work: Two Local ansatz, where qubits are coupled in pairs, the Real Amplitude ansatz, which assumes real-valued amplitude for each base element of the wave-function, and the Pauli Two ansatz, used mainly in quantum machine learning for the mitigation of barren plateu46. Although other ansatzes are provided in Qiskit, they are generally unsuitable for a PO problem. For instance, the Excitation preserving ansatz preserves the ratio between basis vector components, hence does not allow, in principle, any weight imbalance in the output distribution while moving towards the solution of the problem. For all the ansatzes considered, the convergence of four different possible assumptions on the entanglement structure of the wave-function is checked, namely the full entanglement, the linear entanglement, the circular and the pairwise entanglement. The former modifies the ansatz such that any qubit is entangled with all the others pairwisely. In the linear case, the entanglement is built between consecutive pairs of qubits. The circular case is equivalent to the linear entanglement but with an additional entanglement layer connecting the first and the last qubit before the linear sector. Finally, in the pairwise entanglement construction, in one layer, the ith qubit is entangled with qubit (i+1) for all even i, and in a second layer, qubit i is entangled with qubit (i+1), for odd values of i. Once the ansatz is defined, its parameters must be optimized classically until convergence is reached. The choice of the optimizer is crucial because it impacts the number of measurements that are necessary to complete the optimization cycle since, when properly chosen, it can mitigate the barren plateau problem and minimize the number of iterations required to reach convergence. In this work, dealing with the PO problem, different optimizers are tested to select which one fulfills its task faster, among those available on Qiskit, i.e., Cobyla, SPSA, and NFT47. The experimental results presented in this work are obtained on real quantum hardware, specifically using the platforms provided by IBM superconducting quantum computers. These quantum machines belong to the class of NISQ devices, which stands for Noisy Intermediate Scale Quantum devices, i.e., a class of hardware with a limited number of qubits and where noise is not suppressed. Noise, in quantum computers, comes from various sources: decoherence, gate fidelities, and measurement calibration. Decoherence is the process that most quantum mechanical systems undergo when interacting with an external environment48. It causes the loss of virtually all the quantum properties of the qubits, which then collapse into classical bits. Gate fidelities measure the ability to implement the desired quantum gates physically: in the IBM superconducting qubits hardware, these are constructed via pulses, which are shaped and designed to control the superconductors. Given the limited ability to strictly control these pulses, a perfect gate implementation is highly non-trivial and subject to imperfections. Last, measurement errors are caused by the limits of the measurement apparatus, improper calibration, and imperfect readout techniques. Hence, NISQ devices do not always provide reliable results due to the lack of fault tolerance. However, they provide a good benchmark for testing the possibilities of quantum computing. Furthermore, ongoing research is on the possibility of using NISQ in practical applications, such as machine learning and optimization problems. In this work, both simulators and real quantum computers are used. Even though error mitigation techniques49 can be applied, the main goal of this paper is to test the performances of the quantum computers on a QUBO problem, such as PO, without error mitigation, with the binary encoding strategies and the budget constraints as described in the previous sections. Therefore, in all computations, there is no error mitigation, aiming to build an indirect but comprehensive analysis of the hardware limitations and to improve the quality of the results offered by a proper selection of the hyperparameters. This will provide a solid benchmark for the following experimental stages, which will be enabled in the coming years by large and nearly fault-tolerant quantum computers. Hence, the experiments run on simulators (without noise) are also executed by adding noise mimicking real hardware: this operation can be readily implemented on Qiskit by inserting a noise model containing the decoherence parameters and the gate error rate from real quantum hardware. Moreover, experiments are run on IBM NISQ devices with up to 25 qubits. Specifically, a substantial subset of the available quantum computers in the IBM quantum experience was employed: IBM Guadalupe, Toronto, Geneva, Cairo, Auckland, Montreal, Mumbai, Kolkata, and Hanoi. These machines have either 16 or 27 qubits, but they have different quantum volumes (QV) and Circuit Layer Operations Per Second (CLOPS). QV and CLOPS are useful metrics to define the performances of a quantum computation pipeline50. Generally, a bigger QV means that the hardware can sustain deeper circuits with a relatively small price on the performance. At the same time, the CLOPS quantifies the number of operations that can be handled by the hardware per unit of time. Hence, altogether, they qualify the quality and speed of quantum computation. Link:
Best practices for portfolio optimization by quantum computing ... - Nature.com