QuantumCircuitOpt: An Open-source Framework for Provably Optimal Quantum Circuit Design

Mixed-integer programming (MIP), a special case falling under an umbrella of mathematical programming, has been an integral part of important decision making problems such as the classic traveling salesman problem, airline scheduling and energy systems optimization, to name a few [12]. The mathematical theory of MIP, shown in (1), is concerned with finding amongst all solutions of a system of linear algebraic equalities and inequalities, written as functions of decision variables, a variable assignment that minimizes a linear function of those variables. The variables in an MIP can be a mix of continuous and integer values. Model (1) represents a canonical form of an MIP, where $x$ and $z$ are the continuous and integer variables, respectively. Solving MIPs, while NP-Hard in the worst case, has made dramatic strides in practical applications via state-of-the-art commercial solvers such as CPLEX and Gurobi, where decades of research and development (R&D) have yielded algorithmic performance improvements that even outpaced the hardware improvements from Moore’s Law. For example, CPLEX, since its release in 1991, has seen $\approx$200 billion factor speedup in solving MIPs [13].

At their core, MIPs are declarative, they specify what optimization problem to solve and not how to solve it. Hence, a performant mathematical programming language is required to access MIP solvers and to enable a separation from the problem-specific complexities and the underlying algorithms employed by such solvers. For this purpose, JuMP is one of the most widely adopted packages for mathematical programming and provides a declarative domain specific language (DSL) in the Julia programming language for modeling mathematical programs [14], such as MIPs. By taking a “leap of faith” that the problem of optimal quantum circuit design can also be posed as an MIP like shown in (1) (details in section III), we now present the framework of the QuantumCircuitOpt software package to formulate and solve such problems.

Building on the recent success of Julia, JuMP and mixed-integer programming, in this work we introduce QuantumCircuitOpt (or QCOpt), a free and open-source toolkit for quantum circuit design. As illustrated in Figure 1, QCOpt is written in Julia, a relatively new and fast dynamic programming language used for technical computing with support for extensible type system and meta-programming [15]. At a high level, QCOpt provides an abstraction layer to achieve two primary goals: (a) to capture user-specified inputs, such as a desired quantum computation and the available hardware gates, and build a JuMP model of an MIP formulation and (b) to extract, analyze and post-process the solution from the JuMP model to provide exact and approximate circuit decompositions, up to a global phase and machine precision. In the remainder of this section, we will highlight a few important features of QCOpt in detail.

User inputs
QCOpt’s user inputs is a dictionary, with a collection of key-value pairs, where the keys are the options to be set and values are the user inputs. These key-value pairs need not be of the same data type, as seen below. While the below list is not comprehensive, it provides sufficient information to start exploring the package:

1. 1.

   "num_qubits": Total number of qubits of the circuit.

2. 2.

   "maximum_depth": Maximum allowable depth for the circuit’s decomposition ($\geq$ 2). From QCOpt’s perspective, maximum depth is also the total number of input elementary_gates (see below) allowed in the decomposition.

3. 3.

   "elementary_gates": A vector of all one- and two-qubit elementary gates, native to the hardware. For example, on a three qubits circuit, here is an allowable input: ["H_1", "H_2", "T_2", "Tdagger_3", "CNot_1_2", "CNot_3_2", "CNot_3_1", "Identity"]. Notice that an one-qubit gate is associated with it’s qubit location and a two-qubit gate is associated with it’s control and target qubit locations. For example, inputs "CNot_1_2" and "CNot_3_2" implement the regular $\mathrm{CNOT}$ and the reverse $\mathrm{CNOT}$ gates in corresponding qubit locations.

4. 4.

   "target_gate": A target unitary quantum gate which needs to be compiled/decomposed using the above-mentioned native "elementary_gates".

5. 5.

   "objective": The options are (a) "minimize_depth" to minimize the total number of one- and two-qubit gates in the decomposition, and (b) "minimize_cnot" to minimize the total number of $\mathrm{CNOT}$ gates in the decomposition.

6. 6.

   "decomposition_type": The options are: (a) "exact": finds an exact, optimal circuit decomposition up to a global phase and machine precision, provided it exists, (b) "approximate": Finds an approximate decomposition with best gate fidelity found within the limits of run time.

