Various methods of optimizing control pulses for quantum systems with decoherence

We study three methods of obtaining an approximation of unitary evolution of a quantum system under decoherence. We use three methods of optimizing the control pulses: genetic optimization, approximate evolution method and approximate gradient method. To model the noise in the system we use the Lindblad equaiton. We obtain results showing that genetic optimization may give a better approximation of a unitary evolution in the case of high noise.


I. INTRODUCTION
One of the fundamental issues of quantum information science is the ability to manipulate the dynamics of a given complex quantum system. Since the beginning of quantum mechanics, controlling a quantum system has been an implicit goal of quantum physics, chemistry and implementations of quantum information processing [1].
If a given quantum system is controllable, i.e. it is possible to drive it into a previously fixed state, it is desirable to develop a control strategy to accomplish the required control task. In the case of finite dimensional quantum systems the criteria for controllability can be expressed in terms of Lie-algebraic concepts [2][3][4]. These concepts provide a mathematical tool, in the case of closed quantum systems, i.e. systems without external influences.
It is an important question whether the system is controllable when the control is performed only on a subsystem. This kind of approach is called a local-controllability and can be considered only in the case when the subsystems of a given system interact. As examples may serve coupled spin chains or spin networks [4][5][6][7]. Local-control has a practical importance in proposed quantum computer architectures, as its implementation is simpler and the effect of decoherence is reduced by decreased number of control actuators [8,9].
A widely used method for manipulating a quantum system is a coherent control strategy, where the manipulation of the quantum states is achieved by applying semi-classical potentials in a fashion that preserves quantum coherence. In the case when a system is controllable it is a point of interest what actions must be performed to control a system most efficiently, bearing in mind limitations imposed by practical restrictions [10][11][12][13][14]. The always present noise in the quantum system may be considered as a such constrain [15][16][17][18][19][20][21][22]. Therefore it is necessary to study methods of obtaining piecewise constant control pulses which implement the desired quantum operation on a noisy system.
In this paper we present a general method of obtaining a piecewise-constant controls, which is robust with respect to noise in the quantum system. This means, we wish to perform a unitary evolution on a greater system than the target one, and discard the ancilla. This paper is organized as follows. In Section II we introduce the model of the studied quantum system. Section III shows different approaches to solving the Lindblad equation. Next, in Section III C we introduce genetic programming. Detailed description of the optimization procedure and studied cases can be found in Section III C 3. In Section IV we show results of the numerical simulations. Finally, we summarize this work in Section V.

II. MODEL OF THE QUANTUM SYSTEM
Our goal is to implement the unitary operations U NOT = σ x ⊗ 1l and U SWAP = 1l ⊗ ij |i j| ⊗ |j i| on a quantum system modeled as a an isotropic Heisenberg spin-1/2 chain of a finite length N . We will study two and three qubit systems. The total Hamiltonian of the aforementioned quantum control system is given by where is a drift part given by the Heisenberg Hamiltonian. The control is performed only on the n th spin and is Zeeman-like, i.e.
In the above S i k denotes k th Pauli matrix acting on the spin i. Time dependent control parameters h x (t) and h y (t) are chosen to be piecewise constant. For notational convenience, we set = 1 and after this rescaling frequencies and control-field amplitudes can be expressed in units of the coupling strength J, and on the other hand all times can be expressed in units of 1/J [23].
We model the noisy quantum system using the Markovian approximation with the master equation in the Kossakowski-Lindlbad form where L j are the Lindblad operators, representing the environment influence on the system [24] and ρ is the state of the system. The main goal of this paper is to compare various methods for optimizing control pulses h x (t), h y (t) for the model introduced above. In the next section we present three methods for this purpose. The comparison of the control pulses obtained by different methods is done by applying these pulses into the above model and analysis of the obtained results.

III. VARIOUS APPROACHES TO FIDELITY MAXIMIZATION
In this Section we describe three methods we used to Obtain optimal control pulses. The first method is an approximate method for obtaining a mapping which is close to a unitary one. The second uses an approximate derivative of the mapping with respect to control pulses. Both of these methods allows to compute approximate gradient of the fitness function. In our numerical research we perform optimization with use of the L-BFGS-B optimization algorithm [25]. This algorithm requires efficient gradient computation. Its main advantages lies in harnessing the approximation of Hessian of fitness function (we refer to the section 3 of [26] for details). Finally, we use genetic programming to optimize control pulses without the need for computing gradient of the fidelity function.

