Distributed Quantum-Enhanced Optimization: A Topographical Preconditioning Approach for High-Dimensional Search Dominik Soós∗ , Marc Paterno∗† , John Stenger∗‡ , Nikos Chrisochoides∗§ ∗ Department of Computer Science, Old Dominion University, Norfolk, VA, USA † Computational Science and Artificial Intelligence Directorate, Fermi National Accelerator Laboratory, Batavia, Illinois, USA
arXiv:2604.20639v1 [quant-ph] 22 Apr 2026
‡ Chemistry Division, Naval Research Laboratory, Washington, D.C., USA § Department of Physics, Old Dominion University, Norfolk, VA, USA
Abstract—Optimization problems become fundamentally challenging as the number of variables increases. Because the volume of the search space grows exponentially, classical algorithms frequently fail to locate the global minimum of non-convex functions. While quantum optimization offers a potential alternative, mapping continuous problems onto near-term quantum hardware introduces severe scaling limits and barren plateaus. To bridge this gap, we propose the Distributed Quantum-Enhanced Optimization (D-QEO) framework. Instead of forcing the quantum processor to find the exact minimum, we use it simply as a topographical preconditioner. The QPU maps the landscape to locate the most promising basin of attraction, generating highquality seed points for a classical GPU-accelerated solver to refine. To make this approach viable for utility-scale problems, we exploit the mathematical structure of separable functions. This allows us to cut a 50-qubit (i.e., 250 ) global search space into independent and manageable sub-spaces using 5-qubit subcircuits. By executing these fragments concurrently with CUDA-Q, we completely bypass the overhead of cross-register entanglement and classical tensor knitting for separable functions. Benchmarks on the 10-dimensional Rastrigin and Ackley functions show that D-QEO prevents the exponential failure rates observed in purely classical algorithms. Furthermore, this quantum warmstart significantly reduces the number of classical BFGS iterations required to converge, providing a highly practical blueprint for utilizing near-term quantum resources in complex global search. Index Terms—quantum computing, quantum optimization algorithms,
I. I NTRODUCTION Global mathematical optimization is a vital component of modern scientific and industrial applications, ranging from machine learning and financial modeling to the simulation of complex physical systems. Despite significant progress in classical algorithmic design, the “curse of dimensionality” continues to impose a fundamental limit on the scalability of global search strategies. As the number of dimensions increases, the volume of the search space grows exponentially, making the number of sample points required to guarantee convergence to the global minimum prohibitively large. Recent advancements in classical optimization, such as the Z EUS framework [1], have sought to mitigate these challenges
by leveraging the massive parallelism of modern Graphical Processing Units (GPUs). By combining stochastic Particle Swarm Optimization (PSO) for global exploration [2] with the quasi-Newton Broyden-Fletcher-Goldfarb-Shanno (BFGS) method for local refinement [3]–[6] and utilizing forwardmode Automatic Differentiation (AD) for accurate gradient calculation [7], ZEUS provides a high-throughput pipeline for non-convex landscapes. However, even with such highperformance classical methods, the exponential scaling of the search space remains an intractable barrier (see Figure 1). Empirical studies on the multimodal Rastrigin function [8] demonstrate that as dimensionality grows, the probability of successful convergence to the global minimum basin decays exponentially. This failure mode, where the number of correct solutions reaches effectively zero by just 10 dimensions, suggests that purely classical approaches are fundamentally limited even for moderately dimensional spaces. In real-world applications, one often does not know the number of local minima present in the function being minimized. For such problems it is often necessary to try many different starting points for the minimization in order to achieve sufficient confidence that one has obtained the correct global minimum. One such example in high-energy physics is the recent neutrino mixing analyses in the NOvA experiment [9], [10]. This analysis relies on the profiled Feldman-Cousins (FC) technique [11]. In the NOvA analysis, calculating confidence intervals at the 3-sigma level required more than 20 million core-hours to generate a single two-dimensional contour plot [12]. As future high-precision experiments like DUNE [13] target 4-sigma confidence intervals, the computational effort of generating, fitting, and searching millions of pseudoexperiments across the parameter space is projected to increase by a factor of more than 40. Currently, tackling this exascale-level challenge requires developing parallel GPU-accelerated algorithms tailored for next-generation supercomputers. By utilizing the Quantum Processing Unit (QPU) as a topographical preconditioner to map global basins and warm-start classical GPU swarms, our hybrid approach directly addresses the structural
inefficiencies of these massive parameter searches, paving the way to make 4-sigma analyses and other high-demand scientific applications computationally feasible down the road. The emergence of hybrid quantum-classical algorithms offers a promising alternative to bypass these classical bottlenecks. Quantum computing introduces a Hilbert space that can inherently represent exponentially large states, providing a platform for global operations on the entire optimization landscape. However, naive translations of classical optimization to the quantum regime often struggle with the “barren plateau” problem, where gradients vanish exponentially in high-dimensional Hilbert spaces, and the high sampling costs associated with unstructured search. To achieve a true quantum advantage, it is necessary to move beyond unstructured exploration and adopt a mapping strategy that reflects the problem’s dimensionality while leveraging the mathematical structure of the objective function [14]. In this paper, we propose a “quantum-aware” framework, Distributed Quantum-Enhanced Optimization (D-QEO), which achieves scalability by utilizing a register-based dimensional mapping rather than a particle-based mapping. By assigning each of the 10 dimensions of the Rastrigin function to its own dedicated 5-qubit register, the framework utilizes a total of 50 qubits. This collective tensor-product space allows us to represent the entire 250 discretized search space in a single quantum wave-function probe. The primary novelty of our approach lies in formulating the hybrid quantum-classical workflow as a 2-level data- and circuit-decomposition [15], [16] as a preconditioner followed by a classical HPC search [1]. Rather than executing a monolithic search, we employ a search-andprune technique: the QPU is deployed strictly to execute a coarse-grained global search that prunes subspaces and isolates promising basins of attraction. This topographical decomposition drastically reduces the search volume, delegating only the localized, high-resolution refinement to classical GPUs. We demonstrate that by identifying basins of attraction through a variational quantum probe, we can substantially reduce the required number of classical iterations. Furthermore, we provide evidence for the conjecture that increasing qubit resolution (i.e., finer data decomposition and thus increasing the importance of cutting/knitting [17]) directly reduces the energy footprint of global search by replacing dissipative classical iterations with energy-preserving unitary quantum operations [18], [19]. This integration ensures that the convergence rate remains robust even as problem dimensions scale into regimes that have remained inaccessible to traditional (i.e., classical) optimization theory. Our contributions are as follows. • Topographical Preconditioning: We introduce the use of QPU as a preconditioner, identifying global basins of attraction to “warm-start” classical GPU optimizer. • Exponential Volume Reduction: We provide evidence that our framework reduces the search volume exponentially (see Table I). • One-level decomposition (dimension-wise): We partition the problem space into independent quantum sub-
registers (for separable functions) as opposed to particle (data) decomposition. • We show that our implementation outperforms the stateof-the-art HPC classical algorithm [1] in terms of the number of correct solution for both the Rastrigin and Ackley functions and provide evidence that this approach is promising for non-separable functions such as the Himmelblau function.
II. BACKGROUND A. Curse of Dimensionality The primary challenge in global mathematical optimization remains the “curse of dimensionality”, where the volume of the search space grows exponentially with each added dimension. In non-convex landscapes with multiple local minima, classical search algorithms must navigate a solution space that becomes increasingly sparse as the number of dimensions d increases. This scalability bottleneck is most evident in multimodal test functions like the Rastrigin function, which has 11d local minima within the standard range [−5.12, 5.12]. We consider the global optimization of a non-convex function f (x) where x ∈ Rd with d is the dimensionality. A canonical benchmark is the Rastrigin function:
f (x) = A × d +
d X 2 xi − A × cos(2πxi )
(1)
i=1
For A = 10, this function possesses a global minimum at x = 0, surrounded by a dense field of local minima. Assuming each dimension xi is bounded by an interval of length L, the total volume of the search space V increases exponentially as Ld . Consequently, for a constant particle density ρ, the number of required classical samples N = ρV grows exponentially. Experimental data shows that for d = 10, classical swarm success rates drop significantly (d ≥ 5) as particles fail to initialize within the global basin of attraction [1]. While modern GPU-accelerated methods like Z EUS can process thousands of starting points simultaneously using PSO and BFGS methods, they remain fundamentally constrained by the probability of landing in the basin of the global minimum. As shown in Figure 1, when the number of dimensions is increased from 2 to 10 while maintaining a constant number of particles (105 ), the count of correct solutions (Ncorrect ) follows a clear exponential decay. This rapid decay highlights that for a 10-dimensional problem (containing approximately 26 billion local minima), the number of particles that successfully locate the true global minimum basin is effectively zero. This exponential failure rate suggests that simply increasing classical computational resources or particle counts is an insufficient strategy for high-dimensional global search.
A. Quantum Approximate Optimization Algorithm QAOA is an iterative method designed primarily for combinatorial optimizations such as Max-Cut [21]–[23]. It applies a sequence of cost and mixing Hamiltonians to find the ground state of the problem. While theoretically sound for discrete structures, NISQ-era implementations of QAOA have shown only marginal improvements over classical baselines and are highly sensitive to parameter settings and the structural noise of the hardware [23]–[25]. Furthermore, QAOA is difficult to adapt to continuous high-dimensional optimization without significant discretization overhead that often leads to deep, decoherence-prone circuits [18], [26], [27]. D-QEO addresses this gap by utilizing a register-based dimensional mapping.
10000
Ncorrect
1000
100
10
1 2
3
4
5
6
7
8
9
10
Dimensionality
Fig. 1. Box and whisker plot showing performance degrades drastically for a classical solver for the Rastrigin function as the dimensionality of the problem increases when using the same number of particles
B. Non-differentiable Landscapes: The Ackley Function To test the resilience of D-QEO, we also evaluate a separable variant of the widely benchmarked Ackley function [20]: f (x) =
d X
(−20 exp (−0.2|xi |) − exp (cos(2πxi )) + 20 + e)
i=1
(2) Unlike the standard Ackley function that introduces coupling with a global radial ℓ2 -norm in the first term and an averaged exponential term, this modified variant is a sum of independent one-dimensional Ackley slices. This separability allows compatibility with current decoupled1 distributed circuit decomposition methodology for now, allowing d independent quantum circuits to process the landscape without requiring cross-circuit entanglement. The separable Ackley function preserves the challenging geometry of its native counterpart. It still features a deep central funnel leading to the global minimum, surrounded by a nearly flat outer plateau with many shallow local minima. We show that this topography is still difficult for classical exploration (see Figure 5 and Figure 6). Furthermore, the absolute value component |xi | makes the functions non-differentiable at the origin, which is where the global minimum is located. This point of non-differentiability restricts purely gradient-based classical approaches when they approach convergence. III. L ITERATURE R EVIEW A comprehensive review of existing quantum and quantuminspired algorithms reveals critical gaps that this paper aims to fill. We examine the Quantum Approximate Optimization Algorithms (QAOA), Variational Quantum Eigensolvers (VQE), and Quantum Particle Swarm Optimization (QPSO). 1 Future work will address weakly and partially coupled paradigms, outside the scope of this paper.
B. Variational Quantum Eigensolver VQE is a prominent primitive for near-term quantum advantage, particularly in quantum chemistry and molecular groundstate approximations [28]–[30]. It prepares a parameterized quantum state and measures the expectation value of a Hamiltonian, using a classical optimizer to tune the parameters. Although VQE displays resilience against certain noise types, its application to continuous mathematical function optimization remains largely under-explored compared to its use in electronic structure calculations [14], [29], [31]. Furthermore, existing studies rarely scale beyond 20-qubits [29], [30], [32]. D-QEO fills this gap by leveraging the mathematical separability of high-dimensional functions to enable a 50-qubit implementation through scalable circuit cutting [33], [34]. C. Quantum-Inspired and Hybrid PSOs Quantum-behaved PSO (QPSO) and other quantum-inspired variants simulate quantum dynamics on classical hardware to improve global search [35]–[37]. While these methods improve upon standard PSO [2] by introducing mechanisms such as probabilistic tunneling and improved population diversity mechanisms [38]–[40]. However, as they do not utilize true quantum superposition or entanglement, they remain fundamentally bound by classical computational scaling limits [41]– [43]. D-QEO bridges this by integrating a true QPU probe as a hardware-level preconditioner. The strategic gap identified across all these methods is the lack of a high-dimensional (50-qubit) continuous optimization framework that effectively combines quantum global exploration with classical high-precision refinement while managing hardware constraints through the 2-level decomposition. This gap directly motivates the formulation of our approach. IV. P ROBLEM F ORMULATION AND M OTIVATION A. Constructing the Discretized Hamiltonian: The General Case We translate the continuous classical problem into a discrete quantum space through binary encoding. A step-by-step mapping for the general case, using the non-separable Himmelblau function [44] as an example:
Step 1: Classical Discretization First, we must define the classical boundaries of our search space. Let us define a continuous variable x bounded by a minimum and maximum value: x ∈ [xmin , xmax ]. If we dedicate N qubits to represent this variable, we can encode 2N discrete states. We can represent any classical integer k in the range [0, 2N − 1] using standard binary notation: N −1 X k= qi 2i , qi ∈ {0, 1} (3) i=0
To map this integer k to our continuous range, we use a simple linear interpolation formula with step size ∆ = xmax −xmin : 2N −1 xk = xmin + k∆ (4) If we substitute our binary sum into this equation, we get the exact classical value for any bitstring: x(q) = xmin + ∆
N −1 X
qi 2i
(5)
i=0
Note that our explicit use of the σ̂ z notation alongside operator hats to prevent any notational collision with the PauliX operator, σ̂ x . This operator X̂ is a matrix whose eigenvalues correspond exactly to the discretized coordinates on our search grid and whose eigenvectors are the computational basis states (the bitstrings). Step 4: Constructing the Final Hamiltonian With the quantum representation for individual continuous variables established, we now substitute these operators into our target objective function to generate the final energy landscape. The Himmelblau function relies on two variables: f (x, y) = (x2 + y − 11)2 + (x + y 2 − 7)2
(8)
To evaluate this on a quantum computer, we assign N qubits to represent x, and a separate set of N qubits to represent y. We construct the operator X̂ for the first set of qubits, and the operator Ŷ for the second set. Finally, we substitute these operators directly into the Himmelblau polynomial to create our Hamiltonian, Ĥ: Ĥ = (X̂ 2 + Ŷ − 11I)2 + (X̂ + Ŷ 2 − 7I)2
(9)
Step 2: Formulating the binary number operator Next, we formulate an operator whose eigenvalues are classical bits (qi ∈ 0, 1). In quantum mechanics, the computational basis states |0⟩ and |1⟩ are the eigenstates of the Pauli-σ̂ z operator. The action of the σ̂ z operator on these states yields their corresponding eigenvalues: σ̂ z |0⟩ = |0⟩ (eigenvalue is + 1) σ̂ z |1⟩ = −|1⟩
(eigenvalue is − 1)
We require a mathematical transformation that maps the eigenvalues of the Pauli-σ̂ z operator, {1, −1}, to the classical bit values qi ∈ {0, 1}. We achieve this by defining the number operator, n̂i , with the following linear transformation: I − σ̂iz (6) 2 Where I is the identity matrix, and σ̂iz is the Pauli-Z operator applied to the i-th qubit. We can verify the eigenvalues of this operator: z • If the qubit is in state |0⟩, the σ̂ eigenvalue is 1. The formula yields: (1 − 1)/2 = 0. z • If the qubit is in state |1⟩, the σ̂ eigenvalue is −1. The formula yields: (1 − (−1))/2 = 1. We now have a well-defined quantum operator whose computational basis eigenstates perfectly map to the classical binary values. n̂i =
Step 3: Building the Variable Operator We substitute our number operator (n̂i ) back into Eq. 5. This gives us a new Quantum Operator, X̂, that represents our continuous variable in the Hilbert space: X̂ = xmin I + ∆
N −1 X i=0
2i n̂i
(7)
Expanding this polynomial mathematically yields a sum of tensor products of Z and I operators across all qubits. Because the Hamiltonian is constructed from σ̂ z and identity operators, it is not only composed of commuting terms but is also explicitly diagonal in the computational basis. Therefore, its gound state is guaranteed to be a single basis state, the bitstring representing the global minimum on our discretized search space. B. The Scaling Challenge of Non-separable Functions We evaluated our hybrid algorithm using an asymmetrical Himmelblau function, which features four mathematically identical global minima (E = 0) at different x and y locations. While the classical algorithm distributed its solutions relatively evenly across all four target basins, identifying Global Min 1 204 times, Min 2 280 times, Min 3 241 times, and Min 4 175 times, the 5 qubits per dimension hybrid algorithm exhibited severe discretization-induced bias, collapsing all 900 of its successful runs into Global Min 2. This behavior is not an algorithmic bug but a fundamental artifact of mapping a continuous search space onto a finitedimensional Hilbert space. By encoding the [−50, 50] continuous space using 5-qubits register per dimension, the system is limited to 25 = 32 orthogonal basis states. This discretization forces the quantum state vector to sample the objective function on a coarse grid with ≈ 3.2-unit intervals. Due to the alignment of this specific grid, the discrete coordinates nearest to Global Minima 1,3, and 4 land on the steep walls of those valleys, evaluating to high energy floors (≈ 50.9, ≈ 14.6, and ≈ 58.6, respectively). Consequently, the quantum optimizer misses three of the four mathematically identical global minima in the continuous space. The coarse discretization artificially breaks this equivalence, skewing the energy landscape and creating a single,
Algorithm
Classical
Hybrid (10−Qubit)
Hybrid (5−Qubit)
1000
900
Number of Successes
750
500
280 250
204
310
284 241
220
175 86
identifying the global basin of attraction. These identified regions serve as high-quality seed points for classical highresolution solvers. In this work, we restrict our scope strictly to separable functions—encompassing both differentiable (e.g., the Rastrigin function) and non-differentiable (e.g., the Ackley function) landscapes. Exploiting this separability enables us to scale the problem through circuit cutting (parallelization) while completely bypassing the exponential overhead associated with classical tensor knitting. Despite this structural constraint, our preliminary data (Fig. 5) demonstrates that this approach enhances existing classical PSO methods to a degree that, to the best of our knowledge, remains unachievable without the integration of QPUs.
(+ ,− )
A. Register-based Mapping
G
lo
ba lM
in
4
(− ,− ) G
lo
ba lM
in
3
(− ,+ ) 2 in ba lM lo G
G
lo
ba lM
in
1
(+ ,+ )
0
Identified Basin
Fig. 2. Motivation for cutting and knitting. Increasing from 5 to 10 qubits per dimension removes the grid-alignment bias and restores mathematical symmetry but incurs an exponential scaling cost in qubit requirements.
dominant global minimum at Global Min 2 (≈ 6.9). Visually demonstrated in Figure 3, this wave function collapses into a single discrete basin (bitstring 010110) as opposed to four corresponding global minima. To counter this discretization artifact, we increase the register size to 10 qubits per dimension (1024 states), drastically refining the grid resolution. As shown in Figure 2, providing the quantum phase with enough precision and thus distributing the runs much more naturally: 220 in Min 1, 310 in Min 2, 86 in Min 3, and 284 in Min 4, we were able to find all four global minimum. However, in general increasing the qubit register is not a scalable solution. For high-dimensional problems with certain level of entanglement, increasing the qubit count per dimension scales the state space exponentially, quickly exceeding the coherence and connectivity limits of near-term NISQ hardware as well as the memory limits of classical simulators. This scaling bottleneck directly motivates our future work in quantum circuit cutting and classical tensor network knitting. This exponential scaling bottleneck directly motivates our work in this paper. To achieve utility-scale quantum optimization (e.g., 50 qubits) on near-term infrastructure without the cost of classical tensor network knitting, we restrict our focus to separable landscapes, leveraging their mathematical structure to perfectly decouple the quantum registers. V. M ETHODOLOGY The Distributed Quantum-Enhanced Optimization (D-QEO) framework introduces an architectural role reversal in hybrid computing [36], [37]. Rather than relying on the quantum processor to pinpoint the exact global minimum, we deploy QPU as a topographical pre-processor. Its role is to construct a probability density map of the objective function’s landscape,
Unlike particle-based mappings common in classical metaheuristics, D-QEO maps the continuous space into a discrete Hilbert space using a register-based dimensional encoding. For a D-dimensional optimization problem, we assign a register of K qubits to each dimension xi . Building upon the quantum discretization formulated in Eq. 5, we define the mapping for each dimension xi using its respective K-qubit register: X̂i = xmin I + ∆
K−1 X
2k n̂i,k
(10)
k=0
where n̂i,k represents the number operator for the k-th qubit in the i-th register. For our benchmarks, we utilize K = 5, allowing a relatively small sub-circuit to explore 25 = 32 discrete sub-regions per dimension. B. Quantum Topographical Preconditioning with CVaR The core of our approach is the Quantum Topographical Preconditioner outlined in Algorithm 1. To train the variational quantum circuit, we utilize the Conditional Value-at-Risk (CVaR) as our objective function [45]. For a set of M = 1000 samples with energies {E1 , . . . EM } sorted in ascending order, we define the objective for a confidence level α = 0.1 [45] as: ⌈αM ⌉
CVaRα (θ) =
X 1 Ek ⌈αM ⌉
(11)
k=1
The classical optimizer iteratively updates the variational parameters to minimize this tail energy. As the optimization progresses, the quantum wave-function naturally “bunches” around the lowest energy eigenstates, which correspond to the most promising basins of attraction in the classical landscape. To prepare this state, we employ a Hardware-Efficient Ansatz (HEA) [45], [46] (and references within) that has an initial layer of Hadamard gates to create a uniform superposition, followed by parameterized Ry rotation layers (l = 3) and a cyclic CNOT entangling ring. The ansatz for a single dimension and K = 5 qubit register is illustrated in Figure 4.
1250
Measurement Count
1000
Energy
750
750
500
500 250
250
000000 000001 000010 000011 000100 000101 000110 000111 001000 001001 001010 001011 001101 001110 001111 010000 010001 010011 010100 010101 010110 010111 011000 011001 011010 011011 011100 011101 100000 100001 100010 100011 100100 100101 100110 100111 101000 101001 101010 101011 101100 101101 101110 101111 110000 110001 110010 110011 110100 110101 110110 110111 111000 111001 111010 111100 111101 111110 111111
0
Quantum Register State
Fig. 3. Quantum measurement distribution for the 2D Himmelblau function at low spatial resolution (K = 3 qubits per dimension). Although the landscape features four global minima (E = 0), the coarse 8 × 8 grid artificially breaks this energy degeneracy. The grid intersection closest to one specific minimum possesses a lower discrete energy than the intersections near the other three. The quantum optimizer correctly identifies this discrete anomaly, resulting in a pronounced wave function collapse into a single basin. Higher qubit resolutions are required to mitigate this grid-induced degeneracy breaking.
q0
H
Ry (θ0 )
Ry (θ5 )
...
q1
H
Ry (θ1 )
Ry (θ6 )
...
q2
H
Ry (θ2 )
Ry (θ7 )
...
q3
H
Ry (θ3 )
Ry (θ8 )
...
q4
H
Ry (θ4 )
Ry (θ9 )
...
ematical decomposition alongside circuit-level cutting using CUDA-Q [47]. Because the quantum fragments are perfectly decoupled, the memory requirements of each independent Kqubit subcircuit fit within the VRAM of a single GPU. This allows execution of the global evaluation by concurrently processing the subcircuit fragments on a single GPU. By doing so, we maximize classical simulation throughput and demonstrate that high-dimensional quantum optimization can be achieved without requiring massive supercomputing. D. Hybrid VQE and Classical Refinement
Fig. 4. Hardware-Efficient Ansatz utilized for a single K = 5 qubit dimension register. The circuit builds a parameterized probability distribution over the discrete search space.
C. Distributed Circuit Execution for Separable Landscapes A defining feature of this methodology is the exploitation of functional PD separability. For a separable objective function f (x) = i=1 fi (xi ), the corresponding global Hamiltonian decomposes perfectly into D independent PD partial Hamiltonians (Ĥi ) and thus the total is: Ĥtotal = i=1 Ĥi . Consequently, the global quantum state does not require cross-register entanglement. We can mathematically and operationally cut the D × K qubit circuit into D independent K-qubit fragments. This guarantees the distribution and the evaluation of each dimensional fragment asynchronously across multiple QPUs + GPUs (or just GPUs in our CUDA-Q simulation) with O(c) quantum knitting overhead, essentially bypassing the exponential complexity of high-dimensional quantum optimization. To simulate utility-scale environments (e.g., 50 qubits) in our experiments, we employ this math-
The transition from the quantum environment to the classical HPC solver requires translating the localized quantum wave-functions back into continuous floating-point variables. Since we cannot pass the raw measurement bitstrings directly to the classical optimizer, we extract the statistical properties of the optimized probability distribution to define a highly bounded, continuous search space. Once the variational training is complete, we perform a global concatenation to reconstruct the seed coordinates for classical refinement. We extract the absolute best coordinates found xbest , and calculate the weighted centroid of the samples in the CVaR tail, xcvar . The final seed coordinate is computed as a weighted average: xseed = β × xbest + (1 − β) × xcvar
(12)
Inspired by PSO [2], β serves as a hyperparameter, controlling the trust placed in the single best measurement versus the distributional centroid. Through experimentation, we picked β = 0.7. Future thorough investigation will be performed. To construct the bounding box for the classical solver, we define a dynamic search radius δ. Drawing on the principles of
step-size adaptation in Covariance Matrix Adaptation (CMAES) [48] and stochastic trust-region methods [49], δ scales with the root-mean-square (RMS) spatial deviation: δ = δbase + γ × RM S
(13)
We set the minimum algorithmic trust-region δbase = 0.5 to match the ±0.5 boundaries of our target global basins (§ VI-A1). By setting the scaling multiplier γ = 2, the expansion term acts as a 2σ spatial confidence interval. This dynamically expands the classical solver’s bounds strictly when the quantum sub-routines exhibit high spatial uncertainty. These bounded coordinates initialize the classical phase [1]. Algorithm 1 D-QEO: Quantum Topographical Preconditioner Require: Dimensions D, Qubits per dimension P K, Bounds [Xmin , Xmax ], separable objective f (x) = fi (xi ) Phase 1: Initialization & Distribution Initialize D independent quantum registers of size K Map each fi (xi ) to local Hamiltonian Ĥi using Eq. (10) Phase 2: Parallel Quantum Preconditioning for i = 1 to D do in parallel Initialize Ansatz |Ψi (θi )⟩ (Fig. 4) while convergence criteria not met do Sample 1000 shots from |Ψi (θi )⟩ once per iteration Calculate CVaR0.1 from lowest 100 energy samples Update θi using COBYLA to minimize CVaR end while Extract xbest,i and calculate tail centroid xcvar,i end for Phase 3: Reconstruction & Classical Refinement Construct Xseeds and calculate search radius δ with variance Decode B into continuous coordinates Xseeds if f (x) is differentiable then Warm-start PSO+BFGS at Xseeds within radius δ else Warm-start PSO around Xseeds within radius δ end if return Global Minimum found by classical solver [1] VI. E XPERIMENTS To evaluate the performance of the proposed Hybrid DQEO architecture, we benchmark our solver against classical baselines [1] using a set of highly non-convex, symmetric optimization functions. These landscapes are chosen for the dense population of local minima, which are designed to trap standard optimization algorithms. A. Experimental Setup and Hardware Configuration To ensure the reproducibility of our results, all quantum circuit simulations were executed using NVIDIA’s CUDA-Q framework. The simulations were performed in state-vector mode, accelerated by a single NVIDIA A100 Tensor Core GPU (80 GB). This hardware configuration allowed for the efficient classical simulation of our 50-qubit space by executing the independent 5-qubit chunks.
1) Evaluation Metrics and Statistical Trials: Because stochastic optimization algorithms exhibit variance, each experiment was repeated for 100 independent trials. The primary metric for success, Ncorrect , is defined as the number of trials that successfully converge to the true global minimum. A trial is classified as “correct” if the final optimized coordinates fall strictly within the basin of attraction of the global minimum (e.g., bounded within ±0.5 of the origin across all dimensions [1]). 2) Search Space Reduction Metrics: To rigorously quantify the topographical advantage of the quantum preconditioner, we track the continuous bounds generated by the quantum state. Across the 100 independent trials for a given dimension D, we extract the widest (worst-case) bounding box [lb, ub] that successfully captured the global basin. From this, QD we calculate the preconditioned search volume as Vpre = i=1 (ubi − lbi ). Furthermore, because landscapes like Rastrigin and Ackley are highly symmetric, with local minima occurring at approximate integer intervals, we calculate the total number of local minima within these reduced bounds. For a given dimension i, the number of trapped minima is determined by the number of integers bounded within [lbi , ubi ]. The total preconditioned minima count is the product of these surviving traps across all D dimensions (see Table I). B. Differentiable landscapes For the differentiable continuous benchmark, we test the algorithm on the N -dimensional Rastrigin function. Rastrigin is difficult for classical solvers due to its cosine-modulated landscape, which produces thousands of deep local minima surrounding a single global minimum. We demonstrate that utilizing the quantum phase as a preconditioner fundamentally alters the search dynamics. By collapsing the probability wave into the bounding box containing the true global minimum, the quantum phase bypasses the high-frequency local minima. As we increase the number of dimensions, the classical baseline severely drop in the number of correct solutions. In contrast, the hybrid algorithm maintains a higher percentage of correct solutions. To analyze the impact of the quantum optimizer’s computational budget on the final convergence variance, we employ the gradient-free COBYLA [50] optimizer to train the variational circuit during Phase 1. It is well-documented in optimization literature that the performance of gradient-free algorithms like COBYLA is heavily dependent on the available iteration budget, often requiring a substantial number of functional evaluations to yield high-quality results [50]–[52]. To quantify this trade-off between the depth of the quantum search and the stability of the final classical convergence, we systematically limit the quantum preconditioning phase to maximum functional budgets of Neval ∈ {200, 2000, 8000}. C. Non-differentiable landscapes In this setting, we explore the algorithmic trade-off of optimizing highly non-convex landscapes where the classical gradient is undefined. Purely classical optimization on
Classical
D−QEO (200 Evals)
D−QEO (2000 Evals)
D−QEO (8000 Evals)
100,000
Ncorrect
such functions requires gradient-free approximation or purely heuristic methods, which suffer from severe computational overhead and slow convergence rates in high dimensions. By utilizing our quantum preconditioner, which evaluates the landscape with discrete state sampling rather than relying on local gradients, we can effectively isolate the target basin without requiring objective function differentiability. We also demonstrate that feeding this tightly constrained bounding box to a classical gradient-free solver reduces the search volume exponentially. The following section presents the empirical outcomes of these trials, specifically focusing on dimensional scaling, computational effort, and topological volume reduction.
1,000
VII. R ESULTS A. Convergence Stability Across Dimensions We first evaluate the algorithm’s ability to locate the true global minimum basin as the dimensionality of the problem scales. Figure 5 illustrates the number of successful solutions (Ncorrect ) achieved by the classical baseline compared to the D-QEO on the Rastrigin function. As is consistent with established literature regarding the curse of dimensionality [1] (and references within), the purely classical algorithm exhibits an exponential decay in success rate, failing almost entirely as the search space reaches 10 dimensions. Conversely, the quantum preconditioner drastically stabilizes the swarm. By collapsing the probability wave into the bounding box containing the true global minimum, the hybrid algorithm maintains a remarkably high percentage of correct solutions regardless of the ambient dimensional scaling. Even at the restricted quantum evaluation budget of Neval = 200, the hybrid algorithm maintains a reliable success rate in 10 dimensions. Increasing the variational budget to Neval = 8000 further tightens the probability distribution of the quantum state, resulting in a near-perfect success rate across all evaluated dimensions. This behavior directly aligns with established classical findings that the solution quality of gradientfree algorithms like COBYLA scales heavily with the available evaluation budget [50]–[52]. However, our results demonstrate a key quantum advantage: even an “early,” truncated COBYLA run provides sufficient topographical data to the quantum state to effectively warm-start and rescue the classical solver. To demonstrate the versatility of this preconditioning approach, we also evaluated the algorithms on the separable Ackley function shown in Figure 6. The Ackley landscape is highly difficult for quasi-Newton methods because its nondifferentiability at the origin. Because the quantum preconditioner evaluates the landscape with discrete state sampling rather than local gradients, it effectively bypasses this nondifferentiability. As a result, the D-QEO warm-started classical solver consistently achieves a higher yield of correct solutions in high dimensions compared to the purely classical baseline. B. Reduction in Classical Computational Effort In addition to improving accuracy, D-QEO must demonstrate a tangible reduction in the classical workload. Figure 7
10
2
3
4
5
6
7
8
9
10
Dimensionality Fig. 5. Number of correct solutions (Ncorrect ) for the Rastrigin function. The D-QEO preconditioner preserves convergence in high dimensions where the classical baseline experiences exponential failure. The box plots display the distribution of successful runs. The central line indicates the median, the box edges represent the 25th and 75th percentiles (Interquartile Range), the whiskers extend to 1.5 × IQR, and the discrete dots represent outlier runs.
tracks the classical computational effort—measured in total BFGS iterations—required to achieve convergence. It is important to note that because these experiments are executed on classical GPU-based quantum simulators, we do not claim a reduction in total wall-clock execution time, as simulating 50 qubits inherently introduces massive classical overhead. Rather, these results demonstrate a fundamental reduction in algorithmic complexity and classical pathlength. By utilizing the QPU to map the topography and identify the global basin of attraction, the D-QEO framework effectively skips the classical exploration phase. Consequently, the classical BFGS optimizer requires fewer iterations to descend to the exact minimum compared to an un-preconditioned swarm. This validates the conjecture that increasing qubit resolution replaces dissipative, energyintensive classical iterations with energy-preserving quantum unitary operations, thereby fundamentally reducing the energy footprint required for complex global optimization. C. Volume Reduction Beyond classical iteration counts, the mathematical advantage of the D-QEO preconditioner is best understood through the reduction of the search space and the systematic elimination of local minima. Table I compares the original classical search space against the worst-case localized bounding boxes extracted from the quantum distribution. The results demonstrate staggering efficiency gains that fundamentally disrupt the “curse of dimensionality”. For instance,
Classical
D−QEO (200 Evals)
Classical
D−QEO (200 Evals)
D−QEO (2000 Evals)
D−QEO (8000 Evals)
D−QEO (2000 Evals)
D−QEO (8000 Evals)
100
Classical BFGS iterations
Ncorrect
100,000
1,000
10
30
10
3 2
3
4
5
6
7
8
9
10
Dimensionality Fig. 6. Number of correct solutions (Ncorrect ) for the separable Ackley function. The quantum sampling effectively bypasses the classical gradient failures caused by the function’s non-differentiability at the global optimum. Plotting conventions follow those defined in Figure 5.
in the 10-dimensional Rastrigin landscape, the ambient space spans [−5.12, 5.12]D , containing 1110 ≈ 2.59 × 1010 (over 25 billion) local minima. However, the quantum preconditioner bounds the search so tightly around the global basin that in the worst-case trial, 7 of the 10 dimensions were restricted to a width of less than 1.0, trapping only the true minimum coordinate (0). The remaining 3 dimensions exhibited slightly more variance, trapping two minima each. Consequently, the quantum warm-start collapsed a 25billion minima problem into a localized continuous subspace containing only 17 ×23 = 8 total local minima. By feeding this tightly constrained bounding box to the classical solver, the quantum phase mathematically eliminates the vast majority of the landscape’s non-convexity, trapping the classical optimizer within a continuous subspace where failure is highly improbable. Similar exponential reductions are observed in the nondifferentiable Ackley landscape, validating the framework’s robustness across varied topographies. VIII. D ISCUSSION AND F UTURE W ORK A. Algorithmic Complexity vs. Hardware Overhead The results demonstrate a clear reduction in the classical computational effort required to solve high-dimensional separable landscapes. However, it is critical to contextualize these findings within the current state of quantum computing. Currently, the total wall-clock execution time of the D-QEO framework is dominated by the overhead of classically simulating utility-scale quantum circuits.
2
3
4
5
6
7
8
9
10
Dimensionality Fig. 7. Classical optimization effort required to achieve convergence. The quantum topographical warm-start drastically reduces the number of classical BFGS iterations needed to traverse the landscape, indicating a shift from dissipative classical computation to efficient quantum preconditioning. Plotting conventions follow those defined in Figure 5.
Therefore, this study explicitly evaluates algorithmic complexity and classical iteration reduction rather than raw temporal speedup. The framework is designed with the foresight that as physical QPUs mature—specifically regarding gate execution times, active qubit reset capabilities, and QPUCPU network bandwidth—the physical wall-clock execution will naturally align with the theoretical algorithmic efficiency demonstrated in this study. B. NISQ Constraints and Hardware Noise It is also worth noting that the primary scaling analysis presented in this work is using CUDA-Q. By the time of the presentation, we anticipate that we will have data from a subset of IBM Nighthawk or Heron processors. Because D-QEO utilizes the QPU strictly as a coarsegrained topographical preconditioner rather than a highprecision solver, we hypothesize it inherently possesses a high tolerance for such errors. C. Extending to Non-Separable Landscapes While this work successfully circumvents the exponential spatial bottleneck by restricting its scope to separable functions, many high-value industrial [13] and basic research optimization problems are inherently non-separable, featuring heavily coupled variables (as demonstrated by the Himmelblau scaling failure in Section IV-B). Future work will focus on extending the D-QEO architecture to navigate these highly coupled landscapes. To achieve this, we plan to explore advanced dynamic circuit cutting and
TABLE I T OPOGRAPHICAL P RECONDITIONING I MPACT: C OMPARING THE ORIGINAL CLASSICAL SEARCH SPACE AGAINST THE CONSTRAINED BOUNDING BOXES ([lb, ub]) EXTRACTED FROM THE QUANTUM DISTRIBUTION . T HE REDUCTION FACTOR ILLUSTRATES THE EXPONENTIAL ADVANTAGE OF ISOLATING THE GLOBAL BASIN . Landscape
D
Orig. Vol. (Vorig )
Precond. Vol. (Vpre )
Reduction Factor
Orig. Minima
Precond. Minima
Rastrigin
2 3 4 5 6 7 8 9 10
104.86 1, 073.74 1.10 × 104 1.13 × 105 1.15 × 106 1.18 × 107 1.21 × 108 1.24 × 109 1.27 × 1010
2.80 6.90 11.08 19.09 41.75 108.64 76.60 228.67 178.41
37.46 155.53 992.10 5, 896.43 2.76 × 104 1.09 × 105 1.58 × 106 5.41 × 106 7.11 × 107
121 1, 331 1.46 × 104 1.61 × 105 1.77 × 106 1.95 × 107 2.14 × 108 2.36 × 109 2.59 × 1010
2 8 4 8 8 64 8 16 8
Ackley
2 3 4 5 6 7 8 9 10
4, 294.97 2.81 × 105 1.84 × 107 1.21 × 109 7.92 × 1010 5.19 × 1012 3.40 × 1014 2.23 × 1016 1.46 × 1018
28.15 162.13 3, 398.70 5, 394.95 2.44 × 104 1.39 × 105 5.73 × 105 4.73 × 106 3.63 × 107
152.57 1, 736.10 5, 427.60 2.24 × 105 3.24 × 106 3.75 × 107 5.94 × 108 4.71 × 109 4.03 × 1010
4, 225 2.75 × 105 1.79 × 107 1.16 × 109 7.54 × 1010 4.90 × 1012 3.19 × 1014 2.07 × 1016 1.35 × 1018
30 150 3, 584 5, 400 1.88 × 104 1.62 × 105 9.33 × 105 4.86 × 106 4.03 × 107
classical tensor network knitting techniques. By strategically severing the cross-register entanglement operations, we aim to reconstruct the high-resolution probability amplitudes to break the symmetry of complex, non-separable functions without exceeding the strict coherence and connectivity limits of NISQera architectures. Furthermore, future iterations will prioritize the integration of Quantum Error Mitigation (QEM) protocols to facilitate direct execution on physical quantum hardware. IX. C ONCLUSION In this paper, we presented the Distributed QuantumEnhanced Optimization framework, a novel hybrid architecture designed to mitigate the “curse of dimensionality” in continuous global search. By strategically reversing the traditional roles of hybrid computing, we deployed the QPU as a topographical preconditioner to construct probability density maps of the objective landscape, leaving high-precision continuous refinement to classical GPU solvers. To achieve utility-scale execution, we exploited the mathematical structure of separable continuous functions, allowing a monolithic 50-qubit search space to be mathematically severed and distributed into independent subcircuits without incurring exponential tensorknitting overhead. Our experiments on the highly non-convex Rastrigin and non-differentiable Ackley functions demonstrate that D-QEO successfully neutralizes the exponential failure rates inherent to purely classical swarms. By collapsing the quantum wave function into the bounding box of the true global minimum, the quantum warm-start dramatically reduces the required classical pathlength and BFGS iteration count. Ultimately, this framework provides a highly scalable, distributed blueprint for translating near-term quantum resources into tangible computational reductions for high-dimensional optimization.
ACKNOWLEDGMENT This research was supported in part by the Fermi National Accelerator Laboratory (FERMILAB-CONF-26-0212CSAID) graduate Summer Internship Program (DS & MP), the Office of Naval Research (ONR) through the U.S. Naval Research Laboratory (JS), and the Richard T. Cheng Endowment (NC) at ODU. The authors thank Min Dong at ITS in ODU for his help with CUDA-Q at ODU. The authors would like to acknowledge the use of Google’s Gemini AI during the preparation of this manuscript. In accordance with IEEE policy, we note that this generative AI system was utilized strictly as an advanced copyediting and formatting assistant. All scientific concepts, experimental data, algorithms, and core analytical arguments remain the entirely original, humangenerated work of the authors. This work was performed using computational facilities at ODU enabled by grants from the National Science Foundation (MRI grant no. CNS-1828593) and Virginia’s Commonwealth Technology Research Fund and Google Cloud Platform through ODU’s Monarch Sphere initiative. Any subjective views or opinions expressed in this paper do not necessarily represent the views of the national labs, the NSF, or the United States Government. R EFERENCES [1] D. Soós, M. Paterno, D. Ranjan, and M. Zubair, “Zeus: An efficient gpu optimization method integrating pso, bfgs, and automatic differentiation,” in 2025 IEEE 32nd International Conference on High Performance Computing, Data, and Analytics (HiPC). IEEE, 2025, pp. 225–235. [2] J. Kennedy and R. Eberhart, “Particle swarm optimization,” in Proceedings of ICNN’95-international conference on neural networks, vol. 4. ieee, 1995, pp. 1942–1948. [3] C. G. Broyden, “The convergence of a class of double-rank minimization algorithms 1. general considerations,” IMA Journal of Applied Mathematics, vol. 6, no. 1, pp. 76–90, 1970. [4] R. Fletcher, “A new approach to variable metric algorithms,” The computer journal, vol. 13, no. 3, pp. 317–322, 1970.
[5] D. Goldfarb, “A family of variable-metric methods derived by variational means,” Mathematics of computation, vol. 24, no. 109, pp. 23–26, 1970. [6] D. F. Shanno, “Conditioning of quasi-newton methods for function minimization,” Mathematics of computation, vol. 24, no. 111, pp. 647– 656, 1970. [7] L. B. Rall and G. F. Corliss, “An introduction to automatic differentiation,” Computational Differentiation: Techniques, Applications, and Tools, vol. 89, pp. 1–18, 1996. [8] L. A. Rastrigin, “Systems of extremal control,” Nauka, 1974. [9] J. Bian, “The nova experiment: overview and status,” arXiv preprint arXiv:1309.7898, 2013. [10] M. A. Acero et al., “Improved measurement of neutrino oscillation parameters by the nova experiment,” Phys. Rev. D, vol. 106, p. 032004, 2022. [11] ——, “Monte Carlo method for constructing confidence intervals with unconstrained and constrained nuisance parameters in the NOvA experiment,” JINST, vol. 20, no. 02, p. T02001, 2025. [12] N. Buchanan, S. Calvez, D. Doyle, V. Hewes, A. Himmel, J. Kowalkowski, A. Norman, M. Paterno, T. Peterka, S. Sehrish, A. Sousa, T. Thakore, and O. Yildiz, “Analyzing nova neutrino data with the perlmutter supercomputer,” Poster presented at the International Conference for High Performance Computing, Networking, Storage and Analysis (SC22), Dallas, TX, USA, Nov. 2022, poster. [Online]. Available: https://sc22.supercomputing.org/presentation/index-610.htm [13] B. Abi, R. Acciarri, M. A. Acero, G. Adamov, D. Adams, M. Adinolfi, Z. Ahmad, J. Ahmed, T. Alion, S. A. Monsalve et al., “Volume i. introduction to dune,” Journal of instrumentation, vol. 15, no. 08, pp. T08 008–T08 008, 2020. [14] M. Haidar, O. Adjoua, S. Badreddine, A. Peruzzo, and J.-P. Piquemal, “Non-iterative disentangled unitary coupled-cluster based on liealgebraic structure,” Quantum Science and Technology, vol. 10, no. 2, p. 025031, 2025. [15] A. Maciejunes, J. Stenger, D. Gunlycke, and N. Chrisochoides, “Solving large-scale vehicle routing problems with hybrid quantum-classical decomposition,” 2025. [Online]. Available: https://arxiv.org/abs/2507.05373 [16] E. Billias and N. Chrisochoides, “Towards a utility-scale quantum edge detection for real-world medical image data,” 2025. [Online]. Available: https://arxiv.org/abs/2507.10939 [17] W. Tang, “Enabling large-scale quantum computing via distributed and hybrid architectures,” Ph.D. dissertation, Princeton University, 2025. [18] K. Morimoto, Y. Takase, K. Mitarai, and K. Fujii, “Continuous optimization by quantum adaptive distribution search,” Physical Review Research, vol. 6, no. 2, p. 023191, 2024. [19] T. Lubowe and S. Morino, “Best-in-class quantum circuit simulation at scale with nvidia cuquantum appliance,” 2022. [20] D. Ackley, A connectionist machine for genetic hillclimbing. Springer science & business media, 2012. [21] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028, 2014. [22] G. G. Guerreschi and A. Y. Matsuura, “Qaoa for max-cut requires hundreds of qubits for quantum speed-up,” Scientific reports, vol. 9, no. 1, p. 6903, 2019. [23] R. Shaydulin and Y. Alexeev, “Evaluating quantum approximate optimization algorithm: A case study,” in 2019 tenth international green and sustainable computing conference (IGSC). IEEE, 2019, pp. 1–6. [24] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, “Barren plateaus in quantum neural network training landscapes,” Nature communications, vol. 9, no. 1, p. 4812, 2018. [25] M. V. S. Cerezo de la Roca, A. T. Arrasmith, P. J. Czarnik, L. Cincio, and P. J. Coles, “Effect of barren plateaus on gradient-free optimization,” Quantum, vol. 5, no. LA-UR–20-29699, 2021. [26] G. Verdon, J. M. Arrazola, K. Brádler, and N. Killoran, “A quantum approximate optimization algorithm for continuous problems,” arXiv preprint arXiv:1902.00409, 2019. [27] M. Luna, V. Patare, G. Aksoy, and G. Cattan, “Implementation of a quantum approximate optimization algorithm for continuous variables with qiskit,” HAL, 2025. [28] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature communications, vol. 5, no. 1, p. 4213, 2014. [29] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth et al., “The variational quantum
eigensolver: a review of methods and best practices,” Physics Reports, vol. 986, pp. 1–128, 2022. [30] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” nature, vol. 549, no. 7671, pp. 242–246, 2017. [31] G. Intoccia, U. Chirico, V. Schiano Di Cola, G. P. Pepe, and S. Cuomo, “Quantum adaptive search: a hybrid quantum-classical algorithm for global optimization of multivariate functions,” Frontiers in Applied Mathematics and Statistics, vol. 11, p. 1662682, 2025. [32] D. Gunlycke, C. S. Hellberg, and J. P. Stenger, “Cascaded variational quantum eigensolver algorithm,” Physical Review Research, vol. 6, no. 1, p. 013238, 2024. [33] S. Bravyi, G. Smith, and J. A. Smolin, “Trading classical and quantum computational resources,” Physical Review X, vol. 6, no. 2, p. 021043, 2016. [34] T. Peng, A. W. Harrow, M. Ozols, and X. Wu, “Simulating large quantum circuits on a small quantum computer,” Physical review letters, vol. 125, no. 15, p. 150504, 2020. [35] J. Sun, B. Feng, and W. Xu, “Particle swarm optimization with particles having quantum behavior,” in Proceedings of the 2004 congress on evolutionary computation (IEEE Cat. No. 04TH8753), vol. 1. IEEE, 2004, pp. 325–331. [36] W. Fang, J. Sun, Y. Ding, X. Wu, and W. Xu, “A review of quantumbehaved particle swarm optimization,” IETE Technical Review, vol. 27, no. 4, pp. 336–348, 2010. [37] J. Sun, W. Fang, X. Wu, V. Palade, and W. Xu, “Quantum-behaved particle swarm optimization: analysis of individual particle behavior and parameter selection,” Evolutionary computation, vol. 20, no. 3, pp. 349– 393, 2012. [38] S. M. Mikki and A. A. Kishk, “Quantum particle swarm optimization for electromagnetics,” IEEE transactions on antennas and propagation, vol. 54, no. 10, pp. 2764–2775, 2006. [39] L. d. S. Coelho, “Gaussian quantum-behaved particle swarm optimization approaches for constrained engineering design problems,” Expert Systems with Applications: An International Journal, vol. 37, no. 2, pp. 1676–1683, 2010. [40] Y. Fan, D. Tian, Q. Xu, J. Sun, Q. Xu, and Z. Shi, “Particle swarm optimization based on k-means clustering and adaptive dual-groups strategy,” Swarm and Evolutionary Computation, vol. 100, p. 102226, 2026. [41] S. Li, M. Tan, I. W. Tsang, and J. T.-Y. Kwok, “A hybrid psobfgs strategy for global optimization of multimodal functions,” IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), vol. 41, no. 4, pp. 1003–1014, 2011. [42] J.-J. Liang and P. N. Suganthan, “Dynamic multi-swarm particle swarm optimizer,” in Proceedings 2005 IEEE Swarm Intelligence Symposium, 2005. SIS 2005. IEEE, 2005, pp. 124–129. [43] O. U. Rehman, S. Yang, and S. U. Khan, “A modified quantumbased particle swarm optimization for engineering inverse problem,” COMPEL-The international journal for computation and mathematics in electrical and electronic engineering, vol. 36, no. 1, pp. 168–187, 2017. [44] D. Himmelblau, “Reduction methods in nonlinear programming,” 1982. [45] P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, and S. Woerner, “Improving variational quantum optimization using cvar,” Quantum, vol. 4, p. 256, 2020. [46] J. Stenger, C. S. Hellberg, and D. Gunlycke, “Hybrid vqe-cvqe algorithm using diabatic state preparation,” arXiv preprint arXiv:2512.04801, 2025. [47] J.-S. Kim, A. McCaskey, B. Heim, M. Modani, S. Stanwyck, and T. Costa, “Cuda quantum: The platform for integrated quantum-classical computing,” in 2023 60th ACM/IEEE Design Automation Conference (DAC). IEEE, 2023, pp. 1–4. [48] N. Hansen, “The cma evolution strategy: A tutorial,” arXiv preprint arXiv:1604.00772, 2016. [49] A. R. Conn, N. I. Gould, and P. L. Toint, Trust region methods. SIAM, 2000. [50] M. J. Powell, “A direct search optimization method that models the objective and constraint functions by linear interpolation,” in Advances in optimization and numerical analysis. Springer, 1994, pp. 51–67. [51] S. Bagheri, W. Konen, M. Emmerich, and T. Bäck, “Self-adjusting parameter control for surrogate-assisted constrained optimization under limited budgets,” Applied Soft Computing, vol. 61, pp. 377–393, 2017.
[52] L. M. Rios and N. V. Sahinidis, “Derivative-free optimization: a review of algorithms and comparison of software implementations,” Journal of Global Optimization, vol. 56, no. 3, pp. 1247–1293, 2013.