Optimal control of large quantum systems: assessing memory and runtime performance of GRAPE

The rapidly evolving field of quantum information processing has generated much interest in quantum optimal control for tasks such as state preparation [1, 2], error correction [3, 4] and realization of logical gates [5–8]. This methodology has been implemented in various quantum systems [9], such as nuclear magnetic resonance [10–15], trapped ions [16, 17], nitrogen-vacancy centers in diamonds [18–24], Bose–Einstein condensation [25–27] and neutral atoms [28, 29].

One of the prominent techniques for implementing quantum optimal control is Gradient Ascent Pulse Engineering (GRAPE). The central goal of optimal control is to adjust control parameters in such a way that the infidelity of the desired state transfer or gate operation is minimized. In GRAPE, the minimization is based on gradient descent and thus generally requires the evaluation of gradients. Automatic differentiation (AD) [30] is a convenient tool that has also made an impact on quantum optimal control with GRAPE [31–34]. Internally, AD builds a computational graph and applies the chain rule to automatically compute the gradient of a given function.

However, the convenience of AD comes at the cost of memory required for storing the full computational graph. The use of semi-automatic differentiation (semi-AD) [35] partially reduces the memory overhead compared to full automatic differentiation (full-AD). This reduction can still be insufficient to tackle optimal control for quantum systems of rapidly growing size [36, 37]. This motivates us to revisit the original GRAPE scheme [12] which is based on hard-coded gradients (HG), and has the smallest possible memory cost. We perform a scaling analysis of memory usage and runtime, and thereby develop a decision tree guiding the optimal choice of strategy among HG, semi-AD and full-AD. We exemplify the scaling behavior with benchmarks for concrete optimization tasks.

The outline of our paper is as follows. Section 2 briefly reviews GRAPE and introduces necessary notation. In section 3, we illustrate derivation and memory-efficient implementation of the gradients for the most commonly used cost function contributions. Our numerical implementation of HG further improves the efficiency of state propagation. In section 4, we analyze and compare the scaling of runtime and memory usage scaling for HG, AD and semi-AD, and present results from concrete benchmark studies. In section 5, we formulate a decision tree that facilitates the choice of optimal numerical strategy for GRAPE-based quantum optimal control. We summarize our findings and present our conclusions in section 6.

We briefly sketch the basics of GRAPE [12], mainly to establish the notion used in subsequent sections.

Consider a system described by the Hamiltonian

$$\begin{eqnarray}&&H(t)={H}_{s}+a(t){h}_{c}.\end{eqnarray} \tag{ 1 }$$

Here, H_s denotes the static system Hamiltonian and h_c is an operator coupling the classical control field a(t) to the system ¹ . Quantum optimal control aims to adjust a(t) such that a desired unitary gate or state transfer is performed within a given time interval 0 ≤ t ≤ T. To facilitate optimization, the total control time T is divided into N intervals of duration Δt = T/N, and the control field is specified by the amplitudes at the discrete times t_n = nΔt (n = 1, 2,...,N). The discretization interval Δt is chosen such that control amplitudes are approximately constant within each Δt. This treatment is actually quite close to modeling the actual drive output of an arbitrary waveform generator, where Δt can be chosen to align with the available resolution. We denote the time-discretized values by

$$\begin{eqnarray}&&{a}_{n}:=a({t}_{n}),\end{eqnarray} \tag{ 2 }$$

and group these to form the vector $a\in {{\mathbb{R}}}^{N}$. Under the influence of this piecewise-constant control field, the system's quantum state ∣ψ_n 〉:=∣ψ(t_n )〉 evolves by a single time step according to

$$\begin{eqnarray}&&| {\psi }_{n}\rangle ={U}_{n}| {\psi }_{n-1}\rangle .\end{eqnarray} \tag{ 3 }$$

Here, U_n is the short-time propagator

$$\begin{eqnarray}&&{U}_{n}:=U({t}_{n},{t}_{n-1})={e}^{-{{iH}}_{n}{\rm{\Delta }}t}\qquad ({\hslash }=1),\end{eqnarray} \tag{ 4 }$$

where H_n = H_s + a_n h_c . Optimization of the control field proceeds via minimization of a cost function C. Typically, C includes the state-transfer or gate infidelity, along with additional cost contributions employed to moderate pulse characteristics such as maximal power or bandwidth. To this end a composite cost function is constructed, $C(a)={\sum }_{\nu }{\alpha }_{\nu }{C}_{\nu }\left(a\right)$, where α_ν are empirically chosen weight factors and C_ν are individual cost contributions. GRAPE utilizes the method of steepest descent to minimize C(a) by updating the control amplitudes according to the cost function gradient:

$$\begin{eqnarray}&&{{\bf{a}}}^{(j+1)}={{\bf{a}}}^{(j)}-{\eta }_{j}\,{{\rm{\nabla }}}_{{\bf{a}}}C({{\bf{a}}}^{(j)}).\end{eqnarray} \tag{ 5 }$$

Here, j enumerates steepest-descent iterations and η_j > 0 is the learning rate governing the step size.

A critical ingredient for GRAPE is the computation of cost function gradients. One popular option is to leverage automatic differentiation for this purpose, such as implemented in Tensorflow [38] and Autograd [39]. A clear advantage of this strategy is the ease by which complicated cost function gradients are evaluated automatically at runtime. This convenience, however, is bought at the cost of significant memory consumption which can be prohibitive for large quantum systems. This waste of memory resources is prohibitive when optimal control is applied to large quantum systems. Therefore, we are motivated to revisit hard-coded gradients in a computationally efficient scheme for key cost contributions.

3.1. Analytical gradients

We provide a discussion of the analytical gradients of two key cost contributions in this subsection. Expressions for gradients of other cost contributions are shown in appendix A.

In state-transfer problems, a central goal of quantum optimal control is to minimize the state-transfer infidelity given by

$$\begin{eqnarray}&&{C}_{1}^{{st}}=1-{\left|\left\langle {\phi }_{T}| {\psi }_{N}\right\rangle \right|}^{2},\end{eqnarray} \tag{ 6 }$$

where ∣ϕ_T 〉 is the target state and ∣ψ_N 〉 is a state realized at final time t_N . The gradient of ${C}_{1}^{\mathrm{st}}$ takes the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {C}_{1}^{{st}}}{\partial {a}_{n}}=-\langle {\phi }_{T}| {U}_{N}\cdots \displaystyle \frac{\partial {U}_{n}}{\partial {a}_{n}}\cdots {U}_{1}| {\psi }_{0}\rangle \langle {\psi }_{N}| {\phi }_{T}\rangle -{\rm{c}}.{\rm{c}}.=-\langle {\phi }_{n}| \displaystyle \frac{\partial {{\rm{U}}}_{n}}{\partial {{\rm{a}}}_{n}}| {\psi }_{n-1}\rangle \langle {\psi }_{N}| {\phi }_{T}\rangle -{\rm{c}}.{\rm{c}}.,\end{eqnarray} \tag{ 7 }$$

where the state ∣ϕ_n 〉 obeys the recursive relation

$$\begin{eqnarray}&&| {\phi }_{n}\rangle ={U}_{n+1}^{\dagger }| {\phi }_{n+1}\rangle .\end{eqnarray} \tag{ 8 }$$

This can be interpreted as backward propagation.

In addition to state transfer, unitary gates are another operation of interest. Gate optimization can be achieved by minimizing the gate infidelity ${C}_{1}^{g}=1-{\left|\mathrm{tr}\left({U}_{T}^{\dagger }{U}_{R}\right)/d\right|}^{2}$. Here, U_T is the target unitary gate and U_R = U_N ⋯ U₁ is the actually realized gate. For a given orthonormal basis $\{| {\psi }_{0}^{h}\rangle \}{}_{h}$ (h enumerates the basis states), the gradient of ${C}_{1}^{g}$ can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {C}_{1}^{g}}{\partial {a}_{n}}=-\displaystyle \frac{1}{{d}^{2}}\displaystyle \sum _{h}\langle {\phi }_{n}^{h}| \displaystyle \frac{\partial {U}_{n}}{\partial {a}_{n}}| {\psi }_{n-1}^{h}\rangle \displaystyle \sum _{h^{\prime} }\langle {\psi }_{N}^{h^{\prime} }| {\phi }_{N}^{h^{\prime} }\rangle -{\rm{c}}.{\rm{c}}.,\end{eqnarray} \tag{ 9 }$$

