Quantum annealing for systems of polynomial equations

Numerous scientific and engineering applications require numerically solving systems of equations. Classically solving a general set of polynomial equations requires iterative solvers, while linear equations may be solved either by direct matrix inversion or iteratively with judicious preconditioning. However, the convergence of iterative algorithms is highly variable and depends, in part, on the condition number. We present a direct method for solving general systems of polynomial equations based on quantum annealing, and we validate this method using a system of second-order polynomial equations solved on a commercially available quantum annealer. We then demonstrate applications for linear regression, and discuss in more detail the scaling behavior for general systems of linear equations with respect to problem size, condition number, and search precision. Finally, we define an iterative annealing process and demonstrate its efficacy in solving a linear system to a tolerance of 10−8.


INTRODUCTION.
Many problems in science, engineering, and mathematics can be reduced to solving systems of equations with notable examples in modeling and simulation of physical systems and the verification and validation of engineering designs.Conventional methods for solving linear systems range from exact methods, such as matrix diagonalization, to iterative methods, such as fixed-point solvers.The advent of quantum computing has opened up the possibility of new methods for solving these challenging problems.For example, a quantum algorithm for solving systems of linear equations was established for gate-based quantum computers [1] and demonstrated with small-scale problem instances [2].Additionally, an algorithm for solving linear systems within the adiabatic quantum computing model [3] was experimentally demonstrated [4], followed by a more recent proposal [5].
In this Letter, we present an approach for solving a more general system of polynomial equations using quantum annealing.Akin to adiabatic quantum optimization [6], quantum annealing prepares a quantum statistical distribution that approximates a sought-after solution by applying a slowly changing, time-dependent Hamiltonian [7], where measurements drawn from the distribution represent candidate solutions.Unlike adiabatic quantum computing, quantum annealing permits non-adiabatic dynamics at non-zero temperature, making this approach easier to realize experimentally but also more challenging to distinguish quantum mechanically [7][8][9][10][11][12][13][14].While examples of non-trivial advantages have been observed for fixed-size problem instances [15][16][17][18][19], more general statements about computational complexity remain unresolved [20].
Based on the principles of quantum annealing, we demonstrate a novel probabilistic solver for ill-conditioned systems of linear equations.We use the discretized Dirac equation Dφ = χ from lattice quantum chromodynamics (QCD) as a motivating example, in which the Dirac operator D encodes the fundamental interactions between quarks and gluons that lead to formation of hadrons.The solution to the discretized equation is currently the only approach to evaluating QCD non-perturbatively, and yields the quark propagator φ for a given source χ.However, well-known numerical challenges slow convergence with conventional solvers [21,22].We use quantum annealing to solve this system and characterize the performance from experimental demonstrations with a commercial quantum annealer.Our implementation recovers the correct solution with high precision and accuracy for problem sizes that fit within the available hardware.

Polynomial Systems of Equations.
We consider the system of N polynomial equations where i ∈ {1, . . ., N }, and P (n) is a rank n + 1 tensor of known real-valued coefficients for the polynomial of order n, and the real-valued vector x denotes the solution.Truncating to first order recovers a linear system of equations, i.e., P i .Existing approaches for solving Eq. ( 1) include direct diagonalization using Gauss-Jordan elimination or iterative methods such as conjugate-gradient.In practice, direct diagonalization is limited in computational efficiency, as those methods scale sharply with the size of the matrix.By contrast, iterative methods may have greater computational efficiency but the performance and stability are often sensitive to the input matrix.For example, the Dirac operator D becomes more singular as the quark mass is lowered to its physical value.The condition number of the resulting D matrix slows the convergence of the linear solver and causes low quark mass calculations to be prohibitively expensive.
Preconditioning improves convergence by transforming the input as M −1 Ax = M −1 b, where the preconditioner M must be inexpensive to invert and M −1 should be "close" to A −1 , so that M −1 A resembles a matrix close to unity.Identifying an effective preconditioner plays an important role in numerical convergence of iterative methods [23][24][25][26][27].For lattice QCD applications, the low-lying spectrum of the Dirac operator slows iterative convergence and preconditioning has been used to project out these low-lying modes.Acquiring the low-lying eigenpairs or singular triplets of D is in general computationally expensive because it requires the use of additional iterative methods that also suffer from critical slowing down and poor scaling with the physical volume.Solutions to address this issue include EigCG [21,22], inexact deflation [23], and adaptive multigrid [24][25][26][27].

