Control optimization¶

The Hamiltonian of a controlled system can be written as \begin{align} H = H_0(\textbf{x})+\sum_{k=1}^K u_k(t) H_k, \end{align}

where \(H_0(\textbf{x})\) is the free evolution Hamiltonian with unknown parameters \(\textbf{x}\) and \(H_k\) represents the \(k\)th control Hamiltonian with \(u_k\) the corresponding control coefficients. In QuanEstimation, different algorithms are invoked to update the control coefficients. The control optimization algorithms are the gradient ascent pulse engineering (GRAPE) [1,2,3], GRAPE algorithm based on the automatic differentiation (auto-GRAPE) [4], particle swarm optimization (PSO) [5], and differential evolution (DE) [6]. The codes for control optimization are

PythonJulia

control = ControlOpt(savefile=False, method="auto-GRAPE", **kwargs)
control.dynamics(tspan, rho0, H0, dH, Hc, decay=[], ctrl_bound=[],
                 dyn_method="expm")

QFIMCFIMHCRB

control.QFIM(W=[], LDtype="SLD")

control.CFIM(M=[], W=[])

control.HCRB(W=[])

In QuanEstimation, the optimization is executed in Julia and the results are saved to files. The variable savefile indicates whether to save all the control coefficients and its default value is False which means the control coefficients for the final episode and the values of the objective function in all episodes will be saved in "controls.npy" and "f.csv", respectively. If set True then the control coefficients and the values of the objective function in all episodes will be saved during the training. The package contains four control optimization algorithms which can be set via method, **kwargs is the corresponding keyword and default value.

After calling control = ControlOpt(), the dynamics parameters should be input. Here tspan is the time length for the evolution and rho0 represents the density matrix of the initial state. H0 and dH are the free Hamiltonian and its derivatives with respect to the unknown parameters to be estimated. H0 is a matrix when the free Hamiltonian is time-independent and a list of matrices with the length equal to tspan when it is time-dependent. dH should be input as \([\partial_a{H_0}, \partial_b{H_0}, \cdots]\). Hc is a list representing the control Hamiltonians. decay contains decay operators \((\Gamma_1, \Gamma_2, \cdots)\) and the corresponding decay rates \((\gamma_1, \gamma_2, \cdots)\), its input rule is decay=[[\(\Gamma_1\), \(\gamma_1\)], [\(\Gamma_2\), \(\gamma_2\)],...]. The default value for decay is empty which means the dynamics is unitary. The package can be used to optimize bounded control problems by setting lower and upper bounds of the control coefficients via ctrl_bound, which is an array with two elements representing the lower and upper bound of each control coefficient, respectively. The default value of ctrl_bound=[] which means the control coefficients are in the regime \([-\infty,\infty]\). dyn_method="expm" represents the method for solving the dynamics is matrix exponential, it can also be set as dyn_method="ode" which means the dynamics (differential equation) is directly solved with the ODE solvers.

The objective functions for control optimization can be set as QFI \(\left[\mathrm{Tr}(W \mathcal{F}^{-1})\right]\), CFI \(\left[\mathrm{Tr}(W\mathcal{I}^{-1})\right]\), and HCRB, the corresponding codes for them are control.QFIM() (default), control.CFIM(), and control.HCRB(). Here \(\mathcal{F}\) and \(\mathcal{I}\) are the QFIM and CFIM, \(W\) corresponds to W is the weight matrix which defaults to the identity matrix. If the users call control.HCRB() for a single parameter, the program prints "Program terminated. In the single-parameter scenario, HCRB is equivalent to QFI. Please choose QFIM as the objective function."; use control.QFIM() instead. HCRB is only available with PSO and DE, and control.HCRB() raises ValueError for GRAPE and auto-GRAPE. LDtype in control.QFIM() represents the types of the QFIM, it can be set as LDtype="SLD" (default), LDtype="RLD", and LDtype="LLD". M in control.CFIM() represents a set of positive operator-valued measure (POVM) with default value [] which means a set of rank-one symmetric informationally complete POVM (SIC-POVM) is used in the calculation.

