Escaping barren plateaus in approximate quantum compiling

Niall F. Robertson$^{1}$, Albert Akhriev$^{1}$, Jiri Vala$^{2,3}$ and Sergiy Zhuk$^{1}$

$^1$ IBM Quantum, IBM Research Europe - Dublin, IBM Technology Campus, Dublin 15, Ireland

$^2$ Maynooth University, Maynooth, Ireland

$^3$ Tyndall National Institute, Cork, Ireland

Quantum compilation provides a method to translate quantum algorithms at a high level of abstraction into their implementations as quantum circuits on real hardware. One approach to quantum compiling is to design a parameterised circuit and to use techniques from optimisation to find the parameters that minimise the distance between the parameterised circuit and the target circuit of interest. While promising, such an approach typically runs into the obstacle of barren plateaus - i.e. large regions of parameter space in which the gradient vanishes. A number of recent works focusing on so-called quantum assisted quantum compiling have developed new techniques to induce gradients in some particular cases. Here we develop and implement a set of related techniques such that they can be applied to classically assisted quantum compiling. We consider both approximate state preparation and approximate circuit preparation and show that, in both cases, we can significantly improve convergence with the approach developed in this work.

§ INTRODUCTION

Quantum compiling is a set of tools to translate a quantum circuit into a set of gates that can be efficiently implemented on quantum hardware. To do so, one must take into account hardware constraints such as connectivity and maximum circuit depth. Limited hardware connectivity between qubits requires additional operations such as $SWAP$ gates to be included into the compilation process [1].  Similarly, the physical implementation of entangling two-qubit operations such as a controlled-NOT ($CNOT$) operation, illustrated in Fig. <ref>, may require that only one of the two qubits can play the role of the control qubit. Despite recent progress on error mitigation techniques [2, 3, 4], noise on present day devices significantly limits the depth of the circuit that can be implemented - the benefit of optimally decomposing a quantum circuit into the hardware's native gate alphabet and connectivity constraints is hence self-evident. We can classify quantum compiling techniques into two broad categories: classically assisted and quantum assisted. In the former case, the idea is to use classical methods to compile sub-circuits that form part of a larger quantum algorithm. The main advantage of such an approach is that one can avoid the issues arising from noise on a quantum device, the disadvantage being that, due to the exponential scaling of the Hilbert space, one is limited to sub circuits that act on a fixed number of qubits. The quantum assisted approach however is scalable and can hence be applied to the entire quantum circuit of interest. However, one encounters quantum noise which, even in cases where the compilation scheme is somewhat noise resilient [5], can cause significant convergence issues [6]. Furthermore, the circuits that have been proposed to implement a quantum assisted quantum compilation procedure require long range qubit connectivity on the quantum device of interest [7].

Here we focus on the classically assisted approach to quantum compilation. Classically assisted variational quantum algorithms have been recently studied in [8, 9]. In these works, Tensor Networks are used to generate a shallow circuit with weak entanglement. The circuit is subsequently deepened on a quantum device, and the entanglement is thus increased beyond what classical techniques can store. These recent promising results provide significant motivation to improve the state of the art of classical optimisation techniques applied to quantum circuits which are currently limited by convergence issues such as barren plateaus, as will be discussed below.

We consider the problem of approximate quantum compiling as studied in [10]. This approach involves the design of a parametric quantum circuit with fixed depth - the parameters are then adjusted to bring it as “close" as possible to the target, where “close" is defined via some carefully chosen metric (see section <ref>). In [10], it was shown that such an approach could be used to derive a 14-$CNOT$ decomposition of a 4-qubit Toffoli gate from scratch. It was furthermore shown that the Quantum Shannon Decomposition, introduced by Shende, Bullock and Markov [11], can be compressed by a factor of two without practical loss of fidelity. In this work, our main contribution is a new classical algorithm for quantum compilation which significantly improves convergence - see below and main text for details.

In addition to the sub-circuit problem discussed above, we consider state preparation as an example where classical techniques can play a significant role. For a large class of quantum algorithms, the entanglement of the wavefunction is weakly entangled at the beginning of the circuit and grows with the depth. There exist widely studied classical techniques (i.e. Tensor Networks) to store the wavefunction of a large number of qubits when the system is only weakly entangled. Thus, this opens the possibility of using classical Tensor Network techniques for state preparation to optimally compile the early parts of a quantum circuit, in a similar way to [8, 9]. One of the biggest obstacles however to any quantum compilation technique, whether it is quantum assisted or not, is the existence of barren plateaus in the parameter landscape, i.e. large regions where the gradient of the cost function vanishes exponentially in the number of qubits, thus rendering it extremely difficult to train. A proposed solution for the quantum assisted quantum compiling approach is to use shallow circuit Ansatz' and so-caled “local cost functions". In short, a cost function calculated via a quantum circuit is local if it involves the measurement of only a small number of qubits - a global cost function involves the measurement of all qubits. Here we design a new classical algorithm inspired by the distinction between local and global cost functions made in [7, 12] for quantum compilation. In particular, we derive an expression inspired by the quantum approach that can be easily implemented on a classical computer and that allows us to escape the barren plateaus. This new cost function introduces a number of coefficients that we must update throughout the optimsation procedure. Our algorithm implements a scheme to do so effectively, such that short depth parametric quantum circuits can be found that closely approximate the target circuit of interest.

Our main contributions in this work are as follows: we derive an explicit expression for the local cost functions considered in [7, 12] in equation (<ref>), such that they obtain a clear and intuitive operational meaning. We use this to develop a new cost function - see equation (<ref>) - whose gradient can be efficiently calculated classically. We consider a particular case where our cost function can be studied analytically and we show how the variance of its gradient behaves - see equation (<ref>) - which allows us to understand how and why our classical compilation scheme allows for the escape from barren plateaus. We develop an algorithm to update the coefficients appearing in our new cost function - see section <ref>, appendix <ref> and appendix <ref>. Finally, we demonstrate numerically in Figures <ref> and <ref> that the combination of our cost function and our algorithm to update the coefficients significantly improves the convergence of certain problems in quantum compilation.

The outline of this article is as follows: in section <ref>, we review approximate quantum compiling and we show how it can be formulated as an optimisation problem. In section <ref>, we discuss the convergence issues that are commonly encountered in quantum compiling problems - we discuss in particular the appearance of barren plateaus in the parameter landscape. We introduce our classical version of the local cost function and we show how, for a particular example where the cost function can be studied analytically, the variance of the gradient of our cost function behaves as a function of the number of qubits. In section <ref>, we describe our scheme to update the coefficients that appear in our classical local cost function over the course of the optimisation procedure. In section <ref>, we present some numerical results on simulations for up to 24 qubits showing that our classical local cost function significantly improves convergence as compared to the global one. We conclude in section <ref>.

§ SETUP

