Quantum optimization algorithm based on multistep quantum computation

We present a quantum algorithm for finding the minimum of a function based on multistep quantum computation and apply it for optimization problems with continuous variables, in which the variables of the problem are discretized to form the state space of the problem. Usually the cost for solving the problem increases dramatically with the size of the problem. In this algorithm, the dimension of the search space of the problem can be reduced exponentially step by step. We construct a sequence of Hamiltonians such that the search space of a Hamiltonian is nested in that of the previous one. By applying a multistep quantum computation process, the optimal vector is finally located in a small state space and can be determined efficiently. One of the most difficult problems in optimization is that a trial vector is trapped in a deep local minimum while the global minimum is missed, this problem can be alleviated in our algorithm and the runtime is proportional to the number of the steps of the algorithm, provided certain conditions are satisfied. We have tested the algorithm for some continuous test functions.


I. INTRODUCTION
Optimization problem is one of the most important problems in science and engineering.
It includes a wide class of problems ranging from molecular modeling, quantum mechanical calculations, machine learning, to combinatorial optimization.These problems can be classified into different categories, e.g., continuous or discrete optimization, constrained or unconstrained optimization, convex or nonconvex optimization, differentiable or nondifferentiable optimization, deterministic or stochastic optimization [1][2][3][4][5], etc.There is no universal optimization algorithm.Most classical optimization algorithms start with a trial vector that is varied by using different techniques to find the optimum of an objective function.The cost of the algorithms can become very expensive due to the increase of the dimension of the state space of the problem, which is known as "the curse of dimension".Another problem that often happens for optimization algorithms is that the trial vector is trapped in a deep local minimum, while missing the global minimum of the objective function.
Optimization has also been studied in the framework of quantum computation.Adiabatic quantum computing (AQC) is designed for solving combinatorial optimization problems [6], in which starting with the ground state of a simple initial Hamiltonian, the system is evolved adiabatically to a final Hamiltonian whose ground state encodes the solution to the optimization problem.Despite the theoretical guarantee of the adiabatic theorem, the condition of adiabaticity in AQC is difficult to maintain in practice, since the allowed rate of evolution is determined by the minimum energy gap between the ground and the first excited states of the adiabatic evolution Hamiltonian, which is not known a priori.Quantum annealing is a heuristic quantum optimization algorithm [7][8][9][10][11] that can be viewed as a relaxation of AQC, where the conditions of adiabaticity are not met and the evolution time from an initial Hamiltonian to the final Hamiltonian is determined heuristically.Whether or not quantum annealing can provide quantum speed-up over classical heuristic algorithms is still not clear.Variational quantum algorithms such as quantum approximate optimization algorithm (QAOA) [12] are hybrid quantum-classical algorithms designed for near-term noisy intermediate-scale quantum computers [13] without performance guarantees.It is known that in the infinite depth limit, the QAOA recovers adiabatic evolution and would converge to the optimal solution.The gradient decent methods are used for optimization problems with continuous variables.The methods find local minima of a smooth function by moving along the direction of the steepest descent.Quantum algorithm provides an efficient way in calculating numerical gradients [14], and has been used in iterative algorithms for polynomial optimization [15].Optimization algorithms based on gradient decent require that the objective function to be smooth, and they have the problem of being trapped in a local minimum and missing the global minimum.Besides, as the dimensionality of the problem increases, the search of the phase space becomes more and more complicated, and the complexity of the algorithm increases.Another approach for continuous optimization is by using Grover's search algorithm [16].Continuous optimization problems can be discretized and mapped to a search problem, thereby solved by using Grover's algorithm.The Grover adaptive search algorithms iteratively apply Grover search to find the optimum value of an objective function [17][18][19][20][21][22][23], and can achieve quadratic speedup over classical search algorithms.However, these brute force methods are prohibitively expensive due to the large search space of the problems.
In a recent work [24], we proposed an efficient quantum algorithm for solving a search problem with nested structure through multistep quantum computation.The problem can be decomposed and the search space of the problem can be reduced in a polynomial rate.
The runtime of the algorithm is proportional to the number of steps of the algorithm.In this work, we generalize this algorithm for optimization problems with continuous variables.
The nested structured search problem [24] is a search problem that contains N items with one target item, and can be decomposed by using m [O(log N)] oracles to construct m Hamiltonians, respectively, as and where the set Π i contains N i marked items in the N items and |η i are the marked states associated with the marked items, and |η is the target state that defines the problem Hamiltonian of the search problem.These sets are nested as , respectively.The ratio N i−1 /N i are polynomial large, and The goal is to find the the target state |η that is associated with the target item in the set Π m .
Our algorithm solves the nested structured search problem by finding the ground state of the problem Hamiltonian H P via a multistep quantum computation process, which is realized through quantum resonant transition (QRT) [25,26].In this algorithm, a probe qubit is coupled to an n-qubit register R that represents the problem.We construct a sequence of intermediate Hamiltonians to form a Hamiltonian evolution path to the problem Hamiltonian as where Then we start from the ground state of the initial Hamiltonian, and evolve it through the ground states of the intermediate Hamiltonians sequentially through QRT to reach the ground state of the problem Hamiltonian.The ground state of an intermediate Hamiltonian is protected in an entangled state of the probe qubit and the register R, such that it can be used repeatedly without making copies.Therefore the algorithm circumvents the restriction of the no-cloning theorem [27,28] and realizes the multistep quantum computation.The algorithm can be run efficiently provided that: (i) the energy gap between the ground and the first excited states of each Hamiltonian H i and, (ii) the overlaps between ground states of any two adjacent Hamiltonians are not exponentially small.For the nested structured search problem, the conditions of the algorithm are satisfied since the ratio N i−1 /N i are polynomial large, therefore it can be solved efficiently, and the conditions for efficiently running our algorithm are not equivalent to those of the AQC algorithms [24].
In this algorithm, by using the Hamiltonians H P i sequentially in each step, the dimension of the search space of the problem is reduced in a polynomial rate, the solution state to the problem Hamiltonian is obtained step by step.The idea of reducing the search space in a polynomial rate step by step in our algorithm has a classical analogue as follows: suppose there are 80 balls, all of them have equal weights except one that is lighter than the others.
How to find the lighter ball?If we randomly pick up a ball and compare its weight with the other balls, this will take about 40 trials on average.If we have a balance, then how many times do we have to use the balance to find the lighter ball?According to information theory, the number of times the balance has to be used is log 80/ log 3 ≈ 4. The procedure is as follows: we divide all the 80 balls into 3 groups, each group has 27, 27 and 26 balls, respectively; then pick up the two groups that both have 27 balls, and use the balance to determine if they have equal weights.If the answer is positive, pick the group with 26 balls and divide it into 3 groups again: 9, 9, 8; otherwise, take the group that is lighter and divide it into three new groups: 9, 9, 9.This process can be repeated until the lighter ball is found.
In this example, we can see that the problem is divided into a series of nested sub-problems and the size of the search space is reduced in a rate about 1/3 by using a balance.The target ball is found through an iterative procedure and the cost is reduced exponentially.
By using a different oracle in each step, the QRT procedure in our algorithm emulates the usage of the balance in solving the nested structured search problem.
The procedure for solving the nested structured search problem can be applied for optimization problems that are transformed to finding the ground state of a problem Hamiltonian in quantum computation.Here, we propose a quantum algorithm based on multistep quantum computation for optimization problems with continuous variables.We first discretize the variables of the objective function to construct the state space of the problem.Then we construct a sequence of intermediate Hamiltonians to reach the problem Hamiltonian by decomposing the problem using a set of threshold values, and apply a multistep quantum computation process to reduce the search space of the problem step by step.The solution vector to the optimization problem is narrowed in a small state space and can be determined efficiently through measurements.If the search spaces of the Hamiltonians are reduced in a polynomial rate by using an appropriate set of threshold values, then the optimum of the function can be obtained efficiently.Meanwhile if the global minimum of the optimization problem is in the state space of the problem, then it can be obtained efficiently.The problem in many optimization algorithms where the trial vector is trapped in a deep local minimum and missing the global minimum can be avoided in our algorithm, provided the above conditions are satisfied.In quantum computing, the dimension of the Hilbert space of the qubits increases exponentially with the number of qubits, it is more efficient to represent a large state space on a quantum computer than on a classical computer, therefore increasing the probability of finding the global minimum of the problem.This paper is organized as follows: in Sec.II, we describe the quantum algorithm for optimization problems with continuous variables based on multistep quantum computation; in Sec.III, we apply the algorithm for some test optimization problems, and we close with a discussion.