where $| {\psi }_{n}^{h}\rangle ={U}_{n}\cdots {U}_{1}| {\psi }_{0}^{h}\rangle$. The state $| {\phi }_{n}^{h}\rangle$ is obtained as follows: apply the target unitary U_T to basis state $\{| {\psi }_{0}^{h}\rangle \}{}_{h}$ and backward propagate the result to time t_n ,

$$\begin{eqnarray}&&| {\phi }_{n}^{h}\rangle ={U}_{n+1}^{\dagger }\cdots {U}_{N}^{\dagger }{U}_{T}| {\psi }_{0}^{h}\rangle .\end{eqnarray} \tag{ 10 }$$

Equations (7) and (9) provide the analytical expressions needed for gradient descent. In the next two subsections, we will discuss how to numerically evaluate these expressions in a memory-efficient way.

3.2. Numerical implementation of hard-coded gradient

The calculation of gradients in equations (7) and (9) requires numerical evaluation of expressions of the form U_n ∣Ψ〉 and $\tfrac{\partial {U}_{n}}{\partial {a}_{n}}| {\rm{\Psi }}\rangle$. Here, the short-time propagator U_n is given by a matrix exponential e^A and A = − iH_n Δt.

Numerically evaluating e^A ∣Ψ〉 is not without challenge. As pointed out by Moler and Van Loan [40], several methods exist, yet none of them is truly optimal. Three methods ranked highly in reference [40] are based on solving ordinary differential equations, matrix diagonalization and scaling and squaring (combined with series expansion). ODE solving, in this context, tends to be most expensive in runtime [35]. We mainly on focus on scaling and squaring but comment on situations where matrix diagonalization is preferable in section 5.

Scaling and squaring [41] which is based on the identity

$$\begin{eqnarray}&&{e}^{A}={\left({e}^{A/s}\right)}^{s},\end{eqnarray} \tag{ 11 }$$

where the right hand side matrix exponential now involves the matrix B = A/s which has a norm that is reduced compared to $\left|\left|A\right|\right|$. This generally allows for a lower truncation order m of the exponential series ² ,

$$\begin{eqnarray}&&{e}^{B}\approx {e}_{m}^{B}:= \displaystyle \sum _{q=0}^{m}\displaystyle \frac{{B}^{q}}{q!}.\end{eqnarray} \tag{ 12 }$$

Typically, m is smaller than the one required for the original e^A . Our goal is to approximate the state vector ${e}^{A}| {\rm{\Psi }}\rangle \approx {\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle$ where

$$\begin{eqnarray}&&{e}_{m}^{B}| {\rm{\Psi }}\rangle =\left(1+B+\displaystyle \frac{{B}^{2}}{2!}+\cdots +\displaystyle \frac{{B}^{m}}{m!}\right)| {\rm{\Psi }}\rangle =| {\rm{\Psi }}\rangle +B| {\rm{\Psi }}\rangle +\displaystyle \frac{B\left(B| {\rm{\Psi }}\rangle \right)}{2}+\cdots +\displaystyle \frac{B\left(\cdots \left(B| {\rm{\Psi }}\rangle \right)\right)}{m!}.\end{eqnarray} \tag{ 13 }$$

According to the last expression, approximation of e^A ∣Ψ〉 can thus proceed by repeated application of B = A/s to a vector, without ever invoking matrix-matrix multiplication. This boosts runtime performance from O(d³) to O(d²).

To proceed with the evaluation of e^A ∣Ψ〉, the truncation order m and scaling factor s are determined in such a way that the error remains below the error tolerance τ,

$$\begin{eqnarray}&&{\varepsilon }_{A,| {\rm{\Psi }}\rangle }(m,s):= \displaystyle \frac{\left|\left|{e}^{A}| {\rm{\Psi }}\rangle -[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle ]\kern-2pt] \right|\right|}{\left|\left|{e}^{A}| {\rm{\Psi }}\rangle \right|\right|}\lt \tau .\end{eqnarray} \tag{ 14 }$$

Here, 〚•〛 denotes the floating point representation of •. Since this error ε_A,∣Ψ〉 is difficult and costly to evaluate, we instead derive an upper bound E_A (m, s) [see appendix C] and resort to the inequality

$$\begin{eqnarray}&&{\varepsilon }_{A,| {\rm{\Psi }}\rangle }(m,s)\lt {E}_{A}(m,s)\leqslant \tau .\end{eqnarray} \tag{ 15 }$$

For a given A, this inequality generally has multiple solutions for m and s. As shown in appendix C, the evaluation requires ms matrix-vector multiplications, which are observed to dominate the runtime. Since systematic minimization of ms requires significant runtime itself, we instead: (1) select solutions to equation (15) with the smallest s which is denoted by ${\mathbb{s}}$, and (2) among these pairs, choose the one with smallest m (which we denote by ${\mathbb{m}}$).

Thus, calculating e^A ∣Ψ〉 requires

$$\begin{eqnarray}&&\mu := {\mathbb{ms}}\end{eqnarray} \tag{ 16 }$$

matrix-vector multiplications. As shown in appendix D for two concrete examples, our method for evaluating e^A ∣Ψ〉 can achieve better runtime performance and accuracy than Scipy's expm multiply function [42] .

Having addressed the challenges associated with numerically evaluating e^A ∣Ψ〉, we now turn to the calculation of $\tfrac{\partial {U}_{n}}{\partial {a}_{n}}| {\rm{\Psi }}\rangle$, which is crucial for computing cost function gradients as in equations (7) and (9). We base this evaluation on the approximation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {U}_{n}}{\partial {a}_{n}}| {\rm{\Psi }}\rangle \approx \displaystyle \frac{\partial {\left({e}_{m}^{{B}_{n}}\right)}^{s}}{\partial {a}_{n}}| {\rm{\Psi }}\rangle ,\end{eqnarray} \tag{ 17 }$$

where B_n = − iH_n Δt/s. We proceed by using the auxiliary matrix method based on the identity [43, 44]

$$\begin{eqnarray}&&{\left({e}_{m}^{{{ \mathcal A }}_{n}/s}\right)}^{s}\left(\begin{array}{c}0\\ | {\rm{\Psi }}\rangle \end{array}\right)=\left(\begin{array}{c}\displaystyle \frac{\partial {\left({e}_{m}^{{B}_{n}}\right)}^{s}}{\partial {a}_{n}}| {\rm{\Psi }}\rangle \\ {\left({e}_{m}^{{B}_{n}}\right)}^{s}| {\rm{\Psi }}\rangle .\end{array}\right)\end{eqnarray} \tag{ 18 }$$

Here, ${{ \mathcal A }}_{n}$ is defined as

$$\begin{eqnarray}{{ \mathcal A }}_{n}:= \left(\begin{array}{cc}-{{iH}}_{n}{\rm{\Delta }}t & -{{ih}}_{c}{\rm{\Delta }}t\\ 0 & -{{iH}}_{n}{\rm{\Delta }}t\end{array}\right),\end{eqnarray} \tag{ 19 }$$

which is a 2 × 2 matrix of operators each acting on our Hilbert space with dimension d. Once a basis is chosen, ${{ \mathcal A }}_{n}$ can be represented as a 2d × 2d matrix of complex numbers. Numerical evaluation of the left hand side of the identity (18) gives a 2d dimensional vector. The first d entries of the vector correspond to $\partial {\left({e}_{m}^{{B}_{n}}\right)}^{s}/\partial {a}_{n}| {\rm{\Psi }}\rangle$, which is an approximation of $\tfrac{\partial {U}_{n}}{\partial {a}_{n}}| {\rm{\Psi }}\rangle$.