opt = ControlOpt(ctrl=ctrl, ctrl_bound=ctrl_bound, seed=1234)
alg = autoGRAPE(kwargs...)
dynamics = Lindblad(H0, dH, tspan, Hc; ctrl=ctrl, dyn_method=:Expm)

QFIMCFIMHCRB

obj = QFIM_obj(W=nothing, LDtype=:SLD)

obj = CFIM_obj(M=nothing, W=nothing)

obj = HCRB_obj(W=nothing)

scheme = GeneralScheme(probe=rho0, param=dynamics)
optimize!(scheme, opt; algorithm=alg, objective=obj, savefile=false)

In control optimization, a set of control coefficients (optimization variable) and its boundary value can be input via ctrl=ctrl and ctrl_bound=ctrl_bound in ControlOpt(). ctrl is a list of arrays with the length equal to control Hamiltonians, ctrl_bound is an array with two elements representing the lower and upper bound of the control coefficients, respectively. The default value of ctrl_bound=[-Inf, Inf] which means the control coefficients are in the regime \([-\infty,\infty]\). seed is the random seed which can ensure the reproducibility of results.

The package can deal with the parameterization process in the form of master equation, the dynamics parameters should be input via Lindblad(). Here tspan is the time length for the evolution and rho0 represents the density matrix of the initial state. H0 and dH are the free Hamiltonian and its derivatives with respect to the unknown parameters to be estimated. H0 is a matrix when the free Hamiltonian is time-independent and a list of matrices with the length equal to tspan when it is time-dependent. dH should be input as \([\partial_a{H_0}, \partial_b{H_0}, \cdots]\). Hc is a list representing the control Hamiltonians. decay contains decay operators \((\Gamma_1, \Gamma_2, \cdots)\) and the corresponding decay rates \((\gamma_1, \gamma_2, \cdots)\) with the input rule decay=[[\(\Gamma_1\), \(\gamma_1\)], [\(\Gamma_2\), \(\gamma_2\)],...]. The default value decay is nothing which means the dynamics is unitary. dyn_method=:Expm represents the method for solving the dynamics is matrix exponential, it can also be set as dyn_method=:Ode which means the dynamics (differential equation) is directly solved with the ODE solvers.

The objective functions for control optimization can be set as QFI \(\left[\mathrm{Tr}(W \mathcal{F}^{-1})\right]\), CFI \(\left[\mathrm{Tr}(W\mathcal{I}^{-1})\right]\), and HCRB, the corresponding codes for them are QFIM_obj() (default), CFIM_obj(), and HCRB_obj(). Here \(\mathcal{F}\) and \(\mathcal{I}\) are the QFIM and CFIM, \(W\) corresponds to W is the weight matrix which defaults to the identity matrix. If the users call HCRB_obj() for single parameter scenario, the program will exit and print "Program terminated. In the single-parameter scenario, the HCRB is equivalent to the QFI. Please choose 'QFIM_obj' as the objective function". LDtype in QFIM_obj() represents the types of the QFIM, it can be set as LDtype=:SLD (default), LDtype=:RLD, and LDtype=:LLD. M in CFIM_obj() represents a set of positive operator-valued measure (POVM) with default value nothing which means a set of rank-one symmetric informationally complete POVM (SIC-POVM) is used in the calculation.

The variable savefile indicates whether to save all the control coefficients and its default value is false which means the control coefficients for the final episode and the values of the objective function in all episodes will be saved in "controls.dat" (JLD2, under the key "controls") and "f.csv", respectively. If set true then the control coefficients and the values of the objective function in all episodes will be saved during the training. The algorithms used for optimizing the control coefficients in QuanEstimation are GRAPE, auto-GRAPE, PSO, and DE. kwargs... contains the keywords and defaults value corresponding to the optimization algorithm which will be introduced in detail below.

GRAPE and auto-GRAPE¶

The codes for control optimization with GRAPE and auto-GRAPE are as follows

PythonJulia

control = ControlOpt(method="GRAPE", **kwargs)

control = ControlOpt(method="auto-GRAPE", **kwargs)

where kwargs is of the form

