A Stable fixed point topology of PMP-iterative propagator control algorithms

Search complexity and resource scaling for the quantum optimal control of unitary transformations


The optimal control of unitary transformations is a fundamental problem in quantum control theory and quantum information processing. The feasibility of performing such optimizations is determined by the computational and control resources required, particularly for systems with large Hilbert spaces. Prior work on unitary transformation control indicates that (i) for controllable systems, local extrema in the search landscape for optimal control of quantum gates have null measure, facilitating the convergence of local search algorithms; but (ii) the required time for convergence to optimal controls can scale exponentially with Hilbert space dimension. Depending on the control system Hamiltonian, the landscape structure and scaling may vary. This work introduces methods for quantifying Hamiltonian-dependent and kinematic effects on control optimization dynamics in order to classify quantum systems according to the search effort and control resources required to implement arbitrary unitary transformations.


1 in1 in1 in1 in

1 Introduction

The methodology of optimal control theory (OCT) has been applied to achieve various dynamical objectives in quantum systems by manipulating constructive quantum wave interference to maximize the likelihood of attaining desired target states [1]. Three classes of problems - state control, observable control, and unitary transformation or gate control [2] - have received the most attention in the quantum control community to date. The generation of targeted unitary transformations is of both fundamental interest and has direct applications to quantum information sciences since the quantum logic gates required to carry out quantum computation are represented by unitary transformations [3]. Since (and ) are compact Lie groups, it is possible to generate any through sequential application of elements of a complete set of generators for , i.e., . This strategy (uniform finite generation) is now commonly applied in gate decomposition strategies wherein the unitary gate expressed in terms of qubits (i.e., with a corresponding dimensional Hilbert space) is constructed through sequential application of various which each act on only 1-2 qubits [4, 3]. However, provided the system is controllable, it is also possible to generate any by shaping time-dependent control functions over a time interval . Each control function is coupled to a corresponding control Hamiltonian , which are all simultaneously applied in order to generate the desired at time . This is the method of optimal control theory. Typically, uniform finite generation requires a greater total evolution time than OCT methods based on pulse shaping [4].

OCT has been applied to unitary transformation control for the purposes of quantum computation in a variety of quantum systems. Kosloff and coworkers studied the implementation of quantum gates based on vibrational eigenstates of the molecular sodium ion on the ground electronic surface [5, 6]. Similar studies were carried out similar studies in the acetylene molecule using the asymmetric C-H stretching and bending modes by de Vivie-Riedle and coworkers [7]. Gate control on spin-system dynamics with optimally designed NMR pulses has been performed by several groups [8, 9, 10, 11]. Deutsch and coworkers implemented unitary maps on the magnetic sublevels of the ground electronic state of cesium [12]. The computational studies using OCT can be extended to the laboratory by shaping ultrafast laser fields using Optimal Control Experiment (OCE) [13] to generate control functions .

An important issue in determining the feasibility of optimally constructing unitary transformations is how the required search effort scales with the Hilbert space dimension of the system under control. Whereas state control and maximization of observable expectation values have met with widespread success in experimental and computational incarnations [2, 14], with search effort generally invariant to the Hilbert space dimension [15, 16], the achievement of high fidelity unitary transformations has proved more challenging [5, 6, 15], especially for large systems.

Recently, a series of fundamental studies have been carried out on the underlying properties of quantum control landscapes, defined as the map between the external control field and the objective fidelity [8, 17, 18, 19, 20]. These studies revealed that under reasonable assumptions and in the absence of auxiliary costs (e.g., on the field fluence) or constraints, the control landscape contains no sub-optimal extrema, or “traps” that can hinder a gradient-based local algorithm for finding an optimal control field. The importance of the landscape topology to determining the feasibility of quantum control is beginning to be more widely recognized [12]. The topology of quantum control landscapes is dominated by so-called kinematic extremals, which are determined by the cost function alone and independent of the Hamiltonian used [17, 19, 20]. Although the trap-free landscape topology ensures convergence of unconstrained gradient-based algorithms given sufficient effort, the topology does not specify the convergence rate of such algorithms, which additionally depends on the local landscape structure (e.g., slope and curvature). In this paper we examine the effects of control landscape features on the convergence rate of first-order algorithms for the optimal generation of unitary transformations. The primary goals are (i) to assign different Hamiltonians to classes exhibiting exponential or sub-exponential scaling with the system size, and (ii) to quantify the effects of the local landscape structure on the convergence rate.

