globalization_mechanism globalization_strategy step_computation

# Binary Control Pulse Optimization for Quantum Systems

Xinyu Fei  Department of Industrial and Operations Engineering, University of
Michigan at Ann Arbor<EMAIL_ADDRESS>Lucas T. Brady  Joint Center for
Quantum Information and Computer Science, NIST/University of Maryland
<EMAIL_ADDRESS>Jeffrey Larson  Mathematics and Computer Science
Division, Argonne National Laboratory<EMAIL_ADDRESS>Sven Leyffer
Mathematics and Computer Science Division, Argonne National Laboratory
<EMAIL_ADDRESS>Siqian Shen  Department of Industrial and Operations
Engineering, University of Michigan at Ann Arbor<EMAIL_ADDRESS>
###### Abstract

Quantum control aims to manipulate quantum systems toward specific quantum
states or desired operations. Designing highly accurate and effective control
steps is vitally important to various quantum applications, including energy
minimization and circuit compilation. In this paper we focus on discrete
binary quantum control problems and apply different optimization algorithms
and techniques to improve computational efficiency and solution quality.
Specifically, we develop a generic model and extend it in several ways. We
introduce a squared $L_{2}$-penalty function to handle additional side
constraints, to model requirements such as allowing at most one control to be
active. We introduce a total variation (TV) regularizer to reduce the number
of switches in the control. We modify the popular gradient ascent pulse
engineering (GRAPE) algorithm, develop a new alternating direction method of
multipliers (ADMM) algorithm to solve the continuous relaxation of the
penalized model, and then apply rounding techniques to obtain binary control
solutions. We propose a modified trust-region method to further improve the
solutions. Our algorithms can obtain high-quality control results, as
demonstrated by numerical studies on diverse quantum control examples.

##  1 Introduction

Quantum control [1, 2, 3] is a rich field that started with applications in
quantum chemistry [4, 5, 6, 7, 8] but has since expanded to other subfields
such as atomic, molecular, and optical physics [9, 10, 11]. Quantum control
theory has been used in the quantum information community for a long time,
mostly for designing pulse controls and gates in quantum devices [12, 13, 14,
15, 16, 17, 18, 19, 20]; but more recently, researchers have considered using
control theories and optimization for the high-level design of quantum
algorithms [21, 22, 23, 24, 25, 26, 27]. In this paper, we focus primarily on
the quantum information applications of control theory, but our methods are
generalizable to other settings. We do not focus on the analytic aspects of
quantum control but on the more mechanistic aspects of how to apply
optimization techniques to the problem of designing and finding solutions to a
variety of quantum control problems.

Quantum computing architecture falls into different categories, such as
quantum circuits, which utilize discrete unitary gates to control a system,
and quantum annealing [28, 29], which uses smooth Hamiltonian evolution
instead. Quantum control is already familiar within quantum circuit
architecture design as a method for designing quantum gates, and several of
our examples described in Section 2 belong to this category. Within the
framework of quantum algorithms, quantum control is often used in variational
quantum algorithms, which can straddle the line between discrete gates and
smooth evolution. This paper focuses mainly on discrete binary control, which
is more applicable to a quantum gate architecture. Variational techniques have
already been implemented experimentally (see, e.g., [30, 31]), and our
techniques can potentially enhance the performance of the classical
variational loop on such near-term devices.

Most quantum optimal control algorithms are based on gradient descent for
better convergence than using gradient-free algorithms [32]. Khaneja et al.
[9] proposed the gradient ascent pulse engineering (GRAPE) algorithm for
designing pulse sequences in nuclear magnetic resonance. The authors
approximated the control function by a piecewise-constant function and
evaluated the explicit derivative. Based on the GRAPE algorithm, Larocca and
Wisniacki [33] proposed a K-GRAPE algorithm by using a Krylov subspace to
estimate quantum states during the process of time evolution. Brady et al.
[25] applied an analytical framework based on gradient descent to discuss the
optimal control procedure of a special control problem that minimizes the
energy of a quantum state with a combination of two Hamiltonians. Another
class of quantum optimal control algorithms is based on the chopped random
basis (CRAB) optimal control technique, which describes the control space by a
series of basis functions and optimizes the corresponding coefficients [34,
35]. Sørensen et al. [36] combined the GRAPE algorithm and CRAB algorithm to
achieve better results and faster convergence. However, all these algorithms
are designed for unconstrained continuous quantum control problems.

In this work, we focus on a quantum pulse optimization problem having binary
control variables and restricted feasible regions derived by linear
constraints. These constraints describe a so-called bang-bang control in a
model that corresponds to a quantum circuit design similar to that in the
quantum approximate optimization algorithm (QAOA) [37] and other variational
quantum algorithms [38, 39]. Some literature formulates the QAOA algorithm
into bang-bang control problems with a fixed evolution time and investigates
the performance of multiple methods including Stochastic Descent and
Pontryagin’s minimum principle for optimizing control models [40, 41].
However, they only consider quantum systems with two controllers, while we
propose a more general solution framework for quantum systems with multiple
controllers. In this paper, we introduce four quantum control examples: (i)
energy minimization problem, (ii) controlled NOT (CNOT) gate estimation
problem, (iii) NOT gate estimation problem with the energy leakage, and (iv)
circuit compilation problem. The nonconvexity, binary variables, and
restricted feasible sets in these examples lead to extreme challenges and
difficulties in solving the related binary quantum control problems.

Methods for binary control mainly include genetic algorithms [42], branch-and-
bound [43, 44], and local search [45]. To the best of our knowledge, genetic
algorithms have not been applied to binary quantum control, but Zahedinejad et
al. [46] showed that such algorithms fail to obtain high-quality solutions
even to continuous quantum control problems. Branch-and-bound can find high-
quality solutions for binary optimal control, but the long computational time
makes the algorithm intractable and hard to use for large-scale practical
problems. Vogt and Petersson [45] introduced a trust-region method [32] to
solve the single flux quantum control problem with binary controls. We extend
their approach by allowing the addition of constraints, such as min-up-time
and max-switching constraints to control the number of switches.

The main contributions of this paper are as follows. First, we propose a
generic model including continuous and discretized versions for binary quantum
control problems with linear constraints indicating that at each time step
there can be only one active control. Second, we introduce an exact penalty
function for linear constraints and develop an algorithmic framework combining
the GRAPE algorithm and rounding techniques to solve it. Third, to prevent
chattering on controls, we propose a new model that includes a total variation
(TV) regularizer, and then we propose a new approach based on the alternating
direction method of multipliers (ADMM) to solve this model. Fourth, we
introduce a modified trust-region method to improve the solutions obtained
from these algorithms. Compared with other methods, our algorithmic framework
can obtain high-quality binary control sequences and prevent frequent
switches. We demonstrate the performance on multiple quantum control examples
with various objective functions and controllers.

The paper is organized as follows. In Section 2 we construct generic models of
quantum binary control problems using continuous and discretized formulations.
We also introduce the corresponding specific formulations for four quantum
control examples. In Section 3 we review the GRAPE algorithm, utilize a
penalty function to relax certain linear constraints, and introduce our
rounding techniques. In Section 4 we propose and solve the penalty model with
the TV regularizer and the corresponding ADMM algorithm. In Section 5 we
derive trust-region subproblems and propose an approximate local-branching
method based on the trust-region algorithm. In Section 6 we present numerical
simulations for multiple instances of the three quantum control examples in
Section 2 to demonstrate the computational efficacy and solution performance
of our models and algorithms. In Section 7 we summarize our work and propose
future research directions.

##  2 Formulations for Quantum Control

Consider a quantum system with $q$ qubits and a time horizon $[0,t_{f}]$,
where $t_{f}>0$ is defined as the evolution time. Let
$H^{(0)}\in\mathbb{C}^{2^{q}\times 2^{q}}$ be the intrinsic Hamiltonian. Let
$N$ be the number of control operators and $H^{(j)}\in\mathbb{C}^{2^{q}\times
2^{q}},\ j=1,\ldots,N$ be the given control Hamiltonians. Let
$X_{\textrm{init}},\ X_{\textrm{targ}}\in\mathbb{C}^{2^{q}\times 2^{q}}$ be
the initial and target unitary operators of the quantum system, respectively.
Define variables $u_{j}(t)\in\left\\{0,1\right\\},\ j=1,\ldots,N,\ \forall t$
as the real-valued control functions at time $t$, and use $u$ to represent the
corresponding vector form. Define variables $H(t)\in\mathbb{C}^{2^{q}\times
2^{q}}$ as the time-dependent control Hamiltonian function and
$X(t)\in\mathbb{C}^{2^{q}\times 2^{q}}$ as the time-dependent unitary quantum
operator function. A generic quantum control problem for optimizing a function
of quantum operators is formulated as

$\displaystyle\min\quad$ $\displaystyle F(X(t_{f}))$ (1a) s.t. $\displaystyle
H(t)=H^{(0)}+\sum_{j=1}^{N}u_{j}(t)H^{(j)},\ \forall t\in[0,t_{f}]$ (1b)
$\displaystyle\frac{d}{dt}X(t)=-iH(t)X(t),\ \forall t\in[0,t_{f}]$ (1c)
$\displaystyle X(0)=X_{\textrm{init}}$ (1d)
$\displaystyle\sum_{j=1}^{N}u_{j}(t)=1,\ \forall t\in[0,t_{f}]$ (1e)
$\displaystyle u_{j}(t)\in\left\\{0,1\right\\},\ j=1,\ldots,N.$ (1f)

The objective function (1a) is a general function with respect to the final
unitary operator $X(t_{f})$. We assume that $F$ is a differentiable function.
In Sections 2.1–2.4, we provide the specific formulations of the objective
function for four quantum control examples. Constraint (1b) describes the
computation of the time-dependent Hamiltonian function based on the intrinsic
Hamiltonian, control Hamiltonians, and control functions. Constraint (1c) is
the ordinary differential equation describing the time evolution of the
operator of the quantum system with the same units for energy and frequency,
setting $\hbar=1$. Constraint (1d) is the initial condition of the operator.
Constraint (1e) indicates that the summation of all the control functions
should be one at all times, which is described by the Special Ordered Set of
Type 1 (SOS1) property in optimal control theory [47]. In binary control, this
SOS1 property guarantees that only one control field will be active at a time
that mimics the bang-bang nature of some quantum applications such as the
quantum approximate optimization algorithm and the variational quantum
eigensolver (VQE). Constraints (1f) require all time-based values of the
control functions to be feasible.

Following the time discretization process in [9], we explicitly integrate
constraint (1c). In particular, we divide the time horizon $[0,t_{f}]$ into
$T$ time intervals $[t_{k-1},t_{k}),\ k=1,\ldots,T$. We consider time
intervals with an equal length represented by $\Delta t$, but our work readily
extends to nonuniform discretizations. For each controller $j=1,\ldots,N$ and
each time step $k=1,\ldots,T$, we define discretized binary control variables
as $u_{jk}\in\left\\{0,1\right\\}$. For each time step $k=1,\ldots,T$, we
define the discretized time-dependent Hamiltonians as
$H_{k}\in\mathbb{C}^{2^{q}\times 2^{q}}$ and the quantum operators as
$X_{k}\in\mathbb{C}^{2^{q}\times 2^{q}}$. The differential equation (1c) is
thus approximated by

$\displaystyle\frac{d}{dt}X(t)=-iH_{k}X(t),\ \forall t\in[t_{k-1},t_{k}),\
k=1,\ldots,T.$ (2)

For each $k=1,\ldots,T$, the linear differential equation (2) is a Schrödinger
equation of a unitary operator, and we obtain an explicit solution as

$X(t)=\exp\left\\{-iH_{k}(t-t_{k-1})\right\\}X(t_{k-1}),\ \forall
t\in[t_{k-1},t_{k}),$

because $H_{k}$ is time independent. We then obtain a discretized quantum
control problem (DQCP) as a discretization of the differential control model
(1) using explicit solutions on each interval $[t_{k-1},t_{k})$ for unitary
operators as

$\displaystyle(DQCP)\quad\min\quad$ $\displaystyle F(X_{T})$ (3a) s.t.
$\displaystyle H_{k}=H^{(0)}+\sum_{j=1}^{N}u_{jk}H^{(j)},\ k=1,\ldots,T$ (3b)
$\displaystyle X_{k}=e^{-iH_{k}\Delta t}X_{k-1},\ k=1,\ldots,T$ (3c)
$\displaystyle X_{0}=X_{\textrm{init}}$ (3d)
$\displaystyle\sum_{j=1}^{N}u_{jk}=1,\ k=1,\ldots,T$ (3e) $\displaystyle
u_{jk}\in\left\\{0,1\right\\},\ j=1,\ldots,N,\ k=1,\ldots,T.$ (3f)

The objective function (3a) is the objective function (1a) evaluated at the
approximated final operator $X_{T}\approx X(t_{f})$. Constraints (3b)–(3f) are
the discretized formulations of constraints (1b)–(1f). In the following
sections we use this discretized formulation to develop our algorithms. We
present the following discussion on the relationship between the continuous
and discretized formulation.

######  Remark 1.

Problem (3) is equivalent to (1) for piecewise-constant controls $u_{j},\
j=1,\ldots,N$. In addition, we have

$\displaystyle\int_{0}^{t_{k}}u_{j}(\tau)d\tau=\sum_{\tau=1}^{k}u_{j\tau}\Delta
t,\ j=1,\ldots,N,\ k=1,\ldots,T.$ (4)

######  Remark 2.

The model can be generalized to instances with multiple active controllers in
any time interval by considering each possible combination of the controllers.
For an instance with $L$ control Hamiltonians $H^{(1)},\ldots,H^{(L)}$, we
consider $N=2^{L}$ combinations $H^{(l_{1},\ldots,l_{k})},\
\forall\left\\{l_{1},\ldots,l_{k}\right\\}\subseteq\left\\{1,\ldots,L\right\\}$
and define control functions for each combination. Thus, the original problem
is converted to a new problem with only one active controller in any time
interval.

In our paper we consider four examples with different objective functions and
control Hamiltonians. In the following sections we introduce the continuous
model of four example problems in quantum control. They can be formulated as
discretized models (3) following the above discretization process.

###  2.1 Energy Minimization Problem

Consider a quantum spin system consisting of $q$ qubits, no intrinsic
Hamiltonian, and two control Hamiltonians $H^{(1)},\ H^{(2)}$. The initial
operator $X_{\textrm{init}}$ is a $2^{q}$-dimensional identity matrix. Let
$\ket{\psi_{0}}$ be the ground state of the first control Hamiltonian
$H^{(1)}$. The ground state energy of this quantum system is the minimum
eigenvalue of the control Hamiltonian $H^{(2)}$, represented by
$E_{\textrm{min}}$ in our model. Since the ground state energy differs across
problem instances, it is important to compare performance consistently across
these instances. Therefore, we maximize the approximation ratio, defined as
the achieved energy divided by the true ground state energy. The problems we
study are well defined so that the initial state of the system has zero energy
with respect to $H^{(2)}$; furthermore, the ground state energy will always be
a large negative number. Therefore, the approximation ratio reflects the
energy level achieved by our procedure from the initial energy of zero to the
true ground state energy. Because our procedure can improve the state relative
only to its initial configuration, a negative approximation ratio is not
possible.

To ensure consistency with the objective functions of our other two examples,
we minimize one minus the approximation ratio. Because the true ground state
energy is negative, this formula is equivalent to minimizing the achieved
energy. The continuous-time formulation is

$\displaystyle\min\quad\displaystyle$ $\displaystyle
1-\left\langle\psi_{0}\right|X(t_{f})^{\dagger}H^{(2)}X(t_{f})\left|\psi_{0}\right\rangle/E_{\textrm{min}}$
(5a) s.t. $\displaystyle H(t)=u_{1}(t)H^{(1)}+u_{2}(t)H^{(2)},\ \forall
t\in[0,t_{f}]$ (5b) $\displaystyle H^{(1)}=-\sum_{i=1}^{q}\sigma_{i}^{x},\
H^{(2)}=\sum_{ij}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$ (5c) Constraints
(1c)–(1e) $\displaystyle u_{1}(t),\ u_{2}(t)\in\left\\{0,1\right\\},$ (5d)

where $\cdot^{\dagger}$ represents the conjugate transpose of a matrix and
$\langle\cdot|$ represents the conjugate transpose of a quantum state vector
$\ket{\cdot}$. The matrix $[J_{ij}],\ i,j=1,\ldots,q$ is the adjacency matrix
of a randomly generated graph with $q$ nodes, and $\sigma_{i}^{x},\
\sigma_{i}^{z}$ are Pauli matrices of qubit $i$.

###  2.2 CNOT Gate Estimation Problem

Our second example is defined on an isotropic Heisenberg spin-1/2 chain with
length 2 according to [48]. For our study we consider just the unitary
evolution with no energy leakage. The quantum system includes 2 qubits, an
intrinsic Hamiltonian $H^{(0)}$, and two control Hamiltonians
$H^{(1)},H^{(2)}$. The initial operator $X_{\textrm{init}}$ is a
$4$-dimensional identity matrix and the target operator is the matrix
formulation of the CNOT gate represented by
$X_{\textrm{CNOT}}\in\mathbb{C}^{4\times 4}$. This problem’s continuous-time
formulation is

