Benchmarking Quantum Solvers in Noisy Digital Simulations for Financial Portfolio Optimization

TL;DR

This paper compares quantum algorithms QITE and QAOA for portfolio optimization, analyzing their performance and robustness on noisy quantum hardware, highlighting trade-offs between scalability and noise tolerance.

Contribution

It provides the first systematic benchmarking of QITE and QAOA on a practical financial problem under realistic noisy quantum conditions.

Findings

QAOA converges well in noiseless settings.

QITE is more robust to noise but computationally intensive.

QAOA scales better and remains effective with noise mitigation.

Abstract

In this work, we benchmark two prominent quantum algorithms: Quantum Imaginary-Time Evolution (QITE) and the Quantum Approximate Optimization Algorithm (QAOA) for obtaining the ground state of Ising-type Hamiltonians. Specifically, we apply them to the Markowitz portfolio optimization problem in quantitative finance, on both digital quantum computers and local quantum simulators with controllable two-qubit errors (noise). In noiseless settings, we find that QAOA achieves excellent convergence to the optimal results. Under noisy conditions, the QITE method exhibits greater robustness and stability, though it incurs substantially more classical numerical cost. In contrast, we demonstrate that QAOA offers better scalability and can still yield robust results if the noise can be effectively mitigated. Our findings provide valuable insights into the trade-offs between scalability and noise tolerance and demonstrate the practical potential of quantum algorithms for solving real-world optimization problems on near-term quantum devices.

Tables

Table 1. Table 1: Comparison between the two approaches: QAOA and QITE for portfolio optimization and other problems based on finding Ising ground states.

Methods	QAOA	QITE
Noiseless Performance (FIG. 5)	Excellent convergence to exact ground states	Slight discrepancies due to imperfect circuit optimization
Noisy Performance (FIG. 6)       Two-qubit gate error rate $\sim {10}^{-2}$	Unstable; fails to converge reliably	Robust; optimal state still identifiable
Resource Requirements	High quantum resources (extensive circuit evaluations); light classical cost	Moderate quantum resources; heavy classical optimization required (Appendix. S4)
Scalability on NISQ Devices (FIG. 4)	Better scalability if noise is mitigated	Limited by classical optimization cost

Equations

$$F(z)={\theta }_1\sum _{u=1}^m{r}_u{z}_u-{\theta }_2\sum _{u,v=1}^m{c}_{u,v}{z}_u{z}_v-{\theta }_3\sum _{u=1}^m(z_u{p}_{w}b-b{)}^2.$$

$$r_u=\frac{b{p}_w}{a_{u,N_f}}\cdot \sum _{l=1}^{N_f}\frac{a_{u,l+1}-a_{u,l}}{N_f-1},$$

$$c_{u,v}=\frac{b^2{p}_w^2}{a_{u,N_f}{a}_{v,N_f}}\cdot \frac{{\sum }_{l=1}^{N_f}(a_{u,l}-{\bar{a}}_u)(a_{v,l}-{\bar{a}}_v)}{N_f-1}.$$

$$z_u=\sum _{k=1}^w{2}^{k-1}{x}_{i(u,k)},$$

$$H=\sum _i{h}_i{\sigma }_i^z+\sum _{i,j}{J}_{i,j}{\sigma }_i^z{\sigma }_j^z+\delta$$

$$C^0(\mathbf{s})=⟨\Psi (\mathbf{s})|H_p|\Psi (\mathbf{s})⟩,$$

$$C(\mathbf{s})=|E_g-E(\mathbf{s})|.$$

$${|\psi ⟩}_{\mathbf{ground}}\approx e^{-\beta H}|\psi (0)⟩/\sqrt{||e^{-\beta H}|\psi (0)⟩||},$$

$$U_{3}(\theta,\phi,\lambda)=\left[\begin{array}[]{cc}\cos\left(\frac{\theta}{2}\right)&-e^{i\lambda}\sin\left(\frac{\theta}{2}\right)\\ e^{i\phi}\sin\left(\frac{\theta}{2}\right)&e^{i(\phi+\lambda)}\cos\left(\frac{\theta}{2}\right)\end{array}\right],$$

$$C^{\text{QITE}}(\mathbf{s})=1-\mathrm{Tr}[V^{†}(\mathbf{s})U_{\mathrm{circuit}}]/2^n$$

$$=1-⟨{\psi }_h|V^{†}(\mathbf{s})U_{\mathrm{circuit}}|{\psi }_h⟩$$

$$C^{\text{QAOA}}(\mathbf{s})=|E_g-⟨\psi (\mathbf{s})|H|\psi (\mathbf{s})⟩|,$$

$$ECR=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cccc}0&1&0&i\\ 1&0&-i&0\\ 0&i&0&1\\ -i&0&1&0\end{array}\right).$$

$$F^{\mathrm{Error}}=F^{\mathrm{ideal}}-F^{\mathrm{circuit}},$$

$$z_u=\sum _{k=1}^w{2}^{k-1}{x}_{i(u,k)},$$

$$\sum _{u=1}^m{r}_u{z}_u=\sum _{u=1}^m{r}_u\left(\sum _{k=1}^w{2}^{k-1}{x}_{i(u,k)}\right)=\sum _{u=1}^m\sum _{k=1}^w{2}^{k-1}{r}_u{x}_{i(u,k)}.$$

$$-\sum _{u,v}^m{c}_{u,v}{z}_u{z}_v=-\sum _{u,v}^m\sum _{k=1}^w\sum _{k^{\prime }=1}^w{2}^{k-1}{2}^{k^{\prime }-1}{x}_{i(u,k)}{x}_{j(v,k^{\prime })}.$$