A. Approximate method
Assuming piecewise constant control pulses the Hamiltonian in Eq. (4) becomes independent of time during the duration of the pulse. This allows us to simplify the master equation.
For notational convenience, let us write the decoherence part of the Eq. (4) in the following form [27] − The second term in Eq. 4 can be written as This observations allow us to write a mapping representing the evolution of a system in an initial state ρ under Eq. (4) for time t as where −F = −G − iH. The final state of the evolution is where res(·) is a linear mapping defined as res(|φ ψ|) = |φ |ψ .
We can approximate the superoperator A as Note that, only the C(t) term depends on the control pulses. Assuming piecewise constant control pulses, we can write the resulting superoperator as where n is the total number of control pulses and ∆t i is the length of the time interval in which control pulse h i is applied to the system. The derivative of the superoperator with respect to a control pulse h l is We use the fidelity as the figure of merit where N is the number of qubits in the system and A T is the target superoperator. The derivative of the fidelity is given by B. Approximate gradient method In this section we follow the results by Machnes et. al. [28]. In order to introduce the approximate gradient method, we introduce the following notationĤ This allows us to write Eq. (4) in the form The evolution of a quantum map under these equation is given by In order to perform numerical simulations, Eq. (17) needs to be discretized. Given a total evolution time T , we divide it into M small intervals, each of length ∆t = T /M . Hence, the quantum map in the k th time interval is given by We utilize the trace fidelity as the figure of merit for this optimization problem where . This allows us to write the derivative of the fidelity with respect to the control pulses as Since L and iĤ need not commute, we can not calculate the derivative ∂X k ∂hj (t k ) using exact methods. The best approach is to use the following approximation for the gradient This approximation is valid provided that ∆t 1 ||iĤ + L|| 2 .

C. Genetic programming
Genetic programming (GP) is a numerical method based on the evolutionary mechanisms [29,30]. There are two main reasons for using GP for finding optimal control pulses. First of all it enables to perform numerical search in complicated, mathematically untraceable space. On the other hand, one should note that the values of control pulses in different time intervals can be set independently. Thus the idea of genetic code fits well as a model for a control setting. Thanks to such an representation, genetic programming enables to exchange values of control pulses in some fixed intervals between control settings that results with maximally accurate approximation of the desired evolution.

General GP algorithm
Genetic programming belongs to the family of search heuristics inspired by the mechanism of natural evolution. Each element of a search space being candidate for a solution is identified with a representative of a population. Every member of a population has its unique genetic code, which is its representation in optimization algorithm. In most of the cases genetic code is a sequence of values from a fixed set Σ of possible values of all the features that characterize a potential solution x ∈ Σ n in the search space. Searching for the optimal solution is done by the systematic modification and evaluation of genetic codes of population members due to the rules of the evolution such as mutation, selection, crossover and inheritance.
Mutators are functions that change single elements of a genetic code randomly. A basic example of a mutator is a function that randomly changes values of a representative x at all positions with some non-zero probability Crossovers implement the mechanism of inheritance. This function divides given parental genetic code and create a new genetic code. Commonly two new codes are created at the same time from two parental codes. An example of such crossover is so called two point cut, where both parental codes (x i , y i ) are cut into three regions and the middle segments are interchanged where c 1 < c 2 are randomly chosen indices. In every iteration of the algorithm all members of the population are evaluated using fitness function f : Σ n → which enables elements ordering. Then, using a selector function, the set of the best members is obtained and used to create a new generation of the population using mutation and crossover functions. There is a number of strategies for defining selector function -from completely random choices to the deterministic choice of best representatives. Strategy based on evolution mechanism makes genetic programming especially usable when parts of genetic code represent features of elements of a search space that can be interchanged between elements independently. In such case GA is expected to find the features that occur in well fitted representatives and mix them in order to find the best possible combination. Pseudo code representing this approach is presented in Listing 1. CrossOver work as defined in Section III C.
While the customization of population representation and fitness function unavoidably relies on the optimization problem, other parameters of genetic programming such as crossover and mutation methods are universal.

Customization of the GP
Using GP schema requires obtaining proper representation of a problem. First of all one need to model the space of possible solutions as a set of genomes, usually by representing each of unique and independent features of a solution as one gene. Secondly, it is necessary to define a fitness function that allows to estimate genomes in a way that is consistent with the optimization problem. Using this function one need to determine selection method. The last step is to define methods for modifying modeled genomes. It may be methods based on mutation, crossing-over or both.
In this work we investigate methods for optimizing the sequence of control pulses in order to perform given unitary evolution. As we assume that control pulses in each time step are independent, it is a natural to define each of subsequent control pulses as gene. In this case genome is modeled as sequence od real numbers. It is a very convenient method, because there are many already developed mutation and crossing-over methods for such genome model that have been successfully applied to the problems of searching for the optimal evolution [31,32]. Methods of selection do not depend on genome model and can be based on a variety of already existing ones as well.
Since GP schema does not require the ability to compute gradient of considered fitness function one can use any method for genome estimation. In this work we use control pulses represented by given genome to perform simulation of the system evolution and obtain resulting state that is compared with the one resulting from target evolution. For comparison purposes we apply functions described in section III C 3.