3.3. Efficient evaluation of hard-coded gradients

We consider the following three methods to compute gradients within GRAPE: evaluation of hard-coded gradients, full automatic differentiation [32], semi-automatic differentiation [35]. To determine the method most appropriate for a given problem, we inspect the memory usage and runtime scaling of these three approaches. Following this analysis, we present benchmark studies of specific state-transfer and gate-operation problems to illustrate the scaling behavior.

4.1. Scaling analysis of runtime and memory usage

4.2. Benchmark studies

Our first benchmark example studies the following state-transfer problem. Consider a setup in which a driven flux-tunable transmon is capacitively coupled to a superconducting cavity, and perform a state transfer from the cavity vacuum state to the 20-photon Fock state. The Hamiltonian in the frame co-rotating with the transmon drive is

$$\begin{eqnarray}&&{H}_{1}(t)={\rm{\Delta }}{b}^{\dagger }b+\displaystyle \frac{1}{2}\alpha {b}^{\dagger }b({b}^{\dagger }b-1)+g({{cb}}^{\dagger }+{c}^{\dagger }b)+{a}_{z}(t){b}^{\dagger }b+{a}_{x}(t)(b+{b}^{\dagger }).\end{eqnarray} \tag{ 20 }$$

Here, the drive is both longitudinal and transverse, c and b are the annihilation operators for cavity photons and transmon excitations respectively. α, g and Δ denote anharmonicity, interaction strength and difference between the transmon and cavity frequency.

Optimization proceeds by minimizing the infidelity ${C}_{1}^{{st}}$ and limiting the occupation of higher transmon levels by including the cost contribution ${C}_{2}^{{st}}$ [equation (A.1)]. The latter has the additional benefit of enabling the truncation to only a few transmon levels. We measure runtime and peak memory usage of full-AD ⁴ , semi-AD and HG for a single optimization iteration. The runtime is measured by the package Temci [45] and the peak memory usage is monitored by the memory profiler package [46]. The optimization is performed on an AMD Ryzen Threadripper 3970X 32-Core CPU with 3.7 GHz frequency.

Benchmark results for the state-transfer optimization are shown in figure 2. We consider a single time step (N = 1) and measure the runtime versus d. The results in figure 2(a) confirm that full-AD, semi-AD and HG have the same scaling of runtime per time step with respect to d (table 1). The observed runtime scaling ${\rm{\Theta }}({d}^{\tfrac{5}{2}})$ is consistent with κ ∼ Θ(d²) due to dense matrix storage for H₁, and $\mu \sim {\rm{\Theta }}({d}^{\tfrac{1}{2}})$ (see appendix E). The scaling of the required memory resources is illustrated in figure 2(b): memory usage scales with dimension d as Θ(d²) for full-AD, semi-AD and HG, consistent with κ ∼ Θ(d²) being the dominant contribution. As anticipated, the scaling of memory with N differs between full-AD, semi-AD and HG [figure 2(c)]. For full-AD and semi-AD, the scaling is linear with N, while it is independent of N for HG. The slope of full-AD is much larger than the one of semi-AD, since the former has extra dependence on μ (table 1), which is a constant ≈100 in this benchmark. The results confirm the favorable memory efficiency of HG and semi-AD compared to full-AD in state-transfer problems with large number of time steps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next turn to an example for optimizing a unitary gate. This benchmark is performed for three driven transmons with nearest neighbor coupling. The target gate consists of simultaneous Hadamard operations on each of the three transmons. For concreteness, we consider the case of transmons with the same frequency, driven resonantly. The full Hamiltonian in the rotating frame of the drive is

$$\begin{eqnarray}&&{H}_{2}(t)=\displaystyle \sum _{\nu =1}^{3}\left[\displaystyle \frac{1}{2}\alpha {b}_{\nu }^{\dagger }{b}_{\nu }({b}_{\nu }^{\dagger }{b}_{\nu }-1)+{a}^{(\nu )}(t)({b}_{\nu }+{b}_{\nu }^{\dagger })\right]+g({b}_{1}{b}_{2}^{\dagger }+{b}_{2}{b}_{3}^{\dagger }+{\rm{h}}.{\rm{c}}.).\end{eqnarray} \tag{ 21 }$$

The goal of the optimization is to minimize the infidelity ${C}_{2}^{g}$. The setup for this benchmark is the same as in the previous example.

Figure 3 records the benchmarks based on a single time step (N = 1). The observed runtime scaling with dimension d is ${\rm{\Theta }}({d}^{\tfrac{11}{3}})$ for full-AD, semi-AD and HG, see figure 3(a). This power law arises from Θ(μ κ d) for κ ∼ Θ(d²) and $\mu \sim {\rm{\Theta }}({d}^{\tfrac{2}{3}})$ [see appendix E]. The memory usage scaling with d is illustrated in figure 3(b), confirming that full-AD has an unfavorable extra $\mu \sim {\rm{\Theta }}({d}^{\tfrac{2}{3}})$ dependence compared to HG and semi-AD (table 1). The memory usage scaling with N is as expected [figure 3(c)]: for semi-AD, the scaling is linear with N; for HG, it is independent of N. (The memory usage of full-AD in this case exceeds our resources, and thus is not shown.) Throughout our benchmark, this linear dependence with N causes semi-AD to consume about 10 times more memory than HG. This underlines that HG is preferable whenever memory constraints become a limitation in optimization problems with large Hilbert space dimension and large number of time steps. We emphasize that this improvement is achieved without worsening the runtime scaling.

Zoom In Zoom Out Reset image size

There are several considerations which determine whether full-AD, HG or semi-AD is the optimal strategy for a given optimization problem, see the decision tree figure 4. Whenever the Hilbert space dimension d and the number of time steps N are sufficiently small, then memory usage may not be a bottleneck. In this case, full-AD may be preferable, since it eliminates the need for hard-coded gradients. If the memory cost of full-AD is not affordable but analytical gradients are tedious to compute, semi-AD may be a viable compromise, though still less memory efficient than HG. In all other cases, HG is the method of choice, since it can handle optimization with larger N and d compared to full-AD and semi-AD.

Zoom In Zoom Out Reset image size

The choice of implementing HG via scaling and squaring or matrix diagonalization should be informed by whether the Hamiltonian is sparse and κ scales better than Θ(d²). In the latter case, scaling and squaring is preferable and leads to an overall memory usage ∼Θ(d + κ), which is better than the memory scaling O(d²) of matrix diagonalization. By contrast, for Hamiltonians with κ ∼ Θ(d²), these two methods are equivalent in memory usage scaling and we are free to select the one with better runtime for a specific optimization task. While the runtime for optimization using matrix diagonalization scales as O(d³) [appendix F], the runtime for optimization based on scaling and squaring differs between state-transfer and gate optimizations: (1) for state transfer, the runtime scales as Θ(μ d²); hence, whenever μ scales better than Θ(d), scaling and squaring wins over matrix diagonalization. (2) for gate optimization, the runtime scales as Θ(μ d³); hence matrix diagonalization is preferable.

Motivated by the memory bottleneck problem in quantum optimal control of large systems, we have revisited and compared the memory usage and runtime scaling of multiple flavors of GRAPE: hard-coded gradients (HG), full automatic differentiation (full-AD) and semi-automatic differentiation (semi-AD). These three methods are equivalent in runtime scaling but differ characteristically in scaling of memory usage. The ranking among them is as follows, going from the most to least efficient: HG, semi-AD and full-AD. Thus HG is the preferred option when facing a memory bottleneck problem.