The paper is organized as follows. In Section 2, the theoretical formulation of the control problem is presented, including a summary of the associated landscape topology. Section 3 provides a framework that unifies first-order OCT algorithms for gate control, demonstrating that the convergence of all these algorithms is governed by the same underlying landscape topology. Section 4 presents metrics for landscape slope and curvature, including kinematic bounds, and dynamical metrics that quantify the effect of the control system Hamiltonian. Section 5 defines the control systems and target propagators used in the simulations. Section 6 presents numerical results on control optimization search effort and resource scaling with respect to Hilbert space dimension and identifies classes of Hamiltonians exhibiting exponential and sub-exponential scaling, while Section 7 relates search effort and resource scaling to local landscape structure. Finally, in Section 8, we draw conclusions from the findings.

2 Optimal Control Theory

2.a Dynamical Formulation of the Control Objective

Consider an -dimensional isolated quantum system whose dynamics are governed by the time-dependent Schrödinger equation,


where is the time-dependent Hamiltonian whose control variables are denoted as . In atomic units, the unitary propagator at some final time is


where T is the time-ordering operator and is implicitly understood to be a function of , which in this work is represented by an external control field . We consider an isolated quantum system satisfying the dipole formulation where is the field-free (drift) Hamiltonian and is the dipole or control Hamiltonian operator.

This work concerns the class of control objective functionals


The endpoint control objective may be defined as guiding the system’s unitary propagator to match a pre-specified unitary matrix . A convenient cost function is to minimize the Hilbert-Schmidt distance ,


where the desired minimum2 is achieved when , and the global maximum corresponds to . The quantum system is assumed to be controllable such that any desired can be generated by some choice of at time . (see Section 4.E).

In this work, the objective of Eq. (4) is optimized using dynamical controls present in an external electric field . We consider a controllable quantum system with levels in . To determine an optimal control that maximizes or minimizes Eq. (3), it is useful to define a Lagrangian functional that directly imposes the dynamical constraint in Eq. (1):


where is a Lagrange multiplier matrix function and . Denoting by the Hilbert-Schmidt inner product , the first integrand term

is the PMP (Pontryagin maximum principle) Hamiltonian function [21, 22]. A necessary condition for maximizing or minimizing Eq. (3) subject to the dynamical constraint is satisfaction of the first-order conditions (Euler-Lagrange equations) for the Lagrangian [22]. The first Euler-Lagrange equation is simply the Schrödinger equation (1). The second Euler-Lagrange equation of (5) is


where satisfies the boundary condition . For given by Eq. (4), we have [23]


The third Euler-Lagrange equation (critical condition) is . For a control system satisfying Eq. (1),


Consider the control-propagator map and the composition of maps . Then, the functional derivative evaluated at time is denoted as . For simplicity of exposition, we use the symbols interchangeably with , respectively, and refer simply to the derivative . Control fields that satisfy Eqs. (1), (6), and (8) constitute the critical points of the control landscape . Since the problem with J given by (3) is underdetermined, these critical points lie on critical submanifolds, each consisting of an infinite number of fields which produce the same value of .

2.b Critical Topology of

The matrix which appears in equation (8) is the functional derivative . A simplifying condition that facilitates extraction of the critical topology of is that the Hermitian Gramian matrix


is full rank [24]. Here, denotes the “vectorization” of an complex matrix into an -component complex vector and . Note . Satisfaction of the full-rank condition ensures that (8), which may be written , implies (see Section 4.B). The condition is verified numerically for diverse classes of quantum control systems in Section 7.

When this condition is satisfied, the critical topology (number of local optima and their optimality status) of the control landscape is equivalent to that of in Eq. (4) [19, 25]. The control-propagator map associates with each critical submanifold of a critical submanifold of whose number of positive and negative Hessian eigenvalues is identical, but which has an infinite-dimensional nullspace [20]. The topology may thus be analyzed by considering the degrees of freedom, using for example, the matrix elements of as controls. It has been shown that there are distinct critical values of , =0, 4, 8, … corresponding to critical points where [19, 25]. The topology of these critical points can be determined by considering the Hessian operator of ,


where the are a suitable set of local kinematic coordinates around the critical point . Of interest is the number of the positive, negative, and zero eigenvalues of the Hessian at , which correspond to the number of upward, downward, and flat directions at that point.

