# An improvement of D-MORPH method for finding quantum optimal
control††thanks: Punlished in ”Mezdunarodnyj naucno-issledovatel’skij zurnal”
[”International Research Journal”], 2016, Vol. 6(48), p. 94–99. DOI:
10.18454/IRJ.2016.48.057 [in Russian]

Konstantin Zhdanov Department of applied mathematics and control processes,
Saint-Petersburg state university, Russia. ORCID: 0000-0003-2290-6324

(June, 2016)

###### Abstract

The paper examines the prominent algorithm D-MORPH to search for the optimal
control of a quantum system in order to implement desired unitary evolution of
the quantum system at the final time, and reveals new mathematical expressions
for various orders’ corrections to the algorithm, that include information
about the commutators of the system’s Hamiltonian. Inclusion of such
corrections results in faster optimal quantum control’s search with high
precision, i.e. allows saving of computational resources.

Keywords:  quantum systems, optimal control, evolution operator

##  1 Introduction

Many papers of distinguished specialists in the area of quantum control [1, 2,
3, 4, 5] use the D-MORPH method to construct optimal quantum control. The
method is based on introducing a new variable $s$ which represents the
optimization algorithm’s progress towards the minimum of an objective
function. Variable $s$ has an interesting property: while increasing $s$, the
optimization algorithm is constructing new controls that drive the system
closer and closer to the desired state associated with the minimum of the
objective function.

This optimization algorithm is approximate since it’s based on solving a
system of ordinary differential equations numerically with the help of the
popular MATLAB’s ode45 method, i. e. a variable-step Runge–Kutta’s method of
order four. As it’s been demonstrated many times [1, 2, 3, 4, 5], this
optimization algorithm allows construction of controls resulting in close-to-
the-optimal values of the objective function. One weakness of the algorithm is
the requirement to solve the system on a large interval of integration $[0,S]$
in order to achieve a high accuracy, which can be very computationally
demanding.

This paper introduces a way to improve the aforementioned algorithm — new
formulas for corrections of different orders to the original method are
derived so that they allow high-accuracy reaching of the minimum of an
objective function, which represents the distance between the actual evolution
of a quantum system and the desired evolution, with fewer integration steps
than the original method requires and so spending less time. These corrections
allow searching for optimal quantum controls to be effectively carried out on
common computers with moderate performance.

##  2 The original method and corrections to it

Let’s consider an $N$-level quantum system subjected to $n$ control fields and
formulate a task: find such control fields $\epsilon_{k}(t)$,
$k=\overline{1,n}$, that the quantum system at a moment $T$ implements the
desired evolution operator $U_{D}$ as accurately as possible. The accuracy is
expressed in terms of an objective function $J$, which represents the distance
between the actual evolution $U(T,0)$ of the system and the desired evolution
$U_{D}$.

The original algorithm is formulated as follows. A new variable $s$ is
introduced such that objective $J(s)$ and controls $\epsilon_{k}(s,t)$ both
depend on it. Then we find the derivative

$\frac{\mathrm{d}J}{\mathrm{d}s}=\sum_{k=1}^{n}\int_{0}^{T}\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)}\frac{\mathrm{d}\epsilon_{k}(s,t)}{\mathrm{d}s}\mathrm{d}t.$

To find the minimum of the objective function $J$, it’s sufficient to require
that $\nicefrac{{\mathrm{d}J}}{{\mathrm{d}s}}\leq 0$ holds. In that case, if
we increase $s$, $J$ at least won’t increase.

Let’s divide the interval $[0,T]$ into $L$ subintervals of the same length
$\Delta t$ and treat controls as piecewise constant with respect to $t$
functions

$\epsilon_{k}(s,t)=\sum_{l=1}^{L}\epsilon_{k}^{l}(s)\chi_{[t_{l-1},t_{l}]}(t),$