We have illustrated the scaling of HG and full-AD by benchmark studies for concrete state-transfer and gate-operation optimizations, respectively. As part of this, we have presented improvements in the implementation of scaling and squaring, used to evaluate propagator-state products numerically. In the example studied, we observed a speedup compared to Scipy's expm multiply function. The decision tree shown in figure 4 summarizes the choice of appropriate numerical strategy for a given optimization problem.

We thank Thomas Propson for fruitful discussions. This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Superconducting Quantum Materials and Systems Center (SQMS) under contract number DE-AC02-07CH11359. We are grateful for the following open-source packages used in our work: matplotlib [47], numpy [48], scipy [49], sympy [50], scqubits [51], and qoc [31].

All data that support the findings of this study are included within the article (and any supplementary files).

In the main text, we exclusively focused on minimizing state-transfer or gate infidelity. In this appendix, we provide the analytical gradients for other cost contributions of interest.

The cost contribution

$$\begin{eqnarray}&&{C}_{2}^{{st}}=\displaystyle \frac{1}{N}\displaystyle \sum _{n^{\prime} =1}^{N}\langle {\psi }_{n^{\prime} }| {\rm{\Omega }}| {\psi }_{n^{\prime} }\rangle \end{eqnarray} \tag{ A.1 }$$

penalizes the integrated expectation value of the Hermitian operator Ω. The gradient of this cost contribution can be written as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {C}_{2}^{{st}}}{\partial {a}_{n}} & = & \displaystyle \frac{1}{N}\displaystyle \sum _{n^{\prime} =1}^{N}\langle {\psi }_{n^{\prime} }| {\rm{\Omega }}\displaystyle \frac{\partial }{\partial {a}_{n}}\displaystyle \prod _{n^{\prime\prime} =1}^{n^{\prime} }{U}_{n^{\prime\prime} }| {\psi }_{0}\rangle +c.c.=\displaystyle \frac{1}{N}\left(\displaystyle \sum _{n^{\prime} =n}^{N}\displaystyle \prod _{n^{\prime\prime} =n+1}^{n^{\prime} }\langle {\psi }_{n^{\prime} }| {\rm{\Omega }}{U}_{n^{\prime\prime} }\displaystyle \frac{\partial }{\partial {a}_{n}}{U}_{n}| {\psi }_{n-1}\rangle \right)+c.c.\\ & = & \displaystyle \frac{2}{N}\mathrm{Re}(\langle {{\rm{\Phi }}}_{n}| \displaystyle \frac{\partial }{\partial {{\rm{a}}}_{n}}{{\rm{U}}}_{n}| {\psi }_{n-1}\rangle ),\end{array}\end{eqnarray} \tag{ A.2 }$$

with

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{n}\rangle := \displaystyle \sum _{n^{\prime} =n}^{N}\,{U}_{n+1}^{\dagger }\cdots {U}_{n^{\prime} }^{\dagger }{\rm{\Omega }}| {\psi }_{n^{\prime} }\rangle \end{eqnarray} \tag{ A.3 }$$

The vector ∣Φ_n 〉 obeys the recursion relation

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{n}\rangle ={U}_{n+1}^{\dagger }| {{\rm{\Phi }}}_{n+1}\rangle +{\rm{\Omega }}| {\psi }_{n}\rangle .\end{eqnarray} \tag{ A.4 }$$

Due to this relation, the expression in equation (A.2) can be evaluated with a forward-backward scheme analogous to the one discussed in the section 3.3.

The cost contribution ${C}_{3}^{{st}}=1-\tfrac{1}{N}{\sum }_{n^{\prime} =1}^{N}{\left|\left\langle {\phi }_{T}| {\psi }_{n^{\prime} }\right\rangle \right|}^{2}$ sums up infidelities from all time steps (as opposed to recording the final state infidelity only, as in ${C}_{1}^{{st}}$). This steers the system towards its target state more rapidly, thus helping to reduce the duration of state-transfer. Similar to equation (A.2), the gradient of ${C}_{3}^{{st}}$ is:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {C}_{3}^{{st}}}{\partial {a}_{n}}=-\displaystyle \frac{2}{N}\mathrm{Re}\left(\left\langle {{\rm{\Phi }}}_{T,n}| \displaystyle \frac{\partial }{\partial {{\rm{a}}}_{n}}{{\rm{U}}}_{n}| {\psi }_{n-1}\right\rangle \right),\end{eqnarray} \tag{ A.5 }$$

with

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{T,n}\rangle := \displaystyle \sum _{n^{\prime} =n}^{N}{U}_{n+1}^{\dagger }\cdots {U}_{n^{\prime} }^{\dagger }| {\phi }_{T}\rangle \langle {\phi }_{T}| {\psi }_{n^{\prime} }\rangle .\end{eqnarray} \tag{ A.6 }$$

The vector obeys the recursive relation

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{T,n}\rangle ={U}_{n+1}^{\dagger }| {{\rm{\Phi }}}_{T,n+1}\rangle +| {\phi }_{T}\rangle \langle {\phi }_{T}| {\psi }_{n}\rangle .\end{eqnarray} \tag{ A.7 }$$

This allows us to use a forward-backward scheme to evaluate the gradient expression in equation (A.5).

For gate operations, the cost contribution analogous to ${C}_{3}^{{st}}$ is given by ${C}_{3}^{g}=1-{\sum }_{n^{\prime} =1}^{N}{\left|\mathrm{tr}\left({U}_{T}^{\dagger }{U}_{n^{\prime} ,1}\right)\right|}^{2}/{{Nd}}^{2}$, where ${U}_{n^{\prime} ,1}={U}_{n^{\prime} }\cdots {U}_{1}$ describes the evolution from time step 1 to $n^{\prime}$. We obtain the gradient

$$\begin{eqnarray}&&\displaystyle \frac{\partial {C}_{3}^{g}}{\partial {a}_{n}}=-\displaystyle \frac{2}{{{Nd}}^{2}}\displaystyle \sum _{h=1}^{d}\mathrm{Re}\left(\left\langle {{\rm{\Phi }}}_{T,n}^{h}| \displaystyle \frac{\partial }{\partial {{\rm{a}}}_{n}}{{\rm{U}}}_{n}| {\psi }_{n-1}^{h}\right\rangle \right),\end{eqnarray} \tag{ A.8 }$$

with

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{T,n}^{h}\rangle := \displaystyle \sum _{n^{\prime} =n}^{N}{U}_{n+1}^{\dagger }\cdots {U}_{n^{\prime} }^{\dagger }{U}_{T}| {\psi }_{0}^{h}\rangle \mathrm{tr}({U}_{n^{\prime} ,1}{U}_{T}^{\dagger }).\end{eqnarray} \tag{ A.9 }$$

Here, $\{| {\psi }_{0}^{h}\rangle \}$ forms an arbitrary orthonormal basis and h = 1, 2,...,d enumerates the basis states. The vector $| {{\rm{\Phi }}}_{T,n}^{h}\rangle$ obeys the recursion relation

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{T,n}^{h}\rangle ={U}_{n+1}^{\dagger }| {{\rm{\Phi }}}_{T,n+1}^{h}\rangle +{U}_{T}| {\psi }_{0}^{h}\rangle \mathrm{tr}({U}_{n,1}{U}_{T}^{\dagger }),\end{eqnarray} \tag{ A.10 }$$

again enabling the evaluation of the gradient expression with a forward-backward propagation scheme.

For the cost contributions discussed in appendix A, we now analyze the memory usage and runtime scaling of full-AD, semi-AD and HG.

B.1. Scaling analysis: hard-coded gradient

B.2. Scaling analysis: full-automatic differentiation

The approximation of e^A ∣Ψ〉 in section 3.2 required the determination of an upper bound E_A (m, s) for the error ε_A,∣Ψ〉(m, s) [equation (14)]. In this appendix, we show how to acquire such bound.

Considering the expression for the error [equation (14)], we note that the denominator is simply 1. The error ε_A,∣Ψ〉 is bounded as follows:

$$\begin{eqnarray}\begin{array}{rcl} & & {\varepsilon }_{A,| {\rm{\Psi }}\rangle }={\left|\left|{e}^{A}| {\rm{\Psi }}\rangle -{\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle +{\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle -[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle ]\kern-2pt] \right|\right|}_{2}\\ & & \quad \leqslant {\underbrace{{\left|\left|{e}^{A}| {\rm{\Psi }}\rangle -{\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle \right|\right|}_{2}}}_{{\varepsilon }_{t}}+{\underbrace{{\left|\left|{\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle -[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{s}| {\rm{\Psi }}\rangle ]\kern-2pt] \right|\right|}_{2}}}_{{\varepsilon }_{r}}.\end{array}\end{eqnarray} \tag{ C.1 }$$

The error ε_t arises from the truncation of the Taylor series and rounding error ε_r is due to finite-precision arithmetics. We proceed to derive the upper bounds for these two errors separately.

C.1. Upper bound for the truncation error ε_t

We rewrite the truncation error as

$$\begin{eqnarray}&&{\varepsilon }_{t}={\left|\left|{e}^{A}| {\rm{\Psi }}\rangle -{\left({e}^{B}-{R}_{m}(B)\right)}^{s}| {\rm{\Psi }}\rangle \right|\right|}_{2}={\left|\left|\displaystyle \sum _{i=1}^{s}\displaystyle \frac{s!}{i!(s-i)!}{\left(-{R}_{m}(B)\right)}^{i}{\left({e}^{B}\right)}^{s-i}| {\rm{\Psi }}\rangle \right|\right|}_{2},\end{eqnarray} \tag{ C.2 }$$

where ${R}_{m}(B)={\sum }_{q\,=\,m\,+\,1}^{\infty }{B}^{q}/q!$ denotes the matrix-valued error of the truncated exponential series. Employing the triangle inequality and submultiplicativity [52], we find the bound

$$\begin{eqnarray}&&{\varepsilon }_{t}\leqslant \displaystyle \sum _{i=1}^{s}\displaystyle \frac{s!}{i!(s-i)!}{\left({R}_{m}({\left|\left|B\right|\right|}_{2})\right)}^{i}{\left|\left|{e}^{B}\right|\right|}_{2}^{s-i}{\left|\left|| {\rm{\Psi }}\rangle \right|\right|}_{2}=\displaystyle \sum _{i=1}^{s}\displaystyle \frac{s!}{i!(s-i)!}{\left({R}_{m}({\left|\left|B\right|\right|}_{2})\right)}^{i}.\end{eqnarray} \tag{ C.3 }$$

Here the second line utilizes the fact that the state is normalized and e^B is unitary.

This upper bound must be evaluated numerically to determine m and s. However, computing the matrix 2-norm is inconvenient since it requires the use of an eigensolver. To mitigate this, we switch to the matrix 1-norm by exploiting ⁵

$$\begin{eqnarray}&&{\left|\left|B\right|\right|}_{2}\leqslant {\left|\left|B\right|\right|}_{1}.\end{eqnarray} \tag{ C.4 }$$

Monotonicity of R_m leads us to

$$\begin{eqnarray}&&{\varepsilon }_{t}\leqslant \displaystyle \sum _{i=1}^{s}\displaystyle \frac{s!}{i!(s-i)!}{\left({R}_{m}({\left|\left|B\right|\right|}_{1})\right)}^{i}\lt {{sR}}_{m}({\left|\left|B\right|\right|}_{1})\displaystyle \frac{1-{\left({{sR}}_{m}({\left|\left|B\right|\right|}_{1})\right)}^{s}}{1-{{sR}}_{m}({\left|\left|B\right|\right|}_{1})}.\end{eqnarray} \tag{ C.5 }$$

C.2. Upper bound for the rounding error ε_r

The rounding error ε_r originates from the finite precision of the floating-point representation of real numbers as well as the finite accuracy of basic arithmetic operations. The floating-point representation 〚ξ〛 for the real number ξ can be written as 〚ξ〛 = ξ(1 + δ) where δ = δ(ξ) is the relative deviation [52]. Its magnitude is always smaller than machine precision, i.e. ∣δ∣ < u. Similarly, the complex number z obeys

$$\begin{eqnarray}&&[\kern-2pt[ z]\kern-2pt] =z(1+\delta )\end{eqnarray} \tag{ C.6 }$$

with ∣δ∣ ≤ u. Arithmetic operations such as addition and multiplication with two complex numbers z₁, z₂ incur a relative deviation δ = δ(z₁, z₂) such that

$$\begin{eqnarray}&&[\kern-2pt[ {z}_{1}+{z}_{2}]\kern-2pt] =([\kern-2pt[ {z}_{1}]\kern-2pt] +[\kern-2pt[ {z}_{2}]\kern-2pt] )(1+\delta ),| \delta | \leqslant u,\quad [\kern-2pt[ {z}_{1}{z}_{2}]\kern-2pt] =[\kern-2pt[ {z}_{1}]\kern-2pt] [\kern-2pt[ {z}_{2}]\kern-2pt] (1+\delta ),| \delta | \leqslant u^{\prime} ,\end{eqnarray} \tag{ C.7 }$$

where $u^{\prime} := 2\sqrt{2}u/(1-2u)$ [52]. The relative deviation shown in equations (C.6), (C.7) generally depend on both numbers involved as well as the operation in question. In the following, we will distinguish different δ with subscript ν.

Using the equations (C.6), (C.7), we derive an upper bound for the rounding error in the evaluation of the dot product which is the basic arithmetic operation in evaluating ${e}_{m}^{B}| {\rm{\Psi }}\rangle$ according to equation (13).

C.2.1. Upper bound for the rounding error associated with evaluating the dot product

The rounding error incurred when evaluating the dot product is given by ∣Σ_d − 〚Σ_d 〛∣ with Σ_d = w ^T z (${\boldsymbol{w}},{\boldsymbol{z}}\in {{\mathbb{C}}}^{d}$). Here, w ^T is a row of B = i δ tH/s. We illustrate how to bound the rounding error for d=1,2, and give the general result subsequently. Based on equations (C.6) and (C.7), we have

$$\begin{eqnarray}&&| [\kern-2pt[ {{\rm{\Sigma }}}_{1}]\kern-2pt] -{{\rm{\Sigma }}}_{1}| =| [\kern-2pt[ {w}_{1}{z}_{1}]\kern-2pt] -{w}_{1}{z}_{1}| \lt | {w}_{1}| | {z}_{1}| \left|\displaystyle \prod _{\nu =1}^{3}(1+{\delta }_{\nu })-1\right|\lt {\gamma }_{3}| {w}_{1}| | {z}_{1}| .\end{eqnarray} \tag{ C.8 }$$

In the last step, we have used $\left|{\prod }_{\nu =1}^{n}(1+{\delta }_{\nu })-1\right|\lt \tfrac{{nu}^{\prime} }{1-{nu}^{\prime} }$ [52], and defined ${\gamma }_{n}:= \tfrac{{nu}^{\prime} }{1-{nu}^{\prime} }$. For d = 2, the upper bound of the error is

$$\begin{eqnarray}&&| [\kern-2pt[ {{\rm{\Sigma }}}_{2}]\kern-2pt] -{{\rm{\Sigma }}}_{2}| =| ([\kern-2pt[ {{\rm{\Sigma }}}_{1}]\kern-2pt] +[\kern-2pt[ {w}_{2}{z}_{2}]\kern-2pt] )(1+{\delta }_{1})-{{\rm{\Sigma }}}_{2}| \lt {\gamma }_{4}\displaystyle \sum _{i=1}^{2}| {w}_{i}| | {z}_{i}| .\end{eqnarray} \tag{ C.9 }$$

The general upper bound for arbitrary d is given by:

$$\begin{eqnarray}&&| [\kern-2pt[ {{\rm{\Sigma }}}_{d}]\kern-2pt] -{{\rm{\Sigma }}}_{d}| \lt {\gamma }_{d+2}\displaystyle \sum _{j=1}^{d}| {w}_{j}| | {z}_{j}| .\end{eqnarray} \tag{ C.10 }$$

Since w ^T is proportional to a row of the Hamiltonian matrix, sparsity of H implies sparsity of w ^T. Vanishing entries in w ^T is special because they do not incur rounding errors. Thus, one can obtain a stronger upper bound

$$\begin{eqnarray}&&| [\kern-2pt[ {{\rm{\Sigma }}}_{d}]\kern-2pt] -{{\rm{\Sigma }}}_{d}| \lt {\gamma }_{\sigma +2}\displaystyle \sum _{j=1}^{d}| {w}_{j}| | {z}_{j}| \end{eqnarray} \tag{ C.11 }$$

where σ is the number of non-zero elements in w ^T.

C.2.2. Upper bound for the rounding error of ${e}_{m}^{B}| {\rm{\Psi }}\rangle$

The vector ${e}_{m}^{B}| {\rm{\Psi }}\rangle$ is evaluated recursively via

$$\begin{eqnarray}&&{b}_{j}=\displaystyle \frac{B}{j}{b}_{j-1},\qquad {e}_{j}^{B}{b}_{0}={b}_{j-1}+{b}_{j}\end{eqnarray} \tag{ C.12 }$$

where j = 1, 2, ⋯ ,m and b₀ is a vector representation of ∣Ψ〉 under a choice of orthonormal basis. Before evaluating the upper bound for the error ${\left|\left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] -{e}_{m}^{B}{b}_{0}\right|\right|}_{2}$, we first derive an upper bound for the error of an vector component $\left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}{]\kern-2pt] }^{(i)}-{\left({e}_{m}^{B}{b}_{0}\right)}^{(i)}\right|$. We illustrate how to bound this error for m = 0,1, and then give the general result.