$$-(\sum _{u=1}^m{p}_{w}b{z}_u-b{)}^2=-(\sum _{u=1}^m\sum _{k=1}^w{2}^{k-1}{p}_{w}b{x}_{i(u,k)}-b{)}^2,$$

$$\underset{z}{\max }[{\theta }_1\sum _{i=1}^n{r}_i{x}_i-{\theta }_2\sum _{i,j=1}^n{c}_{i,j}{x}_i{x}_j-{\theta }_3(\sum _{i=1}^n{2}^{k-1}b{p}_w{x}_i-b{)}^2],$$

$$n$$

$$r_i$$

$$c_{i,j}$$

$$\mathop{min}\limits_{x}\Big{(}\sum_{i}^{n}q_{i}x_{i}+\sum_{i,j}^{n}Q_{i,j}x_{i}x_{j}+\gamma\Big{)},$$

$$q_i$$

$$Q_{i,j}$$

$$\gamma$$

$$H=\sum _i{h}_i{\sigma }_i^z+\sum _{i,j}{J}_{i,j}{\sigma }_i^z{\sigma }_j^z+\delta ,$$

$$J_{i,j}$$

$$h_i$$

$$\delta$$

$$U=\left[\begin{array}[]{cc}uU_{\rm non}&B\\ C&D\end{array}\right].$$

$$U^{\prime}=\left[\begin{array}[]{cc}uU_{\rm non}&I\\ C&I\end{array}\right],$$

$$U^{\prime }=UM,$$

$$U\big{(}\ket{\psi}\otimes\ket{0}\big{)}=\big{(}u\,U_{\mathrm{non}}\ket{\psi}\big{)}\otimes\ket{0}+\big{(}C\ket{\psi}\big{)}\otimes\ket{1},$$

Full text

I Introduction

Leveraging on quantum computing for real-world applications represents a promising yet challenging frontier [1, 2, 3, 4, 5, 6, 7, 8]. With its ability to process and manipulate information with qualitatively heightened efficiency, quantum processors have significant potential to outperform classical resources in solving complex problems across domains from physics [9, 10, 11, 12, 13, 14, 15, 5, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] to practical applications such as cryptography and drug discovery [29, 30, 31, 32, 33, 34]. More recently, quantum computing has demonstrated substantial promise in tackling computationally demanding optimization problems [35, 36, 37, 38, 39, 40]. A particularly relevant example with practical importance is portfolio optimization in quantitative finance [41, 42, 43, 44]. This task involves constructing a portfolio that maximizes expected returns, a problem traditionally addressed using classical optimization algorithms and machine learning techniques [45, 43, 44, 46]. As the complexity and dimensionality of modern financial portfolios increase, quantum hardware is emerging as a compelling tool for managing sophisticated optimization models.

To solve such optimization problems on digital quantum computers, one first needs to formulate them within a quantum Hamiltonian framework, such that the optimal solution corresponds to the ground state of a spin-1/2 system. A widely used method is Quantum Approximate Optimization Algorithm (QAOA) [47, 35, 48, 36, 49, 50]. This approach is implemented through optimizing a parameterized circuit and minimizing the expectation value of a cost function. Beyond these methods, Quantum Imaginary-Time Evolution (QITE) provides a deterministic route to ground-state preparation by evolving a quantum state along imaginary time [51, 52, 53].

While previous works have explored quantum solvers for portfolio optimization problems [47, 54, 44, 55], most remain constrained by limitations inherent to current quantum hardware [56, 57, 58, 59, 60], and few explicitly examine the impact of noise. To assess the practical performance of these methods in noisy simulations, we conduct a comprehensive benchmarking study of these quantum algorithms applied to the Markowitz portfolio optimization problem. Our results show that QAOA performs exceptionally well under ideal, noiseless conditions; however, it exhibits poor robustness in the presence of realistic noise [61, 62, 63]. In contrast, QITE exhibits greater robustness to noise, but requires extensive training of parameterized circuits to accurately approximate the desired imaginary-time evolution. Overall, our findings suggest that current quantum optimization strategies are still fundamentally limited by challenges related to noise resilience and scalability. While both methods face limitations on noisy devices, QAOA shows strong potential for large-scale applications if hardware noise can be significantly mitigated through advanced error mitigation techniques or future fault-tolerant architectures [56, 64, 65, 66, 67, 68].

This work is organized as follows: In Section II, we introduce the problem of financial portfolio optimization and the potential advantages of tackling it through a quantum algorithm. Following this, in Section III, we introduce the Markowitz portfolio optimization problem and detail its mapping to a quantum Hamiltonian. In Section IV, we describe the quantum optimization algorithms employed in this study, including QAOA and QITE. Then, in Section V, we present simulation results and compare the performance of QAOA and QITE under both noiseless and noisy conditions. Finally, we summarize our benchmarking results of the two quantum algorithms and discuss applications for future research in Section VI.

II Portfolio optimization problem

In the financial sector, there exist a variety of optimization problems [69, 70, 71, 72]. These include forecasting tasks such as pricing, risk assessment, anomaly detection in transactions, and identifying customer preferences, as well as optimization tasks like selecting optimal portfolios, developing effective trading strategies, and designing robust hedging techniques. Portfolio optimization specifically refers to the problem of identifying the optimal allocation of a given set of assets [44, 73, 43]. In this section, we provide an introduction to the problem of Markowitz portfolio selection.

II.1 Markowitz portfolio selection

