Quantum annealing for systems of polynomial equations

TL;DR

This paper introduces a quantum annealing-based direct method for solving systems of polynomial equations, demonstrating its application to linear regression and analyzing its scaling behavior and iterative refinement capabilities.

Contribution

It presents a novel quantum annealing approach for directly solving polynomial systems, validated on real hardware, with analysis of scaling and iterative improvements.

Findings

Successfully solved second-order polynomial equations on a quantum annealer.

Demonstrated application to linear regression problems.

Showed iterative annealing can achieve high-precision solutions.

Abstract

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}$.

Equations

$$F_i=P_i^{(0)}+\sum _j{P}_{ij}^{(1)}{x}_j+\sum _{jk}{P}_{ijk}^{(2)}{x}_j{x}_k+\cdots =0$$

$$x_j=a_j\sum _{r=0}^{R-1}{2}^r{\psi }_{rj}+b_j.$$

$$\mathcal{A}$$

$$\mathcal{B}$$

$$\mathcal{R}$$

$$\mathcal{X}$$

$$\Psi$$

$${\chi }^2=$$

$$\equiv$$

$$Q=$$

$$\mathrm{independent}\{x_i:i\in \{1,...,X\}\}$$

$$⟨y_i⟩=\frac{1}{N}\sum _{g=1}^N{y}_{i;g}$$

$$\sum _{ij}{\left[F(x,p)-⟨y⟩\right]}_i{S}_{ij}^{-1}{\left[F(x,p)-⟨y⟩\right]}_j$$

$$F(x_i,p)=\sum _{n=1}^P{p}_n{f}_n(x_i)$$

$$\sum _{ij}\left[\sum _n{p}_n{f}_n(x_i)-y_i\right]S_{ij}^{-1}\left[\sum _m{p}_m{f}_m(x_j)-y_j\right].$$

$$P^{(1)}=\left(\begin{array}{ccc}{f}_0(x{)}^T{S}^{-1}{f}_0(x)& \dots & f_0(x{)}^T{S}^{-1}{f}_P(x)\\ ⋮& \ddots & ⋮\\ f_P(x{)}^T{S}^{-1}{f}_0(x)& \dots & f_P(x{)}^T{S}^{-1}{f}_P(x)\end{array}\right)$$

$$0=2x_0^2+3x_0{x}_1+x_1^2+2x_0+4x_1-51,$$

$$(x_0,x_1)=(2,3),$$

$$P^{(0)}=\left(\begin{array}{c}-51\\ -46\end{array}\right),$$

$$\Psi =\left(\begin{array}{cccc}0& 1& 1& 1\end{array}\right)\to \mathcal{X}=\left(\begin{array}{cc}2& 3\end{array}\right)$$

$$C({\psi }_i{\psi }_j-2{\psi }_i{\psi }_{ij}^a-2{\psi }_j{\psi }_{ij}^a+3{\psi }_{ij}^a)$$

$$\Psi =$$

$$=$$

$$\mathrm{E}\left[D(x)\right]=8+4x+7x^2$$

$$F(x,A)=A_0+A_{1}x+A_2{x}^2,$$

$$A_i={\psi }_{1i}+2{\psi }_{2i}+4{\psi }_{3i}+8{\psi }_{4i}$$

$${\Psi }_0=$$

$$A_0=\left(0001\right)\to 8$$

$$\kappa (P^{(1)})=\frac{{\sigma }_{\max }(P^{(1)})}{{\sigma }_{\min }(P^{(1)})}$$

$${\chi }_{\text{rel.}}^2=\frac{(P^{(1)}{x}_{\mathrm{approx}}+P^{(0)}{)}^2}{P^{(0{)}^2}}=\frac{E_0}{P^{(0{)}^2}},$$

Full text

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, while polynomial systems are typically solved iteratively with homotopy continuation. 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 work, we present an approach for solving a general system of $n^{\textrm{th}}$-order polynomial equations based on the principles of quantum annealing, followed by a demonstration of the algorithm for a system of second-order polynomial equations on commercially available quantum annealers. We then narrow the scope to examples of linear equations by first demonstrating an application to linear regression, before elucidating results on ill-conditioned linear systems motivated by the discretized Dirac equation $D\phi=\chi$ from lattice quantum chromodynamics (QCD). The solution to the discretized Dirac equation is currently the only approach for evaluating non-perturbative QCD. However, well-known numerical challenges slow convergence with conventional solvers [6, 7]. We end by using quantum annealing to solve a similar system and characterize the performance from experimental demonstrations with a commercial quantum annealer.

Polynomial Systems of Equations.

We consider the system of $N$ polynomial equations

[TABLE]

