Control Methods
All optimizations in the QuantumControl
package are done by calling QuantumControl.optimize
, or preferably the high-level wrapper @optimize_or_load
. The actual control methods are implemented in separate packages. The module implementing a particular method should be passed to optimize
as the method
keyword argument.
QuantumControl.optimize
— MethodOptimize a quantum control problem.
result = optimize(
problem;
method, # mandatory keyword argument
check=true,
callback=nothing,
print_iters=true,
kwargs...
)
optimizes towards a solution of given problem
with the given method
, which should be a Module
implementing the method, e.g.,
using Krotov
result = optimize(problem; method=Krotov)
If check
is true (default), the initial_state
and generator
of each trajectory is checked with check_state
and check_generator
. Any other keyword argument temporarily overrides the corresponding keyword argument in problem
. These arguments are available to the optimizer, see each optimization package's documentation for details.
The callback
can be given as a function to be called after each iteration in order to analyze the progress of the optimization or to modify the state of the optimizer or the current controls. The signature of callback
is method-specific, but callbacks should receive a workspace objects as the first parameter as the first argument, the iteration number as the second parameter, and then additional method-specific parameters.
The callback
function may return a tuple of values, and an optimization method should store these values fore each iteration in a records
field in their Result
object. The callback
should be called once with an iteration number of 0
before the first iteration. The callback
can also be given as a tuple of vector of functions, which are automatically combined via chain_callbacks
.
If print_iters
is true
(default), an automatic callback
is created via the method-specific make_print_iters
to print the progress of the optimization after each iteration. This automatic callback runs after any manually given callback
.
All remaining keyword argument are method-specific. To obtain the documentation for which options a particular method uses, run, e.g.,
? optimize(problem, ::Val{:Krotov})
where :Krotov
is the name of the module implementing the method. The above is also the method signature that a Module
wishing to implement a control method must define.
The returned result
object is specific to the optimization method, but should be a subtype of QuantumControl.AbstractOptimizationResult
.
The following methods of optimal control are implemented by packages in the JuliaQuantumControl organization:
Krotov's Method
See the documentation of the Krotov
package for more details.
QuantumControl.optimize
— Methodusing Krotov
result = optimize(problem; method=Krotov, kwargs...)
optimizes the given control problem
using Krotov's method, by minimizing the functional
\[J(\{ϵ_l(t)\}) = J_T(\{|Ψ_k(T)⟩\}) + ∑_l \int_{0}^{T} \frac{λ_{a,l}}{S_l(t)} [ϵ_l(t) - ϵ_l^{(0)}(t)]^2 \, dt\,,\]
cf. the general form of a quantum control functional. The "reference field" $ϵ_l^{(0)}(t)$ is the guess control for that particular iteration. The above functional implies a first-order update equation
\[Δϵ_l(t) = \frac{S_l(t)}{λ_{a,l}} \Im ∑_k \left[ \Big\langle \chi_k^{(0)}(t) \Big\vert \frac{\partial \hat{H}_k}{\partial ϵ_l(t)} \Big\vert \Psi_k(t) \Big\rangle \right]\,,\]
where $|\chi^{(0)}_k(t)⟩$ is the state backward-propagated under $Ĥ_k^{\dagger}(\{ϵ_l^{(0)}(t)\})$ with the boundary condition $|\chi_k(T)⟩ = \partial J_T / \partial ⟨Ψ_k^{(0)}(T)|$ and $Ĥ_k$ is the generator
of the $k$'th trajectory.
Note that the particular control-dependent running cost in the above functional is required to obtain the given Krotov update equation. Other running costs, or state-dependent running costs are not supported in this implementation of Krotov's method (even though some running costs are mathematically compatible with Krotov's method).
Returns a KrotovResult
.
Keyword arguments that control the optimization are taken from the keyword arguments used in the instantiation of problem
; any of these can be overridden with explicit keyword arguments to optimize
.
Required problem keyword arguments
J_T
: A functionJ_T(Ψ, trajectories)
that evaluates the final time functional from a listΨ
of forward-propagated states andproblem.trajectories
. The functionJ_T
may also take a keyword argumenttau
. If it does, a vector containing the complex overlaps of the target states (target_state
property of each trajectory inproblem.trajectories
) with the propagated states will be passed toJ_T
.
Recommended problem keyword arguments
lambda_a=1.0
: The inverse Krotov step width λₐ for every pulse.update_shape=(t->1.0)
: A functionS(t)
for the "update shape" that scales the update for every pulse.
If different controls require different lambda_a
or update_shape
, a dict pulse_options
must be given instead of a global lambda_a
and update_shape
; see below.
Optional problem keyword arguments
The following keyword arguments are supported (with default values):
pulse_options
: A dictionary that maps every control (as obtained byget_controls
from theproblem.trajectories
) to the following dict::lambda_a
: The value for inverse Krotov step width λₐ.:update_shape
: A functionS(t)
for the "update shape" that scales the Krotov pulse update.
This overrides the global
lambda_a
andupdate_shape
arguments.chi
: A functionchi(Ψ, trajectories)
that receives a listΨ
of the forward propagated states and returns a vector of states $|χₖ⟩ = -∂J_T/∂⟨Ψₖ|$. If not given, it will be automatically determined fromJ_T
viamake_chi
with the default parameters. Similarly toJ_T
, ifchi
accepts a keyword argumenttau
, it will be passed a vector of complex overlaps.sigma=nothing
: A function that calculates the second-order contribution. If not given, the first-order Krotov method is used.iter_start=0
: The initial iteration number.iter_stop=5000
: The maximum iteration number.prop_method
: The propagation method to use for each trajectory; see below.print_iters=true
: Whether to print information after each iteration.store_iter_info=Set()
: Which fields fromprint_iters
to store inresult.records
. A subset ofSet(["iter.", "J_T", "∫gₐ(t)dt", "J", "ΔJ_T", "ΔJ", "secs"])
.callback
: A function (or tuple of functions) that receives the Krotov workspace, the iteration number, the list of updated pulses, and the list of guess pulses as positional arguments. The function may return a tuple of values which are stored in theKrotovResult
objectresult.records
. The function can also mutate any of its arguments, in particular the updated pulses. This may be used, e.g., to apply a spectral filter to the updated pulses or to perform similar manipulations. Note thatprint_iters=true
(default) adds an automatic callback to print information after each iteration. Withstore_iter_info
, that callback automatically stores a subset of the printed information.check_convergence
: A function to check whether convergence has been reached. Receives aKrotovResult
objectresult
, and should setresult.converged
totrue
andresult.message
to an appropriate string in case of convergence. Multiple convergence checks can be performed by chaining functions with∘
. The convergence check is performed after anycallback
.verbose=false
: Iftrue
, print information during initialization.rethrow_exceptions
: By default, any exception ends the optimization but still returns aKrotovResult
that captures the message associated with the exception. This is to avoid losing results from a long-running optimization when an exception occurs in a later iteration. Ifrethrow_exceptions=true
, instead of capturing the exception, it will be thrown normally.
Trajectory propagation
Krotov's method involves the forward and backward propagation for every Trajectory
in the problem
. The keyword arguments for each propagation (see propagate
) are determined from any properties of each Trajectory
that have a prop_
prefix, cf. init_prop_trajectory
.
In situations where different parameters are required for the forward and backward propagation, instead of the prop_
prefix, the fw_prop_
and bw_prop_
prefixes can be used, respectively. These override any setting with the prop_
prefix. This applies both to the properties of each Trajectory
and the problem keyword arguments.
Note that the propagation method for each propagation must be specified. In most cases, it is sufficient (and recommended) to pass a global prop_method
problem keyword argument.
GRAPE
The Gradient Ascent Pulse Engineering (GRAPE) method is implemented in the GRAPE
package. See the GRAPE
documentation for details.
QuantumControl.optimize
— Methodusing GRAPE
result = optimize(problem; method=GRAPE, kwargs...)
optimizes the given control problem
via the GRAPE method, by minimizing the functional
\[J(\{ϵ_{nl}\}) = J_T(\{|ϕ_k(T)⟩\}) + λ_a J_a(\{ϵ_{nl}\})\]
where the final time functional $J_T$ depends explicitly on the forward-propagated states and the running cost $J_a$ depends explicitly on pulse values $ϵ_{nl}$ of the l'th control discretized on the n'th interval of the time grid.
Returns a GrapeResult
.
Keyword arguments that control the optimization are taken from the keyword arguments used in the instantiation of problem
; any of these can be overridden with explicit keyword arguments to optimize
.
Required problem keyword arguments
J_T
: A functionJ_T(Ψ, trajectories)
that evaluates the final time functional from a listΨ
of forward-propagated states andproblem.trajectories
. The functionJ_T
may also take a keyword argumenttau
. If it does, a vector containing the complex overlaps of the target states (target_state
property of each trajectory inproblem.trajectories
) with the propagated states will be passed toJ_T
.
Optional problem keyword arguments
chi
: A functionchi(Ψ, trajectories)
that receives a listΨ
of the forward propagated states and returns a vector of states $|χₖ⟩ = -∂J_T/∂⟨Ψₖ|$. If not given, it will be automatically determined fromJ_T
viamake_chi
with the default parameters. Similarly toJ_T
, ifchi
accepts a keyword argumenttau
, it will be passed a vector of complex overlaps.J_a
: A functionJ_a(pulsevals, tlist)
that evaluates running costs over the pulse values, wherepulsevals
are the vectorized values $ϵ_{nl}$, wheren
are in indices of the time intervals andl
are the indices over the controls, i.e.,[ϵ₁₁, ϵ₂₁, …, ϵ₁₂, ϵ₂₂, …]
(the pulse values for each control are contiguous). If not given, the optimization will not include a running cost.gradient_method=:gradgen
: One of:gradgen
(default) or:taylor
. Withgradient_method=:gradgen
, the gradient is calculated using QuantumGradientGenerators. Withgradient_method=:taylor
, it is evaluated via a Taylor series, see Eq. (20) in Kuprov and Rogers, J. Chem. Phys. 131, 234108 (2009) [17].taylor_grad_max_order=100
: If given withgradient_method=:taylor
, the maximum number of terms in the Taylor series. Iftaylor_grad_check_convergence=true
(default), if the Taylor series does not convergence within the given number of terms, throw an an error. Withtaylor_grad_check_convergence=true
, this is the exact order of the Taylor series.taylor_grad_tolerance=1e-16
: If given withgradient_method=:taylor
andtaylor_grad_check_convergence=true
, stop the Taylor series when the norm of the term falls below the given tolerance. Ignored iftaylor_grad_check_convergence=false
.taylor_grad_check_convergence=true
: If given astrue
(default), check the convergence after each term in the Taylor series an stop as soon as the norm of the term drops below the given number. Iffalse
, stop after exactlytaylor_grad_max_order
terms.lambda_a=1
: A weight for the running costJ_a
.grad_J_a
: A function to calculate the gradient ofJ_a
. If not given, it will be automatically determined. Seemake_grad_J_a
for the required interface.upper_bound
: An upper bound for the value of any optimized control. Time-dependent upper bounds can be specified viapulse_options
.lower_bound
: A lower bound for the value of any optimized control. Time-dependent lower bounds can be specified viapulse_options
.pulse_options
: A dictionary that maps every control (as obtained byget_controls
from theproblem.trajectories
) to a dict with the following possible keys::upper_bounds
: A vector of upper bound values, one for each intervals of the time grid. Values ofInf
indicate an unconstrained upper bound for that time interval, respectively the globalupper_bound
, if given.:lower_bounds
: A vector of lower bound values. Values of-Inf
indicate an unconstrained lower bound for that time interval,
print_iters=true
: Whether to print information after each iteration.print_iter_info=["iter.", "J_T", "|∇J|", "|Δϵ|", "ΔJ", "FG(F)", "secs"]
: Which fields to print ifprint_iters=true
. If given, must be a list of header labels (strings), which can be any of the following:"iter."
: The iteration number"J_T"
: The value of the final-time functional for the dynamics under the optimized pulses"J_a"
: The value of the pulse-dependent running cost for the optimized pulses"λ_a⋅J_a"
: The total contribution ofJ_a
to the full functionalJ
"J"
: The value of the optimization functional for the optimized pulses"ǁ∇J_Tǁ"
: The ℓ²-norm of the current gradient of the final-time functional. Note that this is usually the gradient of the optimize pulse, not the guess pulse."ǁ∇J_aǁ"
: The ℓ²-norm of the the current gradient of the pulse-dependent running cost. For comparison with"ǁ∇J_Tǁ"
."λ_aǁ∇J_aǁ"
: The ℓ²-norm of the the current gradient of the complete pulse-dependent running cost term. For comparison with"ǁ∇J_Tǁ"
."ǁ∇Jǁ"
: The norm of the guess pulse gradient. Note that the guess pulse gradient is not the same the current gradient."ǁΔϵǁ"
: The ℓ²-norm of the pulse update"ǁϵǁ"
: The ℓ²-norm of optimized pulse values"max|Δϵ|"
The maximum value of the pulse update (infinity norm)"max|ϵ|"
: The maximum value of the pulse values (infinity norm)"ǁΔϵǁ/ǁϵǁ"
: The ratio of the pulse update tothe optimized pulse values"∫Δϵ²dt"
: The L²-norm of the pulse update, summed over all pulses. A convergence measure comparable (proportional) to the running cost in Krotov's method"ǁsǁ"
: The norm of the search direction. Should beǁΔϵǁ
scaled by the step withα
."∠°"
: The angle (in degrees) between the negative gradient-∇J
and the search directions
."α"
: The step width as determined by the line search (Δϵ = α⋅s
)"ΔJ_T"
: The change in the final time functional relative to the previous iteration"ΔJ_a"
: The change in the control-dependent running cost relative to the previous iteration"λ_a⋅ΔJ_a"
: The change in the control-dependent running cost term relative to the previous iteration."ΔJ"
: The change in the total optimization functional relative to the previous iteration."FG(F)"
: The number of functional/gradient evaluation (FG), or pure functional (F) evaluations"secs"
: The number of seconds of wallclock time spent on the iteration.store_iter_info=[]
: Which fields to store inresult.records
, given as
a list of header labels, see
print_iter_info
.callback
: A function (or tuple of functions) that receives the GRAPE workspace and the iteration number. The function may return a tuple of values which are stored in theGrapeResult
objectresult.records
. The function can also mutate the workspace, in particular the updatedpulsevals
. This may be used, e.g., to apply a spectral filter to the updated pulses or to perform similar manipulations. Note thatprint_iters=true
(default) adds an automatic callback to print information after each iteration. Withstore_iter_info
, that callback automatically stores a subset of the available information.check_convergence
: A function to check whether convergence has been reached. Receives aGrapeResult
objectresult
, and should setresult.converged
totrue
andresult.message
to an appropriate string in case of convergence. Multiple convergence checks can be performed by chaining functions with∘
. The convergence check is performed after anycallback
.x_tol
: Parameter for Optim.jlf_tol
: Parameter for Optim.jlg_tol
: Parameter for Optim.jlshow_trace
: Parameter for Optim.jlextended_trace
: Parameter for Optim.jlshow_every
: Parameter for Optim.jlallow_f_increases
: Parameter for Optim.jloptimizer
: An optional Optim.jl optimizer (Optim.AbstractOptimizer
instance). If not given, an L-BFGS-B optimizer will be used.prop_method
: The propagation method to use for each trajectory, see below.verbose=false
: Iftrue
, print information during initializationrethrow_exceptions
: By default, any exception ends the optimization, but still returns aGrapeResult
that captures the message associated with the exception. This is to avoid losing results from a long-running optimization when an exception occurs in a later iteration. Ifrethrow_exceptions=true
, instead of capturing the exception, it will be thrown normally.
Trajectory propagation
GRAPE may involve three types of propagation:
- A forward propagation for every
Trajectory
in theproblem
- A backward propagation for every trajectory
- A backward propagation of a gradient generator for every trajectory.
The keyword arguments for each propagation (see propagate
) are determined from any properties of each Trajectory
that have a prop_
prefix, cf. init_prop_trajectory
.
In situations where different parameters are required for the forward and backward propagation, instead of the prop_
prefix, the fw_prop_
and bw_prop_
prefix can be used, respectively. These override any setting with the prop_
prefix. Similarly, properties for the backward propagation of the gradient generators can be set with properties that have a grad_prop_
prefix. These prefixes apply both to the properties of each Trajectory
and the problem keyword arguments.
Note that the propagation method for each propagation must be specified. In most cases, it is sufficient (and recommended) to pass a global prop_method
problem keyword argument.