GRAPE.jl: Quantum Optimal Control in Julia

• GRAPE.jl is a Julia package for quantum optimal control that uses gradient ascent pulse engineering with piecewise-constant controls to steer quantum systems.
• It employs a semi-automatic differentiation architecture that analytically computes gradients for control updates while applying AD only to final-time functionals, ensuring efficiency and accuracy.
• Built within the JuliaQuantumControl ecosystem, its modular design leverages multiple dispatch and JIT compilation to achieve performance comparable to low-level languages.

GRAPE.jl is a Julia implementation of Gradient Ascent Pulse Engineering (GRAPE) for quantum optimal control, built on top of QuantumControl.jl and QuantumPropagators.jl and designed to be both high-performance and highly extensible across a broad class of closed and open quantum control problems (Goerz et al., 3 Nov 2025). Its purpose is to find controls that steer a quantum system in a particular way, including unitary gate implementation, state transfer, robust ensemble control, and optimization under dissipative dynamics. Within the JuliaQuantumControl ecosystem, it is also the concrete realization of the “semi-automatic differentiation” formulation developed for generalized GRAPE, in which analytic propagation derivatives are combined with automatic differentiation applied only to final-time functionals (Goerz et al., 2022).

1. Conceptual basis in quantum optimal control

In quantum optimal control, the central task is to choose control fields $\epsilon_l(t)$ so that a driven quantum system realizes a specified objective. GRAPE.jl adopts the standard formulation in which the dynamics are governed by a Schrödinger-type equation,

$$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$

or, more generally, by a generator such as a Liouvillian super-operator for open-system dynamics (Goerz et al., 3 Nov 2025). The optimization target is encoded in an objective functional $J$, which may depend on the final state, on multiple propagated trajectories, and on explicit control penalties.

The defining numerical idea of GRAPE is the use of piecewise-constant controls over a time grid. The interval $[0,T]$ is divided into $N$ steps, and each control amplitude $\epsilon_{nl}$ is held constant on one interval. This yields stepwise-constant generators $\hat{H}_{k,n}$ and interval propagators

$$\hat{U}_{k,n} = \exp\!\left[ -\frac{i}{\hbar} \hat{H}_{k,n} \Delta t_n \right],$$

so that the total propagator is an ordered product of step propagators (Goerz et al., 3 Nov 2025). Because the controls are discretized in this way, GRAPE can derive analytic expressions for the gradient $\partial J / \partial \epsilon_{nl}$ and evaluate the full gradient by combining a forward propagation with a backward propagation or adjoint construction, at a computational cost that is essentially only a constant-factor overhead over a single forward evolution.

GRAPE.jl follows the generalized GRAPE scheme of Goerz et al. described in “Quantum Optimal Control via Semi-Automatic Differentiation” (Goerz et al., 2022). This generalized formulation is explicitly not restricted to the traditional linear-control Hamiltonian $\hat{H} = \hat{H}_0 + \sum_l \epsilon_l \hat{H}_l$, nor to overlap-only objectives. Instead, it supports arbitrary dependence $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$0, arbitrary sets of trajectories $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$1, and general functionals $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$2 provided that derivatives with respect to final states can be obtained analytically or via automatic differentiation (Goerz et al., 3 Nov 2025).

2. Semi-automatic differentiation and gradient evaluation

A defining feature of GRAPE.jl is its “semi-automatic differentiation” architecture. The method separates the dependence of the objective on propagated final-time quantities from the dependence of those propagated quantities on the controls. In the formulation of (Goerz et al., 2022), the functional is rewritten in terms of intermediary variables such as final states $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$3, overlaps $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$4, or a projected logical gate matrix $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$5. The chain rule then splits the gradient into two parts: $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$6 is handled analytically by GRAPE-style forward/backward propagation, while $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$7 is evaluated by automatic differentiation on a small final-time computational graph (Goerz et al., 2022).