Here, we first introduce the mathematical model for Markowitz portfolio selection [74, 75, 76, 77]. In this framework, a fixed budget $b$ can be allocated among $m$ assets. Each asset $u$ is characterized by its expected return $r_{u}$. The model also incorporates a covariance matrix $c_{u,v}$, which quantifies the historical price correlation between assets $u$ and $v$. Here, the task of Markowitz portfolio selection is to construct an optimal portfolio that maximizes expected return while minimizing overall portfolio variance. This formulation is based on Markowitz portfolio theory, which aims to achieve an efficient trade-off between risk and return [43].

Next, we introduce the quantized minimum investment fraction $p_{w}=1/2^{w-1}$, where $w\in\mathbb{N}$ parametrizes the number of budget slices $2^{w-1}$. The investment fraction for the $u$-th asset can be described as an integer multiple $z_{u}p_{w}$, where the integer $z_{u}$ indicates how many quantized units are allocated to asset $u$. We then impose the constraint that the total investment must match the full budget: $\sum_{u}z_{u}p_{w}=1$. Within this framework, the Markowitz portfolio selection problem is reformulated as an optimization task: identifying the optimal set of $z_{u}$ values that maximizes the following objective function [43]

[TABLE]

Here, $b$ represents the real amount of budget, while the product $z_{u}p_{w}$ ($0\leq z_{u}p_{w}\leq 1$) denotes the fraction of the budget allocated to asset $u$. The parameters $\theta_{1,2,3}$ represent the investor’s preferences with respect to expected return, portfolio risk, and tolerance for budget deviation, respectively. Specifically, the $\theta_{1}$ term rewards high expected return, the $\theta_{2}$ term penalizes risk arising from asset correlations, and the $\theta_{3}$ term imposes a penalty for deviation from the total budget [43]. The inclusion of the factor $b$ in the $\theta_{3}$ term aligns it with the budget dependence of $r_{u}$ [Eq. 2] and $c_{u,v}$ [Eq. 3] in the first two terms. This ensures that the budget deviation penalty remains comparable in magnitude to the return and risk terms, even when $b$ is large.. In our implementations, we set $\theta_{1}=0.8$, $\theta_{2}=0.1$, $\theta_{3}=0.1$, and $b=10$. More details are presented in Appendix. S2.

II.2 Correspondence between Financial and Quantum Models

To solve this optimization problem described by Eq. 1 on a quantum platform, we first map it onto the ground-state problem of a corresponding quantum Hamiltonian. To achieve this, we rescale the expected return based on historical price data and define

[TABLE]

where $a_{u,l}$ is the $l$-th history of price data for the $u$-th asset and $N_{f}$ denotes the number of historical assets price points [78, 79]. In this work, these historical prices are randomly generated, and details are provided in Appendix. S2. The corresponding covariance matrix $c_{u,v}$ is computed as

[TABLE]

where $\bar{a}$ denotes the mean price of assets. To enable encoding these data on a qubit system, we represent each integer $z_{u}$ in binary form:

[TABLE]

where $i(u,k)=(u-1)w+k$ and $x_{i}\in\{0,1\}$ is a binary variable indicating the presence or absence of investment in a specific quantized unit.

For instance, the optimal bitstring shown in FIG. 6 is $010100100$ ($101011011$ under little-endian qubit ordering). This solution encodes an investment configuration for a portfolio of three assets with three binary slices per asset. The actual investments are $z_{0}/p_{w}=1/2$, $z_{1}/p_{w}=1/4$, and $z_{2}/p_{w}=1/4$.

Under this representation, the original optimization problem is reduced to finding the ground state of an Ising Hamiltonian:

[TABLE]

where classical binary quantities are mapped to quantum operators via the transformation $x_{i}\rightarrow\frac{1}{2}(1+\sigma_{i}^{z})$, with the Pauli operator $\sigma^{z}_{i}$. The coefficients $h_{i}$, $J_{i,j}$ and $\delta$ canbe derived from the parameters appearing in Eq. 1, which are computed from the asset prices, as detailed in Appendix S1.

There exists previous works that demonstrated how this portfolio optimization problem can be addressed using quantum methods [80, 81, 43, 82, 47, 54, 44, 55, 40]. However, most focus on noiseless settings, which do not accurately reflect the limitations of noise in current platforms, including gate errors, decoherence, and imperfect measurements. Given these, we implement and benchmark different quantum algorithms on noisy simulators and demonstrate whether these methods can identify optimal results under certain controlled noise levels. In the following section, we will provide a comprehensive discussion of these quantum algorithms and their effectiveness in optimizing financial portfolios.

III Methods

In this section, we outline the methods employed to solve the optimization problem described above using digital quantum processors and their simulators.

III.1 Quantum optimization algorithms

The Quantum Approximate Optimization Algorithm (QAOA) is a variational algorithm that can be employed to approximate the ground state of a given Hamiltonian [80, 35]. This method operates by minimizing the expectation of the target Hamiltonian $H_{p}$ through a parameterized quantum circuit. Generally, this is achieved by minimizing the cost function

[TABLE]

but we later introduce a modification to this form in order to accelerate convergence. Here, $\bm{s}$ is a generalized vector containing the entire set of variational parameters $(s_{0},s_{1},...)$, where each $s_{i}$ denotes the angles in each rotational gate, as defined in Eq. 9. These variational angles are optimized by minimizing the difference between the exact ground-state energy $E_{g}$ and the estimated energy $E(\bm{s})$ i.e. by minimizing the cost function

[TABLE]