7. 7.

   "set_cnot_lower_bound": While this is an optional input, it aids in reducing the size of the search space of MIP and also helps in exploiting theoretical bounds on $\rm CNOT$ gate costs from the literature [16].

Note that if a gate with continuous angles is in the input elementary_gates, then discretization of these angles have to be specified in the inputs. For example, see figure 2 for specifying the angles for the $\rm U_{3}$ gate.

Also, it is important to keep in mind that the optimization tasks handled by QCOpt are NP-hard to compute (in worst-case) in the size of the following set of inputs: num_qubits, maximum_depth and elementary_gates.

Currently, QCOpt supports the gates of the Clifford group ($\rm H,S,S^{\dagger},X,\sqrt{X},\sqrt{X}^{\dagger},Y,Z,CNOT$), and the $\pi/8$ phase gate, $\rm T$, and it’s conjugate transpose, $\rm T^{\dagger}$, which together are sufficient to arbitrarily approximate any target gate with a finite set of operations. However, for enhanced generality, users can also input universal single-qubit operators such as the rotation gates, $\rm(R_{x}(\theta),R_{y}(\theta),R_{z}(\theta)$) and the universal gate with three Euler angles, $\rm U_{3}(\theta,\phi,\lambda)$, along with the corresponding angle discretizations. QCOpt also supports various two-qubit gates such as the $\rm Swap$, $\rm iSwap$, $\rm Magic$-basis, $\rm Grover$’s diffusion operator, $\rm Sycamore$, $\rm QFT2$ and controlled versions of all the single-qubit gates in the Clifford group, $\rm T,T^{\dagger},V,V^{\dagger},R_{x},R_{y},R_{z}$ and $\rm U_{3}$ gates. QCOpt also supports Kronecker products of one- and two-qubit gates as a native gate in the following form: for e.g., on three-qubit circuits, ($\rm H_{1}\otimes T_{2}\otimes T_{3}$), ($\rm T_{1}\otimes CNOT_{2,3}$).

Figure 2 shows a sample implementation for decomposing a two-qubit controlled-$\rm Z$ gate. The native gates set contains discretized $\rm U_{3}$ gates which can be located on both the qubits, $\rm CNOT$ gate and the $\rm I$ gate. The maximum allowable depth for the decomposition is four and Gurobi is specified as an MIP solver. Finally, the results from the optimization model, including an optimal decomposition, run time and solver status, can be accessed from the results dictionary. QCOpt solves this optimization model to optimality, with 72 native gates (after discretization) per depth, in less than 4 seconds on a commodity computing hardware.

Pre-processing and transformation of mathematical programs before solving, known as pre-solving, is a well established procedure for improving run time of mathematical optimization algorithms [17]. To that end, before solving the main optimization model (MIP) to optimally, QCOpt preprocesses the user-defined inputs to make them more amenable to the MIP solver’s algorithm. Here we present a few such pre-solving techniques:

1. 1.

   Elimination of identical gates in the "elementary_gates" input, including the gates obtained via discretization of angle parameters.

2. 2.

   Implementation a compact, smaller-sized MIP formulation if both the input "elementary_gates" and "target_gate" contain purely real-valued elements.

3. 3.

   An optional user-given "input_circuit" is preprocessed to map this readily available circuit to warm-start the MIP solver with a feasible decomposition, as this can improve solver run times. During this pre-processing step, if QCOpt recognizes the input circuit as a feasible solution, then it activates the option to emphasize on proving optimality than finding feasibility for the MIP solver. For example, in the Gurobi solver, the MIPFocus parameter is set to “2” to emphasize optimality.

Currently QCOpt implements two separate objective functions relevant in the field of quantum circuit design. However, the software can easily be extended to support other objectives as needed, such as minimizing the cross-talk noise on quantum processors [25]. To be consistent with the previous derivations, we now present these functions on the two-qubit circuit:

a) To minimize the total depth, which is also the total number of one- and two-qubit gates admitted in the decomposition, as this serves as a proxy for the execution time of the circuit. This objective is modeled as a linear function as follows:

b) To minimize the total number of $\mathrm{CNOT}$ gates in the decomposition, modeled as a linear function as follows:

Note that it is important to have an $\rm Identity$ gate ($\rm I$) in the native gates set in order to obtain a compact circuit decomposition while minimizing the above mentioned objectives. This ensures an optimal decomposition of a depth lesser than the given value of $D$, if such a decomposition is feasible to the optimization model.

To summarize, the optimization model, with the linear objective function and the set of linearized constraints as described earlier in this section, is a mixed-integer program (MIP) which can be solved to global optimality using state-of-the-art MIP solvers. We now provide a brief discussion on the optimality guarantees of these solvers.

State-of-the-art solvers such as CPLEX and Gurobi are very well developed technologies with a solid mathematical footing and decades of R&D towards solving generic MIPs with convex constraints. These solvers employ a combination of branch-and-bound (B&B) algorithm and cutting plane generation to solve a MIP to global optimality [26, 12]. In its basic form, a B&B is a divide-and-conquer approach that partitions an exponential search space sequentially into smaller sub-problems, also referred to as nodes of the B&B tree, which are intended to be easier to solve. At each node of this B&B tree, a linear programming (LP) relaxation is solved by relaxing the integrality of the binary variables. During this process of “branching”, the LP relaxation’s objective serves as a lower bound and an integral solution’s objective, obtained using heuristics, serves an upper bound to the optimal solution of the MIP. If at any node, the LP relaxation value is greater than the best integral solution encountered so far, then the search can safely be stopped at that node as none of its children will yield a better solution than the current best solution. Termination criterion for the B&B algorithm is when the obtained LP relaxation is also integral, that is, when the lower and the upper bounds are close enough to each other, thus guaranteeing the global optimality of this integer solution. Despite the theoretical complexity of this algorithm being exponential in the worst-case, it can be accelerated by orders-of-magnitude by incorporation of effective and valid constraints (cutting planes or cuts) into this scheme. State-of-the-art MIP solvers incorporate an array of such generic cuts, like Gomory, mixed-integer rounding and flow-cover cuts, to name a few [12]. While these cuts perform well for generic MIPs, there still remain several technical insights and careful analysis which can lead to problem-specific cuts, which MIP solvers cannot generate automatically. Hence, in the forthcoming section, we summarize a few effective cuts which QCOpt implements to accelerate the convergence to global optimality.

QCOpt implements various symmetry-breaking, valid constraints in form of inequalities which have been observed to be very critical in reducing the search space or the number of nodes explored in the B&B tree. For example, these mathematical insights are shown in Section IV-A to provide up to a 180 times improvement in runtime of QCOpt even on two-qubit circuits. Here, we summarize these constraints:

1. 1.

   Commuting gates: Given any two unitary native gates $\mathrm{U}_{1}$ and $\mathrm{U}_{2}$, they are said to commute if their commutator, $[\mathrm{U}_{1},\mathrm{U}_{2}]=\mathrm{U}_{1}\cdot\mathrm{U}_{2}-\mathrm{U}_{2}\cdot\mathrm{U}_{1}$ is equal to zero. Thus including the following inequality in a pairwise fashion will forbid the symmetric feasible solution $U_{2}U_{1}$, while not changing the optimal circuit:

   $$z^{u_{2}}_{d}+z^{u_{1}}_{d+1}\leq 1\quad\forall d=1,\ldots,D$$

   Commuting gate constraints can also be further generalized to deeper sub-circuits of equivalent patterns, which can be found extensively in the quantum literature [9, 27, 28]. For example, while the sequence “$\mathrm{CNOT}_{1,2}\cdot\mathrm{H}_{2}\cdot\mathrm{CNOT}_{2,3}\cdot\mathrm{H}_{2}$” is a possible feasible sub-circuit, an equivalent sub-circuit “$\mathrm{H}_{2}\cdot\mathrm{CNOT}_{2,3}\cdot\mathrm{H}_{2}\cdot\mathrm{CNOT}_{1,2}$” can be eliminated from the search space via a valid inequality. QCOpt is currently equipped to identify such classic patterns and apply constraints to forbid equivalent sub-circuits.