$\displaystyle\min_{u(t),X(t),H(t)}\quad$ $\displaystyle
1-\frac{1}{4}\left|\operatorname{tr}\left\\{X_{\textrm{CNOT}}^{\dagger}X(t_{f})\right\\}\right|$
(6a) s.t. $\displaystyle H(t)=H^{(0)}+u_{1}(t)H^{(1)}+u_{2}(t)H^{(2)},\
\forall t\in[0,t_{f}]$ (6b) $\displaystyle
H^{(0)}=\sigma^{x}_{1}\sigma^{x}_{2}+\sigma^{y}_{1}\sigma^{y}_{2}+\sigma^{z}_{1}\sigma^{z}_{2}$
(6c) $\displaystyle H^{(1)}=\sigma^{x}_{1},\ H^{(2)}=\sigma^{y}_{1}$ (6d)
Constraints (1c)–(1d) $\displaystyle u_{1}(t),\
u_{2}(t)\in\left\\{0,1\right\\},$ (6e)

where $\sigma_{i}^{x},\ \sigma_{i}^{y},\ \sigma_{i}^{z},\ i=1,2$ are Pauli
matrices of two qubits. The goal is to minimize the objective function,
defined as the infidelity between the final operator and the target operator.

###  2.3 NOT Gate Estimation with Energy Leakage

We consider an anharmonic multilevel system that is widely used in quantum
experiments to implement physical qubits. Following [49], we consider the
lowest three levels of an energy spectrum with only nearest-level coupling for
a single qubit $q=1$. The first two levels comprise the qubit and the third
level models the energy leakage. Under this setting, the system contains an
intrinsic Hamiltonian $H^{(0)}$ and two control Hamiltonians $H^{(1)},\
H^{(2)}$. The initial operator $X_{\textrm{init}}$ is a 3-dimensional identity
matrix, and the target operator is the matrix formulation of the NOT gate
corresponding to the lowest three levels of the energy spectrum represented by
$X_{\textrm{NOT}}\in\mathbb{C}^{3\times 3}$, which is

$\displaystyle X_{\textrm{NOT}}=\begin{pmatrix}0&1&0\\\ 1&0&0\\\ 0&0&0\\\
\end{pmatrix}.$ (7)

The continuous-time formulation of this problem is

$\displaystyle\min_{u(t),X(t),H(t)}\quad$ $\displaystyle
1-\frac{1}{2}\left|\operatorname{tr}\left\\{X_{\textrm{NOT}}^{\dagger}X(t_{f})\right\\}\right|$
(8a) s.t. $\displaystyle H(t)=H^{(0)}+u_{1}(t)H^{(1)}+u_{2}(t)H^{(2)},\
\forall t\in[0,t_{f}]$ (8b) $\displaystyle H^{(0)}=\mu_{1}|1\rangle\langle
1|+\mu_{2}|2\rangle\langle 2|$ (8c) $\displaystyle
H^{(1)}={\omega_{1}}(|0\rangle\langle 1|+|1\rangle\langle
0|)+{\omega_{2}}(|1\rangle\langle 2|+|2\rangle\langle 1|)$ (8d) $\displaystyle
H^{(2)}={\omega_{1}}(i|0\rangle\langle 1|-i|1\rangle\langle
0|)+{\omega_{2}}(i|1\rangle\langle 2|-i|2\rangle\langle 1|)$ (8e) Constraints
(1c)–(1d) $\displaystyle u_{1}(t),\ u_{2}(t)\in\left\\{0,1\right\\},$ (8f)

where $|j\rangle$ represents a vector with the $j$-th element as 1 and other
elements as 0. The parameters $\mu_{1},\ \mu_{2}$ weigh the relative strength
of transitions and the parameters $\omega_{1},\ \omega_{2}$ correspond to the
drive frequency. We present the specific values of the parameters used in
numerical simulations in Section 6.2.

###  2.4 Circuit Compilation Problem

For a general quantum algorithm, each quantum circuit needs to be compiled to
a representation imposed by specific controllers and constraints in order to
execute the algorithm on specific quantum devices. We take quantum circuits
generated by the unitary coupled-cluster single-double (UCCSD) method [50, 51]
for estimating the ground state energy of molecules in quantum chemistry as
examples. We consider the Hamiltonian controllers in the gmon-qubit system
[52] with $q$ qubits. Each qubit has a flux-drive controller and a charge-
drive controller. In addition, there is a rectangular-grid topology with
nearest-neighbor connectivity denoted by $E$, and each pair of connected
qubits $(j_{1},j_{2})\in E$ has a controller. Define $w_{j_{1}j_{2}}(t)$ as
control functions for the connected qubits controllers and use $w(t)$ to
represent the corresponding vector form of $w_{j_{1}j_{2}}(t)$. The initial
operator $X_{\textrm{init}}$ is a $2^{q}$-dimensional identity matrix. The
target operator is the matrix formulation of the quantum circuit generated by
UCCSD for specific molecules represented by
$X_{\textrm{UCCSD}}\in\mathbb{C}^{2^{q}\times 2^{q}}$. We refer interested
readers to Gokhale et al. [53] for more details. The continuous-time
formulation of this problem is

$\displaystyle\min\quad$ $\displaystyle
1-\frac{1}{2^{q}}\left|\operatorname{tr}\left\\{X_{\textrm{UCCSD}}^{\dagger}X(t_{f})\right\\}\right|$
(9a) s.t. $\displaystyle
H(t)=\sum_{j=1}^{2q}u_{j}(t)H^{(j)}+\sum_{(j_{1},j_{2})\in
E}w_{j_{1}j_{2}}(t)H^{(j_{1}j_{2})},\ \forall t\in[0,t_{f}]$ (9b)
$\displaystyle H^{(2j-1)}=J_{c}\sigma^{x}_{j},\ j=1,\ldots,q$ (9c)
$\displaystyle H^{(2j)}=J_{f}\left(\begin{array}[]{ll}0&0\\\
0&1\end{array}\right),\ j=1,\ldots,q$ (9f) $\displaystyle
H^{(j_{1}j_{2})}=J_{e}\sigma^{x}_{j_{1}}\sigma^{x}_{j_{2}},\ (j_{1},j_{2})\in
E$ (9g) Constraints (1c)–(1e) $\displaystyle u_{j}(t)\in\left\\{0,1\right\\},\
w_{j_{1}j_{2}}(t)\in\left\\{0,1\right\\},\ j=1,\ldots,2q,\ (j_{1},j_{2})\in
E,$ (9h)

where $J_{c},\ J_{f},\ J_{e}$ are parameters corresponding to the quantum
system and $\sigma_{j}^{x}$ is the Pauli matrix of qubit $j$. We finish this
section by showing a general property of the objective function in optimal
control problems.

######  Theorem 1.

For a quantum system with $q$ qubits, the infidelity
$1-\frac{1}{2^{q}}\left|\operatorname{tr}\left\\{X_{\text{targ
}}^{\dagger}X(t_{f})\right\\}\right|\in[0,1]$. The objective functions of the
optimal control models for CNOT gate (6) and circuit compilation problem (9)
are bounded between zero and one.

###### Proof.

By definition, the objective function value cannot be larger than $1$. From
the definition of our quantum control problem, the quantum operators
$X_{\textrm{targ}}$ and $X(t_{f})$ are both unitary matrices. Therefore, we
can write $X(t_{f})$ as $X(t_{f})=\sum_{i=1}^{2^{q}}x_{i}e_{i}^{T}$, where
$x_{i}$ are the orthonormal column vectors of $X(t_{f})$ and $e_{i}$ are
$2^{q}\times 1$ vectors with the $i$th component as $1$ and all the other
components as $0$. Thus we have

$\displaystyle\left|\operatorname{tr}\left\\{X_{\text{targ
}}^{\dagger}X(t_{f})\right\\}\right|=\left|\sum_{i=1}^{2^{q}}\operatorname{tr}\left\\{X_{\text{targ
}}^{\dagger}x_{i}e_{i}^{T}\right\\}\right|\leq\sum_{i=1}^{2^{q}}\left\|X_{\text{targ
}}^{\dagger}x_{i}\right\|_{\infty}\leq 2^{q}.$ (10)

The first inequality follows from the fact that the trace of a product is
invariant under cyclic permutations of the factors. The second inequality
follows from $\left\|X_{\text{targ }}^{\dagger}x_{i}\right\|_{\infty}\leq 1$.
As a result, the objective function is no less than $0$, and this completes
the proof. ∎

######  Remark 3.

All the problems in this section fall into a class known as right-invariant
systems [54] that consists of control problems on Lie groups. Specifically,
this class consists of all control problems described by a Hamiltonian of the
form $H(t)=H^{(0)}+\sum_{j}u_{j}(t)H^{(j)}$. Controllability results on this
class of control problems are well known [54, 55]. A test of controllability
of these problems consists of examining the dynamical Lie algebra generated by
the Hamiltonians, $H^{(j)}$, and their nested commutators and comparing that
Lie algebra with the Lie algebra of the subgroup defined by the symmetries in
a given problem. For the systems considered in this paper, there are indeed
enough control Hamiltonians to reach the target state by using piecewise
continuous control functions, $u_{j}(t)$.

##  3 Binary Quantum Control Algorithms

In this section we propose models and algorithms to obtain binary quantum
control results for the discretized model (3). In Section 3.1 we review the
adjoint approximation method of gradients used in the GRAPE algorithm for
solving the continuous relaxation of the optimal control problem (3) and
derive the gradient update rule for the aforementioned four specific examples.
In Section 3.2 we propose an exact penalty function for the SOS1 property that
allows us to add constraints to GRAPE. In Section 3.3 we introduce rounding
techniques to obtain binary controls. These algorithms form the basis of our
discrete framework discussed in Sections 4–5.

The overall framework is outlined in Figure 1, which presents the methods in
Sections 3–5 and the corresponding tables and results. We summarize the
algorithms of our framework in Table 1. In the figure, the rectangles
represent the three formulations of our model, and the ellipses represent the
algorithms. The parallelograms show the result tables corresponding to
algorithms, where the dashed parallelograms refer to the results of the
continuous models and the solid parallelograms refer to the results of the
discrete models.

Unconstrained model (continuous relaxation of (3) without (3e)) Model with
penalized squared $L_{2}$-norm (continuous relaxation of (3.2)) Model with TV
regularizer (continuous relaxation of (DQCP-TV)) GRAPE (Section 3.1) Penalized
GRAPE (Section 3.2) ADMM algorithm (Section 4) Trust-region method (Section 5)
Table 4 Sum-up-rounding (no constraints, Section 3.3.1) Integral minimization
with min-up-time constraints (Section 3.3.2) Integral minimization with max-
switching constraints (Section 3.3.2) Table 6 Approximate Local-branching
(Section 5) Table 8

Figure 1: Overview and results of the algorithmic framework. Table 1: List of
acronyms of algorithms in numerical simulations. The pGRAPE and GRAPE
algorithms are identical when the SOS1 constraint is absent.

Name Algorithms pGRAPE Solving the continuous relaxation with squared
$L_{2}$-penalty function by penalized GRAPE. TR Solving the continuous
relaxation with TV regularizer and SOS1 property by a trust-region method.
ADMM Solving the continuous relaxation with TV regularizer and SOS1 property
by ADMM. SUR Rounding the continuous solutions without hard control
constraints. MT Rounding the continuous solutions with min-up-time
constraints. MS Rounding the continuous solutions with max-switching
constraints. ALB Improving the binary solutions by the approximate local-
branching method.

###  3.1 GRAPE Approach for Solving Continuous Relaxations

The GRAPE algorithm [9] is a gradient descent algorithm for solving the
unconstrained continuous control problem. We derive a gradient descent
algorithm based on the adjoint approach used in the GRAPE algorithm to solve
the continuous relaxation of our discretized model (3) without the SOS1
property enforced by constraint (3e). We eliminate variables $X$ and $H$ by
converting them into implicit functions of the control variables $u$ using
constraints (3b)–(3d). Then the minimization problem of $F(X_{T})$ on $u,\ X,\
H$ is converted to the minimization problem of $F(X_{T}(u))$ over the
variables $u$. The goal of the adjoint approach is to approximate the gradient
of the objective function $F$ with respect to control variables $u_{jk},\
j=1,\ldots,N,\ k=1,\ldots,T$ based on the discretized time horizon. For
simplicity, we define propagators for each time step $k=1,\ldots,T$ as
$U_{k}=\exp\left\\{-iH_{k}\Delta t\right\\}$. For small $\Delta t$, the
gradient of $U_{k}$ corresponding to $u_{jk}$ is estimated as

$\displaystyle\frac{\partial U_{k}}{\partial u_{jk}}=-i\Delta tH^{(j)}U_{k},\
j=1,\ldots,N,\ k=1,\ldots,T.$ (11)

The gradient of the objective function with respect to the propagators depends
on its specific formulation. We discuss gradients of two specific functions
for the examples in Section 2.

##### Energy function

For generality, we use $\bar{H}$ to represent the Hamiltonian in the energy
objective function, which means that the goal is to minimize the energy
corresponding to $\bar{H}$. Specifically, we have $\bar{H}=H^{(2)}$ in the
example in Section 2.1. The general energy function can be expressed as

$\displaystyle\langle\psi_{0}|X_{T}^{\dagger}\bar{H}X_{T}\ket{\psi_{0}}=\langle\psi_{0}|U_{1}^{\dagger}\cdots
U_{T}^{\dagger}\bar{H}U_{T}\cdots U_{1}\ket{\psi_{0}}.$ (12)

For each time step $k=1,\ldots,T-1$, define variables
$\ket{\kappa_{k}}=U_{k+1}^{\dagger}\cdots
U_{T}^{\dagger}\bar{H}X_{T}\ket{\psi_{0}}$. For the last time step $T$, we
define a variable $\ket{\kappa_{T}}=\bar{H}X_{T}\ket{\psi_{0}}$. Then the
gradient with respect to $u_{jk}$ is computed as

$\displaystyle\frac{\partial F}{\partial
u_{jk}}=\frac{2}{E_{\textrm{min}}}\textrm{Re}\left[i\Delta
t\langle\kappa_{k}|H^{(j)}X_{k}\ket{\psi_{0}}\right],\ j=1,\ldots,N,\
k=1,\ldots,T.$ (13)

##### Infidelity function

The infidelity function can be expressed as

$\displaystyle
1-\frac{1}{2^{q}}\left|\operatorname{tr}\left\\{X_{\textrm{targ}}^{\dagger}X_{T}\right\\}\right|=1-\frac{1}{2^{q}}\left|\operatorname{tr}\left\\{X_{\textrm{targ}}^{\dagger}U_{T}\cdots
U_{1}X_{0}\right\\}\right|.$ (14)

For each time step $k=1,\ldots,T-1$, define variables
$\lambda_{k}=U_{k+1}^{\dagger}\cdots U_{T}^{\dagger}X_{\textrm{targ}}$. For
the last time step $T$, we define a variable $\lambda_{T}=X_{\textrm{targ}}$.
Then the gradient of the trace with respect to $u_{jk}$ is computed as

$\displaystyle\frac{\partial\operatorname{tr}\left\\{X_{\textrm{targ}}^{\dagger}X_{T}\right\\}}{\partial
u_{jk}}=-i\Delta t\lambda_{k}^{\dagger}H^{(j)}X_{k},\ j=1,\ldots,N,\
k=1,\ldots,T.$ (15)

Using the definition of objective function $F$, we compute the gradient as

$\displaystyle\frac{\partial F}{\partial
u_{jk}}=\frac{1}{2^{q}}\textrm{Re}\left[i\Delta
t\,\textrm{tr}\left\\{\lambda_{k}^{\dagger}H^{(j)}X_{k}\right\\}e^{-i\textrm{arg}\left(\operatorname{tr}\left\\{X_{\textrm{targ}}^{\dagger}X_{T}\right\\}\right)}\right],\
j=1,\ldots,N,\ k=1,\ldots,T,$ (16)

where $\textrm{arg}(\cdot)$ represents the argument of a complex number.

With these computed gradients, we minimize the reduced function $F(X_{T}(u))$
with respect to control variables $u\in[0,1]^{N\cdot T}$ by L-BFGS-B [56], a
well-known quasi-Newton algorithm for solving bound-constrained optimization
problems.

######  Remark 4.

Our implementation of GRAPE extends the classical gradient-ascend scheme (see,
e.g., [9, 33]) in a number of ways: We use a bound-constrained quasi-Newton
method to solve problems with bounded control variables, and we introduce a
penalized GRAPE (denoted by pGRAPE) that adds a quadratic penalty of the SOS1
property in Section 3.2.

We take the circuit compilation problem (Section 2.4) on the molecule H2
(dihydrogen) as an example to show the control results, and we present the
continuous control results without the SOS1 property obtained by the GRAPE
algorithm in Figure 2a. We define the absolute violation of the SOS1 property
at each time step $k=1,\ldots,T$ as $|\sum_{j=1}^{N}u_{jk}-1|$ and present the
value at each time step in Figure 2b. The figure shows that the results
violate the SOS1 property at all time steps, with the maximum absolute
violation being $2.404$. Hence we introduce a penalized function in Section
3.2 to impose the SOS1 property.

(a) Control results of $5$ controllers

(b) Absolute violation at each step