This decomposition is motivated by the prohibitive cost of full reverse-mode AD through the entire time evolution. In the benchmark analysis of (Goerz et al., 2022), naive full AD of propagation incurs storage scaling of roughly $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$8 for matrix-exponential propagators, or $$i\hbar \frac{\partial}{\partial t} \vert\Psi(t)\rangle = \hat{H}(\epsilon(t)) \vert\Psi(t)\rangle,$$9 for polynomial propagators, whereas semi-automatic differentiation stores only forward-propagated states, with $J$0 memory scaling, and confines AD to a tiny graph defined on final-time variables. The practical implication is that GRAPE.jl retains the runtime and RAM behavior of a conventional analytic GRAPE implementation while allowing essentially any computable final-time functional (Goerz et al., 2022).

For the propagation derivative itself, GRAPE.jl uses the auxiliary-matrix formalism of Goodwin and Kuprov to obtain derivatives of the step propagator with respect to controls to machine precision (Goerz et al., 3 Nov 2025). In the theoretical presentation of (Goerz et al., 2022), this appears as a gradient-generator or Van Loan construction in which an extended linear operator simultaneously propagates the state and the derivatives of the propagator. GRAPE.jl implements this without constructing the full block matrix explicitly; instead, it applies an abstract linear operator to a block vector. This preserves efficiency and permits in-place sparse complex linear algebra in the propagation layer, while AD is used only for the much smaller functional layer.

The same framework extends beyond simple overlap functionals. In gate optimization, the co-states at final time can be constructed from AD of a logical gate matrix $J$1; in overlap-based optimization, from AD of the overlaps $J$2; and in state-based objectives, directly from AD of the final states (Goerz et al., 2022). This is what allows direct optimization of non-analytic functionals such as perfect-entangler objectives and gate concurrence, a capability highlighted as the first direct optimization of the non-analytic gate concurrence (Goerz et al., 2022).

3. Software architecture and Julia-specific design

GRAPE.jl is part of the JuliaQuantumControl ecosystem and builds tightly on QuantumControl.jl and QuantumPropagators.jl (Goerz et al., 3 Nov 2025). QuantumControl.jl provides the abstract model for trajectories, controls, functionals, and propagation options, while GRAPE.jl implements one specific optimization method: generalized GRAPE with semi-automatic differentiation for piecewise-constant controls. The package uses these abstractions to represent the dynamical model $J$3 or $J$4 for each trajectory, the mapping from discrete control amplitudes to generators, the objective functional, and the derivative of that functional with respect to final states.

The architecture is modular. A problem is specified in terms of states, Hamiltonians or Liouvillians, controls, trajectories, and a functional; GRAPE.jl then performs the optimization by dispatching to suitable propagators and AD tools (Goerz et al., 3 Nov 2025). State representations can be wavefunctions or vectorized density matrices, and the types are user-extensible. The same is true for Hamiltonian or generator types, which may encode arbitrary control dependence, including non-linear couplings.

Propagation is delegated either to QuantumPropagators.jl, for efficient piecewise-constant propagation and auxiliary-matrix derivatives, or to DifferentialEquations.jl for more general dynamics (Goerz et al., 3 Nov 2025). Optimization itself is decoupled from propagation and gradient generation. In original GRAPE, updates were fixed-step gradient steps, but in GRAPE.jl the default is to pass gradients to a quasi-Newton optimizer, most notably L-BFGS-B via the LBFGSB.jl wrapper. This supports bound constraints on amplitudes and, according to the package description, dramatically improves convergence speed and robustness compared to fixed-step gradient ascent, consistent with the findings of Fouquières et al. (Goerz et al., 3 Nov 2025).

The package is explicitly designed around Julia’s language features. Multiple dispatch is used to select specialized propagation and gradient routines based on the types of states, generators, and controls; JIT compilation produces specialized kernels once types and control structures are fixed; and the type system plus metaprogramming are used to encode structural information such as matrix size, sparsity, or symmetry (Goerz et al., 3 Nov 2025). The consequence is a framework intended to combine high-level composability with numerical performance comparable to low-level Fortran or C implementations.

