4.3. Optimization on a closed system: creation of the CNOT gate in a two-qubit system

This tutorial demonstrates the design of pulses for the creation of a CNOT gate in a two-qubit system. This is done using both qocttools and, for the sake of comparison, the GRAPE implementation of QuTiP.

The model for the 2-qubit system that we will use is the following: the static Hamiltonian is given by: \begin{equation} H_0 = \frac{1}{2}\omega_1 \sigma_z \otimes I_2 + \frac{1}{2}\omega_2 I_2 \otimes \sigma_z\,, \end{equation} whereas we will assume that one can use the two following perturbations, that permit to control the system and implement the CNOT gate: \begin{align} V_1 &= \sigma_x \otimes \sigma_x\,, \\ V_2 &= I_2 \otimes \sigma_x + \sigma_x \otimes I_2\,. \end{align}

import numpy as np
import matplotlib
from copy import deepcopy
import matplotlib.pyplot as plt
import qutip as qt
# WARNING: The use of the quantum control implementation in QuTiP is commented out
# in this tutorial, as there seems to be a problem in the prerelease version of QuTiP 5.
#import qutip_qtrl.pulseoptim as cpo
import qutip_qip.operations as operations
import datetime
import nlopt

import qocttools
import qocttools.hamiltonians as hamiltonians
import qocttools.target as target
import qocttools.pulses as pulses
import qocttools.qoct as qoct
import qocttools.solvers as solvers

qocttools.about()

No git dir found

Running qocttools version 0.0.14

data = []

We start by defining the static or “drift” Hamiltonian H0. We will use a qubit with frequency \(\omega_1 = 1.0\), and another one with different frequency \(\omega_2 = 2.0\).

w1 = 1.0
w2 = 2.0
h01 = 0.5 * w1 * qt.tensor(qt.sigmaz(), qt.qeye(2))
h02 = 0.5 * w2 * qt.tensor(qt.qeye(2), qt.sigmaz())
H0 = h01 + h02
# To use in qocttools, we need to eliminate the internal tensor structure of the QuTiP object,
# and use simple 1-1 tensors as operators.
H0 = qt.Qobj(H0, dims = [[4], [4]])

Now we define the two perturbations that we will use.

V1 = qt.tensor(qt.sigmax(), qt.sigmax())
V2 = qt.tensor(qt.sigmax(), qt.qeye(2)) + qt.tensor(qt.qeye(2), qt.sigmax())
V1 = qt.Qobj(V1, dims = [[4], [4]])
V2 = qt.Qobj(V2, dims = [[4], [4]])
V = [V1, V2]

The CNOT gate is defined in QuTiP, although we have to reshape the object to be used in qocttools

UCNOT = qt.Qobj(operations.cnot(), dims = [[4], [4]])

For the total propagation time \(T\), we will use five times the period associated to the frequency of the first qubit. We will also define \(\omega_0 = \frac{2\pi}{T}\) as the base frequency for the Fourier expansion.

omega = w1
tau = 2.0*np.pi/omega
ncycles = 5
T = ncycles * tau
# omega0 will be the base frequency for the Fourier expansion
omega0 = 2.0*np.pi/T
print("omega0 = {:.4f}".format(omega0))

omega0 = 0.2000

4.3.1. Optimization with qocttools

Let us perform the calculations with qocttools. We start by building the hamiltonian object:

H = hamiltonians.hamiltonian(H0, V)

Next, we create the Target. The type of target in this case is a “gate” or evolution operator:

tg = target.Target("evolutionoperator", Utarget = UCNOT)

Now, we must define the time-dependent functions that will be used as controls; the total Hamiltonian will be given by: \begin{equation} H(u, t) = H_0 + f_1(u^{(1)}, t)V_1 + f_2(u^{(2)}, t) V_2\,. \end{equation} where \(u\) represents all the control parameters, and is composed of two sets of parameters, \(u=(u^{(1)}, u^{(2)})\), one for each control function. We will use a normal Fourier parametrization for the control functions, and therefore those parameters are the coefficients of the Fourier expansions. See the documentation on the pulse class for details.