Figure 2: (a) Control results of the continuous relaxation of the discretized
model (3) without the SOS1 property constraint (3e) of the circuit compilation
example on the molecule H2. (b) Absolute violation at each time step.

###  3.2 Penalty Function for SOS1 Property

Our GRAPE algorithm is defined only for unconstrained or bound-constrained
problems [9]. Hence, we convert the SOS1 property (3e) into a squared
$L_{2}$-penalty term and add it to the objective function. Define the squared
$L_{2}$ function of the constraint violation as

$\displaystyle l(u,T)=\sum_{k=1}^{T}\left(\sum_{j=1}^{N}u_{jk}-1\right)^{2}.$
(17)

We let a constant $\rho>0$ be the penalty parameter. Then the model with the
penalty function is

$\displaystyle\min\ F(X_{T})+\rho l(u,T)\quad\textrm{s.t.}\
\textrm{Constraints \eqref{eq:model-d-1-cons-h}--\eqref{eq:model-d-1-cons-u-
sum}}.$

We use $l(u^{*}_{\rho},T)$ to denote the optimized value of the squared
$L_{2}$ function $l(u,T)$ under the penalty parameter $\rho$. The following
theorem discusses the exactness of the penalized objective function.

######  Theorem 2.

[Exactness of Squared $L_{2}$-Penalty] Let $z_{1}$ and $z_{2}(\rho)$ be the
optimal values of (3) and (3.2), respectively. There exists
$\tilde{\rho}<\infty$ such that $z_{1}=z_{2}(\rho)$ for all
$\rho>\tilde{\rho}$.

A more general version of Theorem 2 is proved in [57]. Next we discuss the
value of the optimized squared $L_{2}$-penalty term with respect to the
penalty parameter $\rho$ for the continuous relaxation of (3.2).

######  Theorem 3.

For the continuous relaxation of (3.2), if the original objective function $F$
is upper bounded by a constant $C_{F}$ and (3) is feasible, then for any $T$
the optimized squared $L_{2}$-penalty term uniformly holds that
$|l(u^{*}_{\rho},T)|\leq 2C_{F}/\rho$.

A detailed proof is provided in Appendix A.1.

######  Remark 5.

The GRAPE algorithm is still able to solve the unconstrained continuous
relaxation with the penalty function (3.2). We propose a penalized version
(pGRAPE) in Section 3.1 by adding a term
$2\rho\left(\sum_{j=1}^{N}u_{jk}-1\right)$ to the approximated gradient of the
objective function.

######  Remark 6.

For the quantum control problem with two controllers, the SOS1 can be enforced
directly by substituting

$\displaystyle u_{2k}=1-u_{1k},\ k=1,\ldots,T$ (18)

in the model instead of using the penalty function.

We show the control results obtained from the continuous relaxation with the
squared $L_{2}$-penalty function and penalty parameter $\rho=1.0$ in Figure 3.
The value of the squared penalty term is $5.55\times 10^{-7}$, leading to a
small error between continuous and binary results with the SOS1 property
obtained by rounding techniques. However, we still need to provide rounding
techniques to obtain binary results.

Figure 3: Control results of the continuous relaxation of the discretized
model with a squared $L_{2}$-penalty function (3.2) of the circuit compilation
example on the molecule H2 (objective value: 4.37E$-$07).

###  3.3 Rounding Techniques and Optimality Guarantees

In this section we introduce rounding techniques to obtain binary control
results and investigate their optimality guarantees. In Section 3.3.1 we
review the sum-up rounding technique designed for binary controls without
constraints and discuss the difference compared with the continuous results.
In Section 3.3.2 we introduce the integral minimization problem for rounding
to obtain restricted binary controls.

####  3.3.1 Sum-Up Rounding

The sum-Up rounding (SUR) strategy proposed by Sager et al. [47] is a well-
known method to obtain integer controls from continuous, relaxed ones in
optimal control theory; see, for example, [58, 59]. One concern when applying
SUR is that the solution of the relaxation does not satisfy the SOS1
constraints because the exactness of the squared $L_{2}$-penalty is only for
binary controls. We show that as long as the violation of the SOS1 constraint
is small, the solution that is constructed via SUR still satisfies the strong
convergence properties of SUR. We define the discretized continuous control
function $u^{c}$ as $u^{c}_{jk}\in[0,1],\ j=1,\ldots,N,\ k=1,\ldots,T$. We
define the discretized binary control function $u^{b}$ as
$u^{b}_{jk}\in\left\\{0,1\right\\},\ j=1,\ldots,N,\ k=1,\ldots,T$. For each
time step $k=1,\ldots,T$, we denote the vector form of control variables as
$u^{c}_{k}$ and $u^{b}_{k}$. We define the current cumulative deviation
between continuous and binary controls as $\hat{p}_{jk},\ j=1,\ldots,N,\
k=1,\ldots,T$. The SUR approach to obtain binary controls is described in
Algorithm 1.

Input: Continuous control $u^{c}$ on uniform discretization.

/* Iterate over each time step */

for  _$k=1,\ldots,T$_ do

/* Iterate over each controller */

for  _$j=1,\ldots,N$_ do

Compute cumulative deviation as
$\displaystyle\hat{p}_{jk}=\sum_{\tau=1}^{k}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{k-1}u^{b}_{j\tau}\Delta t$.

Choose controller $j^{*}=\mathop{\rm arg\,max}_{j=1,\ldots,N}\hat{p}_{jk}$. If
there is a tie, we break the tie by choosing the smallest index.

Set binary control $u^{b}_{j^{*}k}=1$ and $u^{b}_{jk}=0,\ \forall j\neq
j^{*}$.

Output: Binary control $u^{b}$.

Algorithm 1 Sum-Up Rounding with the SOS1 Property.

The construction of $u_{jk}^{b}$ ensures that the binary control satisfies the
SOS1 property. Sager et al. [47] proved that if the continuous control $u^{c}$
satisfies the SOS1 property, then the rounded control $u^{b}$ will converge to
$u^{c}$ when the length of a time interval $\Delta t$ converges to zero. This
does not hold for arbitrary continuous controls without the SOS1 property,
however, as the next proposition shows.

######  Proposition 1.

Let continuous control $u^{c}$ and binary control $u^{b}$ defined in Algorithm
1 be given for $N\geq 2$. Define

$\displaystyle\epsilon(\Delta
t)\coloneqq\max_{k=1,\ldots,T}\left|\sum_{\tau=1}^{k}\left(\sum_{j=1}^{N}u^{c}_{j\tau}-1\right)\Delta
t\right|.$ (19)

If $t_{f}<\infty$, then $u^{b}$ satisfies the SOS1 property
$\sum_{j=1}^{N}u^{b}_{jk}=1$ at any time step $k=1,\ldots,T$, and we have

$\displaystyle\max_{k=1,\ldots,T}\left\|\sum_{\tau=1}^{k}\left(u^{c}_{\tau}-u^{b}_{\tau}\right)\Delta
t\right\|_{\infty}\geq\frac{1}{N}\epsilon(\Delta t).$ (20)

The proof is provided in Appendix A.2. We have the following estimate of the
difference between the continuous control without the SOS1 property and the
rounded control.

######  Theorem 4.

With the same definition of $\epsilon(\Delta t)$ in Proposition 1, it holds
that for any time step $k=1,\ldots,T$,

$\displaystyle\left\|\sum_{\tau=1}^{k}\left(u^{c}_{\tau}-u^{b}_{\tau}\right)\Delta
t\right\|_{\infty}\leq\left(N-1\right)\Delta t+\frac{2N-1}{N}\epsilon(\Delta
t).$ (21)

Because $u^{b}$ and $u^{c}$ are bounded, we have $\epsilon(\Delta t)<\infty$
if $t_{f}<\infty$. The proof of the inequality (21) extends the proof for
Theorem 5 in [47]. We modify the right-hand side of the inequality by adding a
term $\frac{2N-1}{N}\epsilon(\Delta t)$. We provide the detailed proof in
Appendix A.2. In the following corollary we show that for our bounded
discretized quantum control problem, $\epsilon(\Delta t)$ converges to zero as
$\Delta t$ converges to zero, ensuring the convergence of the rounded control
to the continuous control.

######  Corollary 1.

Let $l(u^{c},T)$ be the value of the optimized squared $L_{2}$ term in the
continuous relaxation of the model (3.2). Then it holds that

$\displaystyle\epsilon(\Delta t)\leq\sqrt{t_{f}{l(u^{c},T)}\Delta t},$ (22)

where $t_{f}$ is the evolution time. Furthermore, if the original objective
function $F$ is bounded, the rounded control will converge to the continuous
control when $\Delta t$ is small enough.

The proof is provided in Appendix A.2. Based on the convergence of binary
results to continuous results, we present the following proposition to
guarantee the optimality of binary results.

######  Proposition 2.

Under the discretized setting of the quantum control problem, let $u^{c}$ be
the continuous control and $u^{b}$ be the binary control obtained by Algorithm
1. Then the state of time evolution with binary control $X^{b}$ converges to
the state with continuous control $X^{c}$ at each time step, namely,

$\displaystyle\lim_{\Delta t\rightarrow 0}X^{b}_{k}=\lim_{\Delta t\rightarrow
0}X^{c}_{k},\ k=1,\ldots,T,$ (23)

leading to $\displaystyle\lim_{\Delta t\rightarrow 0}F(X_{T}^{b})=\lim_{\Delta
t\rightarrow 0}F(X_{T}^{c})$ for a continuous objective function.

The proof is provided in Appendix A.2.

####  3.3.2 Combinatorial Integral Approximation

To obtain discretized binary controls, Sager [60] proposed a more general
rounding technique by minimizing the integral difference between continuous
and binary controls with certain additional constraints on the binary controls
called combinatorial integral approximation (CIA). Let
$\mathcal{U}_{B}\subseteq\left\\{0,1\right\\}^{N\cdot T}$ be the feasible
region for binary controls. For each controller $j=1,\ldots,N$ and each time
step $k=1,\ldots,T$, define $u^{c}_{jk}$ and $u^{b}_{jk}$ as the discretized
continuous and binary controls. We use $u^{c}$ and $u^{b}$ to represent the
corresponding vector forms. The optimization problem for rounding is
formulated as a mixed-integer problem

$\displaystyle\min_{u^{b}}\quad$
$\displaystyle\max_{j=1,\ldots,N}\max_{k=1,\ldots,T}\left|\sum_{\tau=1}^{k}(u^{c}_{j\tau}-u^{b}_{j\tau})\Delta
t\right|$ (24a) s.t. $\displaystyle u^{b}\in\mathcal{U}_{B}.$ (24b)

The rounding optimization problem can be solved by diverse algorithms for
solving integer programs. In this paper we choose a branch-and-cut algorithm
[61]. Furthermore, the rounding result obtained from SUR is an optimal
solution of model (24) with $\mathcal{U}_{B}=\left\\{0,1\right\\}^{N\cdot T}$
[60]. In practice, researchers and engineers have proposed diverse
restrictions on binary controls to avoid frequent switches. We formulate two
main types of constraints as linear constraints.

##### Min-Up-Time Constraints.

Let $T_{\textrm{minup}}$ be the minimum number of steps $\Delta t$ that the
controller is active. The min-up-time constraints enforce that each controller
is active for at least $T_{\textrm{minup}}$ time steps. The restricted
feasible region is formulated by the following linear constraints:

$\displaystyle\mathcal{U}_{B}=\Bigg{\\{}(u,v):$ $\displaystyle-v_{jk}\leq
u_{jk}-u_{jk+1}\leq v_{jk},\ j=1,\ldots,N,\ k=1,\ldots,T-1$
$\displaystyle\sum_{k=t}^{t+T_{\textrm{minup}}-1}v_{jk}\leq 1,\ j=1,\ldots,N,\
t=1,\ldots,T-T_{\textrm{minup}}\Bigg{\\}}.$ (25)

##### Max-Switching Constraints.

Let $S$ be the maximum number of switches during the evolution time horizon
$[0,t_{f}]$. The max-switching constraints enforce an upper bound $S$ on the
total number of switches for each controller. We describe the restricted
feasible region by linear constraints as

$\displaystyle\mathcal{U}_{B}=\Bigg{\\{}(u,v):$ $\displaystyle-v_{jk}\leq
u_{jk}-u_{jk+1}\leq v_{jk},\ j=1,\ldots,N,\ k=1,\ldots,T-1$
$\displaystyle\sum_{k=1}^{T}v_{jk}\leq S,\ j=1,\ldots,N\Bigg{\\}}.$ (26)

We present the binary results obtained from SUR and CIA with min-up-time
constraints in Figure 4. We show that SUR leads to chattering on singular
arcs. Although CIA can prevent frequent switches by setting hard constraints
on the controls, it leads to a serious objective value increase. In Section 4
and Section 5 we propose models and algorithms to reduce switches when taking
the original objective function into account.

(a) SUR results (objective value $0.021$)

(b) CIA results with min-up time constraints (objective value $0.600$)

Figure 4: Binary control results obtained by conducting SUR and CIA on
continuous results of the circuit compilation example using the molecule H2.
The min-up time constraints reduce the switches but increase the objective
value.

##  4 Model with the TV Regularizer

Binary controls that avoid chattering on singular arcs can be obtained by
solving the rounding model (24). However, the objective function of the
rounding model is not relevant to the original objective function $F$, leading
to a significant increase in the value of $F$, as illustrated in Figure 4.
Here we investigate a TV regularizer as an alternative strategy to obtain
solutions with fewer switches.

The TV of a function is defined as the integral of the absolute change of the
function over the entire space. Since Rudin et al. [62] first introduced the
TV regularizer by considering functions describing images, it has become a
popular method in image noise reduction; see, for example, [63, 64]. We refer
interested readers to [65] for a detailed review of applying a TV regularizer
to denoise images. Stella et al. [66] propose an algorithm to solve the
continuous control problem with the TV regularizer based on the forward-
backward envelope [67]. Sager and Zeile [68] considered the TV of integer
control functions to reduce the absolute change of controls. The authors added
a constraint imposing that the TV regularizer should be no more than a
threshold. However, they considered the TV regularizer only when rounding the
control results of continuous relaxation, which leads to an increase in the
original objective function. Leyffer and Manns [69] extended the trust-region
method to solve the integer optimal control problem with the TV regularizer in
the objective function, but it is still highly dependent on the initial
values.

We use the TV regularizer to penalize the absolute change between control
variables of two consecutive time steps, which is defined as

$\displaystyle TV(u)=\sum_{j=1}^{N}\sum_{k=1}^{T-1}|u_{jk}-u_{jk+1}|$ (27)

for a control variable vector $u$. Let $\alpha>0$ be the parameter for the TV
regularizer. Then the continuous relaxation model with the TV regularizer is
formulated as

$\displaystyle\min\ $ $\displaystyle F(X_{T})+\rho l(u,T)+\alpha
TV(u)\quad\textrm{s.t.}\ \textrm{Constraints
\eqref{eq:model-d-1-cons-h}--\eqref{eq:model-d-1-cons-i}},\ u\in[0,1]^{N\cdot
T}.$ (DQCP-TV)

The first term is the original objective function, the second term is the
squared $L_{2}$-penalty function for the SOS1 property, and the third term is
the TV regularizer function. Because of the $L_{1}$ term, we cannot use GRAPE;
instead, we propose ADMM to solve the continuous relaxation. For each
controller $j=1,\ldots,N$ and time step $k=1,\ldots,T-1$, we introduce
auxiliary variables $v_{jk}=u_{jk}-u_{jk+1}$ to describe the change between
control variables of two consecutive time steps. We reformulate the model as
follows:

$\displaystyle\min\quad$ $\displaystyle F(X_{T})+\rho
l(u,T)+\alpha\sum_{j=1}^{N}\sum_{k=1}^{T-1}|v_{jk}|$ (28a) s.t. $\displaystyle
v_{jk}=u_{jk}-u_{jk+1},\ j=1,\ldots,N,\ k=1,\ldots,T-1$ (28b) Constraints
(3b)–(3d) $\displaystyle u\in[0,1]^{N\cdot T}.$ (28c)

Define $\mu_{jk},\ j=1,\ldots,N,\ k=1,\ldots,T-1$ as the dual variables
corresponding to constraints (28b). We use $u,\ v,\ \mu$ to denote the
corresponding vector forms. Let fixed parameters $\beta>0,\delta>0$ be the
Lagrangian penalty parameter and stopping criterion threshold. The update
procedure of ADMM consists of three steps: (i) updating variables $u$, (ii)
updating variables $v$, and (iii) updating dual variables $\mu$. We solve the
minimization problem for updating variables $u$ by the modified GRAPE
algorithm with a squared $L_{2}$-penalty function and Lagrangian penalty
function. We derive an exact form for the update of variables $v$. We use
gradient descent to update the dual variables $\mu$. The specific procedure
for updating is presented in Algorithm 2.

Input: Initial values of variables $u^{0},\ v^{0},\ \mu^{0}$, and number of
ADMM iterations $L$.

Initialize iteration $l=1$.

while  _$l\leq L$ and the algorithm does not converge_ do