where $\chi_{[t_{l-1},t_{l}]}(t)$ — the indicator function of interval
$[t_{l-1},t_{l}]$,

$\chi_{[t_{l-1},t_{l}]}(t)=\begin{cases}1,&t\in[t_{l-1},t_{l}]\\\
0,&t\notin[t_{l-1},t_{l}]\end{cases}.$

Then inequality $\nicefrac{{\mathrm{d}J}}{{\mathrm{d}s}}\leq 0$ can be written
as

$\frac{\mathrm{d}J}{\mathrm{d}s}=\sum_{k=1}^{n}\sum_{l=1}^{L}\int_{t_{l-1}}^{t_{l}}\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)}\frac{\mathrm{d}\epsilon_{k}^{l}(s)}{\mathrm{d}s}\mathrm{d}t\leq
0.$

The inequality is satisfied if we let

$\frac{\mathrm{d}\epsilon_{k}^{l}(s)}{\mathrm{d}s}=-\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)},\quad t\in[t_{l-1},t_{l}],\quad
l=\overline{1,L},\quad k=\overline{1,n}.$ (1)

Therefore the original method solves the system of differential equations (1),
where $\epsilon_{k}^{l}(s)$ are unknown functions to be found. As it’s assumed
that a functional form of $J$ with respect to $\epsilon_{k}(s,t)$ is known in
advance, the right-hand side of the expression above is assumed to be known as
well. It’s worth noting that the right-hand side depends on all
$\epsilon_{k}^{l}(s)$.

In papers [1, 2, 3, 4, 5] such systems of differential equations are solved
with a variable-step Runge–Kutta’s method of order four (MATLAB’s ode45) on
$[0,S]$, where $S$ — a large number. Then an approximately optimal control is
retrieved at $s=S$ owing to the method guarantees that the objective function
at least doesn’t increase while increasing $s$. In this paper the
aforementioned original method is modified so as to get approximate optimal
control with given accuracy by solving the system on smaller interval $[0,S]$
than used in the original method. A decrease in the length of the interval of
integration should diminish the time spent by ode45 method, which should be
beneficial when using computers with moderate performance.

Since the right-hand side of system (1) depends on time while the left-hand
side doesn’t, we integrate each equation of the system with respect to $t$
from $t_{l-1}$ to $t_{l}$.

$\int_{t_{l-1}}^{t_{l}}\frac{\mathrm{d}\epsilon_{k}^{l}(s)}{\mathrm{d}s}\mathrm{d}t=-\int_{t_{l-1}}^{t_{l}}\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)}\mathrm{d}t,\quad l=\overline{1,L},\quad
k=\overline{1,n},$
$\frac{\mathrm{d}\epsilon_{k}^{l}(s)}{\mathrm{d}s}=-\frac{1}{\Delta
t}\int_{t_{l-1}}^{t_{l}}\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)}\mathrm{d}t,\quad l=\overline{1,L},\quad
k=\overline{1,n}.$

In what follows, we are going to expand the integral on the right-hand side of
the expression above. For a task of implementing the desired evolution
operator, the objective function in the original method is taken as [1, 3]

$J=\frac{1}{2}-\frac{1}{2N}\operatorname{Re}\operatorname{Tr}\left[U_{D}^{*}U(T,0)\right].$

This objective function represents the distance between the desired evolution
operator $U_{D}$ and the actual evolution operator $U(T,0)$ at time $T$. The
derivative of $J$ with respect to controls is written as [3]

$\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)}=-\frac{1}{2N}\operatorname{Im}\operatorname{Tr}\left[U_{D}^{*}U(T,t)H_{k}U(t,0)\right],$

while the integral of this expression is equal to