The Hessian eigenvalue enumeration may be obtained from the Hessian quadratic form (HQF) of :


where is an infinitesimal Hermitian matrix [25]. At a critical point, it may be shown [25] that


for some unitary and where . Evaluating the HQF explicitly at a critical point and writing the elements of as yields terms ( terms in the first sum and terms in the second sum) corresponding to the eigenvalues of the Hessian [25],


The sign of each term in Eq. (13) corresponds to the sign of each Hessian eigenvalue. It has been shown [25, 26] that for a critical point with an objective value of for , the number of positive (), negative () and zero () type of eigenvalue is


For =0, =0, there are positive eigenvalues, indicating that the optimum is an isolated point. Similarly, for =4, where , there are negative eigenvalues. For intermediate values of , there are a mixture of positive, negative, and zero eigenvalues, indicating that all intermediate critical points have a saddle topology. For example, at =4, =1, and , , and . Assuming that minimization of is desired, this saddle point may be expected to pose a hindrance to search effort, as there is only one direction out of leading down to the global minimum. The higher saddles contain more negative eigenvalues and thus are expected to pose less of a hindrance in search effort. This matter will be examined in Section 6.A.

3 Optimization Methods

For unitary transformation control, deterministic first-order algorithms are typically used for control optimization. In this Section we compare these first-order algorithms and demonstrate that they share a common fixed point topology. In Section 4 we extend these results to demonstrate that the algorithms share common bounds on their convergence rates to identify optimal controls.

The simplest first-order algorithm is the gradient flow of the objective function. Using the variable to index the search path, the gradient flow trajectory is the solution to the initial value problem


for a specified initial guess for the control , where is an adaptive step size. Associated with the control field trajectory is a trajectory for in for the final dynamical propagator, which is induced by the control-propagator map . In the numerical simulations in this work, Eq. (15) will be solved using a variable step size fourth order Runge-Kutta integrator built into MATLAB [27]. A primary concern in this paper is the convergence rate of such algorithms, whose fixed points include all such that in Eq. (12).

As , converges toward stable fixed points that are critical points of . These points are neutrally stable, i.e., within any neighborhood of consisting of controls such that

there exists a subneighborhood such that if , for all .

The neutrally stable solutions are the global optima of , which can be seen from the corresponding trajectory induced by the map . converges to asymptotically stable fixed points that are optima of - points such that

for some that is equal to the radius of the attracting region of the fixed point. The latter are the critical points identified in Eq. (12) with positive definite Hessian (13) - and according to (14), the only critical point satisfying this criterion is the unique global optimum . Due to the asymptotic stability of , any such that is within the attracting region of will converge to a neutrally stable fixed point that lies on the global optimum submanifold of . The instability of fixed points lying on other critical submanifolds with follows from the indefiniteness of the HQF at .

Many unitary control studies use so-called PMP-iterative algorithms [28], which can be formulated only in discrete time. These algorithms iteratively integrate equations (6) and (1) at each step , using control fields , , respectively, where are scalars. The fixed points of these algorithms are points on the control landscape where or . In Appendix A we show that under appropriate regularity conditions the only neutrally stable lie on the global optimum submanifold where . A third class of gate control optimization algorithms consists of first-order tracking algorithms which follow a prescribed path in the space of propagators to the target gate; these have been shown to be capable of achieving gate fidelities approaching machine precision [23]. In this work we focus on the application of steepest descent algorithms (15), due to the mathematical convenience of formulating convergence to a stable fixed point for these continuous time algorithms, but our conclusions on convergence efficiency are applicable to PMP-iterative and tracking algorithms as well.

4 Landscape Structure Metrics

The landscape topology summarized in Section 2.B suggests that an optimal field to achieve a desired unitary transformation may be readily found because no suboptimal extrema exist on the landscape. This attractive behavior does not preclude the possibility that complicated landscape features, including strong influence by saddle regions, may impede optimal searches. Thus, an understanding of the effects of the local landscape structure on optimal searches is necessary in order to explain and predict the scaling of search effort with system size . We introduce landsape metrics and show that the same landscape local structural features govern the convergence of PMP-iterative and gradient-based algorithms.

4.a First-order metrics

The local structure metrics of the landscape are based on a Taylor expansion of the cost functional with respect to . The slope metric at a point on the landscape is defined as