Update variables $u$ as  $\displaystyle u^{l}=$ $\displaystyle\mathop{\rm
arg\,min}_{u\in\mathcal{U}_{C}}\ F(X_{T})+\rho
l(u,T)+\frac{\beta}{2}\sum_{j=1}^{N}\sum_{k=1}^{T-1}\|u_{jk}-u_{jk+1}-v_{jk}^{l-1}+\mu_{jk}^{l-1}\|^{2}$
$\displaystyle\textrm{s.t. Constraints
\eqref{eq:model-d-1-cons-h}--\eqref{eq:model-d-1-cons-i}}.$ for
_$k=1,\ldots,T$_ do

for  _$j=1,\ldots,N$_ do

Update variables $v_{jk}$ as
$v^{l}_{jk}=\begin{cases}u^{l}_{jk}-u^{l}_{jk+1}+\mu^{l-1}_{jk}-\alpha/\beta,&\text{if
}u^{l}_{jk}-u^{l}_{jk+1}+\mu^{l-1}_{jk}>\alpha/\beta,\\\
u^{l}_{jk}-u^{l}_{jk+1}+\mu^{l-1}_{jk}+\alpha/\beta,&\text{if
}u^{l}_{jk}-u^{l}_{jk+1}+\mu^{l-1}_{jk}<-\alpha/\beta,\\\
0,&\text{otherwise}.\end{cases}$ Update dual variables $\mu_{jk}$ as
$\displaystyle\mu_{jk}^{l}=\mu_{jk}^{l-1}+(u_{jk}^{l}-u_{jk+1}^{l}-v_{jk}^{l}).$

if
_$\sum_{j=1}^{N}\sum_{k=1}^{T-1}\|u^{l}_{jk}-u^{l}_{jk+1}-v^{l}_{jk}\|^{2}\leq\delta$_
then

Break the loop.

Update $l\leftarrow l+1$.

Output: Final solutions of variables $u^{l},\ v^{l},\ \mu^{l}$.

Algorithm 2 ADMM Algorithm for Solving Continuous Relaxation of (DQCP-TV).

We present the continuous control results obtained from the model (DQCP-TV) in
Figure 5 for the circuit compilation example on the molecule H2. We test the
parameter of the TV regularizer $\alpha=10^{-5},\ 10^{-4},\ 10^{-3},\
10^{-2},\ 10^{-1}$ and the Lagrangian parameter $\beta=0.1,\ 0.5,\ 1.0$. We
choose $\alpha=10^{-3}$ and $\beta=0.5$, which has the smallest objective
value after rounding with min-up-time constraints. The number of switches in
control results decreases significantly as compared with the results in Figure
3, showing the benefits of the TV regularizer.

Figure 5: Control results of the continuous relaxation of the discretized
model with the TV regularizer (DQCP-TV) of the circuit compilation example on
the molecule H2 (objective value: 1.33E$-$05).

##  5 Improvement Heuristic: Approximate Local-Branching Method

Because the ADMM algorithm does not ensure global optimality for nonlinear
objective functions and the objective value increases after rounding, the
control results can be improved by applying an approximate local-branching
(ALB) heuristic. The trust-region method [32] is a widely used local search
method in optimization but highly depends on the initial values. In this
section we propose a trust-region subproblem for the quantum control problem
and then modify the trust-region method to solve the problem starting from the
solutions obtained in Sections 3–4 to improve the quantum controls.

We introduce our improvement heuristic based on the approach of [69] to solve
the binary model with the TV regularizer. Given a feasible point $\hat{u}$, we
use the first-order gradient to approximate the objective value around the
point $\hat{u}$ and define the trust-region subproblem with radius $R$ as

$\displaystyle\min_{u}\quad$ $\displaystyle\langle\nabla
F(\hat{u}),u-\hat{u}\rangle_{L_{2}}+\alpha TV(u)-\alpha TV(\hat{u})$ (29a)
s.t. $\displaystyle\|u-\hat{u}\|_{1}\leq R$ (29b)
$\displaystyle\textrm{Constraints \eqref{eq:model-d-1-cons-u-
sum}--\eqref{eq:model-d-1-cons-u}}.$

We note that (29a) approximates the objective value $F(u)$ but uses an exact
form of the TV regularizer. Constraint (29b) indicates that we consider only
the points with $L_{1}$ distance to $\hat{u}$ no more than $R$. For each
controller $j=1,\ldots,N$ and time step $k=1,\ldots,T-1$, define variables
$v_{jk}$ as the upper bound of the absolute change between the control values
of two consecutive time steps. The trust-region subproblem is reformulated as
the following mixed-integer linear program:

$\displaystyle\min_{u}\quad$ $\displaystyle\langle
F(\hat{u}),u-\hat{u}\rangle_{L_{2}}+\alpha\sum_{k=1}^{T-1}\sum_{j=1}^{N}v_{jk}-\alpha
TV(\hat{u})$ (30a) s.t.
$\displaystyle\sum_{k=1}^{T}\left(\sum_{j:\hat{u}_{jk}=0}u_{jk}+\sum_{j:\hat{u}_{jk}=1}(1-u_{jk})\right)\leq
R$ (30b) $\displaystyle-v_{jk}\leq u_{jk}-u_{jk+1}\leq v_{jk},\ j=1,\ldots,N,\
k=1,\ldots,T-1$ (30c) $\displaystyle\textrm{Constraints
\eqref{eq:model-d-1-cons-u-sum}--\eqref{eq:model-d-1-cons-u}}.$

The objective function (30a) is a reformulation of (29a) with variables $v$.
Constraint (30b) is the trust-region constraint. Constraints (30c) ensure that
for each controller $j=1,\ldots,N$, $v_{jk}\geq|u_{jk}-u_{jk+1}|,\
k=1,\ldots,T-1$. Let $\bar{u}$ be an optimal solution of the model (30).
Define $\Delta F_{p}(\hat{u},\bar{u})$ and $\Delta F_{a}(\hat{u},\bar{u})$
respectively as the predictive and actual decrease of the objective function
with the following formulations:

$\displaystyle\Delta F_{p}(\hat{u},\bar{u})$ $\displaystyle=\langle\nabla
F(\hat{u}),\hat{u}-\bar{u}\rangle_{L_{2}}+\alpha TV(\hat{u}_{jk})-\alpha
TV(\bar{u}_{jk}),$ (31a) $\displaystyle\Delta F_{a}(\hat{u},\bar{u})$
$\displaystyle=F(\hat{u})-F(\bar{u})+\alpha TV(\hat{u}_{jk})-\alpha
TV(\bar{u}_{jk}).$ (31b)

The trust-region algorithm consists of an inner loop and an outer loop. In the
inner loop we solve a sequence of trust-region subproblems a with
monotonically decreasing radius until the ratio between the actual decrease
and predictive decrease is large enough. To obtain a balance between the
computational cost and the searched area, we set a threshold $\bar{R}$ for the
radius. When the radius is greater than the threshold, we decrease it
according to a geometric sequence [69]; otherwise, we decrease it by an
arithmetic sequence. In the outer loop we repeat the inner loop for each
updated point until the predictive decrease is nonpositive. This procedure is
described in Algorithm 3.

Input: Starting radius $R^{0}>0$, threshold of radius $\bar{R}>0$, initial
feasible point $u^{0}$, and threshold of decrease $\eta>0$.

Initialize predictive decrease $\Delta F_{p}(u^{0},\bar{u})=\infty$ and actual
decrease $\Delta F_{a}(u^{0},\bar{u})=-\infty$. Set the number of outer
iterations $l\leftarrow 0$.

while  _$\Delta F_{p}(u^{l},\bar{u}) >0$_ do

Initialize radius $R=R^{0}$.

while  _$\Delta F_{a}(u^{l},\bar{u}) <\eta\Delta F_{p}(u^{l},\bar{u})$_ do

Solve model (30) with $\hat{u}=u^{l}$ to obtain the minimizer $\bar{u}$.

Compute predictive and actual decrease $\Delta F_{p}(u^{l},\bar{u})$ and
$\Delta F_{a}(u^{l},\bar{u})$ by (31).

if  _$R >\bar{R}$_ then

$R\leftarrow\max\left\\{\lfloor R/2\rfloor,\bar{R}\right\\}$.

else

$R\leftarrow R-1$.

Set $l\leftarrow l+1$. Update the central point $u^{l}\leftarrow\bar{u}$.

Output: Control results $u^{l}$.

Algorithm 3 Trust-Region Method for Quantum Control.

######  Remark 7.

If we relax the feasible region of each controller to $[0,1]$, the trust-
region method can solve the continuous relaxation of the model with the TV
regularizer. Trust-region methods show convergence to stationary points on
continuous relaxation if we allow the radius $R$ to be real and adjusted
accordingly (Corollary 3.9 in [70]).

The trust-region approach can also improve solutions obtained by CIA
approaches from Section 3.3.2. Based on the restricted feasible region
$\mathcal{U}_{B}$ instead of the TV regularizer, we propose the trust-region
subproblem with additional linear constraints as follows:

$\displaystyle\min_{u}\quad$ $\displaystyle\langle\nabla
F(\hat{u}),u-\hat{u}\rangle_{L_{2}}$ (32a) s.t. Constraints (3e)–(3f), (30b)
$\displaystyle(u,v)\in\mathcal{U}_{B}.$ (32b)

The objective function (32a) is a first-order gradient approximation of the
original objective function $F(u)$, and constraint (32b) restricts feasible
regions for controls. The computations of the predictive decrease and actual
decrease are modified as

$\displaystyle\Delta F_{p}(\hat{u},\bar{u})$ $\displaystyle=\langle\nabla
F(\hat{u}),\hat{u}-\bar{u}\rangle_{L_{2}},$ (33a) $\displaystyle\Delta
F_{a}(\hat{u},\bar{u})$ $\displaystyle=F(\hat{u})-F(\bar{u}).$ (33b)

For the trust-region algorithm (Algorithm 3), we replace solving model (30)
with solving model (32) and compute the decrease by (33). The other parts
remain the same.

In theory, all time-evolution processes in our algorithms can be conducted on
quantum computers. The inputs for the quantum computers are the control
sequences, and we require the output of the objective function and its
(approximate) gradient. The updates of variables are still computed on
classical computers.

##  6 Numerical Results

We apply our algorithmic framework proposed in Sections 3–5 to diverse
instances of the four examples introduced in Section 2. In Section 6.1 we show
the results from state-of-the-art optimization solvers as the baseline. In
Section 6.2 we introduce the design of numerical instances and parameter
settings. In Sections 6.3–6.5 we present the numerical results, including the
continuous relaxation results, binary results by combinatorial integral
approximation, and binary results after an improvement heuristic.

###  6.1 State-of-the-Art Optimization Solvers

The NEOS server is a frequently used internet-based service containing several
state-of-the-art solvers for numerical optimization problems [71, 72, 73]. For
nonlinear constrained continuous problems, SNOPT [74] and IPOPT [75] are two
widely used solvers. For mixed-integer nonlinear constrained problems, BARON
[76] and Couenne [77] are commonly used for obtaining global optimal
solutions, while MINLP [78, 79], SCIP [80], and Bonmin [81] are effective
solvers with good performance for finding local optima.

As a baseline we first apply the optimization solvers from the NEOS server to
solve the energy minimization example problem in Section 2.1. All the
experiments are conducted on the online server. We set the number of qubits
$q=2$ and $q=6$. We set the evolution time $t_{f}=2$ and the number of time
steps $T=40$. We notice that for all the evolution times in Section 6, we use
dimensionless units, with $\hbar=1$. We apply the solver MUSCOD-II (6.0) [82]
to solve the continuous formulation with the differential equations (1), and
we set the maximum iterations to $100$ for $q=2$ and $700$ for $q=6$. To apply
the state-of-the-art nonlinear optimization solvers to a discretized model, we
derive a discretized formulation of the ordinary differential equation by the
implicit Euler method (see, e.g., [83]) as

$\displaystyle X_{k}-X_{k-1}=-iH_{k}X_{k}\Delta t,\ k=1,\ldots,T.$ (34)

We apply the solvers SNOPT (7.6.1) and IPOPT (3.13.4) to solve the continuous
relaxation model and apply the solvers MINLP, Bonmin (1.8.8), Couenne (0.5.8),
and SCIP (7.0.3.5) to solve the binary model. We set the time limit to $5$
minutes when $q=2$ and $80$ minutes when $q=6$. In Table 2 we present the
objective values, TV-norm values, CPU times, and explored nodes (only for
binary solvers).

Table 2 shows that MUSCOD-II is good for solving the continuous relaxation for
the instance with $q=2$, but it takes a long time to solve the large instance
with $q=6$. SNOPT and IPOPT perform worse than the MUSCOD-II solver especially
when the quantum systems include more qubits. For the binary control problem,
all the methods reach the time limit. Bonmin obtains the best solution, but
the gaps between their results and true energy are all large. For the instance
with $q=6$, some optimization solvers such as MINLP and Couenne run out of
memory. The CPU time increases and the number of explored nodes decreases
significantly when the size of the quantum system increases.

Table 2: Results of solvers on energy minimization example. The results are marked by “OOM” if a solver runs out of memory and “LIMIT” if a solver reaches the time or iteration limit. The explored nodes of continuous solvers are marked by “-” because the node exploration process is conducted only by binary solvers. Solver | $q=2$ | $q=6$
---|---|---
Obj | Time (s) | Nodes | Obj | Time (s) | Nodes
SNOPT | 0.084 | 0.18 | - | 0.736 | 199.64 | -
IPOPT | 0.084 | 0.10 | - | 0.736 | 46.85 | -
MUSCOD-II |  $5.90$E$-09$ | 1.67 | - | 0.154 | 4399.54 | -
MINLP | 0.150 | LIMIT | 9056 | OOM | OOM | OOM
Bonmin | 0.148 | LIMIT | 11672 | 0.793 | LIMIT | 355
Couenne | 0.471 | LIMIT | 6932 | OOM | OOM | OOM
SCIP | 0.149 | LIMIT | 25679 | 0.000 | LIMIT | 2480

This experiment shows that existing standard nonlinear programming and mixed-
integer nonlinear programming solvers cannot solve discretized optimal control
problems. In the following sections we present the numerical results of our
proposed algorithms on multiple instances. Our pGRAPE+SUR method obtains
binary solutions with an objective value $4.22$E$-04$ in $0.14$ seconds for
$q=2$ and binary solutions with an objective value $0.157$ in $27.95$ seconds
for $q=6$. We show that our methods can obtain better results and with shorter
computational time than do the state-of-art solvers.

###  6.2 Experimental Design and Parameter Settings

When the problems have no SOS1 property, the squared $L_{2}$ term in the
penalized model (3.2) is eliminated, so solving (3.2) by pGRAPE is equivalent
to solving the original model (3) by GRAPE. Because the energy minimization
problem has only two controllers, we eliminate the SOS1 property directly by
substituting $u_{2}(t)=1-u_{1}(t)$. We remove the SOS1 property in the CNOT
gate estimation problem to show the generality of our models. We solve the
model with a squared $L_{2}$-penalty function for the circuit compilation
problem.

For all the examples, we apply the pGRAPE algorithm to solve the continuous
relaxation. We also employ the TR and ADMM algorithms to solve the continuous
relaxation with the TV regularizer. Then we obtain the binary results without
constraints by SUR, the binary results with min-up-time constraints, and the
binary results with max-switching constraints by integral minimization.
Furthermore, we apply the approximate local-branching methods to improve the
binary solutions with the TV regularizer and hard control constraints.

For the energy minimization problem, when $q=2$, the matrix $J$ is two-
dimensional and includes only one independent element because it is symmetric.
Hence, we set diagonals in $J$ to $0$ and other elements to $1$. When $q=4,\
6$, we generate $5$ instances where each instance has a random symmetric
matrix $J$ with zero diagonals and elements within a range $[-1,1]$. We will
present averaged objective value results for these problems. The instances are
represented by Energy2, Energy4, and Energy6 corresponding to the number of
qubits $q=2,\ 4,\ 6$ in our following results. For the CNOT problem, we
conduct experiments with evolution times $t_{f}=5,\ 10,\ 15,\ 20$ and
represent the corresponding instances as CNOT5, CNOT10, CNOT15, CNOT20,
respectively. For the NOT gate estimation problem, we follow the parameter
settings in [49] and set $\mu_{1}=0,\ \mu_{2}=2\pi,\ \omega_{1}=1,\
\omega_{2}=\sqrt{2}$. We conduct numerical simulations with evolution times
$t_{f}=2,\ 6,\ 10$ and represent the instances as NOT2, NOT6, and NOT10,
respectively. For the circuit compilation problem, we test two molecules, H2
(dihydrogen) and LiH (lithium hydride), and generate the UCCSD circuits with
minimum energy by the VQE algorithm in the Python package Qiskit [84]. The
instances are represented by CircuitH2 and CircuitLiH. We set the parameters
corresponding to quantum systems as $J_{c}=0.2\pi,\ J_{f}=3\pi,\ J=0.1\pi$. We
test values of the penalty parameter for the squared $L_{2}$-penalty function
$\rho=10^{-6},\ 10^{-5},\ 10^{-4},\ 10^{-3},\ 10^{-2},\ 10^{-1},\ 1,\ 10$, and
we choose the one with the smallest objective value among the parameters with
a squared penalized term less than $10^{-6}$. We set the penalty parameter to
be $\rho=1.0$ and $\rho=0.1$ for two the instances, respectively. We test
values for the TV parameter $\alpha=10^{-5},\ 10^{-4},\ 10^{-3},\ 10^{-2}$ and
choose $\alpha$ with the smallest rounding objective value with min-up-time
constraints. The settings of $\alpha$ and other parameters are presented in
Table 3. All the numerical simulations were conducted in Python 3.8 on a MacOS
computer with 8 cores, 16 GB of RAM, and a 3.20 GHz processor. All the
computational time results are the times on classical computers.

