Numerical Optimal Control and Gradient Checks

Numerical optimal control adjusts a time-dependent Hamiltonian to achieve a target state or gate under explicit constraints. An optimizer can exploit errors in the propagator or an unrealistically flexible pulse model. The central computational task is therefore to optimize a well-defined discrete objective and then verify the resulting control with a finer, independent simulation.

Helpful background. Quantum control objectives provides the operational setting; matrix exponentials provides finite-dimensional propagation.

Piecewise controls and objective functions

Let

$$H(t)=H_0+\sum_a u_a(t)H_a,$$

with Hermitian drift and control operators. Divide the duration into $N$ slices of width $\Delta t$, using constant amplitudes $u_{a,j}$ on each slice. Then $U_j=e^{-iH_j\Delta t/\hbar}$ and $|\psi_N\rangle=U_N\cdots U_1|\psi_0\rangle$.

For a normalized target state, set $A=\langle\psi_{\rm tar}|\psi_N\rangle$ and maximize $F=|A|^2$. A gate objective such as $|\operatorname{Tr}(U_{\rm tar}^\dagger U)|^2/d^2$ tests an entire $d$-dimensional unitary rather than one input state. State transfer alone does not certify a gate.

Include amplitude bounds, bandwidth or slew limits, and any fluence penalty in the actual optimization. For example minimize

$$J=1-F+\lambda\Delta t\sum_{a,j}u_{a,j}^2.$$

The coefficient $\lambda$ has units fixed by the amplitude convention. A post hoc smoothing operation changes the simulated pulse and must be followed by a fresh fidelity calculation.

Exact slice derivatives and adjoint reuse

Define forward states $|\psi_{j-1}\rangle=U_{j-1}\cdots U_1|\psi_0\rangle$ and backward bras $\langle\chi_j|=\langle\psi_{\rm tar}|U_N\cdots U_{j+1}$. Then

$$\frac{\partial F}{\partial u_{a,j}} =2\operatorname{Re}\left[ A^*\langle\chi_j| \frac{\partial U_j}{\partial u_{a,j}} |\psi_{j-1}\rangle\right].$$

Forward and backward sweeps reuse propagation instead of separately rerunning the entire trajectory for each parameter. The exact matrix derivative is

$$\frac{\partial U_j}{\partial u_{a,j}} =-\frac{i}{\hbar}\int_0^{\Delta t} e^{-iH_j(\Delta t-s)/\hbar} H_a e^{-iH_js/\hbar}\,ds.$$

Replacing this by $-(i\Delta t/\hbar)H_aU_j$ is exact only when $H_a$ commutes with $H_j$ and otherwise has a slice-dependent error. Fréchet derivatives or divided differences in an eigenbasis avoid that approximation. Check gradients with centered finite differences over several perturbation sizes before optimizing.

An exact rotation benchmark

For $H(t)=\hbar u(t)\sigma_x/2$, all time slices commute. Starting in $|0\rangle$, the probability of $|1\rangle$ is

$$F=\sin^2(\Theta/2),\qquad \Theta=\int_0^T u(t)\,dt.$$

Thus $\partial F/\partial u_j=(\Delta t/2)\sin\Theta$ for constant-width slices. A pulse area $\Theta=\pi$ gives unit transfer. If $|u|\leq u_{\max}$, this model requires at least $T=\pi/u_{\max}$ for the direct rotation. This is a bound for this specified control model, not a general quantum speed-limit claim.

The benchmark detects missing factors of two, inconsistent angular-frequency units, and incorrect propagation order. Adding a detuning introduces noncommuting dynamics and provides a stronger gradient test.

Robustness and propagation accuracy

Optimize an ensemble of detunings or amplitude errors when robustness is part of the goal, then evaluate a denser validation ensemble not used in training. Report nominal, average, and worst-case infidelity separately. A high mean fidelity can conceal a narrow failure region.

Refine time slices, tighten exponential-action tolerances, and propagate the final control with a different method. Include the transfer function of the intended control hardware when it is part of the model. Distinguish local optimizer convergence from global optimality; repeated starts and different parameterizations explore the landscape without proving a global optimum.

Exercises

At $\Theta=\pi$, the exact first derivative vanishes. Does this mean a small pulse-area error has no effect?

Solution

For $\Theta=\pi+\delta$, $F=\cos^2(\delta/2)=1-\delta^2/4+O(\delta^4)$. The leading error is quadratic, not zero. A vanishing gradient at a perfect transfer point does not guarantee broad robustness.

References

• Khaneja, N., T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, “Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms,” Journal of Magnetic Resonance 172, 296–305 (2005), doi:10.1016/j.jmr.2004.11.004 — forward/backward gradient optimization.
• Nocedal, J., and S. J. Wright, Numerical Optimization, 2nd ed., Springer, 2006, doi:10.1007/978-0-387-40065-5 — constrained objectives and convergence tests.