{'fraction_non_alphanumeric': 0.030644202180376612, 'fraction_numerical': 0.009791873141724479, 'mean_word_length': 4.516291274874262, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A hybrid classical-quantum approach to solve the heat equation using quantum annealers

The numerical solution of partial differential equations by discretization techniques is ubiquitous in computational physics. In this work we benchmark this approach in the quantum realm by solving the heat equation for a square plate subject to fixed temperatures at the edges and random heat sources and sinks within the domain. The hybrid classical-quantum approach consists in the solution on a quantum computer of the coupled linear system of equations that result from the discretization step. Owing to the limitations in the number of qubits and their connectivity, we use the Gauss-Seidel method to divide the full system of linear equations into subsystems, which are solved iteratively in block fashion. Each of the linear subsystems were solved using 2000Q and Advantage quantum computers developed by D-Wave Systems Inc. By comparing classical numerical and quantum solutions, we observe that the errors and chain break fraction are, on average, greater on the 2000Q system. Unlike the classical Gauss-Seidel method, the errors of the quantum solutions level off after a few iterations of our algorithm. This is partly a result of the span of the real number line available from the mapping of the chosen size of the set of qubit states. We verified this by using techniques to progressively shrink the range mapped by the set of qubit states at each iteration (increasing floating-point accuracy). As a result, no leveling off is observed. However, an increase in qubits does not translate to an overall lower error. This is believed to be indicative of the increasing length of chains required for the mapping to real numbers and the ensuing limitations of hardware.

The numerical solution of partial differential equations by discretization techniques is ubiquitous in computational physics. In this work we benchmark this approach in the quantum realm by solving the heat equation for a square plate subject to fixed temperatures at the edges and random heat sources and sinks within the domain. The hybrid classical-quantum approach consists in the solution on a quantum computer of the coupled linear system of equations that result from the discretization step. Owing to the limitations in the number of qubits and their connectivity, we use the Gauss-Seidel method to divide the full system of linear equations into subsystems, which are solved iteratively in block fashion. Each of the linear subsystems were solved using 2000Q and Advantage quantum computers developed by D-Wave Systems Inc. By comparing classical numerical and quantum solutions, we observe that the errors and chain break fraction are, on average, greater on the 2000Q system. Unlike the classical Gauss-Seidel method, the errors of the quantum solutions level off after a few iterations of our algorithm. This is partly a result of the span of the real number line available from the mapping of the chosen size of the set of qubit states. We verified this by using techniques to progressively shrink the range mapped by the set of qubit states at each iteration (increasing floating-point accuracy). As a result, no leveling off is observed. However, an increase in qubits does not translate to an overall lower error. This is believed to be indicative of the increasing length of chains required for the mapping to real numbers and the ensuing limitations of hardware.

I. INTRODUCTION
In fields of computational physics, the mathematical description of the problem of interest relies on the solution of partial differential equations restricted to boundary conditions and initial conditions. Most of the techniques are based on discretization methods, where the physical domain is divided into cells or volumes. When implicit methods are chosen, the discrete set of coupled equations to be solved are cast as a system of linear equations of the form Ax = b, where A is generally a N × N sparse matrix and the vector b is known, while the vector x is the sought solution to the problem. The largest fraction of computational time is spent on the solution of the linear systems.
After the advent of the Harrow, Hassidim, and Lloyd (HHL) algorithm [1] to solve the linear quantum problem, subsequent work has used it as subroutine to solve ordinary or partial differential equations [2][3][4][5][6][7][8][9]. However, the solutions for the differential equations were presented in form of a vector state, restricted to the calculation of global properties of the system of interest, otherwise the cost to obtain each amplitude of the quantum state would render the quantum approach too expensive. Clader et al. [3] argue that an exponential speedup can be achieved by using the HHL algorithm to solve the electromagnetic scattering cross-section problem via the finite element method. However, Montanaro and Pallister [6] showed that the relation between the size of the system of linear equations and the accuracy of the solution was missed, which implies that the speedup can be reduced or even eliminated. A thorough analysis of the resources needed to solve the electromagnetic scattering cross-section of a 2D target was made by Scherer et al. [10], concluding that the number of resources must be reduced by many orders of magnitude for the algorithm to become practical. There are other issues concerning quantum state preparation and the necessity of an oracle to encode the coefficients of the matrix A. For reasons like these, it is not yet clear if exponential speedups have been achieved by quantum algorithms using the linear quantum problem as a subroutine. Recently, Linden et al. [11] showed that HHL-like algorithms used to solve the heat equation are never faster than the best classical algorithms. However, a quantum alternative which uses amplitude amplification can accelerate the process, achieving at most a quadratic speedup. Although this a very encouraging result, the gate model of quantum computing needs advance to approach large problems. Following this reasoning, quantum annealers can be more attractive at the current stage of quantum computing due to their large number of qubits and increase connectivity. For instance, a graph-coloring-based methodology to solve differential equations through quantum annealers was proposed in [12].
An alternative method for solving the linear system of equations by quantum annealing was presented in Refs. [13] and Rogers and Singleton Jr [14]. An analysis of the resources employed for this task was made in Ref. [15], showing that in some circumstances it is possible to have advantage when compared to traditional algorithms to solve least squares problems. The linear system is encoded in a QUBO (quadratic unconstrained binary optimization) problem [16] through a well-defined protocol, whose ground state is the solution to our problem of interest. As opposed to a HHL-type algorithm, the final solution of the linear system of equations is a classical state codifying the vector x. In this work we solve the steadystate 2D heat equation by a hybrid classical-quantum algorithm. Applying the finite difference method to the partial differential equation yields a system of coupled linear equations. In order to enable the solution of larger linear systems of equations, we use the block Gauss-Seidel method and solve the full system iteratively, where the smaller subsystems of equations are solved in the quantum processing units (QPU). The quantum algorithm is implemented in the 2000Q and Advantage systems of D-Wave Systems Inc.
The paper is divided as follows. In §II we introduce the D-Wave quantum annealing methodology, review the approach of Rogers and Singleton Jr [14] to solve the linear system of equations, present the block Gauss-Seidel method, and elaborate about quantum resources required to implement the algorithm. The results are presented in §III and are followed by conclusions and perspectives in §IV.

II. METHODS
As described above, the basic idea of our hybrid approach is to replace a classical subroutine to solve the system of linear equations that result from the discretization of partial differential equations by a quantum algorithm. In the proposed method the system of linear equations is solved using quantum annealing, where the final Hamiltonian encodes a QUBO problem [14]. In the following we briefly review the quantum annealing followed by the QUBO problem.
A. D-Wave quantum annealing D-Wave quantum annealing is based on the adiabatic quantum evolution of the transverse Ising Hamiltonian with the system in thermal equilibrium at a very low temperature [17]. The goal of the computation is to minimize an objective function (energy) for which the solution to the desired problem is approximately encoded in its ground state [18]. In order to achieve the ground state of the system with a high success rate, the quantum fluctuations must be dominant over thermal fluctuations [17].
The idea behind this model is the preparation of the system as the ground state of an easy-to-prepare Hamiltonian, for instance ( ω = 1), where σ i x is the Pauli matrix in the x direction for the i-particle and its ground state is |ψ(0) = |+ ⊗n , which contains the superposition of all states in the computational basis equally distributed. By choosing the final Hamiltonian H(1) as the operator whose ground state contains the solution to the desired problem, it is possible to interpolate between these two Hamiltonians with a schedule function f (s), subject to f (0) = 0 and f (1) = 1, resulting in, where s = t/T (0 ≤ s ≤ 1) is the dimensionless time, t is the current time and T is the total time of the computation. The final Hamiltonian implemented by D-Wave is the Ising Hamiltonian, with the parameters α i and β ij being chosen to encode the solution to the problem and σ i z is the Pauli matrix in the z direction for the particle i.
According to the adiabatic theorem [19], if the evolution is performed adiabatically, then the ground state of the problem Hamiltonian H(1) is found with high probability, provided that the evolution time T from the initial Hamiltonian H(0) to H(1) is proportional to g −2 min , where g min = E 1 − E 0 is the minimum energy gap between the ground and first excited states [19,20]. It is important to observe that a suitable choice for the schedule function may result in totally different evolution times [20][21][22]. Therefore, the final Hamiltonian H(1) must contain the solution to the linear system of equations [14]. Next, we show how to build H(1), thus casting the system of linear equations into a QUBO problem.

B. QUBO algorithm to solve linear systems of equations
O'Malley and Vesselinov [13] and Rogers and Singleton Jr [14] showed how to map the problem of solving a linear system of equations into a QUBO problem [16]. Also, in Ref. [23], the linear system of equations can be obtained as a particular case of a nonlinear system of equations. Given a linear system Ax = b, the solution x can be found by minimizing the function We can approximate each component x i (i = 1 . . . N ) of the vector x by its binary representation in a given , with R bits, as: where each q i r ∈ {0, 1} is the r-th bit in the binary representation of x i . If d i > 0 and c i > d i /2, the domain of x i contains positive and negative real values. This can easily be generalized to complex variables [14]. After substituting (5) into (4) and dropping the positive constant term, the equation becomes, where the coefficients are given by, Here, A ij are coefficients of the matrix A and b i are components of the vector b. It is important not to confuse b i with b ij rs . In quantum annealers, the physical qubits are indexed by a 1-dimensional linear index (see Eq. (3)), while the logical qubits are indexed by 2-dimensional indices i and r, where i = 0, 1, ..., N − 1 and r = 0, 1, ..., R − 1, as shown in Eq. (6). The mapping between the logical an physical indices is given by, such that the inverse mapping is, The Finally, the quantization of the QUBO problem is made throughq |q = q |q , where the operators are such thatq 2 =q , with eigenvalues q ∈ {0, 1} and |q = |q 0 ⊗ ... ⊗ |q N R−1 . In order to solve the system of linear equations, the coefficients a i r → a and b ij rs → b κ , described in Eqs. (7) and (8), must be provided to the D-Wave Ocean Software [24], the interface that communicates with the quantum hardware.
The readout of the quantum algorithm is a bit string q 0 q 1 ...q N R−1 of N.R qbits that returns to the classical algorithm for further processing. In principle, the qubits must be fully connected. However, in situations where the matrix A is sparse, several coefficients b ij rs are set to zero, hence not all connections between the qubits are necessary. With knowledge of the above, if the full linear system of equations is solved directly, the number of required qubits N.R for equivalent floating-point accuracy of classical algorithms makes this implementation prohibitive at the current stage of technological development. However, the linear system can be divided into subsystems, and the full linear system solved iteratively, in a hybrid quantum-classical approach, as detailed next.

C. Iterative QUBO method for solving linear systems
The block Gauss-Seidel method consists in solving the linear system Ax = b iteratively, by applying the Gauss-Seidel method on the original system divided in blocks [25]. Suppose, for simplicity, that the size of the system is an even number N and divide it into two subsystems, This translates to In the iterative method we make an initial guess for x 2 , such as x  (12) and (13) allow us to set up an iterative process where the initial guess is updated by solving two linear systems of size N/2.
At the k-th step, If the sequences x x (k−1) 2 and the original linear system is satisfied. Supposing A 11 and A 22 are invertible, it can be shown that, Hence, a sufficient condition for the convergence of this method is: The generalization of the method is possible by dividing the original linear system into several blocks whose size is adjustable to the particular interest. If the system is divided into D blocks, each smaller system has size N/D × N/D. The number of required iterations is not known beforehand, but defined by a convergence criterion. A typical criterion depends on the normalized residual r (k) , defined by where x (k) is the solution of the system of linear equations after k iteration steps. A solution is considered converged when r (k) meets a specified tolerance τ , such that r (k) ≤ τ . Here, we instead monitor the relative error of the solution vector at the k-th iteration x (k) , given as, It can be shown that e (k) and r (k) are related according to, where κ(A) is the condition number of the coefficient matrix A. Preconditioning techniques should be used when A is ill-conditioned. In this work we do not pursue this issue further, since it would be relevant for both classical and quantum algorithms.

Scaling of the method
The representation of a vector x with N components requires N · R qubits, where R is the number of qubits that gives the desired accuracy of floating-point numerical representation, , of each variable x i . However, as stated above, it is possible to solve smaller subsystems of size N/D by partitioning the full system in D parts. In this case the number of qubits scales as (N/D) · R. In order to reduce the number of qubits required in the computation for a given floating-point accuracy of x i , one possibility is to shrink the interval [−d i , 2c i − d i ] at each iteration, which we adopt here. Any reduction in the number of required qubits is significant in the era of noisy intermediate-scale quantum (NISQ) computers, where at most few thousands of qubits are available. Even more dramatic is the number of connections necessary to solve the QUBO problems. In the worst case, every logical qubits must be fully connected (see the coefficient b ij rs in Eq. (7)). The block Gauss-Seidel method can reduce the number of connections among qubits by a factor N 2 R 2 1 − 1/D 2 . Additionally, when A is sparse, connectivity is reduced considerably.
If is a bound for the accuracy of the floating-point numerical representation of a physical variable x i , for a solution represented with R qubits, then where the solution is represented and 2 R is the number of different points that can be represented with R bits inside the aforementioned interval. From the expression above, the number of qubits required to represent one variable with floating-point accuracy scales as Since the total computing time T depends on the energy gap between the ground and first excited states of the system, it is hard to obtain precise estimates of T . As noted in [26], the task of analytically extracting the minimum gap scaling has been extremely difficult in practice, except for a few special cases, such as for Deutsch-Jozsa [27] and Grover algorithms [21,28].
In Table I we summarize the consuming quantum resources used in the iterative and non-iterative QUBO algorithms to solve the linear system of equations. In order to assess the viability and performance of the proposed method, we solve the heat equation for a square plate of length L subject to fixed temperature at its edges and multiple heat sources and sinks within the domain. The steady state solution T (x, y) is then determined with the following PDE and boundary conditions,   T (x, 0) = 0 • C , x ∈ [0, L] The problem is solved numerically by the finite difference method. The domain is divided into m segments of equal length, with nodes located at positions x i = (L/m) · i and y j = (L/m) · j, where i, j = 0 . . . m. The points (x i , y j ) of the grid corresponding to i = 0, i = m, j = 0 or j = m form the boundary of the plate and are determined by the boundary conditions (26). The points inside the grid are the unknowns, and the partial derivatives can be written as for some Ignoring terms proportional to (L/m) 2 or higher orders, the PDE becomes, where i, j = 1, 2, . . . , m−1. Equation (28) can be written in the form of a linear system Ax = b, where A is a sparse matrix with few non-zero entries and b is defined by the boundary conditions and eventual source terms. The spatial discretization forms a 9 × 9 matrix (m = 11), resulting in a system of coupled linear equations with 81 unknown variables, which is solved by the iterative QUBO algorithm on the Advantage and D-Wave 2000Q systems [29]. A sample solution of the problem is shown in Fig. 1.
This approach enables us to solve larger linear systems of equations, circumventing the limited number of qubits in the D-Wave quantum processing units. For sake of comparison, in this work we use at most 7 precision qubits, which demands at most 56 logical qubits to describe all unknown variables in the block Gauss-Seidel approach, while 567 would be required to solve the full system of linear equations non-iteratively. We will see that the number of logical qubits is much smaller than the number of physical qubits in both D-Wave systems, because of the limited connectivity between physical qubits.
In Fig. 2 the dependence of the relative error defined in Eq. (23) is shown as a function of the number of iterations of the block Gauss-Seidel method, for different numbers of qubits used in the numerical representation of the unknown variables. In this example, the smallest error is found using three qubits in the Advantage system. We also notice that the error does not continually decrease, as would be expected in a classical algorithm. This occurs in part because of the limited floating-point accuracy of the variables x i , such that after a certain number of iterations, an improvement of the solution is not possible. Another observation is that the increase in qubits does not necessarily translate into smaller errors. This result was not expected. However, it can be related to the appearance of more noise in the system, since there are more qubits involved in the numerical representation of a given variable x i , which in turn increases the the probability of chain breaks. This subject will be discussed in the following. In general, the Advantage system performs marginally better than 2000Q.
In order to create a logical qubit, generally more than one physical qubit is necessary. This is a consequence of the limited number of couplers (connections) between qubits. In the D-Wave quantum processing units, a set of physical qubits that describe a vertex of a graph form a chain. For this chain to behave as a basic quantum processing unit, all of its physical qubits must be in the same state. When this does not occur, the chain is broken and the state is determined by a majority vote. In Fig.  (3) we plot the average chain break fraction as function of the number of iterations and the number of qubits used to represent each value of temperature. The average chain break fraction for each number of precision qubits is greater in the 2000Q than in the Advantage system and is on average independent on the iteration number. This result is expected, since in the Chimera topology of the 2000Q system each qubit is connected to 6 other qubits, while in the Pegasus topology of the Advantage system, each qubit connects to 15 other qubits. Furthermore, the Advantage system has 5000+ qubits compared to 2000+ qubits in the 2000Q system. These features allow for more compact chains that are less susceptible to noise [30]. A significant improvement in the relative error of Fig. 2 can be obtained by reducing the length of the interval 2c i represented by a given number of qubits at each iteration. The iterative QUBO algorithm is modified by shrinking the interval by factor of γ = 0.8 at each iteration, i.e., the interval for estimating the solution at the k-th iteration is given by is the estimate for x i at the (k − 1)-th iteration. This procedure increases the capacity of the QUBO solver to represent the solution of the system with a fixed number R of qubits, as can be seen in Figure 4. The implementation of this updated search interval length brings the quantum solvers on par with the error obtained on classical hardware. It must be noted, however, that the update of the rhs of equations (16) and (17) during the Gauss-Seidel iterations is performed on a classical computer, with double precision, even in the quantum algorithm.
A comparison of Figs. 2 and 4 shows that the 2000Q performs slightly better than the Advantage system when the update of the search interval length is incorporated into the algorithm. There is not a clear interpretation why this occurs, since the Advantage system is a improved version of the 2000Q system, where the chains are more compact and less susceptible to noise [30]. The results obtained from both figures also suggest that a greater number of precision qubits does not necessarily translate to a smaller error, as already noticed earlier.

ACKNOWLEDGMENTS
The authors acknowledge the financial support from the Brazilian funding agencies Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -CAPES, Conselho Nacional de Desenvolvimento Científico e Tecnológico -CNPq, Fundação de Amparoà Pesquisa e Inovação do Estado de Santa Catarina -FAPESC, and the Instituto Nacional de Ciência e Tecnologia de Informação Quântica -INCT-IQ.

IV. CONCLUSIONS AND PERSPECTIVES
The hybrid classical-quantum algorithm we have developed solves partial differential equations by the finite difference method, yielding a system of coupled linear equations. This technique was tested with success for the 2D steady-state heat equation, but in principle can be applied to any partial differential equation, encompassing a large array of applications for solution by quantum annealers. The method with updated search interval length allows us to obtain comparable accuracy to a classical algorithm with double precision using 3 qubits to numerically represent each variable. The main limitation of the quantum hardware is the reduced connectivity between qubits, albeit great strides have been made in order to increase the number of couplers for each qubit. The iterative algorithm allows for a significant reduction in the number of logical qubits, from N.R to N.R/D, at the expense of an increase in the number of required iterations for a given accuracy.
It is our understanding that there are good prospects for the iterative QUBO algorithm and that there is room for further improvements, such as robust heuristics for update of the search interval length. Another improvement can be obtained by using an application-oriented embedding, as many linear systems typical of finite difference and finite volume methods are characterized by sparseness of the coefficient matrix A. This has the advantage of requiring less couplers than a fully connected problem. Possibly a specific embedding can be found for this class of problems, allowing for less qubits and enabling the solution of larger systems of coupled linear equations. Owing to our limited access to the D-Wave QPUs, we believe that more tests are required to fully