kwargs = {"Adam":True, "ctrl0":[], "max_episode":300, "epsilon":0.01, \
            "beta1":0.90, "beta2":0.99, "seed":1234}

The keywords and the default values of GRAPE and auto-GRAPE can be seen in the following table

\(~~~~~~~~~~\)**kwargs\(~~~~~~~~~~\)	\(~~~~\)default values\(~~~~\)
"Adam"	True
"ctrl0"	[ ]
"max_episode"	300
"epsilon"	0.01
"beta1"	0.90
"beta2"	0.99
"seed"	1234

Adam algorithm can be introduced to update the control coefficients when using GRAPE and auto-GRAPE for control optimization, which can be realized by setting Adam=True. In this case, the Adam parameters include learning rate, the exponential decay rate for the first moment estimates and the second moment estimates can be set by the users via epsilon, beta1, and beta2, respectively. If Adam=False, the control coefficients will update according to the learning rate "epsilon". ctrl0 is a list representing the initial guesses of control coefficients and max_episode is the number of episodes.

alg = GRAPE(Adam=true, max_episode=300, epsilon=0.01, beta1=0.90,
            beta2=0.99)

alg = autoGRAPE(Adam=true, max_episode=300, epsilon=0.01, beta1=0.90,
                beta2=0.99)

The keywords and the default values of GRAPE and auto-GRAPE can be seen in the following table

\(~~~~~~~~~~\)keywords\(~~~~~~~~~~\)	\(~~~~\)default values\(~~~~\)
"Adam"	true
"max_episode"	300
"epsilon"	0.01
"beta1"	0.90
"beta2"	0.99

Adam algorithm can be introduced to update the control coefficients when using GRAPE and auto-GRAPE for control optimization, which can be realized by setting Adam=true. In this case, the Adam parameters include learning rate, the exponential decay rate for the first moment estimates and the second moment estimates can be set by the users via epsilon, beta1, and beta2, respectively. If Adam=false, the control coefficients will update according to the learning rate "epsilon".

PSO¶

The code for control optimization with PSO is as follows

PythonJulia

control = ControlOpt(method="PSO", **kwargs)

where kwargs is of the form

kwargs = {"p_num":10, "ctrl0":[], "max_episode":[1000,100], \
          "c0":1.0, "c1":2.0, "c2":2.0, "seed":1234}

The keywords and the default values of PSO can be seen in the following table

\(~~~~~~~~~~\)**kwargs\(~~~~~~~~~~\)	\(~~~~\)default values\(~~~~\)
"p_num"	10
"ctrl0"	[ ]
"max_episode"	[1000,100]
"c0"	1.0
"c1"	2.0
"c2"	2.0
"seed"	1234

Here p_num is the number of particles. c0, c1, and c2 are the PSO parameters representing the inertia weight, cognitive learning factor, and social learning factor, respectively. max_episode accepts both integer and array with two elements. If it is an integer, for example max_episode=1000, it means the program will continuously run 1000 episodes. However, if it is an array, for example max_episode=[1000,100], the program will run 1000 episodes in total but replace control coefficients of all the particles with global best every 100 episodes. seed is the random seed which can ensure the reproducibility of the results.

alg = PSO(p_num=10, ini_particle=nothing, max_episode=[1000,100],
          c0=1.0, c1=2.0, c2=2.0)

The keywords and the default values of PSO can be seen in the following table

\(~~~~~~~~~~\)keywords\(~~~~~~~~~~\)	\(~~~~\)default values\(~~~~\)
"p_num"	10
"ini_particle"	nothing
"max_episode"	[1000,100]
"c0"	1.0
"c1"	2.0
"c2"	2.0

Here p_num is the number of particles. c0, c1, and c2 are the PSO parameters representing the inertia weight, cognitive learning factor, and social learning factor, respectively. max_episode accepts both integer and array with two elements. If it is an integer, for example max_episode=1000, it means the program will continuously run 1000 episodes. However, if it is an array, for example max_episode=[1000,100], the program will run 1000 episodes in total but replace control coefficients of all the particles with global best every 100 episodes. ini_particle is a tuple representing the initial guesses of control coefficients.