First, we will define some initial values for the control parameters. The first step is to set the cutoff frequency; the Fourier expansion will be made up of the frequencies \(\omega_i = i\omega_0\) for \(i=0,1,2,\dots M\). Therefore, one has to choose a value for \(M\).

Then, in order to choose the initial control parameter values, one often wants to use random ones. Here, in order to ensure the reproducibility of the tutorial, we will reuse some values that were previously randomly generated. Note that, in order to generate those random numbers, one needs to set an amplitude bound, i.e. \(\vert u_k\vert \le \kappa\).

M = 10 # The cutoff index of the Fourier expansion
kappa = 0.1 # The amplitude bound for the initial random values.

# The following code sets random numbers for the control parameters
#bound = kappa
#a = -bound
#b = bound
#u1 = (b-a) * np.random.rand(2*M+1) + a
#u2 = (b-a) * np.random.rand(2*M+1) + a
#print(repr(u1))
#print(repr(u2))

# But instead, we will hardcode here values that were generated before.
u1 = np.array([-0.03203772,  0.08398405, -0.07384907, -0.06540511,  0.04327422,
       -0.03974848, -0.06372874,  0.01909529,  0.09770391, -0.00169402,
       -0.0471008 , -0.05903078, -0.09347964, -0.02234591, -0.0223768 ,
        0.09680869,  0.04258134, -0.05684655, -0.01742488, -0.00879573,
       -0.01601734])
u2 = np.array([ 0.05399459,  0.08808636, -0.00589051, -0.00989415,  0.03005704,
        0.07721242, -0.0536454 ,  0.03674837, -0.03142538,  0.07709192,
       -0.03038883, -0.06144648,  0.0870775 , -0.01474247,  0.06243688,
        0.01285235, -0.08496945, -0.08987118,  0.02714249, -0.01303071,
        0.00023745])

We can now create the two pulse objects, that we place into a list:

f = [pulses.pulse("fourier", T, u1), pulses.pulse("fourier", T, u2)]
f0 = deepcopy(f) # We create a copy to store the initial values
# u0 will store the whole collection of initial parameters of the two pulses
u0 = pulses.pulse_collection_get_parameters(f)

Now we create the Qobj object. We need to specify the number of steps in the time discretization that is used to do the propagations.

ntsteps = ncycles * 20
U0 = qt.identity(4)
opt = qoct.Qoct(H, T, ntsteps, tg, f, U0,
                interaction_picture = False)

Before doing the optimization itself, let us do an initial propagation to see how the initial guess pulses perform. We compute the value of the target function along the propagation. In this case where the target is the creation of a gate, the definition of the target function is: \begin{equation} F(U) = \frac{1}{{\rm dim}^2} \vert {\rm Tr} U^\dagger U_{\rm target} \vert^2 \end{equation} In our case, the target gate is the CNOT gate, and \({\rm dim}=4\).

def F(U):
    return (1/4**2) * np.abs((UCNOT.dag() * res[j]).tr())**2

ts = np.linspace(0, T, ntsteps)
res = solvers.solve('cfmagnus4', H, f, qt.qeye(4), ts)
Ftinitial = np.zeros(ntsteps)
for j in range(ntsteps):
    Ftinitial[j] = F(res[j])

Now we can launch the optimization. In this case, we will use the LBFGS algorithm from the nlopt library (in principle, one can use any of the algorithms defined there). We will set a maximum of 100 function evaluations as a limit to stop the run, and we will also ask the algorithm to stop whenever a stopval value is reached, that we set to 0.99 (ideally, a full overlap with the target gate is achieved when the target function is one). See the documentation for the maximize method for more details and options.

# Before the optimization, one could do a check on the quality of the gradient with this call:
#derqoct, dernum, error = opt.check_grad(u0)
#print(derqoct, dernum, error)

x, optval, result = opt.maximize(algorithm = nlopt.LD_LBFGS,
                                 maxeval = 100,
                                 stopval = 0.99,
                                 verbose = True)
data.append(optval)