Quantum Annealing for Polynomial Solvers
Quantum annealing offers an alternative approach to solving ill-condition systems of equations.By recasting Eq. ( 1) as the scalar minimization problem arg min we reduce finding the solution to minimizing an x-variable objective function.We establish the connection to quantum annealing by first approximating the objective function using an R-bit representation for each real-valued variable x j as xj = a j R−1 r=0 2 r ψ r,j + b j . ( where ψ r,j ∈ {0, 1} and a i and b i control the range and numerical shift, respectively, such that xj ∈ [b j , b j + 2 R−1 a j ] in increments of a j .Direct substitution yields the multibody objective function where its ground state solves a system of polynomial equations.For current commercial quantum annealers, ancillary binary variables are required to decompose multi-linear terms into two-body interactions.This quadratization process introduces more variables into the problem with additional penalty terms [28][29][30]. Quantum Annealing for Linear Solvers Restricting the transformation of Eq. (3) to a system of linear equations simplifies Eq. ( 4) by involving only terms up to two-body interactions.As a result the explicit binary quadratic optimization problem (QUBO) is . . .a 1 a N P (1) 1N . . . . . . . . . where n,i .Note that the matrix Q is real-valued and symmetric.In addition, constant terms that arise from the substitution of Eq. (3) are omitted, as these simplifications shift the energy spectrum by a constant but leave the minimizing solution unchanged.
Linear Least Squares -Least squares minimization is a standard approach for regression analysis of an overdetermined system [31].In general, the unknown parameters of interest, dependent variables and independent variables can be continuous and correlated.Generalized least squares offers a strategy for defining the best-fit such that the parameters in the regression function best describe the given data.
Given a set of N identical and independent observations of the independent {x i : i ∈ {1, ..., X}} dependent {y i;n : i ∈ {1, ..., X}, g ∈ {1, ..., N }} variable, the mean and covariance of y i follows where the angle brackets denote the expectation value over N observations.A fitting function F (x i , p) may be defined with respect to the set of P unknown parameters p = {p n : n ∈ {1, ..., P }}, and a corresponding objective function for generalized least squares may be defined as where the optimal value for the set p is determined by minimizing the value of χ 2 .Restriction to linear least squares demands that the fitting function is linear in the unknown parameters, and therefore may be written in the form where f n (x i ) can be any function.The solution for linear regression is obtained by expanding Eq. ( 6) with Eq. ( 7) and yields The extrema of the objective function can be determined by taking the derivative of Eq. ( 8) with respect to p n yielding a matrix equation of the form A (1) analogous to Eq. ( 1) where The solution to least squares minimization can then be mapped to a QUBO problem following Eq.( 5), and amenable to methods of quantum annealing.

DISCUSSION Linear Least Squares Example
As an example, consider the following artificially generated data where x ∈ Z : x ∈ [0, 49].The Toeplitz correlation matrix [32] is chosen to simulate a correlated time-series dataset, where the correlations decay exponentially as a function of x.Following the notation in Eq. ( 7), we assume a linear fit and we estimate the parameters A n given the data D(x).Using Eq. ( 3), we express each parameter A i as a 4-bit unsigned integer and construct the problem Hamiltonian following Eqs.( 5) and ( 9).The required 12 logical spins (3 parameters × 4-bit representation) support a total of 4096 possible solutions.Explicit evaluation finds the true ground-state to have energy E 0 = −1.418and eigenstate which corresponds to the parameter values These correct coefficients for the generating function in Eq. ( 10) verify the design of the algorithm.
We next test the algorithm by solving the objective function using quantum annealing.The target Hamiltonian of Eq. ( 10) is mapped into the D-Wave 2000Q, as discussed above.Results for 100,000 independent evaluations are acquired using an annealing schedule with T = 200 µs.The histogram (blue) and cumulative distribution (red) of solutions are presented in Fig. 1, with 0.5% of the results reproducing the ground state given by Eq. ( 15).In addition, the lowest 0.8% of the true eigenvalue spectrum is obtained by 10% of these solutions with the overall results biased towards the lower-lying energy eigenspectrum.

Conditioned Systems of Linear Equations
In this section we show results and scaling of a classical method and the quantum annealer.One of the criteria for categorizing the "difficulty" of a linear system is condition number.The condition number of a matrix is defined as the ratio of maximum and minimum singular values.
In the case of symmetric matrices, this is equivalent to the ratio of largest and smallest eigenvalues.All of our numerical experiments in the subsequent sub-sections are performed with symmetric matrices.We vary our test matrices in two ways: 1) vary the rank of the matrix while holding the condition number fixed, 2) the dimension is held constant with varying condition number.A common metric to assess how accurately the linear system has been solved is the L 2 -norm of the relative residual.
For conjugate gradient, a tolerance for Eq. ( 20) is utilized as a terminating criterion and the number of iterations when this point is reached is recorded.For quantum annealing the role of the relative residual is more subtle.The annealer is run many times and the lowest energy eigenpair is returned.The eigenvector from this set is substituted for x approx , allowing a relative residual to be defined for the total anneal.Classical Solutions -For the examples with a classical linear solver, conjugate gradient is used on a matrix of rank 12, with varying condition number.Although conjugate gradient is not the optimal choice for classically solving such systems, the scaling comparison in condition number with the quantum algorithm is informative.Fig. 2 shows slightly worse than square root scaling of conjugate gradient with condition number.
The matrices from this result are constructed by creating a random unitary matrix of rank 12, denoted as U .A linear spacing of real eigenvalues is chosen with the maximum equaling 5 and the minimum selected precisely to give the condition number that is desired.The matrices are trivially formed as A = U ΛU † .A common right-hand side is taken for all A: a vector of length 12 with linearly spaced decimals between 1 and -1.
Quantum annealing -In the following section, we demonstrate the scaling of the annealing algorithm under varying problem size, condition number, and precision of the search space.We conclude by applying the algorithm iteratively on a fixed problem and study the convergence of the relative residual.
Problem size -We observe that as the problem size increases linearly, the relative residual remains relatively constant.This is consistent with the notion that at fixed condition number, the range of numerical values in the result remain unchanged, and the algorithm provides constant precision under such conditions.However, for problem sizes beyond 16 linear equations, hardware limitations devoid us of the possibility of sampling the correct ground state, as demonstrated by the divergence of the annealed result versus the theoretical solution.In the on-set of this limitation, the residual grows more rapidly as expected.
In Fig. 3a, we study the scaling behavior for a fixed condition number of 1.1 using a 2-bit representation for each of n parameters.Due to prior knowledge of the conjugategradient solution, the search space for all n parameters are fixed for the set of problems, and are known to encompass the minimum and maximum results of the solution vector x.Additionally, knowledge of the result allows us to identify the ground-state QUBO solution by minimizing the difference between the conjugate-gradient and QUBO results (the forward error), and studies the theoretical scaling of the algorithm absent of current hardware limitations.Due to the small condition number of this study, minimizing the forward error is equivalent to minimizing the backwards error.
With increased problem size, we observe that the percentage of annealed solutions which return the ground state decreases exponentially.This indicates the solution for a dense matrix may require exponentially more evaluations to obtain for current quantum annealers.The observed scaling is consistent with the assumption that the energy gap exponentially vanishes with increasing size for a dense Hamiltonian.In particular, beyond n = 16, only one out of 100,000 evaluations yield the resulting annealed solution, demonstrating that the real ground-state is well beyond the reach of the available statistics.
Condition number -Fig.3b demonstrates the scaling of the algorithm with respect to changing condition number.The problem size is fixed to a system of 12 linear equations, and a 2-bit precision search is employed.The condition number affects the solution vector x, and therefore for this study we restrict the search range to span exactly the minimum and maximum values of x.The chosen search range keeps the resulting relative residual approximately constant under varying condition number.For linear systems with larger condition numbers, minimizing the forward error is no longer a reliable estimate of the residual of the backwards error, and is therefore dropped from this study.
With increasing condition number, we observe that the percentage of solutions that converge to the lowest-lying state is a relatively constant value as demonstrated by a less than one order-of-magnitude change between the different examples.This behavior is in stark contrast with the scaling observed in Fig. 3a, and suggests that with increasing condition number, the ground state is exponentially easier to identify.This is in amazing contrast to the classical result from Fig. 2, in which convergence to the solution decreases as condition number is raised.
Precision of search -Fig.3c explores the behavior of the algorithm as the search precision is increased for a system of four linear equations and a condition number of 1.1.We observe that the relative residual exponentially decreases, as expected due to sampling an exponential number of solutions.In contrast however, increasing the search precision linearly increases the size of the target Hamiltonian, and results in requiring exponentially more evaluations from the annealer in order to resolve the ground state.Similarly to Fig. 3a, we observe that the forward error for problem sizes beyond 5-bits of search precision starts to deviate from the backward error, an indication that the limits of hardware control have been reached.
Iterative approach -Finally, we explore the possibility of iteratively applying the algorithm in order to decrease the relative residual of the final solution.We demonstrate this technique on a system of 4 linear equations, with a condition number of 1.1, and perform a 4-bit search for each unknown variable.For this study, we initiate the search space ranging from -1 to 1.At the end of each iteration, we narrow the optimization to two neighboring values of the result allowed by the search space.Fig. 4 shows how the search space is refined with each iteration of the algorithm and converges to the conjugate gradient solution.Fig. 3d shows that the relative residual exponentially decreases with the application of each iteration, while the number of anneals required to sample the ground state stays relatively constant.The solution from quantum annealing at the final (ninth) iteration agrees with conjugate-gradient at single precision accuracy.