This setup is schematically shown in FIG. 1 (a). Note that $E_{g}$ is the exact ground level, which is a relative energy shift. In general, this term is not strictly necessary, since QAOA can be implemented by directly minimizing the cost function defined in Eq. 6. However, incorporating the relative shift often accelerates convergence, thereby reducing the number of iterations and lowering the overall simulation cost.

To simulate this optimization process on quantum circuits, we first evaluate the expectation value $E(\bm{s})$, obtained by measuring the expectation of each local term in the Hamiltonian. We then iteratively adjust the parameters in $\bm{s}$ to minimize the deviation $|E_{g}-E(\bm{s})|$. After each update, the new parameters are re-uploaded to the quantum circuit, and the cycle of measurement and optimization is repeated until the cost function reaches convergence within a reasonable tolerance. For the classical optimization step, we primarily employ the COBYLA (Constrained Optimization BY Linear Approximations) optimizer to achieve optimization [83]. As an alternative, we also test the L-BFGS-B (Limited-memory Broyden–Fletcher–Goldfarb–Shanno with Box constraints) algorithm [83]. IWhile both methods are viable, we consistently observe that COBYLA achieves faster and more reliable convergence in practice.

Through an alternating evolution process, QAOA efficiently explores the Hilbert space to approximate the ground state. The specific structure of the variational circuit used in this work is detailed in section III.3. While QAOA has demonstrated strong effectiveness and versatility across a variety of optimization tasks, its performance is not robust in the presence of noise in tasks that require high-fidelity identification of optimal quantum states, as we will demonstrate below. In the following section, we conduct a detailed evaluation of QAOA performance in solving this problem, with a particular focus on its robustness in noisy environments and its effectiveness in approximating optimal quantum states.

III.2 Quantum imaginary-time evolution

We next discuss the Quantum Imaginary Time Evolution (QITE) method. This approach obtains the ground state of a quantum system via the following evolution:

[TABLE]

where $H$ is the system Hamiltonian, and $\beta$ is the imaginary time parameter [51, 52]. As $\beta$ increases, excited states are suppressed, and the evolution converges towards the ground state of the system. Since this process is intrinsically non-unitary, its implementation on quantum hardware requires an effective embedding strategy: the non-unitary operator $e^{-\beta H}$ is realized through an extended unitary operation which we denote as $U_{\text{circuit}}$, as schematically illustrated in FIG. 1 (b). In this setup, an ancilla qubit (shown in green) is coupled to the physical qubits to implement the QITE process. Then, post-selection on an ancilla qubit enables the extraction of the desired state $e^{-\beta H}\ket{\psi(0)}/\sqrt{||e^{-\beta H}\ket{\psi(0)}||}$. Further technical details are provided in Appendix S4.

III.3 Variational circuits for implementing QAOA and QITE

To implement QAOA and QITE on an actual quantum circuit, we employ variational circuits i.e. parametrized quantum circuits to be optimized. The structure of a variational circuit $V$ is shown in FIG. 2, where the variational parameters enter each of the $U_{3}$ gates

[TABLE]

which is parameterized by three tunable angles $\theta,\phi,\lambda$. The entire variational circuit $V$ is thus parameterized through the collection of all the vectors $s_{i}=(\theta_{i},\phi_{i},\lambda_{i})$.

Variational circuits are used differently in QITE and QAOA. Exactly representing the target unitary $U_{\rm circuit}$ used in the QITE process (corresponding to the circuit shown in FIG. 1 (b)) as a quantum circuit is typically infeasible, as it may involve circuit depths that are orders of magnitude greater than what is allowed by current NISQ noise thresholds. However, a properly optimized variational circuit $V$ can approximate $U_{\rm circuit}$ rather accurately through the minimization of the following cost function

[TABLE]

where $n$ is the number of qubits in the entire circuit, and $\ket{\psi_{h}}=\text{Hadamard}^{\otimes n}\ket{000...}$. In the well-optimized limit where $C^{{\text{QITE}}}(\bm{s})\approx 0$, we have $V\approx U_{\rm circuit}$.

For the QAOA algorithm, the variational circuit $V$ is optimized to minimize the discrepancy between the energy expectation and the ground state energy, instead of maximizing its overlap with a desired evolution operator. We write the QAOA cost function Eq. 7 as

[TABLE]

where the expectation of the Hamiltonian $H$ is taken with respect to ${\ket{\psi(\bm{s})}=}V(\bm{s}))\ket{\psi_{0}}$, the state evolved from $\ket{\psi_{0}}=\ket{000...}$ under the variationally optimized operator $V$.

From FIG. 2, the variational circuit also consists of (non-parametrized) two-qubit entangling gates known as Echoed Cross Resonance (ECR) gates, which are natively supported on the IBM Quantum processors that we use. Functionally equivalent to CX gates, they are defined by

[TABLE]

However, in our local noisy simulations, it is more straightforward to include noise effects in CX gates through imperfect X operations. Therefore, we use CX gates for all local simulations in this work, with each ECR gate decomposed as $(H\otimes I)CX(I\otimes S^{\dagger})CX(H\otimes X)$, with $S^{\dagger}$ denoting the S-adjoint gate.

In our study, the implementation of QITE requires deeper circuits with a larger number of trainable parameters. The 9-qubit instance shown later requires only 2 layers in the QAOA circuit, whereas the corresponding QITE circuit demands between 4 and 8 layers to achieve convergence. Thus, QITE in turn results in a substantially higher classical optimization cost

IV Benchmarking Results

IV.1 Benchmarking simulations of QAOA