0 0.191753 (0.045546 s)
1 0.072449 (0.044287 s)
2 0.093749 (0.044134 s)
3 0.249927 (0.034458 s)
4 0.287447 (0.040124 s)
5 0.242607 (0.044513 s)
6 0.380359 (0.034123 s)
7 0.398719 (0.040273 s)
8 0.550735 (0.044356 s)
9 0.557701 (0.043438 s)
10 0.725678 (0.034986 s)
11 0.743634 (0.042686 s)
12 0.772435 (0.034957 s)
13 0.820350 (0.044193 s)
14 0.855783 (0.044028 s)
15 0.891240 (0.035857 s)
16 0.909609 (0.033904 s)
17 0.924110 (0.044283 s)
18 0.946126 (0.044127 s)
19 0.960249 (0.034228 s)
20 0.970578 (0.040936 s)
21 0.976559 (0.035011 s)
22 0.984004 (0.044072 s)
23 0.993447 (0.034598 s)
Successfull termination with result code = 2
("Stop value reached")

Let us check how the optimal pulse does. The optimal parameters are stored in the x array. The following updates the list of pulses f with those values:

pulses.pulse_collection_set_parameters(f, x)

Now we can re-propagate with the optimized pulses:

res = solvers.solve('cfmagnus4', H, f, qt.qeye(4), ts)
Ft = np.zeros(ntsteps)
for j in range(ntsteps):
    Ft[j] = F(res[j])

And finally, we plot the results. In the top panel, the value of the target function as it evolves in time, for the initial guess and for the optimal pulses. In the bottom panel, those pulses.

plt.rcParams["font.size"] = "8"
plt.rcParams["figure.autolayout"] = False

fig, ax = plt.subplots(2, 1, figsize=(5, 5), sharex = True)
fig.subplots_adjust(hspace=0)

ax[0].plot(ts/tau, Ftinitial, label = r"$F(U_{\rm initial}(t))$")
ax[0].plot(ts/tau, Ft, label = r"$F(U(t))$")
#ax[0].set_xlabel(r"$t/\tau$")
ax[0].set_xlim(0, T/tau)
ax[0].legend(bbox_to_anchor = (1.0, 1.0), loc = 'upper left')

ts = np.linspace(0, T, ntsteps)

ax[1].plot(ts/tau, f0[0].fu(ts), label = r"$g_x^{\rm ini}(t)$")
ax[1].plot(ts/tau, f0[1].fu(ts), label = r"$g_y^{\rm ini}(t)$")
ax[1].plot(ts/tau, f[0].fu(ts), label = r"$g_x^{\rm opt}(t)$")
ax[1].plot(ts/tau, f[1].fu(ts), label = r"$g_y^{\rm opt}(t)$")
ax[1].set_xlabel(r"$t/\tau$")
ax[1].set_xlim(0, T/tau)
#ax[1].label_outer()
ax[1].legend(bbox_to_anchor = (1.0, 1.0), loc = 'upper left')

#fig.savefig("closed-gate.pdf", bbox_inches = 'tight')
plt.show()

4.3.2. Calculations with the GRAPE algorithm as implemented in QuTiP

4.3.3. (WARNING: TEMPORARILY DISABLED)

We will now approach the same problem, but using the GRAPE algorithm as implemented in QuTiP. This means that we can no longer use the Fourier parametrization of the pulses; instead, the basic assumption is that the pulses are piece-wise constant functions. The control parameters are the values of the function on each of those steps. This is also often called a real-time parametrization of the function.

See the QuTiP documentation and tutorials for an explanation of the following steps.

#logger = qt.logging_utils.get_logger()
#log_level = qt.logging_utils.INFO

## Number of time slots; we will use the same number of steps than we did with qocttools.
#n_ts = ntsteps -1 # 10
## Time allowed for the evolution
#evo_time = T # 10

## Fidelity error target
#fid_err_targ = 1e-10
## Maximum iterations for the optisation algorithm
#max_iter = 200
## Maximum (elapsed) time allowed in seconds
#max_wall_time = 120
## Minimum gradient (sum of gradients squared)
## as this tends to 0 -> local minima has been found
#min_grad = 1e-20

## pulse type alternatives: RND|ZERO|LIN|SINE|SQUARE|SAW|TRIANGLE|
#p_type = 'SAW'