Quantum compilation is an active research area with a wide spectrum of research directions. One such direction is to find the the best representation - in terms of some distance or structural properties - of a generic unitary in a certain canonical basis (e.g. [13, 14, 15, 16, 17]). One can also consider an approach based on exact qubit routing, i.e. mapping logical qubits to physical qubits so that the resulting unitary coincides with the original one (up to a permutation) while allowing for compatibility with the hardware topology (e.g. [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]). As discussed in the introduction, the approach to quantum compilation that we focus on here is to design a parametric quantum circuit with fixed depth and to adjust the parameters of this circuit to bring it as “close" as possible to the target, where “close" is defined via some carefully chosen metric. While one can in principle use universal parametric circuits which provide a theoretical guarantee that there exist parameters such that any target circuit can be exactly compiled, such circuit templates result in impractically deep circuits for present day quantum devices. Approximate quantum compiling can be seen as a compilation approach that trades off the circuit accuracy for its depth. It was shown in [32] that generic unitary operations over $n$ qubits are computationally hard even to approximate. More recently it has been shown that quantum circuit compilation is an NP-complete problem [33]. However, the fact that there is a subset of unitary operations that can be realized in quantum circuits of lower complexity, such as those in the Shor algorithm [34] for example, suggests that quantum compilation may offer both an interesting perspective on complexity of quantum circuits and perhaps a path to new quantum algorithms that would provide quantum advantage.

[scale = 1.5]
[gray, thick] (-0.3,1) – (0.3,1);
[gray, thick] (0,1) – (0,0);
[black] (0,1) circle (2pt);
[black] (0,0.12) circle (4pt);
[gray, thick] (-0.3,0.13) – (0.3,0.13);

[gray, thick] (4,1) – (4,0);
[gray, thick] (3.7,0) – (4.3,0);
[black] (4,0) circle (2pt);
[black] (4,0.88) circle (4pt);
[gray, thick] (3.7,0.87) – (4.3,0.87);

Controlled-NOT ($CNOT$) operations with the role of a control qubit played by the first and second qubit respectively. $CNOT$ transforms each standard computational basis element in a quantum superposition by adding a value of the control qubit to that of the target qubit modulo 2: $\ket{a}_c \otimes \ket{b}_t \rightarrow \ket{a}_c \otimes \ket{b \oplus a\mod2}_t$ where $a,b \in \{0, 1\}$.

§.§ Parametric quantum circuits and ansatz

As discussed in [10], a natural circuit Ansatz is one based on so-called $CNOT$ blocks. A $CNOT$ block is a $CNOT$ gate followed by single qubit rotations (see Figure <ref>). Single-qubit rotation gates are defined as complex exponential functions of the Pauli matrices:

R_x(θ)=exp(-i σ_x θ/ 2)=[
cos(θ/ 2)     -i sin(θ/ 2)

-i sin(θ/ 2)     cos(θ/ 2)

R_y(θ)=exp(-i σ_y θ/ 2)=[
cos(θ/ 2)     -sin(θ/ 2)

sin(θ/ 2)     cos(θ/ 2)

R_z(θ)=exp(-i σ_z θ/ 2)=[
exp(-i θ/ 2)     0

0     exp(i θ/ 2)

[scale = 1.8]
[gray, thick] (0,1) – (0,0);
[black] (0,1) circle (2pt);
[black] (0,0.12) circle (4pt);
[gray, thick] (-0.14,0.13) – (0.14,0.13);

[gray, thick] (-1.14,1) – (1.1,1);
[gray, thick] (-1.14,0.13) – (1.1,0.13);

[gray, thick] (1.1,0.7) – (1.1,1.3);
[gray, thick] (1.9,0.7) – (1.9,1.3);
[gray, thick] (1.1,1.3) – (1.9,1.3);
[gray, thick] (1.1,0.7) – (1.9,0.7);

[gray, thick] (1.1,-0.17) – (1.1, 0.43);
[gray, thick] (1.9,-0.17) – (1.9, 0.43);
[gray, thick] (1.1, 0.43) – (1.9, 0.43);
[gray, thick] (1.1,-0.17) – (1.9,-0.17);

[gray, thick] (1.9,1) – (2.3,1);
[gray, thick] (1.9,0.13) – (2.3,0.13);

[gray, thick] (2.3,-0.17) – (2.3, 0.43);
[gray, thick] (3.1,-0.17) – (3.1, 0.43);
[gray, thick] (2.3, 0.43) – (3.1, 0.43);
[gray, thick] (2.3,-0.17) – (3.1,-0.17);

[gray, thick] (2.3,0.7) – (2.3,1.3);
[gray, thick] (3.1,0.7) – (3.1,1.3);
[gray, thick] (2.3,1.3) – (3.1,1.3);
[gray, thick] (2.3,0.7) – (3.1,0.7);

[gray, thick] (3.1,1) – (3.5,1);
[gray, thick] (3.1,0.13) – (3.5,0.13);

at (-1.3,1) j;
at (-1.3,0.13) k;

at (1.5,1) $R_y(\theta_1)$;
at (2.7,1) $R_z(\theta_2)$;
at (1.5,0.13) $R_y(\theta_3)$;
at (2.7,0.13) $R_x(\theta_4)$;

$CNOT$ block forms the basic building block of our circuit ansatz.

R_z(θ_1)     R_y(θ_2)     R_z(θ_3)     [wires=2]    [wires=2]   [wires=2]                 [wires=2]     [wires=2]     [wires=2]

R_z(θ_4)     R_y(θ_5)     R_z(θ_6)                [wires=2]    [wires=2]    [wires=2]                [wires=2]     [wires=2]     [wires=2]

R_z(θ_7)     R_y(θ_8)     R_z(θ_9)     [wires=2]    [wires=2]   [wires=2]                [wires=2]     [wires=2]     [wires=2]

R_z(θ_10)     R_y(θ_11)     R_z(θ_12)

Parameterised circuit with $n=4$ qubits, $l=2$ layers and $b=3$ CNOT blocks in each layer.

A block with a $CNOT$ gate acting on a “control" qubit $j$ and “target" qubit $k$ is written as $\mathrm{CU}_{jk}(\theta_1, \theta_2, \theta_3, \theta_4)$. For a given hardware connectivity, one can then write down a full parameterised circuit as:

V_ct(θ)=    CU_ct(L)(θ_3 n+4 L-3, …, θ_3 n+4 L) ⋯CU_ct(1)(θ_3 n+1, …, θ_3 n+4)

[R_z(θ_1) R_y(θ_2) R_z(θ_3)] ⊗⋯⊗[R_z(θ_3 n-2) R_y(θ_3 n-1) R_z(θ_3 n)]

The position of the $CNOT$ blocks in the parameterised circuit can be customised to suit the particular target circuit that one is interested in. We will use a particular structure of the ansatz which is illustrated in Figure <ref>. In this circuit Ansatz, we have 3 rotation angles for each qubit at the beginning, we then have $l$ layers in which we repeat each $CNOT$ block $b$ times. In Figure <ref>, we have $n=4$ qubits, $l=2$ and $b=3$. We note that for $b=3$, each block of $CNOT$s acting on $2$ qubits constitutes a universal two-qubit quantum circuit which represents an arbitrary unitary rotation of a two-qubit quantum state.

For an $n$-qubit circuit, we typically take the initial state to be $\ket{0}^{\otimes n}$ which, in the illustrated example, is first transformed by a sequence of three single-qubit rotations $R_z(\theta_1)$, $R_y(\theta_2)$, $R_z(\theta_3)$. Strictly speaking, only two inequivalent single-qubit rotations are required if the initial state is the fiducial initial state $\ket{\boldsymbol{0}}$ but we use three for implementation in order to consider more general situations [10].

§.§ Optimization problem

The compilation now becomes an optimisation problem whereby one aims to tune the parameters of the rotation gates to minimise the distance between the parametric circuit and the target. Given a target unitary matrix $U \in \mbox{U}(2^n)$ and a set of constraints, such as connectivity and length , the approximate quantum compiling problem is to find the unitary matrix $V(\theta) \in \mbox{U}(2^n)$, that can be implemented with parametrized quantum circuit and that is closest distance to the target $U$. One can define this distance via a number of metrics such as the Frobenius norm

min_V ∈U(2^n) f(V):=1/2V(θ)-U_F^2 = 1 - 1/2^n  [(V^†(θ) U) ]

Similarly, one can also use the Hilbert-Schmidt test:

C_HS = 1 - 1/d^2|(V^†(θ)U)|^2

where $d$ is the dimension of the Hilbert space. The Hilbert-Schmidt test (<ref>) is closely related to the Frobenius norm in equation (<ref>) and used as the cost function in [7]. We have:

C_FN≡1/2V(θ)-U_F^2 = d - Re [(V^†(θ)U)]

From (<ref>) and (<ref>), we get:

1-C_HS = (1-C_FN/d)^2 + 1/d^2(Im ((V^†U)))^2

and hence:

1-C_HS ≥(1-C_FN/d)^2

We can thus conclude that minimising $C_{FN}$ minimises $C_{HS}$ automatically. We will show in this work that one can improve the convergence issues caused by barren plateaus in approximate quantum compiling by minimising a cost function that is related to (<ref>) above but also includes some additional perturbation terms that “localise" the original cost function - a concept that will be precisely defined below.

§ BARREN PLATEAUS

The Frobenius norm in equation (<ref>) provides a natural cost function for Approximate Quantum Compiling. In particular, one can indeed see that this quantity is $0$ if and only if $V(\theta) = U$. However, this cost function is problematic in practice and is only feasible for a small number of qubits $n$ due to the appearance of barren plateaus, i.e. large regions of parameter space where the gradient of the cost function vanishes exponentially, thus rendering it impossible to train the cost function. Not only is this observed in practice, but a number of recent works [12, 6, 35] have even provided strict proofs that guarantee the existence of these barren plateaus in a variety of circumstances, where the definition of a barren plateau for these proofs is that the variance of the gradient vanishes exponentially with the number of qubits $n$. We remark that these works were focused on quantum assisted quantum compiling in which case an exponentially vanishing gradient requires an exponential number of measurements to estimate the gradient. In the classical implementation of quantum compiling that we consider here, we calculate the gradient directly without taking measurements and thus in principle an exponentially vanishing gradient is not as significant of an issue. However, in practice we observe that even in the classical case, the appearance of these barren plateaus cause the training times to increase beyond practical limits. Furthermore, we observe that the barren plateaus appear even in systems with small numbers of qubits $n$. We hence require new tools to overcome these issues.

In [7, 12], it was shown that one can escape the barren plateau with so-called local cost functions as opposed to global cost functions - to be defined below - as long as the circuit ansatz $V(\theta)$ is sufficiently shallow. In particular, as shown in [12], if the depth of $V(\theta)$ is of order $\mathcal{O}(\log n)$ then the variance of the gradient of a local cost function vanishes only polynomially instead of exponentially. Broadly speaking, a global cost function is one that is evaluated by measuring all qubits at once, whereas a $k-$local cost function is evaluated by measuring only $k$-qubits at once. More specifically, in the context of state preparation circuits as studied in [12], a global cost function is defined as one of the form:

C_G=Tr[O_G V(θ)|ψ_0⟩⟨ψ_0| V(θ)^†]

where $O_{\mathrm{G}}$ is an operator that acts non-trivially on all qubits and $\ket{\psi_0}$ is the target state. An example would be:

A local cost function, on the other hand, is defined as:
C_L=Tr[O_L V(θ)|ψ_0⟩⟨ψ_0| V(θ)^†]
where $O_{\mathrm{L}}$ is a sum of operators which act non-trivially on a local subset of qubits. An example would be:
O_L=1-1/n ∑_j=1^n|0⟩⟨0|_j ⊗1_j̅.,

where $\bar{j}$ is the complement of the $j$-th qubit. In [7] the evaluation of the cost function in (<ref>) on a quantum computer was discussed. In particular, it was shown that for a target circuit $U$ and approximate circuit $V(\theta)$ acting on $n$ qubits, one can calculate the quantity $\frac{1}{d^2}|\Tr(V^{\dagger}(\theta)U)|^2$ in equation (<ref>) by running the circuit in Figure <ref> - which acts on $2n$ qubits - and measuring the probability of observing all zeros at the end of the circuit. Since this circuit involves the measurement of all qubits the cost function is classified as a global cost function and hence suffers from the barren plateau problems discussed above.

In [7], the authors describe a natural way to create a local version of the Hilbert-Schmidt cost function. One starts with the circuit in Figure <ref> and replaces the measurements on all qubits with a measurement of qubits $j$ and $j+n$, as shown in Figure <ref>, and measures the probability $p_0^j$ of observing zeros on both qubits. The local cost function is then obtained by taking the average of all the $p_0^j$'s. In what follows, we will explicitly calculate an expression for this quantity such that it has an intuitive operational meaning.

As discussed in the introduction, another case for which a classical implementation of approximate quantum compiling is useful is when the target circuit appears at the beginning of the circuit. In this case, the problem essentially becomes one of state preparation where the target state is given by $\ket{\psi_0} = U\ket{0}$. One can use classical techniques to store the wavefunction $\ket{\psi_0}$ of the system even for large numbers of qubits if the entanglement is low enough [8, 9]. For some quantum circuits, such as those that simulate the time evolution of quantum spin chains after a quantum quench [36], the initial state is weakly entangled and grows linearly with the depth of the circuit. For such circuits, our classical implementation of approximate quantum compiling can be applied to a sub-circuit at the beginning of the algorithm when the entanglement is low. The Hilbert-Schmidt test in (<ref>) applied to state preparation becomes:

C_hs^state = 1 -  |⟨0|V^†(θ)|ψ_0⟩|^2 .

The cost function in (<ref>) is now equal to the probability of observing of all zeros after preparing the state $\ket{\psi_0}$ on the quantum computer followed by the application of $V^{\dagger}$. The local version of (<ref>) becomes:

C_local = 1-1/n∑_j=1^np_0^j

where $p_0^j$ is the probability of observing zero on the $j$-th qubit. Here we explicitly compute the expression in (<ref>) to obtain a quantity with a clear operational meaning. The final expression can be written as a sum of terms with “bit-flips”:

C_lhs^state = 1 -  |⟨0|V^†(θ)|ψ_0⟩|^2 -    (n-1/n)∑_j=1^n| ⟨0|X_jV^†(θ)|ψ_0⟩ | ^2-
(n-2/n)∑_j<k|⟨0|X_j X_kV^†(θ)|ψ_0⟩ |^2

- ... - 1/n∑_j<k<l<...|⟨0|X_j X_k X_l...V^†(θ)|ψ_0⟩|^2

Global cost function

[wires = 4]A
$\ket{0}$     H     4              [wires=4, nwires=3]U                 4     H

$\ket{0}$     H         4                      4         H

$\ket{0}$     H             4            4             H

[wires=4]B         $\ket{0}$                       [wires=4, nwires=3]V^*

The quantum circuit discussed in [7] that is used to evaluate the expression in (<ref>).

[wires = 4]A         $\ket{0}$     H     4              [wires=4, nwires=3]U                 4     H

$\ket{0}$     H         4

$\ket{0}$     H             4

[wires=4]B         $\ket{0}$                       [wires=4, nwires=3]V^*

The quantum circuit discussed in [7] that is used to evaluate the local version of the Hilbert-Schmidt test.

Local cost function

The explicit expression for $C_{lhs}^{\text{state}}$ in (<ref>) allows us to use a local cost function on a classical computer without running the circuit several times and taking shots. It was shown in [12] that if the circuit Ansatz $V(\theta)$ is sufficiently shallow then the variance of the gradient of the local cost function vanishes polynomially in the number of qubits, instead of exponentially. One can therefore expect that using the cost function in (<ref>) will lead to faster convergence when applied to classical compilation of quantum circuits. However, the large number of terms appearing in (<ref>) increases the difficulty of calculating its gradient, as compared with the expression in (<ref>) - note that this is not the case for quantum compilation schemes where one runs the circuit in Figure <ref> instead of calculating each term in (<ref>) explicitly. In what follows, we will introduce a scheme to truncate the expression in (<ref>) such that its gradient can be calculated more easily. In particular, we will take only the first $k$ terms in equation (<ref>). This opens the question as to whether the truncated expression has a barren plateau. Below, we will take a particular example such that the cost function can be studied analytically and show exactly how each "bit-flip" term in (<ref>) affects the variance of the gradient.

The specific example that we use to demonstrate the effect of the “bit-flip" terms was also studied in [12]. We take

V(θ) = ⊗_j=1^ne^-iθ_j σ_x/2

Equations (<ref>) and (<ref>) become respectively:

C_hs^state = 1 - ∏_j=1^ncos^2θ_j/2

C_lhs^state = 1 - 1/n∑_j=1^ncos^2θ_j/2

As discussed in [12], the variance of the expression in (<ref>) is given by:

Var[∂C_hs^state/∂θ_j ] = 1/8(3/8)^n-1

whereas for the expression in (<ref>) we get:

Var[∂C_lhs^state/∂θ_j] = 1/8n^2

Let's define a new truncated cost function as follows:

C_L^(k)(α_1,...,α_k) =    1 - |⟨0|V^†(θ)|ψ_0⟩|^2 - α_1∑_j=1^n|⟨0|X_jV^†(θ)|ψ_0⟩|^2

The goal is now to tune $\alpha_1,...,\alpha_k$ to increase the variance of this truncated expression. Note that if we take $k = n-1$ and $\alpha_1,...,\alpha_k = \frac{n-1}{n},...,\frac{1}{n}$ we recover the original expression in (<ref>). We can understand the intuitive meaning behind (<ref>) in terms of "bit-flips"; first note that the second term in (<ref>) is equal to exactly $1$ at the optimal values of $\theta$ and all subsequent terms are zero. We can therefore understand this expression in the following way: the second term encourages the parameters $\theta$ to tend towards those that correspond to the target state $\ket{\psi_0}$, while all other terms bias $\theta$ towards states with increasing numbers of "bit-flip" errors, i.e. $X_j\ket{\psi_0}$ and $X_kX_j\ket{\psi_0}$, which act so as to increase the size of the target in regions of barren parameter space.

It is enlightening to write the local cost function explicitly in terms of bit-flips as in equation (<ref>), in which case one obtains the expression:

C_L = 1 -    ∏_j=1^ncos^2θ_j/2 - n-1/n∑_i_1sin^2θ_i_1∏_j ≠i_1cos^2θ_j - n-2/n∑_i_1 < i_2sin^2θ_i_1sin^2θ_i_2∏_j ≠i_1, i_2cos^2θ_j

-...-1/n∑_i_1cos^2θ_i_1∏_j ≠i_1sin^2θ_j

We can rearrange (<ref>) into an expression of the form:

C_L = 1 - a_n∏_j=1^ncos^2θ_j - a_n-1∑_i_1∏_j≠i_1^ncos^2 θ_j - a_n-2∑_i_1<i_2∏_j≠i_1, i_2^ncos^2 θ_j  -...+1/n∑_i_1 cos^2 θ_j

We recover equation (<ref>) by observing that all the $a_k$ vanish. We have, for example:

a_n = 1 - (n-1/n) n 1 + (n-2/n) n 2-...+(1/n) n n-1 = 0

and similarly for the other coefficients $a_{n-1},...,a_1$.

Let's now take $k=1$, in which case equation (<ref>) becomes

C_L^(1)(α_1) ≡1 - ∏_j=1^ncos^2θ_j/2 - α_1∑_i_1sin^2θ_i_1∏_j ≠i_1cos^2θ_j

By setting $\alpha_1 = \frac{1}{n}$ we get:
C_L^(1)(α= 1/n) = 1 - 1/n∑_i∏_j≠i cos^2 θ_j

Similarly, we can define:
C_L^(2)(α_1, α_2) ≡1 - ∏_j=1^ncos^2θ_j/2 - α_1∑_i_1sin^2θ_i_1∏_j ≠i_1cos^2θ_j - α_2∑_i_1<i_2sin^2θ_i_1sin^2θ_i_2∏_j≠i_1, i_2cos^2θ_j

Now by taking $\alpha_1$ and $\alpha_2$ to be the solution to the set of linear equations:

n α_1 -α_2 n 2    = 1

-α_1 + α_2(n-1)    = 0

we get
C_L^(2) = 1 - 2/n(n-1)∑_i_1, i_2∏_j≠i_1, i_2 cos^2 θ_j

Now, we can see that each truncated local expression reduces the power of the variance in (<ref>) by 1. In particular, we get:

Var[∂C_L^(1)/∂θ_j ] ∝(3/8)^n-2

Var[∂C_L^(2)/∂θ_j ] ∝(3/8)^n-3

In the general case $C_L^{(k)}$, we have $k$ independent variables $\alpha_1,...,\alpha_k$ and we will obtain $k$ linear equations similar to (<ref>) such that we get:

Var[∂C_L^(k)/∂θ_j ] ∝(3/8)^n-k-1

We have thus shown that there exist coefficients $\alpha_1,...,\alpha_k$ such that the exponent that appears in the variance of the gradient of the truncated cost function in (<ref>) is increased by $k$ with respect to the global cost function. This has practical implications as we will demonstrate in section <ref>, where it will be shown that the convergence of (<ref>) with $k=1$ is signicantly better than (<ref>). We can understand this behaviour from the analytical results in equations (<ref>) to (<ref>). We conjecture that such results will hold for more general Ansatz and target states - see section <ref>. In section <ref> we will present our method to update the $\alpha_1,...,\alpha_k$ over the course of the optimisation.

We note that we need to extend this analysis further in order to apply it to full approximate quantum compiling, i.e. where the target is a sub-quantum circuit that does not necessarily appear at the beginning of the algorithm, in which case the target to approximate is a matrix instead of a vector. In this case, a natural way to extend the cost function in (<ref>) is:

C_lhs^circuit = 1 - 1/d^2 [     |(V^†U)|^2 + (n-1/n)∑_j=1^n |(X_j V^†U)|^2 + (n-2/n)∑_j<k|(X_j X_k V^†U)|^2

+ ... + 1/n∑_j<k<l<...|(X_j X_k X_l... V^†U)|^2 ]

In section <ref>, we will present results showing how (<ref>) improves the convergence as compared to the Hilbert-Schmidt test in equation (<ref>).

Escape from barren plateau

§ ALGORITHM

We now discuss the scheme we have developed to update the coefficients $\alpha$ appearing in equation (<ref>) during the optimisation. Consider the cost function $C_L^{(1)}(\alpha)$:

C_L^(1)(α) = 1 - |⟨0|V^†(θ)|ψ_0⟩|^2 -
α∑_j=1^n| ⟨0|X_j V^†(θ)|ψ_0⟩|^2

The parameter $\alpha$ is a function of the number of iterations, i.e. $\alpha \equiv \alpha(k)$. The first, most straightforward scheme that we implement is to gradually decrease $\alpha$ with increasing $k$ and hence reduce the influence of the final term in (<ref>) as the cost function approaches zero. In particular, we take $\alpha(k) = \sqrt{C_L^{(1)}(k)}$ and $\alpha(0) = 1$. Note, it can be shown that $C_L^{(1)}(\alpha)$ is always non-negative, even when $\alpha=1$, because $\ket{0}$ and all “bit-flipped" states constitute an orthogonal basis. Intuitively, the weighting scheme acts such as to tighten a grip around the state $\ket{0}$ as $C_L^{(1)}(\alpha)$ approaches zero.  Eventually, when $C_L^{(1)}(\alpha) \approx 0$, all the flip terms are switched off and we arrive at the desired solution $V(\theta)\ket{0} \approx \ket{\psi_0}$. In <ref> we discuss this weighting scheme in more detail, as well as our implementation of a fast-gradient approach based on [37] which allows us to do experiments with up to 24 qubits. We also developed a number of techniques to improve performance, in particular we use a “surrogate model" which speeds up the calculation of the gradient - we discuss these aspects of the algorithm in <ref>.

§ RESULTS

§.§ State Preparation

Influence of the flip terms on convergence rate. Not only does the local cost function $C_L^{(1)}$ converge at a much faster rate than $C_G$, but so too does $C_G$ when evaluated at the parameters $\theta$ that are found by optimising for $C_L^{(1)}$ - see the green curve “$C_G$ via $C_L^{(1)}$”.

We now consider the results of our experiments on approximate state preparation using the localised cost functions discussed in previous sections. We take the circuit Ansatz to be that discussed in section <ref>. In the notation defined in Figure <ref> and the surrounding text, we consider the case with $n=24$ qubits, $b=1$ repetitions and $l=1$ layers. Note that in this notation, Figure <ref> corresponds to $n=4$, $b=3$ and $l=2$. In our case, there are therefore $24$ qubits and $24$ connecting every couple of adjacent qubits. The target state $\ket{\psi_{0}}$ was generated from a circuit with the same $CNOT$ structure but parameterized by random angles. Figure <ref> clearly demonstrates that the objective function without the “bit-flip" terms $C_G$ performs much worse than $C_L^{(1)}$ which takes into account the bit-flip terms.

It was observed that second-order, quasi-Newton  method terminates prematurely if it does not detect any progress after a small number of iterations. In order to circumvent this problem, we employ two-stage optimization. Firstly, the first-order gradient-based  optimizer reduces the objective function from $1$ to a value below $0.9$, then  does the rest of the job.

§.§.§ Trotterized ansatz

R_y(0)     R_x(0)
1     R_y(0)     R_z(θ_1)
R_y(0)     R_x(0)     R_z(π/2)

-1     R_y(θ_2)              R_z(0)
R_y(θ_3)              R_x(0)
-1     R_y(0)     R_z(0)

Trotterized ansatz and its initial angular values.

In the universal ansatz, presented on Fig. <ref>, are organized in layers, where each layer contains $b \cdot (n-1)$ blocks and every block connects a pair of adjacent qubits except the first and last ones, where $b$ is the block multiplicity, $1 \le b \le 3$, and $n$ is the number of qubits. In the example depicted on Fig. <ref>, $b = 3$. Usually in our experiments, have the same “orientation”, that is, index of control qubits is smaller than index of target ones. While such an ansatz has a quite good expressibility, when we initialize parametric circuit with random initial angles (initial guess), the optimization process often converges to a local minimum with poor fidelity value. This happens because the problem in hand is highly non-linear and high-dimensional.

However, for the concrete class of problems, a better initial guess, which is close to a global minimum, could be obtained. Since we are considering time evolution approximated by Trotterization, it makes sense to change slightly ansatz layout so that it resembles a Trotter circuit.

§.§ Circuit preparation

[scale = 1.5]
[gray, thick] (0,1.13) – (0,0.13);
[black] (0,0.13) circle (2pt);
[black] (0,1) circle (4pt);
[gray, thick] (-0.14, 1) – (0.14,1);

[gray, thick] (-1.5,1) – (1.1,1);
[gray, thick] (-0.7,0.13) – (1.1,0.13);

at (-1.1,0.13) $R_z(-\frac{\pi}{2})$;

[gray, thick] (-1.5,-0.17) – (-1.5, 0.43);
[gray, thick] (-0.7,-0.17) – (-0.7, 0.43);
[gray, thick] (-1.5, 0.43) – (-0.7, 0.43);
[gray, thick] (-1.5,-0.17) – (-0.7,-0.17);

[gray, thick] (1.1,0.7) – (1.1,1.3);
[gray, thick] (1.9,0.7) – (1.9,1.3);
[gray, thick] (1.1,1.3) – (1.9,1.3);
[gray, thick] (1.1,0.7) – (1.9,0.7);

[gray, thick] (1.1,-0.17) – (1.1, 0.43);
[gray, thick] (1.9,-0.17) – (1.9, 0.43);
[gray, thick] (1.1, 0.43) – (1.9, 0.43);
[gray, thick] (1.1,-0.17) – (1.9,-0.17);

[gray, thick] (1.9,1) – (6.5,1);
[gray, thick] (1.9,0.13) – (4,0.13);

at (1.5,1) $R_y(\theta_1)$;
at (1.5,0.13) $R_y(\theta_3)$;

[gray, thick] (3,1) – (3,0);
[black] (3,1) circle (2pt);
[black] (3,0.13) circle (4pt);

at (4.4,0.13) $R_y(\theta_3)$;

[gray, thick] (4,-0.17) – (4, 0.43);
[gray, thick] (4.8,-0.17) – (4.8, 0.43);
[gray, thick] (4, 0.43) – (4.8, 0.43);
[gray, thick] (4,-0.17) – (4.8,-0.17);

[gray, thick] (4.8,0.13) – (7.8,0.13);

[gray, thick] (5.8,1.13) – (5.8,0.13);
[black] (5.8,0.13) circle (2pt);
[black] (5.8,1) circle (4pt);
[gray, thick] (5.66, 1) – (5.94,1);

[gray, thick] (6.5,0.7) – (6.5,1.3);
[gray, thick] (7.3,0.7) – (7.3,1.3);
[gray, thick] (6.5,1.3) – (7.3,1.3);
[gray, thick] (6.5,0.7) – (7.3,0.7);

[gray, thick] (7.3,1) – (7.8,1);

at (6.9,1) $R_z(\frac{\pi}{2})$;

Implementation of two site operator $e^{i(\alpha \sigma^x \otimes \sigma^x + \beta \sigma^y \otimes \sigma^y + \gamma \sigma^z \otimes \sigma^z)}$ as a quantum circuit. We have the correspondences $\theta = \frac{\pi}{2}-2\gamma$, $\phi = 2\alpha-\frac{\pi}{2}$ and $\lambda=\frac{\pi}{2}-2\beta$. The Hamiltonian in (<ref>) corresponds to the case $\alpha = \beta = \gamma = dt$

We now consider the results for full approximate quantum compiling for which the target is a matrix instead of a vector. As discussed in section <ref>, the truncated local cost function expression permits a natural extension from (<ref>) to (<ref>). In [10], full approximate compiling was applied to random unitaries. Here we consider the approximate compilation of circuits used to simulate time evolution of quantum systems. The motivation for considering these circuits is that they are approximations of the exact time evolution operator - to be discussed in detail below - and thus provide a useful benchmark against which we can compare the accuracy of our compilation. Time evolution of quantum systems is governed by the Schrödinger equation:

|ψ(t)⟩= e^-iHt |ψ(0)⟩

where $\ket{\psi(0)}$ is the wavefunction at time $t=0$, and $H$ is the Hamiltonian of the system under study. We will consider the XXX Hamiltonian - a paradigmatic model for quantum magnetism - written explicitly in terms of Pauli matrices $\sigma^x, \sigma^y$ and $\sigma^z$ as:

\begin{equation}\label{hxxz}
H_{XXZ} = -\sum_{\ell=0}^{L-1}\left( \sigma^x_\ell \sigma^x_{\ell+1} + \sigma^y_\ell \sigma^y_{\ell+1} + \sigma^z_\ell \sigma^z_{\ell+1}  \right),
\end{equation}

The time evolution operator $U(t)\equiv e^{-iHt}$ can be executed as a quantum circuit via Trotterisation which we now explain. We define $h_{\ell, \ell+1} = \left(\sigma^x_\ell \sigma^x_{\ell+1} + \sigma^y_\ell \sigma^y_{\ell+1} + \sigma^z_\ell \sigma^z_{\ell+1}\right)$ such that the Hamiltonian in (<ref>) becomes $H_{XXZ} = H_1 + H_2$ where $H_1 = -\sum\limits_{\ell\text{ odd}} h_{\ell, \ell+1}$ and $H_2 = \sum\limits_{\ell\text{ even}} h_{\ell, \ell+1}$. Note that all operators in a given sum commute with all other operators in their respective sums. We then define the Trotterised time evolution operator $\mathcal{U}_{\text{trot}}(dt)$ in the following way:

𝒰_trot(dt) = ∏_j=0^L/2-1U_2j, 2j+1(dt)∏_j=1^L/2-1U_2j-1, 2j(dt) = e^-iH_XXZdt + O(dt^2)
where $U_{j k}(dt) = e^{-ih_{j k}dt}$. The exact time evolution operator $U(t)$ is thus approximated by $m$ repeated applications of $\mathcal{U}_{\text{trot}}(dt=\frac{t}{m})$, i.e. $U(t) \approx \mathcal{U}^m_{\text{trot}}(dt=\frac{t}{m})$. As discussed in [36], each $U_{j k}(dt)$ appearing in (<ref>) can be implemented by the quantum circuit in Figure <ref>. Equation (<ref>) is a first order expression, one can also consider second and higher order Trotter operators [38] but with increased circuit depth. Note that $\mathcal{U}_{\text{trot}}(dt)$ is exact if one considers a time dependent Hamiltonian in the Floquet formalism, thus allowing for highly accurate quantum simulations [39].

We now compare the fidelity of the operator we find via approximate quantum compiling with the exact time evolution operator vs the fidelity of the Trotter operator vs the time evolution operator. We note that a number of other techniques have been developed recently that aim to find short depth quantum circuits that implement time evolution [40, 41, 42, 43]. The fidelity of two matrix operators $U$ and $V$ of dimension $d \times d$ is defined as:

f(U,V) = 1 + |V^† U|^2/d/d+1

Approximate compilation of time evolution operator on 8 qubits - the horizontal dotted line corresponds to the fidelity that one obtains with the standard first order Trotter decomposition with two Trotter steps and with $dt=0.2$.

In Figure <ref> we plot the fidelity of the approximate quantum compiling operator with the exact time evolution operator vs the number of iterations in the optimisation procedure. We can see that both the global cost function $C_G$ defined in (<ref>) and the partially localised cost function defined in (<ref>) with $k=1$ reach a barren plateau and converge to an operator that has poor fidelity with the target. However, we can see that $C_L^{(2)}$ performs much better and not only does it escape from the barren plateau, but achieves a fidelity that is higher than that of the Trotter circuit.

§ DISCUSSION

In this work we have developed a set of techniques that significantly improve the convergence of classically assisted quantum compilation algorithms. Just like the local cost functions considered in [7, 12], we have demonstrated in practice that our truncated local cost functions and their gradients (see e.g. equations (<ref>) and (<ref>)) can be feasibly calculated on a classical computer while also being powerful enough to escape barren plateaus that appear in both approximate state preparation and full circuit approximation - see Figures <ref> and <ref>. As discussed in the introduction, there are a number of circumstances where this classically assisted approach can be useful. One might for example be interested in compiling sub-circuits that appear in a given quantum algorithm, or alternatively one might be interested in finding the most efficient way to generate a weakly entangled state on a quantum computer [9]; one can then use this state generated by a shallow circuit as the input for a quantum algorithm [8]. In the latter case, one can use Tensor Networks to represent the state classically and thus one can apply the techniques developed here on large numbers of qubits - see our upcoming work.

[1]
Giacomo Nannicini, Lev S Bishop, Oktay Günlük, and Petar Jurcevic.
Optimal qubit assignment and routing via integer programming.
ACM Transactions on Quantum Computing, 2021.

[2]
Kristan Temme, Sergey Bravyi, and Jay M Gambetta.
Error mitigation for short-depth quantum circuits.
Physical review letters, 119(18):180509, 2017.

[3]
Youngseok Kim, Christopher J Wood, Theodore J Yoder, Seth T Merkel, Jay M
Gambetta, Kristan Temme, and Abhinav Kandala.
Scalable error mitigation for noisy quantum circuits produces
competitive expectation values.
arXiv preprint arXiv:2108.09197, 2021.

[4]
Tudor Giurgica-Tiron, Yousef Hindy, Ryan LaRose, Andrea Mari, and William J
Digital zero noise extrapolation for quantum error mitigation.
In 2020 IEEE International Conference on Quantum Computing and
Engineering (QCE), pages 306–316. IEEE, 2020.

[5]
Kunal Sharma, Sumeet Khatri, Marco Cerezo, and Patrick J Coles.
Noise resilience of variational quantum compiling.
New Journal of Physics, 22(4):043006, 2020.

[6]
Samson Wang, Enrico Fontana, Marco Cerezo, Kunal Sharma, Akira Sone, Lukasz
Cincio, and Patrick J Coles.
Noise-induced barren plateaus in variational quantum algorithms.
Nature communications, 12(1):1–11, 2021.

[7]
Sumeet Khatri, Ryan LaRose, Alexander Poremba, Lukasz Cincio, Andrew T
Sornborger, and Patrick J Coles.
Quantum-assisted quantum compiling.
Quantum, 3:140, 2019.

[8]
Manuel S Rudolph, Jacob Miller, Jing Chen, Atithi Acharya, and Alejandro
Synergy between quantum circuits and tensor networks: Short-cutting
the race to practical quantum advantage.
arXiv preprint arXiv:2208.13673, 2022.

[9]
Manuel S Rudolph, Jing Chen, Jacob Miller, Atithi Acharya, and Alejandro
Decomposition of matrix product states into shallow quantum circuits.
arXiv preprint arXiv:2209.00595, 2022.

[10]
Liam Madden and Andrea Simonetto.
Best approximate quantum compiling problems.
ACM Transactions on Quantum Computing, 3(2):1–29, 2022.

[11]
Vivek V. Shende, Stephen S. Bullock, and Igor Markov.
Synthesis of quantum logic circuits.
IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst., 25:1000,

[12]
Marco Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J Coles.
Cost function dependent barren plateaus in shallow parametrized
quantum circuits.
Nature communications, 12(1):1–12, 2021.

[13]
Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca.
Fast and efficient exact synthesis of single qubit unitaries
generated by clifford and t gates.
arXiv preprint arXiv:1206.5236, 2012.

[14]
Neil J Ross and Peter Selinger.
Optimal ancilla-free Clifford+ T approximation of z-rotations.
Quantum Information & Computation, 16(11-12):901–953, 2016.

[15]
Peter Sarnak.
Letter to scott aaronson and andy pollington on the solovay-kitaev
theorem and golden gates.
Google Scholar Google Scholar Reference, 2015.

[16]
Peter Selinger.
Generators and relations for n-qubit Clifford operators.
Logical Methods in Computer Science, 11, 2015.

[17]
Jun Zhang, Jiri Vala, Shankar Sastry, and K. Birgitta Whaley.
Minimum construction of two-qubit quantum operations.
Phys. Rev. Lett., 93:020502, Jul 2004.

[18]
Dmitri Maslov, Gerhard W Dueck, D Michael Miller, and Camille Negrevergne.
Quantum circuit simplification and level compaction.
IEEE Transactions on Computer-Aided Design of Integrated
Circuits and Systems, 27(3):436–444, 2008.

[19]
Matthew Amy, Dmitri Maslov, Michele Mosca, and Martin Roetteler.
A meet-in-the-middle algorithm for fast synthesis of depth-optimal
quantum circuits.
IEEE Transactions on Computer-Aided Design of Integrated
Circuits and Systems, 32(6):818–830, 2013.

[20]
Yunseong Nam, Neil J. Ross, Yuan Su, Andrew M. Childs, and Dmitri Maslov.
Automated optimization of large quantum circuits with continuous
npj Quantum Information, 4(1):23, 2018.

[21]
Dmitri Maslov, Sean M Falconer, and Michele Mosca.
Quantum circuit placement.
IEEE Transactions on Computer-Aided Design of Integrated
Circuits and Systems, 27(4):752–763, 2008.

[22]
Debjyoti Bhattacharjee and Anupam Chattopadhyay.
Depth-optimal quantum circuit placement for arbitrary topologies.
arXiv preprint arXiv:1703.08540, 2017.

[23]
Angelo Oddi and Riccardo Rasconi.
Greedy randomized search for scalable compilation of quantum
In Integration of Constraint Programming, Artificial
Intelligence, and Operations Research, pages 446–461, Cham, 2018. Springer
International Publishing.

[24]
Kyle E. C. Booth, Minh Do, J Christopher Beck, Eleanor Rieffel, Davide
Venturelli, and Jeremy Frank.
Comparing and integrating constraint programming and temporal
planning for quantum circuit compilation.
In Proceedings of the International Conference on Automated
Planning and Scheduling (ICAPS2018), pages 366–374, 2018.

[25]
Alwin Zulehner and Robert Wille.
Compiling SU(4) Quantum Circuits to IBM QX Architectures.
In Proceedings of the 24th Asia and South Pacific Design
Automation Conference, page 185–190, 2019.

[26]
D. Bhattacharjee, A. A. Saki, M. Alam, A. Chattopadhyay, and
S. Ghosh.
MUQUT: Multi-Constraint Quantum Circuit Mapping on NISQ Computers:
Invited Paper.
In 2019 IEEE/ACM International Conference on Computer-Aided
Design (ICCAD), pages 1–7, 2019.

[27]
A. Zulehner, A. Paler, and R. Wille.
An Efficient Methodology for Mapping Quantum Circuits to the IBM QX
IEEE Transactions on Computer-Aided Design of Integrated
Circuits and Systems, 38(7):1226–1236, 2019.

[28]
Alexander Cowtan, Silas Dilkes, Ross Duncan, Alexandre Krajenbrink, Will
Simmons, and Seyon Sivarajah.
On the Qubit Routing Problem.
In 14th Conference on the Theory of Quantum Computation,
Communication and Cryptography (TQC 2019), volume 135, pages 5:1–5:32,

[29]
Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington,
and Ross Duncan.
t$\vert$ket${\rangle}$: a retargetable compiler for NISQ devices.
Quantum Science and Technology, 6(1):014003, nov 2020.

[30]
Toshinari Itoko, Rudy Raymond, Takashi Imamichi, and Atsushi Matsuo.
Optimization of quantum circuit mapping using gate transformation and
Integration, 70:43–50, 2020.

[31]
Bochen Tan and Jason Cong.
Optimality study of existing quantum computing layout synthesis
IEEE Transactions on Computers, 70(9):1363–1373, 2020.

[32]
Emanuel Knill.
Approximation by quantum circuits.
arXiv preprint quant-ph/9508006, 1995.

[33]
Adi Botea, Akihiro Kishimoto, and Radu Marinescu.
On the complexity of quantum circuit compilation.
In Eleventh annual symposium on combinatorial search, 2018.

[34]
Peter W Shor.
Polynomial-time algorithms for prime factorization and discrete
logarithms on a quantum computer.
SIAM review, 41(2):303–332, 1999.

[35]
Zoë Holmes, Kunal Sharma, Marco Cerezo, and Patrick J Coles.
Connecting ansatz expressibility to gradient magnitudes and barren
PRX Quantum, 3(1):010313, 2022.

[36]
Adam Smith, MS Kim, Frank Pollmann, and Johannes Knolle.
Simulating quantum many-body dynamics on a current digital quantum
npj Quantum Information, 5(1):1–13, 2019.

[37]
Liam Madden, Albert Akhriev, and Andrea Simonetto.
Sketching the best approximate quantum compiling problem.
arXiv preprint arXiv:2205.04025, 2022.

[38]
David Layden.
First-order trotter error from a second-order perspective.
Physical Review Letters, 128(21):210501, 2022.

[39]
Nathan Keenan, Niall Robertson, Tara Murphy, Sergiy Zhuk, and John Goold.
Evidence of kardar-parisi-zhang scaling on a digital quantum
arXiv preprint arXiv:2208.12243, 2022.

[40]
Joe Gibbs, Kaitlin Gili, Zoë Holmes, Benjamin Commeau, Andrew Arrasmith,
Lukasz Cincio, Patrick J Coles, and Andrew Sornborger.
Long-time simulations with high fidelity on quantum hardware.
arXiv preprint arXiv:2102.04313, 2021.

[41]
Cristina Cirstoiu, Zoe Holmes, Joseph Iosue, Lukasz Cincio, Patrick J Coles,
and Andrew Sornborger.
Variational fast forwarding for quantum simulation beyond the
coherence time.
npj Quantum Information, 6(1):1–10, 2020.

[42]
Christa Zoufal, David Sutter, and Stefan Woerner.
Error bounds for variational quantum time evolution.
arXiv preprint arXiv:2108.00022, 2021.

[43]
Kishor Bharti and Tobias Haug.
Quantum-assisted simulator.
Physical Review A, 104(4):042418, 2021.

§ WEIGHTING SCHEME

As discussed in section <ref>, the parameter $\alpha$ in equation (<ref>) is gradually decreased over the optimisation, enforcing $V^{\dagger}(\theta)\ket{\psi_0}$ to have a large projection on $\ket{0}$ while progressing mostly inside the sub-space spanned by states $\ket{0}$ and $\{X_j \ket{0} | 1 \le j \le n\}$. One such approach - as mentioned in the main text - is to take $\alpha(k) = \sqrt{C_L^{(1)}(k)}$ and $\alpha(0) = 1$. We find that slightly modifying the expression for $\alpha(k)$ to e.g.:

\begin{equation}\label{eq:weight}
w^{k} = 0.1 \cdot \sqrt{C_{lhs}^{k-1}} + 0.9 \cdot w^{k-1}, \quad w^0 = 1.
\end{equation}

can enhance the optimisation procedure due to the reduced noise in the gradient profile. To increase speed, we apply a fast gradient approach based on [37] by using a more economical matrix-vector multiplication. However this approach must be modified in order to apply it to our case; the main bottleneck is that the gradient of every term in the sum in equation (<ref>) should be computed separately and then all of them summed up. This causes issues with speed when considering sizes of up to 24 qubits. To handle this, we designed and implemented a surrogate model approach which we describe in <ref>.

§ SURROGATE MODEL

In the second formulation, instead of summing up squared projections onto individual flip-states, we minimize the maximal projection of $V^{\dagger}(\theta)\ket{\psi_0}$ onto the sub-space spanned by states $\ket{0}$ and $\{X_j \ket{0} | 1 \le j \le n\}$, which we call a surrogate model:

\begin{equation}\label{eq:fobj-sur}
\begin{aligned}
C_{lhs}^{surr}(\theta; \alpha) &= \min_{{\boldsymbol\beta}, \|{\boldsymbol\beta}\|^2=1} \Bigg\{
1 - (1 - \alpha) |\bra{0}V^{\dagger}(\theta)\ket{\psi_0}|^2 - \\
& \qquad\qquad\qquad - \alpha\left|\left(\beta_0^* \bra{0} + \sum\limits_{j=1}^n \beta_j^* \bra{0}X_j\right)
\Bigg\},
\end{aligned}
\end{equation}

where on every iteration the objective is minimized over the parameters ${\boldsymbol\beta}$ by taking the maximal projection onto the sub-space of states $\ket{0}$ and $\{X_j \ket{0} | 1 \le j \le n\}$ (which is the same as maximizing the third term) before optimization over the angular parameters $\theta$. It is straightforward to show that the following normalized $\beta$'s, $\|{\boldsymbol\beta}\|^2 = 1$, maximize the projection: $\beta_0 \sim \bra{0}V^{\dagger}(\theta)\ket{\psi_0}$ and $\beta_j \sim \bra{0}X_jV^{\dagger}(\theta)\ket{\psi_0}$[Really, let us define $y_0$ = $\bra{0}V^{\dagger}(\theta)\ket{\psi_0}$ and $y_j$ = $\bra{0}X_jV^{\dagger}(\theta)\ket{\psi_0}$. Then the last term in (<ref>) takes the form of squared module of a dot product: $|\sum\nolimits_{i=0}^n \beta_i^* y_i|^2$, which is maximized when $\beta_i \sim y_i$.].

On every iteration, we compute the latter dot products for $\beta$'s and normalize the vector ${\boldsymbol\beta}$, then substitute $\beta$-parameters in the last term in (<ref>) and compute the composite state $\beta_0^* \bra{0} + \sum\nolimits_{j=1}^n \beta_j^* \bra{0}X_j$. Afterwards, we proceed with optimization over the parameters $\theta$ as usually. The computational gain is obvious — the gradient is computed for just two dot products instead of $n+1$ ones, and the gradient calculation is known to be the most time-consuming operation.

For the additional speed up we compute the state $V^{\dagger}(\theta)\ket{\psi_0}$ at the beginning of each optimization iteration given the current parameters $\theta$ and use it for computation of both, the objective function and its gradient, see [37] for details. This substantially reduces the processing time.

Despite the advantage of surrogate model, adding up individual states to get a composite one, $\beta_0^* \bra{0} + \sum\nolimits_{j=1}^n \beta_j^* \bra{0}X_j$, can be challenging if we want to generalize the framework from using an ordinary vectors to matrix-product state (MPS) format for scalability. In view of future generalization, we devised even more simple surrogate model, where instead of composite state we pick up one of flipped states or $\ket{0}$, whichever gives the maximal projection onto $V^{\dagger}(\theta)\ket{\psi_0}$, and solely use it in the last term:

\begin{equation}\label{eq:fobj-sur-max}
C_{lhs}^{max}(\theta; \alpha) =
1 - (1 - \alpha) |\bra{0}V^{\dagger}(\theta)\ket{\psi_0}|^2 -
\alpha\left|\bra{j_{max}}V^{\dagger}(\theta)\ket{\psi_0}\right|^2,
\end{equation}

where $j_{max} = 0$, if $|\bra{0}V^{\dagger}\ket{\psi_0}|$ gives the largest projection, or $j_{max} = j^{\prime}$, $1 \le j^{\prime} \le n$, if $|\bra{0}X_{j^{\prime}}V^{\dagger}\ket{\psi_0}|$ has the greatest magnitude.

One more feature is required to make the objective (<ref>) viable. If component of the largest projection is altered too often, this could contribute a strong disturbance into the optimization process. In order to alleviate that adverse effect, we employ the idea of hysteresis. Namely, the former leading state $\ket{j_{max}}$ only gives up if its projection onto $V^{\dagger}(\theta)\ket{\psi_0}$ is $10\%$ lower than the projection of the currently leading state, otherwise optimization proceeds with the former leader.