Conjugate Gradient
The method of conjugate gradient [33] is an iterative algorithm that builds an optimal polynomial for the solution within an order-q Krylov subspace.
The first basis vector in this procedure p 0 is chosen to be the negative gradient of 1 2 x † Ax − x † b, where † refers to the Hermitian conjugation.All other basis vectors p i =0 will be conjugate to p 0 , thereby satisfying p † i Ap 0 = 0, thus the name Conjugate Gradient.This may be used to solve positive-definite symmetric systems or to solve the normal equations of a non-symmetric system (A † A = A † b).Deploying conjugate gradient on the normal equations was used efficaciously in [34] to compute the weak axial coupling of the nucleon with the Mobius Domain Wall action [35], whose spectrum is indefinite.

Quantum Annealing
We may reduce the QUBO problem of Eq. ( 5) to finding the ground-state eigenvalue and eigenvector of an encoded Hamiltonian using quantum annealing.Quantum annealing uses adiabatic evolution under a time-dependent Hamiltonian, for example, the linear interpolation over s = t/T as where the QUBO function is expressed by the target Hamiltonian H problem , while the quantum state of system is initialized as the previously well-characterized ground state of the Hamiltonian H init .Under idealized conditions, purestate evolution from time s = 0 to s = 1 transforms the initial quantum state to the sought-after ground state of the problem Hamiltonian.The minimum time T required to evolve the time-dependent Hamiltonian for values of s ∈ [0, 1] which yields a vector Ψ that is close to the ground state of H( 1) is given by the adiabatic condition [36] T where the numerator is given by the Frobenius norm, and ∆(H(s)) is the energy gap between the instantaneous ground-state and first excited-state eigenvalues.As a result, the time complexity of adiabatic quantum computing depends on the scaling of the eigenvalue spectrum of H(s).However, the time-dependent spectrum for a generic H(s) is a priori unknown and experimental implementations must use a value for T that meets desired tolerances.Because of the absence of detailed Hamiltonian information, quantum annealing is assumed to operate under a relaxation of the condition for adiabatic time-evolution.This relaxation permits coupling of the lowest-lying energy level with higher-lying states, which compromises the promise of terminating in the ground state.In addition, quantum annealing operates at non-zero temperature and may be influenced by external perturbations such as fluctuating control fields and poorly-characterized Hamiltonian interactions.As a result, the ground state is typically prepared with non-unit probability by quantum annealing, and the process must be repeated multiple times so that a distribution of solutions representing possible eigenvalues and eigenstates is generated.Notably, QA often resolves the low-lying spectrum of H problem , and for a generic problem, the interesting statistic is to measure the percentage of measurements which resolve the ground-state, and serves as the QA analogue of Eq. (23).
We use quantum annealing to solve systems of equations by mapping the corresponding QUBO to an Ising model Hamiltonian.Due to hardware topology, we restrict our applications to only linear systems.We validate the algorithm using the D-Wave 2000Q commercial quantum annealer.This hardware is based on cryogenically cooled superconducting electronic elements that implement a programmable Ising model.Each quantum register element expresses a single Ising spin variable, but the D-Wave 2000Q supports only a limited connectivity between these elements.The underlying connectivity pattern is represented by the Chimera graph shown in Fig. 5 [37].In particular, the i-th spin variable may be assigned a bias h i = Q ii and can be coupled to a unique set of six neighboring registers through the coupling J ij = Q ij .A densely connected Hamiltonian can be embedded into the sparsely connected Chimera topology by using secondary constraints to build chains of strongly correlated elements in which J constraint ij J Ising ij .This coupling constraint favors chains of spin elements which behave as a single spin variable [38].Previous studies have identified the optimal mappings for a fully-connected graph into the Chimera structure [39,40].For the D-Wave 2000Q quantum annealer, approximately 64 logical spin variables may be represented within the 2048 physical spin elements.Our subsequent examples use the dwave-sapi2 Python library [41], which is a software tool kit that facilitates cloud access to the 2000Q and supports a heuristic embedding method for the available hardware.