In this section, we first evaluate the performance of the QAOA method, with the structure of the variational circuits detailed in FIG. 2. During the optimization process, the cost function, defined by the energy difference in Eq. (7), is minimized to approximate the ground state of the target Hamiltonian.

IV.1.1 Noiseless QAOA simulations

As depicted in FIG.3 (a), we implement the QAOA procedure using a noiseless local simulator. Here, in each optimization iteration, we perform a full round of circuit simulation (distinct from an actual circuit run in Qiskit) and subsequently update the trainable parameters based on the feedback. For these noiseless simulations, the optimization exhibits near-perfect convergence towards $E_{g}$ after approximately $600$-$1000$ iterations. However, this high accuracy comes at the cost of significant quantum resources, requiring over $600$ circuit evaluations on quantum hardware. The lowest discrepancy ${\rm Min}|E_{g}-E_{i}|$ from the true ground state energy $E_{g}$ achieved throughout the QAOA optimization process in each instance are presented in FIG. 3 (b). Here, the optimized circuits converge to $C^{{\text{QAOA}}}(\bm{s})<2.5$, indicating good approximations to the ground state, as further detailed in Appendix S6. Note that different problem instances correspond to Hamiltonians with distinct parameter sets derived from independently generated asset data, as detailed in Appendix S2. Consequently, their ground states are associated with different optimal bitstrings. That said, even in the absence of noise, the finite number of circuit executions introduces intrinsic sampling errors, preventing convergence to the exact ground state.

IV.1.2 QAOA simulations with controlled noise

We then conduct controlled noisy simulations of QAOA on a local simulator, and the effect of CX gate noise on the trajectory of this optimization process is presented in FIG. 3 (c). Under increasing CX noise (error rates), the optimization becomes highly unstable and exhibits no clear convergence, highlighting the significant impact of noise on QAOA performance. This trend is further illustrated in Fig.3(d), where it is evident that the minimum energy achieved during optimization increases rapidly with growing CX gate error rates. Even for a $0.001$ CX error rate, which is comparable to that in contemporary IBM quantum processors, $E_{g}-E_{i}$ is already far larger than that from noiseless simulations.

In this work, QAOA optimization is performed using the COBYLA optimizer [83], while the corresponding quantum state measurements are collected from a noisy simulator. While COBYLA is designed to handle constrained problems efficiently, it relies on iterative evaluations of the cost function. When those evaluations are affected by hardware noise, the optimizer struggles to navigate the resulting rugged energy landscape, often leading to poor convergence behavior.

IV.1.3 QAOA simulations – scaling to larger systems

To further evaluate the scalability of QAOA, we simulate its performance on 20- and 30-qubit systems using a matrix product state (MPS) simulator. Here, we present the results of two instances in FIG. 4. For 20-qubit instances, the QAOA optimization process exhibits relatively stable convergence, and the optimal outcome $|E_{g}-E_{i}|\approx 5$ is very small compared to the overall energy scale. Despite notable fluctuations during the optimization process of the 30-qubit instances, the final stage can reach a similar energy level of $|E_{g}-E_{i}|\approx 10$, comparable with the results in smaller systems. These results demonstrate that QAOA can effectively approach optimal solutions even at intermediate system sizes. Since poor convergence under noise has already been observed in smaller systems, and noisy simulations at the 20-qubit scale lead to substantial computational cost, we do not present noisy simulations for these larger instances.

IV.2 Comparison between QAOA and QITE

IV.2.1 Benchmarking from predicted returns

In this section, we compare the performance of QAOA and QITE in terms of their predicted returns $F$, as defined and computed by Eq. 1. Using the same representative instances, we compute the return error between the circuit outputs and the exact (ideal) results, defined as

[TABLE]

where $F^{\rm ideal}\approx 0$, and $F^{\rm circuit}<0$ denotes loss in portfolio value due to suboptimal outcomes. FIG. 5 illustrates the simulated results of such returns. FIG. 5 (a) demonstrates that, remarkably, the QITE algorithm executed on IBM quantum hardware (with $\approx 0.007$ ECR/CX error rate, as obtained from Fig. S1 from Appendix. S3) predicts return values that exhibit very low $F^{\rm Error}<1$, almost comparable to that from the QAOA instances obtained in ideal (noiseless) settings. Note that because the implementation of QAOA requires a large number of circuit runs, so to save cost on IBM Quantum hardware, we only conduct QAOA simulations on local simulators which can be calibrated to realistic noise levels of actual quantum hardware (see Appendix. S3).

To put the excellent performance of the QITE hardware runs into perspective, we compare the return error $F^{\rm Error}$ from QAOA with simulated noise of similar magnitude [FIG. 5 (b)]. Saliently, even the smallest levels of noise significantly affect QAOA performance, with return errors exceeding $F^{\rm Error}>20$, and increasing further as noise levels rise. Notably, the optimal results obtained from QAOA are not stable: under a moderately elevated noise level ($0.011$), the first two instances exhibit a sharp increase in return error. This instability arises because noise can introduce temporal fluctuations in the measured cost function, leading to a non-smooth landscape in optimization. As a result, the performance of QAOA becomes highly sensitive to both the specific problem instance and the stochastic behavior of the circuit runs.