2. 2.

   Involutory gates: Any unitary native gate, $\mathrm{U}$ is said to be involutory, if and only if $\mathrm{U}^{2}$ is an Identity ($\mathrm{I}$) gate. Assuming that $\mathrm{I}$ is already in the native gates set, such consecutive pairs of $\rm U$ gate can be eliminated using an inequality similar as above.

3. 3.

   Idempotent gates: Any unitary native gate, $\mathrm{U}$ is said to be idempotent, if and only if $\mathrm{U}^{2}=\mathrm{U}$. Again, such consecutive pairs of $\rm U$ gate can be eliminated using an inequality.

4. 4.

   Redundant gate-pairs: Let the set of native gates be $\mathcal{U}=\{\rm U_{1},U_{2},\ldots,U_{m}\}$. Any gate-pair index which belongs to the set $\{(i,j)|(i\neq j)\wedge(\mathrm{U}_{i}\cdot\mathrm{U}_{j}\in\mathcal{U}\ \forall i,j=1,\ldots,m)\}$ is considered as a redundant pair, and such pairs can be eliminated via inequalities as described above.

In the process of compiling any arbitrary entangled state into a circuit that is supported by a particular quantum processor, analyzing it’s feasibility with minimal measurement errors itself becomes a critical task. Moreover, obtaining insights on cost of noisy $\rm CNOT$ gates with feasibility certificates for an entangled target gate can also be very useful for quantum processor designer. To this end, QCOpt can be an incredibly useful diagnostic tool for proving infeasibility of decomposing any arbitrary unitary using hardware-based elementary gates. Although, proving infeasibility of an MIP can be NP-hard in worst case, there are numerous sophisticated techniques in MIP solvers such as conflict graph analysis, dual proof analysis, bound-tightening presolves and isolation of an Irreducible Infeasible Subset of constraints (IIS) from a large set of constraints [29]. These aforementioned methods can also get very efficient by further reducing the search space of feasibility via valid constraints as described in section III-C. Thus, the design philosophy of QCOpt makes it easy to perform this kind of infeasibility analysis, about which we present a result on the classic Toffoli gate in section IV-E.

Table I summarizes the run times (in seconds) of QCOpt’s implementation for decomposing various target gates as listed in the table. Note that the “Magic basis” and “Grover” gates in this table correspond to the circuits shown in Figures 4(b) and 5(b), respectively. The objective function for all these test gates was to minimize the total number of one and two-qubit gates in the resulting decomposition.

It is clear from Table I that the symmetry-breaking valid constraints (see section III-C) provide huge speedups in runtimes, with up to 180 times, in comparison with the basic MIP formulation derived in section III. However, note that the optimal solutions obtained by both these methods are identical, because the valid constraints are only speeding up the MIP solver by reducing the feasible search space. We also observed in our experiments that the basic MIP formulation (w/o VCs) stalled at much larger optimality gaps (without convergence) for circuits with $N\geq 3$. Thus, the MIP models and the valid constraints in combination with the pre-solving techniques discussed in section II-C make QCOpt an efficient tool for synthesizing up to medium-scale quantum circuits.

We now present new, compact circuit realizations provided by QCOpt on two-qubit circuits for target gates, “magic basis” and “Grover diffusion operator”. The circuits provided here are exact (up to global phase and machine precision) for each of these target gates. Run times for obtaining these realizations are provided in Table I (column 5). In this section, “compressed depth” of a circuit is the depth obtained by compressing the one-qubit gates on adjacent-depths which are acting on non-identical qubits into a single depth.

Figure 4 presents four equivalent circuit representations for the magic basis gate with different native gate sets. In the circuit in Figure 4(b), obtained using $\rm\{U_{3},CNOT_{2,1}\}$ native gates, we observe $\approx$33.3% reduction in the compressed depth in comparison with the circuit implemented in [7] (see Figure 4(a)).