Table 3: Parameter settings of examples. The parameters include the number of qubits ($q$), number of controllers ($N$), evolution time ($t_{f}$), number of time steps ($T$), TV parameter ($\alpha$), minimum uptime steps ($T_{\textrm{minup}}$), and maximum switches ($S$). Instance | $q$ | $N$ | $t_{f}$ | $T$ | $\alpha$ | $T_{\textrm{minup}}$ | $S$
---|---|---|---|---|---|---|---
Energy2 | 2 | 2 | 2 | 40 | 0.01 | 10 | 5
Energy4 | 4 | 2 | 2 | 40 | 0.01 | 10 | 5
Energy6 | 6 | 2 | 2 | 40 | 0.01 | 10 | 5
CNOT5 | 2 | 2 | 5 | 100 | 0.01 | 10 | 20
CNOT10 | 2 | 2 | 10 | 200 | 0.001 | 10 | 20
CNOT15 | 2 | 2 | 15 | 300 | 0.0001 | 10 | 20
CNOT20 | 2 | 2 | 20 | 400 | 0.0001 | 10 | 20
NOT2 | 1 | 2 | 2 | 20 | 0.001 | 5 | 4
NOT6 | 1 | 2 | 6 | 60 | 0.001 | 5 | 12
NOT10 | 1 | 2 | 10 | 100 | 0.001 | 5 | 20
CircuitH2 | 2 | 5 | 4 | 80 | 0.001 | 10 | 8
CircuitLiH | 4 | 12 | 20 | 200 | 0.001 | 5 | 40

In Sections 6.3–6.5, for brevity we select four different instances, Energy6,
CNOT20, NOT10, and CircuitLiH, to analyze the results of objective values and
CPU time. Detailed results of all the methods and instances are presented in
Tables 4, 5, 6, 7, 8, and 9 in Appendix B. Our full code and results are
available on our Github repository [85].

###  6.3 Results of Continuous Relaxations

In Figure 6 we show how the common logarithm of the squared $L_{2}$-norm value
varies with the penalty parameter represented by the lines. We also present an
asymptotic bound $\log_{10}c/\rho$ represented by the dashed lines, where $c$
is a constant marking that a selected data point lies on the dashed lines in
the simulation results. Notably, the figure confirms that
$l(u^{*}_{\rho},T)\sim O(1/\rho)$, as shown in Theorem 3.

(a) Instance CircuitH2

(b) Instance CircuitLiH

Figure 6: Common logarithm of squared $L_{2}$-norm varying with
$\log_{10}\rho$, where $\rho$ is the squared $L_{2}$ penalty parameter. Blue
lines are the common logarithm of the squared $L_{2}$-norm. Orange dashed
lines are functions $\log_{10}c-\log_{10}\rho$, where $c$ is a constant. The
figures show that the squared $L_{2}$-norm $l(u^{*}_{\rho},T)$ decreases as
the penalty parameter $\rho$ increases with a rate $O(1/\rho)$.

We show that pGRAPE always obtains the lowest objective value but the highest
TV regularizer values because it solves the model without the TV regularizer.
ADMM has the best performance for reducing the TV regularizer values, showing
the benefits of introducing a TV regularizer in Section 4. The detailed
results are presented in Appendix B.

We measure the problem size of each instance by $2^{q}\times N\times T$, where
$q$, $N$, and $T$ represent the number of qubits, number of controllers, and
number of time steps. For four groups of instances, Energy, CNOT, NOT, and
Circuit, we present in Figure 7 common logarithm-logarithm figures to show how
the CPU time and the number of iterations vary with the problem size. The CPU
time increases significantly with an increase in the number of qubits because
the dimension of simulation Hamiltonian matrices increases exponentially. The
CPU time also increases when the number of controllers and time steps
increases. Compared with pGRAPE, TR and ADMM spend more time solving the model
because it contains the TV regularizer. The CPU time of ADMM is more than that
of TR for Energy instances, while TR spends more time on CNOT, NOT, and
Circuit instances than ADMM spends.

Figure 7: Common logarithm-logarithm figure of CPU time and iterations of
continuous relaxation for all the instances and number of time steps. Blue,
orange, green, and red colors represent instances Energy, CNOT, NOT, and
Circuit. Circles, triangles, and stars represent obtaining continuous results
by pGRAPE, TR, and ADMM, respectively. Dots with increasing size and
transparency represent results by MT and MS, respectively.

###  6.4 Results of Rounding Techniques

The results for Sections 6.4 and 6.5 are summarized in Figure 12, which shows
the objective values and TV regularizer values of the binary results obtained
by SUR as well as CIA with min-up-time constraints and max-switching
constraints of selected instances. Compared with pGRAPE+MT/MS, TR+MT/MS or
ADMM+MT/MS obtains lower objective values and TV regularizer values for the
binary results with min-up-time constraints and max-switching constraints
because they solve the model with the TV regularizer.

From the results of pGRAPE we show that the squared $L_{2}$-penalty function
ensures that the difference between binary controls of SUR and continuous
results is small. To further demonstrate the convergence of pGRAPE+SUR in
Section 3.3.1, we present how the maximum absolute integral error, which is
the left-hand side of equation (21), varies with the number of time steps in
Figure 8 represented by the lines. The figures show that the error converges
to zero with the increase in the number of time steps, as is claimed in
Proposition 2. We also present the theoretical upper bound of the integral
error, which is the right-hand side of (21), represented by the dashed lines.
We observe that the integral error is always less than the upper bound,
demonstrating the conclusion in Theorem 4. In addition, we present the figures
of how $\epsilon(\Delta t)$ and the upper bound $\sqrt{t_{f}{l(u^{c},T)}\Delta
t}$ vary with the number of time steps for instances CircuitH2 and CircuitLiH
in Figure 9. We show that $\epsilon(\Delta t)$ is always smaller than the
upper bound, demonstrating the conclusion (22) in Corollary 1.

(a) Instance CicruitH2

(b) Instance CircuitLiH

Figure 8: Binary logarithm of maximum absolute integral error and upper bound
of continuous results and SUR binary results. Blue lines marked by circles
represent the maximum absolute integral error, and orange lines marked by
stars represent the upper bound. The error is always smaller than the upper
bound and converges to zero as the number of time steps increases,
demonstrating the conclusion (21).

(a) Instance CicruitH2

(b) Instance CircuitLiH

Figure 9: Binary logarithm of the maximum SOS1 difference of continuous
control $\epsilon(\Delta t)$ and its upper bound in (22). Blue lines represent
$\epsilon(\Delta t)$ and orange dashed lines represent the upper bound. We
show that $\epsilon(\Delta t)$ converges to $0$ with the increase in time
steps and is always smaller than the upper bound, demonstrating Corollary 1.

We set the time limit to $60$ seconds for solving the combinatorial integral
approximation with min-up-time and max-switching constraints. In Figure 10 we
present the CPU time and iterations of all the instances. All the CPU times of
SUR and Energy instances are less than 0.01 seconds so we do not show them in
the figure. The CPU times increase significantly with the number of variables.
For the instance CircuitLiH, all the methods exceed the computational time
limit. In most cases, iterations and CPU times of MS are more than MT because
the problem with MS has a larger integrality gap.

Figure 10: Common logarithm-logarithm figure of CPU time and iterations of CIA
for all the instances. Blue, orange, and green colors represent instances of
CNOT, NOT, and Circuit, respectively. Circles, triangles, and stars represent
obtaining continuous results by pGRAPE, TR, and ADMM, respectively. Dots with
increasing size and transparency represent results by SUR, MT, and MS,
respectively.

###  6.5 Results of Improvement Heuristic

The objective values and TV regularizer values after the improvement heuristic
with the control results of combinatorial integral approximation as initial
points for selected instances are represented by opaque small dots in Figure
12. Because there are multiple local minima for these quantum control
problems, each method may obtain a different performance after the improvement
heuristic. For systems with two controllers, all the methods obtain
performance at a similar level. The performance of different methods varies
more on the circuit compilation problem because it contains more local optima.

We show in Figure 11 the CPU time and iterations for the improvement heuristic
of selected instances. For each group of instances Energy, CNOT, NOT, and
Circuit, the CPU time increases with the problem size. The number of
iterations highly depends on the initial points.

In addition, we observe that for both continuous relaxation and ALB, the CPU
time is exponentially increasing but the number of iterations increases
slightly. We note that the main increase of the computational time comes from
the time-evolution process for simulating the quantum system because the
dimension of the Hamiltonian matrices increases exponentially with the number
of qubits $q$. In our algorithm, we could directly conduct the time evolution
on quantum computers with a control pulse, and then only use the final state
to compute the objective value and approximate the gradient by the finite
difference method. Hence, we do not need to evaluate the state at each time
step, indicating the scalability of our approach for handling large-scale
instances.

Figure 11: Common logarithm-logarithm figure of CPU time and iterations of ALB
for all the instances. Blue, orange, green, and red colors represent groups of
instances Energy, CNOT, NOT, and Circuit. Circles, triangles, and stars
represent obtaining continuous results by pGRAPE, TR, and ADMM, respectively.
Dots with increasing size and transparency represent results by SUR, MT, and
MS, respectively.

In Figure 12 we show the objective values and TV regularizer values of the
combinatorial integral approximation and the improvement heuristic. We show
that the objective value of SUR keeps the same or slightly increases,
demonstrating the theorems and propositions in Section 3.3.1. Adding min-up-
time and max-switching constraints leads to an increase in objective values
because of adding additional constraints stated in Section 3.3.2 to the
feasible set. The increase from adding min-up-time constraints is higher
because the restricted feasible region contains fewer feasible solutions. From
the aspect of the chattering measured by the TV regularizer value, we show
that SUR leads to more chattering compared with continuous solutions because
of the rounding process. As discussed in Section 3.3.2, imposing min-up-time
and max-switching constraints reduces the chattering, and adding min-up-time
constraints reduces it more significantly because the constraints are
stricter. The improvement heuristic reduces the objective values and
chattering remarkably for the results obtained from all the methods with
rounding techniques, especially pGRAPE, which shows the benefits of our
improvement algorithms in Section 5. In Figure 12, the dot closest to the
lower-left corner indicates that it obtains the best balance between improving
objective values and reducing chattering. If the original objective function
matters more, one can choose methods that consist of SUR and the improvement
heuristic. If one focuses more on reducing the switches, the methods with min-
up-time or max-switching constraints and the improvement heuristic are better.

Figure 12: Objective values and TV regularizer values of binary results of
selected instances. Blue, orange, and green dots represent solving continuous
relaxation by pGRAPE, TR, and ADMM, respectively. Circles, triangles, and
stars represent obtaining binary results by SUR, MT, and MS, respectively. Big
dots represent results before ALB, and opaque small dots represent results
after ALB. We annotate the method represented by the closest point to the
lower-left corner.

##  7 Conclusion

In this paper we built a generic control model with both continuous and
discretized formulations for the quantum control problem. We introduced a
penalized squared $L_{2}$ function into the model to ensure that only one
control is active at all times. In addition, we proposed a model with the TV
regularizer aiming to reduce the absolute change of controls.

We developed an algorithmic framework combining the GRAPE approach,
combinatorial integral approximation, and local-branching improvement
heuristic. With numerical simulations on multiple examples with various
objective functions and controllers, we demonstrated the feasibility of the
discrete optimal quantum control problem and illustrated that our algorithms
obtained trade-off controls between high quality and low absolute changes
within a short computational time. Specifically, the performance of different
relaxation models varies among instances, and therefore, testing all three
methods (pGRAPE, TR, and ADMM) is helpful to select the best control. If one
is more interested in obtaining controls with a lower energy or infidelity, we
recommend methods using the sum-up-rounding technique. On the other hand, if
one aims to reduce switches, it is better to choose methods with min-up-time
constraints added. Both cases require the use of the approximate local
branching algorithm to improve final solutions. In practice, the performance
can be improved further by fine-tuning the corresponding parameters.

All the numerical simulations in our paper were conducted on classical
computers , and the CPU time mainly comes from the time-evolution process of
quantum states. Because our algorithm only requires the final state to compute
the objective value and approximate the gradient, running the time-evolution
process in our algorithms on quantum computers is easy in practice and could
help eliminate a significant amount of the exponential slow-down as we observe
in our simulations. Furthermore, with continuous control results, we can
extract appropriate controller sequences and optimize the location of
switching points to obtain high-quality controls. Switching-time optimization
using continuous control results is future research.

## Acknowledgment

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, under contract number DE-AC02-06CH11357. This work was
supported by the U.S. Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research, Accelerated Research for Quantum
Computing program. The authors gratefully acknowledge the support from the
Givens Associate program in Argonne National Laboratory. Xinyu Fei and Siqian
Shen are thankful for the support from U.S. National Science Foundation (NSF)
grants #CMMI-1727618 and Department of Energy (DOE) grant #DE-SC0018018.

Certain commercial equipment, instruments, or materials are identified in this
paper in order to specify the experimental procedure adequately. Such
identification is not intended to imply recommendation or endorsement by NIST,
nor is it intended to imply that the materials or equipment identified are
necessarily the best available for the purpose.