To further contextualize our $F^{\rm Error}$ results, we also compute returns from reference states with uncorrelated uniformly random coefficients, which represent completely unoptimized outcomes 111These random states are not produced by random circuits; instead, we generate them numerically as $\ket{\psi}=\sum_{i}a_{i}\ket{\psi_{i}}/||\sum_{i}a_{i}\ket{\psi_{i}}||$, which are normalized. Here, the coefficients $a_{i}$ are randomly sampled from a uniform distribution within $[0,1]$, and $\ket{\psi_{i}}$ range over the computational basis states.. As shown in FIG.5(b) (cyan bars), these unoptimized reference random states exhibit large errors of $F^{\rm Error}\approx 38$, which substantially exceed those in FIG.5(a), and are also larger than the low-noise QAOA results also plotted in FIG. 5(b). However, under significant but realistic noise levels, such as a CX error rate of $0.011$, QAOA performance error approaches the level of these fully unoptimized outcomes, offering no optimization value at all.

IV.2.2 Benchmarking from measured optimal state distribution

We now examine the actual solution states obtained from both QITE and QAOA, which represents the ideal portfolio composition. In principle, the ideal optimal solution corresponds to a single bitstring with the highest return. Thus, to assess how the outputs of QAOA and QITE deviate from the ideal solution, we examine a representative problem instance and analyze the full probability distributions of measured quantum states.

As shown in FIG. 6, the QITE results (green bars) obtained from simulations on IBM hardware exhibit a reasonably prominent peak at the bitstring 101011011 corresponding to the exact optimal solution. Although the obtained state probability distribution is broadened due to the presence of noise, the optimal bitstring remains clearly identifiable. For comparison, the blue bars in FIG. 6 display the probability distribution obtained from the noiseless QAOA simulation with the same instance. In this ideal noiseless setting, QAOA yields a nearly perfectly localized state on this optimal bitstring (see FIG. S6) In contrast, QAOA results (red bars) conducted on a local simulator under realistic errors (see Appendix. S3), fail to reveal any distinct peak corresponding to the optimal solution. This behavior highlights the instability under noise, consistent with the unstable convergence trends observed in FIG. 3.

To further validate the robustness of QITE, we simulate the QITE circuit for additional 50 portfolio optimization instances on the IBM quantum computer, also with $L=9$ qubits. These simulations are performed on both IBM Quantum hardware. As shown in FIG. 7, the return errors $F^{\text{Error}}$ range from about 1 to 6, which is much smaller than those from QAOA with similar simulated CX gate error rate of $0.007$ (see Appendix. S3), as previously shown in FIG. 5. Note that the return error range can arise from variations in the optimization quality of circuits across different instances, as well as from noise fluctuations in simulations. FIGs.7(b,c) present the measured solution state distributions for two representative instances drawn from FIG.7(a). Despite higher return errors of $2.66$ and $6.17$, the optimal bitstring still remain prominent.

V Conclusion

In this work, we benchmarked two quantum algorithms, QAOA and QITE, which are widely used in digital quantum simulations for solving not just the Markowitz portfolio optimization problem, but also other diverse problems with quantum Ising formulations [85, 86, 87, 88]. Our results show that the QAOA method achieves excellent convergence to the exact ground state in noiseless simulations when hardware noise is negligible or extensively mitigated. However, the QAOA method is substantially sensitive to noise, particularly due to its reliance on deep parameterized circuits and repeated quantum evaluations, which are highly sensitive to gate and readout errors. In contrast, the QITE method exhibits much stronger robustness to realistic noise levels, achieving comparably optimal solutions on actual quantum hardware runs as noiseless QAOA simulated runs.

From a resource perspective, as touched on earlier in Section. III.3, QAOA demands extensive quantum resources, typically requiring hundreds to thousands of circuit executions, but relatively low classical computational cost. However, QITE, on the other hand, requires less quantum cost but substantial classical computational cost due to the need for extensive circuit training. As illustrated in FIG. 3 and FIG. 4, the 9-qubit instance of QAOA requires a few hundred iterations for reasonable convergence, which increases to around 1000 iterations for the 20-qubit case. For QITE, only a single round of circuit simulation is required, equivalent to the cost of one QAOA iteration, resulting in markedly lower quantum resource requirements. For problems of the same scale, QAOA requires hundreds of times more quantum resources than QITE.

On quantum hardware, qubit measurements are straightforward and can be performed on systems with over 100 qubits. As a result, in terms of scalability, QAOA holds greater theoretical promise for tackling larger problem instances, if hardware noise can be effectively mitigated. A detailed side-by-side comparison of both methods is summarized in Table 1.

Overall, our comparative analysis of the noise resilience and scalability in quantum optimization strategies suggests that the choice between QAOA and QITE should be guided by the specific noise characteristics of the quantum device, the availability of classical computational resources, and the complexity and scale of the optimization task [37, 89, 90, 55, 91, 92, 93].

VI Acknowledgements

We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Quantum team. All data and codes of this work are available from the corresponding authors upon reasonable request. This work is supported by the Singapore Ministry of Education Academic Research Fund Tier-I preparatory grant (WBS no. A-8002656-00-00) and MOE’s Tier-II grant (WBS no: A-8003505-00-00).

Appendix S1 Correspondence between the Markowitz Portfolio model and the Ising Hamiltonian

In this section, we detail the construction of the Ising Hamiltonian for the Markowitz portfolio selection problem described by Eq. 1. We begin by reformulating the total expected return using binary variables, a necessary step for quantum computing applications [43].

We use $z_{u}$ to represent the investment in asset $u$, quantified as an integer number of quantized minimum investment fractions. We can express each integer $z_{u}$ using a standard binary representation with $w$ binary variables, $x_{i(u,k)}\in\{0,1\}$:

[TABLE]

