An integrated programming and development environment for adiabatic quantum optimization

The physical basis for the AQC model was first established in terms of quantum annealing by Kadowaki and Nishimori [38]. Farhi et al as well as others later formalized these ideas as a means of universal computation [2, 7, 39]. Several efforts have since shown the equivalence between the AQC model and other quantum computing models [40, 41]. In a generalized AQC algorithm [2, 3, 7], a quantum physical system of n qubits is evolved under the Schrödinger equation

$$\begin{eqnarray}&&i\hbar \frac{\partial \left| \psi \left( t \right) \right\rangle }{\partial t}=H\left( t \right)\left| \psi \left( t \right) \right\rangle \end{eqnarray} \tag{ 1 }$$

according to a time-dependent Hamiltonian

$$\begin{eqnarray}&&H\left( t \right)=A\left( t \right){{H}_{I}}+B\left( t \right){{H}_{P}}\end{eqnarray} \tag{ 2 }$$

that interpolates between the initial Hamiltonian ${{H}_{I}}$ and the final (problem) Hamiltonian ${{H}_{P}}$ from an initial time $t={{t}_{0}}$ to a final time t = T. We shall assume ${{t}_{0}}=0$. In equation (2), the schedules A(t) and B(t) satisfy the boundary conditions $A\left( 0 \right)\gg B\left( 0 \right)$ and $A\left( T \right)\ll B\left( T \right)$, while the quantum system is initially prepared in the lowest-energy eigenstate of ${{H}_{I}}$. Given the instantaneous eigenvalue equations

$$\begin{eqnarray}&&H\left( t \right)\left| {{\tilde{\phi }}_{j}}\left( t \right) \right\rangle ={{E}_{j}}\left( t \right)\left| {{\tilde{\phi }}_{j}}\left( t \right) \right\rangle ,\end{eqnarray} \tag{ 3 }$$

with $j=0,1,\ldots {{2}^{n}}-1$ labeling states of monotonically increasing energy, the initial state condition implies $\left| \psi \left( 0 \right) \right\rangle =\left| {{\tilde{\phi }}_{0}}\left( 0 \right) \right\rangle$.

We define the energy gap between the ground and first excited state as

$$\begin{eqnarray}&&\Delta \left( t \right)={{E}_{1}}\left( t \right)-{{E}_{0}}\left( t \right),\end{eqnarray} \tag{ 4 }$$

in which we neglect possible ground state degeneracy for simplicity. If the energy gap is always strictly greater than zero, i.e., $\forall t:\Delta \left( t \right)>0$, then the state $\left| \psi \left( t \right) \right\rangle$ will remain in the instantaneous ground state with high probability provided certain bounds on the rate of change of the Hamiltonian are satisfied [2]. Consequently, evolution under equation (1) to the time T prepares the final state $\left| \psi \left( T \right) \right\rangle$ in the lowest energy eigenstate of ${{H}_{P}}$. By making a judicious choice of the final Hamiltonian ${{H}_{P}}$, the prepared final state may encode the solution to a computation. In order to ensure the computation is correct, the adiabatic condition must be satisfied. In its simplest interpretation, this implies the global time T must be much larger than the inverse of the minimum spectral gap of H(t) [2]. More sophisticated analysis shows that better results may be obtained by adjusting the evolution schedule according to the local energy gap [30]. In either case, failure to ensure the adiabatic condition risks the possibility that the final state will not belong to the ground state manifold of ${{H}_{P}}$ but rather to an excited state and an error in the computation. It is notable that the spectral gap depends not only on the problem to be solved, but also on how the problem is implemented as a quantum program. Understanding input influence on a program run time and error rates is an open question in quantum computer science.

In specializing to the AQO algorithm, we require a quantum logical system of n qubits with an initial Hamiltonian

$$\begin{eqnarray}&&{{H}_{I}}=-\mathop{\sum }\limits_{i\in {{V}_{P}}}{{X}_{i}}\end{eqnarray} \tag{ 5 }$$

and final Ising Hamiltonian

$$\begin{eqnarray}&&{{H}_{P}}=-\mathop{\sum }\limits_{i\in {{V}_{P}}}{{\alpha }_{i}}{{Z}_{i}}-\mathop{\sum }\limits_{\left( i,j \right)\in {{E}_{P}}}{{\beta }_{i,j}}{{Z}_{i}}{{Z}_{j}},\end{eqnarray} \tag{ 6 }$$

