Sweeping quantum annealing algorithm for constrained optimization problems


 As a typical quantum computing algorithm, quantum annealing is widely used in the optimization of glass-like problems to find the best solution. However, the optimization problems in constrained complex systems usually involve topological structures, and the performance of the quantum annealing algorithm is still largely unknown. Here, we take an Ising system as a typical example with local constraints accompanied by intrinsic topological properties that can be implemented on quantum computing platforms such as the D-wave machine, and study the effectiveness of the quantum annealing algorithm in its optimization and compare it with that of the thermal annealing. We find that although conventional quantum annealing is difficult for the optimization of constrained topological problems, a generalized algorithm --- the sweeping quantum annealing method --- can be designed and solve the problem with better efficiency than both conventional quantum and thermal annealing algorithms. The sweeping quantum annealing algorithm, therefore, opens up a promising avenue for quantum computing of constrained problems and can be readily employed on the optimizations in quantum material design, engineering, and even social sciences.

S Olutions to complex optimization problems often have direct connections with the finding of ground state of many-body systems. For example, many optimization problems can be translated into Ising Hamiltonian with random couplings, and the minimum energy configuration of the Hamiltonian corresponds to their optimal solutions [1][2][3][4][5]. In classical many-body systems, the real ground state is found by traversing different states through thermal fluctuations of thermal annealing (TA). For quantum systems, superposition and quantum tunneling can search for multiple classical configurations at the same time, and find the final ground state through quantum annealing (QA) [6].
In experiments, some confirmations of the superiority of quantum annealing have been obtained, for example, results on the spin-1/2 disordered Ising ferromagnet LiHo 0.44 Y 0.56 F 4 have suggested that quantum annealing procedure, which attempts to follow the adiabatic quantum algorithm [7][8][9], works considerably well than classical thermal annealing [10,11].
Following the recent technological advancements in manufacturing coupled qubit systems, the idea of building dedicated quantum annealing equipment to solve optimization problems has attracted much attention, and such prototype device (quantum computer) has been implemented on superconducting flux qubits [12][13][14]. Currently, quantum annealing computers have also been commercialized, such as the Dwave machine. Although it is not a general-purpose quantum computer, it does show higher efficiency than classical computers in certain optimization problems [15][16][17][18]. However, the optimization problems in the real world are usually more complicated than those mentioned above. For example, besides the nearly extensive degeneracy that is shared by all the glassy Ising type Hamiltonians, the real optimization problems are often mapped to Hamiltonian with local constrained conditions, which greatly increase the difficulty of optimization. From a physical perspective, the difficulty of such problems lies in the fact that the lowenergy Hilbert space not only has many nearly degenerate local minima but is also restricted, and one needs to search for the optimal solution therein.
To understand the constrained optimization problem better, one can borrow daily experience: assuming a kindergarten in the present COVID-19 pandemic, we hope that every child can sit in a row separated by a seat to keep the social distance, which is a local constraint. In the state shown in Fig. 1(a), most of the children have met this condition except one on the left edge of the row, creating a "topological defect" that needs to be optimized. This can be achieved by a global operation of every child moving one position at the same time, rather than repeated simple local operations. This oversimplified example shows that the local constraints usually introduce a topological perspective into the problem and would require global movement (children change seats together) that could change the topological sectors to properly find the minimal. These constrained optimization problems have important applications in many fields, such as structural optimization and material design in engineering [23,24], ethnic differentiation, and interpersonal relationships in social sciences [25][26][27], but the analysis of the annealing algorithms (thermal and quantum) of constrained systems is still rare.
In this article, we construct an Ising model with local constraints and focus on the comparison of the effectiveness of different annealing schemes to find the ground state of the problem. We found that although the normal quantum annealing is difficult (the thermal annealing is slightly better but also difficult) for the optimization of such model, a generalized algorithm -the sweeping quantum annealing method -can be designed and solve the problem with better efficiency than both simple quantum and thermal annealing algorithms. Our sweeping quantum annealing algorithm can also be easily implemented on quantum computing platforms such as D-wave machine and therefore opens up a promising avenue for the application of quantum computing on realistic constrained problems.