where $i(u,k)=(u-1)w+k,k=1,2,...,w.$ In this expression, the variables $\{x_{i(u,1)},\dots,x_{i(u,w)}\}$ are the binary digits of $z_{u}$. For example, if we use $w=3$ bits to represent the investment in the $u-$th asset, an amount of $z_{u}=5$ (which is $101_{2}$ in binary) would be encoded by setting $x_{i(u,1)}=1$, $x_{i(u,2)}=0$, and $x_{i(u,3)}=1$, since $5=(1\cdot 2^{0})+(0\cdot 2^{1})+(1\cdot 2^{2})$. Substituting this binary expansion of $z_{u}$ [Eq. S1] into the formula for total expected return in Eq. 1 yields :

[TABLE]

Here, $m$ is the total number of assets, $w$ is the number of binary variables (or ”qubits”) allocated to each asset, and $r_{u}$ is the expected return for the $u$-th asset. The function $i(u,k)=(u-1)w+k$ simply assigns a unique index to each binary variable across all assets.

The corresponding covariance term is rewritten as follows

[TABLE]

The allocation constraint is given by:

[TABLE]

where $b$ is the budget and $p_{w}=1/2^{w-1}$ is the quantized minimum investment fraction.

Based on such a binary representation, Eq. 1 from the main text can be rewritten as:

[TABLE]

with the following rescaling parameters:

[TABLE]

The above optimization is equivalent to a Quadratic Unconstrained Binary Optimization (QUBO) problem, expressed as:

[TABLE]

where the parameters are defined as:

[TABLE]

By transforming $x_{i}=\frac{1}{2}(1+\sigma_{i}^{z})$, the above optimization can be reduced to finding the ground state of the following Ising Hamiltonian:

[TABLE]

with the corresponding parameters:

[TABLE]

In our implementation, we discard the constant term $\delta$, as it does not affect the optimization process.

Appendix S2 Dataset Preparation

In this section, we outline the method used to prepare our dataset, which simulates the historical prices for financial assets. The asset price dataset is structured as a matrix of dimensions $(N_{f},m)$, where $N_{f}$ represents the number of historical price points per asset, and $m$ is the number of assets under consideration. Each column in this matrix corresponds to the price history of a distinct asset.

The dataset preparation process begins by setting the price generator $a_{u,0}$, for the $u$-th asset. This price generator is randomly selected from the range [$b/10,b$], where $b$ is the total budget. To generate brownian-like temporal price fluctuations, each subsequent historical price point, $a_{u,l}$, is then generated within a fluctuation range of $75\%$ to $125\%$ of the price generator, specifically, $a_{u,l}=(1+\alpha)a_{u,0}$. This procedure generates the $u$-th column of the price data matrix. To construct the full matrix, this process is repeated across all assets. Based on the generated price data, the expected return $r_{u}$ and covariance matrix $c_{u,v}$ are calculated using Eq. 2. and Eq. 3 respectively.

With the calculated price data, we proceed to generate the Ising Hamiltonian problem instances. This is achieved by selecting one combination of the parameters $\theta_{1},\theta_{2},\theta_{3}$ subject to the constraint $\sum_{i}\theta_{i}=1$, where different choices of $\bm{\theta}$ correspond to distinct investment preferences. These parameters, which represent different investor preferences, are integrated with the calculated $r_{u}$ and $c_{u,v}$. Subsequently, the rescaling process described in Eq. S6 is applied. A detailed workflow (9-qubit instance in the main text) is demonstrated as follows, and the corresponding code and data are available at [94]:

Input parameters: The number of assets is $m=3$, with each asset having $N_{f}=100$ price points. Each asset is represented using $w=3$ binary bits (slices), with $p_{w}=1/2^{w-1}$. A total of 100 problem instances are generated, with a total budget $b=10$. Additional parameters are set to $\theta_{1}=0.8$, $\theta_{2}=0.1$, and $\theta_{3}=0.1$. 2. 2.

Price generation: For each asset $u$, randomly initialize the price generator $a_{u,0}$ within the range $[b/10,b]$. Generate historical prices using $a_{u,l}=(1+\alpha)a_{u,0}$, where $\alpha$ is a random variable uniformly sampled from $[-0.25,0.25]$. 3. 3.

Covariance matrix and expected return: For each instance:

• •

Compute the covariance matrix $c_{uv}$ as defined in Eq. 3 of the main text:

$c_{u,v}=\frac{b^{2}p_{w}^{2}}{a_{u,N_{f}}a_{v,N_{f}}}\cdot\frac{\sum_{l=1}^{N_{f}}(a_{u,l}-\bar{a}_{u})(a_{v,l}-\bar{a}_{v})}{N_{f}-1},$ where $\bar{a}$ denotes the mean price of the assets and $p_{w}=2^{1-w}$.

• •

Compute daily returns $\Delta a_{u,l}=a_{u,l+1}-a_{u,l}$ and the expected return $r_{u}$ from Eq. 2 of the main text: $r_{u}=\frac{bp_{w}}{a_{u,N_{f}}}\cdot\sum_{l=1}^{N_{f}}\frac{a_{u,l+1}-a_{u,l}}{N_{f}-1}.$ 4. 4.

Rescaling: Rescale $c_{uv}$ and $r_{u}$ to $Q_{i,j}$ and $q_{i}$ (parameters for QUBO problems) according to Eq. S6 and Eq. S8. 5. 5.

Ising coupling and field parameters: The total number of binary variables is $m\cdot w=9$. For binary indices $i,j\in\{1,\dots,9\}$, the corresponding Ising coupling $J_{ij}$ and Ising field $h_{i}$ are computed based on Eq. S10: $J_{i,j}=\frac{1}{4}Q_{i,j},h_{i}=\frac{1}{2}q_{i}+\frac{1}{2}\sum_{j}Q_{i,j}$.