II. QUANTUM OPTIMIZATION ALGORITHM BASED ON MULTISTEP QUAN-TUM COMPUTATION
Let S be the domain of x, an optimization problem can be formulated as a minimization problem: where F is a real-valued objective function and x is the vector of the variables.Here we focus on optimization problems with continuous variables, which can be described as follows: for a real-valued function of r variables, F (x 1 , x 2 , • • • , x r ), find a vector of the variables such that the function has the minimum value.In the following, we present a quantum optimization algorithm based on multistep quantum computation for this problem.
We discretize the continuous variables in the function domain into intervals of same length for all the variables, and map the problem on a quantum computer.For simplicity, suppose each variable is discretized into l elements in its definition domain, the dimension of the state space of the function is l r .We prepare r quantum registers and each register contains ⌈log 2 l⌉ qubits that represents the elements of the variable.Therefore n = r⌈log 2 l⌉ qubits form the register R that represents the problem with state space of size N = 2 n on a quantum computer.A vector of the discretized variables x where Then we construct m Hamiltonians as: where and H Pm = H m = H P .This can be achieved by using an oracle that recognizes whether where , and M i is an approximate estimation of N i .We have demonstrated that as M i = N i , the conditions for efficiently running the algorithm are satisfied provided that the ratio N i−1 /N i are polynomial large [24].The parameters M i can be estimated efficiently by using the Monte Carlo sampling method [31], and we can adjust the threshold values such that the ratio N i−1 /N i are polynomial large.Detailed analysis of the effect of M i on the efficiency of the algorithm is presented in the appendix.The ground state of H P can be obtained through the following multistep quantum computation process based on QRT in m steps.We use the ith step of the algorithm to illustrate the procedures.
In the ith step, given the Hamiltonian H i−1 , its ground state eigenvalue E of H i by using the QRT method.The algorithm requires (n + 1) qubits with a probe qubit coupling to the n-qubit register R. The algorithm Hamiltonian of the ith step is constructed as where I N is the N-dimensional identity operator, and σ x and σ z are the Pauli matrices.The first term in Eq. ( 9) is the Hamiltonian of the probe qubit, the second term contains the Hamiltonian of the register R and describes the interaction between the probe qubit and R, and the third term is a perturbation with c ≪ 1.The parameter α i is used to re-scale the energy levels of H i−1 , and the ground state energy of α i H i−1 is used as a reference energy level to the ground state eigenvalue 0 of H i .The initial state of the (n + 1) qubits is set as |1 |ϕ . First we obtain the eigenvalue E (i) 0 of H i by using the QRT method through varying the frequency of the probe qubit as shown in Ref. [24].Then we set = ω for resonant transition between the probe qubit and the transition between states |ϕ (i−1) 0 and |ϕ (i) 0 is satisfied.When obtaining the eigenvalue E (i) 0 of H i , we can also obtain the overlap g between the ground states of H i−1 and H i through the Rabi's formula [32].Then we can set the optimal runtime t i = π/(2cg (i) 0 ) at which the probability for the system to be evolved to the state |0 |ϕ (i) 0 reaches its maximum.
The procedures for obtaining the ground state of H i are summarized as follows: (i) Initialize the probe qubit to its excited state |1 and the register R in state |ϕ (i−1) 0