DE¶

The code for control optimization with DE is as follows

PythonJulia

control = ControlOpt(method="DE", **kwargs)

where kwargs is of the form

kwargs = {"p_num":10, "ctrl0":[], "max_episode":1000, "c":1.0, \
            "cr":0.5, "seed":1234}

The keywords and the default values of DE can be seen in the following table

\(~~~~~~~~~~\)**kwargs\(~~~~~~~~~~\)	\(~~~~\)default values\(~~~~\)
"p_num"	10
"ctrl0"	[ ]
"max_episode"	1000
"c"	1.0
"cr"	0.5
"seed"	1234

p_num and max_episode represent the number of populations and episodes. c and cr are the mutation constant and the crossover constant.

alg = DE(p_num=10, ini_population=nothing, max_episode=1000,
         c=1.0, cr=0.5)

The keywords and the default values of DE can be seen in the following table

\(~~~~~~~~~~\)keywords\(~~~~~~~~~~\)	\(~~~~\)default values\(~~~~\)
"p_num"	10
"ini_population"	nothing
"max_episode"	1000
"c"	1.0
"cr"	0.5

ini_population is a tuple representing the initial guesses of the control coefficients, p_num is the number of populations, and max_episode is the number of episodes. c and cr are the mutation constant and the crossover constant.

Example 5.1 In this example, the free evolution Hamiltonian of a single qubit system is \(H_0=\frac{1}{2}\omega \sigma_3\) with \(\omega\) the frequency and \(\sigma_3\) a Pauli matrix. The dynamics of the system is governed by

\[\begin{align} \partial_t\rho=-i[H_0, \rho]+ \gamma_{+}\left(\sigma_{+}\rho\sigma_{-}-\frac{1}{2}\{\sigma_{-} \sigma_{+},\rho\}\right)+ \gamma_{-}\left(\sigma_{-}\rho\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-}, \rho\}\right), \end{align}\]

where \(\gamma_{+}\), \(\gamma_{-}\) are decay rates and \(\sigma_{\pm}=(\sigma_1 \pm \sigma_2)/2\). The control Hamiltonian \begin{align} H_\mathrm{c}=u_1(t)\sigma_1+u_2(t)\sigma_2+u_3(t)\sigma_3. \end{align}

Here \(\sigma_{1}\), \(\sigma_{2}\) are also Pauli matrices. The probe state is taken as \(|+\rangle\) and the measurement for CFI is \(\{|+\rangle\langle+|, |-\rangle\langle-|\}\) with \(|\pm\rangle:=\frac{1}{\sqrt{2}}(|0\rangle\pm|1\rangle)\). \(|0\rangle\) \((|1\rangle)\) is the eigenstate of \(\sigma_3\) with respect to the eigenvalue \(1\) \((-1)\).

Example 5.2 In the multiparameter scenario, the dynamics of electron and nuclear coupling in NV\(^{-}\) can be expressed as \begin{align} \partial_t\rho=-i[H_0,\rho]+\frac{\gamma}{2}(S_3\rho S_3-S^2_3\rho-\rho S^2_3) \end{align}

with \(\gamma\) the dephasing rate. And \begin{align} H_0/\hbar=DS^2_3+g_{\mathrm{S}}\vec{B}\cdot\vec{S}+g_{\mathrm{I}}\vec{B}\cdot\vec{I}+\vec{S}^ {\,\mathrm{T}}\mathcal{A}\vec{I} \end{align}

is the free evolution Hamiltonian, where \(\vec{S}=(S_1,S_2,S_3)^{\mathrm{T}}\) and \(\vec{I}=(I_1,I_2,I_3)^{\mathrm{T}}\) with \(S_i=s_i\otimes I\) and \(I_i=I\otimes \sigma_i\) \((i=1,2,3)\) the electron and nuclear operators. \(s_1, s_2, and s_3\) are spin-1 operators with

\[\begin{eqnarray} s_1 = \frac{1}{\sqrt{2}}\left(\begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{array}\right), s_2 = \frac{1}{\sqrt{2}}\left(\begin{array}{ccc} 0 & -i & 0\\ i & 0 & -i\\ 0 & i & 0 \end{array}\right), \nonumber \end{eqnarray}\]