In this work, for the 9-qubit problem, we generate 100 distinct instances, each associated with a specific Hamiltonian. Consequently, the optimization landscape varies across instances, and the optimal solution bitstring is not universally fixed.

Appendix S3 IBM quantum device used

For simulations in this work, we implement jobs on the IBM ”Brisbane” device, and details of the gate errors are shown in FIG. S1. Here, we select the qubit chain $[0,1,2,3,4,5,6,7,8,9,10,11,12]$ to be used for our computations. Note that simulations shown in FIG. 6 are based on a local noise model calibrated using the error parameters extracted from the “Brisbane” device [95].

Appendix S4 Quantum imaginary-time evolution

In this section, we describe our approach for simulating the non-unitary evolution operator $U_{\rm non}=e^{-\beta H}$. Since quantum circuits are inherently unitary, this non-unitary operation must be embedded into a larger unitary operator that can be physically implemented. This embedding can be constructed in the following form:

[TABLE]

where the scaling factor $u^{-2}$ denotes the maximum eigenvalue of $U^{\dagger}_{\rm non}U_{\rm non}$ and the lower-left block $C$ is numerically solved, $C=A\sqrt{I-u^{2}\Sigma^{2}}E^{\dagger}$ with identity matrix $I$. The matrix elements $A$, $E$, and $\Sigma$ are obtained from the singular value decomposition (SVD) of the operator $U_{\rm non}$: $U_{\rm non}=A\Sigma E^{\dagger}$ [23, 96], where $\Sigma$ is a diagonal matrix containing the singular values, and $A$, $E$ are unitary matrices.

To construct the full unitary operator $U$, we first perform QR decomposition on an associated ansatz matrix:

[TABLE]

where the right-side blocks are simply taken to be identity matrices. To recover $U$, we can apply the QR decomposition to $U^{\prime}$, such as to obtain

[TABLE]

where $U$ ($UU^{\dagger}=I$) is the desired unitary operator, and $M$ is an upper triangular matrix.

Here, we present a concrete example to demonstrate how $U_{\text{non}}$ can be executed from the $U$ circuit. We can consider an initial composite state $\ket{\psi}\otimes\ket{0}$, where $\ket{\psi}$ is the input state for physical problems, and $\ket{0}$ is for the ancilla qubit. The action of the unitary operator $U$ on this state can be expressed as

[TABLE]

where $U_{\mathrm{non}}$ denotes the desired non-unitary evolution. Then, we perform a projective measurement on the ancilla qubit. The post-selecting outcome with $\ket{0}$ is recorded, and the undesired branch associated with $\ket{1}$ is discarded. The resulting normalized state of system qubits is

[TABLE]

The implementation of this unitary matrix $U$ can be achieved through variational circuits. Note that the Hamiltonian defined by Eq. S9 includes long-range couplings, so the optimization of variational circuits requires deep layers. As a result, the QITE method in this work is primarily applied to small-scale systems, where circuit depth and computational cost remain manageable. For example, the 9-qubit instance presented in the main text requires only 2 layers in the QAOA circuit, whereas the corresponding QITE circuit demands between 4 and 8 layers to achieve convergence. Despite the scaling challenges, larger instances, such as those with 20 and 30 qubits, can still be addressed using shallow QAOA circuits with just 2 layers.

For the optimization of QITE circuits, we set a convergence threshold of $C^{\text{QITE}}<0.1$ [Eq. 10]. This ensures that the evolved state closely approximates the ground state, and also keeps the circuit depth moderate and computationally manageable. As shown in FIG. 7(a), such imperfect convergence can give rise to a broad spread in the return error.

Appendix S5 Quantum annealing and adiabatic evolution

Quantum annealing is another approach designed to solve optimization and minimization problems [97, 98, 99]. This method is based on the adiabatic evolution of a quantum system governed by a time-dependent Hamiltonian that interpolates between an initial Hamiltonian $H_{i}$ and the target Hamiltonian $H_{p}$. This evolution is described by: $H(\lambda(t))=(1-\lambda(t))H_{i}+\lambda(t)H_{p},$ where $\lambda(t)$ is a time-dependent annealing parameter satisfying the boundary conditions $\lambda(0)=0,\lambda(T)=1$. The system is initially prepared in the ground state of $H_{i}$ and evolves under $H_{ad}(\lambda(t))$, with the target Hamiltonian gradually turned on. According to the quantum adiabatic theorem, if the evolution proceeds slowly enough, the system remains in its instantaneous ground state and eventually transitions to the ground state of $H_{p}$.

In general, quantum annealing can be naturally implemented on analog quantum simulators, where the evolution on such quantum platforms closely follows a continuous-time annealing schedule. However, the simulator employed in this work is digital, and emulating quantum annealing dynamics on such a platform requires extremely deep quantum circuits to Trotterize the continuous evolution. These deep circuits are highly susceptible to noise and decoherence. Given this, we do not include quantum annealing in our performance benchmarks.

Appendix S6 Additional results

Here, we present additional results to demonstrate that the noiseless optimization of QAOA can reach reasonable convergence. These results are presented in FIG. S2, for the four instances shown in FIG. 5 (main text). These results clearly demonstrate that, in the absence of noise, QAOA can converge toward the ground state, successfully identifying the optimal bitstring for each instance. Note that for different instances, the corresponding Hamiltonian might also contain different parameters corresponding to different instances of generated asset data, leading to the convergence to distinctive bitstrings.