MICROSCOPIC SPIN MODEL
Similar to the usual annealing problems, we can translate the constrained optimization problem into a problem of finding the ground state of a Hamiltonian. To facilitate the implementation of our algorithm on quantum computers such as D-wave machine, we study a sufficiently simple constrained system -the antiferromagnetic (AF) Ising model on a triangular lattice. As shown in Fig. 1(b), the phase diagram of the model with transverse field (to introduce quantum fluctuations) has been mapped out from previous theoretical and numerical works [19][20][21][22]. The constrained problem we would like to solve is seeking the ground state at the limits where h and T are approaching zero. Due to the antiferromagnetic interaction, every triangle must be composed of two parallel spins and the other one opposite to them, and it can be taken as the constraint condition dubbed "triangle rule".
It is straightforward to see that, exactly at the classical Ising limits of h = T = 0, the ground state has extensive degeneracy: any configuration that satisfies the triangle rule is a legitimate ground state. Then, we introduce a small anisotropy for the nearest-neighbor (NN) AF Ising model to partially lift the degeneracy, where σ is the Pauli matrix, i j x and i j ∧ represent the nearest-neighbor sites on the horizontal bonds and the interchain bonds, respectively. If we set J = 1 and J x = 0.9, the ground state of H A is a stripe-ordered phase as shown in the left panel of Fig. 1(c) with the ground state energy exactly computable (for a 6 × 6 lattice, the energy for stripe state is E g = −39.6). However, it is not a simple task to anneal the system to the stripe phase when taking a random configuration as the initial state at high temperature or in a large transverse field. This is because flipping a spin arbitrarily will result in violation of the triangle rule with energy penalty in the order of J, i.e., it is difficult to evolve from arbitrary state to ground state through local changes. Thus, the annealing process towards the ground state of H A provides an ideal object to test the performance of different annealing methods. By comparing thermal annealing and various quantum annealing methods, we finally invent the sweeping quantum annealing algorithm which exhibits the best performance.

THERMAL ANNEALING VS QUANTUM ANNEALING
The thermal annealing (TA) scheme is to slowly reduce the temperature of H A from a region where local minima can be easily escaped from through tunneling to nearly zero. We simulate the thermal annealing process by using the Metropolis Monte Carlo [28] with local spin flips. The Boltzmann constant is set to k B = 1 for convenience, and we perform thermal annealing through linearly decreasing temperature from T = 5.00 to T = 0.05 with cooling interval ∆T = 0.05 ( Fig. 1(b), red arrow). In contrast, quantum annealing (QA) requires introducing the quantum fluctuation, so that the Hamiltonian becomes slowly reduced to zero such that the original Ising Hamiltonian H A is recovered at the end of the process. The quantum annealing is implemented by utilizing the quantum Monte Carlo (QMC) simulation [29][30][31], and we adopt a continuous-time QMC -stochastic series expansion (SSE) method [32][33][34] to avoid the Trotter discretization error [29,35]. The QMC schemes have been successfully employed in many previous works (including ours) on quantum Ising models [21,22,36].
Similar to thermal annealing, we carry out quantum annealing with linearly decreasing transverse field from h = 5.00 to 0.00 with annealing interval ∆h = 0.05 ( Fig. 1(b), blue arrow), and the temperature T = 0.05 is low enough to ignore the thermal effect. In order to compare thermal and quantum annealing, we study how the energy E evolves under the same annealing time, in the unit of Monte Carlo step (MCS). We have simulated the system with sizes 6 × 6 with periodic boundary conditions and parallelized the simulations on 64 independent Markov chains to achieve better statistics.
As shown in Fig. 2, the QA turns out to have poor performance with energy difference ∆E = E − E g ≈ 5. This can be understood from the phase diagram. The annealing path of TA is mostly inside the disordered phase whereas the path of QA experiences the disordered, clock-ordered, and critical phases.
The phase transitions between these phases obviously disturb the effectiveness of the simulated annealing, in particular, the quantum phase transition where the three phases meet at finite h c is known to have emergent continuous U(1) symmetry [20,21] and the significant critical slowing down of the quantum annealing is thus expected [37,38].
To overcome such difficulty in QA, we have also tried the quantum annealing with inhomogeneous field It has been suggested that disordered magnetic field would weaken the effect of phase transition and improve the efficiency of quantum annealing [39][40][41][42][43]. At initial time t = 0, we randomly generate values from 0 to 10 to for each h i . Then all h i are linearly reduced to 0 in the QA process. We denote such annealing as "QA-h" in Fig. 2. However, it turns out that our constrained optimization problem is more complicated than those overcome by random quantum fields [39][40][41][42][43], the result of QA-h looks similar to QA without obvious improvement.