Implementing an entangling gate like $\rm CNOT_{2,1}$ (like in Figure 4(a) and (b)) itself can be expensive, particularly in IBM-type hardware [31]. It usually necessitates five separate gates using a Hadamard gate on both the qubits and a $\rm CNOT_{1,2}$ gate. Hence, we implemented two equivalent circuits for the magic basis using $\rm\{S,H,CNOT_{1,2}\}$ and $\rm\{U_{3},CNOT_{1,2}\}$ native gates, respectively. Clearly, by optimizing using $\rm\{U_{3},CNOT_{1,2}\}$ gates, as seen in Figure 4(d), QCOpt produces a circuit with one less depth, which also is the most compact way to implement a magic basis gate. To summarize, it is worth to notice that this seemingly simple result in two qubits can provide a substantial improvement in the theoretical bounds for circuit complexity, as derived in [7], when the magic basis is repeatedly applied for realizing other two-qubit gates.

In Figure 5(b), the circuit is optimal with respect to the native gates in $\rm\{U_{3},CNOT_{1,2}\}$. To the best of our knowledge, this circuit with three gates is the most compact representation for the Grover operator in the literature. This realization has $\approx$57.1% lesser depth in comparison with the circuit provided in [32]. Again, to summarize, compact representations for fundamental two-qubit gates can be very beneficial, particularly when a gate like a Grover operator is applied repeatedly for larger problems as discussed in [32].

Figure 6 represents QCOpt’s optimal circuit realizations of target unitarites on three and four qubits, using only two-qubit native gates. We present these gates in comparison with [36], where the authors implement a brute-force-type enumeration algorithm, based on so-called “symbolic reachability analysis”, to realize optimal circuits. Thus the circuit realizations in Figure 6 obtained by QCOpt match the ones in [36]. Although, the run times mentioned in [36] were based on older computing platforms, we observe speedups in QCOpt by orders-of-magnitude, due to both improvement in computing platforms and the algorithmic enhancements implemented in QCOpt and state-of-the-art MIP solvers. To be clear, in Figure 6, QCOpt’s run times for decomposing Toffoli, Fredkin, double-Toffoli and quantum full-adder gates are 3.7 s., 72.1 s., 397.4 s., 348.7 s., respectively.

It is also noteworthy to mention that authors in [35] present an implementation for Fredkin gate and mention the following: “However, since the numerical search often gets stuck in local minima, even in cases where it eventually finds a solution, the fact that we were unable to find a smaller implementation of the Fredkin gate is not a proof that one does not exist.”. Using QCOpt, we prove that the decomposition in [35] for the Fredkin gate is indeed globally optimum (in $\approx$75 s.), as shown in Figure 6(b).

The naïve approach to decomposition of quantum circuits with optimality guarantees is a brute-force enumeration. This approach tries every possible combination of gates and therefore arrives at the optimal decomposition. To verify every possible circuit realization scales $O({N_{g}}^{D})$ for the native gate set of size $N_{g}$ and a depth of $D$. Additionally, the evaluation of each circuit scales with the complexity of matrix multiplication which grows exponentially with the number of qubits in the circuit. Thus, we implemented an approach which goes through all circuit combinations randomly up to a prescribed time limit. For this purpose, we set a time limit of 500 s.

For the Fredkin gate, the above-described approach could not produce any feasible decomposition (fidelity =1) within the time limit for all three randomized trials. However, for the quantum full-adder gate, out of three randomized trials, it produced an optimal decomposition in just one trial (within $\approx$435 s.), and was infeasible for the remaining two trials. For the quantum full-adder gate, as mentioned in section IV-C, QCOpt produces an optimal decomposition in less than $\approx$385 seconds. To summarize, in comparison with the above-described naïve approach, rigorous mathematical programming-based methods in QCOpt can be incredibly reliable and robust for quantum circuit applications.

As discussed in section III-D, QCOpt can not only can prove the existence of an optimal circuit, but can also be used as a tool to prove infeasibility of an existence of an implementation for given target gate. For this purpose, we tested the following: Given a three qubit Toffoli [16] as the target gate, with $\rm\{T_{1},T_{2},T_{3},T^{\dagger}_{1},T^{\dagger}_{2},T^{\dagger}_{3},CNOT_{1,2},CNOT_{2,3},CNOT_{1,3},I\}$ as the elementary gates, and a maximum depth for an allowable decomposition ($D$) equal to $10$, QCOpt provided a proof on infeasibility within $\approx$9 minutes that there does not exist an exact decomposition of depth lesser than or equal to $D$. This again shows the efficacy of mathematical models implemented in QCOpt.