$\int_{t_{l-1}}^{t_{l}}\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)}\mathrm{d}t=-\frac{1}{2N}\operatorname{Im}\operatorname{Tr}\left[U_{D}^{*}\int_{t_{l-1}}^{t_{l}}U(T,t)H_{k}U(t,0)\mathrm{d}t\right]=\\\
=-\frac{1}{2N}\operatorname{Im}\operatorname{Tr}\left[U_{D}^{*}U(T,0)\int_{t_{l-1}}^{t_{l}}U^{*}(t,0)H_{k}U(t,0)\mathrm{d}t\right].$

In general, the evolution operator of a quantum system with Hamiltonian
$H(s,t)=H_{0}+\sum_{k=1}^{n}\epsilon_{k}(s,t)H_{k}$ is expressed as
($\hbar=1$) [1]

$U(T,0)=\operatorname{T}\exp\left\\{-\imath\int_{0}^{T}H(s,\tau)\mathrm{d}\tau\right\\},$

where $\operatorname{T}$ — time-ordering operator. Since the controls are
piecewise constant, the evolution operator turns into a composition of matrix
exponentials

$U(T,0)=\exp\left\\{-\imath\int_{t_{L-1}}^{T}H(s,\tau)\mathrm{d}\tau\right\\}\exp\left\\{-\imath\int_{t_{L-2}}^{t_{L-1}}H(s,\tau)\mathrm{d}\tau\right\\}\cdots\\\
\cdots\exp\left\\{-\imath\int_{0}^{t_{1}}H(s,\tau)\mathrm{d}\tau\right\\}.$

Using the following properties of the evolution operator

$U(t,0)=U(t,t_{l-1})\cdots U(t_{1},0)=U(t,t_{l-1})U(t_{l-1},0),\quad
t\in[t_{l-1},t_{l}],$ $U^{*}(t,0)=U^{*}(t_{1},0)\cdots
U^{*}(t,t_{l-1})=U^{*}(t_{l-1},0)U^{*}(t,t_{l-1}),\quad t\in[t_{l-1},t_{l}]$

we get

$\int_{t_{l-1}}^{t_{l}}U^{*}(t,0)H_{k}U(t,0)\mathrm{d}t=U^{*}(t_{l-1},0)\int_{t_{l-1}}^{t_{l}}U^{*}(t,t_{l-1})H_{k}U(t,t_{l-1})\mathrm{d}tU(t_{l-1},0).$

Let’s express the evolution operator as the matrix exponential.

$\int_{t_{l-1}}^{t_{l}}U^{*}(t,t_{l-1})H_{k}U(t,t_{l-1})\mathrm{d}t=\\\
=\int_{t_{l-1}}^{t_{l}}\exp\left[(t-t_{l-1})X\right]H_{k}\exp\left[-(t-t_{l-1})X\right]\mathrm{d}t=\int_{t_{l-1}}^{t_{l}}\operatorname{Ad}_{\exp[(t-t_{l-1})X]}H_{k}\mathrm{d}t,$

where $X=\imath\left(H_{0}+\sum_{k=1}^{n}\epsilon_{k}^{l}H_{k}\right)$,
$\operatorname{Ad}_{A}B=ABA^{-1}$. It’s well known in the general theory of
Lie algebras that

$\operatorname{Ad}_{\exp(A)}B=\exp(\operatorname{ad}_{A})\circ B,$

where $\operatorname{ad}_{A}=AB-BA$. Using this expression we get

$\int_{t_{l-1}}^{t_{l}}\operatorname{Ad}_{\exp[(t-t_{l-1})X]}H_{k}\mathrm{d}t=\int_{t_{l-1}}^{t_{l}}\exp\left(\operatorname{ad}_{(t-t_{l-1})X}\right)\circ
H_{k}\mathrm{d}t$

Let’s expand operator $\exp\left(\operatorname{ad}_{(t-t_{l-1})X}\right)$ in a
Taylor series and use identity
$\operatorname{ad}_{(t-t_{l-1})X}=(t-t_{l-1})\operatorname{ad}_{X}$

