Every optimal-control algorithm spends most of its time propagating: computing the time evolution of states or operators under the current controls — and, for gradient methods, its derivatives — over and over again. How fast and how accurately that is done sets the size of the systems and the number of iterations that are practical. This page collects the numerical methods behind the optimizers.

Propagating one time step

Over a short step \(\Delta t\) in which the Hamiltonian is (approximately) constant, one needs \(e^{-iH\Delta t}\) or its action on a state.

  • Dense matrix exponentials. For small systems, \(e^{-iH\Delta t}\) is computed directly. Scaling-and-squaring with Padé approximants is the standard general-purpose method. Since \(H\) is Hermitian, an eigendecomposition is also convenient, and it gives exact derivatives (below). Many methods exist and some fail badly in practice [1].
  • Krylov subspace (Lanczos). For large sparse \(H\), the action \(e^{-iH\Delta t}|\psi\rangle\) is approximated in the small subspace spanned by \(|\psi\rangle, H|\psi\rangle, H^2|\psi\rangle, \dots\), using only matrix–vector products [2].
  • Chebyshev expansion. The exponential is expanded in Chebyshev polynomials of the spectrally rescaled Hamiltonian. The error is distributed uniformly over the spectrum and the method is very efficient for long steps [3]. Variants with iterative time ordering handle explicitly time-dependent Hamiltonians [4].
  • Split-operator. For particles on spatial grids, kinetic and potential terms are exponentiated separately, in momentum and position space, using fast Fourier transforms [5].

Kosloff’s review compares these approaches for molecular dynamics [6]. Their trade-offs carry over to control.

Time-dependent Hamiltonians

Control Hamiltonians change in time. Options include:

  • piecewise-constant approximation — natural for GRAPE, with a time step small compared with the inverse spectral width of \(H(t)\);
  • Magnus integrators, which exponentiate a truncated Magnus series and stay exactly unitary, with higher-order accuracy for smooth controls [7];
  • general ODE solvers (e.g. adaptive Runge–Kutta), which are flexible and error-controlled but not exactly unitary — used, for example, by GOAT.

Derivatives of propagators

Gradient methods need \(\partial e^{-iH\Delta t}/\partial u\) where \(H = H_0 + u H_k\). Three exact approaches are common:

Eigenbasis formula. With \(H|a\rangle = \lambda_a|a\rangle\), the matrix elements of the derivative are

$$\Big\langle a\Big|\frac{\partial e^{-iH\Delta t}}{\partial u}\Big|b\Big\rangle = \langle a|H_k|b\rangle\, \frac{e^{-i\lambda_a \Delta t} - e^{-i\lambda_b \Delta t}}{\lambda_a - \lambda_b},$$

with the limit \(-i\Delta t\, e^{-i\lambda_a\Delta t}\) when \(\lambda_a = \lambda_b\).

Augmented (auxiliary) matrix. The exponential of a block-triangular matrix contains the derivative in its off-diagonal block [8]:

$$\exp\!\left[\Delta t\begin{pmatrix} A & B \\ 0 & A\end{pmatrix}\right] = \begin{pmatrix} e^{A\Delta t} & \displaystyle\int_0^{\Delta t} e^{A(\Delta t-s)}\,B\,e^{As}\,ds \\ 0 & e^{A\Delta t}\end{pmatrix},$$

and with \(A = -iH\), \(B = -iH_k\) the upper-right block is exactly \(\partial e^{-iH\Delta t}/\partial u\). This has been developed into a general tool for optimal control and spin dynamics [9].

Fréchet-derivative algorithms compute the derivative of the matrix exponential along a direction with the same scaling-and-squaring machinery as the exponential itself [10].

Automatic differentiation obtains the same derivatives automatically, provided the framework differentiates matrix functions correctly.

Scaling up

  • Open systems. Density matrices live in a space of dimension \(d^2\), so Liouville-space propagators are expensive. Stochastic quantum trajectories trade this for averaging over many state vectors [11]; see Control of open quantum systems.
  • Many-body systems. Tensor-network methods (TEBD, tDMRG, TDVP) propagate states with limited entanglement in very large Hilbert spaces, and have been coupled to optimal control from the start [12]; see Many-body optimal control.
  • Hardware. GPUs accelerate both the dense linear algebra and automatic differentiation [13].
  • Approximations as a feature. Truncated Magnus and Chebyshev expansions also turn the control problem into a polynomial one that can be solved globally (QCPOP).