and \(s_3=\mathrm{diag}(1,0,-1)\) and \(\sigma_i (i=1,2,3)\) is Pauli matrix. \(\mathcal{A}=\mathrm{diag} (A_1,A_1,A_2)\) is the hyperfine tensor with \(A_1\) and \(A_2\) the axial and transverse magnetic hyperfine coupling coefficients. The coefficients \(g_{\mathrm{S}}=g_\mathrm{e}\mu_\mathrm{B}/\hbar\) and \(g_{\mathrm{I}}=g_\mathrm{n}\mu_\mathrm{n}/\hbar\), where \(g_\mathrm{e}\) (\(g_\mathrm{n}\)) is the \(g\) factor of the electron (nuclear), \(\mu_\mathrm{B}\) (\(\mu_\mathrm{n}\)) is the Bohr (nuclear) magneton, and \(\hbar\) is the Planck's constant. \(\vec{B}\) is the magnetic field which is to be estimated.

In this case, the initial state is taken as \(\frac{1}{\sqrt{2}}(|1\rangle+|\!-\!1\rangle)\otimes |\!\!\uparrow\rangle\), where \(\frac{1}{\sqrt{2}}(|1\rangle+|\!-\!1\rangle)\) is an electron state with \(|1\rangle\) \(\left(|-\!1\rangle\right)\) the eigenstate of \(s_3\) with respect to the eigenvalue \(1\) (\(-1\)). \(|\uparrow\rangle\) is a nuclear state and the eigenstate of \(\sigma_3\) with respect to the eigenvalue 1. \(W\) is set to be identity.

In this example, auto-GRAPE, PSO, and DE are used to search the optimal controls; GRAPE is also available.

PythonJulia

Minimum-time optimization¶

Search for the minimum time to reach a given value of the objective function.

PythonJulia

control.mintime(f, W=[], M=[], method="binary", target="QFIM",
                LDtype="SLD")

mintime(f, opt, alg, obj, dynamics; method="binary", savefile=false)

f is the given value of the objective function. Two search methods are provided: binary search and forward search from the beginning of time. They can be selected with method="binary" (default) and method="forward". target represents the objective function for searching the minimum time; the users can choose QFI \(\left[\mathrm{Tr}(WF^{-1})\right]\), CFI \(\left[\mathrm{Tr}(WI^{-1})\right]\), and HCRB. HCRB is not available as a target for GRAPE and auto-GRAPE. If target="QFIM", the types of logarithmic derivatives can be set via LDtype. M is the POVM used for the CFI objective; when M=[], the QFIM is constructed instead, so a POVM must be supplied explicitly to search on the CFI. The minimum-time search is implemented for Lindblad dynamics.

Bibliography¶

[1] N. Khaneja, T. Reiss, C. Hehlet, T. Schulte-Herbruggen, and S. J. Glaser, Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms, J. Magn. Reson. 172, 296 (2005).

[2] J. Liu and H. Yuan, Quantum parameter estimation with optimal control, Phys. Rev. A 96, 012117 (2017).

[3] J. Liu and H. Yuan, Control-enhanced multiparameter quantum estimation, Phys. Rev. A 96, 042114 (2017).

[4] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: a survey, J. Mach. Learn. Res. 18, 1-43 (2018).

[5] J. Kennedy and R. Eberhar, Particle swarm optimization, Proc. 1995 IEEE International Conference on Neural Networks 4, 1942-1948 (1995).

[6] R. Storn and K. Price, Differential Evolution-A Simple and Efficient Heuristic for global Optimization over Continuous Spaces, J. Global Optim. 11, 341 (1997).

See also: Parameterization process — construct the controlled Hamiltonian and define dynamics parameters. Quantum metrological tools — the objective functions being optimized. Output files — how to load the saved controls.npy results.

Control optimization¶

GRAPE and auto-GRAPE¶

PSO¶

DE¶

input the dynamics data¶

build the scheme and run the control optimization problem¶

Minimum-time optimization¶

Bibliography¶