where $i\in\{1,\dots,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^{(0)}_{i}+\sum_{j}P^{(1)}_{ij}x_{j}=0$.

Prior to this work, there exists no direct methods for solving a general $n^{\textrm{th}}$-order polynomial system. For linear systems, existing approaches 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. Preconditioning improves convergence of linear systems by transforming the input as $M^{-1}P^{(1)}x=M^{-1}b$, where the preconditioner $M$ must be inexpensive to invert and $M^{-1}$ should be “close" to $P^{(1)^{-1}}$, so that $M^{-1}P^{(1)}$ resembles a matrix close to unity. Identifying an effective preconditioner plays an important role in numerical convergence of iterative methods [8, 9, 10, 11, 12]. For lattice QCD applications [13], 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 and requires the use of additional iterative methods that also suffer from critical slowing down. Solutions to address this issue include EigCG [6, 7], inexact deflation [8], and adaptive multigrid [9, 10, 11, 12].

RESULTS

Quantum Annealing for Polynomial Solvers

Quantum annealing offers an alternative approach to solving a general system of equations. We map each variable $x_{j}$ using $R$ number of qubits such that

[TABLE]

where $\psi_{rj}\in\{0,1\}$, $a_{i}\in\mathbb{R}$ and $b_{i}\in\mathbb{R}$ such that $x_{j}\in\{b_{j}+2^{r-1}a_{j}|r\in\mathbb{Z}^{\geq}<R\}$. Defining the vectors

[TABLE]

where $|V|$ is the cardinality operator yielding the number of elements in a generic vector $V$. The objective function $\chi^{2}$ which solves Eq. (1) is given by minimizing the residual sum of squares in the qubit-basis

[TABLE]

where $\cdot$ is the dot product, $\circ$ is the Hadamard product, and $\otimes$ is the tensor product. In particular, $(\circ\mathcal{A}\otimes\mathcal{R})^{n}\equiv\circ\mathcal{A}^{\otimes n}\otimes\mathcal{R}^{\otimes n}$, where $V^{\otimes n}$ is a repeated $n$ sequence of tensor products. The ground state of Eq. (3) solves a system of polynomial equations. For current commercial quantum annealers, auxiliary qubits are required to reduce multi-linear terms down to bilinear interactions through quadratization [14, 15, 16, 17, 18, 19, 20, 21]. We provide the details of quadratization through reduction-by-substitution on a system of second-order polynomials in Methods.

Finally, we note that resulting energy at the end of the optimization corresponds to exactly the residual sum of squares if the constant terms in $\chi^{2}$ are correctly accounted for. It follows that the entire energy spectrum is positive, and if the exact solution is recovered, then the ground state energy must be zero.

Quantum Annealing for Linear Solvers

A system of linear equations simplifies Eq. (3) to involve only bilinear terms without quadratization, and reduces to a quadratic unconstrained binary optimization (QUBO) problem where $H^{\mathrm{QUBO}}(\psi)=\sum_{ij}\psi_{i}Q_{ij}\psi_{j}$ with

[TABLE]

where $P^{\prime}_{n}=P^{(0)}_{n}+\sum_{i}b_{i}P^{(1)}_{ni}$. In addition, constant terms that arise from the substitution of Eq. (2) are omitted for simplicity and leaves the solution vector $\Psi$ unchanged, but should be included when interpreting the energy as the residual sum of squares.

Application to Linear Regression

Given a set of $N$ identical and independent observations of the

[TABLE]

variable, the mean and covariance of $y_{i}$ follows

[TABLE]

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\in\{1,...,P\}\}$, and a corresponding objective function for generalized least squares may be defined as

[TABLE]

where the optimal value for the set $p$ is determined by minimizing Eq. (6).

Restriction to linear least squares demands that the fitting function is linear in the unknown parameters, and therefore may be written in the form

[TABLE]

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

[TABLE]

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 $\sum_{j}A^{(1)}_{ij}p_{j}=A^{(0)}_{i}$ analogous to Eq. (1) where

[TABLE]

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

System of Second Order Polynomial Equations

We demonstrate the validity of the algorithm on a system of two second order polynomial equations. The problem is chosen to be small such that the solution can be confirmed by an explicit search over the entire Hilbert space, and evaluated onto a D-Wave annealer. Consider the following system of equations,

[TABLE]

with four real solutions at

[TABLE]

For the sake of discussion, we set up Eq. (3) to solve for the solution at $(2,3)$ by choosing $\mathcal{A}=\begin{pmatrix}1&1\end{pmatrix}$, $\mathcal{B}=\begin{pmatrix}0&0\end{pmatrix}$, and $\mathcal{R}=\begin{pmatrix}2^{0},2^{1}\end{pmatrix}$. The tensors $P^{(n)}$ are obtained by inspection,

[TABLE]

After transforming to the qubit-basis, direct search of the ground state of the 4-body Hamiltonian yields