For the case m = 0, we have the bound

$$\begin{eqnarray}&&\left|[\kern-2pt[ {e}_{0}^{B}{b}_{0}{]\kern-2pt] }^{(i)}-{\left({e}_{0}^{B}{b}_{0}\right)}^{(i)}\right|=\left|[\kern-2pt[ {b}_{0}{]\kern-2pt] }^{(i)}-{b}_{0}^{(i)}\right|\lt {\gamma }_{2}| {b}_{0}^{(i)}| .\end{eqnarray} \tag{ C.13 }$$

We then obtain the upper bound of the error for the example of m = 1:

$$\begin{eqnarray}\begin{array}{rcl} & & | [\kern-2pt[ {e}_{1}^{B}{b}_{0}{]\kern-2pt] }^{(i)}-{\left({e}_{1}^{B}{b}_{0}\right)}^{(i)}| =\left|\left([\kern-2pt[ {b}_{0}{]\kern-2pt] }^{(i)}+[\kern-2pt[ {b}_{1}{]\kern-2pt] }^{(i)}\right)(1+{\delta }_{1})-{\left({e}_{1}^{B}{b}_{0}\right)}^{(i)}\right|\\ & & \quad =\left|[\kern-2pt[ {b}_{0}{]\kern-2pt] }^{(i)}(1+{\delta }_{1})+{\displaystyle \frac{[\kern-2pt[ {{Ab}}_{0}]\kern-2pt] }{s}}^{(i)}\displaystyle \prod _{\nu =1}^{3}{\left(1+{\delta }_{\nu }\right)}^{3}-{\left({e}_{1}^{B}{b}_{0}\right)}^{(i)}\right|\mathop{\lt }\limits^{C.11}{\gamma }_{3}| {b}_{0}^{(i)}| +{\gamma }_{{\sigma }_{i}+5}\displaystyle \sum _{l}| {B}^{({il})}| | {b}_{0}^{(l)}| ,\end{array}\end{eqnarray} \tag{ C.14 }$$

where σ_i is number of non-zero elements in ith row of the matrix B. Here, the third line arises from the division error between a complex number and a real number, i.e. 〚z/ξ〛 = z/ξ(1 + δ), ∣δ∣ < u [52]. The general upper bound for arbitrary m is given by:

$$\begin{eqnarray}&&\left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}{]\kern-2pt] }^{(i)}-{\left({e}_{m}^{B}{b}_{0}\right)}^{(i)}\right|\lt \displaystyle \sum _{k=0}^{m}{\gamma }_{k({\sigma }_{i}+2)+m+2}{\left(\displaystyle \frac{| B{| }^{k}}{k!}| {b}_{0}| \right)}^{i},\end{eqnarray} \tag{ C.15 }$$

where ∣B∣, ∣b₀∣ denote the matrix and the vector with elements ∣B^il ∣, $| {b}_{0}^{i}|$, respectively. Defining

$$\begin{eqnarray}&&\sigma ^{\prime} := \mathop{\max }\limits_{i}{\sigma }_{i}\ \end{eqnarray} \tag{ C.16 }$$

and using equation (C.15), we obtain

$$\begin{eqnarray}\begin{array}{rcl} & & {\left|\left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] -({e}_{m}^{B}{b}_{0})\right|\right|}_{2}={\parallel \left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] -({e}_{m}^{B}{b}_{0})\right|\parallel }_{2}\lt \displaystyle \sum _{k=0}^{m}{\gamma }_{k(\sigma ^{\prime} +2)+m+2}{\parallel \displaystyle \frac{{\left|B\right|}^{k}}{k!}| {b}_{0}| \parallel }_{2}\\ & & \quad \lt \displaystyle \sum _{k=0}^{m}{\gamma }_{k(\sigma ^{\prime} +2)+m+2}\displaystyle \frac{{\parallel \left|B\right|\parallel }_{2}^{k}}{k!}{\parallel \left|{b}_{0}\right|\parallel }_{2}\mathop{\lt }\limits^{({\rm{C}}.4)}\displaystyle \sum _{k=0}^{m}{\gamma }_{k(\sigma ^{\prime} +2)+m+2}\displaystyle \frac{{\parallel \left|B\right|\parallel }_{1}^{k}}{k!}{\parallel \left|{b}_{0}\right|\parallel }_{2}=\beta {\parallel {b}_{0}\parallel }_{2}\end{array}\end{eqnarray} \tag{ C.17 }$$

with $\beta := {\sum }_{k=0}^{m}{\gamma }_{k(\sigma ^{\prime} +2)+m\,+\,2}\tfrac{| | B| {| }_{1}^{k}}{k!}$. Here, β∣∣b₀∣∣₂ the desired upper bound for the rounding error of ${e}_{m}^{B}| {\rm{\Psi }}\rangle$.

C.2.3. Upper bound for the error ε_r

The vector ${\left({e}_{m}^{B}\right)}^{s}{b}_{0}$ is evaluated by the recursive relation.

$$\begin{eqnarray}&&{c}_{l}={e}_{m}^{B}{c}_{l-1},\qquad l=1,2,\cdots ,s,\end{eqnarray} \tag{ C.18 }$$

where c₀ = b₀. We illustrate how to bound the error ${\varepsilon }_{r}={\left|\left|[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{s}{b}_{0}]\kern-2pt] -{\left({e}_{m}^{B}\right)}^{s}{b}_{0}\right|\right|}_{2}$ for s = 1, 2, and give the general result subsequently.

For s = 1, we have

$$\begin{eqnarray}&&{\left|\left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] -{e}_{m}^{B}{b}_{0}\right|\right|}_{2}\mathop{\lt }\limits^{(C.17)}\beta | | {b}_{0}| {| }_{2}=\beta .\end{eqnarray} \tag{ C.19 }$$