;
(ii) Implement the unitary evolution operator U(t i ) = exp −iH (i) t i ; (iii) Read out the state of the probe qubit.
The system is approximately in state 1 − p as the resonant transition occurs, where p is the decay probability of the probe qubit of the ith step.The state |ϕ (i−1) 0 from the previous step is protected in this entangled state.By performing a measurement on the probe qubit, if the probe decays to its ground state |0 , it indicates that the resonant transition occurs and the system evolves from the state |1 |ϕ  , then we repeat procedures ii)-iii) until the probe qubit decays to its ground state |0 .Therefore we can obtain the ground state |ϕ By protecting the state |ϕ (i−1) 0 through entanglement, the state can be used repeatedly without copying it, such that the algorithm realizes multistep quantum computation.The runtime of the algorithm is proportional to the number of steps of the algorithm, and the success probability of the algorithm is polynomial large by setting the coupling coefficient c appropriately [24].After running the algorithm for m steps, we obtain the ground state of the problem Hamiltonian H P , which is a superposition state of a few CBS with eigenvalues below the threshold value d m .Then we can perform measurement on the state and find the CBS that has the minimum function value, therefore solving the optimization problem.We can run the algorithm for a few rounds by discretizing the variables in the neighborhood of the optimized vector to improve the precision of the solution to the optimization problem.
The algorithm can be run efficiently if both the energy gap between the ground and the first excited states of each Hamiltonian H i and the overlap between ground states of any two adjacent Hamiltonians g (i) 0 are not exponentially small.By solving the eigen-problem of the Hamiltonian H i , these conditions can be satisfied if the ratio N i−1 /N i are polynomial large, and the parameters M i are set such that: the point (N i /N, M i /N) is far away from the neighborhood of the point (0, 1/2), and 2M i (N − N i ) < N 2 .Detailed analysis is shown in the appendix.

