You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Piece-wise constant control functions and total-variation penalty term (#27)
* Implemented piece-wise constant B-spline ctrl functions. New configuration option for control types:
control_segment0 = spline0, <nsplines>
* Add total variation regularization to prevent noisy. New configuration option:
optim_penalty_variation = <value> for the c++ code and
gamma_variation = <value> for the python interface
* Documenting the zeroth order B-spline basis function in the user's guide.
* Setting default carrier wave frequencies to zero, if piecewise constant controls are used through the python interface
// Define the controllable segments for each oscillator and the type of parameterization. Multiple segments can be listed behind each other, with corresponding starting and finish times.
48
48
// Format: <controltype>, <number of basis functions> [, <tstart>, <tstop>]
49
+
// Available control types: "spline" for 2nd order Bspline basis functions (recommended), "spline0" for piecewise constant control parameterization (aka 0th order Bspline basis functions)
The time-dependent rotating-frame control functions $d^k(\vec{\alpha}^k,t)$ are parameterized using $N_s^k$piecewise quadratic B-spline basis functions $B_s(t)$ acting as envelope for $N_f^k$ carrier waves:
219
-
\begin{align}
218
+
The time-dependent rotating-frame control functions $d^k(\vec{\alpha}^k,t)$ are parameterized using $N_s^k$ basis functions $B_s(t)$ acting as envelope for $N_f^k$ carrier waves:
By default, the basis functions are piecewise quadratic (2nd order) B-spline polynomials, centered on an equally spaced grid in time. To instead use a piecewise constant (0th order) B-spline basis, see Section~\ref{subsec:bspline-0}.
222
223
The amplitudes $\alpha_{s,f}^{k(1)}, \alpha_{s,f}^{k(2)} \in\R$ are the control
223
224
parameters (\textit{design} variables) that Quandary can optimize in order to realize a
224
225
desired system behavior, giving a total number of $2\sum_kN_s^kN_f^k$ real-valued optimization variables. (Note that the number of carrier wave frequencies $N_f^k$ as well as the number of spline basis functions $N_s^k$ can be different for each subsystem $k$.) $\Omega_k^f \in\R$ denote the carrier wave frequencies in the rotating frame which can be choosen to trigger certain system frequencies. The corresponding Lab-frame carrier frequencies become $\omega_k^r + \Omega_k^f$. Those frequencies can be chosen to match the transition frequencies in the lab-frame system Hamiltonian. For example, when $\xi_{kl} << \xi_k$, the transition frequencies satisfy $\omega_k - n\xi_k$. Thus by choosing $\Omega_k^f = \omega_k-\omega_k^r - n \xi_k$, one triggers transition between energy levels $n$ and $n+1$ in subsystem $k$. Choosing effective carrier wave frequencies is quite important for optimization performance. We recommend to have a look at \cite{petersson2021optimal} for details on how to choose them.
@@ -285,6 +286,17 @@ \subsubsection{Storage of the control parameters}
285
286
q(t) = \sum_f \sum_s \alpha_{f,s} \sin(\Omega_f t + b_f) B_s(t)
286
287
\end{align}
287
288
289
+
\subsubsection{Zeroth order B-spline basis functions (piecewise constant controls)}\label{subsec:bspline-0}
290
+
A piecewise continuous envelope function can be generated by using zeroth order B-spline basis functions. When the carrier wave frequency is set to zero, this results in a control function that is piecewise constant in the rotating frame. For example, to use the zeroth order basis functions for controlling sub-system number 0 with 50 constant control segments, use the configuration option:
291
+
\begin{verbatim}
292
+
control_segments0 = spline0, 50
293
+
\end{verbatim}
294
+
When optimizing with zeroth order B-spline control functions, strong variations between consequtive control amplitudes can be avoided by enabling the total variation penalty term through the command
295
+
\begin{verbatim}
296
+
optim_penalty_variation= 1.0
297
+
\end{verbatim}
298
+
Compare Section \ref{sec:penalty}.
299
+
288
300
289
301
\subsection{Interfacing to Python environemnt}
290
302
You can use the Python interface for Quandary to simulate and optimize from within a python environment (version $\geq$ 3). It eases the use of Quandary, and adds some additional functionality, such as automatic computation of the required number of time steps, automatic choice of the carrier frequencies, and it allows for custom Hamiltonian models to be used (system and control Hamiltonian operators $H_d$ and $H_c$). A good place to start is to have a look into the example \texttt{example\_swap02.py}. This test case optimizes for a 3-level SWAP02 gate that swaps the state of the zero and the second energy level of a 3-level qudit.
@@ -310,10 +322,10 @@ \section{The Optimal Control Problem} \label{sec:optim}
310
322
In the most general form, Quandary can solve the following optimization problem
where $\rho_i(T)$ denotes one or multiple quantum states evaluated at a final time $T>0$, which solve either Lindblad's master equation \eqref{mastereq} or Schroedinger's equation \eqref{eq:schroedinger} in the rotating frame for initial conditions $\rho_i(0)$, as specified in Section \ref{subsec:initcond}, $i=1,\dots, n_{init}$. The first term in \eqref{eq:minproblem} minimizes an objective function $J$ (see Section \ref{sec:objectivefunctionals}) that quantifies the discrepancy between the realized states $\rho_i(T)$ at final time $T$ driven by the current control $\boldsymbol{\alpha}$ and the desired target $\rho^{target}_i$, see Section \ref{sec:targets}.
316
-
The second term is a standart Tikhonov regularization that can be added with parameter $\gamma_1\geq0$ in order to regularize the optimization problem (stabilize optimization convergence) by favoring solutions with small norm. Additional penalty terms, e.g. for leakage prevention, can be added, compare \ref{sec:penalty}.
328
+
The remaining terms are regularization and penalty terms that can be added to stabilize convergence, or prevent leakage, compare Section \ref{sec:penalty}
317
329
318
330
\subsection{Fidelity}\label{sec:fidelity}
319
331
@@ -459,7 +471,7 @@ \subsubsection{Basis states}
459
471
460
472
461
473
462
-
\subsubsection{Diagonal density matrices (aka all pure states)}
474
+
\subsubsection{Only the diagonal density basis matrices}
463
475
464
476
For density matrices (Lindblad solver), one can choose to propagate only those basis states that correspond to pure states of the form $\boldsymbol{e}_k\boldsymbol{e}_k^\dagger$, i.e. propagating only the $B^{kk}$ in \eqref{eq:basismats} for $k=0,\dots, N-1$, and then $n_{init}=N$. For the Schroedinger solver, this is equivalent to all basis states.
465
477
@@ -514,22 +526,36 @@ \subsubsection{Reading an initial state from file}
514
526
515
527
516
528
517
-
\subsection{Penalty terms and leakage prevention}\label{sec:penalty}
529
+
\subsection{Tikhonov regularization, penalty terms, and leakage prevention}\label{sec:penalty}
518
530
519
-
Three optional penalty terms can be added to the objective function:
531
+
In order to regularize the optimization problem (stabilize optimization convergence), a standard Tikhonov regularization term can added to the objective function.
The first penalty term can be added with $\gamma_2\geq0$ to achieve a desired behavior over the entire time-domain $0\leq t\leq T$, rather than only at the final time. If extra (non-essential levels are considered through the optimization for at least one oscillator ($n_k^e < n_k$ for at least $k$, compare Sec. \ref{sec:essentiallevels}), then this term can be used to prevent leakage to higher energy levels that are not modelled. In particular, in that case the occupation of all \textit{guard levels} are penalized through
535
+
By adding this term with a small parameter $\gamma_1 > 0$, the optimization problem will favor optimal control vectors that have a small norm. It regularizes the optimization problem since it adds a small but positive identity matrix to the Hessian of the objective function, hence ``convexifying'' the problem.
536
+
537
+
In addition to the Tikhonov regularization term, four additional penalty terms can optionally be added to the objective function if desired:
538
+
\begin{align*}
539
+
Penalty &= \frac{\gamma_2}{T} \int_0^T P\left(\{\rho_i(t)\}\right) \,\mathrm{d} t \hspace{3cm} \rightarrow\text{Leakage prevention}\\
The first penalty term can be added with $\gamma_2 > 0$ to drive the quantum system towards the desired state over the entire time-domain $0\leq t\leq T$, rather than only at the final time. If extra (non-essential levels are considered through the optimization for at least one oscillator ($n_k^e < n_k$ for at least $k$, compare Sec.~\ref{sec:essentiallevels}), then this term can be used to prevent leakage to higher energy levels that are not modelled. In particular, in that case, the occupation of all \textit{guard levels} are penalized with
524
545
\begin{align}\label{eq:leakprevention}
525
546
P(\rho(t)) = \sum_{r} \|\rho(t)_{rr} \|^2_2
526
547
\end{align}
527
-
where $r$ iterates over all indices that correspond to a guard level (i.e. last non-essential energy level) of at least one of the subsystems, and $\rho(t)_{rr}$ denotes their corresponding population.
548
+
where $r$ iterates over all indices that correspond to a guard level (i.e., the final (highest) non-essential energy level) of at least one of the subsystems, and $\rho(t)_{rr}$ denotes their corresponding population.
528
549
550
+
The second penalty term can be added with parameter $\gamma_3 > 0$ to encourage solutions whose populations vary slowly in time by penalizing the second derivative of the populations of the state.
529
551
530
-
The second penalty term can be added with parameter $\gamma_3\geq0$ to encourage solutions whose populations vary slowly in time by penalizing the second derivative of the populations of the state.
552
+
The third penalty term can be added with parameter $\gamma_4 >0$ to encourage small control pulse amplitudes by penalizing the control pulse energy. This term can be useful if hardware bounds are given for the control pulse amplitudes: Rather than include amplitude bounds on control pulse directly, which often leads to more non-convex optimization problems and convergence deterioration, one can utilize this penalty term to favor short control pulses with small amplitudes. Compare also \cite{gunther2023practical} for its usage to determine minimal gate durations.
531
553
532
-
The third penalty term can be added with parameter $\gamma_4\geq0$ to encourage small control pulse amplitudes by penalizing the control pulse energy. This term can be useful if hardware bounds are given for the control pulse amplitudes: Rather than include amplitude bounds on control pulse directly, which often leads to more non-convex optimization problems and convergence deterioration, one can utilize this penalty term to favor short control pulses with small amplitudes. Compare also \cite{gunther2023practical} for its usage to determine minimal gate durations.
554
+
The last penalty term, activated by setting $\gamma_5>0$, is used to penalize variation in control strength between consequtive B-spline coefficients. It is currently only implemented for piecewise zeroth order spline functions, see Section \ref{subsec:bspline-0}, where it is useful to prevent noisy control pulses. Referring to the control function representation in \eqref{eq:spline-ctrl}, this penalty function takes the form:
in terms of the complex-valued control parameters $\alpha_{s,f}^k = \alpha_{s,f}^{k(1)} + i \alpha_{s,f}^{k(2)}$. Penalizing the variance can significantly reduce the noise level in the optimized control functions.
533
559
534
560
Note: All regularization and penalty coefficients $\gamma_i$ should be chosen small enough so that they do not dominate the final-time objective function $J$. This might require some fine-tuning. It is recommended to always add $\gamma_1>0$, e.g. $\gamma_1 = 10^{-4}$, and add other penalties only if needed.
# Further, make sure that your quandary executable is in your $PATH variable. E.g. with
4
+
# > export PATH=/path/to/quandary/:$PATH
5
+
fromquandaryimport*
6
+
7
+
## Two qubit test case, demonstrating the use of piecewise constant control functions with total variation penalty term. CNOT gate, two levels each, no guard levels, dipole-dipole coupling 5KHz ##
Jkl= [0.005] # Dipole-Dipole coupling of qubit 0<->1
13
+
# Frequency of rotations for computational frame [GHz] per oscillator
14
+
favg=sum(freq01)/len(freq01)
15
+
rotfreq=favg*np.ones(len(freq01))
16
+
17
+
# Set the pulse duration (ns)
18
+
T=200.0
19
+
20
+
# Set up the CNOT target gate
21
+
unitary=np.identity(4)
22
+
unitary[2,2] =0.0
23
+
unitary[3,3] =0.0
24
+
unitary[2,3] =1.0
25
+
unitary[3,2] =1.0
26
+
# print("Target gate: ", unitary)
27
+
28
+
# Flag for printing out more information to screen
29
+
verbose=False
30
+
31
+
# For reproducability: Random number generator seed
32
+
rand_seed=1234
33
+
34
+
# For piecewise constant control functions, choose spline order of 0 (Default spline order would be 2, being 2nd order Bsplines). Note, the number spline basis functions for piecewise constant controls has to be much larger than if you use 2nd order Bsplines. Also note that if the spline order is zero, it is recommended not to use any carrier frequencies, which is already the default.
35
+
spline_order=0
36
+
nsplines=1000
37
+
38
+
# In order get less noisy control functions, activate the penalty term for variation of the control parameters
39
+
gamma_variation=1.0
40
+
41
+
# Optionally: let controls functions start and end near zero
42
+
control_enforce_BC=True
43
+
44
+
# Set up the Quandary configuration for this test case. Make sure to pass all of the above to the corresponding fields, compare help(Quandary)!
0 commit comments