where the unit gradient vector is . The metric is thus equivalent to the magnitude of the gradient on the landscape at the th point. Beginning from the expression in Eq. 16, at any point is bounded by


Above, the Cauchy-Schwarz inequality is used twice. A greater value of at results in a locally faster descent.

We can also establish a bound on in PMP-iterative algorithms (see Appendix A):


In PMP-iterative algorithms, the bound on the increment in the field for infinitesimally small step length is equivalent to that for steepest descent with a Mayer cost.

4.b Second-order metrics

From the second variation of the objective functional , we may derive the Hessian kernel; the elements of the Hessian are given by [23]


At a critical point, the last term of Eq. (18) drops out. The relationship between the HQF expression (11) and (18) is described in Appendix B.

In steepest descent algorithms, according to the gradient expression (8), is composed of linear combinations of the real and imaginary components of . The eigenvalues of the Hessian (18) specify the rates at which new frequency modes required for optimal control (contained within the eigenfunctions of ) can be added to . Thus, several bounds on the Hessian are given below.

The Hessian trace or mean curvature is given by


At a critical point,


At =0, , and Eq. (20) becomes


where the second step uses the cyclic permutation trace rule. Similarly, at the maximum =4, the trace is given by the negative value of Eq. (21).

We can also calculate a bound on Hessian mean curvature away from a critical point:

Finally, we consider the local curvature, or the projection of the Hessian matrix on to the normalized gradient vector ,


The curvature near the optimum may influence the required search effort by determining the ease of convergence to the optimum. Note that since the gradient and Hessian can be expressed in terms of the same basis functions of time, only inner products of components of the time-evolved dipole operator contribute to the local curvature of the control landscape.

4.c Distance metrics

On a search trajectory, the field starts out at algorithmic index =0 with and progresses in steps (i.e., ) until the trajectory ends at an optimal field, at . The complexity of the search may be characterized by the ratio of the trajectory path length to the Euclidian distance between initial and final control fields ,


The closer is to unity, then the more direct the path, i.e., the closer the path is to a straight line in search space. This metric will be used to assess the complexity of the search trajectories followed during optimizations.

Since the presence of saddle manifolds on the landscape may influence the efficiency of an optimal search, the distance of points on the search trajectory to the nearest saddle also provides important structural information. This distance may be measured by examining the eigenvalues of the matrix , since these eigenvalues are all at any saddle (all at =0 and all at =). If are the eigenvalues of , a convenient metric to express the distance to the nearest saddle is


where the normalization factor is if and if , which makes the maximal allowed value of equal to 1 for any value. Finding that is close to zero near the saddle values (=4, 8, 12, …) indicates that the search trajectory encounters the saddle manifold at this value. The effects of search trajectories approaching saddle manifolds will be examined in detail in Section 6.A.

4.d Gramian matrix

Unlike the Hessian, which depends on as well as the system Hamiltonian, the Gramian matrix (9) provides a means of characterizing purely dynamical effects on the optimization trajectory. Consider the equation for the controlled propagator corresponding to the gradient flow (15). Denoting by the vectorization of the propagator , we have [24]:


from which it can be seen that is a linear map between vectors in and all system-dependent effects. The eigenvectors of the Gramian are orthogonal directions in the tangent space to the unitary group . The magnitudes of the eigenvalues of represent the dynamical contributions of the corresponding orthogonal directions in to the propagator variation induced by the gradient flow control variation .

If the eigenvalues of are always sufficiently far from zero, corresponding to a well-conditioned matrix, then can be infinitesimally small only near the global optimum submanifold, irrespective of the direction of , which facilitates convergence. This can be seen as follows: the critical point condition implies

then if at a critical point, then . Since the only neutrally stable fixed points of flow (15) lie on the global optimum manifold, the claim holds.

More generally, control systems for which the expected values of the condition number of are low near the global optima typically exhibit faster convergence, as shown in Section 7. A special case is where is degenerate; then, the dynamical contribution of each eigenvector in contributes equally when the control variation is integrated over time and convergence is governed by the kinematic gradient3.

4.e Higher-order analysis

The above metrics and associated bounds characterize the local properties of the control landscape. Properties of the global solution to the gate control problem can be analyzed using the Dyson series expansion for the controlled unitary propagator in the interaction picture [29]:


where . The -th term in the Dyson expansion corresponds physically to the set of possible -photon transition pathways between eigenstates over time . Note that each successive term in (26) contains higher order products of and than the previous terms. For any bounded field fluence and tolerance , the series converges to at some finite order [29].

The matrices and , which define the control system, also fully determine the minimal order in series (26) required to produce any given . Let denote the dynamical Lie algebra of the quantum control system - i.e., the Lie algebra spanned by repeated commutators of and :

Beyond some critical (called the depth of ), which depends on the control system, the dynamical Lie algebra saturates [4]. For bilinear quantum control systems, a sufficient condition for full controllability - i.e., the existence of a control field such that the corresponding induced by the Schrödinger equation can be any unitary matrix - is that [4, 30].

Higher-order terms in series (26) correspond to higher order commutators, and hence generally require higher field fluences (higher amplitudes of the corresponding Fourier modes). If these amplitudes in are below certain minimal values required for the corresponding Dyson series terms to be nonnegligible, it is not possible for the field to be a solution to the gate control problem. In general, systems with control Hamiltonian (dipole) operators that systematically exclude distant transitions require higher order terms in the Dyson series to reach any arbitrary gate , which increases the nonlinearity of the optimization problem for these control systems with greater Lie algebra depth . This results in the optimal controls containing more frequency modes, corresponding to more complex control mechanisms [29]. In Sections 6 and 7, we demonstrate that optimizing with control Hamiltonians with weak or forbidden transitions between distant quantum states yields higher-fluence and more complex optimal controls , which require a greater algorithmic search effort to find.

Equation (26) also determines when the Gramian matrix (9) may be nonsingular. Only for controls where all the terms in the Dyson series required for full controllability of a given and have become populated can the Gramian matrix be well-conditioned. This is quantified by the relation4

In Section 7, the effects of control system Hamiltonians on Gramian condition number are studied numerically.

5 Control systems

An exhaustive sampling of system structures for the drift Hamiltonian and the control Hamiltonian is impractical. Here, we choose two and four structures motivated by common propagator control systems, which differ qualitatively in their Lie algebra depth, as described in Section 4.E. The goal is to provide an overview of search behavior that might be expected without making any specific predictions for the behavior of any particular quantum system. In order to delimit the space of drift and control Hamiltonians studied, only diagonal structures are considered, with structural variation restricted to the control Hamiltonians . Although control systems requiring multiple fields to ensure controllability (e.g., coupled spin systems) may be used to realize quantum computation [11], we consider only systems controllable by a single field here so that the search behavior across different systems can be directly compared.

We consider an -level quantum system in arbitrary dimensionless units. Two model systems with in its diagonal basis are considered. First is that of a rigid rotor,


and second is that of an anharmonic oscillator,


with and .

Four physically relevant real matrix control Hamiltonian structures will be considered, paired with one of the two structures above. Control Hamiltonians of different matrix distributions of off-diagonal elements were chosen because of their different Lie algebra depths and optimal control mechanisms, as discussed in Section 4.E. All of the control Hamiltonian structures used here have nonzero trace in order to make the systems controllable on and not only on . For many physical systems the coupling between states decreases as the difference between the quantum numbers of the states increases, and the first choice of takes this property into account, with the following structure


where , is the coupling parameter, and all elements of have a random phase of with the restriction that remains symmetric. This “” structure qualitatively corresponds to diatomic molecules and other anharmonic vibrational systems. The second control Hamiltonian structure examined is the related “banded” structure where a fixed number of rows nearest to the diagonal have elements of with the remaining rows having elements of zero; the extreme of having only one row with allowed transitions qualitatively corresponds to a harmonic oscillator or rigid rotor. Third, we consider a “sparse” structure with 50 of the off-diagonal elements randomly chosen as and the remaining 50 of the off-diagonal elements being zero, while maintaining as symmetric. Sparse control Hamiltonians with fewer than 50 allowed couplings are examined as well in Section 6.C. Control Hamiltonian structures containing some allowed and some forbidden transitions qualitatively correspond to coupled-spin system qubit structures commonly used in quantum computation, although only certain specified distributions of couplings are allowed for qubit systems. In order to investigate the search effort for systems with such control operators, we consider a “tensor product” control Hamiltonian on qubits,


where the diagonal matrix is added to make the system controllable on and is

We consider pairings of this with the diagonal operators above 5.