where ${{X}_{i}}$ and ${{Z}_{i}}$ are the Pauli operators for the ith qubit [1], ${{\alpha }_{i}}$ is the bias on the ith qubit, and ${{\beta }_{i,j}}$ is the coupling between qubits i and j. As shown in the section below, the graph ${{G}_{P}}=\left( {{V}_{P}},{{E}_{P}} \right)$ with vertex set ${{V}_{P}}$ $\left( \left| {{V}_{P}} \right|\equiv n \right)$ and edge set ${{E}_{P}}$ defines an input optimization problem, cf the weight matrix ${\bf P}$ in equation (7). The final Hamiltonian is diagonal in the basis defined by the tensor products of the $\pm 1$ eigenstates of the ${{Z}_{i}}$ operators. This basis will also serve as the computational basis. For comparison, the ground state of the initial Hamiltonian (and initial state of the AQO algorithm) is the symmetric superposition of these computational basis states and has an eigenvalue $-n$.

An important consequence arising from the choice of the initial and final Hamiltonians, respectively, equations (5) and (6), is that the time-dependent Hamiltonian H(t) of equation (2) is not capable of universal adiabatic quantum computation. Extending the form of the Hamiltonian beyond the Ising model, for example, to the 2-local ZZXX Hamiltonian of Biamonte and Love [42], would support universal computation, but we do not consider that case here.

Any binary optimization problem (BOP) can be mapped into the form of the final Hamiltonian in equation (6). In doing so, we define the classical input to the AQO algorithm as a QUBO problem. This is because non-binary as well as constrained optimization problems can be reduced to QUBO [43], with multiple methods for performing the reduction available [44]. The QUBO problem is to find

$$\begin{eqnarray}&&{\rm arg} \mathop{{\rm min} }\limits_{{\bf x}\in {{{\bf B}}^{m}}}{{{\bf x}}^{T}}{\bf Px},\end{eqnarray} \tag{ 7 }$$

where ${\bf x}$ is a vector of n binary variables with ${{x}_{i}}\in \left\{ 0,1 \right\}$ and ${\bf P}$ is an n-by-n symmetric real-valued cost matrix. We use the weight matrix ${\bf P}$ to define the weighted adjacency matrix of the input (problem) graph ${{G}_{P}}$ introduced in equation (6). The graph ${{G}_{P}}$ has a vertex set ${{V}_{P}}$ of size $\left| {{V}_{P}} \right|\equiv n$ and an edge set ${{E}_{P}}$ defined as $\left( i,j \right)\in {{E}_{P}}$ iff ${{P}_{i,j}}\ne 0$. From this point of view, programming the AQO algorithm requires mapping the matrix ${\bf P}$ to the biases and couplings of the Ising Hamiltonian. It has been shown previously by Choi that parameterization of the logical Ising Hamiltonian in equation (6) may be given in terms of the QUBO problem as [24]

$$\begin{eqnarray}&&{{\alpha }_{i}}=\frac{1}{2}{{P}_{i,i}}+\frac{1}{4}\mathop{\sum }\limits_{j=1}^{n}{{P}_{i,j}}\quad {\rm for}\;\;i=1\;\;{\rm to}\;\;n,\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{{\beta }_{i,j}}=\frac{1}{4}{{P}_{i,j}}\quad {\rm for}\;\;i<j=1\;\;{\rm to}\;\;n.\end{eqnarray} \tag{ 9 }$$

We may also add an energy shift to the Ising Hamiltonian in equation (6) of the form

$$\begin{eqnarray}&&\gamma =\frac{1}{4}\mathop{\sum }\limits_{i,j=1}^{n}{{P}_{i,j}}+\frac{1}{2}\mathop{\sum }\limits_{i=1}^{n}{{P}_{i,i}}\end{eqnarray} \tag{ 10 }$$

in order to match the energies of the solution state. Although this shift does not affect the solution obtained using AQC, it must be accounted for in reporting the minimal value in equation (7).

Whether or not the logical Hamiltonian in equation (6) is supported directly on a given hardware depends on the available connectivity of that hardware. We express the connectivity of a targeted processor in terms of its hardware graph ${{G}_{H}}=\left( {{V}_{H}},{{E}_{H}} \right)$. When any vertex can be coupled to any other vertex and $\left| {{V}_{H}} \right|\geqslant \left| {{V}_{P}} \right|$, then it is possible to support all possible input problems using a one-to-one mapping between the logical and physical qubits and the biases and couplings of the physical Hamiltonian. However, when ${{G}_{H}}$ is less than fully connected, then there are certain input problems that will not map directly into hardware. In such circumstances, it may be possible to embed the problem graph ${{G}_{P}}$ into the hardware graph ${{G}_{H}}$ via graph minor embedding [25, 45].