For s = 2, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\left|\left|[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{2}{b}_{0}]\kern-2pt] -{\left({e}_{m}^{B}\right)}^{2}{b}_{0}\right|\right|}_{2} & & \leqslant {\left|\left|[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{2}{b}_{0}]\kern-2pt] -{e}_{m}^{B}[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] \right|\right|}_{2}+{\left|\left|{e}_{m}^{B}[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] -{\left({e}_{m}^{B}\right)}^{2}{b}_{0}\right|\right|}_{2}\\ & & \lt {\left|\left|[\kern-2pt[ {e}_{m}^{B}[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] ]\kern-2pt] -{e}_{m}^{B}[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] \right|\right|}_{2}+\beta {\left|\left|{e}_{m}^{B}\right|\right|}_{2},\end{array}\end{eqnarray} \tag{ C.20 }$$

where the last step is enabled by equation (C.19). Finally, we bound the matrix norm in the second term by

$$\begin{eqnarray}&&{\left|\left|{e}_{m}^{B}\right|\right|}_{2}={\left|\left|{e}^{B}-{R}_{m}(B)\right|\right|}_{2}\leqslant 1+{R}_{m}({\left|\left|B\right|\right|}_{2})\mathop{\leqslant }\limits^{({\rm{C}}.4)}1+{R}_{m}({\left|\left|B\right|\right|}_{1})=: \alpha .\end{eqnarray} \tag{ C.21 }$$

Since equation (C.17) holds for any complex vector b₀, the error in equation (C.20) is further bounded by:

$$\begin{eqnarray}&&\begin{array}{l}{\left|\left|[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{2}{b}_{0}]\kern-2pt] -{\left({e}_{m}^{B}\right)}^{2}{b}_{0}\right|\right|}_{2}\mathop{\lt }\limits^{({\rm{C}}.17)}\beta {\left|\left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] \right|\right|}_{2}+\beta \alpha \\ =\beta {\left|\left|[\kern-2pt[ {e}_{m}^{B}{b}_{0}]\kern-2pt] -{e}_{m}^{B}{b}_{0}+{e}_{m}^{B}{b}_{0}\right|\right|}_{2}+\beta \alpha \quad \mathop{\lt }\limits^{({\rm{C}}.19)}\beta (\beta +\alpha )+\beta \alpha .\end{array}\end{eqnarray} \tag{ C.22 }$$

Generalizing this strategy further, we obtain the expression for the upper bound of ε_r :

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{r}={\left|\left|[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{s}{b}_{0}]\kern-2pt] -{\left({e}_{m}^{B}\right)}^{s}{b}_{0}\right|\right|}_{2} & \leqslant & \displaystyle \sum _{l=0}^{s-1}{\left|\left|{\left({e}_{m}^{B}\right)}^{l}[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{s-l}{b}_{0}]\kern-2pt] -{\left({e}_{m}^{B}\right)}^{l+1}[\kern-2pt[ {\left({e}_{m}^{B}\right)}^{s-l-1}{b}_{0}]\kern-2pt] \right|\right|}_{2}\\ & \lt & \beta \displaystyle \sum _{l=0}^{s-1}{\left(\alpha +\beta \right)}^{s-l-1}{\alpha }^{l}={\left(\alpha +\beta \right)}^{s}-{\alpha }^{s}.\end{array}\end{eqnarray} \tag{ C.23 }$$

Combining the upper bounds for ε_r and ε_t , we return to equation (C.1) and obtain

$$\begin{eqnarray*}&&{\varepsilon }_{A,| {\rm{\Psi }}\rangle }\lt {E}_{A}(m,s)={{sR}}_{m}({\left|\left|B\right|\right|}_{1})\displaystyle \frac{1-{\left({{sR}}_{m}({\left|\left|B\right|\right|}_{1})\right)}^{s}}{1-{{sR}}_{m}({\left|\left|B\right|\right|}_{1})}+{\left(\alpha +\beta \right)}^{s}-{\alpha }^{s},\end{eqnarray*}$$

where ${\left|\left|B\right|\right|}_{1}={\left|\left|A\right|\right|}_{1}/s$, $\alpha =\alpha ({\left|\left|A\right|\right|}_{1},m,s)$ and $\beta =\beta (\sigma ^{\prime} ,{\left|\left|A\right|\right|}_{1},m,s)$.

Scaling and squaring is used to evaluate the propagator-state product, e^−iHΔt ∣Ψ〉. We compare relative error and runtime between our implementation of scaling and squaring and Scipy's implementation, expm multiply [42]. The relative errors for our implementation and Scipy's implementation follow the general expression in equation (14), but differ in specific choice of m and s. As a proxy for the exact propagator-state product, we apply Taylor expansion to matrix exponential with 400 digits precision by Sympy and the result of the product converges at 30 digits after decimal point. The runtime is measured with the help of the Temci package [45].

As the basis for this comparison, we use the examples of Hamiltonians (20) and (21) and implement sparse storage for them. The relative errors are shown in figures D1(a) and (c) as a function versus the matrix norm ${\left|\left|-{iH}{\rm{\Delta }}t\right|\right|}_{2}$ respectively. In both examples, the relative errors for our implementation are bounded by the error tolerance τ, while the relative errors for Scipy's implementation tend to exceed the tolerance for large norms There are two reasons for this tendency. First, rounding errors are not considered in the derivation of the upper bound for the errors [42]; second the truncation of e^B _m [equation (12)] is merely based on a heuristic criterion.

Zoom In Zoom Out Reset image size

Even though the relative errors of Scipy's implementation are worse, it gains smaller number of matrix-vector multiplications, which is denoted by μ. This is illustrated in figure D1(b) and (d), that the ratio of μ between our implementation and Scipy's implementation is always larger than 1. Although Scipy's implementation spends less time on matrix-vector multiplications, the total runtime of Scipy's implementation is larger than the one of our implementation, see figure D1(b) and (d). This is due to the larger runtime overhead of Scipy's implementation for determining (m, s).

In conclusion, our implementation has the better numerical stability and runtime than Scipy's implementation in our numerical experiments.

The asymptotic behavior of runtime, discussed in the main text, is governed by the scaling of μ (16) with dimension d of the Hilbert space. The dependence μ(d) is specific to the matrix A to be exponentiated (which is proportional to the Hamiltonian). This prevents us from stating a general scaling relation. We thus turn to a discussion of several concrete examples.

We first consider the example of the transmon-cavity system, see equation (20), where we increase d via the cavity dimension d_c while leaving the transmon cutoff fixed. The corresponding Hamiltonian H₁, written in the bare product basis, then has a constant maximum number of non-zero elements in each row. (I.e., the parameter $\sigma ^{\prime}$ (C.16) is independent of d.) Numerically, we find in this case that μ scales linearly with ${\left|\left|A\right|\right|}_{1}$ ($\mu \sim {\rm{\Theta }}({\left|\left|A\right|\right|}_{1})$), and this linearity holds for a wide range of τ, see figure E1(a). In this example, two ingredients lead to the scaling ${\left|\left|A\right|\right|}_{1}\sim {\rm{\Theta }}({d}^{\tfrac{1}{2}})$ : (1) d is a constant multiple of d_c , and (2) the cavity ladder operators in H₁ which has one-norm scaling ${\rm{\Theta }}({d}_{c}^{\tfrac{1}{2}})$. Thus, we obtain the overall scaling $\mu \sim {\rm{\Theta }}({d}^{\tfrac{1}{2}})$, see figure E1(e). In the second example, we consider the system three coupled transmons, see equation (21). The Hamiltonian H₂ is written in the harmonic oscillator basis and has a constant $\sigma ^{\prime}$, which leads to $\mu \sim {\rm{\Theta }}({\left|\left|A\right|\right|}_{1})$ [figure E1(b)]. One-norm of operator ${\left({b}^{\dagger }b\right)}^{2}$ in H₂ has quadratic scaling with each transmon dimension. Since we increase d by enlarging the dimensions of all three transmons simultaneously, we obtain ${\left|\left|A\right|\right|}_{1}\sim {\rm{\Theta }}({d}^{\tfrac{2}{3}})$. This scaling results in $\mu \sim {\rm{\Theta }}({d}^{\tfrac{2}{3}})$, see figure E1 (f).