4. Supported problem classes and optimization workflows

GRAPE.jl inherits the generality of QuantumControl.jl’s trajectory framework and therefore spans closed and open quantum systems (Goerz et al., 3 Nov 2025). In the abstract formulation used by the package, $J$5 may be a pure state, a vectorized density matrix, or another vector representation, and the generator may be a Hamiltonian or a Liouvillian. This is sufficient for Lindblad-type or more general Markovian master equations, so the package supports optimization of unitary operations under dissipative dynamics (Goerz et al., 3 Nov 2025).

For quantum-computing applications, the package supports both gate optimization and state transfer. In gate optimization, rather than propagating the full unitary, GRAPE.jl can propagate multiple trajectories corresponding to logical basis states and compare the induced basis-state mapping to the target gate; this is described as more memory- and time-efficient than propagating the full propagator (Goerz et al., 3 Nov 2025). For state transfer, the objective may be an overlap-based functional such as

$J$6

or a more general state-map objective (Goerz et al., 3 Nov 2025).

A significant generalization over many classical GRAPE implementations is support for arbitrary dependence of the Hamiltonian on the controls. The package description explicitly states that $J$7 can be any function of the controls, not just a linear form, which is important for realistic device models and auxiliary model constructions (Goerz et al., 3 Nov 2025). Ensemble optimization is implemented through multiple trajectories with different initial states, parameter values, or dynamical models, and this is used to build robustness against frequency detunings, noise in couplings, and calibration uncertainties (Goerz et al., 3 Nov 2025).

The practical workflow described for GRAPE.jl proceeds through a standard sequence. One defines the system, including Hilbert-space dimension, basis, drift Hamiltonians, and control-dependent terms; chooses the controls and time grid; defines trajectories for gate optimization or ensemble robustness; specifies the objective functional; runs the optimization with a chosen propagator, optimizer, and control bounds; and then performs post-processing to examine fidelities and possibly refine the discretization or robustness specification (Goerz et al., 3 Nov 2025). The stated range of applications includes NMR spin systems, superconducting qubits, trapped atoms, and systems with spatial degrees of freedom, including the tractor atom interferometry problem in Dash et al. (Goerz et al., 3 Nov 2025).

5. Functionals, benchmarks, and comparative position

GRAPE.jl is designed to accommodate a broad family of functionals. In the theoretical framework of (Goerz et al., 2022), these include state-based functionals, overlap-based functionals, gate-based functionals, state-dependent running costs, and non-analytic functionals derived from gate geometry. Concrete examples include the standard square-modulus gate functional,

$J$8

the perfect-entangler functional in the Weyl chamber built from Makhlin invariants and a population-loss penalty, and a direct concurrence-based perfect-entangler objective (Goerz et al., 2022). The importance of this design is that a user may write a Julia function that computes the functional from $J$9, overlaps, or final states, and GRAPE.jl can differentiate that top layer with Zygote while leaving the propagation gradient analytic.

In the software-announcement paper, GRAPE.jl is positioned relative to a wide range of existing implementations: SIMPSON, Spinach, pulse-finder, QuTiP, C3, QuOCS, QuanEstimation, QDYN, and full-AD frameworks based on TensorFlow, PyTorch, or JAX (Goerz et al., 3 Nov 2025). The stated comparison emphasizes four points. First, Julia plus optimized numerical libraries can deliver performance comparable to Fortran/C codes such as QDYN. Second, the use of QuantumPropagators.jl and auxiliary-matrix methods yields machine-precision gradients with no numerical overhead beyond the traditional GRAPE scheme. Third, semi-automatic differentiation avoids the runtime and memory penalties of full AD through time evolution. Fourth, the architecture is parallelizable over trajectories because different trajectory indices $[0,T]$0 correspond to independent solves (Goerz et al., 3 Nov 2025).