#result = cpo.optimize_pulse_unitary(H0, V, U0, UCNOT, n_ts, evo_time,
#                fid_err_targ=fid_err_targ, min_grad=min_grad,
#                max_iter=max_iter, max_wall_time=max_wall_time,
#                out_file_ext=None, init_pulse_type=p_type,
#                log_level=log_level, gen_stats=True)

#result.stats.report()
#print("Final evolution\n{}\n".format(result.evo_full_final))
#print("********* Summary *****************")
#print("Final fidelity error {}".format(result.fid_err))
#print("Final gradient normal {}".format(result.grad_norm_final))
#print("Terminated due to {}".format(result.termination_reason))
#print("Number of iterations {}".format(result.num_iter))
#print("Completed in {} HH:MM:SS.US".format(
#        datetime.timedelta(seconds=result.wall_time)))

The calculations are very fast, and the final achieved fidelity is very high. Let us now propagate the evolution operator, both for the initial and for the optimized pulses, to double-check that everything worked, and also to be able to plot later the evolution of the target function in time.

#dt = result.time[1]
#UTinitial = qt.identity(4)
#UToptimized = qt.identity(4)
#Ftinitial = np.zeros(ntsteps)
#Ft = np.zeros(ntsteps)
#Ftinitial[0] = (1/4**2) * np.abs((UCNOT.dag() * UTinitial).tr())**2
#Ft[0] = (1/4**2) * np.abs((UCNOT.dag() * UToptimized).tr())**2
#for j in range(n_ts):
#    Hoptimized = H0 + result.final_amps[j, 0] * V1 + result.final_amps[j, 1] * V2
#    Uoptimizeddt = (-1j * dt * Hoptimized).expm()
#    UToptimized = Uoptimizeddt*UToptimized
#    Ft[j+1] = (1/4**2) * np.abs((UCNOT.dag() * UToptimized).tr())**2
#    Hinitial = H0 + result.initial_amps[j, 0] * V1 + result.initial_amps[j, 1] * V2
#    Uinitialdt = (-1j * dt * Hinitial).expm()
#    UTinitial = Uinitialdt*UTinitial
#    Ftinitial[j+1] = (1/4**2) * np.abs((UCNOT.dag() * UTinitial).tr())**2

#plt.rcParams["font.size"] = "8"
#plt.rcParams["figure.autolayout"] = False

#fig, ax = plt.subplots(2, 1, figsize=(5, 5), sharex = True)
#fig.subplots_adjust(hspace=0)

#ax[0].plot(ts/tau, Ftinitial, label = r"$F(U_{\rm initial}(t))$")
#ax[0].plot(ts/tau, Ft, label = r"$F(U(t))$")
##ax[0].set_xlabel(r"$t/\tau$")
#ax[0].set_xlim(0.0, T/tau)
#ax[0].legend(bbox_to_anchor = (1.0, 1.0), loc = 'upper left')

#ts = np.linspace(0, T, ntsteps)

#ax[1].step(result.time/tau,
#           np.hstack((result.initial_amps[:, 0], result.initial_amps[-1, 0])),
#           where='post', label = r"$g_x^{\rm ini}(t)$")
#ax[1].step(result.time/tau,
#           np.hstack((result.initial_amps[:, 1], result.initial_amps[-1, 1])),
#           where='post', label = r"$g_y^{\rm ini}(t)$")
#ax[1].step(result.time/tau,
#           np.hstack((result.final_amps[:, 0], result.final_amps[-1, 0])),
#           where='post', label = r"$g_x^{\rm opt}(t)$")
#ax[1].step(result.time/tau,
#           np.hstack((result.final_amps[:, 1], result.final_amps[-1, 1])),
#           where='post', label = r"$g_y^{\rm opt}(t)$")

#ax[1].set_xlabel(r"$t/\tau$")
#ax[1].set_xlim(0.0, T/tau)
##ax[1].label_outer()
#ax[1].legend(bbox_to_anchor = (1.0, 1.0), loc = 'upper left')

#plt.show()

# This file is used by the testing script of the code.
with open("data", "w") as f:
    for i in data:
        f.write("{:.14e}\n".format(i))