Zoom In Zoom Out Reset image size

For other Hamiltonian matrices where $\sigma ^{\prime}$ depends on d, the scaling $\mu \sim {\rm{\Theta }}({\left|\left|A\right|\right|}_{1})$ may still hold. In the following, we show two concrete examples that one where linearity holds, and another where it breaks. We first consider a driven system of N_q qubits with nearest-neighbor σ_z coupling. In the rotating frame, the system is described by

$$\begin{eqnarray*}&&{H}_{3}(t)=\displaystyle \sum _{\nu =1}^{{N}_{q}}[{a}_{x}^{(\nu )}(t){\sigma }_{x}^{(\nu )}+{a}_{y}^{(\nu )}(t){\sigma }_{y}^{(\nu )}]+\displaystyle \sum _{\nu =1}^{{N}_{q}-1}g{\sigma }_{z}^{(\nu )}{\sigma }_{z}^{(\nu +1)}.\end{eqnarray*}$$

Each qubit is controlled via drive fields ${a}_{x}^{(n)}$ and ${a}_{y}^{(n)}$. For this example, the numerical results show that scaling $\mu \sim {\rm{\Theta }}({\left|\left|A\right|\right|}_{1})$ still holds, see figure E1(c). Since we increase d by enlarging N_q , we obtain ${\left|\left|A\right|\right|}_{1}\sim {\rm{\Theta }}({N}_{q})\sim {\rm{\Theta }}({\mathrm{log}}_{2}d)$. This scaling results in $\mu \sim {\rm{\Theta }}({\mathrm{log}}_{2}d)$, see figure E1(g). For the second example, we consider two capacitively coupled fluxonium qubits with Hamiltonian

$$\begin{eqnarray}&&\begin{array}{l}{H}_{4}(t)=\displaystyle \sum _{i=1}^{2}\left(4{E}_{{Ci}}{n}_{i}^{2}-{E}_{{Ji}}\cos ({\varphi }_{i})+\displaystyle \frac{{E}_{{Li}}}{2}{\left[{\varphi }_{i}-{\phi }_{i}(t)\right]}^{2}\right)+{{gn}}_{1}{n}_{2}.\end{array}\end{eqnarray} \tag{ E.1 }$$

Here, E_Ci , E_Ji , E_Li denote the charging, inductive and Josephson energies of two fluxonium qubits, and g is the interaction strength. φ and n are the phase and charge number operators, defined in the usual way. We represent the Hamiltonian in the harmonic oscillator basis. Figure E1(d) shows that scaling of μ is superlinear with ${\left|\left|A\right|\right|}_{1}$ for τ = 10⁻⁸. We thus settle for a numerical determination of the scaling relation between μ and d for τ = 10⁻⁸ by increasing the dimensional cutoff for both fluxonium qubits simultaneously. The observed scaling is μ ∼ Θ(d^0.65), see figure E1(h).

The numerical implementation of the hard-coded gradient scheme requires evaluation of the short-time propagator and its derivative. Under certain conditions, it turns out to be advantageous to numerically evaluate the short-time propagator by employing matrix diagonalization. In this appendix, we first discuss the method to compute propagator derivative and then analyze the scaling of HG.

F.1. Evaluation of propagator derivative

The operator −iHΔt is anti-Hermitian, and hence can be diagonalized and exponentiated according to

$$\begin{eqnarray}&&A={{SDS}}^{\dagger },\end{eqnarray} \tag{ F.1 }$$

$$\begin{eqnarray}&&{e}^{A}={{Se}}^{D}{S}^{\dagger }.\end{eqnarray} \tag{ F.2 }$$

Here, S is unitary and D diagonal. As a result, the propagator derivative can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}^{A}}{{da}}={{Se}}^{D}\displaystyle \frac{{dD}}{{da}}{S}^{\dagger }+\displaystyle \frac{{dS}}{{da}}{e}^{D}{S}^{\dagger }+{{Se}}^{D}\displaystyle \frac{{{dS}}^{\dagger }}{{da}}.\end{eqnarray} \tag{ F.3 }$$

We apply the inverse unitary transformation to $\tfrac{{{de}}^{A}}{{da}}$ (F.3) and obtain

$$\begin{eqnarray}&&{S}^{\dagger }\displaystyle \frac{{{de}}^{A}}{{da}}S={e}^{D}\displaystyle \frac{{dD}}{{da}}+{S}^{\dagger }\displaystyle \frac{{dS}}{{da}}{e}^{D}+{e}^{D}\displaystyle \frac{{{dS}}^{\dagger }}{{da}}S.\end{eqnarray} \tag{ F.4 }$$

Based on $\tfrac{{{dS}}^{\dagger }}{{da}}S+{S}^{\dagger }\tfrac{{dS}}{{da}}=0$, equation (F.4) can be rewritten as

$$\begin{eqnarray}&&{S}^{\dagger }\displaystyle \frac{{{de}}^{A}}{{da}}S={e}^{D}\displaystyle \frac{{dD}}{{da}}+E\circ \left({S}^{\dagger }\displaystyle \frac{{dS}}{{da}}\right)\end{eqnarray} \tag{ F.5 }$$

with ${E}_{{ij}}:= {e}^{{D}_{{jj}}}-{e}^{{D}_{{ii}}}$. Here, ◦denotes the Hadamard product (element-wise multiplication). To obtain dD/da and dS/da, we take the derivative of equation (F.1) and repeat the same steps as for equations (F.4) and (F.5), leading to

$$\begin{eqnarray}&&{S}^{\dagger }\displaystyle \frac{{dA}}{{da}}S=\displaystyle \frac{{dD}}{{da}}+E^{\prime} \circ \left({S}^{\dagger }\displaystyle \frac{{dS}}{{da}}\right),\end{eqnarray} \tag{ F.6 }$$

where $E{{\prime} }_{{ij}}:= {D}_{{jj}}-{D}_{{ii}}$. Since the diagonal elements of $E^{\prime}$ are zero, it follows that

$$\begin{eqnarray}&&\displaystyle \frac{{dD}}{{da}}=I\circ \left({S}^{\dagger }\displaystyle \frac{{dA}}{{da}}S\right).\end{eqnarray} \tag{ F.7 }$$

For off-diagonal elements of the operator in equation (F.6), we obtain

$$\begin{eqnarray}&&F\circ \left({S}^{\dagger }\displaystyle \frac{{dA}}{{da}}S\right)+G={S}^{\dagger }\displaystyle \frac{{dS}}{{da}},\end{eqnarray} \tag{ F.8 }$$

where we define

$$\begin{eqnarray}\begin{array}{l}{F}_{{ij}}:= \left\{\begin{array}{ll}\displaystyle \frac{1}{{D}_{{jj}}-{D}_{{ii}}}, & \mathrm{if}\,{D}_{{jj}}-{D}_{{ii}}\ne 0.\\ 0, & \mathrm{otherwise}.\end{array}\right.\\ {G}_{{ij}}:= \left\{\begin{array}{ll}0, & \mathrm{if}\,{D}_{{jj}}-{D}_{{ii}}\ne 0.\\ {\left({S}^{\dagger }\displaystyle \frac{{dS}}{{da}}\right)}_{{ij}}, & \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray} \tag{ F.9 }$$

Inserting equations (F.7) and (F.8) into equation (F.5) gives

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}^{A}}{{da}}=S\left(\left({e}^{D}+E\circ F)\circ ({S}^{\dagger }\displaystyle \frac{{dA}}{{da}}S\right)\right){S}^{\dagger },\end{eqnarray} \tag{ F.10 }$$

where we use E ◦ G = 0.

F.2. Scaling analysis