III. APPLICATION OF THE ALGORITHM FOR SOME TEST FUNCTIONS
We now apply the quantum optimization algorithm described above for some test functions of optimization problem: the Damavandi function, the Griewank function and the Price function.

A. The Damavandi function
The two dimensional Damavandi function is defined as and the graph of this function is shown in  We can see that the dimension of the state space is reduced to a few CBS after a number of steps.Therefore the final state that encodes the solution to the optimization problem can be readout and checked efficiently to find the global optimum of the function.

B. The Griewank function
The Griewank function has the form Fig. 3 shows the second-order Griewank function with two variables, we can see that the function has many local minima.For classical optimization algorithms, it is very easy for a trial vector to be trapped in one of the local minima, while missing the global minimum of the function.This situation can be avoided in our algorithm.
We  the problem is reduced smoothly in each step of the algorithm in rate of {0.31, 0.39, 0.45, 0.55, 0.44, 0.25, 0.36, 0.44, 0.24, 0.30, 0.22, 0.20}, respectively.After running the algorithm for a number of steps, the state space of the problem is reduced to a very small space and can be readout to calculate the corresponding function value and find the global minimum.

C. The Price function
The Price01 function can be written in form of where A m are determined by Eqs. ( 6) and (7), and the sizes of the In the following we solve the eigen-problem of the Hamiltonian H i .Let where Then in basis {|q i } q i ∈A i , |q ⊥ i , we have and Then Define ẽ to be an N i × 1 vector of all ones and we can rewrite H i as (i) Define the vector space V = {e, e n } of dimension 2. Then from Eq. (A12), ∀x ∈ V ⊥ , we have H i x = − 1 − M i N x, therefore the eigenvalues are − 1 − M i N , corresponding to N i − 1 eigenvectors.
(ii) The vector space of V can be spanned by vectors It is easy to check that and Then we can verify that Solving the eigen-problem of the above matrix, we can obtain the eigenvalues: as: Correspondingly, the state |V (i−1) − can be written as: Thus the overlap between the ground states |V ) r is represented by state |ij • • • k , where x (j) s represents the jth element of the variable x s .The states |i , |j , • • • , |k are binary representation of the elements x (i) 1 , x (j) 2 , . .., x (k) r on the quantum registers.These vectors form the computational basis states (CBS) of r quantum registers of dimension N as |J = |i |j . . .|k , J = 0, 1, . .., N −1, and the corresponding function value is F x qr) r is the minimum of the function F .By using an oracle O F where O F |J |0 = |J |F (J) , the Hamiltonian of the optimization problem can be constructed as previous step, we are to prepare the ground state |ϕ (i) 0 if the probe qubit stays in state |1 , it means that the register R remains in state |ϕ (i−1) 0