We formally define the minor embedding of a graph ${{G}_{P}}$ into a graph ${{G}_{H}}$ as a mapping $\phi :{{V}_{P}}\to {{V}_{H}}$ such that:

• (i)  
  each vertex i in ${{V}_{P}}$ is mapped to the vertex set of a connected subgraph ${{T}_{i}}$ of ${{G}_{H}}$.
• (ii)  
  if $\left( i,j \right)\in {{E}_{P}}$, then there exist ${{\tau }_{i}},{{\tau }_{j}}\in {{V}_{H}}$ such that ${{\tau }_{i}}\in {{T}_{i}}$, ${{\tau }_{j}}\in {{T}_{j}}$, and $\left( {{\tau }_{i}},{{\tau }_{j}} \right)\in {{E}_{H}}$.

If such a mapping ϕ exists, then ${{G}_{P}}$ is minor-embeddable in ${{G}_{H}}$, or ${{G}_{P}}$ is a minor of ${{G}_{H}}$. In subsequent discussions, we simply use the term embedding as a reference to minor embedding.

In adiabatic quantum programming, the vertices of the input graph ${{G}_{P}}$ represent the bits of a candidate solution to the QUBO problem, while the edges represent the presence of non-zero coupling coefficients, as defined in equations (8) and (9), respectively. The vertices of the hardware graph ${{G}_{H}}$ represent the physical qubits and the edges represent the couplings between qubits that are available in the hardware. An embedding maps each vertex in ${{V}_{P}}$ to a subset of ${{V}_{H}}$ and each edge in ${{E}_{P}}$ to edges between these subsets. When an embedding exists, then the resulting subgraph ${{G}^{*}}=\left( {{V}^{*}},{{E}^{*}} \right)$ of the hardware graph defines the physical Ising model

$$\begin{eqnarray}&&{{H}_{{{G}^{*}}}}=-\mathop{\sum }\limits_{k\in {{V}^{*}}}\alpha _{k}^{*}{{Z}_{k}}-\mathop{\sum }\limits_{\left( k,\ell \right)\in {{E}^{*}}}\beta _{k,\ell }^{*}{{Z}_{k}}{{Z}_{\ell }}.\end{eqnarray} \tag{ 11 }$$

The bias and coupling coefficients $\alpha _{k}^{*}$ and $\beta _{k,\ell }^{*}$ depend on the selected embedding ϕ per the requirements (i) and (ii) listed above. The physical Ising coefficients are defined as [45]

$$\begin{eqnarray}&&\alpha _{k}^{*}={{\alpha }_{i}}/\left| {{T}_{i}} \right|\quad {\rm for}\;{\rm each}\;\;k\in {{V}_{{{T}_{i}}}}\end{eqnarray} \tag{ 12 }$$

and for $k\ne \ell$