## References

  * Rabitz et al. [2000] Herschel Rabitz, Regina De Vivie-Riedle, Marcus Motzkus, and Karl Kompa.  Whither the future of controlling quantum phenomena?  _Science_ , 288(5467):824–828, 2000.  doi: 10.1126/science.288.5467.824.
  * Werschnik and Gross [2007] J. Werschnik and E. K. U. Gross.  Quantum optimal control theory.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ , 40(18):R175–R211, 2007.  doi: 10.1088/0953-4075/40/18/r01.
  * Brif et al. [2010] Constantin Brif, Raj Chakrabarti, and Herschel Rabitz.  Control of quantum phenomena: Past, present and future.  _New Journal of Physics_ , 12:075008, 2010.  doi: 10.1088/1367-2630/12/7/075008.
  * Shi et al. [1988] Shenghua Shi, Andrea Woody, and Herschel Rabitz.  Optimal control of selective vibrational excitation in harmonic linear chain molecules.  _Journal of Chemical Physics_ , 88(11):6870–6883, 1988.  doi: 10.1063/1.454384.
  * Peirce et al. [1988] Anthony P. Peirce, Mohammed A. Dahleh, and Herschel Rabitz.  Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications.  _Physical Review A_ , 37(12):4950–4964, 1988\.  doi: 10.1103/PhysRevA.37.4950.
  * Shi and Rabitz [1989] Shenghua Shi and Herschel Rabitz.  Selective excitation in harmonic molecular systems by optimally designed fields.  _Chemical Physics_ , 139(1):185–199, 1989.  doi: 10.1016/0301-0104(89)90011-6.
  * Kosloff et al. [1989] R. Kosloff, S. A. Rice, P. Gaspard, S. Tersigni, and D. J. Tannor.  Wavepacket dancing: Achieving chemical selectivity by shaping light pulses.  _Chemical Physics_ , 139(1):201–220, 1989.  doi: 10.1016/0301-0104(89)90012-8.
  * Jakubetz et al. [1990] W. Jakubetz, J. Manz, and H. J. Schreier.  Theory of optimal laser pulses for selective transitions between molecular eigenstates.  _Chemical Physics_ , 165(1):100–106, 1990.  doi: 10.1016/0009-2614(90)87018-M.
  * Khaneja et al. [2005] Navin Khaneja, Timo Reiss, Cindie Kehlet, Thomas Schulte-Herbrüggen, and Steffen J. Glaser.  Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms.  _Journal of Magnetic Resonance_ , 172(2):296–305, 2005.  doi: 10.1016/j.jmr.2004.11.004.
  * Gorshkov et al. [2008] Alexey V. Gorshkov, Tommaso Calarco, Mikhail D. Lukin, and Anders S. Sørensen.  Photon storage in $\Lambda$-type optically dense atomic media, IV: Optimal control using gradient ascent.  _Physical Review A_ , 77:043806, 2008.  doi: 10.1103/physreva.77.043806.
  * van Bijnen and Pohl [2015] R. M. W. van Bijnen and T. Pohl.  Quantum magnetism and topological ordering via Rydberg dressing near Förster resonances.  _Physical Review Letters_ , 114(24):243002, 2015\.  doi: 10.1103/physrevlett.114.243002.
  * Palao and Kosloff [2002] José P. Palao and Ronnie Kosloff.  Quantum computing by an optimal control algorithm for unitary transformations.  _Physical Review Letters_ , 89(18):188301, 2002\.  doi: 10.1103/PhysRevLett.89.188301.
  * Palao and Kosloff [2003] José P. Palao and Ronnie Kosloff.  Optimal control theory for unitary transformations.  _Physical Review A_ , 68(6):062308, 2003.  doi: 10.1103/PhysRevA.68.062308.
  * Montangero et al. [2007] Simone Montangero, Tommaso Calarco, and Rosario Fazio.  Robust optimal quantum gates for Josephson charge qubits.  _Physical Review Letters_ , 99(17):170501, 2007\.  doi: 10.1103/PhysRevLett.99.170501.
  * Grace et al. [2007] Matthew Grace, Constantin Brif, Herschel Rabitz, Ian A. Walmsley, Robert L. Kosut, and Daniel A. Lidar.  Optimal control of quantum gates and suppression of decoherence in a system of interacting two-level particles.  _Journal of Physics B_ , 40(9):S103–S125, 2007\.  doi: 10.1088/0953-4075/40/9/s06.
  * Waldherr et al. [2014] G. Waldherr, Y. Wang, S. Zaiser, M. Jamali, T. Schulte-Herbrüggen, H. Abe, T. Ohshima, J. Isoya, J. F. Du, P. Neumann, and J. Wrachtrup.  Quantum error correction in a solid-state hybrid spin register.  _Nature_ , 506:204–207, 2014.  doi: 10.1038/nature12919.
  * Dolde et al. [2014] Florian Dolde, Ville Bergholm, Ya Wang, Ingmar Jakobi, Boris Naydenov, Sébastien Pezzagna, Jan Meijer, Fedor Jelezko, Philipp Neumann, Thomas Schulte-Herbrüggen, Jacob Biamonte, and Jörg Wrachtrup.  High-fidelity spin entanglement using optimal control.  _Nature Communications_ , 5(3371), 2014.  doi: 10.1038/ncomms4371.
  * Venturelli et al. [2018] Davide Venturelli, Minh Do, Eleanor Rieffel, and Jeremy Frank.  Compiling quantum circuits to realistic hardware architectures using temporal planners.  _Quantum Science and Technology_ , 3(2):025004, 2018.  doi: 10.1088/2058-9565/aaa331.
  * Omran et al. [2019] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin.  Generation and manipulation of Schrödinger cat states in Rydberg atom arrays.  _Science_ , 365(6453):570–574, 2019.  doi: 10.1126/science.aax9743.
  * Khatri et al. [2019] Sumeet Khatri, Ryan LaRose, Alexander Poremba, Lukasz Cincio, Andrew T. Sornborger, and Patrick J. Coles.  Quantum-assisted quantum compiling.  _Quantum_ , 3:140, 2019.  doi: 10.22331/q-2019-05-13-140.
  * Yang et al. [2017] Zhi-Cheng Yang, Armin Rahmani, Alireza Shabani, Hartmut Neven, and Claudio Chamon.  Optimizing variational quantum algorithms using Pontryagin’s minimum principle.  _Physical Review X_ , 7:021027, 2017.  doi: 10.1103/PhysRevX.7.021027.
  * Bapat and Jordan [2019] Aniruddha Bapat and Stephen Jordan.  Bang-bang control as a design principle for classical and quantum optimization algorithms.  _Quantum Information & Computation_, 19:424–446, 2019\.  doi: 10.26421/QIC19.5-6-4.
  * Mbeng et al. [2019] Glen Bigan Mbeng, Rosario Fazio, and Giuseppe Santoro.  Quantum annealing: A journey through digitalization, control, and hybrid quantum variational schemes.  _arXiv:1906.08948_ , 2019.  doi: 10.48550/arXiv.1906.08948.
  * Lin et al. [2019] Chungwei Lin, Yebin Wang, Grigory Kolesov, and Uroš Kalabić.  Application of Pontryagin’s minimum principle to Grover’s quantum search problem.  _Physical Review A_ , 100:022327, 2019.  doi: 10.1103/PhysRevA.100.022327.
  * Brady et al. [2021a] Lucas T Brady, Christopher L Baldwin, Aniruddha Bapat, Yaroslav Kharkov, and Alexey V Gorshkov.  Optimal protocols in quantum annealing and QAOA problems.  _Physical Review Letters_ , 126:070505, 2021a.  doi: 10.1103/PhysRevLett.126.070505.
  * Brady et al. [2021b] Lucas T. Brady, Lucas Kocia, Przemyslaw Bienias, Yaroslav Kharkov Aniruddha Bapat, and Alexey V. Gorshkov.  Behavior of analog quantum algorithms.  _arXiv:2107.01218_ , 2021b.  doi: 10.48550/arXiv.2107.01218.
  * Venuti et al. [2021] Lorenzo Campos Venuti, Domenico D’Alessandro, and Daniel A. Lidar.  Optimal control for quantum optimization of closed and open systems.  _Physical Review Applied_ , 16(5), 2021.  doi: 10.1103/physrevapplied.16.054023.
  * Kadowaki and Nishimori [1998] Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse Ising model.  _Physical Review E_ , 58:5355, 1998.  doi: 10.1103/PhysRevE.58.5355.
  * Farhi et al. [2000] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser.  Quantum computation by adiabatic evolution.  _arXiv:quant-ph/0001106_ , 2000.  doi: 10.48550/arXiv.quant-ph/0001106.
  * Pagano et al. [2020] Guido Pagano, Aniruddha Bapat, Patrick Becker, Katherine S. Collins, Arinjoy De, Paul W. Hess, Harvey B. Kaplan, Antonis Kyprianidis, Wen Lin Tan, Christopher Baldwin, Lucas T. Brady, Abhinav Deshpande, Fangli Liu, Stephen Jordan, Alexey V. Gorshkov, and Christopher Monroe.  Quantum approximate optimization of the long-range Ising model with a trapped-ion quantum simulator.  _PNAS_ , 117(41):25396–25401, 2020.  doi: 10.1073/pnas.2006373117.
  * Harrigan et al. [2021] Matthew P. Harrigan, Kevin J. Sung, Matthew Neeley, Kevin J. Satzinger, Frank Arute, Kunal Arya, Juan Atalaya, Joseph C. Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B. Buckley, David A. Buell, Brian Burkett, Nicholas Bushnell, Yu Chen, Zijun Chen, Ben Chiaro, Roberto Collins, William Courtney, Sean Demura, Andrew Dunsworth, Daniel Eppens, Austin Fowler, Brooks Foxen, Craig Gidney, Marissa Giustina, Rob Graff, Steve Habegger, Alan Ho, Sabrina Hong, Trent Huang, L. B. Ioffe, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Cody Jones, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Seon Kim, Paul V. Klimov, Alexander N. Korotkov, Fedor Kostritsa, David Landhuis, Pavel Laptev, Mike Lindmark, Martin Leib, Orion Martin, John M. Martinis, Jarrod R. McClean, Matt McEwen, Anthony Megrant, Xiao Mi, Masoud Mohseni, Wojciech Mruczkiewicz, Josh Mutus, Ofer Naaman, Charles Neill, Florian Neukart, Murphy Yuezhen Niu, Thomas E. O’Brien, Bryan O’Gorman, Eric Ostby, Andre Petukhov, Harald Putterman, Chris Quintana, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Andrea Skolik, Vadim Smelyanskiy, Doug Strain, Michael Streif, Marco Szalay, Amit Vainsencher, Theodore White, Z. Jamie Yao, Ping Yeh, Adam Zalcman, Leo Zhou, Hartmut Neven, Dave Bacon, Erik Lucero, Edward Farhi, and Ryan Babbush.  Quantum approximate optimization of non-planar graph problems on a planar superconducting processor.  _Nature Physics_ , 17:332–336, 2021.  doi: 10.1038/s41567-020-01105-y.
  * Nocedal and Wright [2006] Jorge Nocedal and Stephen Wright.  _Numerical Optimization_.  Springer Science & Business Media, 2006.  doi: 10.1007/978-0-387-40065-5.
  * Larocca and Wisniacki [2021] Martín Larocca and Diego Wisniacki.  Krylov-subspace approach for the efficient control of quantum many-body dynamics.  _Physical Review A_ , 103(2), 2021.  doi: 10.1103/physreva.103.023107.
  * Doria et al. [2011] Patrick Doria, Tommaso Calarco, and Simone Montangero.  Optimal control technique for many-body quantum dynamics.  _Physical Review Letters_ , 106(19):190501, 2011\.  doi: 10.1103/PhysRevLett.106.190501.
  * Caneva et al. [2011] Tommaso Caneva, Tommaso Calarco, and Simone Montangero.  Chopped random-basis quantum optimization.  _Physical Review A_ , 84(2):022326, 2011.  doi: 10.1103/physreva.84.022326.
  * Sørensen et al. [2018] J. J. W. H. Sørensen, M. O. Aranburu, T. Heinzel, and J. F. Sherson.  Quantum optimal control in a chopped basis: Applications in control of Bose-Einstein condensates.  _Physical Review A_ , 98(2):022119, 2018.  doi: 10.1103/PhysRevA.98.022119.
  * Farhi et al. [2014] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  _arXiv:1411.4028_ , 2014.  doi: 10.48550/arXiv.1411.4028.
  * Bharti et al. [2022] Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik.  Noisy intermediate-scale quantum algorithms.  _Reviews of Modern Physics_ , 94(1), 2022.  doi: 10.1103/revmodphys.94.015004.
  * Cerezo et al. [2021] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles.  Variational quantum algorithms.  _Nature Reviews Physics_ , 3(9):625–644, 2021\.  doi: 10.1038/s42254-021-00348-9.
  * Liang et al. [2020] Daniel Liang, Li Li, and Stefan Leichenauer.  Investigating quantum approximate optimization algorithms under bang-bang protocols.  _Physical Review Research_ , 2(3):033402, 2020\.  doi: 10.1103/physrevresearch.2.033402.
  * Bao et al. [2018] Seraph Bao, Silken Kleer, Ruoyu Wang, and Armin Rahmani.  Optimal control of superconducting gmon qubits using Pontryagin’s minimum principle: Preparing a maximally entangled state with singular bang-bang protocols.  _Physical Review A_ , 97(6):062343, 2018.  doi: 10.1103/physreva.97.062343.
  * Mühlenbein et al. [1988] Heinz Mühlenbein, Martina Gorges-Schleuter, and Ottmar Krämer.  Evolution algorithms in combinatorial optimization.  _Parallel Computing_ , 7(1):65–85, 1988.  doi: 10.1016/0167-8191(88)90098-1.
  * Lawler and Wood [1966] Eugene L Lawler and David E Wood.  Branch-and-bound methods: A survey.  _Operations Research_ , 14(4):699–719, 1966.  doi: 10.1287/opre.14.4.699.
  * Leyffer [2001] Sven Leyffer.  Integrating SQP and branch-and-bound for mixed integer nonlinear programming.  _Computational Optimization and Applications_ , 18(3):295–309, 2001.  doi: 10.1023/A:1011241421041.
  * Vogt and Petersson [2022] Ryan H. Vogt and N. Anders Petersson.  Binary optimal control of single-flux-quantum pulse sequences.  _SIAM Journal on Control and Optimization_ , 60(6):3217–3236, 2022.  doi: 10.1137/21m142808x.
  * Zahedinejad et al. [2014] Ehsan Zahedinejad, Sophie Schirmer, and Barry C Sanders.  Evolutionary algorithms for hard quantum control.  _Physical Review A_ , 90(3):032310, 2014.  doi: 10.1103/PhysRevA.90.032310.
  * Sager et al. [2012] Sebastian Sager, Hans Georg Bock, and Moritz Diehl.  The integer approximation error in mixed-integer optimal control.  _Mathematical Programming_ , 133(1):1–23, 2012\.  doi: 10.1007/s10107-010-0405-3.
  * Pawela and Sadowski [2016] Łukasz Pawela and Przemysław Sadowski.  Various methods of optimizing control pulses for quantum systems with decoherence.  _Quantum Information Processing_ , 15(5):1937–1953, 2016.  doi: 10.1007/s11128-016-1242-y.
  * Motzoi et al. [2009] F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm.  Simple pulses for elimination of leakage in weakly nonlinear qubits.  _Physical Review Letters_ , 103(11), 2009.  doi: 10.1103/physrevlett.103.110501.
  * Bartlett and Musiał [2007] Rodney J. Bartlett and Monika Musiał.  Coupled-cluster theory in quantum chemistry.  _Reviews of Modern Physics_ , 79(1):291, 2007\.  doi: 10.1103/RevModPhys.79.291.
  * Romero et al. [2018] Jonathan Romero, Ryan Babbush, Jarrod R McClean, Cornelius Hempel, Peter J. Love, and Alán Aspuru-Guzik.  Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz.  _Quantum Science and Technology_ , 4(1):014008, 2018.  doi: 10.1088/2058-9565/aad3e4.
  * Chen et al. [2014] Yu Chen, C Neill, P Roushan, N Leung, M Fang, R Barends, J Kelly, B Campbell, Z Chen, B Chiaro, et al.  Qubit architecture with high coherence and fast tunable coupling.  _Physical Review Letters_ , 113(22):220502, 2014\.  doi: 10.1103/PhysRevLett.113.220502.
  * Gokhale et al. [2019] Pranav Gokhale, Yongshan Ding, Thomas Propson, Christopher Winkler, Nelson Leung, Yunong Shi, David I. Schuster, Henry Hoffmann, and Frederic T Chong.  Partial compilation of variational algorithms for noisy intermediate-scale quantum machines.  In _Proceedings of the 52nd Annual IEEE/ACM International Symposium on Microarchitecture_ , pages 266–278, 2019.  doi: 10.1145/3352460.3358313.
  * Jurdjevic and Sussmann [1972] Velimir Jurdjevic and Héctor J Sussmann.  Control systems on Lie groups.  _Journal of Differential Equations_ , 12(2):313–329, 1972.  doi: 10.1016/0022-0396(72)90035-6.
  * Ramakrishna et al. [1995] Viswanath Ramakrishna, Murti V. Salapaka, Mohammed Dahleh, Herschel Rabitz, and Anthony Peirce.  Controllability of molecular systems.  _Physical Review A_ , 51(2):960, 1995.  doi: 10.1103/PhysRevA.51.960.
  * Byrd et al. [1995] Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A limited memory algorithm for bound constrained optimization.  _SIAM Journal on Scientific Computing_ , 16(5):1190–1208, 1995.  doi: 10.1137/0916069.
  * Sinclair [1986] Marius Sinclair.  An exact penalty function approach for nonlinear integer programming problems.  _European Journal of Operational Research_ , 27(1):50–56, 1986.  doi: 10.1016/S0377-2217(86)80006-6.
  * You and Leyffer [2011] Fengqi You and Sven Leyffer.  Mixed-integer dynamic optimization for oil-spill response planning with integration of a dynamic oil weathering model.  _AIChE Journal_ , 57(12):3555–3564, 2011.  doi: 10.1002/aic.12536.
  * Manns and Kirches [2020] Paul Manns and Christian Kirches.  Multidimensional sum-up rounding for elliptic control systems.  _SIAM Journal on Numerical Analysis_ , 58(6):3427–3447, 2020.  doi: 10.1137/19M1260682.
  * Sager [2005] Sebastian Sager.  _Numerical Methods for Mixed-Integer Optimal Control Problems_.  PhD thesis, 2005.
  * Wolsey [2020] Laurence A Wolsey.  _Integer Programming_.  John Wiley & Sons, 2020.  doi: 10.1002/9781119606475.
  * Rudin et al. [1992] Leonid I Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  _Physica D: Nonlinear Phenomena_ , 60(1-4):259–268, 1992.  doi: 10.1016/0167-2789(92)90242-F.
  * Condat [2013] Laurent Condat.  A direct algorithm for 1-D total variation denoising.  _IEEE Signal Processing Letters_ , 20(11):1054–1057, 2013.  doi: 10.1109/LSP.2013.2278339.
  * Kunisch and Hintermüller [2004] Karl Kunisch and Michael Hintermüller.  Total bounded variation regularization as a bilaterally constrained optimization problem.  _SIAM Journal on Applied Mathematics_ , 64(4):1311–1333, 2004.  doi: 10.1137/S0036139903422784.
  * Rodríguez [2013] Paul Rodríguez.  Total variation regularization algorithms for images corrupted with different noise models: A review.  _Journal of Electrical and Computer Engineering_ , 2013, 2013.  doi: 10.1155/2013/217021.
  * Stella et al. [2017] Lorenzo Stella, Andreas Themelis, Pantelis Sopasakis, and Panagiotis Patrinos.  A simple and efficient algorithm for nonlinear model predictive control.  In _56th Annual Conference on Decision and Control_ , pages 1939–1944. IEEE, 2017.  doi: 10.1109/CDC.2017.8263933.
  * Themelis et al. [2018] Andreas Themelis, Lorenzo Stella, and Panagiotis Patrinos.  Forward-backward envelope for the sum of two nonconvex functions: Further properties and nonmonotone linesearch algorithms.  _SIAM Journal on Optimization_ , 28(3):2274–2303, 2018.  doi: 10.1137/16M1080240.
  * Sager and Zeile [2021] Sebastian Sager and Clemens Zeile.  On mixed-integer optimal control with constrained total variation of the integer control.  _Computational Optimization and Applications_ , 78(2):575–623, 2021.  doi: 10.1007/s10589-020-00244-5.
  * Leyffer and Manns [2021] Sven Leyffer and Paul Manns.  Sequential linear integer programming for integer optimal control with total variation regularization.  _arXiv:2106.13453_ , 2021.  doi: 10.48550/arXiv.2106.13453.
  * Aravkin et al. [2022] Aleksandr Y. Aravkin, Robert Baraldi, and Dominique Orban.  A proximal quasi-Newton trust-region method for nonsmooth regularized optimization.  _SIAM Journal on Optimization_ , 32(2):900–929, 2022.  doi: 10.1137/21m1409536.
  * Czyzyk et al. [1998] Joseph Czyzyk, Michael P. Mesnier, and Jorge J. Moré.  The NEOS server.  _IEEE Journal on Computational Science and Engineering_ , 5(3):68–75, 1998.  doi: 10.1109/99.714603.
  * Dolan [2001] Elizabeth D. Dolan.  The NEOS server 4.0 administrative guide.  Technical Memorandum ANL/MCS-TM-250, Mathematics and Computer Science Division, Argonne National Laboratory, 2001.
  * Gropp and Moré [1997] William Gropp and Jorge J. Moré.  Optimization environments and the NEOS server.  In Martin D. Buhman and Arieh Iserles, editors, _Approximation Theory and Optimization_ , pages 167–182. Cambridge University Press, 1997.
  * Andrei [2017] Neculai Andrei.  A SQP algorithm for large-scale constrained optimization: SNOPT.  In _Continuous nonlinear optimization for engineering applications in GAMS technology_ , pages 317–330. Springer, 2017.  doi: 10.1007/978-3-319-58356-3.
  * Wächter and Biegler [2006] Andreas Wächter and Lorenz T Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  _Mathematical Programming_ , 106(1):25–57, 2006\.  doi: 10.1007/s10107-004-0559-y.
  * Sahinidis [1996] Nikolaos V Sahinidis.  BARON: A general purpose global optimization software package.  _Journal of Global Optimization_ , 8(2):201–205, 1996.  doi: 10.1007/bf00138693.
  * Belotti [2009] Pietro Belotti.  Couenne: A user’s manual.  Technical report, Lehigh University, 2009.
  * Belotti et al. [2013] Pietro Belotti, Christian Kirches, Sven Leyffer, Jeff Linderoth, James Luedtke, and Ashutosh Mahajan.  Mixed-integer nonlinear optimization.  _Acta Numerica_ , 22:1–131, 2013.  doi: 10.1017/S0962492913000032.
  * Leyffer and Mahajan [2011] Sven Leyffer and Ashutosh Mahajan.  Software for nonlinearly constrained optimization.  In James J. Cochran, Louis A. Cox, Pinar Keskinocak, Jeffrey P. Kharoufeh, and J. Cole Smith, editors, _Wiley Encyclopedia of Operations Research and Management Science_. John Wiley & Sons, Inc., 2011.  doi: 10.1002/9780470400531.eorms0570.
  * Gamrath et al. [2020] Gerald Gamrath, Daniel Anderson, Ksenia Bestuzheva, Wei-Kun Chen, Leon Eifler, Maxime Gasse, Patrick Gemander, Ambros Gleixner, Leona Gottwald, Katrin Halbig, Gregor Hendel, Christopher Hojny, Thorsten Koch, Pierre Le Bodic, Stephen J. Maher, Frederic Matter, Matthias Miltenberger, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Christine Tawfik, Stefan Vigerske, Fabian Wegscheider, Dieter Weninger, and Jakob Witzig.  The SCIP Optimization Suite 7.0.  ZIB-Report 20-10, Zuse Institute Berlin, 2020.
  * Bonami et al. [2008] Pierre Bonami, Lorenz T. Biegler, Andrew R. Conn, Gérard Cornuéjols, Ignacio E. Grossmann, Carl D. Laird, Jon Lee, Andrea Lodi, François Margot, Nicolas Sawaya, and Andreas Wächter.  An algorithmic framework for convex mixed integer nonlinear programs.  _Discrete optimization_ , 5(2):186–204, 2008\.  doi: 10.1016/j.disopt.2006.10.011.
  * Kirches and Leyffer [2013] Christian Kirches and Sven Leyffer.  TACO: A toolkit for AMPL control optimization.  _Mathematical Programming Computation_ , 5(3):227–265, 2013.  doi: 10.1007/s12532-013-0054-7.
  * Butcher [2016] John Charles Butcher.  _Numerical Methods for Ordinary Differential Equations_.  John Wiley & Sons, 2016.  doi: 10.1002/9781119121534.
  * Aleksandrowicz et al. [2021] Gadi Aleksandrowicz, Thomas Alexander, Panagiotis Barkoutsos, Luciano Bello, Yael Ben-Haim, David Bucher, Francisco Jose Cabrera-Hernández, Jorge Carballo-Franquis, Adrian Chen, Chun-Fu Chen, et al.  Qiskit: An open-source framework for quantum computing.  2021\.  doi: 10.5281/ZENODO.2562110.
  * Fei [2022] Xinyu Fei.  Code and Results: Binary Control Pulse Optimization for Quantum Systems.  https://github.com/xinyufei/Quantum-Control-qutip, 2022.
  * Rebentrost and Wilhelm [2009] Patrick Rebentrost and Frank K Wilhelm.  Optimal control of a leaking qubit.  _Physical Review B_ , 79(6):060507, 2009.  doi: 10.1103/physrevb.79.060507.