$\exp\left(\operatorname{ad}_{(t-t_{l-1})X}\right)=\sum_{j=0}^{\infty}\frac{(t-t_{l-1})^{j}}{j!}\operatorname{ad}_{X}^{j}.$

Then we integrate it

$\int_{t_{l-1}}^{t_{l}}\exp\left(\operatorname{ad}_{(t-t_{l-1})X}\right)\circ
H_{k}\mathrm{d}t=\Delta t\left(H_{k}+\frac{\Delta t}{2}[X,H_{k}]+\frac{(\Delta
t)^{2}}{3!}[X,[X,H_{k}]]+\cdots\right).$

Finally we get

$\int_{t_{l-1}}^{t_{l}}\frac{\partial
J(s)}{\partial\epsilon_{k}(s,t)}\mathrm{d}t=-\frac{\Delta
t}{2N}\operatorname{Im}\operatorname{Tr}\left[U_{D}^{*}U(T,t_{l-1})\left(H_{k}+\frac{\Delta
t}{2}[\imath H(t_{l}),H_{k}]+\right.\right.\\\ \left.\left.+\frac{(\Delta
t)^{2}}{3!}[\imath H(t_{l}),[\imath
H(t_{l}),H_{k}]]+\cdots\right)U(t_{l-1},0)\right].$

Then, system (1) can be written in its final form as

$\frac{\mathrm{d}\epsilon_{k}^{l}}{\mathrm{d}s}=\frac{1}{2N}\operatorname{Im}\operatorname{Tr}\left[U_{D}^{*}U(T,t_{l-1})\left(H_{k}+\frac{\Delta
t}{2}[\imath H(t_{l}),H_{k}]+\right.\right.\\\ \left.\left.+\frac{(\Delta
t)^{2}}{3!}[\imath H(t_{l}),[\imath
H(t_{l}),H_{k}]]+\cdots\right)U(t_{l-1},0)\right],\quad k=\overline{1,n},\quad
l=\overline{1,L}.$

Numerically solving this system on interval $[0,S]$ with the MATLAB’s ode45
should give controls that implement the desired evolution operator more
accurately than the original method. This improvement could be accounted for
by the use of information about the Hamiltonian’s commutators.

It’s expected that the more terms of the series on the right-hand side of the
expression above we use during computations, the more accurate we get
approximately optimal control. Of course, when using terms of order $(\Delta
t)^{j}$, we have to compute repeated commutators of order $j$ and this can be
very computationally demanding in general. For this reason, a maximum order of
terms to use on the right-hand side should be chosen such that the speed-up of
the new algorithm beats a slow-down induced by the commutators’ computation.

##  3 Numerical experiment

To confirm the expected speed-up provided by the new formula, a quantum system
composed of two spin-$\nicefrac{{1}}{{2}}$ particles and described by the
following dimensionless Hamiltonian ($\hbar=1$) was taken.

$H=\sum_{i=1}^{2}S_{z}^{i}\omega_{i}+\sum_{k=1}^{2}\epsilon_{k}S_{x}^{k}+C_{x}^{(12)}S_{x}^{1}S_{x}^{2}+C_{y}^{(12)}S_{y}^{1}S_{y}^{2}+C_{z}^{(12)}S_{z}^{1}S_{z}^{2},$

where $\omega_{1}=20$, $\omega_{2}=30$, $C_{x}^{(12)}=110$,
$C_{y}^{(12)}=120$, $C_{z}^{(12)}=130$, $S_{i}^{1}=S_{i}\otimes I$,
$S_{i}^{2}=I\otimes S_{i}$, $i=x,y,z$,

$S_{x}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1\\\
1&0\end{pmatrix},S_{y}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&-\imath\\\
\imath&0\end{pmatrix},S_{z}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},I=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}.$

The task was to implement the following operators as accurate as possible