SWEEPING QUANTUM ANNEALING ALGORITHM
As discussed in the introduction, the constrained optimization problems usually involve non-local excitations and require global operations, e.g. Fig. 1(a). However, the system can easily reach the global  Fig. 1(b). The energy is computed for 6 × 6 system and the true ground state energy E g = −39.6 is denoted by the grey dashed line. The expectation value of E and its error bars are obtained from averaging over 64 independent annealing runs.
The straightforward QA and QA-h have the worst annealing performance, whereas TA, SQA-1, SQA-2, SQA-3 have better performance and the best among these is the SQA-6. minimum, if we can change the boundary condition, such as gluing the first and last seats together to construct the periodical boundary condition and then cutting the connection between first and second seats. It inspires us to invent the sweeping quantum annealing algorithm (SQA). In practice, as shown in the right panel of Fig. 1(c), we can cut the system to artificially create an open boundary. The corresponding energy shift of one edge spin (e.g., the collapsed spin-up and spin-down in Fig. 1(c), right panel) flipping process is 2 × (2J e − J x ) and positive at J e = J = 0.9J x , where J e is the boundary interaction. It implies that the boundaries are so "hard" for the topological defect passing through [44][45][46]. Therefore, in order to soften the boundary, J e should be set to J x /2 so that the edge spins can be flipped freely. Then the edges need to be glued to restore the original boundary conditions when utilizing the quantum annealing.
After that, we move the process of cutting and gluing to another position. If n positions are selected in the whole system, they will be located with equal distance for keeping translational symmetry, and we label such annealing scheme with SQA-n. Fig. 2 shows the results obtained from different SQA-n, and it clearly demonstrates SQA-6 is of the best performance, as it sweeps over all the edges of 6 × 6 lattice. Meanwhile, SQA-n outperforms TA and QA at large annealing time and its performance is positively associated with the number of cutting positions. The implementation of SQA-n on H A can be presented with the pseudo-code shown in Table I.
One annealing step with QMC

10.
Set J x− = J x− + ∆J x ; J e = J e + ∆J e ! smoothly restoring the boundary condition 11. End

12.
Move to next cutting position