The most explicit quantitative evidence is provided in (Goerz et al., 2022) for two-transmon gate optimization. There, semi-AD using a Chebyshev propagator is reported to have indistinguishable runtimes from direct analytic GRAPE for the square-modulus functional. Peak RAM for semi-AD is reported as approximately 100 MB and essentially independent of Hilbert-space size and time steps in the tested range, while full AD with a generic ODE solver is reported at about 700 MB, and full AD with a Chebyshev propagator shows RAM increasing super-linearly with $[0,T]$1, reaching 8 GB at $[0,T]$2 (Goerz et al., 2022). Runtime-wise, the full-AD Chebyshev implementation is described as about $[0,T]$3 slower per time step than semi-AD, while semi-AD adds effectively no overhead over purely analytic GRAPE (Goerz et al., 2022).

These comparisons also clarify the package’s intended niche. GRAPE.jl is not presented merely as a generic AD-based control framework; rather, it is a specialized GRAPE implementation that preserves the classical efficiency profile while extending the admissible functionals. The package announcement states that this addresses common shortcomings in existing implementations, notably rigid functional forms, lack of robust ensemble support, overhead from generic AD tools, and limited extensibility to custom data structures and dynamics (Goerz et al., 3 Nov 2025).

6. Practical considerations, limitations, and broader role

Several practical constraints follow directly from the numerical design. The time step must be small enough to resolve the fastest relevant dynamics and control features, but not smaller than necessary, because the piecewise-constant approximation increases the control dimension with the number of time slices (Goerz et al., 3 Nov 2025). The package description notes that the approximation is accurate when $[0,T]$4 is small compared to dynamical timescales and control bandwidth. A plausible implication is that hardware-facing pulse parameterizations may require additional filtering or reparameterization after optimization when the experimental controls are not naturally piecewise constant.

Initialization and functional design remain consequential. The package description recommends shaped or heuristic initial guesses when possible, since random initial pulses may converge more slowly, and analytic derivatives of the functional with respect to final states are preferable whenever available in order to avoid AD overhead (Goerz et al., 3 Nov 2025). For custom functionals, the semi-automatic differentiation interface should be used together with explicit gradient checks.

The principal scaling limitations are the standard ones for GRAPE. Computational cost grows with Hilbert-space dimension and with the number of trajectories, while open-system formulations based on density matrices scale as $[0,T]$5 with Hilbert-space dimension $[0,T]$6 (Goerz et al., 3 Nov 2025). The package also notes that poorly conditioned or highly non-linear Hamiltonian dependence can challenge both propagation and optimization, especially in the presence of non-smooth behavior or stiff dynamics (Goerz et al., 3 Nov 2025).

Two recurrent misconceptions are explicitly contradicted by the package design. One is that piecewise-constant GRAPE necessarily assumes a linear-control Hamiltonian; GRAPE.jl instead allows arbitrary dependence of $[0,T]$7 on the controls (Goerz et al., 3 Nov 2025). The other is that support for arbitrary functionals requires differentiating the complete propagation with a generic AD system; GRAPE.jl instead differentiates only the functional layer with AD and keeps the control-to-dynamics derivative analytic (Goerz et al., 2022). In that sense, the package sits between hand-coded GRAPE and end-to-end differentiable simulation frameworks.

Future directions are also stated. The current focus is on gradient-based methods with quasi-Newton updates, but extension to full Hessian-based methods, such as modified Newton-Raphson GRAPE with exact Hessians as in Goodwin and Kuprov, is described as being within scope for future versions (Goerz et al., 3 Nov 2025). The broader workflow context is experimental quantum technology: a device-specific model is built, GRAPE.jl is used offline to design pulses, pulses are implemented on hardware, experimental data feed back into model refinement, and the optimization is rerun with updated ensemble constraints or calibration uncertainties (Goerz et al., 3 Nov 2025). This positions GRAPE.jl as a foundational optimization engine for iterative design and calibration workflows in quantum computing, quantum sensing, quantum simulation, and spectroscopy.