##  Appendix A Proofs of Results in Section 3

In this section, we provide proofs of the results in Section 3.

###  A.1 Proofs of Results in Section 3.2

###### Proof of Theorem 3.

Because $F$ is upper bounded by constant $C_{F}$, we have $|F(u)|\leq C_{F}$
for all the $u$. We prove that for any $\rho$ and $T$, $|\rho
l(u^{*}_{\rho},T)|\leq 2C_{F}$. Assume that there exists an optimal solution
$u^{(1)}$ of (3.2) such that $|\rho l(u^{(1)},T)|>2C_{F}$, then we have
$F(u^{(1)})+\rho l(u^{(1)},T)>C_{F}$. This is a contradiction to the
definition of $u^{(1)}$ because we can find a feasible solution $u^{(2)}$ of
(3) such that $F(u^{(2)})+\rho l_{(}u^{(2)},T)=F(u^{(2)})\leq C_{F}$.
Therefore we drive that $|l(u^{*}_{\rho},T)|\leq 2C_{F}/\rho$, which means
that $l(u^{*}_{\rho},T)\sim O(1/\rho)$. ∎

###  A.2 Proofs of Results in Section 3.3

###### Proof of Proposition 1.

Because the SOS1 property holds for binary control $u^{b}$, we have

$\displaystyle\max_{k=1,\ldots,T}\left\|\sum_{\tau=1}^{k}(u^{c}_{\tau}-u^{b}_{\tau})\Delta
t\right\|_{\infty}$
$\displaystyle=\max_{k=1,\ldots,T}\max_{j=1,\ldots,N}\left|\sum_{\tau=1}^{k}(u^{c}_{j\tau}-u^{b}_{j\tau})\Delta
t\right|$
$\displaystyle\geq\max_{k=1,\ldots,T}\frac{1}{N}\left|\sum_{j=1}^{N}\sum_{\tau=1}^{k}(u^{c}_{j\tau}-u^{b}_{j\tau})\Delta
t\right|$
$\displaystyle=\max_{k=1,\ldots,T}\frac{1}{N}\left|\sum_{\tau=1}^{k}\left(\sum_{j=1}^{N}u^{c}_{j\tau}-1\right)\Delta
t\right|=\frac{1}{N}\epsilon(\Delta t).$ (35)

∎

###### Proof of Theorem 4.

For simplicity, we use $\epsilon$ to denote $\epsilon(\Delta t)$. It is
obvious that for any time step $k=1,\ldots,T$, exactly one control can be set
to $1$ by $u^{b}$ by the construction. Hence $u^{b}$ satisfies the SOS1
property. For a specific time step $k=1,\ldots,T$, define

$\displaystyle i_{k}=\mathop{\rm
arg\,max}_{j=1,\ldots,N}\left|\sum_{\tau=1}^{k}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{k}u^{b}_{j\tau}\Delta t\right|.$ (36)

We assume that there exists a time step $r$ such that the claim does not hold.
For simplicity, we use $i$ to represent the corresponding $i_{r}$. We consider
the assumption in two cases.

##### Case 1:

We assume that,

