Box algorithm for the solution of differential equations on a quantum annealer
Siddhartha Srivastava Veera Sundararaghavan








PhD candidate, sidsriva@umich.edu
 

 Associate Professor, veeras@umich.edu

Department of Aerospace Engineering

University of Michigan Ann Arbor MI 48109





             


Differential equations are ubiquitous in models of physical phenomena. Applications like steady-state analysis of heat flow and deflection in elastic bars often admit to a second order differential equation. In this paper, we discuss the use of a quantum annealer to solve such differential equations by recasting a finite element model in the form of an Ising hamiltonian. The discrete variables involved in the Ising model introduce complications when defining differential quantities, for instance, gradients involved in scientific computations of solid and fluid mechanics. To address this issue, a graph coloring based methodology is proposed which searches iteratively for solutions in a subspace of weak solutions defined over a graph, hereafter called as the `box algorithm.' The box algorithm is demonstrated by solving a truss mechanics problem on the D-Wave quantum computer.


Valid PACS appear here
                             

                              
Box algorithm for the solution of differential equations on a quantum annealer
    Veera Sundararaghavan
    December 30, 2023
==============================================================================





















§ INTRODUCTION


Computational methods are rapidly emerging as an essential tool to understand and solve complex engineering problems complementing the traditional means of experimentation and theory.  Richard Feynman's statement <cit.> “with a suitable class of quantum machines you could imitate any quantum system, including the physical world" has driven our vision towards a machine that can solve computational problems inaccessible to classical computers. Early versions of such quantum computers have already appeared. Mirroring gate-based classical computers, gate–based quantum computers with a small number of qubits have been demonstrated and promise an eventual path towards universal quantum computation. However, noise limits the number of gate operations that can be enforced before the quantum states decohere. In parallel, quantum annealers have been developed that provide a significant number of qubits for solving a class of combinatorial optimization problems. In these machines, an Ising hamiltonian is engineered such that the solution to the computational problem is encoded in its ground state. The system evolves adiabatically to the ground state as governed by the Schrodinger equation for the time-dependent hamiltonian. 

The D-wave system is a quantum annealer that currently provides more than 2000 qubits modeling a transverse Ising hamiltonian whose ground state is NP-complete. The  hamiltonian with q_i describing the state of the i^th qubit is given by: 

    E(𝐪) = ∑_i ∈sites H_i q_i + ∑_(i,j) ∈links J_ij q_i q_j

The hamiltonian includes self-interaction (H_i is the on-site energy of qubit q_i) and site-site interaction terms (J_ij are the interaction energies of two connected qubits q_i and q_j), with the qubits connected in a Chimera graph.  The system is first initialized in the ground state of a hamiltonian which is known and easy to prepare. 
Then, the Hamiltonian is changed such that the system equilibrates to the ground state of the Ising Hamiltonian E(𝐪) according to the adiabatic theorem.
NP-hard combinatorial optimization problems can be encoded through the field and site interaction strengths. The system particularly holds promise for solving graph coloring problems with large sizes (N) where classical polynomial time algorithms cannot be devised.   Many engineering problems in airline scheduling, image segmentation, and pattern recognition have been encoded as graph coloring problems solvable on quantum annealers.  

While differential equations are ubiquitous in models of physical phenomena, the use of quantum annealers for scientific computing in solid and fluid mechanics has not yet been explored. Scientific computing mostly involves solving a linear system of equations Ax=b defined on a continuum domain discretized with finite elements. The matrix A, generally being sparse, structured and positive definite matrix obtained by assembling element-level stiffness matrices. In the past, gate–based quantum computing algorithms have been devised to solve the system of linear equations using QLSA algorithms (HHL algorithm <cit.>) and its variants <cit.>. This algorithm, unlike a classical solver, does not give a direct solution x but rather allows sampling from the solution vector. Nevertheless, this has spawned several works in differential equation modeling on quantum computers (<cit.>). The sampling task by itself requires solving  Ax=b. In the classical setting, the complexity scales with the size of the problem and goes as O(Nsklog(1/ϵ)) for conjugate gradient method where N is the number of unknowns, k is the condition number, s is the sparsity of A and ϵ is the precision of the solution. On the other hand, the QLSA <cit.> has a favorable running time of O(log(N)k^2s^2/ϵ) which scales logarithmically with the size of the problem. Quantum annealers are especially attractive for scientific computing with the ability to scale up the simulations to a more significant number of qubits. However, algorithms for the solution of differential equations have not been devised yet on these systems<cit.>. The merit of the paper is the algorithm to solve a differential equation on an annealer. 

Here, we note that the solution to Ax=b can be encoded in an equivalent minimization problem min (1/2 x^TAx - x^Tb) which contains field and interaction terms similar to an Ising model. Thus, in this paper, we explore mapping of this energy to Ising hamiltonian on the D-wave machine recognizing that the graph in the D-wave chip by itself models a finite element mesh-like topology. The element level stiffness and force vectors are then encoded in the Ising hamiltonian as interaction weights and field variables. Dirichlet boundary conditions are enforced by modifying field terms to favor one qubit state over another. An illustration of this procedure is presented in Fig <ref>. A discretized version of the differential equation is solved using energy minimization on a graph. Direct minimization of energy may hold advantages over the conventional finite element approach in systems which lead to a matrix with large condition numbers, zero or negative eigenvalues leading to bifurcation events such as buckling in shells and phase transitions <cit.>.




In the solution to a differential equation, the qubits must encode a rational number. However, the qubit encoding the Ising lattice point carries two discrete levels (up/down spin) in the ground state. In classical computers, with similar binary (0/1) encoding, anywhere from 32 bits (float) to 80 bits (long double) of memory can be used to encode more than 12 million high precision variables in 1 GB memory. In contrast, currently available quantum annealers have a limited number of physical qubits. This restriction makes the representation of solutions of double precision similar to a classical computer extremely expensive.
In (<cit.>), the problem of minimizing ||Ax-b|| in the least squares sense was posed by encoding physical qubits to represent rational numbers using a radix 2 representation. This format requires a significant number of physical qubits and connections to represent positive rational numbers and an additional qubit to represent the sign of the number (<cit.>). 
In comparison, the box algorithm searches within a small discrete set of up/down qubit values with each element of the set mapped to a double precision value, thereby eliminating the need for additional qubits to achieve higher precision. 

In this paper, we consider a self-adjoint form of second order differential equation as the model problem. The problem statement and the relevant mathematical details are presented in Section <ref>. The Graph representation of the problem is formulated in Section <ref>.  The iterative procedure, referred as 'Box algorithm', is presented in Section <ref>. All procedures are accompanied with numerical examples for elucidation. This algorithm is demonstrated by solving a truss mechanics problem on the D-Wave quantum computer in Section <ref>



§ MATHEMATICAL PRELIMINARIES


A self-adjoint form of a second order differential equation on an interval (x_l, x_r) is defined as, 

    -(p(x)u'(x))' + q(x)u(x) = f(x)      x_l<x<x_r

Dirichlet boundary conditions are considered at both ends i.e.
u(x_l) = u_l  and u(x_r) = u_r.  Well-posedness of this problem requires p(x) ≥ p_min > 0  and  q(x) ≥ q_min≥ 0 
Furthermore, for convenience, it is assumed that p,q ∈ C([x_l,x_r]) and f ∈ L^2([x_l,x_r]). These conditions are sufficient to show the existence of a unique solution to the weak form (<cit.>). 



 §.§ Functional minimization


Motivated by the intractability of direct integration of the differential equation (<ref>), it is often convenient to employ functional minimization techniques. Calculus of variations can be used to observe that the minimization of the functional  (<ref>) leads to the strong form described in Eq (<ref>).

    Π[u] = ∫_x_l^x_r( 1/2pu'^2 + 1/2qu^2 -fu )dx

Square integrability of u  and its first derivative are required in this definition of Π[u]. The implication is that the minimizing solution, u lies in the Sobolov space H^1([x_l, x_r]). 
A discrete problem is obtained by using a finite basis for the solution defined in Eq (<ref>) which satisfies the Dirichlet boundary conditions. The admissible choices of 𝐚=(a_0, a_1, ..., a_N) satisfy u_N(x_d)=u_d where x_d is a Dirichlet boundary and u_d is the prescribed value at that point. This approximation reduces the infinite dimensional functional minimization problem to finite dimensions. The approximated functional Π_N is entirely determined by the representation of u in the finite basis as shown in Eq (<ref>). 
It is worth observing that the choice of ϕ_i(x) is such that ϕ_i ∈ H^1([ x_l, x_r  ]) i.e. for any  u_N ∈ V_N = span{ϕ_1, ϕ_2, ...,ϕ_N }⊆ H^1([ x_l, x_r  ]). Additionally the proper inclusion, V_i ⊆ V_i+1, guarantees convergence of the solution with increasing N.

    u_N(x) = ∑_i=0^N a_i ϕ_i(x)




    Π_N[a_0, ...,a_r,.. ,a_N ] = ∫_x_l^x_rp/2(∑_i=0^N a_i ϕ'_i)^2 
    + q/2(∑_i=0^N a_i ϕ_i)^2 -f(∑_i=0^N a_i ϕ_i) dx


As the solution is completely determined by the variable 𝐚, the functional minimization of Eq (<ref>) is reformalized as Eq (<ref>) where 𝐚^b.a. refers to the coefficients of best approximation of solution, u_N, in the subspace V_N


    𝐚^b.a. = _𝐚Π_N(𝐚)




 §.§ Finite Element approximation

Finite element basis provides a popular choice of compactly supported shape functions. For the purpose of simplicity, `tent/hat functions' (defined in Eq (<ref>)) are used in this work. The domain is split into N elements with N+1 nodes. The generalization to higher order families of piecewise-continuous basis functions is immediate but is omitted for brevity. 


    ϕ_i(x) = {[                                0,                        x < x_i-1,;        (x - x_i-1)/(x_i - x_i-1),                   x_i-1≤ x < x_i,; 1 -
         (x - x_i)/(x_i+1 - x_i),                   x_i≤ x < x_i+1,;                                0,                        x≥ x_i+1 . ].


The usage of compact basis further reduces the complexity by reducing the integration over the whole domain to a summation of integration over smaller elements. It is shown in section <ref> that this choice of shape functions lead to a relatively sparse graph. It simplifies the computation by reducing the size of the graph optimization problem. The simplified form of Π specialized for the hat-functions is presented in Eq (<ref>).


    Π_N(𝐚)  
        =∑_i=1^N a_i-1^2 ( ∫_x_i-1^x_ip/2ϕ'^2_i-1 + q/2ϕ^2_i-1 dx ) 
    
        + a_i^2 ( ∫_x_i-1^x_ip/2ϕ'^2_i + q/2ϕ^2_i dx )
    
        + a_i-1a_i(∫_x_i-1^x_i pϕ'_i-1ϕ'_i + qϕ_i-1ϕ_i dx)
    
        - a_i-1(∫_x_i-1^x_i fϕ_i-1dx) -a_i(∫_x_i-1^x_i fϕ_idx)


This form of Π promotes modularity in computation
and allows expressing the functional as 

    Π_N = ∑_i=1^N𝐀_𝐢.𝐒_𝐢

where vectors 𝐀_𝐢≡𝐀_𝐢(a_i-1,a_i) and 𝐒_𝐢≡𝐒_𝐢(p,q,f) are defined for each element in (<ref>). The vector 𝐒_𝐢 is independent of state 𝐚 and is therefore only computed once in the whole procedure. 

    𝐀_𝐢 = [a^2_i-1  ,   a^2_i  ,   a_i-1a_i  ,   a_i-1  ,   a_i]^T 
    𝐒_𝐢 = [  ∫_x_i-1^x_ip/2ϕ'^2_i-1 + q/2ϕ^2_i-1 dx    ,  ∫_x_i-1^x_ip/2ϕ'^2_i + q/2ϕ^2_i dx    , 
    ∫_x_i-1^x_i pϕ'_i-1ϕ'_i + qϕ_i-1ϕ_i dx   ,   -∫_x_i-1^x_i fϕ_i-1dx   ,
       -∫_x_i-1^x_i fϕ_idx ]^T










    Example


Consider the differential equation  with boundary conditions u(0)=0 and u(1) = 1. 

    d^2 u /dx^2 = 0      0<x<1

Functional: 

    Π[u] = 1/2∫_0^1 u'^2 dx

For simplicity, consider a grid with a uniform mesh of 2 elements and 3 nodes: 


Using linear interpolants for the elements,

    u(x) = {
    a_0 (1 - 2x)+ a_1 (2x)     0<x≤ 0.5
    
     
    a_1 (2 - 2x)+ a_2 (2x-1)     0.5<x≤ 1
    .

The functional with the FE discretization: 

    Π_N(a) = (a_0-a_1)^2 + (a_1-a_2)^2

Modular representation of functional (Π_N = 𝐀_1.𝐒_1 + 𝐀_2.𝐒_2): 

    𝐀_1 = [a^2_0  ,   a^2_1  ,   a_0a_1  ,   a_0  ,   a_1]^T 
    𝐀_2 = [a^2_1  ,   a^2_2  ,   a_1a_2  ,   a_1  ,   a_2]^T 
    𝐒_1 = 𝐒_2  =  [1   ,   1   ,   -2   ,   0   ,   0 ]^T









§ GRAPH COLORING PROBLEM


Quantum annealing methods are tailored to find the lowest energy states in an Ising system defined in Eq (<ref>). The Ising hamiltonian defines a binary graph coloring problem with each vertex of graph or qubit labeled as +1 or -1. The value of the qubits determine the free variable, in this case, 𝐚. The parameters H_i and J_ij are defined such that the Ising hamiltonian, for a labeling representing the state, 𝐚, corresponds to the functional Π_N(𝐚). These problems, namely, the representation of state and estimation of parameters are addressed in this section.  



 §.§ Representation of State

Representation of a functional in terms of continuous variables is not feasible on quantum architectures. Due to this limitation, the values of each a_i (i^th component of 𝐚) are chosen from a finite set of values based on the labeling of qubits. The representation presented here permits 3 possible values of a_i at each node. In particular, for each node (indexed `i'), the state (a_i) is exactly determined by the labeling of qubits q^i_1, q^i_2 and q^i_3 with the i^th node taking values in the set { v_i_1, v_i_2, v_i_3}. Eq (<ref>) defines a mapping (q^i_1,q^i_2,q^i_3)→ a_i as tabulated in Table <ref>. It is observed that the mapping results in a_i ∈{v_i_1, v_i_2, v_i_3} when two qubits are labeled -1 and one qubit is labeled +1. Next it shown that the Ising parameters can be manipulated to make these labelings energetically favourable, thereby eliminating the occurrence of undesirable labels. 


    a_i = ∑_j=1^3 v_i_jq^i_j +1/2












    Example (Continued)


In general, the set {v_i_1, v_i_2, v_i_3} is different for each node. However, for simplicity, consider the same set of admissible states for all three nodes given by {v_i_1, v_i_2, v_i_3}≡{0,0.5,1}. Each node is defined by three qubits as follows:




We know that the three qubits each defining the solution at the first and last nodes should take up choices 1 and 3, respectively, due to boundary conditions. The choice for the second node is to be solved.





























 §.§ Parameter Estimation

To promote modularity, the graph representation is decomposed into two component subgraphs, namely, nodal graph and element graph. Each node and element of the FE discretization is endowed with a node graph and element graph, respectively. This allows to refine the mesh by extending the graph. 





  §.§.§ Nodal Graph

The nodal graph is given by a fully connected graph with three vertices representing the three qubits of the FE node. The nodal graph ensures that the energy minimizing states of the Ising hamiltonian corresponds to state 𝐚 with favorable choice of a_i ∈{ v_i_1, v_i_2, v_i_3} with equal probability. As mentioned earlier, the set of favorable labeling of qubits at a node is given by {Q_1, Q_2, Q_3}≡{ (1,-1,-1),(-1,1,-1),(-1,-1,1) }. Since each of the three labelling is equally likely in the absence of any functional minimization, it is expected that the same value of the coupling strength (Ĵ) for each connection and the field strength (H) for each node is used. A choice of Ĵ and H that fulfill these conditions is presented in Fig <ref>. Here, all the field and interaction terms for the nodal graph are given a value of one. The application of the Dirichlet boundary condition is also done by augmenting the field strength of the nodal graph. For example, by switching the field term H corresponding to the second qubit q_2 of a boundary node `b' to -1 forces a lower value of the functional for the boundary node state of ( -1, +1, -1), which corresponds to the solution v_b_2. This allows us to encode the value at the boundary to be v_b_2. 








  §.§.§ Element Graph


The element graph is used to make the energy of minimizing states of graph correspond to the value of functional Π_N of the continuous problem. Each element graph encodes the contribution of the respective element to the functional. Since the contribution of each element is dependent on the values at the nodes of the element, the element graph is constructed by connecting the vertices of neighbouring nodes. In particular, the site-site interaction in the n^th element graph can be estimated as a matrix, J^n, where (J^n)_kl represents the coupling energy between qubits q^i_k (k^th qubit of i^th node) and q^j_l (l^th qubit of j^th node) with i,j being the nodes of n^th element. As shown in the previous section, the contribution of the n^th element towards the functional, based on the choice of a compact basis function, is evaluated as 𝐀_𝐧.𝐒_𝐧. The elements of the vector, 𝐀_𝐧≡𝐀_𝐧(a_i,a_j) can therefore, take nine (3× 3) possible values based on the values of (a_i,a_j). For a particular choice of labeling of qubit the Ising energy of element graph is estimated as E = ∑_k=1^3 ∑_l=1^3 (J^n)_kl q^i_k q^j_l. When the labeling is chosen appropriately (each node has two `-1' and one `+1' label), this energy equals to the value of functional for corresponding state, 𝐚, as shown in Eq (<ref>). This relation can be used to estimate J^n by solving a set of nine independent linear equations presented. It is important observation that the independence of these set of equations relies on the fact that for any node, v_i_k≠ v_i_l for k≠ l. Additionally, the energy of the element graph breaks the symmetry between the states that minimize the energy of nodal graph, however, the values of J^n should be judiciously scaled (uniformly along all elements) such that energy of unfavourable states remain high.  


    ∑_k=1^3 ∑_l=1^3 (J^n)_kl q^i_k q^j_l = 𝐀_𝐧(a_i,a_j).𝐒_𝐧










    Example (Continued)


A single element with two nodes admits to the following connectivity: 

The estimated parameters reflect the contribution of element to the functional for a given choice of labeling: 

Sample 1: 


In the above Figure, both nodes take up choice 1 (a_i=a_j=0). The interaction energy for qubits: E = J_11 - J_12 - J_13 - J_21  +J_22 + J_23 - J_31 + J_32 + J_33 = (a_i-a_j)^2 = 0

Sample 2: 


In the above Figure, node i takes up choice 1 (a_i=0), while node j takes up choice 2 (a_j=0.5). The interaction energy for qubits: E = -J_11 + J_12 + J_13 + J_21  -J_22 - J_23 - J_31 + J_32 + J_33 = (a_i-a_j)^2 = 0.25 


Collectively solving the equation for all 9 such possibilities (as shown in Eq (<ref>)).

    [ +1 -1 -1 -1 +1 +1 -1 +1 +1; -1 +1 +1 +1 -1 -1 -1 +1 +1; -1 +1 +1 -1 +1 +1 +1 -1 -1; -1 +1 -1 +1 -1 +1 +1 -1 +1; +1 -1 +1 -1 +1 -1 +1 -1 +1; +1 -1 +1 +1 -1 +1 -1 +1 -1; -1 -1 +1 +1 +1 -1 +1 +1 -1; +1 +1 -1 -1 -1 +1 +1 +1 -1; +1 +1 -1 +1 +1 -1 -1 -1 +1 ][ J^n_11; J^n_12; J^n_13; J^n_21; J^n_22; J^n_23; J^n_31; J^n_32; J^n_33 ]
    
    = [ (v_i_1-v_j_1)^2; (v_i_2-v_j_1)^2; (v_i_3-v_j_1)^2; (v_i_1-v_j_2)^2; (v_i_2-v_j_2)^2; (v_i_3-v_j_2)^2; (v_i_1-v_j_3)^2; (v_i_2-v_j_3)^2; (v_i_3-v_j_3)^2 ]



    J^1 =  J^2 =   [ 0.1250 0.3750 0.3750; 0.3750 0.5000 0.3750; 0.3750 0.3750 0.1250 ]




The above parameters will exactly reproduce the functional in the interaction term. The boundary conditions are enforced, by setting self interaction term for qubits q^0_1, q^2_3 to H=-1. This locks the state at the 1^st boundary node as a_0 = v_0_1 = 0 and at the 2^nd boundary node as a_2 = v_2_3 = 1. Energy minimization of the resulting Ising hamiltonian gives a_1 = v_1_2= 0.5, which is the exact solution for the discretized problem. 






The process of the graphical representation of the discretized functional using the nodal and element graphs is referred to as “Assembly". Each node and element is endowed with a nodal and element graph, respectively. The effective site-site interaction energy is estimated by summing the nodal coupling strength, Ĵ, over all nodes, and element coupling strength, J, over all elements. Due to the nature of discretization, N element graphs and N+1 nodal graphs are required for representing an N-element discretization of the domain. The assembled graph, from here on, is referred to as the logical graph. Connectivity of logical graphs for one-element and four-element discretization is presented in Fig <ref>.

Two fundamental issues in this approach are addressed next using the box algorithm. Firstly, the choices at a node { v_i_1, v_i_2, v_i_3} were set in stone during initialization. The box algorithm makes this choice flexible. Secondly, as the number of nodes increase, three choices are insufficient. The number of qubits needed at a node must increase to make more choices available. Box algorithm, however, only requires three qubits per node for any level of discretization.







§ BOX ALGORITHM


In this section, an iterative procedure is developed to minimize the functional, Π_N, using the graph coloring representation discussed in the previous section. For a particular choice of {v_i_1,v_i_2,v_i_3}, defined as Eq (<ref>), the possible values of the state a_i at the i^th node are specialized to the set { u_i^c - r , u_i^c , u_i^c + r }, i.e., 


    v_i_j = u_i^c + r(j-2)


The quantities, 𝐮^𝐜= ( u_0^c, u_2^c,..., u_N^c ) and r  are hereafter referred to as box center and the slack variable, respectively. The intention is to approximate functions using the box center while the slack variable provides a bound on this approximation. The precise meaning of this bound is presented later in this section. A linear approximation of f(x) = √(x) using ten nodes is presented in Fig <ref> for different box centers and the slack variable (which can be interpreted as the box size). The function, f(x), is approximated as 𝐮^𝐜 at the nodes with linear interpolation in between the nodes. The blue region describes the possible value of the interpolation if the value at any node is perturbed within the range of ± r. In Fig <ref>(a), an exact approximation of the function at the nodes is presented with a slack variable of 0.2. In (b) the same approximation with a slack variable of 0.02 is presented. The same approximation is given in the two cases, but the bound on nodal values of (b) is tighter than (a). In part (c), the approximation is not exact, however, it lies within the bounds. In part (d), the approximation is neither exact nor within the bound. In the context of the vectorial representation of the coefficients, 𝐚, these bounds are represented as 3^N - 1 points on the surface of a box, defined as, ||𝐚- 𝐮^𝐜||_∞ = r. An illustration for the vectorial representation of a two nodes element is presented in Fig <ref>.  The solution is sampled from a 3×3 grid in the a_1-a_2 vector space.
 


In the discrete setting, the solution to the differential equation can be equivalently reduced to minimization of a function of the form:  𝐚^𝐓 𝐌 𝐚 where 𝐌 is some positive definite matrix. The vector 𝐚 takes value in one of the 3^N possibilities. The minimizer (not necessarily unique) is given by Eq (<ref>). The solution, 𝐚^min, need not coincide with the best approximation solution, 𝐚^b.a. of the continuous problem. In the illustration presented in Fig <ref>, the center is depicted as the solution (𝐚^min = 𝐮^c), the minimum is then contained within the elliptic region of the contour with 𝐚^min on the edge. Geometrically, this gives ||𝐚^min- 𝐚^b.a|| ≤ d ≤√(2)r(1+λ_max/λ_min) where λ_max and λ_min are the maximum and the minimum non–zero eigenvalues of M, respectively. This suggests that as the box size decreases, the corresponding 𝐮^c approaches to the best approximation solution of u. This argument is extended to n dimensions with the bound given in Eq (<ref>). 


    𝐚^min =   _a_i∈{u_i^c - r , u_i^c , u_i^c + r }𝐚^𝐓 𝐌 𝐚



    ||𝐚^min- 𝐚^b.a|| ≤ 2(1+(n-1)λ_max/λ_min)r/√(n)








 §.§ Iterative Procedure




In this section we present the details of the iterative procedure which locates the solution of the discretized problem, 𝐚^min, and updates the box center and slack variable such that 𝐮^C approaches the solution of the continuous problem (in the sense of best approximation).  

The necessary information required for defining the functional is stored in the vector 𝐒_𝐢. It is computed once at the beginning of the procedure as the problem definition stage. The procedure is initiated with a guessed solution of the vector, 𝐚, provided as a box centered at 𝐮^𝐜. The boundary nodes with Dirichlet boundary conditions are assigned the boundary value as the initial guess. The slack variable is initialized with an arbitrary scalar value. A better initial guess for r is the one which bounds the solution in the box defined by 𝐮^𝐜. Such initial guesses require fewer iterations in comparison to arbitrary ones; however, starting with a good guess is not a necessary condition for convergence. The Ising parameters, H, Ĵ and J are estimated as discussed in section <ref>.  


In this work, D-Wave's 2000Q processor is used. This processor has a Chimera-type structure with 2048 qubits and 6016 couplers <cit.>. A direct solution of the optimization problem by re-numeration of qubits is not possible as the assembled graph cannot be found in any subgraph of the physical graph, i.e., the processor. Therefore, it is required that the logical graph is mapped onto the physical graph via the process of embedding. This problem in itself is NP-hard and is not discussed here for brevity. The reader is referred to <cit.> for a discussion on this topic. The topology of the logical graph remains unchanged over the iterations. The search for embedding of a map is only conducted once, and in subsequent iterations, the self-interaction and the site-site interaction parameters are updated for the same embedding. 

The use of three qubits per node in this paper allows the D-wave system to search for a minimum over a space of 3^N solution vectors in a single run. In each iteration, the box center is translated to the energy minimizer, 𝐚^min. This move is referred to as the translation step. In the case where the minimizing state is found at the center, the box size is reduced, and the search is continued with a smaller bound on error. This move is referred to as the contraction step. The complete procedure is presented in Algorithm <ref>.




    ∑_<k,l>J^i_kl q^i-1_k q^i_l = 𝐀_𝐢.𝐒_𝐢

In each iteration, the best solution among 3^N possibilities is obtained using quantum annealing. 
  map all nodes from the graph model onto the physical graph. 
This mapping should preserve the minimum energy states.
The energy minimization G(V,E)  G'(V',E'). T(v_i) ⊂ V' for all v_i ∈ V such that (1) T(v_i)∩ T(v_j) = ϕ for all v_i ≠ v_j, (2) T(v_i) is connected in G', and (3) If (v_i,v_j) ∈ E then there exist an edge between subsets T(v_i) and T(v_j) in E'


𝐯_𝐜 = { v^1_c,...,v^i_c,...,v^n_c }

Moved from conclusions:  These solution vectors are at corners, face and edge centers and center of a box in the N-dimensional solution space. The iteration process shown here either moves along the direction of the minimum energy vertex (excluding center) in the box or shrinks the box towards the center.













    Example (Continued)


In the box algorithm, the set {v_i_1,v_i_2,v_i_3} is constructed using the box center and the slack variable.

With the application of the boundary condition, the favourable labeling of qubits give following three choices in the solution (𝐈, 𝐈𝐈, 𝐈𝐈𝐈). 

The values of  u^c_1 and r are initialized arbitrarily. One of the solutions among 𝐈, 𝐈𝐈 or 𝐈𝐈𝐈 is selected by the annealer. If the minimizer is found to be solution  𝐈𝐈 then the algorithm proceeds with the contraction step by halving the value of r. 
If solutions 𝐈 or 𝐈𝐈𝐈 are selected, then the algorithm proceeds with the translation step by setting the new box center to u_1^c + r or u_1^c - r, respectively.








§ RESULTS

The deformation of a bar under axial loading is modelled using an equation of form (<ref>). In particular, the deflection (u) of a bar is related to the elastic stiffness, (E), cross-sectional area (A), and the applied body force, (f) using Eq (<ref>). The functional (Eq <ref>) is referred to as the potential energy of the system. The corresponding discretized form of the potential energy for piece-wise linear E,A and f is calculated using Eq (<ref>) where E_i,A_i and f_i represents the elastic stiffness, area and applied body force, respectively, at the center of the i^th element.

    (EAu')' + f = 0


    Π_N[u] = ∑_i=1^NN/2 E_iA_i (a_i-a_i-1)^2 - 1/2Nf_i (a_i+a_i-1)



Two test cases are presented in Fig <ref>. In test case (a), a bar with a discontinuity in EA is simulated. The body force is not applied in this case. A four-element discretization is used. The initial guesses are taken as 𝐮^𝐜 = {0,0.25, 0.5, 0.75, 1 } and r=0.2. The numerical solution is observed to approach the exact solution in this case. The convergence in the functional is also evident. 

In the test case (b), a bar with continuously varying EA is simulated. A linearly varying body-force is supplied. A six-element discretization of the bar is used with 𝐮^𝐜 = {0,1/6, 2/6,3/6,4/6,5/6,1 } and r=0.2. Based on the theory of finite element methods, the exact minimization of the energy in discretized space leads to a stiffer solution in comparison to the exact solution. It is observed that the numerical solution approaches the exact solution at nodes which is characteristic of finite element methods. Energy is also observed to be converging towards the finite element solution u^fem in this case. The mismatch of u within the element is expected to decrease with refinement in discretization.

Some implementation details on the D-wave architecture are relevant here. Although the mapping only requires three qubits per node, embedding of this graph into the Chimera graph produced an overhead of 9 qubits per node - constant over a range of discretizations. It is understandable since a complete graph of three qubits used to represent a node is not directly represented on the Chimera graph. In the future, the use of two qubits to represent a node can also be explored. While this ensures we still sample a large enough (2^N vectors) solution space in a single run, we lose out information on the box center energy that is important for reducing the range of the slack variable. However, it is possible to compute the solution at the box center classically. The overhead of performing classical solutions can be offset by the fact that a two-qubit representation has smaller complete sub-graphs and is easier to embed in the physical graph. Another important task in quantum computing is error suppression. Quantum processors, unlike classical computers, do not have parity correction algorithms due to the no-cloning theorem. A compilation of popular methods for quantum error correction is presented in <cit.>. Energy re-scaling is one of the simpler approaches and is employed in this work. Here, in the estimation of J^n, the energy was rescaled to ensure that the energy gap between feasible and unfeasible states is increased while maintaining a similar energy landscape. This step is a heuristic remedy for minimizing noise in quantum computation, and has no bearing in the theoretical convergence of the algorithm.



§ CONCLUSIONS AND FUTURE WORK

Recent rapid developments in quantum annealers warrant further investigation into re-formulation of scientific computation problems as graph coloring problems. The use of quantum computing for solving differential equations has, to date, focused on the use of a gate–computing based QLSA algorithm. This algorithm attempts to sample from the solution space of the linear system of equations Ax=b. In the quantum annealer based algorithm described here, we do not solve the system of equations. We instead map the discretized version of the energy function of the differential equation to an Ising hamiltonian. The solution vector, x, is directly obtained as the ground state of the qubits. The algorithm has low memory requirements since the global matrix is not stored and the local matrices are encoded in the interaction weights of the Ising model. Further, the box algorithm allows mapping of up/down spin states of qubits in the ground state to rational numbers involved in the solution vector. 

Since we primarily solve the Ising model, the cost of computation is tied to the performance of the quantum annealer (<cit.>). Within each iteration however, equation Eq (<ref>) is solved for each element, leading to at least O(n) operations. 

We have shown that the box algorithm indeed guarantees convergence to the best approximation of the solution in the discretized space as the box size goes to zero. However, some improvements could be made to reduce the number of minimization runs. We could use the statistics of solutions that D-wave system returns from a single minimization run to drive the iteration process in an arbitrary direction. This data can also be used to heuristically develop `local' potential energy maps that can be used to identify larger step sizes for faster convergence. With future scaling up of quantum annealers up to millions of qubits, it will be possible to solve challenging engineering solid and fluid mechanics problems using quantum annealers.



§ ACKNOWLEDGMENTS

The authors would like to acknowledge Universities Space Research Association (USRA) for providing access to the use of the D-Wave quantum computer in the USRA-NASA-Google quantum Artificial Intelligence Laboratory at NASA's Ames Research Center. 


[breakable,enhanced] 


    Example (Cont.)




Representation of a functional in terms of continuous variables is not feasible on quantum architectures. Thus it is imperative to further approximate this functional using discrete variables. In this section, a graph coloring based functional representation is formulated. An Ising hamiltonian is used to encode the functional on discrete points in V_N. Quantum annealing methods are tailored to find the lowest energy states in an Ising system defined in Eq (<ref>). The parameters of self-interaction term (H_i) and the site-site interaction term (J_ij) are chosen in a way that the lower energy states correspond to feasible solutions for u({ a_i }). Additionally, it is required that the Ising hamiltonian of the lower energy states should correspond to Π[{ a_i }].
The value of a_i for each node is represented using 3 logical qubits, namely, { q_i^j}. The variable, a_i takes value from a finite set { v_i^j} with cardinality 3. The evaluation of the test function at the i^th node is estimated using Eq (<ref>) where q_i^j takes value in {+1,-1}. 
Different graphs are required to study the different levels of discretization. Two sub-graphs are introduced to maintain generality, namely, Element graph and Nodal graph. All logical graphs can be generated using a combination of these sub-graphs. The contribution of each element to the functional is encoded in the element graph. The hamiltonian is augmented to enforce the preference of feasible states. This additional energy is provided through the Nodal graph. 


 §.§ Nodal Graph

The three qubits assigned to individual nodes are fully connected among themselves. For the purpose of elucidation, the state of the i^th node can be denoted as a triplet { q^1_i,q^2_i,q^3_i }. The states defined in the set Q = {{ 1,-1,-1 },{ -1,1,-1 }, { -1,-1,1 }}≡{{ Q_1 },{ Q_2 }, { Q_3 }} are chosen as allowable states. It is left for the reader to observe that these states result in a_i = {{ v_i^1},{ v_i^2}, { v_i^3}}, respectively. The strength of parameters are decided in a way that the energy landscape of the initial problem is preserved. This requirement is imposed by using the same value of the coupling strength (Ĵ) for each connection and the field strength (H) for each node. A choice of Ĵ and H that fulfill these conditions is presented in Fig <ref>.  

The connection between the qubits of consecutive nodes is determined using the element graph. This graph is constructed in a manner that for a feasible labeling of nodes (i.e. when state of all nodes are in the set Q), the evaluation of the Ising hamiltonian equals the evaluation of functional for corresponding a_i. The site-site interaction for the i^th element can be estimated as a matrix, J^i, where (J^i)_kl represents the coupling energy between qubits q_i-1^k and q_i^l, 

The contribution of the i^th element towards the functional, based on the choice of a compact basis function, is evaluated as 𝐀_𝐢.𝐒_𝐢. The elements of the vector, 𝐀_𝐢 are polynomials in a_i-1 and a_i which are uniquely defined for a feasible state of nodes. This allows determination the matrix J^i  by solving the system of nine linear equations given in Eq (<ref>). The notation 𝐀_𝐢(Q_k,Q_l) refers to the evaluation of the vector 𝐀_𝐢 when the k^th state is selected for the (i-1)^th node and the l^th state is selected for the i^th node. 


    [ ⟵ Q_1 ⟶; ⟵ Q_2 ⟶; ⟵ Q_3 ⟶ ][ J^i_11 J^i_12 J^i_13; J^i_21 J^i_22 J^i_23; J^i_31 J^i_32 J^i_33 ][   ↑   ↑   ↑; Q_1 Q_2 Q_3;   ↓   ↓   ↓ ]
    
    =
    [ 𝐀_𝐢(Q_1,Q_1).𝐒_𝐢 𝐀_𝐢(Q_1,Q_2).𝐒_𝐢 𝐀_𝐢(Q_1,Q_3).𝐒_𝐢; 𝐀_𝐢(Q_2,Q_1).𝐒_𝐢 𝐀_𝐢(Q_2,Q_2).𝐒_𝐢 𝐀_𝐢(Q_2,Q_3).𝐒_𝐢; 𝐀_𝐢(Q_3,Q_1).𝐒_𝐢 𝐀_𝐢(Q_3,Q_2).𝐒_𝐢 𝐀_𝐢(Q_3,Q_3).𝐒_𝐢 ]




    [  1 -1 -1 -1 +1 +1 -1 +1 +1; -1 +1 +1 +1 -1 -1 -1 +1 +1; -1 +1 +1 -1 +1 +1 +1 -1 -1; -1 +1 -1 +1 -1 +1 +1 -1 +1;  1 -1 +1 -1 +1 -1 +1 -1 +1;  1 -1 +1 +1 -1 +1 -1 +1 -1; -1 -1 +1 +1 +1 -1 +1 +1 -1;  1 +1 -1 -1 -1 +1 +1 +1 -1;  1 +1 -1 +1 +1 -1 -1 -1 +1 ][ J^i_11; J^i_12; J^i_13; J^i_21; J^i_22; J^i_23; J^i_31; J^i_32; J^i_33 ]
    
    = [ 𝐀_𝐢 (v_(i-1)_1,v_i_1).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_2,v_i_1).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_3,v_i_1).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_1,v_i_2).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_2,v_i_2).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_3,v_i_2).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_1,v_i_3).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_2,v_i_3).𝐒_𝐢; 𝐀_𝐢 (v_(i-1)_3,v_i_3).𝐒_𝐢 ]



apsrev