[TABLE]

where the consecutive pairs of binary variables maps to the binary representation of $x_{i}$ with little-endianness due to specific choices of $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{R}$. The solution is reproduced when transforming the set of generalized $n$-dimensional $\mathcal{Q}\equiv\{Q^{(0)},Q^{(1)},\dots,Q^{(N)}\}$ matrix as defined in Eq. (4), to upper-triangular tensors, and also reproduced when further reducing the dimensionality of elements in $Q^{(n\geq 2)}$ with repeated indicies, yielding in general the most sparse upper-triangular representation $\mathcal{Q}^{\textrm{sparse}}$.

We quadratize the $\mathcal{Q}^{\textrm{sparse}}$ set of rank 0 to $N$ tensors to the QUBO representation with reduction-by-substitution  [14, 15, 16, 21] by introducing $\frac{1}{2}N(N-1)$ auxiliary qubits $\psi^{a}$ to enforce the following constraint,

[TABLE]

such that the constraint is minimized when $\psi^{a}_{ij}=\psi_{i}\psi_{j}$. The coefficient $C$ should be chosen large enough such that the constraint is satisfied under optimization. Additional details of the quadratization used is given in Methods.

We repeat the exercise of solving the same system of polynomial equations on the D-Wave annealer with the symmetrized and quadratized representation, and successfully reproduce the solution,

[TABLE]

where the last six auxiliary qubits confirm consistency of reducing many body interactions down to bilinear terms. Finally, we mention that with reduction-by-substitution, a system of $m^{\mathrm{th}}$-order polynomials needs to be quadratized $m-1$ times, requiring exponentially more auxiliary qubits.

The software implementation of various many-body and 2-body quadratized representations for a system of second-order polynomial equations, the brute force solver, and D-Wave solver are made publicly available [22].

Linear Regression

As an example, consider the following artificially generated data

[TABLE]

where $x\in\mathbb{Z}:x\in[0,49]$. The Toeplitz correlation matrix [23] 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

[TABLE]

and we estimate the parameters $A_{n}$ given the data $D(x)$. Using Eq. (2), we express each parameter $A_{i}$ as a 4-bit unsigned integer

[TABLE]

and construct the problem Hamiltonian following Eqs. (5) and (9). The required 12 logical qubits (3 parameters $\times$ 4-bit representation) support a total of 4096 possible solutions. Explicit evaluation finds the true ground-state to have energy $E_{0}=-1.418$ and eigenstate

[TABLE]

which corresponds to the parameter values

[TABLE]

These correct coefficients for the generating function in Eq. (11) verify the design of the algorithm.

We next test the algorithm by solving the objective function using quantum annealing. The target Hamiltonian of Eq. (11) is solve with a D-Wave annealer, and results for 100,000 independent evaluations are acquired using an annealing schedule with $T=200\leavevmode\nobreak\ \mu$s. The correct result is reproduced in 0.5% of the solves, while the lowest 0.8% of the eigenvalue spectrum is obtained by 10% of the solves with overall results biased towards the lower-lying 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 $P^{(1)}$ is defined as the ratio of maximum and minimum singular values.

[TABLE]

In the case of symmetric matrices, this is equivalent to the ratio of largest and smallest eigenvalues.

We vary our test matrices in two ways: 1) vary the problem size while holding the condition number fixed, 2) the problem size is held constant with varying condition number. The accurately of the solution will be judged by the relative residual sum of squares,

[TABLE]

