QHDOPT: A Software for Nonlinear Optimization with Quantum Hamiltonian Descent
Samuel Kushnir Jiaqi Leng Yuxiang Peng Lei Fan Xiaodi Wu



1] Samuel Kushnir
2,3,4,†,[Samuel Kushnir and Jiaqi Leng contributed equally to this work. Most of the work was completed at the University of Maryland.]] Jiaqi Leng
1,3] Yuxiang Peng
5,6] Lei Fan
1,3] Xiaodi Wu
[1]Department of Computer Science, University of Maryland
[2]Department of Mathematics, University of Maryland
[3]Joint Center for Quantum Information and Computer Science, University of Maryland
[4]Department of Mathematics and Simons Institute for the Theory of Computing, University of California, Berkeley
[5]Department of Engineering Technology, University of Houston
[6]Department of Electrical and Computer Engineering, University of Houston
[†]mailto:jiaqil@terpmail.umd.edujiaqil@terpmail.umd.edu


QHDOPT: A Software for Nonlinear Optimization with Quantum Hamiltonian Descent
    [
    
==============================================================================



§ ABSTRACT 
We develop an open-source, end-to-end software (named QHDOPT), which can solve nonlinear optimization problems using the quantum Hamiltonian descent (QHD) algorithm. 
QHDOPT offers an accessible interface and automatically maps tasks to various supported quantum backends (i.e., quantum hardware machines). 
These features enable users, even those without prior knowledge or experience in quantum computing, to utilize the power of existing quantum devices for nonlinear and nonconvex optimization tasks.
In its intermediate compilation layer, QHDOPT employs SimuQ, an efficient interface for Hamiltonian-oriented programming, to facilitate multiple algorithmic specifications and ensure compatible cross-hardware deployment.
The detailed documentation of QHDOPT is available at <https://github.com/jiaqileng/QHDOPT>.




§ INTRODUCTION

Nonlinear optimization, also known as nonlinear programming, is a branch of mathematical optimization concerned with solving problems in which the objective function, constraints, or both, exhibit nonlinearity. While nonlinear optimization problems are common in various application fields such as engineering, management, economics, and finance, these problems are in general nonconvex with complicated landscape features like multiple local stationary points, valleys, and plateaus. As the number of variables grows, the complexity of the problem could grow rapidly, posing a significant challenge in obtaining globally optimal solutions.

Several open-source and commercial software packages, including Ipopt <cit.>, Gurobi <cit.>, and CPLEX <cit.>, have been developed to tackle large-scale nonlinear optimization problems. While these optimizers can incorporate powerful heuristics to enhance performance for certain problem instances, there is no polynomial-time guarantee for these optimizers because nonlinear optimization is generally -hard. Often, the problem structure is unknown, and there is no commonly agreed-upon go-to optimizer for nonlinear optimization in practice.

Quantum computers are emerging technologies that can leverage the laws of quantum mechanics to offer theoretical and practical advantages over classical computers in solving large-scale computational problems. Unlike their classical counterparts, quantum computers utilize a unique phenomenon known as quantum tunneling to accelerate the solution of nonconvex optimization problems. Specifically, a quantum particle can pass through a high potential barrier that would be insurmountable classically due to insufficient energy. This exotic behavior enables a quantum computer to bypass sub-optimal solutions, efficiently navigating the complex landscape of nonlinear optimization.

Recently, <cit.> proposes a novel quantum algorithm named Quantum Hamiltonian Descent (QHD). QHD is inspired by the observation that many first-order (i.e., gradient-based) methods can be interpreted as dynamical processes governed by physical laws. For example, it has been shown that the celebrated Nesterov's accelerated gradient descent algorithm can be modeled by a time-dependent Lagrangian mechanical system that would find local minima in the system <cit.>. 
By upgrading the classical Lagrangian mechanics to quantum mechanics, we end up with a minimum-finding quantum process, just like gradient descent. Additionally, this quantum dynamical process demonstrates the quantum tunneling effect, making it a competitive candidate for solving nonconvex optimization problems.
Simulating this quantum dynamical process on a quantum computer gives rise to QHD, a simple but powerful quantum algorithm for continuous optimization, especially nonlinear problems with nonconvex objective functions.
A follow-up work by <cit.> shows that QHD can solve a family of hard optimization instances in polynomial time, while an empirical study suggests that these problem instances are intractable for many classical optimization algorithms such as branch-and-bound, stochastic gradient descent, interior point method, etc.

A key feature of QHD is that it is formulated as a quantum evolution, which can be simulated on both digital and analog quantum computers. This feature allows us to implement QHD to tackle real-world tasks with near-term realizable quantum computers.
Digital quantum computers perform computation by applying a sequence of elementary quantum gates to an initial quantum state. These machines exhibit provable quantum advantages over classical (digital) computers for certain computational tasks, however, they require a large number of digital (i.e., error-corrected) qubits. Although there has recently been a groundbreaking experimental demonstration of early fault tolerance <cit.>, existing digital quantum computers have not yet reached the size necessary to accelerate the solution of real-world problems in application domains such as management, finance, and engineering <cit.>.
Analog quantum computers solve computational tasks by configuring and emulating a real quantum system and then performing quantum measurements. These devices are easier to fabricate, control, and scale <cit.>, while they are unavoidably noisy, and no general error correction technique is currently practical <cit.>.
<cit.> proposed a systematic technique named Hamming encoding that enables us to implement QHD to solve quadratic programming (QP) problems on analog quantum computers with Ising Hamiltonian.
This technique is exemplified in solving 75-dimensional nonconvex QP problems, where the noisy real-machine implementation of QHD outperforms existing open-source nonlinear optimization software like Ipopt.

In this paper, we develop , an end-to-end implementation of QHD for nonlinear optimization.
A notable feature of  is that it supports the deployment of QHD to multiple quantum computing hardware, including gate-based quantum computers such as IonQ, and analog quantum computers such as D-Wave.
 provides a user-friendly interface, with which a nonlinear optimization problem can be specified via either matrix/numeric or symbolic description.
Then, the implementation of QHD is fully automatized and the (approximate) optimal solutions will be returned once the computation is completed.
The mid-level compilation and cross-hardware deployment are achieved by utilizing  for Hamiltonian-oriented programming <cit.>.



  
Organization.
The rest of the paper is organized as follows.[This paper is not intended to be a comprehensive tutorial or documentation on . Instead, we direct the readers to <https://github.com/jiaqileng/QHDOPT> for the source code, examples, tutorials, and documentation.]
In problem-formulation, we explain the general problem formulation for nonlinear optimization problems that can be processed and solved by .
In workflow-main, we discuss the workflow of   including the quantum backend and classical refinement.
In hop, we discuss several unique design features of , especially the multi-backend compatibility achieved by incorporating the Hamiltonian-oriented programming (HOP) framework.
In qhd, we briefly review the QHD algorithm and its implementation on both digital and analog quantum computers.
Then, in workflow, we sketch the workflow of the software, including all major steps in the implementation of QHD and classical post-processing. example provides two worked examples of modeling and solving nonlinear optimization problems.
In state, we review the current state and trend of quantum optimization software.
We conclude with a comparison of  with other available open-source optimizers in comparison.





 §.§ Problem formulation: box-constrained nonlinear optimization

The package  solves nonlinear programming problems of the following form:


    
    min_x     f(x_1, …, x_n) = ∑^n_i=1  g_i(x_i)_univariate part +∑^m_j=1 p_j(x_k_j)q_j(x_ℓ_j)_bivariate part,
    
        s.t.       L_i≤ x_i ≤ U_i, ∀ i ∈{1,…,n},


where x_1,…,x_n are n variables subject to the box constraint x_i ∈ [L_i,U_i] ⊂ℝ for each i = 1,…,n, and the indices k_j, ℓ_j ∈{1,…,n} and k_j ≠ℓ_j for each j = 1,…, m. The functions g_i(x_i), p_j(x_k_j), and q_j(x_ℓ_j) are real univariate differentiable functions defined on ℝ. Note that the univariate part in obj has at most n terms because we can always combine separate univariate functions of a fixed variable x_i into a single one. However, there is no upper bound for the integer m (i.e., the number of bivariate terms).[It is generally impossible to combine a sum of products into a single product form. For example, we can not find two univariate functions p(x) and q(y) such that p(x)q(y) = sin(x)y + xe^y.]

The nonlinear optimization problem primal is in general NP-hard <cit.> and can be used to model several common classes of optimization problems, including linear programming, quadratic programming, and polynomial optimization (with box constraints).
In the following examples, we show how to formulate some standard nonlinear optimization problems in the form of obj.

[Box-constrained Quadratic programming]
    A quadratic programming problem with a box constraint takes the form:
    
        
    min_x    f(x) 1/2x Qx + b x 
    
        s.t.      0 ≤ x ≤ 1,

    
where Q∈^n× n is a symmetric matrix and b is a real-valued vector of dimension n. The objective function can be written as 

    f(x) = ∑^n_i = 1(1/2Q_i,i x^2_i + b_i x_i) + ∑^n_1 ≤ k < ℓ≤ n Q_k, ℓx_k x_ℓ.

This function is represented by obj by choosing 


    g_i(x_i) = 1/2Q_i,i x_i^2 + b_i x_i,  ∀ i = 1,…, n, 
       p_j(x_k_j) = Q_k_j,ℓ_j x_k_j,    q_j(x_ℓ_j) = x_ℓ_j,   ∀ j ∈{1,…, n(n-1)/2}.


Here (k_j, ℓ_j) are the j-th pair in the enumeration {(k, ℓ): 1≤ k<ℓ≤ n}.




While the problem formulation can only handle box constraints, we note that many optimization problems with more sophisticated constraints can be reformulated in the form of primal by adding the constraints as a penalty term in the objective function.

[Spherical constraints]
    Consider the optimization problem with n variables:
    
        
    min_x     f(x) ∑^n_j=1α_j x_j
    
                s.t.     ∑^n_j=1x^2_j = 1,

    
    where α_j are real scalars for all j = 1,…,n. The feasible set of this problem is the n-dimensional sphere with radius 1, which can not be directly recast as a box in the form of box. Meanwhile, we observe that all the variables must take values between 0 and 1 because the unit sphere is contained in the unit (hyper-)cube. Therefore, we can reformulate spherical to a box-constrained optimization problem by the penalty method:
    
        
    min_x     f(x) ∑^n_j=1α_j x_j + λ(∑^n_j=1x^2_j - 1)^2,
    
               s.t.      0 ≤ x ≤ 1.

    
    This new problem can be handled by our software  since the objective function penalty-obj only involves uni- and bi-variate monomials. As the penalty coefficient λ > 0 grows, we can show that the solution to the box-constrained problem penalty will eventually converge to the optimal solution to the original problem spherical. 




It is worth noting that the problem formulation supported by  is restrictive, and there exist many general nonlinear optimization problems that cannot be directly expressed in primal. For example, our formulation cannot deal with objective functions involving trivariate monomials (e.g., xyz). 
While, in theory, QHD can handle box-constrained optimization models given access to ideal quantum hardware, in  we limit the appearance of trivariate parts or higher to cater to the current quantum hardware restrictions.

Additionally, we note that there may be several corner cases that are representable by primal but would require an excessively long time for  to parse and solve. For example, when the objective function involves thousands of bivariate functions, it might take  minutes to compile and implement the automatic differentiation subroutine based on JAX. We advise users to prioritize the use cases with low-degree polynomials, bounded exponential functions, and simple trigonometric functions.





 §.§ Solving problems in 

 utilizes the Quantum Hamiltonian Descent algorithm to facilitate the solution of nonlinear and nonconvex optimization problems.
Theoretically, Quantum Hamiltonian Descent, when running with an ideal fault-tolerant quantum computer, can solve many optimization problems up to global optimality given sufficiently long runtime <cit.>.
However, at the current stage, due to the lack of fault tolerance, we can only implement Quantum Hamiltonian Descent in a low-precision and noisy manner, which significantly reduces the solution quality promised by the theoretical guarantee.
To mitigate the noisy performance of near-term quantum hardware with limited resources, we adopt a hybrid quantum-classical computing workflow in  to achieve optimal performance, as illustrated in workflowB.



  
Pre-processing and problem encoding.
First, we map a box-constrained nonlinear optimization problem to a quantum-mechanical system with finite degrees of freedom. This reduced quantum model can be regarded as a finite-precision approximation of the original QHD model.
Then, this quantum model is embedded into a larger quantum system that is natively executable using one of the supported quantum backends. This process is called Hamiltonian programming. While the quantum hardware only “sees” a reduced version of the original problem, the Hamiltonian embedding technique <cit.> ensures that the spatial structure inherited from the original problem is preserved and naturally encoded in the quantum operator. Therefore,  allows us to run a coarse-grained version of QHD on near-term quantum devices.



  
Deployment and decoding.
Then, the quantum operator that encodes the original nonlinear optimization problem is constructed and executed on a quantum backend. Currently,  supports three backends: the D-Wave quantum computer, the IonQ quantum computer, and a classical simulator based on QuTiP. The measurement results from quantum devices are in 0-1 format (i.e., binaries), which requires a decoder to recover the corresponding solution in the continuous space (e.g., the unit box). 



  
Classical refinement.
Limited by the size and coherent time of current quantum devices, the quantum-generated solutions are of low precision and intrinsically noisy.  relies on classical local search algorithms, such as first- and second-order methods, to improve numerical precision. Currently,  supports two local optimizers: a general-purpose interior point method (Ipopt) and a truncated Newton method (TNC) implemented in SciPy. While we do not include other local search subroutines in , we note that generic local optimizers allowing box constraints should work as well.



Since  leverages local search algorithms as refiners, the output solutions are necessarily locally optimal (i.e., first- or second-order stationary points, depending on the choice of refinement subroutine).
That being said, we would like to note that the Quantum Hamiltonian Descent (QHD) algorithm, when executed on a large fault-tolerant quantum computer, is able to find the global minimum for a large family of nonconvex functions with mild assumptions, provided that the runtime is sufficiently long <cit.>. The performance of  for practical problems, however, heavily depends on the quality of near-term quantum devices, which are often of limited scale and prone to physical noise. Meanwhile, it is also possible to refine the quantum-generated solutions using a global solver (e.g., Gurobi, BARON); in this case, the global optimality is guaranteed, but the post-processing time could be significantly longer. Due to the limited time frame, we leave a global-solver-based refinement as future work.





 §.§ Unique design features




In what follows, we discuss a few unique design features of our software.



  
Hamiltonian-oriented programming (HOP).
 exploits the Quantum Hamiltonian Descent (QHD) algorithm to solve nonlinear optimization problems. This quantum algorithm is formulated as a Hamiltonian simulation (i.e., simulating the evolution of a quantum-mechanical system), encompassing a novel abstraction of computation on quantum devices, which we call Hamiltonian-Oriented Programming (HOP).
In contrast to the conventional circuit-based quantum computation paradigm where theorists describe quantum algorithms in terms of quantum circuits, the HOP paradigm describes quantum algorithms as a single or a sequence of quantum Hamiltonian evolution. 
This new paradigm enables us to build a stack of quantum applications by leveraging the native programmability of quantum hardware in the development of quantum algorithms and software, as illustrated in workflowA. 
The HOP paradigm is empowered by , a recent framework for programming and compiling quantum Hamiltonian systems by <cit.>. 
In , the programming and simulation of quantum Hamiltonian systems are wrapped in user-friendly Python methods. This makes the high-level programming and deployment of Hamiltonian-oriented quantum algorithms accessible to users with little exposure to real-machine engineering and manipulation.
A detailed discussion on the Hamiltonian programming and compilation in  is available in workflow.




  
Multi-backend compatibility.
In , we utilize  as an intermediate layer for the programming of QHD and leverage the  compiler to realize multi-backend compatibility. Through ,  initially constructs a hardware-agnostic Hamiltonian representation of QHD (i.e., Hamiltonian embedding) that can be deployed on various quantum backends, including D-Wave devices, IonQ devices, and classical simulators via QuTiP <cit.>. 




  
Automatic differentiation.
 relies on JAX, a high-performance numerical computing library, to perform automatic differentiation of smooth, nonlinear objective functions. This feature enables  to seamlessly post-process quantum-generated solutions using local search optimizers.





§ QUANTUM HAMILTONIAN DESCENT

In our software, we utilize QHD to solve box-constrained nonlinear optimization problems as described in primal. QHD solves a continuous optimization problem by simulating a quantum dynamical system governed by an evolutionary partial differential equation called Schrödinger equation. Here, we give a high-level review of this quantum algorithm and more details can be found in <cit.>. 



 §.§ Mathematical formulation and interpretation

Consider a nonlinear objective function f(x) with a box constraint Ω= {(x_1,…,x_n)∈^n L_i ≤ x_i ≤ U_i, ∀ i = 1,…, n}. To solve this optimization problem, QHD requires simulating the following Schrödinger equation over the feasible set Ω with Dirichlet boundary condition, i.e., Ψ(t,x) = 0 for x ∈∂Ω,

    i ∂/∂ tΨ(t,x) = [e^φ_t(-1/2Δ) + e^χ_tf(x)]Ψ(t,x),

subject to an initial state Ψ(t,x) = Ψ_0(x). Here, the operator Δ∑^n_i=1∂^2/∂ x^2_i is the Laplacian operator defined in the interior of Ω, and the time-dependent functions e^φ_t and e^χ_t control the total energy distribution of the quantum system. In practice, the initial state Ψ_0(x) is often chosen as a quantum state that is easy to prepare, for example, a Gaussian state or a uniformly random state. For general (nonconvex) optimization problems, it is observed that an inverse polynomially decaying e^φ_t and polynomially increasing e^χ_t (e.g., φ_t = -log(1+γ t^2), χ_t = log(1+γ t^2) with a positive γ) work well for many test problems <cit.>. With a Gaussian initial state and smooth time-dependent functions, the dynamics generated by qhd-pde can be simulated using 𝒪(nT) elementary gates and 𝒪(T) queries to the objective function f <cit.>.


Physically, the equation qhd-pde describes the time evolution of a quantum particle in the box Ω. The time-dependent functions e^φ_t and e^χ_t control the total energy distribution of this quantum particle: when their ratio e^φ_t/e^χ_t is large, the kinetic energy dominates and the particle tends to bounce around; otherwise, the potential energy takes over and the particle tends to stay still. If we choose these functions such that lim_t→∞ e^φ_t/e^χ_t = 0, the kinetic energy of the system is dissipated over time and eventually the quantum particle will take a low-energy configuration. At this point, if we measure this quantum particle, the measured position (which must lie in the feasible set Ω) is likely to give an approximate solution to the problem f(x). In some sense, QHD can be regarded as a quantum version of Polyak's heavy ball method <cit.>.

QHD describes a quantum particle exploring the optimization landscape f(x). When a high-energy barrier emerges, the quantum particle may leverage the quantum tunneling effect to go through the barrier and find a lower local minimum. However, simulating the quantum evolution qhd-pde with a classical computer would require exponential space and time, making this idea impractical as a classical optimization algorithm. 
On the other hand, the evolution qhd-pde can be efficiently simulated using a quantum computer, which makes QHD a genuine quantum algorithm that can leverage the quantum tunneling effect for nonconvex optimization.
Theoretically, it has been shown that QHD can efficiently find the global minimum for certain nonconvex problems with exponentially many local minima, while many classical optimizers such as simulated annealing and SGD appear to require a much longer time to obtain a global solution <cit.>. Numerical experiments also show that QHD outperforms classical first- and second-order methods in a broad class of nonconvex problems with many local stationary points <cit.>.



 §.§ Real-machine implementation


Quantum Hamiltonian Descent is formulated as a Hamiltonian simulation task, i.e., solving a quantum Schrödinger equation as in qhd-pde. While efficient quantum algorithms, such as those proposed by <cit.>, can tackle this simulation task exponentially faster than any known classical algorithms, these quantum simulation algorithms require large fault-tolerant quantum computers. Such ideal quantum computing hardware has not yet been realized due to the immature progress of quantum technology.

To fully exploit the limited programmability of current quantum hardware such as D-Wave and IonQ,  employs a technique named Hamiltonian embedding <cit.> to implement QHD. This technique enables us to map the QHD Hamiltonian to a larger Hamiltonian, and the latter can be natively simulated on existing quantum devices. This real-machine implementation technique is detailed in ham-program-compile.




§ THE WORKFLOW OF 




 §.§ Modeling of nonlinear problems




 offers support for two Python-based input formats: the SymPy format for symbolic input and the QP format for numerical input (i.e., arrays). These two input formats enable users to define their target optimization problems both efficiently and with great flexibility.





  
SymPy format.  <cit.> is a Python package that supports symbolic expression processes. Users can specify f(x) in primal by declaring variables in  and constructing the expression, as in the code snippet in sympy-exp. Here, we import necessary functions like exp from  and QHD from package  in lines 1 and 2. We declare variables x and y in  using symbols in line 3 where the passed string is for  to print the expressions. In line 4, we construct the function f(x), where y**1.5 represents the exponential y^1.5, exp(4*x) represents e^4x, and so on. Lastly, we create a QHD model instance in line 5 and pass f and a symbol list [x, y] to it, informing the QHD model the target optimization function is f with symbols x and y.



  
QP format. For users with specific interests in quadratic programming (QP), we provide a more efficient input model for them. To specify a QP instance with objective function f(x) = 1/2x Q x + b x, we can directly input the matrices Q and b, as in the code snippet in qp-exp. First, we construct Q by a nested Python list or a NumPy array in lines 2 and 3. It is required that Q forms a symmetric square matrix. Then we input the vector b as b. Similar to , we construct the instance by calling the QP method from  and pass Q, b into it. 




 §.§ Hamiltonian programming and compilation

Once a nonlinear optimization problem f(x) is defined using one of the supported input formats,  will form a Hamiltonian description of the corresponding QHD algorithm, as described in qhd-pde. This Hamiltonian description serves as an intermediate layer in the compilation stack and is independent of the choice of the backend (i.e., hardware-agnostic). Although  automates this process, making manual execution unnecessary in most cases, we provide detailed discussions for readers who are interested in gaining a deeper understanding of our software's design.

There are two major steps in the construction of the Hamiltonian description of QHD, namely, spatial discretization and Hamiltonian embedding.



  §.§.§ Spatial discretization

First, we need to perform spatial discretization of the QHD Hamiltonian (which is an unbounded operator) so that it can be described by a finite-dimensional quantum system. For a thorough and mathematically rigorous discussion, readers are encouraged to refer to <cit.>. Given a nonlinear optimization in the form of primal, the QHD Hamiltonian reads the following,

    H(t) = e^φ_t(-1/2Δ) + e^χ_t(∑^n_i=1  g_i(x_i) +∑^m_j=1 p_j(x_k_j)q_j(x_ℓ_j)),

which acts on any L^2-integrable functions over the feasible set Ω = {(x_1,…,x_n)∈^n L_i ≤ x_i ≤U_i, ∀ i = 1,…, n}. Here, for simplicity, we assume the feasible set is the unit box, i.e., L_i = 0 and U_i = 1 for all i = 1,…,n.
We utilize the centered finite difference scheme to discretize this differential operator. Suppose that we divide each dimension of the unit box Ω using N quadrature points {0, h,…, (N-2)h, 1} (where h = 1/(N-1)), the resulting discretized QHD Hamiltonian is an N^n-dimensional operator of the form,

    Ĥ(t) = e^φ_t(-1/2L_d) + e^χ_tF_d,

where (assuming k_j < ℓ_j for all j = 1,…,m)

    L_d = ∑^n_i=1 I⊗…⊗L_the i-th operator⊗… I,


    F_d = ∑^n_i=1 I⊗…⊗D(g_i)_the i-th operator⊗… I + ∑^m_j=1 I⊗…⊗D(p_j)_the k_j-th operator⊗…⊗D(q_j)_the ℓ_j-th operator⊗… I.

Here, I is the N-dimensional identity operator, L and D(g) are N-dimensional matrices given by (g is a differentiable function defined on [0,1] and g_i  g(ih) for i = 0,…,N-1),

    L = 1/h^2[  -2   1        ;   1  -2   1    ; ... ... ... ...;   1  -2   1;       1  -2;  ],   D(g) = [   g_0                  ;         g_1            ;   ...   ...   ...   ...;       g_N-2      ;             g_N-1;  ].

The tridiagonal L matrix corresponds to the finite difference discretization of the second-order differential operator d^2/d x^2, and the diagonal matrix D(g) corresponds to the finite difference discretization of the univariate function g(x). Note that L has a global phase -2/h^2, i.e., L = L' - 2/h^2 with L' only contains the off-diagonal part of L. Since the global phase does not affect the quantum evolution (therefore, the result of the QHD algorithm), we replace L with L' in the rest of the discussion.



  §.§.§ Hamiltonian embedding

The discretized QHD Hamiltonian, as described in discretized_qhd, is a Hermitian matrix with an explicit tensor product decomposition structure. This particular structure allows us to leverage the Hamiltonian embedding technique <cit.> to construct a surrogate Hamiltonian H(t) such that the QHD algorithm (i.e., simulating the Hamiltonian Ĥ(t)) can be executed by simulating H(t). In our case, the surrogate Hamiltonian H(t) is an Ising-type quantum Hamiltonian that involves at most nN qubits and max(n,m) N two-body interaction terms. This means H(t) can be efficiently simulated on current quantum computers, including IonQ's trapped ion systems and D-Wave's quantum annealer.

To construct the Hamiltonian embedding of Ĥ(t), the first step is to build the Hamiltonian embeddings of the N-by-N matrices L' and D(g) (for arbitrary differentiable function g). Both are sparse matrices so we can utilize the embedding schemes provided in <cit.>.  allows users to choose from three embedding schemes: Hamming[The details of Hamming embedding can be found in <cit.>. Note that this embedding scheme is referred to as “Hamming encoding” in <cit.>.], unary, and one-hot[More precisely, the one-hot embedding we implemented in  is referred to as “penalty-free one-hot” embedding in <cit.>.]. In embedding-scheme, we list the details of these embedding schemes when applied to L' and D(g). 
We note that the Hamming embedding scheme only works for quadratic programming, while the other two schemes (unary, one-hot) work for a broader class of nonlinear functions such as exponential functions.
To be consistent with our source code, we adopt the left-to-right 0-indexing system for bits/qubits, e.g., 1_0 0_1 1_2 1_3.



In embedding-scheme, the integer r represents the number of qubits used to embed an N-dimensional matrix. The operators 𝐗_k, 𝐘_k, and 𝐧_k are the Pauli-X, Pauli-Y, and number operator acting at site k, respectively,

    𝐗 = [ 0 1; 1 0 ],  𝐘 = [  0 -i;  i  0 ],  𝐧 = [ 0 0; 0 1 ].

Since the Hamming embedding scheme is only allowed for quadratic programming, we do not consider the Hamming embedding for general nonlinear functions g. Instead, we only consider the embedding of the identity and quadratic functions, i.e., g(x) = x and g(x) = x^2, their corresponding Hamming embeddings are ℰ_1 = 1/r∑^r-1_k=0𝐧_k and ℰ_2 = (ℰ_1)^2, respectively.

Now, we denote ℰ^(i)[A] as a Hamiltonian embedding of an N-by-N Hermitian matrix A acting on sites (i-1)r,(i-1)r+1,…,ir-1, where i = 1,…,n. Then, using the rules of building Hamiltonian embeddings <cit.>, we obtain a nr-qubit Hamiltonian that embeds the discretized QHD Hamiltonian Ĥ(t),

    H(t) = e^φ_t(-1/2∑^n_i=1ℰ^(i)[L']) + e^χ_t(∑^n_i=1ℰ^(i)[D(g_i)] + ∑^m_j=1ℰ^(k_j)[D(p_j)]ℰ^(ℓ_j)[D(q_j)]).


[One-hot embedding for x_1x_2]
    We give a simple example for the Hamiltonian embedding of the discretized QHD Hamiltonian when the objective function is f(x_1,x_2) = x_1x_2. This objective only involves a single bivariate term with p(x) = q(x) = x. We use the one-hot embedding with N = r = 3. Then, the Hamiltonian embeddings of L' and D(x) are,
    
    ℰ[L'] = 1/2h^2(𝐗_0𝐗_1 + 𝐗_1𝐗_2 + 𝐘_0𝐘_1 + 𝐘_1𝐘_2),  ℰ[D(x)] = 𝐧_0 + 1/2𝐧_1,

    respectively. As a result, the full Hamiltonian embedding reads
    
    H(t) =    2e^φ_t/h^2(𝐗_0𝐗_1 + 𝐗_1𝐗_2 + 𝐗_3𝐗_4 + 𝐗_4𝐗_5 + 𝐘_0𝐘_1 + 𝐘_1𝐘_2 + 𝐘_3𝐘_4 + 𝐘_4𝐘_5) + 
       e^χ_t(𝐧_0 + 1/2𝐧_1)(𝐧_3 + 1/2𝐧_4).



In , we use  to construct the Hamiltonian embedding H(t). The users only need to specify the number of qubits r (for each continuous variable), the embedding scheme, and a desired backend in the model.optimize() function, as detailed in the next subsection.




 §.§ Deployment and post-processing


When the Hamiltonian embedding H(t) of a given problem is built, it can be executed on a supported quantum backend by running optimize() (refer to example for sample code). The quantum measurement results are then retrieved from the executing backend in the form of bitstrings. Following this,  implements a series of classical post-processing subroutines. These include decoding the raw measurement results (i.e., bitstrings) into low-resolution solutions and refining them via a classical local solver. The refined solutions are then returned to the users as final results.




  §.§.§ Deployment on quantum devices

Currently,  supports three backend devices for deployment, including classical simulators (e.g., QuTiP), IonQ, and D-Wave. For all three backend devices, the quantum register is initialized to the uniform superposition state. On the IonQ device, the uniform superposition state can be prepared using a single layer of Hadamard gate; on the D-Wave device, the uniform superposition state is the default initial state and it can be prepared in microseconds. When deployed on IonQ,  uses φ_t = -log(1+γ t^2) and χ_t = log(1+γ t^2) for the time-dependent functions (see math-formulation for details).
The time-dependent functions on D-Wave are more restricted and they can only be specified as piece-wise linear functions. We find the default annealing schedule (20 microseconds) provided by the D-Wave device usually works well in practice. We also showcase user-specified time-dependent functions (annealing schedules) in a notebook in the “examples” folder.

Here, we demonstrate the deployment procedures in  using the D-Wave backend, while the same process applies to the other two backends. 
In dwave-execute, a snippet of the source code for the function QHD.dwave_exec() is displayed. After programming the Hamiltonian embedding and the quantum system realizing QHD (lines 2-3), we initiate an abstract D-Wave machine (line 5). Then  employ  to compile the Hamiltonian embedding into low-level device instructions readable by D-Wave (line 7) using 's DWaveProvider(), effectively generating Hamiltonian H_dev(t) on the D-Wave devices. Next, the instructions are sent to the D-Wave quantum computer to execute (line 9), and the raw quantum samples are collected by  as bitstrings (line 12-17).







  §.§.§ Decoding

As we have seen, the real-machine results are in the bitstring format because they are retrieved by computational basis measurements in the quantum computer.
These bitstrings need to be converted to floating-point arrays via the built-in decoder, as presented in decoder. This decoder maps a bitstring to a floating-point array that represents a low-resolution solution to the input optimization problem. For example, if we use the unary embedding for a 2-dimensional problem with resolution parameter r = 4, the decoder will map an 8-bit string to a length-2 array. For example, 00010011 is mapped to [0.25, 0.5]. More details of the embedding schemes and their decoding are available in <cit.>.





  §.§.§ Refinement

Limited by the size of current quantum devices, in most cases, we can only use a small resolution parameter (e.g., r =8) in the real-machine implementation of QHD. Therefore, the retrieved measurement results are merely low-resolution solutions to the specified optimization problem. To improve the precision of the solutions,  then post-processes the measurement results using local search classical optimization methods. 

In principle, any generic local optimizers allowing box constraints should work as well; due to the limited resources, we provide two classical refinement options for the users in , including the truncated Newton method (TNC) using SciPy and the interior point method using Ipopt. TNC is a quasi-Newton method, and the interior point method exploited by Ipopt is a second-order method. 
These classical refiners require the gradient and/or Hessian information of the objective functions. For quadratic programming problems (using the QP input format), the gradient and Hessian can be computed explicitly:

    ∇ f(x) = Q x,   H f(x) = Q.

For more general nonlinear optimization problems specified using the SymPy input format,  computes the gradient and Hessian information by employing Jax <cit.>, a high-performance numerical computing library developed by Google, to perform auto-differentiation.

The post-processed results can be retrieved from model.post_processed_samples. By default, the post-processing subroutine is enabled and automatically executed by running optimize(). However, users can also disable post-processing by specifying optimize(fine_tune=False). In this case, model.post_processed_samples returns None type. 





§ EXAMPLES USING 


In this section, we exhibit two simple examples showcasing the use cases of .



 §.§ Quadratic programming


We first consider a 2-dimensional quadratic programming problem, whose objective function is defined as follows,

    f(x,y) = -x^2 + xy - 1/2y^2 + 3/4x - 1/4y =1/2[ x y ][ -2  1;  1 -1 ][ x; y ]+[  3/4 -1/4 ][ x; y ],

for x, y ∈ [0,1].

In example-1, we exemplify using  to solve this QP problem with all three backends. Note that the API key (not included in ) is required to access cloud-based quantum computers such as D-Wave and IonQ.
By default, QHD.optimize() automatically executes the Scipy TNC method to fine-tune the raw quantum measurement data. To switch to the Ipopt optimizer in the fine-tuning step, one can specify post_processing_method="IPOPT" in the setup, as shown in line 17.







 §.§ Nonlinear optimization involving exponential function


Next, we consider the following nonlinear minimization problem with the objective function:

    f(x,y) = y^3/2 - e^4x(y-3/4),   x,y∈ [0,1].

This objective function f(x,y) is not a polynomial; it involves a fractional power and an exponential function.
In example-2, we illustrate a sample code that runs  to solve the problem defined above. The function is constructed using the SymPy input format, then deployed on the D-Wave quantum computer with resolution parameter r=8.
We may set verbose=1 to print a detailed summary of this execution, including the best-so-far coarse and fine-tuned solutions, as well as a total runtime breakdown, as shown in example2-readout.



§ THE STATE OF SOFTWARE FOR QUANTUM OPTIMIZATION


Software packages are crucial for lowering the barrier to developing and implementing quantum programs across broad user communities. Upon examining the current landscape of quantum software for mathematical optimization, we observe that the majority of the software dedicated to quantum optimization focuses on addressing combinatorial and discrete optimization problems, with limited options available for continuous optimization.

Generally, a combinatorial optimization problem can be reformulated as a Quadratically Unconstrained Binary Optimization (QUBO) problem, the solution of which is believed to be a promising application of quantum computing <cit.>.
There is a rich collection of libraries for quantum and quantum-inspired optimization that can be employed to generate QUBO reformulations, including  <cit.>,  <cit.>,  <cit.>,  <cit.>. 
These QUBO problems can be tackled by several methods, such as quantum annealing, Quantum Approximate Optimization Algorithms (QAOA) <cit.>, and other hybrid approaches <cit.>. 
D-Wave's  <cit.> enables users to interface with their direct QPU (i.e., quantum annealer) and hybrid solvers and retrieve results.
QuEra's  <cit.> is a high-level language for configuring programmable Rydberg atom arrays, which can be used to implement annealing-type quantum algorithms and discrete optimization problems like QUBO <cit.>.
Los Alamos Advanced Network Science Initiative has also released a package named  <cit.> for simulation and execution of quantum annealing. 
Besides, several software packages have been published for programming quantum circuits, including IBM's  <cit.>, Google's  <cit.>, Amazon's  <cit.>, Microsoft's  <cit.>, and Xanadu's  <cit.>. These tools can be used to deploy QAOA on gate-based quantum computers.

While there have been a few proposals for solving continuous optimization problems using quantum or hybrid computing devices such as photonic quantum computers <cit.> and coherent continuous variable machines <cit.>, we are not aware of a software library customized for nonlinear continuous optimization problems. In practice, some nonlinear optimization problems, such as quadratic programming, may be reformulated as QUBO problems and handled by the aforementioned software tools. However, it remains unclear whether this approach could lead to robust quantum advantages.




§ COMPARISON WITH EXISTING TOOLS


As we discussed in workflow,  first obtains some low-resolution solutions by executing the QHD algorithm for a nonlinear optimization problem through Hamiltonian embedding.
Next, the software employs a classical local search strategy for fast post-processing of the raw quantum results. It is of interest to understand to what extent the quantum component (i.e., the noisy implementation of QHD) improves the overall performance of . To this end, we have designed a benchmark test to evaluate the performance of  for nonlinear and nonconvex optimization problems.



 §.§ Test problems

We demonstrate the performance of  using fifteen randomly generated nonlinear optimization instances, all with unit box constraints.
Problem instances 1 - 5 are nonlinear programming (NLP) problems involving two or three continuous variables, as detailed in test-info.
Problem instances 6 - 10 are quadratic programming (QP) problems drawn from the benchmark devised in <cit.>.
Problems instances 11 - 15 are nonlinear programming (NLP) problems involving exponential functions, as specified in the following expression:

    f(x) = 1/2∑^N_i=1∑^N_j=1Q_i,je^x_ie^x_j + ∑^N_i=1 b_i e^-x_i.

The last ten test instances (6 - 15) are intermediate-scale problems with 50 continuous variables, ranging from 0 to 1. To ensure the successful mapping of these test problems to quantum computers with limited connectivity, these test problems are generated in a way such that their Hessians are sparse matrices.[Detailed expressions of these test instances are provided in the software repository, see the “examples” folder.]



These instances were generated in a largely random manner, and each possesses multiple local solutions, making them fairly challenging for classical optimization software. In our experiment, we observed that local solvers, such as Ipopt, cannot find globally optimal (or even approximately optimal) solutions unless a large number of random initial guesses are tried. BARON, a highly optimized commercial solver for global optimization, can find globally optimal solutions to sparse quadratic programming problems in under 1 second, but it takes several minutes to certify global optimality for nonlinear programming problems that involve exponential-type objectives.




 §.§ Experiment setup and results

In this subsection, we discuss the basic setup of the experiment and the numerical results. We test  with two different post-processing optimizers (i.e., Ipopt and Scipy-TNC) using the randomly generated nonlinear programming instances discussed in the previous section. As a comparison, we also run the two classical optimizers on the same test instances using uniformly random initialization.
These two classical optimizers are assessed as baselines to illustrate the quantum advantage introduced by the D-Wave-implemented Quantum Hamiltonian Descent (QHD).
The classical components in both experiments, including the decoding of D-Wave samples and classical refinement, were executed on a 2022 MacBook Pro laptop with an Apple M2 chip.
Our findings assert that QHD, when implemented with D-Wave, brings a significant advantage compared to the standalone use of classical optimizers.



  §.§.§ Experiment setup for 

We evaluate  on this benchmark using the D-Wave Advantage_system6.3 as the quantum backend. For a fair comparison, we use the unary embedding scheme for all instances, including quadratic and non-quadratic problems.
The anneal time is set to be the default value, i.e., 20 microseconds. 
The total quantum runtime per shot (see the “QPU” columns in benchmark) is calculated as the arithmetic mean of the “qpu_access_time” reported by D-Wave, which includes the programming, state preparation, annealing, and decoding. Note that we do not include the transmission time and the task queuing time in our report.
We test  with two post-processing optimizers (i.e., Ipopt and Scipy-TNC), and the classical post-processing (more precisely, classical refinement) time is reported in the “Classical Refine” columns in benchmark. The standard deviation of the classical refinement time is also reported in the parenthesis. Except for the initial guesses, both solvers use the default parameters as provided with their Python API.



  §.§.§ Baseline using classical optimizers

As a comparison, we also test three classical optimizers for the same set of problems: Ipopt, Scipy-TNC, and BARON. The first two optimizers are initialized with 1000 uniformly random guesses in the unit box [0,1]^d (where d is the problem dimension), and the runtime data for the 1000 runs have been collected. For a fair comparison, we use the same random seeds for both methods. BARON is executed to generate global solutions with a 2-minute timeout. Except for the initial guesses, all solvers use the default parameters as provided with their Python API. control-group shows the runtime of the three classical solvers: for Ipopt and TNC, the arithmetic mean (and standard deviation) of the runtime is reported; for BARON, the total runtime is reported. 

Note that for the last five test instances, BARON failed to certify the global optimality of the obtained solutions within the 2-minute timeout window. In baron-sol-quality, we further investigate BARON's solution quality. Our results suggest that, while BARON can find solutions as good as those from the other tested solvers in a comparable timescale, a much longer time is required to prove the global optimality of the obtained solutions. Therefore, we regard the solution returned by BARON as the global minimum.





  §.§.§ Performance metric

In the experiments, we use time-to-solution (TTS) as the key metric to evaluate the performance of various optimization methods. 
TTS is defined as the total runtime required by a method to achieve at least 0.99 success probability. It can be calculated using the following formula,

    TTS = t_0 ×⌈ln(1-0.99)/ln(1-p_s)⌉,

where t_0 is the (average) runtime per shot, and p_s is the success probability.
For all the 15 test instances, the time-to-solution data of four methods (QHD+Ipopt, QHD+TNC, Ipopt, and TNC) are presented in tts-summary.
For the experiments involving QHD (i.e., the results in benchmark), t_0 is calculated as the sum of average QPU time and average classical refinement time; for the experiments that only involve classical optimizers (i.e., the results in control-group), t_0 is equivalent to the (average) classical runtime. 
The success probability p_s is estimated by the fraction of “successful” events in the 1000 samples/trials. Here, a result x' is considered successful if the optimality gap f(x') -f(x^*) is less than 0.001, where x^* is the solution obtained by BARON.












 §.§ Interpretation of the experiment results


Based on the time-to-solution data as reported in tts-summary, we observed that  (QHD + a classical optimizer) always outperforms the standalone use of a classical optimizer for the 15 randomly generated test instances. As for the two local optimizers, we find that Scipy-TNC works better than Ipopt as a post-processing subroutine. While the two optimizers usually return refined samples with comparable success probability, Scipy-TNC always shows a lower TTS due to a notably faster runtime. This is potentially because TNC is a quasi-Newton method that does not need to solve the full Newton linear system in the iterations.


Another interesting finding is that the classical refinement times of quantum-generated samples (see “Classical Refine” in benchmark) in  are significantly shorter than the average runtime of the direct use of local optimizers (see “Avg. Runtime” in control-group). For example, in test instance 1, the average post-processing time (using Ipopt) for a quantum-generated initial guess is 0.2 s, while the average runtime of Ipopt given uniformly random guesses is 13 s.
To further investigate this phenomenon, we plot the distribution of objective function values corresponding to three different sample groups, including (1) randomly generated initial guesses, (2) quantum (D-Wave) generated samples, and (3) TNC refined samples (using quantum-generated samples as initial guesses), as shown in solution-quality.[In solution-quality, we only plot objective function values for high-dimensional problems, i.e., test instances 6 - 15.] 
While the quantum-generated solutions are limited by low precision, it is observed that they are still qualitatively better than random initial guesses. In the subsequent post-processing,  performs a local search subroutine to refine solution quality by improving numerical accuracy. In other words, the quantum sampler in  can be regarded as a fast and efficient warm-start that devises initial guesses of better quality.






§ CONCLUSION AND FUTURE WORK

 is the first open-source software leveraging quantum devices for nonconvex nonlinear optimization problems, providing an accessible interface for domain experts without quantum computing knowledge. Exploiting the idea of Hamiltonian-oriented programming, it efficiently uses quantum devices by implementing the Quantum Hamiltonian Descent algorithm with the  framework. We demonstrated 's effectiveness through examples and benchmarks, showing its advantage over classical solvers, especially in large, complex instances. However, the current limitations of quantum device programmability and scalability constrain our benchmarks' scale. While QHD shows promise in solving complex optimization problems, further empirical studies are needed for real-world performance evaluation.

There are several avenues for future development of . First, it is desired to broaden the problem class that can be handled by .
Currently, due to hardware limitations,  supports only the optimization of box-constrained nonlinear problems defined as a sum of univariate and bivariate functions. We anticipate that, shortly,  can be extended to address more complicated objective functions as quantum technology and quantum algorithm design continue to co-evolve.
Second, while local search algorithms work well to improve the precision of quantum-generated samples, to obtain a global optimality guarantee, it might be promising to replace the refinement/post-processing subroutine in  with global optimizers (for example, those based on branch-and-bound).
Third, with further progress in quantum engineering,  is expected to support more quantum devices from different platforms, including commercial or laboratory devices, which is essential to understanding the advantage of  given different combinations of embedding schemes and quantum devices. Last but not least,
 can be expanded into a plugin for various domain-specific tools, including those in engineering, management, finance, and economics. Adaptions to specific domains are invaluable for users to better utilize quantum devices for their domain problems.  
Our overarching goal is to establish a user-friendly tool, empowering individuals and organizations to harness the power of quantum devices to solve challenging problems in the real world.



§ ACKNOWLEDGMENT

This work was partially funded by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research, Accelerated Research in Quantum Computing under Award Number DE-SC0020273, the Air
Force Office of Scientific Research under Grant No. FA95502110051, the U.S. National Science Foundation grant CCF-1816695, CCF-1942837 (CAREER), ECCS-2045978, a Sloan research fellowship, the Simons Quantum Postdoctoral Fellowship, and a Simons Investigator award through Grant No. 825053. J.L. and Y.P. are also supported by an open-source quantum software grant from the Unitary Fund. 




plainnat


