{'fraction_non_alphanumeric': 0.0491005291005291, 'fraction_numerical': 0.018042328042328044, 'mean_word_length': 4.181195175438597, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tensor product approach to quantum control

In this proof-of-concept paper we show that tensor product approach is efficient for control of large quantum systems, such as Heisenberg spin wires, which are essential for emerging quantum computing technologies. We compute optimal control sequences using GRAPE method, applying the recently developed tAMEn algorithm to calculate evolution of quantum states represented in the tensor train format to reduce storage. Using tensor product algorithms we can overcome the curse of dimensionality and compute the optimal control pulse for a 41 spin system on a single workstation with fully controlled accuracy and huge savings of computational time and memory. The use of tensor product algorithms opens new approaches for development of quantum computers with 50 to 100 qubits.


Introduction
Quantum control plays a central part in both established technologies such as nuclear magnetic resonance (NMR) and emerging technologies such as quantum computing [16]. The 'hardware' of a quantum computer is a quantum system (e.g. a spin system), which we manipulate with time-dependent fields (e.g. laser pulses). Our goal is to design a pulse sequence that steers the system from an initial state |ψ 0 to a desired state |ψ T during time T -a 'firmware', so to say. The main challenge for numerical methods is the notoriously known curse of dimensionality -the dimension of the state space growing exponentially with the number of particles d. Even for simplest spin- 1 2 systems, the wavefunction |ψ , the density matrix ρ, and the HamiltonianĤ require at least O(2 d ) storage. To find the optimal trajectory we need to compute ρ(t) in 100s to 1000s points of time. This leads to gigabyte-scale storage for d ≃ 20 spins and becomes unfeasible for d 30.
The main approaches used currently are the sparse format forĤ, |ψ 's, ρ's (only large entries are stored, zeroes and small elements are discarded), and dimensionality reduction -the equations of system dynamics are projected into a subspace of low-lying eigenstates [44] or its heuristically chosen equivalent [27,4].
Tensor product algorithms are a novel approach to high-dimensional problems. Based on the idea of separation of variables, they approximate high-dimensional data with a tensor product of low-dimensional factors. For example, when interaction of particles is a superposition of a few near-neighbour interactions, the corresponding quantum state is a superposition of a few unentangled states, i.e it has a low Schmidt rank. The corresponding tensor format also has a low rank, providing good compression and low storage costs. By maintaining the compressed format for all operations we can compute the system's trajectory and optimise it without any uncontrollable truncations or heuristic assumptions. Model examples of such systems include Ising and Heisenberg spin chains [35], more realistic applications include quantum wires [3] and backbones of simple protein molecules [34].
In quantum physics a tensor product format was first proposed in 1970s as the renormalization group formalism [43], followed in 1990s by MPS/MPO [14,26] and DMRG [42] algorithms for finding the ground state of a quantum spin chain. Similar formats were applied in statistics and can be traced back to works of Hitchcock in 1930s [22,21] and Tucker in 1960s [40]. Tensor formats in three dimensions were rigorously studied in the numerical linear algebra community as a generalisation of the low-rank decomposition of matrices [30,15,18,33]. Generalisations to higher dimensions eventually led to re-discovery of the MPS/DMRG framework under the name of the tensor train (TT) format [31,29]. In the last decade the TT format and the more general Hierarchical Tucker format [20] have been successfully applied to a variety of problems, such as finding several low-lying states of a quantum system [8], performing superfast Fourier transform [9], solving high-dimensional linear systems [10] and first-order differential equations [7]. First applications of tensor product methods to optimal control problems are considered in [6,1].
The recently proposed tAMEn algorithm [7] allows to compute evolution of a quantum system under a given (time-dependent) Hamiltonian with high accuracy. Although the trajectory lies in extremely high-dimensional state space, tAMEn keeps the computation costs feasible by performing all calculations in the TT format and adapting the TT ranks to meet the required accuracy. In this paper we employ the tAMEn algorithm as a building block to develop a version of a classical GRAPE algorithm [25] and use it to control a spin chain with d = 41 spins. This shows that tensor product algorithms can be used to design control sequences for quantum computers with 50 to 100 qubits, which can be expected within a decade [19].

Dynamic optimisation problem
The dynamics of a quantum system with HamiltonianĤ is described by Liouvillevon Neumann equation in superoperator forṁ where ρ is the quantum density matrix stretched into a vector form (referred to later as quantum state), andĤ TraditionallyĤ(t) is split into a constant drift and a time-dependent control term: Control operatorsĤ k are usually defined by model or instrument used; we manipulate the system by changing control pulse functions c k (t). Now choosing a suitable basis set {ϕ n (t)} N n=1 and dropping the relaxation partR for simplicity, we obtain an ODE for system dynamics: Our goal is to choose control parameters c = [c n,k ] to maximise fidelity (or overlap)

First-order optimisation framework
Classical optimisation methods require a gradient of F w.r.t. c n,k which reads To find the gradient, we differentiate (1) w.r.t. control parameters and obtain ∂ ∂c n,k ∂ρ ∂t = −i ϕ n (t)Ĥ k ρ +Ĥ(t) ∂ρ ∂c n,k .
Changing the order of derivatives in the left hand side, we obtain a system of coupled ODEs for the density matrix and its gradient: To justify the initial condition for the gradient, write a first-order approximation for ρ(δt) with an infinitesimal time δt: and differentiate it w.r.t. c n,k to see ∂ρ ∂c n,k (δt) = O(δt) → 0. The system (4) can be integrated on [0, T ] using a suitable method (e.g. Runge-Kutta scheme [32]) for all k, n to produce both the fidelity (2) and its gradient (3).

GRAPE algorithm
A classical gradient ascend pulse engineering (GRAPE) method [25] is a simple version of the approach described above. Control pulses c k (t) are assumed piecewise-constant on intervals (t n−1 , t n ] between the nodes of the grid {t n } N n=0 with t 0 = 0 and t N = T. We will use a uniform grid t k = τk with step-size τ = T/N. The basis functions ϕ n (t) can be taken as Heaviside functions thus making c n,k the amplitude of control pulse c k (t) on (t n−1 , t n ]. The time-locality of basis functions θ n (t) allows us to compute (3) more efficiently. Note that θ n (t) = 0 for t t n−1 hence ODE (4) for the gradient component simplifies to the homogeneous ODE which of course gives ∂ρ ∂c n,k (t n−1 ) = 0. Therefore, to compute (3) we can propagate the state ρ(t) on (0, t n−1 ] and start solving the coupled system (4) only from t n−1 . When we reach t n and compute ∂ρ ∂c n,k (t n ), we can note that θ n (t) = 0 for t > t n hence (5) holds on t ∈ (t n , T ] as well. SinceĤ(t) is piecewise-constant, and plugging this in (3) gives The state is the solution of ODE with the boundary condition at the right end, colloquially referred to as 'back-propagation': Hence the coupled system (4) needs to be solved only on [t n−1 , t n ] to produce the n-th components of the gradient. By pre-computing all λ(t n ) for n = 1, . . . , N, we can implement the GRAPE iteration with complexity scaling linearly in N.
1. For n = N, . . . , 1, propagate (7) backward: 3. Compute (6) and update the control amplitudes: c n,k −→ c n,k − ǫ ∂F ∂c n,k . Here ǫ is the step size for the gradient ascend. It should be carefully selected, since a too small value will make the algorithm converging very slowly while a very big one will make the algorithm divergent. In our implementation the step size is chosen adaptively: it is increased each time the step is successful (i.e. results in better fidelity) and decreased each time we have an overstep (i.e. fidelity decreases and the step has to be rejected). Note that line-search strategies for this optimisation problem are rather pointless, because the gradient (6) and the fidelity function (2) are computed simultaneously from (4), and have therefore the same complexity.
In the classic GRAPE algorithm, the propagators Q m in (6) are calculated as matrix exponentials [28] and applied to the states; in more advanced versions of the algorithm the action of matrix exponentials on the states is computed on-fly maintaining the sparse format for λ's and ρ's [36].
Instead of computing propagators, we use the tAMEn algorithm to solve the first-order ODEs (4) and (7) with Hamiltonians and states represented in the TT format. The next section focuses on mathematical details of this approach.

Tensor train format and the tAMEn algorithm
Introducing individual state indices for each site, we can treat ρ as a multi-index (high-dimensional) array which is referred to as tensor in numerical linear algebra: We should admit that our use of the term tensor does not imply proper differentiation of upper and lower indices or Einstein summation convention [11]. We use the Tensor Train (TT) decomposition [29] to separate pairs of i k , j k , belonging to different sites, by a product approximation of the form Here, ρ (k) , k = 1, . . . , d, are called TT blocks, and the summation ranges r 1 , . . . , r d are called TT ranks. The grouping of i k and j k is such that unentangled (e.g. pure) states are represented by (8) exactly with elementary TT ranks r 0 = · · · = r d =  (1), (4) and (7) in time. Recall that we solve ODEs on time intervals t ∈ (t n−1 , t n ], where control pulses c k (t) and henceĤ(t) are constant and only ρ(t) depends on time. Without loss of generality we can assume that t ∈ (0, τ] and drop the time index n for the rest of this section. We choose a basis of Lagrange polynomials {L m (t)} M m=1 centered at Chebyshëv nodes {τ m } M m=1 and represent the state as Since the length of the interval τ is usually quite small to provide better resolution of control pulses, the basis size M can be very moderate. We collect the nodal values ρ(τ m ) into a vector ρ = [ρ(τ m )] M m=1 of length 4 d M, which we aim to find. The time derivative operator is discretised as a differentiation matrix S = [L ′ ℓ (τ m )] M m,ℓ=1 in accordance to the spectral approximation theory [39]. Applying this discretisation scheme to ODE (1) we obtain the following linear system Here ρ is the unknown vector of ρ(τ m )'s, 1 M is the vector of size M full of 1's, and I denotes identity matrices of appropriate size. Note that the Kronecker products are implicitly realised by the TT representation (10). Therefore, we never actually compute them in (9), but just collect the corresponding factors into TT blocks. The same applies to the states. In particular, we consider unentangled states ρ 0 and ρ T , and construct their TT blocks explicitly. Simple mixed states can be made by concatenation of TT blocks instead of a direct summation of a large number of density matrices [34]. The TT ranks can be truncated down to optimal values for a desired accuracy using the singular value decomposition (SVD) [29]. Algebraic equations (9) can be solved in the TT approximation by an alternating iteration, such as the Alternating Minimal Energy (AMEn) algorithm [10], which is an enhanced version of Alternating Least Squares (ALS) [23]. We expand the TT representation (8) by an extra block which carries the time index m and write: Now the last block ρ (d+1) = [ρ (d+1) α d ,α d+1 (m)] carries degrees of freedom in time, the rightmost TT rank is r d+1 = 1, and r d can now be larger than 1.
The basic ALS algorithm would solve (9) by sweeping over TT blocks, i.e. representing ρ by (10) and reducing the equation to only one block ρ (k) in each step. This can be written using a frame matrix , which is of size 4 d M × 4r k−1 r k for k d, and of size 4 d M × r k−1 Mr k for k = d + 1.
The frame matrix realises a linear map from the elements of a single TT block into the elements of the whole vector given (10), that is ρ = P =k ρ (k) , assuming that elements of ρ (k) are stretched in a vector. The ALS method solves the reduced Galerkin systems (P † =k AP =k )ρ (k) = P † =k f subsequently for k = 1, · · · , d + 1. This system can be assembled efficiently [23] due to the TT formats of P =k , A and f and solved using standard methods.
The AMEn method [10] improves the convergence of ALS by expanding ρ (k) (and hence all P =q for q > k) with a TT approximation of the residual f − Aρ, where ρ is formed by (10) with an updated block ρ (k) plugged in. This also allows to adapt the TT ranks to ensure a desired accuracy of the TT approximation. The tAMEn (time-dependent AMEn) algorithm [7] utilises the special meaning of the last TT block, carrying the time variable m, in order to preserve conservation laws (e.g. the Frobenius norm), and to estimate the time discretisation error.
For the interval lengths τ chosen in our numerical experiments, we have found that M = 8 Chebyshëv nodes in each interval are sufficient to resolve the time derivative with the relative accuracy of 10 −5 or better. In each step the tAMEn algorithm produces a set of TT blocks ρ (1) , . . . , ρ (d+1) , comprising the discrete-time solution ρ on the interval (0, τ]. After ρ is obtained, we compute ρ(τ) in the TT format by taking the slice of the last TT block ρ (d+1) (M) and hence removing the time index m from consideration. Since t = τ represents the end point of the interval (t n−1 , t n ] the obtained state represents ρ(t n ) and can be used now as an initial state in the next GRAPE step over the time interval (t n , t n+1 ]. The tAMEn algorithm is agnostic to a particular form ofĤ (the only assumption is that it must admit a TT decomposition), and can be applied to any tensorstructured differential equationẋ = Ax such as the auxiliary matrix formalism (4) and the back-propagation (7).

Numerical experiments
We consider a Heisenberg chain of spin-1 2 particles, for which a Hamiltonian is a simple sum of nearest neighbour interactionŝ where the Pauli matricesσ When J x = J y = J z , this system is called the XXX Heisenberg model, when J x = J y = J z it is called the XXZ model, and for J x = J y = J z it is called the XYZ model. It is commonly mentioned that linear Heisenberg chains can be diagonalised exactly using the Bethe ansatz [2], however the eigenstates can only be written via roots of a system of algebraic equations of degree d [37], that are not computable in closed form neither using analytic methods nor numerically with reasonable accuracy for d 20. For simpler XX models (J x = J y , J z = 0) eigenvectors are available in closed form and can be used for dimensionality reduction, allowing chains with d ≃ 200 spins to be controlled [44]. We could not find examples of optimal control pulse computed numerically for XXX, XXZ and XYZ Heisenberg models with the number of spins d 20. We decided therefore to test our tensor product approach for XXX and XXZ Heisenberg spin-1 2 chains with d = 11, d = 21 and d = 41. The initial and target states are taken as ψ 0 = |↑↓↓ · · · ↓↓ , ψ T = |↓↓ · · · ↓↓↑ , so the task is to move the |↑ state from the first to the last position in the chain. The control operator H c (t) = c(t)σ (1) z is the magnetic field acting on the first spin only.
In theory, Heisenberg chains are fully controllable [41], which means the fidelity (2) can be made infinitely close to 1 as T → ∞. In practice, however, the available time T for the pulse is limited, and the final state ρ(T ) will not be fully focused, leaving some infidelity 1 − | ρ T | ρ(T ) | > 0. We optimise the pulse sequence to reduce the infidelity as much as possible.
The results are shown on Fig. 1. Computing a short high-fidelity pulse for large spin system appears to be a challenging task. Using 4 cores on a Xeon E5-2650 CPU and 10 GB of memory, we reduced infidelity to 10 −3 for d = 21 and to 10 −2 for d = 41 within a few days of computation. This seems fast compared to high-fidelity calculations reported in [38], where infidelity 10 −4 was reached after weeks/month of optimisation for a system with d = 4 qubits with 4 levels each.
The efficiency of the method depends primarily on reasonable choice of the pulse length T, which is usually of the order of the inverse of the natural frequency of the system [24]. The time T typically should increase with the number of spins in the chain due to the finite rate of the information exchange among them. We see in our experiments that larger intervals T are required to reduce infidelity for larger chains; also, XXZ chains demand larger T than XXX chains. We keep resolution τ = 10 −2 for the control pulse constant for all experiments, which is the main reason why computational time per iteration grows with T. The TT ranks grow only mildly with d and remain moderate r k 30 for d = 41. This resulted in about 70 · 10 3 unknowns in the TT decomposition of each ρ(t k ) for d = 21 and 140 · 10 3 unknowns for d = 41. For comparison, a single diagonal of ρ(t k ) has 2 21 ≈ 2 · 10 6 entries for d = 21 and 2 41 = 2 · 10 12 entries for d = 41.

Conclusions and future work
Tensor product algorithms open new possibilities to control long spin wires and multi-qubit gates in quantum computers. We explored a proof-of-concept example, in which we considered a simple XXX or XXZ Heisenberg chain, optimised control pulse using a classical GRAPE method [25], while representing all states n the TT format [29] and propagating them using the tAMEn algorithm [7]. This combination of relatively simple algorithms allowed us to reach fidelity of 99% for a chain of d = 41 spins using a single workstation for calculations. We used only controllable approximation techniques and avoided any heuristic or random truncations of the state space. Our algorithm is deterministic and flexible, i.e. can be applied to any linear or quasi-linear quantum system, for which we can expect the states to have moderate entanglement (as measured by Schmidt ranks).
We have not yet reached the fidelity of ∼ 99.99% required for topological error correction [38], mostly because of slow first-order convergence of GRAPE, but we hope to make it possible using a second-order optimisation algorithm, such as Newton-Raphson [17] or a quasi-Newton e.g. BFGS [5, 13]. We have not imposed any constraints on pulse shape, which resulted in a few outbreaks seen on Fig. 1, but implementation of the box constraints |c(t)| < C is a straightforward step which we postpone for future work.