$U_{D}^{1}=\mathrm{e}^{\frac{\imath\pi}{4}}\begin{pmatrix}1&0&0&0\\\
0&1&0&0\\\ 0&0&0&1\\\ 0&0&1&0\end{pmatrix},\quad
U_{D}^{2}=\mathrm{e}^{\frac{\imath\pi}{4}}\begin{pmatrix}1&0&0&0\\\ 0&0&1&0\\\
0&1&0&0\\\ 0&0&0&1\end{pmatrix},$

which are respectively the CNOT quantum gate (controlled NOT) and the SWAP
quantum gate (a swap of qubits’ values).

Two models were compared to each other:

$\frac{\mathrm{d}\epsilon_{k}^{l}}{\mathrm{d}s}=\frac{1}{2N}\operatorname*{Im}\operatorname*{Tr}\left(U_{D}^{*}U(T,t_{l-1})H_{k}U(t_{l-1},0)\right),$
(2)
$\frac{\mathrm{d}\epsilon_{k}^{l}}{\mathrm{d}s}=\frac{1}{2N}\operatorname*{Im}\operatorname*{Tr}\left(U_{D}^{*}U(T,t_{l-1})\left\\{H_{k}+\frac{\imath\Delta
t}{2}\left[H(t_{l}),H_{k}\right]\right\\}U(t_{l-1},0)\right).$ (3)

To numerically solve these systems, the MATLAB’s ode45 method was employed
with the absolute error tolerance set to $10^{-4}$. For the CNOT gate, the
initial controls were set to zero, whereas for the SWAP gate the initial
controls were taken by sampling the function
$10^{-5}\sin(\nicefrac{{t}}{{T}})$ owing to the fact that zero initial
controls didn’t allow the optimization to converge.

The special structure of the given Hamiltonian allows a simplification:

$[H(t),H_{k}]=[H_{0}+\sum_{i=1}^{2}\epsilon_{i}(t)H_{i},H_{k}]=[H_{0},H_{k}].$

The final time $T$ was set to several values: $0.5$, $1$, $5$, $10$. The
number of subintervals $L$ was taken as $150$, $300$. The length of the
interval of integration $S$ was taken as the multiples of $100$: $100$, $200$,
$300$, $400$, … The times required to numerically solve the aforementioned
models on $[0,S]$ with the error in the gates’ implementation not greater than
$10^{-7}$ were compared to each other. Computations were performed on a laptop
with a quad-core Intel Core i7-4702MQ 2.20 GHz processor and 8 Gb RAM. Results
of the comparison for the CNOT gate are in tables 1 and 2, for the SWAP gate —
in tables 3 and 4. The time spent to solve the systems was measured with the
MATLAB’s functions tic() and toc().

Table 1 — Results for the CNOT gate.

| $T=10$ | $T=10$ | $T=5$ | $T=5$
---|---|---|---|---
| $L=300$ | $L=150$ | $L=300$ | $L=150$
Method | (2) | (3) | (2) | (3) | (2) | (3) | (2) | (3)
S | 400 | 100 | x* | 600 | 900 | 200 | 1200 | 400
Time, sec. | 73.76 | 35.76 | x | 46.98 | 97.53 | 39.17 | 29.92 | 20.27
Error ($\times 10^{-8}$) | $0.409$ | $3.94$ | x | $0.412$ | $1.007$ | $1.85$ | $3.54$ | $3.19$

*Method didn’t converge.

Table 2 — Results for the CNOT gate (contrinued)

| $T=1$ | $T=1$ | $T=0.5$ | $T=0.5$
---|---|---|---|---
| $L=300$ | $L=150$ | $L=300$ | $L=150$
Method | (2) | (3) | (2) | (3) | (2) | (3) | (2) | (3)
S | 1000 | 800 | 1000 | 700 | 3900 | 3600 | 4000 | 3600
Time, sec. | 36.68 | 37.66 | 11.66 | 10.39 | 66.22 | 63.21 | 19.52 | 18.6
Error ($\times 10^{-8}$) | $2.2$ | $1.4$ | $3.36$ | $5.1$ | $9.36$ | $9.41$ | $7.80$ | $6.93$