Figure 1
Figure1|Histogram and cumulative distributions for solutions to the linear least squares example problem.The solutions are discrete by construction and no binning is performed to create the histogram.A total of 100,000 measurements are performed on the quantum annealer.

Figure 2
Figure2|The number of conjugate gradient iterations grows slightly worse than κ(A).The stopping criterion is a tolerance of 10 −6 for the norm of the relative residual.All matrices are rank 12, with smaller eigenvalues as κ(A) increases, but identical eigenvectors.The same right-hand side is solved for all cases.

Figure 3 Figure 4
Figure3| a, (Top) The black dashed line shows the theoretical minimum relative residual of the algorithm as predicted by minimizing the forward error, given the search precision and search range used for the test.The red crosses are results of the lowest energy state from 100,000 quantum annealing measurements.Physical measurements deviate from the theoretical minimum as problem size grows.(Bottom) The corresponding percentage of measurements in the minimum energy state are shown in blue.The vertical axis is shown in a logarithmic scale.b, Analogous to Fig.3abut for varying condition number.The forward error prediction (dashed black line) is omitted and not a reliable measure of the relative residual for larger condition numbers.Note that the vertical axis of the bottom plot showing the percentage of measurements observed in the lowest-lying state is on a linear scale.c, This plot is analogous to Fig.3abut for varying search precision.Note that both vertical axes are on a logarithmic scale.d, (Top) The relative residual exponentially decreases with each iteration of the algorithm.By the ninth iteration the result reaches single precision.(Bottom) The percentage of quantum annealing solutions in the lowest-lying state.The algorithm successfully resolves the solution at single precision for this example without issue.

Figure 5 |Example
Figure 5 |Example Chimera graph for a D-Wave 2000Q Quantum Computer solving a system of 12 linear equations.The lattice sites are qubits, the values on the sites set the biases hi, and the connections set the couplers Jij.The figure is generated with the online solver visualizer provided by D-Wave qubist [42].