Fig. 1 .
It has a very sharp global minimum of zero at {x 1 = 2, x 2 = 2}.For classical optimization algorithms based on the gradients methods, it is very easy for a trial vector to be trapped in the bowl-like local minimum, while missing the global minimum.The overall success probability of current global optimization algorithms for finding the global minimum of this function is about 0.25% [33].To apply our algorithm for this optimization problem, the two variables of the Damavandi function are discretized into 281 elements evenly with an interval of 0.05 in the range [0, 14].The dimension of the state space of the function is 78961.By discretizing the value of the function in an interval of 0.01 in the range of [0, 149], and counting the number of states in each interval, we can obtain the distribution of states of the function in each energy interval as shown in Fig. 2. The largest degeneracy is about 56, which is a small number compare to the dimension of the state space.We construct a set of threshold values as {d 1 = 70, d 2 = 30, d 3 = 15, d 4 = 7, d 5 = 4, d 6 = 3, d 7 = 2.5, d 8 = 2.2, d 9 = 2.1, d 10 = 2.02, d 11 = 2.0, discretize the two variables of the Griewank function into 801 elements evenly with interval of 0.1 in the range [−40, 40].The dimension of the state space of the function is 641601.By discretizing the function value in intervals of 0.0001, the distribution of states in each energy interval of the function is shown in Fig. 4. The largest degeneracy is 32.

FIG. 2 .
FIG. 2. Distribution of states by discretizing the value of the Damavandi function in intervals of 0.01.

)FIG. 6 .
FIG. 6. (Color online) Energy gap between the ground and the first excited states of the Hamiltonian H i as a function of a and b.

1 + 1 x
FIG. 8. Ratio x (i) 1 /x (i) 2 between components x (i) 1 and x (i)2 of the ground state of the Hamiltonian larger or less than a threshold value d i .It is a comparison logic circuitry and can be implemented efficiently on a quantum computer[29, p.264] [17, 20, 30].The CBS associated with integers that are less than or equal to d i form a set A i with size N i .They have thenested structure as A m ⊂ A m−1 ⊂ • • • ⊂ A 1 .The ground state of the problem Hamiltonian H P contains CBS in A m with eigenvalues that are below the threshold value d m .We construct a sequence of Hamiltonians that form a Hamiltonian evolution path to the problem Hamiltonian H P as and run the algorithm.The dimension of the corresponding state space in each step of the algorithm are reduced to {56634, 24939, 11573, 4452, 1772, 892, 448, 178, 94, 20, 5, 1}, respectively.The dimension of the state space is reduced smoothly with reduction rates of {0.717, 0.440, 0.464, 0.385, 0.398, 0.503, 0.502, 0.397, 0.528, 0.213, 0.250, 0.200} in each step of the algorithm, respectively.The parameter M i can be estimated through Monte Carlo sampling, if it is set approximately as N i above, the conditions of the algorithm can be satisfied.The ratio M i /N that are closest to 1/2 is 0.316.
respectively, and N 1 > • • • > N m .Let H m = H P = H Pm and the set A m contains the target states |q with size N m .We construct a Hamiltonian evolution path H 0 → H 1 → • • • → H m = H P and start from the ground state |ϕ