SQA AND TOPOLOGY
The high performance of SQA can be understood in the topological aspect. As shown in Fig. 3(b), the constraint (triangle rule) will restrict the Hilbert space so that the TFIM is equivalent to a dual honeycomb lattice dimer model [47][48][49][50][51][52] by mapping two parallel spins to one dimer on the dual lattice. Then, the triangle rule is translated to dimer counterpart that each dual site belongs to one and only one dimer. The honeycomb lattice can be divided into A and B sublattices, then the dimer density n d can be understood as the lattice electric field defined as E = n d − 1/3 on the dual bond from A to B dual-site, the constraint can be written as Gauss law ∇ · E = 0. Therefore, these many-body configurations with constraints can be mapped to lattice electromagnetic fields which are naturally classified to different topological sectors [19,53].
The ground state stripe phase in Fig. 3(c) is composed of horizontal ferromagnetic bonds and the interchain antiferromagnetic bonds. The number of horizontal antiferromagnetic bonds is conserved in each row with keeping constraint, and linking such AF bonds on different chains will construct the one dimensional global (topological) defect -the domain wall (DW) (e.g., the ones shown in Fig. 3(d)). The topological sectors can be distinguished by N D the number of the DWs, and changing it requires generating a pair of spinons (constraint-breaking triangles) connected by DWs as shown in Fig.3(f), letting them go around the cylinder of connected boundary to meet and annihilate, and then the DW number changes by two. However, exciting spinons at low temperatures is extremely hard, so that it is not easy to go across topological sectors.
Both quantum and thermal fluctuations can only deform the DWs but cannot create/annihilate them. In essence, the optimization problem of constraints is an annealing problem of different topological sectors as Fig. 3(a). In this case, the number of states in a topological sector is positively correlated with the number of DWs, so the stripe configuration of the ground state is hard to be reached in Hilbert space as one needs an annealing algorithm that can both tunnel through different topological sectors and reduce the energy.
Therefore, neither TA nor QA can anneal the system to the best solution, nor is QA-h.
In contrast, SQA considers the topological aspect by introducing sweeping. As demonstrated in inefficiencies. Although the TA can also change the topological sectors due to the creation of pairs of spinon excitations in Fig. 3(f), clearly the best of all annealing schemes is the SQA-6 which can straightforwardly change the topology of the system at any position.

CONCLUSION
As mentioned above, in realistic optimization problems, local constraints exist ubiquitously. For such constrained optimization problems, there is still a lack of systematic quantum computing research. In this work, we find that local constraints often bring topological properties, which greatly reduces the efficiency of both conventional thermal and quantum annealing. Borrowing the idea of changing topology by cutting and gluing the system, we invent a generalized algorithm -the sweeping quantum annealing method. After comparing with conventional quantum and thermal annealing algorithm, we find that this algorithm presents high efficiency and validity. Our sweeping quantum annealing algorithm can also be easily implemented on quantum computing platforms such as D-wave machine and therefore opens up an effective and innovative way for the application of quantum computing on realistic constrained problems. In the future, we will check the effectiveness of SQA under noises and dissipations [54,55] which commonly exist in the quantum annealer. The new algorithm can also be implemented for solving more realistic problems, such as opinion formation [56] or elections [57].

METHODS
Classical Monte Carlo for thermal annealing. We use Metropolis Monte Carlo [28] which flip single spin to simulate Eq. (1) thermal annealing. Its general process is as follows: assume we have reached configuration C 1 with energy E 1 . Firstly we choose a spin randomly and flip it to get a new configuration C 2 with energy E 2 . We use the acceptance probability P = min(e −β (E 2 −E 1 ) , 1) to determine whether to update the system to configuration C 2 , otherwise we keep configuration C 1 . For the thermal annealing process, we do the Monte Carlo simulation from T = 5 to T = 0.05, the decreasing interval is ∆T = 0.05, and MCS of every interval is fixed.
Quantum Monte Carlo for quantum annealing. For the quantum annealing in this paper, we use a quantum Monte Carlo (QMC) method with stochastic series expansion (SSE) algorithm [32][33][34]. In this method, the evaluation of partition function Z is done by a Taylor expansion, and the trace is taken by summing over a complete set of suitably chosen basis.
By writing the Hamiltonian as the sum of a set of operators whose matrix elements are easy to calculate H = − ∑ i H i and truncating the Taylor expansion at a sufficiently large cutoff M, we can further obtain To carry out the summation, a Markov chain Monte Carlo procedure can be used to sample the operator sequence {i p } and the trial state α. One step of the update process contains diagonal update and offdiagonal update. In the diagonal update, the diagonal operators are inserted into and removed from the operator sequence. Meanwhile, in the off-diagonal update, the diagonal and off-diagonal operators can be converted into each other. In the quantum annealing process, we decrease the transverse field from h = 5.00 to h = 0.00, with a decreasing interval ∆h = 0.05. A tiny residual field is kept to facilitate the update process.