$\displaystyle\sum_{\tau=1}^{r}u^{c}_{i\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{i\tau}\Delta t<-\left(N-1\right)\Delta
t-\frac{2N-1}{N}\epsilon.$ (37)

Let $\hat{k}$ be the highest index of the time step in which control $i$ is
rounded up,

$\displaystyle\hat{k}=\mathop{\rm arg\,max}_{1\leq k\leq r}\left\\{k\ |\
u^{b}_{ik}=1\right\\}.$ (38)

Then by (37), we have,

$\displaystyle\sum_{\tau=1}^{\hat{k}}u^{c}_{i\tau}\Delta
t+\left(N-1\right)\Delta t+\frac{2N-1}{N}\epsilon\leq$
$\displaystyle\sum_{\tau=1}^{r}u^{c}_{i\tau}\Delta t+\left(N-1\right)\Delta
t+\frac{2N-1}{N}\epsilon$ (39a) $\displaystyle<$
$\displaystyle\sum_{\tau=1}^{r}u^{b}_{i\tau}\Delta
t=\sum_{\tau=1}^{\hat{k}}u^{b}_{i\tau}\Delta t.$ (39b)

Because $u^{b}_{i\hat{k}}=1$, we know that $i$ has the maximum value of (39a)
among $j=1,\ldots,N$. Hence it follows from (39),

$\displaystyle\sum_{\tau=1}^{\hat{k}}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{\hat{k}}u^{b}_{j\tau}\Delta t<-\left(N-1\right)\Delta
t-\frac{2N-1}{N}\epsilon,\ j=1,\ldots,N.$ (40)

Summing up over all the controls $j$,

$\displaystyle-\epsilon\leq\sum_{j=1}^{N}\left(\sum_{\tau=1}^{\hat{k}}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{\hat{k}}u^{b}_{j\tau}\Delta t\right)<-N\left(N-1\right)\Delta
t-(2N-2)\epsilon.$ (41)

The first inequality comes from the definition of $\epsilon$. This leads to
$0<-N(N-1)\Delta t-(N-1)\epsilon$, which is a contradiction because
$\epsilon\geq 0$ and $N\geq 2$.

##### Case 2:

We assume that,

$\displaystyle\sum_{\tau=1}^{r}u^{c}_{i\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{i\tau}\Delta t>\left(N-1\right)\Delta
t+\frac{2N-1}{N}\epsilon.$ (42)

From the definition of $\epsilon$, we have

$\displaystyle\sum_{1=j\neq i}^{N}\left(\sum_{\tau=1}^{r}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{j\tau}\Delta
t\right)+\sum_{\tau=1}^{r}u^{c}_{i\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{i\tau}\Delta t\leq\epsilon$ (43)

and by substituting (42) into (43), it holds that

$\displaystyle\sum_{1=j\neq i}^{N}\left(\sum_{\tau=1}^{r}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{j\tau}\Delta t\right)+(N-1)\Delta
t+\frac{N-1}{N}\epsilon<0.$ (44)

Because the left-hand side can be written as the sum of $N-1$ terms as

$\displaystyle\Delta
t+\frac{1}{N}\epsilon+\sum_{\tau=1}^{r}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{j\tau}\Delta t,$ (45)

at least one of them must be negative, therefore there exists $\hat{j}$ such
that

$\displaystyle\Delta
t+\frac{1}{N}\epsilon+\sum_{\tau=1}^{r}u^{c}_{\hat{j}\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{\hat{j}\tau}\Delta t<0.$ (46)

Let $\hat{k}$ be the highest index of the time step in which control $\hat{j}$
is rounded up,

$\displaystyle\hat{k}=\mathop{\rm arg\,max}_{1\leq k\leq r}\left\\{k\ |\
u^{b}_{\hat{j}k}=1\right\\}.$ (47)

Then we have

$\displaystyle\sum_{\tau=1}^{\hat{k}}u^{c}_{\hat{j}\tau}\Delta
t-\sum_{\tau=1}^{\hat{k}-1}u^{b}_{\hat{j}\tau}\Delta t\leq\Delta
t+\sum_{\tau=1}^{r}u^{c}_{\hat{j}\tau}\Delta
t-\sum_{\tau=1}^{r}u^{b}_{\hat{j}\tau}\Delta t<-\frac{1}{N}\epsilon.$ (48)

The first inequality follows from $u_{\hat{j}\hat{k}}^{b}=1$ and
$u_{\hat{j}k}^{b}=0,\ k\geq\hat{k}$. The second inequality follows from (46).
Because $\hat{j}$ is the control which is rounded up at time step $\hat{k}$,
for any $j=1,\ldots,N$, it holds that,

$\displaystyle\sum_{\tau=1}^{\hat{k}}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{\hat{k}}u^{b}_{j\tau}\Delta
t\leq\sum_{\tau=1}^{\hat{k}}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{\hat{k}-1}u^{b}_{j\tau}\Delta t<-\frac{1}{N}\epsilon.$ (49)

Summing over all the controls, we have

$\displaystyle\sum_{j=1}^{N}\sum_{\tau=1}^{\hat{k}}\left(u^{c}_{j\tau}-u^{b}_{j\tau}\right)\Delta
t<-\epsilon.$ (50)

This contradicts the definition of the parameter $\epsilon$. ∎

###### Proof of Corollary 1.

Based on the formulation of the discretized control in the continuous
relaxation of the model (3.2), we have

$\displaystyle\epsilon(\Delta t)$
$\displaystyle=\max_{k=1,\ldots,T}\left|\sum_{\tau=1}^{k}\left(\sum_{j=1}^{N}u^{c}_{j\tau}-1\right)\Delta
t\right|\leq\sum_{k=1}^{T}\left|\sum_{j=1}^{N}u^{c}_{jk}-1\right|\Delta t$
$\displaystyle\leq\sqrt{T}\sqrt{\sum_{k=1}^{T}\left(\sum_{j=1}^{N}u^{c}_{jk}-1\right)^{2}}\Delta
t=\sqrt{Tl(u^{c},T)}\Delta t=\sqrt{t_{f}{l(u^{c},T)}\Delta t}.$ (51)

From Theorem 4, we directly obtain the statement (22). From Theorem 3, if the
original function $F$ is bounded, the optimized squared $L_{2}$ term
$l(u^{c},T)$ is uniformly bounded over time steps $T$, then the absolute
integral error between continuous and discretized control converges with
$O(\sqrt{\Delta t})$. ∎

###### Proof of Proposition 2.

For any time step $k=1,\ldots,T$, we have

$\displaystyle\lim_{\Delta t\rightarrow 0}X^{c}_{k}$
$\displaystyle=\lim_{\Delta t\rightarrow 0}\prod_{\tau=1}^{k}e^{-i\Delta
tH_{\tau}}X_{0}=\lim_{\Delta t\rightarrow 0}\prod_{\tau=1}^{k}e^{-i\Delta
t\left(H^{(0)}+\sum_{j=1}^{N}u^{c}_{j\tau}H^{(j)}\right)}X_{0}$
$\displaystyle=\lim_{\Delta t\rightarrow 0}e^{-i\Delta tH^{(0)}-i\Delta
t\sum_{j=1}^{N}\sum_{\tau=1}^{k}u^{c}_{j\tau}H^{(j)}}X_{0}.$ (52)

The last equality follows by Trotter expansion that $\displaystyle e^{\Delta
t(A+B)}=e^{\Delta tA}e^{\Delta tB}+O(\Delta t^{2})$ and taking the limit as
$\Delta t$ goes to zero. Similarly, for binary control, we have

$\displaystyle\lim_{\Delta t\rightarrow 0}X^{b}_{k}=\lim_{\Delta t\rightarrow
0}e^{-i\Delta tH^{(0)}-i\Delta
t\sum_{j=1}^{N}\sum_{\tau=1}^{k}u^{b}_{j\tau}H^{(j)}}X_{0}.$ (53)

From Theorem 4, we have for any time step $k=1,\ldots,T$ and controller
$j=1,\ldots,N$,

$\displaystyle\left|\sum_{\tau=1}^{k}u^{c}_{j\tau}\Delta
t-\sum_{\tau=1}^{k}u^{b}_{j\tau}\Delta
t\right|\leq\left\|\sum_{\tau=1}^{k}\left(u^{c}_{\tau}-u^{b}_{\tau}\right)\Delta
t\right\|_{\infty}\leq\left(N-1\right)\Delta t+\frac{2N-1}{N}\epsilon(\Delta
t).$ (54)

Combining with Corollary 1 and taking the limit as $\Delta t$ goes to zero, we
have

$\displaystyle\lim_{\Delta t\rightarrow 0}\sum_{\tau=1}^{k}u^{c}_{j\tau}\Delta
t=\lim_{\Delta t\rightarrow 0}\sum_{\tau=1}^{k}u^{b}_{j\tau}\Delta t.$ (55)

We substitute (55) into the formulation of states (A.2) and (53), then we
obtain the conclusion that

$\displaystyle\lim_{\Delta t\rightarrow 0}X^{b}_{k}=\lim_{\Delta t\rightarrow
0}X^{c}_{k},\ k=1,\ldots,T.$ (56)

Substituting states into the objective function, we prove that

$\displaystyle\lim_{\Delta t\rightarrow 0}F(X_{T}^{b})=\lim_{\Delta
t\rightarrow 0}F(X_{T}^{c}).$ (57)

∎

##  Appendix B Detailed Numerical Results

In this section, we present the results of continuous relaxation,
combinatorial integral approximation, and the improvement heuristic for all
the methods and instances in Table 4–9. We also present the objective values
and TV regularizer values for all the instances and annotate the best method
in Figure 13.

Table 4: Objective value results of continuous relaxation.  Instance | Objective value | TV regularizer
---|---|---
pGRAPE | TR | ADMM | pGRAPE | TR | ADMM
Energy2 | 1.10E$-12$ | 9.11E$-05$ | 8.94E$-05$ | 0.999 | 0.567 | 0.523
Energy4 | 0.154 | 0.155 | 0.168 | 5.064 | 4.114 | 2.752
Energy6 | $0.213$ | $0.213$ | $0.219$ | 5.999 | 4.508 | 3.237
CNOT5 | 0.169 | 0.124 | 0.191 | 16.000 | 9.419 | 6.094
CNOT10 | 1.16E$-09$ | 2.65E$-04$ | 3.21E$-04$ | 20.979 | 15.194 | 11.056
CNOT15 | 1.00E$-10$ | 8.43E$-06$ | 3.65E$-06$ | 29.846 | 24.348 | 16.795
CNOT20 | 5.93E$-10$ | 4.06E$-06$ | 8.07E$-07$ | 26.162 | 23.481 | 15.099
NOT2 | 0.163 | 0.163 | 0.163 | 2.226 | 1.188 | 0.908
NOT6 | 4.28E$-$10 | 6.55E$-$05 | 1.07E$-$04 | 3.251 | 0.652 | 0.578
NOT10 | 6.55E$-$11 | 5.14E$-$05 | 4.79E$-$05 | 2.740 | 1.530 | 1.005
CircuitH2 | 4.37E$-07$ | 1.24E$-04$ | 1.33E$-05$ | 18.720 | 8.421 | 2.744
CircuitLiH | 1.14E$-03$ | 1.45E$-03$ | 1.43E$-03$ | 53.976 | 48.720 | 0.677
Table 5: CPU time and iterations of continuous relaxation.  Instance | CPU time (s) | Iterations
---|---|---
pGRAPE | TR | ADMM | pGRAPE | TR | ADMM
Energy2 | 0.13 | 1.77 | 19.56 | 7 | 132 | 100
Energy4 | 2.89 | 27.99 | 163.08 | 44 | 914 | 100
Energy6 | 27.94 | 341.43 | 1195.41 | 51 | 1578 | 100
CNOT5 | 1.12 | 79.55 | 21.75 | 111 | 3146 | 100
CNOT10 | 0.75 | 192.78 | 70.21 | 31 | 3279 | 100
CNOT15 | 1.26 | 348.96 | 100.59 | 39 | 3900 | 100
CNOT20 | 1.02 | 432.73 | 150.25 | 21 | 3929 | 100
NOT2 | 0.05 | 4.93 | 2.63 | 15 | 1268 | 100
NOT6 | 0.15 | 48.81 | 10.93 | 24 | 3563 | 100
NOT10 | 0.10 | 91.39 | 34.05 | 6 | 3762 | 100
CircuitH2 | 3.01 | 154.05 | 39.86 | 244 | 3832 | 100
CircuitLiH | 663.08 | 1403.03 | 564.95 | 4345 | 4247 | 100

Table 6: Objective value results of combinatorial integral approximation.

Instance Objective value/TV regularizer pGRAPE+SUR TR+SUR ADMM+SUR pGRAPE+MT
TR+MT ADMM+MT pGRAPE+MS TR+MS ADMM+MS Energy2 4.22E$-04$/54 4.91E$-03$/54
4.01E$-04$/48 $0.159$/4 $0.159$/6 $0.154$/6 $0.029$/10 $0.040$/10 $0.028$/10
Energy4 $0.158$/26.8 $0.158$/32.4 $0.170$/44.8 $0.367$/6 $0.317$/6 $0.363$/6
$0.163$/10 $0.162$/10 $0.195$/10 Energy6 $0.216$/38.8 $0.216$/43.6
$0.221$/52.8 $0.684$/6 $0.616$/6 $0.606$/6 $0.247$/10 $0.258$/10 $0.240$/10
CNOT5 0.170/16 0.332/24 0.190/41 0.243/10 0.525/6 0.285/7 0.170/16 0.593/15
0.191/32 CNOT10 6.01E$-04$/116 1.78E$-03$/116 1.68E$-03$/86 0.158/22 0.323/21
0.084/15 0.011/39 0.019/38 0.006/32 CNOT15 1.12E$-03$/266 2.30E$-03$/276
2.90E$-03$/279 0.539/37 0.290/36 0.176/27 0.325/38 0.284/40 0.214/39 CNOT20
1.45E$-03$/491 8.30E$-04$/480 1.46E$-03$/467 0.782/53 0.314/51 0.517/47
0.654/39 0.697/39 0.619/39 NOT2 0.164/3 0.164/1 0.164/1 0.164/1 0.164/1
0.164/1 0.164/3 0.164/1 0.164/1 NOT6 2.38E$-$03/40 3.10E$-$03/46 1.55E$-$03/52
4.05E$-$02/8 3.65E$-$02/9 7.56E$-$02/10 1.38E$-$02/19 1.58E$-$01/14
1.54E$-$01/14 NOT10 4.73E$-$03/126 3.91E$-$04/126 9.22E$-$03/21 2.78E$-$02/24
1.34E$-$01/25 6.38E$-$02/6 6.16E$-$02/40 2.84E$-$01/40 9.22E$-$03/21 CircuitH2
0.027/32 0.027/36 0.006/76 0.600/8 0.591/8 0.063/8 0.026/22 0.038/24 0.008/22
CircuitLiH 0.168/380 0.182/380 0.033/252 0.963/68 0.966/72 0.658/48 0.287/290
0.224/288 0.557/148

Table 7: CPU time and iterations of combinatorial integral approximation.

Instance | CPU time (s)/Iterations
---|---
pGRAPE+MT | TR+MT | ADMM+MT | pGRAPE+MS | TR+MS | ADMM+MS
Energy2 |  $<0.01$/39 |  $<0.01$/39 |  $<0.01$/18 |  $<0.01$/3715 |  $<0.01$/4627 |  $<0.01$/5106
Energy4 |  $<0.01$/15 |  $<0.01$/17 |  $<0.01$/17 |  $<0.01$/759 |  $<0.01$/2578 |  $<0.01$/925
Energy6 |  $<0.01$/8 |  $<0.01$/8 |  $<0.01$/8 |  $<0.01$/2473 |  $<0.01$/3730 |  $<0.01$/2508
CNOT5 | 0.55/5 | 0.55/3 | 0.57/2 | 0.50/2 | 0.52/2 | 0.72/7
CNOT10 | 19.29/36786 | 13.50/5916 | 6.49/4136 | 61.83/25685 | 61.83/27432 | 61.83/54921
CNOT15 | 7.92/555 | 8.89/10 | 25.79/21326 | 64.36/14349 | 64.62/9878 | 64.70/29421
CNOT20 | 28.46/31236 | 22.86/8482 | 38.42/1837 | 67.41/2069 | 67.22/1242 | 67.30/26378
NOT2 | 0.05/1 | 0.05/1 | 0.04/1 | 0.04/1 | 0.06/1 | 0.07/1
NOT6 | 0.80/7403 | 6.62/32640 | 0.26/1 | 41.37/47882 | 0.26/1 | 60.19/90179
NOT10 | 0.97/1255 | 60.51/1815 | 4.26/4136 | 60.50/29462 | 0.66/1 | 0.59/1
CircuitH2 | 0.89/1 | 0.87/1 | 3.83/1514 | 1.34/1 | 1.87/1 | 60.76/24639
CircuitLiH | 70.94/349 | 70.87/889 | 71.01/1570 | 71.45/1545 | 70.83/1595 | 71.23/2137

Table 8: Objective value results of improvement heuristic.

Instance Objective value/TV regularizer Improvement of objective
value/Improvement of TV regularizer pGRAPE+ TR+ ADMM+ pGRAPE+ TR+ ADMM+
pGRAPE+ TR+ ADMM+ SUR+ALB SUR+ALB SUR+ALB MT+ALB MT+ALB MT+ALB MS+ALB MS+ALB
MS+ALB Energy2 $0.003$/10 $0.005$/8 $0.003$/4 $0.003$/4 $0.003$/4 $0.041$/6
$0.001$/10 $0.003$/10 $0.002$/8 $-0.20$%/81.48% 0.00%/85.19% $-0.20$%/91.67%
18.55%/0.00% 18.55%/0.00% 13.36%/0.00% 2.88%/0.00% 3.85%/0.00% 2.67%/20.00%
Energy4 $0.160$/17.2 $0.162$/14 $0.166$/14.4 $0.199$/6 $0.198$/6 $0.218$/5.6
$0.160$/9.2 $0.160$/10 $0.184$/8.4 $-0.14$%/35.82% $-0.42$%/56.79%
$0.53$%/67.86% 26.61%/0.00% 17.34%/0.00% 22.84%/6.67% 0.38%/8.00% 0.32%/16.00%
1.31%/0.00% Energy6 $0.218$/32.4 $0.216$/40.4 $0.219$/50.0 $0.429$/6 $0.429$/6
$0.451$/6 $0.229$/10 $0.230$/10 $0.216$/10 $-0.25$%/16.49% $-0.01$%/7.34%
0.25%/5.30% 80.74%/0.00% 49.09%/0.00% 39.28%/0.00% 2.47%/0.00% 3.77%/0.00%
0.57%/0.00% CNOT5 0.176/9 0.206/7 0.176/9 0.195/9 0.196/10 0.195/9 0.170/16
0.173/20 0.172/16 $-3.53$%/43.75% 37.95%/70.83% 7.37%/78.05% 19.75%/10.00%
62.67%/$-66.67$% 31.58%/$-28.57$% 0.00%/0.00% 70.83%/$-33.33$% 9.95%/$50.00$%
CNOT10 1.58E$-03$/30 2.46E$-03$/24 1.15E$-03$/20 4.06E$-03$/23 9.43E$-03$/16
6.04E$-03$/15 9.80E$-04$/39 1.31E$-03$/38 1.18E$-03$/36 $-66.42$%/74.14%
$-38.20$%/70.69% 0.00%/76.74% 97.47%/$-4.55$% 97.21%/23.81% 92.86%/0.00%
90.91%/0.00% 94.74%/0.00% 83.33%/$-12.50$% CNOT15 5.59E$-04$/262
4.67E$-04$/256 8.51E$-04$/263 6.31E$-03$/33 6.25E$-03$/34 1.63E$-03$/30
1.30E$-03$/39 3.72E$-03$/39 1.91E$-03$/40 44.10%/1.50% 76.65%/7.25%
71.64%/5.73% 98.89%/10.81% 97.93%/5.56% 99.43%/$-11.11$% 99.69%/$-2.63$%
98.94%/2.50% 99.07%/$-2.56$% CNOT20 4.56E$-04$/479 6.13E$-04$/471
5.07E$-04$/441 1.20E$-03$/38 8.22E$-04$/49 1.35E$-03$/48 9.47E$-04$/40
2.81E$-03$/40 7.45E$-04$/40 54.44%/2.44% 38.68%/1.88% 49.27%/5.57%
99.87%/28.30% 99.74%/3.92% 99.81%/$-2.13$% 99.86%/$-2.56$% 99.57%/$-2.56$%
99.88%/$-2.56$% NOT2 0.164/1 0.164/1 0.164/1 0.164/1 0.164/1 0.164/1 0.163/3
0.163/3 0.163/3 0.46%/66.67% 0.36%/0.00% 0.36%/0.00% 0.36%/0.00% 0.36%/0.00%
0.36%/0.00% 0.62%/0.00% 0.52%/$-$200.00% 0.52%/$-$200.00% NOT6 1.42E$-$03/6
7.33E$-$04/10 7.24E$-$04/13 1.45E$-$03/9 1.27E$-$02/11 3.39E$-$03/12
1.06E$-$03/22 8.90E$-$04/18 1.74E$-$03/16 40.24%/85.00% 76.34%/86.96% 53.25%
/75.00% 96.41%/$-$12.50% 65.21%/$-$22.22% 95.52%/$-$20.00% 92.33%/$-$15.79%
99.44%/$-$28.57% 98.87% /$-$14.29% NOT10 5.50E$-$04/32 1.66E$-$03/6
6.16E$-$04/7 7.37E$-$04/24 9.09E$-$04/20 8.90E$-$04/8 5.44E$-$04/38
7.57E$-$04/39 8.90E$-$04/21 88.37%/74.60% $-$76.45%/92.06% 93.32%/66.67%
97.35%/0.00% 99.32%/20.00% 98.60%/$-$33.33% 99.12%/5.00% 99.73%/2.50%
90.35%/0.00% CircuitH2 0.003/24 0.013/32 0.006/76 0.245/12 0.007/2 0.054/10
0.014/18 0.003/22 0.008/22 88.89%/25.00% 51.85%/11.11% 0.00%/0.00%
59.17%/$-50.00$% 98.82%/75.00% 14.29%/$-25.00$% 46.15%/18.18% 92.11%/8.33%
0.00%/0.00% CircuitLiH 0.069/378 0.083/378 0.033/252 0.002/6 0.496/70 0.355/48
0.165/286 0.120/290 0.021/158 58.93%/0.53% 54.40%/0.53% 0.00%/0.00%
99.79%/91.18% 48.65%/2.78% 46.05%/0.00% 42.51%/1.38% 46.43%/$-0.69$%
96.23%/$-6.76$%

Table 9: CPU time and iterations of improvement heuristic.

Instance CPU time (s)/Iterations pGRAPE+ TR+ ADMM+ pGRAPE+ TR+ ADMM+ pGRAPE+
TR+ ADMM+ SUR+ALB SUR+ALB SUR+ALB MT+ALB MT+ALB MT+ALB MS+ALB MS+ALB MS+ALB
Energy2 3.41/255 5.21/379 2.26/262 2.26/170 2.25/170 1.59/150 1.27/93 1.90/154
1.24/91 Energy4 3.36/127 5.53/191 7.90/280 5.12/207 5.13/207 4.94/194 1.21/43
1.23/43 2.14/68 Energy6 16.95/94 11.11/62 9.88/56 41.42/219 39.53/208
33.36/174 18.51/100 24.02/144 10.84/56 CNOT5 12.69/394 16.96/506 30.09/993
19.67/521 20.34/550 21.46/575 0.91/1377 33.99/33 12.58/464 CNOT10 70.36/1245
56.97/1031 56.90/998 74.83/945 84.95/1050 51.37/574 19.89/301 24.04/323
16.65/270 CNOT15 9.12/106 30.48/395 25.57/317 150.25/1238 169.47/1335
107.17/895 113.02/944 129.21/1288 82.44/935 CNOT20 27.82/249 26.22/251
33.77/397 208.79/1246 139.69/1010 157.57/1009 153.63/1035 155.07/1138
157.87/1182 NOT2 0.46/93 0.40/63 0.52/93 0.39/63 0.46/62 0.31/63 0.34/63
0.47/63 0.31/62 NOT6 11.71/725 7.36/306 2.49/188 9.46/592 5.47/231 8.63/582
7.82/458 5.82/236 6.90/470 NOT10 29.69/1181 16.16/409 7.87/298 21.46/824
21.89/558 14.11/541 7.22/256 17.05/465 4.06/161 CircuitH2 4.84/97 3.18/65
4.70/33 5.40/88 9.21/156 2.73/39 4.88/108 3.14/65 1.68/33 CircuitLiH 26.08/108
46.26/108 9.05/36 64.29/186 74.88/153 59.21/135 41.43/143 20.02/71 39.78/158

Figure 13: Objective values and TV regularizer values of binary results of all
the instances.

##  Appendix C Control Results of NOT Estimation Problem

For the NOT gate estimation problem, we present the control obtaining the best
trade-off between objective values and TV regularizer values which are
annotated in Figure 14. All three instances show that controller 1 has a more
significant impact on the final infidelity. When the evolution time is short
($t_{f}=2$), the optimal control sets controller 1 active all the time. When
the evolution time is longer ($t_{f}=6,\ 10$), the active time of controller 2
is shorter. The main reason is that only considering controller 1 can also
obtain low infidelity with a sufficiently long evolution time but the
performance of a single controller is worse than two controllers with the same
evolution time [86, 49].

(a) $t_{f}=2$

(b) $t_{f}=6$

(c) $t_{f}=10$

Figure 14: Control results of NOT gate estimation example with different
evolution times. Blue lines represent controller 1 and orange dashed lines
represent controller 2.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator
of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of
Energy Office of Science laboratory, is operated under Contract No. DE-
AC02-06CH11357. The U.S. Government retains for itself, and others acting on
its behalf, a paid-up nonexclusive, irrevocable worldwide license in said
article to reproduce, prepare derivative works, distribute copies to the
public, and perform publicly and display publicly, by or on behalf of the
Government. The Department of Energy will provide public access to these
results of federally sponsored research in accordance with the DOE Public
Access Plan http://energy.gov/downloads/doe-public-access-plan.