Table 3 — Results for the SWAP gate.

| $T=10$ | $T=10$ | $T=5$ | $T=5$
---|---|---|---|---
| $L=300$ | $L=150$ | $L=300$ | $L=150$
Method | (2) | (3) | (2) | (3) | (2) | (3) | (2) | (3)
S | 900 | 300 | x* | 300 | 3200 | 400 | 1900 | 800
Time, sec. | 143.96 | 85.99 | x | 21.52 | 245.22 | 46.68 | 35.13 | 19.31
Error ($\times 10^{-8}$) | $1.08$ | $0.0001$ | x | $0.0026$ | $4.79$ | $8.14$ | $8.43$ | $0.237$

*Method didn’t converge.

Table 4 — Results for the SWAP gate (continued)

| $T=1$ | $T=1$ | $T=0.5$ | $T=0.5$
---|---|---|---|---
| $L=300$ | $L=150$ | $L=300$ | $L=150$
Method | (2) | (3) | (2) | (3) | (2) | (3) | (2) | (3)
S | 2700 | 2600 | 2900 | 2500 | 3100 | 3200 | 3400 | 3200
Time, sec. | 67.25 | 65.67 | 22.22 | 22.77 | 52.05 | 48.67 | 17.69 | 16.01
Error ($\times 10^{-8}$) | $9.005$ | $7.57$ | $9.57$ | $6.9$ | $9.42$ | $4.24$ | $9.42$ | $4.94$

For $T=0.01$ and $T=0.001$, the both methods didn’t give accurate results,
which perfectly agrees with the fact that in similar tasks there is a lower
bound on time to construct control fields [5].

These results show that the improved method proposed in this paper almost
always achieves the desired accuracy in less time and on a smaller interval
$[0,S]$ than the original method. This improvement is noticeable only for
large $T$ and large $\Delta t$, which agrees with the corrections being
proportional to $\Delta t$ and hence they being small given small $\Delta t$.
Moreover, for some large $\Delta t$ the original method doesn’t converge at
all whereas the improved method achieves the desired accuracy.

##  4 Conclusion

New formulas for constructing optimal controls in a quantum system were
obtained in the form of corrections to the original D-MORPH method [1, 2, 3,
4, 5]. By numerically solving a problem in the area of quantum computations,
it was shown that the new formulas do speed up computations compared to the
original method, i. e. they allow implementation of the desired evolution
operator with high accuracy with fewer ode45’s steps and in less time, even if
only one correcting term proportional to $\Delta t$ is used. Thus, inclusion
of additional information about the system’s Hamiltonian (commutators) had a
positive impact on accuracy achieved and on time spent, which is an advantage
for finding optimal controls on computers with moderate performance.

## References

  * [1] Moore, K. W. Search complexity and resource scaling for the quantum optimal control of unitary transformations / R. Chakrabarti, G. Riviello, H. Rabitz // Phys. Rev. A. 2011. Vol. 83(1).
  * [2] Moore, K. W. Exploring constrained quantum control landscapes / H. Rabitz // The Journal of Chemical Physics. 2012. Vol. 137(13).
  * [3] Moore Tibbetts, K. Exploring the trade-off between fidelity and time optimal control of quantum unitary transformations / C. Brif, M. D. Grace, A. Donovan, et. al. // Phys. Rev. A. 2012. Vol. 86(6).
  * [4] Riviello, G. Searching for quantum optimal controls in the presence of singular critical points / C. Brif, R. Long, R. Wu, K. Moore Tibbetts, T. Ho, H. Rabitz // Phys. Rev. A. 2014. Vol. 90(1).
  * [5] Riviello, G. Searching for quantum optimal controls under sever constraints / K. Moore Tibbetts, C. Brif, R. Long, et. al. // Phys. Rev. A. 2015. Vol. 91(4).