As we will demonstrate in Sections 6 and 7, the distribution of couplings between states is important for assessing the scaling of effort with . With this in mind, comparing the sparse and tensor product structures reveals some important differences, for example at =8.


While each operator has an equal number of allowed transitions, the tensor product structure allows no transitions more than four states apart (note the zeros in the upper-right and lower-left corners of the matrix). The sparse example shown, in contrast, allows transitions between states and , as well as some other transitions five or more states apart. Structural differences in coupling distributions are even more evident at =16 and =32. We thus define two distinct classes of control Hamiltonians: those that allow distant transitions, including the =1.0 and sparse structures, and those that forbid or have very weak distant transitions, including 1.0, banded, and tensor product structures. The differences in coupling distributions between these two classes influence the required search effort, as will be shown in Sections 6 and 7.

The simulations will consider both random Haar-distributed unitary matrices [32] and the Fourier transform quantum gate,


where and denote the matrix elements and run from 1 to .

The initial field at is chosen as


where are the Fourier components of the field, which are selected randomly and bounded by the frequency of the transition in , is a random phase on , and is the field fluence. Prior to multiplication by , the field is normalized to have unit fluence.

6 Control Search Complexity

The search effort required to find an optimal field has important implications for determining the feasibility of controlling the dynamics of complex systems. In Section 6.A, the influence of the saddle point topology of the control landscape is assessed. In Section 6.B, we examine the search effort as a function of for a broad range of choices of , , initial field strength , and . Further exploration of the control Hamiltonian structure’s effect on search effort in Section 6.C identifies control Hamiltonian properties that result in the most efficient searches. Details of the numerical parameters in the simulations are given in Appendix C.

6.a Influence of Landscape Saddle Point Topology

The simulations in this section address how the landscape topology, which is primarily determined by the kinematic cost function , influences the behavior of gradient-based optimizations for systems of dimension up to =8. In particular, we examine the extent to which the saddle manifolds influence the search trajectory and whether the saddle effects are dependent on the choice of or the initial control field. The Hamiltonian is given by Eq. (27) and given by Eq. (29) with =1.0, 0.9, or 0.6.

The trajectories of three searches for =4 are shown in Figure 1. Comparison of the saddle metric (c.f., Eq. 24) in the right panel with the optimization trajectories in the left shows that interaction with saddles retards convergence to the optimum. Examination of the trajectory of the Hessian eigenvalues during the search confirms interaction with a saddle manifold. The Hessian eigenvalues of the search interacting with the =4 saddle are shown versus in Figure 2. At =4 (dotted vertical line), there are nine positive eigenvalues (marked by circles) and one negative eigenvalue (marked by square), in agreement with the Hessian spectrum derived in Eq. (14). Furthermore, at the optimum =0, there are 16 positive Hessian eigenvalues (marked by small circles), in agreement with the maximally allowed positive eigenvalues. The remaining Hessian eigenvalues are null, as predicted [19, 25].

Table I presents statistics on optimizations using a variety of conditions. 1000 searches starting from different initial fields were used to generate the statistics. Shown is the required search effort (defined as the number of algorithmic iterations to reach ) as well as the fraction of searches that interact with saddles. Three degrees of saddle interaction are examined: , , and . The probability of saddle interaction decreases with rising Hilbert space dimension , such that by =8, negligibly few searches have strong interactions with saddles. The decrease in saddle interactions as rises is favorable to performing large-scale unitary transformation optimizations.

6.b Scaling of effort with

Simulations were performed for =2, 4, 8, 16, and 32 with a statistical sample size of 20 (with the exception for some cases of =32, where a single optimization was performed) and a convergence criterion of . The mean search effort with statistical error is shown in Table II for all optimizations. The effort is plotted versus for =1.0, 0.9, sparse, and tensor product structures with rotor in Figure 3. The observed scaling of effort for a fixed structure is similar for both the rotor (Eq. (27)) and oscillator (Eq. (28)) structures, as seen in Table II. The search effort scaling was found to be strongly dependent on the control Hamiltonian structure. For =1.0 or sparse , i.e., class that allows distant transitions between states, the effort scales slowly with . In contrast, for the 1.0 and tensor product structure (Eq. (30)), i.e., the class that forbids or has weak distant transitions, the effort scales exponentially with , as shown by the least-square fit lines on the semi-log plot for =0.9 and tensor product structures in Figure 3.