where $E_{0}$ is the ground-state energy. For conjugate gradient, a tolerance for Eq. (16) 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_{\operatorname{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 $N=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. 1 shows slightly worse than square root scaling of conjugate gradient with condition number.

The matrices from this and subsequent results are constructed by creating a random unitary matrix of rank $N$, denoted as $U$. A diagonal matrix $\Lambda$ is then linear populated by evenly spaced real-valued eigenvalues, such that $\mathrm{max}(\Lambda)/\mathrm{min}(\Lambda)=\kappa$. The matrices are trivially formed as $P^{(1)}=U\Lambda U^{\dagger}$. A common right-hand side is taken for all $P^{(1)}$: a vector of length $N$ 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 – In Fig. 2a, we study the scaling behavior for $\kappa=1.1$ and $R=2$. Due to prior knowledge of the conjugate-gradient solution, the search space for all $N$ parameters are fixed for the set of problems, and 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. 2b demonstrates the scaling of the algorithm with respect to changing condition number. The problem size is fixed to $N=12$, and $R=2$. 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. 2a, 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. 1, in which convergence to the solution decreases as condition number is raised.

Precision of search – Fig. 2c explores the behavior of the algorithm as $R$ is increased for $N=4$ and $\kappa=1.1$. We observe that the relative residual exponentially decreases, as expected due to sampling an exponential number of solutions. However, increasing $R$ also requires exponentially more evaluations from the annealer in order to resolve the ground state. Similarly to Fig. 2a, we observe that the forward error for problem sizes beyond $R=5$ 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 $N=4$, with $\kappa=1.1$, and $R=4$. For this study, we initially set $\textrm{min}(\mathcal{X})=-1$ and $\textrm{max}(\mathcal{X})=1$. With each iteration, we narrow the optimization to two neighboring values of the result allowed by the search space. Fig. 3 shows how the search space is refined with each iteration of the algorithm and converges to the conjugate gradient solution. Fig. 2d 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.

METHODS

Quadratization

Quadratization (i.e. to make quadratic) maps terms in the Hamiltonian that are multi-linear with respect to the binary variables $\mathcal{X}$, to a larger Hilbert space involving only bilinear (i.e. quadratic) contributions. This transformation is required to realize quantum algorithms with non-linear operations on near-term quantum computers. There exists in current literature, a rich selection of methods to perform such a task [14, 15, 16, 17, 18, 19, 20, 21]. For this work, we apply reduction-by-substitution [14, 15, 21], where the constraint equation is given by Eq. (10).

One constraint equation is required to define each auxiliary qubit $\psi_{ij}$. After taking into consideration that in the qubit basis, the Hamiltonian is symmetric under permutations of all indices, $\frac{1}{2}N^{\prime}(N^{\prime}-1)$ auxiliary qubits are needed to account for every unique quadratic combination of the underlying basis of length $N^{\prime}$. In Fig. (4), we provide the smallest non-trivial example which maps a system of two second-order polynomial equations (i.e. $N=2$), with $R=2$ after the $n$-body Hamiltonian has been reduced to the set $\mathcal{Q}^{\textrm{sparse}}$. The constraint equations have an overall coefficient $C$ which needs to be large enough such that the constraints are satisfied under optimization. For visual clarity, the constraints in the second quadrant are entered in the lower triangular section, but in practice should be accumulated with the upper triangular section occupied by $Q^{(2)}$.

Generalization to a system of more equation ($N>2$), or finer searches ($R>2$) is straightforward, and implemented in the accompanying software [22]. Generalizing to higher-order polynomial equations required additional levels of quadratization that was not implemented in this study.

Quantum Annealing

Akin to adiabatic quantum optimization [24], quantum annealing prepares a quantum statistical distribution that approximates the solution by applying a slowly changing, time-dependent Hamiltonian [25, 26, 27], 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 [26, 28, 29, 30, 31, 32, 33, 34]. While examples of non-trivial advantages have been observed for fixed-size problem instances [35, 36, 37, 38, 39], more general statements about computational complexity remain unresolved [40].

We demonstrate the proposed 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. In particular, the $i$-th spin variable may be assigned a bias $Q_{ii}$ and can be coupled to a unique set of six neighboring registers through the coupling $Q_{ij}$. A densely connected Hamiltonian can be embedded into the hardware by using secondary constraints to build chains of strongly correlated elements in which $Q_{ij}^{\mathrm{constraint}}\gg Q_{ij}^{\mathrm{problem}}$. This coupling constraint favors chains of spin elements which behave as a single spin variable [41]. Previous studies have identified optimal mappings of the infinite dimensional to three-dimensional Ising model [42, 43]. For the D-Wave 2000Q, approximately 64 logical spin variables may be represented within the 2048 physical spin elements. Our examples use the dwave-sapi2 Python library [44], which is a software tool kit that facilitates cloud access to the annealer and supports a heuristic embedding method for the available hardware.

DATA AVAILABILITY

Software to reproduce the solutions of the second-order system of polynomial equations discussed in this work is made publicly available at https://github.com/cchang5/quantum_poly_solver.

ACKNOWLEDGMENTS

We thank Hidetoshi Nishimori (Tokyo Tech), Tetsuo Hatsuda (RIKEN), and Chih-Chieh Chen (ITRI) for discussions. This work was performed under the auspices of the U.S. Department of Energy by LLNL under Contract No. DE-AC52-07NA27344 (AG). Access to the D-Wave 2000Q computing system was provided by Oak Ridge National Laboratory. TSH acknowledges support from the Department of Energy, Office of Science, Early Career Research Project and the ORNL Directed Research and Development funds. This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

AUTHOR CONTRIBUTIONS

Initial idea was proposed by C.C.C.. Design of test cases was done by A.G. and C.C.C.. Conjugate-gradient solves was done by A.G.. Calculation done on the DWave 2000Q was performed by C.C.C. and T.S.H.. Cross checks with density matrix renormalization group was performed by S.S. and C.C.C.. All authors contributed to writing and editing of the final manuscript.

ADDITIONAL INFORMATION

Competing Interests: The authors declare no competing interests.