$$\begin{eqnarray}\beta _{k,\ell }^{*}=\left\{ \begin{array}{ccccccccccccccc} {{\beta }_{i,j}}/edges\left( {{T}_{i}},{{T}_{j}} \right) & {\rm for}\;\;k\in {{T}_{i}}\;\;{\rm and}\;\;\ell \in {{T}_{j}}\;\;{\rm and}\;\;i\ne j \\ J & {\rm for}\;\;k\in {{T}_{i}}\;\;{\rm and}\;\;\ell \in {{T}_{j}}\;\;{\rm and}\;\;i=j \\ 0 & {\rm otherwise} \\ \end{array} \right.,\end{eqnarray} \tag{ 13 }$$

where $edges\left( {{T}_{i}},{{T}_{j}} \right)$ is the number of edges between subgraphs ${{T}_{i}}$ and ${{T}_{j}}$. The constant J is chosen sufficiently large to force the qubits in each subgraph to be strongly correlated, with an upper bound on its value given previously by Choi [45]. Setting the embedded Ising model coefficients requires knowledge of the matrix ${\bf P}$ and the selected embedding implied by ${{G}^{*}}$ [25, 45]. The embedding need not be unique and, consequently, different instances of the Hamiltonian in equation (11) may correspond to the same logical problem of equation (7).

A key dependency in finding an embedding is the target hardware graph ${{G}_{H}}$. The hardware graph defines the vertices and connectivity that are available to express the Ising model. An example hardware graph is shown in figure 2. Finding those graphs that can be embedded into a fixed hardware graph is an example of subgraph isomorphism, which is known to be NP-Complete in difficulty [46]. For small hardware graphs, it is tractable to calculate the maximal minors of the graph, i.e., the minors of ${{G}_{H}}$ whose subgraphs represent all other graphs contained in ${{G}_{H}}$ [25]. However, this is a brute force approach and therefore does not scale favorably with hardware size. Similarly, attempts to find complete graphs as minors of an arbitrary hardware graph is known to be NP-Complete [47]. Alternatives to these brute force approaches include heuristic algorithms that incorporate knowledge of ${{G}_{H}}$ or that limit the types of input problems [25].

Zoom In Zoom Out Reset image size

At this point, we emphasize that the role of minor embedding is not as simple as identifying a physical Ising model that is equivalent to the logical Hamiltonian. Indeed, the embedding of a problem into a processor is not unique and, moreover, it is not well understood how different embeddings influence program behavior. There are known tradeoffs in the amount of time spent finding an embedding relative to the size of the embedded problem [25], but it remains unclear how to account for those costs in the benchmarking program performance.

In addition, the current approach to hardware embedding taken by JADE follows the decomposition of a BOP into a QUBO form using quadratization, i.e., decomposing into quadratic form [44]. However, an alternative programming sequence is to map a BOP directly into a multi-linear Ising model that is then decomposed into bilinear form [48]. The latter approach has led to the development of generalized gadgets [49] and, more recently, to resource efficient gadgets that replace multi-linear terms in the Hamiltonian with bilinear ones [50–52]. Gadget decompositions introduce additional ancilla qubits in much the same way that quadratization introduces ancilla bits. Biamonte has presented decompositions minimal in the number of gadget ancilla that would be especially relevant for comparing performance [50]. We have not explored the use of gadgets in the JADE programming model, or compared the overhead using quadratization, but we believe that the impact of this alternative programming model should be investigated.

We restrict our discussion to AQO algorithms that use a time-dependent Hamiltonian fitting the form of equation (2), which interpolates between an initial Hamiltonian ${{H}_{I}}$ and the problem Hamiltonian ${{H}_{P}}$ according to the time-dependent annealing schedules A(t) and B(t). More generally, individual biases and couplings can be time-dependent, e.g., ${{\alpha }_{i}}={{\alpha }_{i}}\left( t \right)$. In either case, the time-dependent schedules specify the rate at which the total Hamiltonian H(t) changes and, consequently, they play an important role in the computational error rates. In particular, the final time T needs to be sufficiently large to ensure the validity of the adiabatic condition, namely,

$$\begin{eqnarray}&&T\gg \frac{\mathcal{E}}{{{\Delta }^{*}}},\end{eqnarray} \tag{ 14 }$$

where ${{\Delta }^{*}}={{{\rm min} }_{t}}\;\Delta \left( t \right)$ is the global minimum of the spectral gap defined in equation (4) and $\mathcal{E}={{{\rm max} }_{t}}\;\left\langle {\rm d}H\left( t \right)/{\rm d}t \right\rangle$ is the maximal rate of change during evolution [2]. In the absence of information about $\Delta \left( t \right)$, it is difficult to ensure the adiabatic condition is satisfied. This uncertainty is one source of the difficulty in benchmarking adiabatic quantum programs. Recent results on amplifying spectral gaps [53] and developing fault-tolerant programs [54] suggest new methods for mitigating this uncertainty.

Although the annealing schedules are sufficient for coarsely specifying program execution, it is ultimately necessary to provide the physical implementation of those schedules in terms of hardware controls. The hardware controls that are available for tuning the biases and coupling of a processor must be capable of expressing programmed schedules. However, available controls are highly dependent on the physics underlying a processor and ensuring the exact implementation of an arbitrary annealing schedule may not be possible. Limitations on annealing schedules arising from constraints and dependencies of control values creates additional uncertainty in the benchmarking effort. Accounting for control constraints and quantification noise is necessary to provide a clear picture of how processor differences impact program behavior. For example, in the case of the family of processors from D-Wave Systems, biases and couplings can be mapped directly to models for the underlying superconductor Josephson-junction. However, the precision of this mapping is limited by the resolution of the on board digital to analog converters [8, 10].

In addition to the constraints expected from hardware design, it is also necessary to anticipate the influence of noise on program behavior. Two types of noise affecting quantum dynamics are classical noise in the controls and quantum noise in the system dynamics. Quantum noise may be modeled as an undesired interaction between computational qubits and non-control elements of the hardware. A specific example is the case of thermal influences on the quantum dynamics, which invalidate the pure state description in section 2 and undermine the adiabatic conditions [55]. Similarly, classical noise in the hardware controls yields a mixed-state description of the quantum dynamics and may bias program execution away from the solution of interest.

Once the time-dependent behavior of the Hamiltonian H(t) has been fully specified, it remains to execute the program. As noted before, the typical sequence begins by initializing the quantum computational register in the ground state of the initial Hamiltonian ${{H}_{I}}$. How initialization is implemented varies with processor and, more important, it may not be implemented perfectly. This additional source of noise must also be accounted for in evaluating program behavior as it is likely to influence the computational result. The remaining step in execution is to carry out the hardware control schedule and, therefore, the programmed computation.

After evolving to the final time T, the state of the computational register is determined using a suitable measurement or readout method. For the case of the AQO algorithm, the ground states at time T are computational eigenstates and, therefore, readout implies a direct measurement in the computational (Z) basis. As with program execution, it is more realistic to describe the readout process in terms of the hardware controls. This description includes capturing any noise or uncertainty in the measurement process.

The bit string generated from computational readout is the result of the quantum annealing process. However, mapping this result back to a solution for the original QUBO problem requires decoding measurements according to the inverse of the embedding map. For those cases where a tree of physical qubits represents a single logical qubit, it is necessary to check the value of all such qubits. In cases where measurement results within a tree disagree, then various strategies can resolve the uncertainty. One simple example is to use a majority vote. After decoding the computational readout, a solution to the original QUBO problem is produced and the program is complete. It may be necessary to repeat the execution of the program, for example, to gather statistics on the readout or solution states, however, the steps performed are similar to those described above.

JADE is designed to provide infrastructure for developing AQO programs and a computational engine for simulating them. This includes functionality for parsing input optimization problems, configuring new quantum hardware, and performing program profiling. Given this broad scope in functionality, JADE was designed for two distinct actors: the Analyst and the Engineer.

An Analyst represents a JADE user whose primary goal is to solve a discrete optimization problem. The Analyst requires a development environment that automates programming choices and execution sequences. In contrast, an Engineer expects to perform additional programming tasks such as customizing low-level Hamiltonian parameters, constructing specialized processor configurations, and defining embedding maps or annealing schedules. As seen in figure 3, this desired JADE functionality is encapsulated by the following use case model.

• Create a Problem—the Analyst constructs a discrete optimization problem as either a BOP or QUBO problem. In the case of the former, JADE converts the BOP to its corresponding QUBO representation. This use case creates a Problem entity.
• Solve a Problem—the Analyst selects a previously created Problem to solve using the AQO algorithm. This use case returns a Solution entity, which is the computed solution to the input problem.
• Create a Processor—the Engineer creates a processor configuration by specifying the number and connectivity of physical qubits. The Engineer may also customize the processor by specifying classical and quantum noise models as well as hardware control constraints. This use case creates a Processor entity.
• Create a Program—the Engineer creates a quantum program that is either a logical program or a physical program. A logical program is synthesized from selected Problem, Processor, and Embedding entities, while a physical program is synthesized only from a Processor. For the physical program, the Engineer sets the parameters of the final Ising Hamiltonian including biases, coupling, and annealing schedules. Both instances of this use case create a Program entity.
• Execute a Program—the Engineer executes a Program. With JADE, the Engineer submits the Program for simulation along with any profiling and simulations options. This use case creates a Result entity that corresponds to the computational readout following program execution. Note that the Result of a Program does not correspond necessarily with the Solution to a Problem, as the Result may require additional processing to generate a Solution.

Zoom In Zoom Out Reset image size

Alongside the use case model, we also present the system context model in figure 4. The system context describes the communication between JADE and its environment as driven by the use case model. The system context details how the Analyst and Engineer interact with the various input–output (I/O) data. As shown in figure 4, the six types of I/O data are: Problem, Processor, Embedding, Program, Result and Solution. These I/O entities are further specified in section 3.3.

Zoom In Zoom Out Reset image size

An Analyst only has access to Problem and Solution entities. However, we anticipate that JADE must synthesize other entities internally, for example, a Program is required to generate a Solution. Consequently, JADE will need private non-interactive methods for internal synthesis of the remaining entities. Although Processor, Embedding, and Program are generated by the system during the Analyst workflow, we do not explicitly model that dependency in figure 4.

JADE comprises three distinct components: JadeD, Sapphire, and NiCE. The JadeD component is responsible for data creation, management, synthesis, and verification, i.e., domain logic. The Sapphire component is responsible for the simulation of quantum programs according to user-defined plug-ins. The NiCE component, a pre-existing open source project, is used to integrate the JadeD and Sapphire components and to manage the computational work flow [59]. Each component provides an independent application programming interface (API).

Figure 5 highlights the interactions between the three components and the associated interfaces. Both JadeD and Sapphire couple to NiCE, which provides a user-driven coupling between program development and program execution. There is a dependency between JadeD and Sapphire due the latterʼs need to parse Program data structures. This dependency is restricted to a very narrow subset of the JadeD functionality and we expect future versions will isolate it in a separate shared library.

Zoom In Zoom Out Reset image size

The JadeD component handles creation and manipulation of quantum programming by exposing a basic create, retrieve, update, and delete interface. This interface enables generation, manipulation, and persistence of Entity data objects, which represent high-level abstractions of the various types of I/O data. The functional scope of JadeD includes parsing user-provided input into verified formats, validating that input, and generating subclasses of Entity tailored to specific input types. We define an IJadeD interface to specify how the JadeD component interacts with clients. By defining a formal interface, we are able to offer the option of supporting multiple JadeD variants.

As shown in figure 6, the IJadeD interface includes a number of methods for creating and storing entity instances. The JadeD class is a realization of this interface that provides a concrete implementation of the defined functionality. The JadeD implementation presented here uses a variety of object-oriented design patterns with the factory design pattern being the most significant [57]. The factory pattern is used to create and modify entities in an abstract manner, which pushes the underlying details of construction to the varying entity subclasses. A registry enhances this factory pattern by permitting the sharing of objects across domain boundaries. The use of factories and a data registry allows future developers to add new entity specializations in an easy and efficient manner. In figure 6, the factory pattern and the corresponding data registry are implemented as EntityFactory and EntityRegistry, respectively.

Zoom In Zoom Out Reset image size

3.4.1. Graph

The Graph data structure represents a set of vertices together with a set of edges coupling those vertices. Graph structures are common to the Problem, Processor, Embedding, and Program entities. The JadeD Graph model shown in figure 7 provides an abstraction of this structure in a way that promotes customization and extensibility with respect to a given entity type.

Zoom In Zoom Out Reset image size

In supporting this versatility, the Graph class utilizes two factory design patterns for generating vertices and edges [56]. This ensures object polymorphism by allowing custom subclasses to inject specialized edges and vertices. For example, this mechanism allows the production of static graphs for Problem, graphs that evolve in time for Program, and graphs that alter their state according to predefined conditions or controls for Processor.

3.4.2. Problem

The Problem class is a subclass of Entity that encapsulates the input data describing a discrete optimization problem. It is created by either an Analyst or Engineer in order to define the logical problem that the system will solve.

The current implementation of JadeD permits users to construct two distinct types of Problem. The first is a weighted or pseudo-Boolean optimization problem. The user inputs an arbitrary number of Boolean clauses in terms of the literals ${{b}_{i}}$, e.g., $\left( \left( {{b}_{1}}\;\;{\rm AND}\;\;{{b}_{2}} \right)\;\;{\rm OR}\;\;{\rm NOT}\;\;{{b}_{3}} \right)$, and each clause also has an associated real-valued weight ${{w}_{i}}$. The pseudo-Boolean function is then cast into an equivalent BOP by converting each Boolean literal to a corresponding binary variable, e.g., ${{b}_{i}}\mapsto {{x}_{i}}$, True $\mapsto 1$ and False $\mapsto 0$. The Boolean clauses are then recast into equivalent binary arithmetic expressions. Denoting the ith binary arithmetic clause as ${{f}_{i}}$ and the corresponding weight as ${{w}_{i}}$, the equivalent BOP over m bits is

$$\begin{eqnarray}&&{\rm arg} \mathop{{\rm min} }\limits_{{\bf x}}\mathop{\sum }\limits_{i}{{w}_{i}}{{f}_{i}}\left( {\bf x} \right),\end{eqnarray} \tag{ 15 }$$

where ${\bf x}\in {{\left\{ 0,1 \right\}}^{m}}$ is an m-bit vector [43]. In JADE, the BOP class stores both the original Boolean clauses and the reductions to algebraic expressions with corresponding weights.

The second type of Problem supported by JadeD is the QUBO problem defined in equation (7). For this type, the input corresponds to the elements of the matrix ${\bf P}$. The matrix ${\bf P}$ is then interpreted as a weighted adjacency matrix and parsed by JadeD into a Graph. Accordingly, the QUBO class is a subclass of Graph. The dependencies between the various Problem subclasses are illustrated in figure 8.

Zoom In Zoom Out Reset image size

As discussed in section 2, a BOP of the form in equation (15) can be reduced to a corresponding QUBO problem of the form in equation (7). The reduction, however, requires introduction of penalty terms to replace multilinear terms with quadratic or linear terms [43]. Expressing these penalties ultimately requires additional ancilla bits which enlarge the binary state space. When JadeD instantiates a BOP, the corresponding QUBO is immediately generated as part of the Problem. JadeD uses a QUBO reduction method that replaces the product of two binary variables by a new binary variables; the process repeats until a quadratic form remains [44]. The relevant BOP information is maintained as part of the Problem in order to facilitate developing the Solution entity returned to the Analyst.

3.4.3. Processor

3.4.4. Program

The Program class is a subclass of Entity that is used to synthesize specific instances of Problem and Processor into an implementation of the AQO algorithm. A Program is the primary input to the Sapphire simulation component and two different types can be constructed, physical or logical. The main difference between these two types for Program is the presence or absence of a high-level logical Problem definition.

Zoom In Zoom Out Reset image size

As shown in figure 10, type-switching is accomplished by composing Program with two classes: Logical Part and Physical Part. The physical part of a Program encapsulates the physical representation of the time-dependent Hamiltonian defined in equation (2). This includes a reference to a Processor and the parameters defining the final Ising Hamiltonian as well as the annealing schedule for each qubit. The logical part of a Program encapsulates a physical program as well as a reference to the specified Problem entity that is being solved. While the physical part of a Program entity is always required, the logical part is not. For Analyst use cases, the Program always has a logical part. In the absence of a logical input, the Program corresponds to an Engineer defined instance of an Ising Hamiltonian.

Zoom In Zoom Out Reset image size

The mapping of the Logical Part into the Physical Part generates an Embedding of the Problem into the Processor. As described in section 2, embedding generates a map between each logical vertex and a connected subgraph in the Processor. Within JadeD, this is accomplished using a subclass of Graph called Embedding. The Embedding class finds an embedding of the Logical Part into the provided Processor and Hardware. The current Embedding class supports the maximal minors methods described by Klymko et al. Its use is limited to a ${{K}_{4,4}}$, but the extensibility of Embedding means that the additional, greedy methods described by Klymko et al can also be incorporated [25].

The NiCE component is responsible for accepting user input, returning JADE output, and managing the computational workflow. It also provides a graphical front end for JADE. NiCE is an existing open-source project that was leveraged for reducing development time and ensuring extensibility. In addition to I/O management, the NiCE component orchestrates the interactions between the JadeD and Sapphire components. It enables users to create input files, launch simulations and examine program metrics.

NiCE is based on a client-server model, where the server handles primary data management and the client acts as the user front end. It is also possible for the server to manage remote workloads including, for example, simulations launched on remote hosts. We use the NiCE server as the primary means for launching and monitoring numerical simulations on both local and remote machines.

We have developed several plug-ins for NiCE that allow direct interaction with the JadeD component for the creation and revision of the Problem, Processor, and Program entities. A screenshot of one such NiCE form is provided in figure 11. NiCE is based on the Open Source Gateway Initiative framework that, among other things, permits dynamic registration of services. We use NiCEʼs implementation of dynamic registration to recognize and load user-defined plug-ins into JADE. This feature permits, for example, user-defined methods for simulation that are developed independently from JADE to be added during runtime. Additional information about NiCE is available from its website [59].

Zoom In Zoom Out Reset image size

Sapphire is the JADE component responsible for profiling Program entities. This includes carrying out numerical simulations of the quantum dynamics as well as other characterizations such as computing the time-dependent energy eigenspectra and computational error rates. While its primary use is to compute the Result of a Program, Sapphire permits a robust set of possible use cases. This is a result of our use of a plug-in architecture to support user-defined extensions to Sapphire. For example, numerical simulation techniques can be tailored to specific questions or physical assumptions. This promotes analysis at any desired fidelity and gives the user the ability to compare different simulation techniques against experimental benchmarks.

The extensibility of Sapphire is achieved through the interplay of a number of abstractions and design patterns, as shown in figure 12. Sapphire only exposes a few methods to external clients through the ISapphire interface. This decoupling between behavioral definition and actual implementation allows Sapphire to take on a number of varied forms. For example, JADE currently provides a Sapphire implementation for multi-threaded, shared memory architecture. We have also implemented SapphireMPI, which uses the Message Passing Interface (MPI) library to execute simulations on distributed architectures. The most significant difference between the two implementations is the MPI dependency and the need to perform unique initialization steps for SapphireMPI prior to beginning the numerical simulation.

Zoom In Zoom Out Reset image size

All implementations of Sapphire must define the method execute. When execute is invoked, Sapphire utilizes the JadeD file-parsing capabilities to construct the Program object defining the parameters of the numerical simulation. Sapphire next parses the simulations options provided by the user to create a Simulation object using the SimulationFactory. The Simulation class is the basis for the extensibility of Sapphire using plug-in libraries. A plug-in is essentially a subclass of Simulation that provides a specialized numerical or algorithmic approach to simulation.

The Simulation class is the primary unit of functionality within Sapphire and it is used to encapsulate a specific mathematical evolution of a quantum state. The factory design pattern allows Sapphire to remain completely agnostic to simulation details. However, there is a specific sequence of execution statements that are a necessary part of Sapphire. Program execution always begins with an initialization statement followed by a loop over a time-dependent solver. Once the exit condition is met, i.e., when t = T, the computational state undergoes readout before the program issues finalization commands. All plug-ins for Sapphire must adhere to the Simulation class functionality defined below.

• Initialize: this method is used primarily to initialize the quantum state of the simulation. Additional tasks include setting up any pre-simulation conditions or parameters.
• Anneal: this method is called every time step by Sapphire to advance the program quantum state. Developers should implement this method to update the state vector with the mathematics inherent to a specific technique for solving the time-dependent Schrödinger equation.
• QueryState: this method is used to query the state of the simulation, including the computational state of the simulated program. The output generated by this method is highly variable and it can include the internal representation of the quantum computational state or the complete eigenenergy spectrum written to an output file. These output files can also be used as checkpoints for restarting the simulation.
• Measure: this method is called after anneal completes and it represents measurement of the final computational state.
• Finalize: this method is used for any final calculations or clean up routines.

Developers of simulation plug-ins must subclass Simulation and implement the purely virtual anneal method. All other methods have default implementations that can be overwritten for specialized functionality. JADE also provides a specialized HamiltonianGenerator abstraction that permits decoupling of numerical dynamics from the actual form of the Hamiltonian describing the system.

3.7.1. Plug-in examples

The Sapphire plug-in architecture maintains extensibility to new simulation methodologies. A plug-in represents a user-created library that implements the Simulation class defined above. JADE users are therefore able to tailor quantum computing simulation techniques to specific problems or metrics of interest. We provide examples of plug-ins that implement Simulation below.

• SimulationZero: this plug-in provides a zero-th order approximation about the state of the computational register. Specifically, this simulation calculates the time-dependent eigenspectrum and instantaneous eigenstates of the time-dependent Hamiltonian defined by a Program. SimulationZero does not provide information about the quantum dynamics but essentially diagonalizes the Hamiltonian at each time step. This analysis provides information about the time-dependent energy gap. Our implementation is modeled in figure 13 and makes use of the Eigen library, which is an open-source C++ template library for linear algebra [60].
• RK4Simulation: this plug-in provides a fourth-order Runge–Kutta solver for the time-dependent Schrödinger equation as in equation (1). RK4Simulation uses two time steps, one for the outer anneal method which updates the Hamiltonian and a second for the inner evolve loop that numerically solves a finite-difference equation. For each evolve time step, the plug-in updates the quantum state and for each anneal it computes the instantaneous eigenspectrum. The plug-in also implements the queryState method to provide a Snapshot output that contains details about the computational state and eigenspectrum. Simulation options include the time steps, number of Snapshot files created, and number of eigenstates reported by queryState. This plug-in also makes use of the linear algebra functionality provided by the Eigen library.
• FOPSimulation: the FOPSimulation plug-in is based on a first-order perturbative solution to the time-dependent Schrödinger equation. It evolves a pure state according to a first-order Magnus expansion for the time-dependent propagation operator. Numerically, the propagation operator is diagonalized by the anneal method and applied successively to the state during the evolve method. This method has an error of $\mathcal{O}\left( \Delta {{t}^{3}} \right)$. Similar to the other simulation methods, Eigen is used to perform the matrix exponential and matrix-vector multiplications.

Zoom In Zoom Out Reset image size

The design and implementation of JADE relies heavily on test-driven development. A formal and rigorous testing model was defined before any actual product code was developed. This has ensured that (1) the functionality of each test unit was defined prior to its implementation and (2) the implementation of each source unit was fully compliant with the predetermined functionality. We employed test-driven development by modeling and designing surrogate classes whose sole purpose was for unit testing critical behavior in actual JADE classes. An example is shown in figure 14, where we test the Simulation class using surrogates for most objects in the Sapphire component. There is a corresponding SimulationTester class. Every class in JADE has a corresponding test class in order to provide the greatest assurance that the code adheres to design requirements.

Zoom In Zoom Out Reset image size