Optimization
In order to optimize a controlled evolution of a system governed by the Lindblad equation, we perform optimization of the average of distances between target state operator and the resulting states for each basis matrix of the space of input states. Our fitness function is defined as where Tr A denotes tracing out the ancila, N is the total number of qubits in the system, ρ i 0 denotes the i th basis matrix, ρ i T = U T Tr A (ρ i 0 )U † T is the target density matrix and Φ is the quantum channel corresponding to the time evolution of ρ i 0 under Eq. (4) for control pulses [c 0 , ..., c M ]. We study two noise models, the amplitude damping and phase damping noise. The former is given by the Lindblad operator L 1 = σ − = |0 1|, while the latter is L 2 = σ z .

IV. RESULTS AND DISCUSSION
In this section we present the results obtained for all of the methods introduced in section III. The comparison of the control pulses obtained by different methods is done by applying these pulses into the model from section II and analysis of the results.
In all of the simulations, we set the number of control pulses to 32 for two-qubit systems and 128 for the three-qubit systems. We limit the strength of the control pulses to h max = 100. In each case we split all of the qubits forming the system into two subsystems: the one that performs some fixed evolution and the auxiliary one.
To find the best spin chain configuration, we study the following systems: 1. One-qubit system with one-qubit ancilla. The control is performed on the target qubit.
3. One-qubit system with one-qubit ancilla. The control is performed on the ancillary qubit.
5. Two-qubit system with one-qubit ancilla. The control is performed on the ancillary qubit.
In our study we choose the phase damping and the amplitude damping channels as noise models. The former is given by the Lindblad operator L = σ z , the latter is given by the operator L = σ − = |0 1|. In order to objectively compare the algorithms, we study two setups: with equal number of steps of the algorithm and with equal computation time. The number of steps and length of these intervals are shown in Table I and Table II Figure 2 shows the results for the phase damping channel. In this case, the methods perform very similarly for the majority of studied cases. The approximate evolution method tends to perform poorly for high noise values. This is due to the fact that in this method, we have periods of coherent evolution followed by decoherence, as Equation (11) states. This is in contrast with genetic optimization, where we make no approximations and use the Lindblad equation. This results in simultaneous decoherence and control. Thanks to this fact, the genetic optimization gives the best Next, in Figure 3 we show the results for the amplitude damping channel. Similarly to the phase damping case, the approximate evolution method performs the worst. The approximate gradient method gives far better results, especially for high values of γ. Again, this may be explained by the fact that in the approximate evolution, we have periods of coherent evolution, followed by decoherence. On the other hand, the genetic optimization performs quite well, even for high noise values we were able to find sets of control parameters which gave fidelity higher than a half. Again, the trade-off was the computation time. Genetic optimization took about an order of magnitude longer compared to other methods.
Unfortunately, all the studied methods appear to fail for γ > 0.01 in the majority of investigated cases. For γ = 0.01 only the genetic optimization gives a high value of fidelity.
The next step of our study focuses on comparing the algorithms when the computation times are on the same order of magnitude. To achieve this, we added more control pulses in the gradient based methods. In the case of both gradient based methods we used 128 pulses for the two-qubit scenario and 512 for the three-qubit scenario. The computation times are summarized in Table III. Results are presented in Figures 4 and 5. Also in this setup the genetic optimization performs better compared to gradient based approaches. In this case the gap between these methods is narrower than in the case with equal number of control pulses.                             We studied different methods of obtaining piecewise constant control pulses that implement an unitary evolution on a system governed by Kossakowski-Lindblad equation. The studied methods included genetic optimization and the BFGS algorithm with the use of fidelity gradient based on an approximate evolution of the quantum system and an approximate gradient method for the exact evolution case. Our results show that, by adding an ancilla, it is possible to implement a unitary evolution on a system under the Markovian approximation. Furthermore, the results heavily depend not only on the size, but also on the location of the ancilla in the spin chain.
What one can notice about the possibility to perform unitary computation in noisy quantum systems is that in majority of the cases it is much better to treat non-target qubits as an ancilla. When comparing systems extended with auxiliary qubits (labeled as a, c, e) and systems with no ancilla (b, d, f ) the difference in possible approximation is emphatic.
The genetic optimization method outperforms gradient based methods in two studied setups: with equal number of control pulses and with equal computation time.