The fluence of the initial field has some effect on the absolute search effort, but not on its scaling with (columns 4 and 5 of Table II). Increasing the fluence cannot overcome exponential scaling for the class of control Hamiltonians that forbids distant transitions. The effort can be reduced to some extent by allowing to scale with (Table II), in agreement with the conclusion that longer control times are needed for systems that have few accessible control pathways [33], but the effort still scales exponentially with . The choice of (i.e. random unitary or FT gate) does not greatly affect the search effort, as shown by comparing the effort to find the FT gate or a random using the sparse structure.

6.c Control Hamiltonian Structure and Search Effort

Of all the search parameters explored above, only the control Hamiltonian structure has a systematic effect on the scaling of the search effort with . In order to determine the effects of the structure of for fixed (here =8), we compared control Hamiltonians with =1.0, 0.9, 0.75, and 0.6, randomly generated sparse structures with 14, 10, or 8 allowed transitions, and banded control Hamiltonian structures where 2, 3, or 4 rows nearest to the diagonal contain allowed transitions. The resulting search effort is plotted versus the norm in Figure 4, which clearly shows that does not determine search effort. Rather, the distribution of strong couplings between states is important.

For a given value of , the banded structure has the greatest search effort, followed by the structure, and the sparse structure has the smallest search effort. For approximately the same value of , the search effort varies by over a factor of 10, from 70 iterations for the sparse structure (14 transitions), through 100 iterations for =0.75 structure, to 1200 iterations for the banded structure with 2 rows of allowed transitions (13 transitions). This indicates that systematic exclusion or suppression of transitions between distant states while allowing only transitions between near states raises the search effort, compared to having an equal number of allowed transitions that include some distant transitions. Thus, the control Hamiltonian class that allows distant transitions is expected to have a lower search effort than the class that forbids distant transitions. This observation is consistent with the observed exponential scaling of the search effort with for 1 and tensor product structures. The reasons behind the dependence of the search effort on the control Hamiltonian structure will be explored in Section 7.

7 Search Effort and Landscape Geometry

Here, we assess the local landscape features in terms of the metrics in Section 4 for unitary propagator control. We first consider the structure of the landscape in terms of the local metrics in Section 7.A. The effects of the landscape structure on the search trajectories, as defined by the directness metric and the Gramian matrix, are examined in Section 7.B.

7.a Local Landscape Structure

The bound on the slope metric derived in Section 4 was found to be conservative. The recorded maximal slope metric was always significantly below this bound and observed to grow linearly with , while the bound grows quadratically with (not shown). The analytically derived Hessian trace at =0 (c.f. Eq. 20) was found to hold; the deviation at the optimum was always less than 0.001 when the convergence criterion was used. The Hessian trace does not predict search effort regardless of where it is measured, since it is only dependent on or , and the effort can vary widely for different structures with similar values of (c.f., Figure 4).

The slope metric and the local curvature were found to correlate with search effort. The statistical distribution of the maximal slope metric over the search samples is plotted versus in Figure 5(a). For =1.0, the maximal slope rises linearly with , but the growth with is slower for 1.0 and tensor product control Hamiltonians. Since the maximal slope metric for any search is often recorded at or near the initial field (depending on the exact choice of ), an estimate of search effort scaling with for any structure can be made simply by measuring the gradient at random initial fields for systems of different with the same type of control Hamiltonian structure: a linearly increasing with indicates minimal scaling with , while sub-linear increase of indicates exponential scaling with . A more accurate prediction of the search effort can be made by measuring the near, but not at, the optimum. Figure 6 shows the absolute search effort for searches using different structures plotted versus the measured value of at =2 and =0.01. Even at =2, the gradient provides a good estimate of the absolute effort. The local curvature at the optimum is also an indicator of search effort scaling, as shown in Figure 5(b). For control Hamiltonians that allow distant transitions and have sub-exponential scaling of effort with (=1.0 and sparse structures), the curvature is flat as increases from 4 to 16. For control Hamiltonians that forbid distant transitions and have exponential scaling (1.0 and the tensor product structures), the curvature decreases with according to a power law (note the log-log plot). Under all circumstances, a smaller value of near the optimum indicates a greater search effort.

To understand the effects of Hessian local curvature on convergence efficiency, note that near a stable fixed point of on the global optimum submanifold, (i.e., for ), the objective function is approximately quadratic in , and we can linearize the differential equation (15) around as