# Long-time simulations with high fidelity on quantum hardware

Joe Gibbs  Theoretical Division, Los Alamos National Laboratory, Los Alamos,
NM, USA.  Kaitlin Gili  Theoretical Division, Los Alamos National Laboratory,
Los Alamos, NM, USA.  Department of Physics, University of Oxford, Clarendon
Laboratory, Oxford, U.K.  Zoë Holmes  Information Sciences, Los Alamos
National Laboratory, Los Alamos, NM, USA.  Benjamin Commeau  Information
Sciences, Los Alamos National Laboratory, Los Alamos, NM, USA.  Department of
Physics, University of Connecticut, Storrs, Connecticut, CT, USA.  Andrew
Arrasmith  Theoretical Division, Los Alamos National Laboratory, Los Alamos,
NM, USA.  Lukasz Cincio  Theoretical Division, Los Alamos National Laboratory,
Los Alamos, NM, USA.  Patrick J. Coles  Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM, USA.  Andrew Sornborger  Information
Sciences, Los Alamos National Laboratory, Los Alamos, NM, USA.

###### Abstract

Moderate-size quantum computers are now publicly accessible over the cloud,
opening the exciting possibility of performing dynamical simulations of
quantum systems. However, while rapidly improving, these devices have short
coherence times, limiting the depth of algorithms that may be successfully
implemented. Here we demonstrate that, despite these limitations, it is
possible to implement long-time, high fidelity simulations on current
hardware. Specifically, we simulate an XY-model spin chain on the Rigetti and
IBM quantum computers, maintaining a fidelity of at least 0.9 for over 600
time steps. This is a factor of 150 longer than is possible using the iterated
Trotter method. Our simulations are performed using a new algorithm that we
call the fixed state Variational Fast Forwarding (fsVFF) algorithm. This
algorithm decreases the circuit depth and width required for a quantum
simulation by finding an approximate diagonalization of a short time evolution
unitary. Crucially, fsVFF only requires finding a diagonalization on the
subspace spanned by the initial state, rather than on the total Hilbert space
as with previous methods, substantially reducing the required resources. We
further demonstrate the viability of fsVFF through large numerical
implementations of the algorithm, as well as an analysis of its noise
resilience and the scaling of simulation errors.

##  I Introduction

The simulation of physical systems is both valuable for basic science and
technological applications across a diverse range of industries, from
materials design to pharmaceutical development. Relative to classical
computers, quantum computers have the potential to provide an exponentially
more efficient means of simulating quantum mechanical systems. Quantum
hardware has progressed substantially in recent years Arute  _et al._ (2019,
2020). However, despite continual progress, we remain in the ‘noisy
intermediate-scale quantum’ (NISQ) era in which the available hardware is
limited to relatively small numbers of qubits and prone to errors. Simulation
algorithms designed for fault-tolerant quantum computers, such as
Trotterization methods Lloyd (1996); Sornborger and Stewart (1999),
qubitization methods Low and Chuang (2019), and Taylor series methods Berry
_et al._ (2015), require deeper circuits than viable given the short coherence
times of current hardware. Thus alternative approaches are needed to
successfully implement useful simulations on NISQ hardware.

Variational quantum algorithms Cerezo  _et al._ (2020a); Endo  _et al._
(2021); Bharti  _et al._ (2021); Peruzzo  _et al._ (2014); Farhi  _et al._
(2014); McClean  _et al._ (2016); Khatri  _et al._ (2019); LaRose  _et al._
(2019); Arrasmith  _et al._ (2019); Cerezo  _et al._ (2020b); Li and Benjamin
(2017); Endo  _et al._ (2020); Yao  _et al._ (2020); Heya  _et al._ (2019);
Cirstoiu  _et al._ (2020); Commeau  _et al._ (2020), where a classical
computer optimizes a cost function measured on a quantum computer, show
promise for NISQ quantum simulations. An early approach introduced an
iterative method, where the state is variationally learned on a step-by-step
basis using action principles Trout  _et al._ (2018); Endo  _et al._ (2020);
Yao  _et al._ (2020); Benedetti  _et al._ (2020). Subsequently, a
generalization of the variational quantum eigensolver Peruzzo  _et al._ (2014)
was developed for simulations in low lying energy subspaces Heya  _et al._
(2019). Very recently, quantum-assisted methods have been proposed that
perform all necessary quantum measurements at the start of the algorithm
instead of employing a classical-quantum feedback loop Bharti and Haug (2020);
Lau  _et al._ (2021); Haug and Bharti (2020).

In this work, we improve upon a recently proposed variational quantum
algorithm known as Variational Fast Forwarding (VFF) Cirstoiu  _et al._
(2020). VFF allows long time simulations to be performed using a fixed depth
circuit, thus enabling a quantum simulation to be ‘fast forwarded’ beyond the
coherence time of noisy hardware. The VFF algorithm requires finding a full
diagonalization of the short time evolution operator $U$ of the system of
interest. Once found, the diagonalization enables any initial state of that
system to be fast forwarded. However, for practical purposes, one is often
interested in studying the evolution of a particular fixed initial state of
interest. In that case a full diagonalization of $U$ is overkill. Instead, it
suffices to find a diagonal compilation of $U$ that captures its action on the
given initial state. Here, we show that focusing on this commonly encountered
but less exacting task can substantially reduce the resources required for the
simulation.

Specifically, we introduce the fixed state VFF algorithm (fsVFF) for fast
forwarding a fixed initial state beyond the coherence time of a quantum
computer. This approach is tailored to making dynamical simulation more
suitable for NISQ hardware in two key ways. First, the cost function requires
half as many qubits as VFF. This not only allows larger scale simulations to
be performed on current resource-limited hardware, but also has the potential
to enable higher fidelity simulations since larger devices tend to be noisier.
Second, fsVFF can utilize simpler ansätze than VFF both in terms of the depth
of the ansatz and the number of parameters that need to be learnt. Thus, fsVFF
can reduce the width, depth, and total number of circuits required to fast
forward quantum simulations, hence increasing the viability of performing
simulations on near-term hardware.

We demonstrate these advantages by implementing long-time high fidelity
quantum simulations of the 2-qubit XY spin chain on Rigetti’s and IBM’s
quantum computers. Specifically, while the iterated Trotter approach has a
fidelity of less than 0.9 after 4 time steps and has completely thermalized by
25 time steps, with fsVFF we achieve a simulation fidelity greater than 0.9
for over 600 time steps. We further support the effectiveness of this approach
for NISQ simulations, with 4 qubit noisy and 8 qubit noiseless numerical
simulations of the XY model and Fermi-Hubbard model respectively.

In our analytical results, we prove the faithfulness of the fsVFF cost
function by utilizing the newly developed No-Free-Lunch theorems for quantum
machine learning Poland  _et al._ (2020); Sharma  _et al._ (2020a). We also
provide a proof of the noise resilience of the fsVFF cost function,
specifically the optimal parameter resilience Sharma  _et al._ (2020b).
Finally, we perform an analysis of simulation errors under fast-forwarding.

The diagonalizations obtained using fsVFF may further be useful for
determining the eigenstates and eigenvalues of the Hamiltonian on the subspace
spanned by the initial state. This can be done using a time series analysis,
by using fsVFF to reduce the depth of the quantum phase estimation (QPE)
algorithm, or using a simple sampling method. We demonstrate on IBM’s quantum
computer that, while standard QPE fails on real hardware, fsVFF can be used to
obtain accurate estimates of the spectrum.

##  II Background

Before presenting our fsVFF algorithm, let us first review the original VFF
algorithm from Ref. Cirstoiu  _et al._ (2020). Consider a Hamiltonian $H$ on a
$d=2^{n}$ dimensional Hilbert space (i.e., on $n$ qubits) evolved for a short
time $\Delta t$ with the simulation unitary $e^{-iH\Delta t}$, and let $T$
(larger than $\Delta t$) denote the desired simulation time. Then the VFF
algorithm consists of the following steps:

  1. 1.

Approximate $e^{-iH\Delta t}$ with a single-timestep Trotterized unitary
denoted $U=U(\Delta t)$.

  2. 2.

Variationally search for an approximate diagonalization of $U$ by compiling it
to a unitary with a structure of the form

$\displaystyle V(\bm{\alpha},\Delta t):=W(\bm{\theta})D(\bm{\gamma},\Delta
t)W(\bm{\theta})^{\dagger}\,,$ (1)

where $\bm{\alpha}=(\bm{\theta},\bm{\gamma})$ is a vector of parameters. Here,
$D(\bm{\gamma},\Delta t)$ is a parameterized unitary that will (after
training) encode the eigenvalues of $U(\Delta t)$, while $W(\bm{\theta})$ is a
parameterized unitary matrix that will consist of the corresponding
eigenvectors Cirstoiu  _et al._ (2020). The compilation is performed using the
local Hilbert-Schmidt test Khatri  _et al._ (2019) to find the parameters
$\bm{\theta}_{\rm opt}$ and $\bm{\gamma}_{\rm opt}$ that minimize the local
Hilbert-Schmidt cost.

  3. 3.

Use the compiled form to simulate for time $T=N\Delta t$ using the circuit

$\displaystyle W(\bm{\theta}_{\rm opt})D(\bm{\gamma}_{\rm opt},N\Delta
t)W(\bm{\theta}_{\rm opt})^{\dagger}\,.$ (2)

VFF has proven effective for providing a fixed quantum circuit structure with
which to fast-forward beyond the coherence time of current noisy quantum
devices. However, the algorithm requires a full diagonalization of $U$ over
the entire Hilbert space. The local Hilbert-Schmidt test used to find this
diagonalization requires $2n$ qubits. Additionally, the ansatz must be
sufficiently expressible to diagonalize the full unitary $U$ to a high degree
of approximation Sim  _et al._ (2019); Nakaji and Yamamoto (2020); Holmes  _et
al._ (2021a). This typically requires a large number of parameters and a
reasonably deep circuit. These overheads limit VFF’s utility on current
hardware.

In what follows, we introduce a more NISQ-friendly refinement to VFF that
reduces these overheads when one is interested in fast-forwarding a fixed
initial state $\ket{\psi_{0}}$, rather than any possible initial state. The
fixed state VFF algorithm (fsVFF) is summarised in Fig. 1.

We note that VFF, like the standard iterated Trotter approach to quantum
simulation, necessarily incurs a Trotter error by approximating $e^{-iH\Delta
t}$ with $U=U(\Delta t)$. This Trotter error may be removed using the
Variational Hamiltonian Diagonalization algorithm (VHD), which directly
diagonalizes the Hamiltonian $H$ Commeau  _et al._ (2020). However, VHD is yet
more resource intensive than VFF on current hardware, so we focus here on
refining VFF.

##  III Fixed State Variational Fast Forwarding Algorithm

Figure 1: The fsVFF Algorithm. (a) An input Hamiltonian and an initial input
state are necessary (b) to create a single time-step Trotterized unitary,
$U(\Delta t)$ and (c) to calculate the number of eigenstates spanned by the
initial state. The value of $n_{\rm eig}$ can be calculated by constructing a
matrix of state overlaps $U^{k}\ket{\psi_{0}}$ and increasing the matrix
dimension until the determinant is zero. (d) The unitary is then variationally
diagonalized into the form, $V(\bm{\alpha},\Delta
t)=W(\bm{\theta})D(\bm{\gamma},\Delta t)W^{\dagger}(\bm{\theta})$. The cost
function $C_{\rm fsVFF}$ is minimized with a classical optimizer (e.g.,
gradient descent), where the parameters $\bm{\theta}$ and $\bm{\gamma}$ are
updated. (e) The optimal parameters $\bm{\theta}_{\rm opt}$ and
$\bm{\gamma}_{\rm opt}$ are then used to implement a fast-forwarded simulation
with the diagonalized unitary form.

###  III.1 Cost function

In fsVFF, instead of searching for a full diagonalization of $U$ over the
entire Hilbert space, we search for a diagonal compilation of $U$ that
captures the action of $U$ on the initial state $\ket{\psi_{0}}$ and its
future evolution, $e^{-iHt}\ket{\psi_{0}}$. Here, we introduce a cost function
tailored to this task.

To make precise what is required of the cost for fsVFF, let us first note that
as the state $\ket{\psi_{0}}$ evolves, it remains within its initial energy
subspace. This can be seen by expanding the initial state in terms of the
energy eigenbasis $\\{\ket{E_{k}}\\}_{k=1}^{2^{n}}$ (the eigenbasis of $H$) as

$\ket{\psi_{0}}=\sum_{k=1}^{n_{\rm eig}}a_{k}\ket{E_{k}},$ (3)

where $a_{k}=\braket{E_{k}}{\psi_{0}}$, and noting that

$e^{-iHt}\ket{\psi_{0}}=\sum_{k=1}^{n_{\rm
eig}}a_{k}e^{-iE_{k}t}\ket{E_{k}}\,.$ (4)

Thus it follows that if $\ket{\psi_{0}}$ spans $n_{\rm eig}$ energy
eigenstates of $H$, so does $e^{-iHt}\ket{\psi_{0}}$ for all future times.
Therefore to find a compilation of $U$ that captures its action on
$e^{-iHt}\ket{\psi_{0}}$ (for all times $t$) it suffices to find a compilation
of $U$ on the $n_{\rm eig}$ dimensional subspace spanned by
$\\{\ket{E_{k}}\\}_{k=1}^{n_{\rm eig}}$. We stress that the eigenstates
$\\{\ket{E_{k}}\\}_{k=1}^{2^{n}}$ need not be ordered, and therefore the
subspace spanned by the subset $\\{\ket{E_{k}}\\}_{k=1}^{n_{\rm eig}}$ is not
necessarily low lying in energy.

A No-Free-Lunch Theorem for quantum machine learning introduced in Ref. Poland
_et al._ (2020) proves that to perfectly learn the action of a unitary on a
$d$-dimensional space requires $d$ training pairs. In the context of fsVFF, we
are interested in learning the action of a unitary on an $n_{\rm
eig}$-dimensional subspace. Since the unitary is block diagonal, one can
directly apply this NFL theorem to the subspace of interest. Therefore $n_{\rm
eig}$ training pairs are required to learn the unitary’s action on this
subspace. (Note, we assume here that the training states are not entangled
with an additional register. It was shown in Ref. Sharma  _et al._ (2020a)
that using entangled training data can reduce the required number of training
states. In fact, this more powerful method is used by the VFF algorithm.
However, producing such entangled training data requires additional qubits and
two-qubit gates and therefore is less NISQ-friendly.)

The No-Free-Lunch theorem therefore implies that $n_{\rm eig}$ states are
required to learn $U$ on $\ket{\psi_{0}}$ (assuming leakage due to Trotter
error is negligible). In general these states may be freely chosen from the
subspace spanned by $\ket{\psi_{0}}$. Here a convenient choice in training
states would be $\ket{\psi_{0}}$ and its Trotter evolutions, that is the set
$\\{U^{k}\ket{\psi_{0}}\\}_{k=1}^{n_{\rm eig}}$. Motivated by these
observations, we define our cost function for fsVFF as

$C_{\rm fsVFF}:=1-\frac{1}{n_{\rm eig}}\sum_{k=1}^{n_{\rm
eig}}|\bra{\psi_{0}}(V^{\dagger})^{k}U^{k}\ket{\psi_{0}}|^{2}\,,$ (5)

where similarly to VFF we use a diagonal ansatz $V(\bm{\alpha},\Delta
t):=W(\bm{\theta})D(\bm{\gamma},\Delta t)W(\bm{\theta})^{\dagger}$. This cost
quantifies the overlap between the initial state evolved under $U$ for $k$
time steps, $U^{k}\ket{\psi_{0}}$, and the initial state evolved under the
trained unitary for $k$ time steps, $WD^{k}W^{\dagger}\ket{\psi_{0}}$,
averaged over $n_{\rm eig}$ time steps. Assuming we have access to the unitary
that prepares the state $\ket{\psi_{0}}$, the state overlaps can be measured
using $n$ qubits, via a circuit that performs a Loschmidt echo Sharma  _et
al._ (2020b). Therefore $C_{\rm fsVFF}$ can be evaluated using only $n$
qubits. This is half as many as standard VFF, opening up the possibility of
performing larger simulations on current hardware.

It is important to note that while the exact time-evolved state
$\exp(-iHt)\ket{\psi_{0}}$ is perfectly confined to the initial subspace, the
approximate evolution induced by $U(\Delta t)$ allows for leakage from the
initial subspace Sahinoglu and Somma (2020). Thus the subspace spanned by
$\\{U^{k}\ket{\psi_{0}}\\}_{k=1}^{n_{\rm eig}}$ in general does not perfectly
overlap with $\\{\ket{E_{k}}\\}_{k=1}^{n_{\rm eig}}$. However, by reducing
$\Delta t$ and considering higher order Trotter approximations Sornborger and
Stewart (1999); Suzuki (1976), this leakage can be made arbitrarily small. In
Appendix A, we prove that in the limit that leakage from the initial subspace
is negligible, $C_{\rm fsVFF}$ is faithful. That is, we show that the cost
vanishes, $C_{\rm fsVFF}=0$, if and only if the fidelity of the fast-forwarded
simulation is perfect,

$F_{\tau}=|\bra{\psi_{0}}W^{\dagger}D^{\tau}WU^{\tau}\ket{\psi_{0}}|^{2}=1\,,$
(6)

for all times $\tau$. Note, that the reverse direction is trivial. If
$F_{\tau}=1$ for all $\tau$, then $C_{\rm fsVFF}=0$.

Similar to the VFF cost, the fsVFF cost is noise resilient in the sense that
incoherent noise should not affect the global optimum of the function. This is
proven for a broad class of incoherent noise models using the results of Ref.
Sharma  _et al._ (2020b) in Appendix B.

Nonetheless, it is only possible to measure $C_{\rm fsVFF}$ if the unitary
$U^{n_{\rm eig}}$ can be implemented comfortably within the coherence time of
the QC. Additionally, the number of circuits required to evaluate $C_{\rm
fsVFF}$ in general scales with $n_{\rm eig}$. Given these two restrictions,
fsVFF is limited to simulating quantum states spanning a non-exponential
number of eigenstates. Consequently, we advocate using fsVFF to simulate
states with $n_{\rm eig}=\text{poly}(n)$. Crucially these states need not be
low lying and therefore our approach is more widely applicable than the
Subspace Variational Quantum Simulator (SVQS) algorithm Heya  _et al._ (2019),
which simulates fixed low energy input states. In Section V.3 we develop
methods for reducing the resources required to evaluate the cost for larger
values of $n_{\rm eig}$.

While $C_{\rm fsVFF}$ was motivated as a natural choice of cost function to
learn the evolution induced by a target unitary on a fixed initial state, it
is a global cost Cerezo  _et al._ (2021) and hence it encounters what is known
as a barren plateau for large simulation sizes McClean  _et al._ (2018);
Cerezo  _et al._ (2021); Cerezo and Coles (2021); Arrasmith  _et al._ (2020a);
Holmes  _et al._ (2020, 2021a); Volkoff and Coles (2021); Sharma  _et al._
(2020c); Pesah  _et al._ (2020); Uvarov and Biamonte (2020); Marrero  _et al._
(2020); Patti  _et al._ (2020). In Appendix C we suggest an alternative local
version of the cost to mitigate such trainability issues.

###  III.2 Calculating $n_{\rm eig}$

In this section, we present an algorithm to calculate $n_{\rm eig}$ and
therefore determine the number of training states required to evaluate $C_{\rm
fsVFF}$. Our proposed algorithm utilizes the fact that the number of energy
eigenstates spanned by $\ket{\psi_{0}}$ is equivalent to the number of
linearly independent states in the set $\mathcal{V}_{\infty}$ where
$\mathcal{V}_{k}:=\\{\ket{\psi_{l}}\\}_{l=0}^{l=k}$ with
$\ket{\psi_{l}}:=U(\Delta t)^{l}\ket{\psi_{0}}$. The subspace
$\mathcal{K}_{k}(U,\psi_{0})$ spanned by $\mathcal{V}_{k}$ is known as the
Krylov subspace associated with the operator $U$ and vector $\ket{\psi_{0}}$
Krylov (1931). Therefore, $n_{\rm eig}$ is equivalently the dimension of the
Krylov subspace $\mathcal{K}_{\infty}(U,\psi_{0})$.

To determine the dimension of $\mathcal{K}_{\infty}(U,\psi_{0})$ we can
utilize the fact that the determinant of the Gramian matrix of a set of
vectors (i.e., the matrix of their overlaps) is zero if and only if the
vectors are linearly dependent. The Gramian corresponding to $\mathcal{V}_{k}$
is given by

$G(k)=\begin{pmatrix}\braket{\psi_{0}}{\psi_{0}}&\braket{\psi_{0}}{\psi_{1}}&\cdots&\braket{\psi_{0}}{\psi_{k}}\\\
\braket{\psi_{1}}{\psi_{0}}&\braket{\psi_{1}}{\psi_{1}}&\cdots&\braket{\psi_{1}}{\psi_{k}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\braket{\psi_{k}}{\psi_{0}}&\braket{\psi_{k}}{\psi_{1}}&\cdots&\braket{\psi_{k}}{\psi_{k}}\end{pmatrix}\,.$
(7)

If $\text{Det}(G(k))\neq 0$, then the vectors in $\mathcal{V}_{k}$ are
linearly independent and therefore span at least a $k+1$ dimensional subspace.
Conversely, if $\text{Det}(G(k))=0$, the set $\mathcal{V}_{k}$ contains linear
dependencies and the subspace they span is less than $k+1$ dimensional.
Therefore, if we can find $k_{\rm min}$, the smallest $k$ such that
$\text{Det}(G(k))=0$, then (noting that $G(k)$ is a $k+1$ dimensional matrix)
we know that $k_{\rm min}$ is the largest number of linearly independent
vectors spanned by $\mathcal{V}_{\infty}$. That is, $k_{\rm min}$ is the
dimension of $\mathcal{K}_{\infty}(U,\psi_{0})$ and so we have that $n_{\rm
eig}=k_{\rm min}$.

The overlaps $\braket{\psi_{l}}{\psi_{l^{\prime}}}$ for any $l$ and
$l^{\prime}$ can be measured using the Hadamard Test, shown in Fig. 1, and
thus the Hadamard test can be used to determine $G(k)$ on quantum hardware.
Since the Gramian here contains two symmetries, hermiticity and the invariance
$\braket{\psi_{l}}{\psi_{l^{\prime}}}=\langle\psi_{0}|U^{-l}U^{l^{\prime}}|\psi_{0}\rangle=\langle\psi_{0}|U^{{l^{\prime}}-l}|\psi_{0}\rangle=\langle\psi_{0}|\psi_{{l^{\prime}}-l}\rangle$,
we only have to calculate the first row of the matrix $G(k)$ on the quantum
computer.

In summary, our proposed algorithm to determine $n_{\rm eig}$ consists of the
following loop. Starting with $k=1$,

  1. 1.

Construct $G(k)$ using the Hadamard test.

  2. 2.

Calculate (classically) $\text{Det}(G(k))$.

If $\text{Det}(G(k))=0$, terminate the loop and conclude that $n_{\rm eig}=k$.

If $\text{Det}(G(k))\neq 0$, increase $k\rightarrow k+1$ and return to step 1.

This is shown schematically in Fig. 1.

We remark that in the presence of degeneracies in the spectrum of $H$, the
eigenvectors corresponding to degenerate eigenvalues are not unique.
Therefore, in this case, the number of states spanned by $\ket{\psi_{0}}$
depends on how the eigenvectors corresponding to degenerate eigenvalues are
chosen. However, as detailed in Appendix A, to learn the action of $U$ on
$\ket{\psi_{0}}$, what matters is the number of eigenstates spanned by
$\ket{\psi_{0}}$ corresponding to unique eigenvalues. This is equivalent to
the dimension of the Krylov subspace $\mathcal{K}_{\infty}(U,\psi_{0})$.
Consequently, the algorithm detailed above can also be used in this case.

While it is beneficial to learn $n_{\rm eig}$ to determine how many training
states are required to perfectly learn the diagonalization on the subspace
spanned by the initial state, we stress that it is not strictly necessary for
the successful implementation of fsVFF. One could always train on an
increasing number of states and study the convergence of an observable of
interest. More concretely, one could train on $k$ states and then use the
resultant diagonalization to compute the evolution of a particular observable
as a function of time. For $k<n_{\rm eig}$ the trajectory of the observable
will alter as $k$ is increased. However, for $k\geqslant n_{\rm eig}$
increasing $k$ further will no longer change the trajectory of the observable
because it will have already converged on the true trajectory. Using this
approach, $n_{\rm eig}$ need not be already known to implement fsVFF.

###  III.3 Ansatz

The fsVFF algorithm, similarly to VFF, employs an ansatz of the form

$\displaystyle V(\bm{\alpha},\Delta t)=W(\bm{\theta})D(\bm{\gamma},\Delta
t)W^{\dagger}(\bm{\theta})\,,$ (8)

to diagonalize the initial Trotter unitary $U(\Delta t)$. Here
$W(\bm{\theta})$ is a quantum circuit that approximately rotates the standard
basis into the eigenbasis of $H$, and $D(\bm{\gamma})$ is a diagonal unitary
that captures the (exponentiated) eigenvalues of $H$. A generic diagonal
operator $D$ can be written in the form

$D(\bm{\gamma},\Delta t)=\prod_{\bm{q}}e^{i\gamma_{\bm{q}}Z^{\bm{q}}\Delta
t}\,,$ (9)

where $\gamma_{\bm{q}}\in\mathbb{R}$ and we use the notation

$Z^{\bm{q}}=Z_{1}^{q_{1}}\otimes\cdots\otimes Z_{n}^{q_{n}}\,,$ (10)

with $Z_{j}$ the Pauli $Z$ operator acting on qubit $j$. While Eq. (9)
provides a general expression for a diagonal unitary, for practical ansätze it
may be desirable to assume that the $Z^{\bm{q}}$ operators are local operators
and the product contains a polynomial number of terms, i.e., is in
$\mathcal{O}(\text{poly}(n))$. There is more flexibility in the construction
of the ansätze for $W$ since these are generic unitary operations. A natural
choice might be to use a hardware-efficient ansatz Kandala  _et al._ (2017) or
an adaptive ansatz Khatri  _et al._ (2019); Bilkis  _et al._ (2021).

One of the main advantages of fsVFF is that diagonalization is only necessary
over the subspace spanned by the initial state, rather than the entire Hilbert
space which will be significantly larger. To outperform standard VFF, it is in
our interest to take advantage of this small subspace to find compact ansätze.

The two main impeding factors we wish to minimize to aid diagonalization are
error rates and optimization time. Therefore, when searching for ansätze, our
priorities are to minimize the number of CNOT gates required (the noisiest
component in the ansätze) and the number of rotation parameters. There is,
however, a trade off between expressibility of the ansatz and its
trainability. There needs to be enough freedom in the unitary to map the
required eigenvectors to the computational basis but generically highly
expressive ansätze exhibit barren plateaus Holmes  _et al._ (2021a).

For systems with symmetries and/or systems that are nearby perturbations of
known diagonalizable systems, it may be possible to find a fully expressive,
compact ansatz by inspection. This is the case for a simple 2-qubit XY
Hamiltonian, as discussed in Section IV.

More generally, it can be challenging to analytically find compact but
sufficiently expressible ansätze. Nonetheless, it is possible to variationally
update the ansatz structure and thereby systematically discover simple
structures. One straightforward approach is to use a layered ansatz where each
layer initializes to the identity gate Grant  _et al._ (2019); Skolik  _et
al._ (2020). The ansatz can be optimized until it plateaus, redundant single
qubit gates removed, then another layer can be appended and the process
repeats. Alternatively, more sophisticated discrete optimization techniques
may be used to variationally search the space of ansätze.

###  III.4 Summary of algorithm

The fixed state Variational Fast Forwarding algorithm (fsVFF) is summarized in
Fig. 1. We start with an initial state $\ket{\psi_{0}}$ that we wish to evolve
under the Hamiltonian $H$.

  1. 1.

The first step is to approximate the short time evolution using a single step
Trotter approximation $U$.

  2. 2.

This Trotter approximation can be used to find an approximation for $n_{\rm
eig}$, the dimension of the energy eigenspace spanned by $\ket{\psi_{0}}$,
using the method outlined in Section III.2.

  3. 3.

Equipped with a value for $n_{\rm eig}$, we then variationally search for a
diagonalization of $U$ over the energy subspace spanned by $\ket{\psi_{0}}$
using $C_{\rm fsVFF}$, Eq. (5). At each iteration step the gradient of the
cost with respect to a parameter $\theta_{i}$ is measured on the quantum
computer for a fixed set of parameters using the analytic expressions for
$\partial_{\theta_{i}}C_{\rm fsVFF}$ provided in Appendix D. These gradients
are used to update the parameters using a classical optimizer, such as those
in Refs. Kübler  _et al._ (2020); Arrasmith  _et al._ (2020b); Sweke  _et al._
(2020). The output of the optimization loop is the set of parameters that
minimize $C_{\rm fsVFF}$,

$\displaystyle\\{\bm{\theta}_{\rm opt},\bm{\gamma}_{\rm
opt}\\}=\underset{\bm{\theta},\bm{\gamma}}{\text{ arg min }}C_{\rm
fsVFF}(\bm{\theta},\bm{\gamma})\,.$ (11)

  4. 4.

Finally, the state $\ket{\psi_{0}}$ can be simulated for time $T=N\Delta t$
using the circuit

$\displaystyle W(\bm{\theta}_{\rm opt})D(\bm{\gamma}_{\rm opt},N\Delta
t)W(\bm{\theta}_{\rm opt})^{\dagger}\,.$ (12)

That is, by simply multiplying the parameters $\bm{\gamma}_{\rm opt}$ in the
diagonalized unitary by a constant number of iterations $N$.

In Appendix E, we show that the total simulation fidelity, in the limit that
leakage is small, is expected to scale sub-quadratically with the number of
fast-forwarding time steps $N$. Thus, if the minimal cost from the
optimization loop is sufficiently small, we expect the fsVFF algorithm to
allow for long, high fidelity simulations.

##  IV Hardware Implementation

In this section we demonstrate that fsVFF can be used to implement long time
simulations on quantum hardware. Specifically, we simulate the XY spin chain,
which has the Hamiltonian

$H_{\rm XY}:=\sum_{j=1}^{n-1}X_{j}X_{j+1}+Y_{j}Y_{j+1}\,,$ (13)

where $X_{j}$ and $Y_{j}$ are Pauli operators on the $j_{\rm th}$ qubit. In
what follows, we first present results showing that we can determine $n_{\rm
eig}$ for an initial state $\ket{\psi_{0}}$ using the method described in
Section III.2. We then demonstrate that the fsVFF cost can be trained to find
an approximate diagonalization of $H_{\rm XY}$ on the subspace spanned by
$\ket{\psi_{0}}$. We finally use this diagonalization to perform a long time
fast forwarded simulation. In all cases we focus on a two qubit chain, i.e.
$n=2$, and we approximate its evolution operator using a first-order Trotter
approximation.

###  IV.1 Determining $n_{\rm eig}$

Figure 2: Gramian Determinant Calculation. Here we plot the determinant of the
Gramian matrix, $\text{Det}(G)$, for $G$ measured on the Honeywell quantum
computer (solid) and simulated classically (dashed) for a 2-qubit XY spin
chain. Specifically we looked at states spanning $k=1$ (blue), $k=2$ (yellow)
and $k=3$ (red) eigenstates. For both sets of data $\text{Det}(G(n_{\rm
eig}))\approx 0$, demonstrating the effectiveness of the method for
determining $n_{\rm eig}$ that we introduce in Section III.2. For the
Honeywell implementation we used 1000 measurement samples per circuit.

The 2-qubit XY Hamiltonian has the eigenvectors
$\\{\ket{00},\frac{1}{\sqrt{2}}(\ket{10}+\ket{01}),\frac{1}{\sqrt{2}}(\ket{10}-\ket{01}),\ket{11}\\}$,
corresponding to the eigenvalues {0, 1, -1, 0}. As proof of principle, we
tested the algorithm for determining $n_{\rm eig}$ on the states $\ket{00}$
(corresponding to $n_{\rm eig}=1$), $\ket{10}$ ($n_{\rm eig}=2$) and
$\frac{1}{\sqrt{2}}(\ket{00}+\ket{10})$ ($n_{\rm eig}=3$). As described in
Section III.2, the $n_{\rm eig}$ of these states can be found by calculating
$\text{Det}(G(k))$ for increasing values of $k$ since, as $k$ is increased,
the determinant first equals 0 when $k=n_{eig}$.

To verify this for the states considered here, we first determine $G$ using a
classical simulator. As seen in Figure 2, in this case $\text{Det}(G(k))$
exactly equals 0 when $k=n_{eig}$. We then measured $G$ on Honeywell’s quantum
computer. Although on the real quantum device gate noise and sampling errors
are introduced, the results reproduce the classical results reasonably well.
Namely, at the correct value of $k$, $\text{Det}(G(k))$ drastically reduces
and approximately equals 0. Thus, we have shown that it is possible to
determine $n_{\rm eig}$ for an initial state by measuring $G$ on quantum
hardware.

###  IV.2 Training

We tested the training step of the algorithm on IBM and Rigetti’s quantum
computers, specifically ibmq_toronto and Aspen-8. For the purpose of
implementing a complete simulation, we chose to focus on simulating the
evolution of the state $|\psi_{0}\rangle=|10\rangle$. As discussed in the
previous section, this state spans $n_{\rm eig}=2$ eigenstates.

To diagonalize $H_{\rm XY}$ on the 2-dimensional subspace spanned by
$|10\rangle$, we used a hybrid quantum-classical optimization loop to minimize
$C_{\rm fsVFF}$. For a state with $n_{\rm eig}=2$ the cost $C_{\rm fsVFF}$,
Eq. (5), uses two training states $\\{U(\Delta
t)^{k}|\psi_{0}\rangle\\}_{k=1,2}$ where $U(t)$ is the first-order Trotter
approximation of $H_{\rm XY}$. On the IBM quantum computer we evaluated the
full cost function for each gradient descent iteration. However, the time
available on the Aspen-8 device was limited, so to speed up the rate of
optimization we evaluated the overlap on just one of the two training states
per iteration, alternating between iterations (instead of evaluating the
overlaps on both training states every iteration). To allow the movement
through parameter space to use information averaged over the two timesteps,
whilst only using a single training state per cost function evaluation,
momentum was added to the gradient updates Defazio (2020).

To take advantage of the fact that more compact ansätze are viable for fsVFF,
we variationally searched for a short depth ansatz, tailored to the target
problem. Specifically, we started training with a general 2-qubit unitary and
then during training the structure was minimised by pruning unnecessary gates.
In Figure 3, we show the circuit for the optimal ansatz obtained using the
method. The ansatz requires one CNOT gate and two single qubit gates for $W$
and only one $R_{z}$ rotation for $D$. This is a substantial compression on
the most general two qubit ansatz for $W$ which requires 3 CNOTs and 15 single
qubit rotations and the most general 2 qubit ansatz for $D$ which requires 2
$R_{z}$ rotations and one 2-qubit ZZ rotation (though in the case of the XY
Hamiltonian this may be simplified to only 2 $R_{z}$ rotations Lieb  _et al._
(1961)).

Figure 3: Ansatz for Hardware Implementation. The ansatz used to diagonalize
the 2-qubit XY Hamiltonian in the subspace of initial state $\ket{10}$ for the
implementation on Rigetti and IBM’s quantum computers. Here
$R_{j}(\theta)=\exp(-i\theta\sigma_{j}/2)$ for $j=x,y,z$. Figure 4: Hardware
Implementation. a) The 2-qubit parameterised quantum circuit shown in Figure 3
was trained to diagonalize $U(\Delta t)$, a first order Trotter expansion of
the 2-qubit XY Hamiltonian with $\Delta t=0.5$, in the subspace spanned by
$|10\rangle$. The dashed line plots the noisy cost as measured on ibmq_toronto
(yellow) and Aspen-8 (red) using 30,000 samples per circuit. The solid line
indicates the equivalent noise-free cost that was calculated on a classical
simulator. b) The initial state $|\psi_{0}\rangle=|10\rangle$ is evolved
forwards in time on the ibmq_rome quantum computer using the iterated Trotter
method (blue) and using fsVFF with the optimum parameters found on
ibmq_toronto (yellow) and Aspen-8 (red). The quality of the simulation is
evaluated by plotting the fidelity $F=\langle\psi|\rho|\psi\rangle$ between
the evolved state and exact evolution. The grey dotted line at $F=0.25$
represents the overlap with the maximally mixed state. The black dotted line
denotes a threshold fidelity at $F=0.9$. The inset shows the fast-forwarding
of the ansatz trained on ibmq_toronto on a longer timescale, where the
fidelity dropped below 0.9 (0.8) at 625 (1275) timesteps. All simulation data
was taken using 1000 samples per circuit.

Figure 4(a) shows the fsVFF cost function versus the number of iterations for
the implementations on ibmq_toronto (yellow) and Aspen-8 (red). The dashed
line indicates the noisy cost value obtained from the quantum computer. To
evaluate the quality of the optimization, we additionally classically compute
the true cost (indicated by the solid lines) using the parameters found on
ibmq_toronto and Aspen-8. While the noisy cost saturates at around $10^{-1}$,
we obtained a minimum noise-free cost of the order $10^{-3}$. The two orders
of magnitude difference between the noisy and the noise-free cost is
experimental evidence that the cost function is noise resilient on extant
quantum hardware.

###  IV.3 Fast forwarding

Finally we took the two sets of parameters found from training on ibmq_toronto
and Aspen-8, and used them to implement a fast-forwarded simulation of the
state $\ket{10}$ on ibmq_rome. To evaluate the quality of the fast forwarding
we calculated the fidelity, $F(N)=\langle\psi(N)|\rho(N)|\psi(N)\rangle$,
between the density matrix of the simulated state, $\rho(N)$, after $N$
timesteps, and the exact time evolved state, $|\psi(N)\rangle$, at time
$T=N\Delta t$. We used Quantum State Tomography to reconstruct the output
density matrix, and then calculated the overlap with the exact state
classically.

As shown by the plots of $F(N)$ in Figure 4(b), fsVFF significantly
outperforms the iterated Trotter method. Let us refer to the time before the
simulation falls below an error threshold $\delta$ as the high fidelity time.
Then the ratio of the high fidelity time for fsVFF ($T^{\rm FF}_{\delta}$) and
for standard Trotterization ($T^{\rm Trot}_{\delta}$) is a convenient measure
of simulation performance,

$R^{\rm FF}_{\delta}=T^{\rm FF}_{\delta}/T^{\rm Trot}_{\delta}\,.$ (14)

A simulation can be said to have been successfully fast-forwarded if $R^{\rm
FF}_{\delta}>1$. The iterated Trotter method dropped below a simulation
infidelity threshold of $\delta=1-F=0.1\ (0.2)$ after 4 (8) timesteps. In
comparison, fsVFF maintained a high fidelity for 625 (1275) timesteps. Thus we
achieved a simulation fast-forwarding ratio of $R_{0.1}^{\rm FF}=156$
($R_{0.2}^{\rm FF}=159$).

##  V Numerical Simulations

###  V.1 Noisy Training

We further validate fsVFF’s performance by testing it on a simulator of a
noisy quantum computer. The noise levels on the simulator are lower than those
experienced on current devices and hence these results are indicative of the
performance of the algorithm in the near future as hardware improves.

For these numerics we diagonalize the evolution of the 4 qubit XY Hamiltonian,
in the subspace spanned by the domain wall state $\ket{\psi_{0}}=\ket{1100}$.
This state spans 5 energy eigenstates of the XY Hamiltonian and so we use the
training states $\\{U(\Delta t)^{k}\ket{\psi_{0}}\\}_{k=1}^{5}$. Here
$U(\Delta t)$ is chosen to be a second-order Trotter-Suzuki decomposition for
the evolution operator under $H_{XY}$ with $\Delta t=0.5$. The noise model
used was based upon the IBM architecture.

To construct the ansatz for the diagonalizing unitary, $W$, we developed an
adaptive technique, similar to that proposed in Khatri  _et al._ (2019);
Bilkis  _et al._ (2021), to evolve the discrete circuit structure, as well as
optimize the rotation parameters using gradient descent. This method tends to
produce shallower circuits than the ones obtained with fixed ansatz
approaches. It is also less prone to get stuck in local minima. Since, the XY
Hamiltonian is particle number conserving we further use only particle number
conserving gates. This reduces the number of parameters in $W$, as well as
minimizing the leakage out of the symmetry sector when the circuit is executed
with a noisy simulator. Additional details on this adaptive learning method
are provided Appendix F. The ansatz for $D$, as in our 2 qubit hardware
implementation, simply consisted of $R_{z}$ rotations on each qubit.

The result of the training is shown in the inset of Fig. 5. The noisy cost was
measured by the noisy quantum simulator, whereas the noise-free cost is
calculated simultaneously but in the absence of any noise. The significant
separation between the noisy and noise-free cost again demonstrates the noise
resilience of the VFF algorithm. After successfully training the cost, the
fast-forward performance was then evaluated. Using the same noise model, the
output density matrix of the Iterated-Trotter state and the fast-forwarded
state was compared against the the Iterated-Trotter state in the absence of
noise, with the fidelity between the two states plotted. As shown in Fig. 5,
the fast-forwarded evolution significantly out performs the Iterated-Trotter
evolution, with the former’s fidelity dropping below 0.8 after 700 timesteps,
compared to only 8 steps of the latter. Thus we achieved a fast forwarding
ratio of $R^{\text{FF}}_{0.8}=87.5$.

Figure 5: Noisy Training and Fast-Forwarding of the 4 qubit XY Hamiltonian.
The inset shows the cost curve as the ansatz is evolved and optimized to
diagonalize the 4 qubit XY Hamiltonian in the 5-dimensional subspace spanned
by initial state $\ket{1100}$. The final circuit found by the learning
algorithm for the diagonalizing unitary, $W$, had 50 CNOT gates. The main plot
evaluates the fast-forwarding performance of the trained ansatz, in comparison
to the Iterated-Trotter evolved state. The fidelity is calculated against the
ideal state found in simulation using the iterated Trotter method in the
absence of noise, $F(N)=\bra{\psi(N)}\rho\ket{\psi(N)}$ with
$\ket{\psi(N)}=U(\Delta t)^{N}\ket{\psi_{0}}$. The black dotted line
highlights a threshold value $F=0.8$. The gray dotted line at $F=1/2^{4}$
represents the overlap with the maximally mixed state.

###  V.2 Fermi-Hubbard model

Finally, to probe the scalability and the breadth of applicability of the
fsVFF algorithm we performed a larger (noiseless) numerical implementation of
the algorithm on the Fermi-Hubbard model. Specifically, we considered the 1D
Fermi-Hubbard Hamiltonian on an $L$-site lattice with open boundary
conditions:

$\begin{split}H_{\mathrm{FH}}=&-J\sum_{j=1}^{L-1}\sum_{\sigma=\uparrow,\downarrow}c_{j,\sigma}^{\dagger}c_{j+1,\sigma}\
+\ \mathrm{h.c.}\\\ &+U\sum_{j=1}^{L}n_{j,\uparrow}n_{j,\downarrow}\
.\end{split}$ (15)

Here, $c_{j,\sigma}$ ($c^{\dagger}_{j,\sigma}$) denotes fermionic creation
(annihilation) operator at site $j$ for each of the two spin states
$\sigma=\uparrow,\downarrow$ and
$n_{j,\sigma}=c^{\dagger}_{j,\sigma}c_{j,\sigma}$ is a particle number
operator. The total number of fermions with a spin $\sigma$ is given by
$N_{\sigma}=\sum_{j}n_{j,\sigma}$. The term with coefficient $J$ in Eq. (15)
represents a single-fermion nearest-neighbor hopping and the term with
coefficient $U$ introduces on-site repulsion. The Hamiltonian preserves
particle numbers $N_{\uparrow}$ and $N_{\downarrow}$.

In our numerical studies, we choose $L=4$ (which requires 8 qubits to
simulate) and $J=1$, $U=2$ as well as $N_{\uparrow}=N_{\downarrow}=2$ (half
filling). The initial state is chosen to be a superposition of
$n_{\mathrm{eig}}=5$ eigenvectors of $H_{\mathrm{FH}}$ in the particle sector
$N_{\uparrow}=N_{\downarrow}=2$. Similar to our noisy simulations of the XY
model, we utilize an adaptive ansatz for $W$ that is made out of gates that
preserve particle number $N_{\uparrow}$ and $N_{\downarrow}$. The ansatz for
$D$ takes the form of Eq. (9), where we only allow for single-$Z$ terms in Eq.
(10). For this numerical result, we trained using the full exponentiation of
the Hamiltonian as the evolution operator, with no Trotter error.

Figure 6: Training and Fast-Forwarding of the 8 qubit Hubbard Model. The inset
shows the cost as it is iteratively minimized using an adaptive ansatz.
Various quality diagonalizations are indicated by the colored circles. In the
main figure, we plot the fidelity between the simulated state and the exact
evolution as a function of time. The red, green and yellow lines denote fsVFF
simulations using the corresponding quality diagonalization shown in the
inset. In blue we plot the fidelity of the Iterated-Trotter simulation,
$F_{\text{exact}}(T)=\bra{\psi(T)}\rho_{\rm trot}\ket{\psi(T)}$ with
$\ket{\psi(T)}=e^{-iHT}\ket{\psi_{0}}$ and $\rho_{\rm trot}$ the simulated
iterated Trotter state.

In the inset of Fig. 6 we show the cost function as it is iteratively
minimized. We then test the performance on a noisy simulator based upon a
fully connected 8-qubit trapped-ion device Trout  _et al._ (2018). As shown in
Fig. 6, small final cost values typically require deeper circuits to achieve,
the optimum diagonalization to use depends on the length of time one wishes to
simulate. At short times, a larger final cost function value performs better
since this corresponds to a shorter ansatz which experiences less noise.
However, to simulate longer times, a higher quality diagonalization is
required, with the additional noise induced by increased circuit depth
resulting in a relatively small decrease in fidelity. As shown in Fig. 6, we
find that the fast forwarding corresponding to an optimized cost of $1.1\times
10^{-5}$ maintained a fidelity of greater than 0.8 for $T<800$. In contrast,
the iterated Trotter method drops below 0.8 for $T>4.6$ and hence we here
achieve a fast forwarding ratio of $R^{FF}_{0.8}=174$.

###  V.3 Randomized Training

While the fsVFF cost as stated in Eq. (5) has $n_{\rm eig}$ terms, this does
not necessarily mean that the number of circuits required to evaluate it also
scales with $n_{\rm eig}$. Analogous to mini-batch gradient descent methods
popular for the training of classical neural networks, we can use only a small
random selection of the total training dataset per gradient evaluation, yet
over the whole optimization the total training set will be fully explored many
times over. Therefore, instead of restricting ourselves to a discrete set of
training states, which requires setting the size of the training set to be
equal to or greater than $n_{eig}$, we can instead randomly select our
training states from a continuum. This has the added advantage that it is then
unnecessary to explicitly compute $n_{\rm eig}$.

This approach results in a modified cost function of the form

$\widetilde{C}_{\rm fsVFF}:=1-\frac{1}{|R|}\sum\limits_{r\in
R}|\bra{\psi_{0}}V(-rt_{\rm max})U(rt_{\rm max})\ket{\psi_{0}}|^{2}$ (16)

where $V(t)=WD(t)W^{\dagger}$ is the fsVFF ansatz, $U(t)$ is a Trotter
approximation for the short time unitary evolution, and the elements of the
set $R$ are randomly generated numbers from the interval $[-1,1]$. The
gradients of the cost function are smaller when the unitary acts close to the
identity operation, so to maintain stronger gradients it is advantageous for
the elements of $R$ to be slightly biased towards the edges of the interval.
Specifically, in our numerics to test this approach, the absolute magnitude of
$r$ was raised to the power of 0.75. Although this approach does not require
an a priori calculation of $n_{\rm eig}$, there is a caveat that $t_{max}$
needs to be large enough to get sufficient separation of the training states
so they are not functionally identical. This alternative training setup is
potentially a yet more NISQ friendly variant, as the unitary does not need to
be decomposed into the form $U(\Delta t)^{n_{eig}}$, as required in Eq. 5, and
therefore allows for shorter depth circuits.

Figure 7: Randomized Training. a) The 5 (6) qubit XY Hamiltonian with initial
state $\ket{11100}(\ket{111000})$ is diagonalized using the cost function Eq.
(16), using only 2 training states per cost function evaluation to learn the
evolution within the 9 (12) dimensional subspace. The final circuit found by
the learning algorithm for the diagonalizing unitary, $W$, had 32 (134) CNOT
gates. b) After completion of the randomised training, the Hamiltonians were
fast-forwarded, with the infidelity evaluated in comparison to the noiseless
Trotter-iterated state $U(\Delta t)^{N}\ket{\psi_{0}}$.

To demonstrate the viability of this batched training method we diagonalized
the 5 (6) qubit XY Hamiltonian with initial state $\ket{11100}(\ket{111000})$,
which has an $n_{eig}=9\ (12)$. For both training curves shown in Fig. 7, only
2 training states per cost function evaluation were used. In both cases, we
trained with the unitary $U(t/6)^{6}$ where $U$ was the second order Trotter-
Suzuki operator, and $t_{max}=1$. The cost was successfully minimised to
$10^{-5}$ in both cases, and a noiseless simulation error of less than
$10^{-2}$ was maintained for over 100 time steps on fast forwarding.

##  VI Energy estimation

Figure 8: Energy estimation circuits a)/b) show circuit diagrams depicting the
enhancement of QPE/QEE using fsVFF. A circuit depth reduction is achieved
through replacing $U(\Delta t)$ with $D(\gamma_{\rm opt},\Delta t)$, and
removing the need to prepare an eigenstate in favour of a computational basis
state, $\ket{v_{k}}=W^{\dagger}\ket{E_{k}}$. QPE relies on implementing
controlled unitaries of the form $U(\Delta t)^{2^{j}}$ and therefore replacing
these with $D(\gamma_{\rm opt},2^{j}\Delta t)$ results in an exponential
reduction in circuit depth.

The diagonalization obtained from the optimization stage of fsVFF,
$W(\bm{\theta}_{\rm opt})D(\bm{\gamma}_{\rm opt},\Delta t)W(\bm{\theta}_{\rm
opt})^{\dagger}$, implicitly contains approximations of the eigenvalues and
eigenvectors of the Hamiltonian of the system of interest. In this section we
discuss methods for extracting the energy eigenstates and eigenvalues from a
successful diagonalization and implement them on quantum hardware.

The energy eigenvectors spanned by the initial state $\ket{\psi_{0}}$ can be
determined by the following simple sampling method. The first step is to apply
$W^{\dagger}$ to the initial state $\ket{\psi_{0}}$. In the limit of perfect
learning and vanishing Trotter error, this gives

$W(\bm{\theta}_{\rm
opt})^{\dagger}\ket{\psi_{0}}=\sum_{k=1}^{n_{eig}}a_{k}\ket{v_{k}}$ (17)

where $a_{k}=\langle E_{k}|\psi_{0}\rangle$ and
$\\{\ket{v_{k}}\\}_{k=1}^{n_{eig}}$ is a set of computational basis states.
The energy eigenstates spanned by $\ket{\psi_{0}}$ are then found by applying
$W(\bm{\theta}_{\rm opt})$ to any of the states obtained from measuring
$W(\bm{\theta}_{\rm opt})^{\dagger}\ket{\psi_{0}}$ in the computational basis,
that is $\\{\ket{E_{k}}\\}_{k=1}^{n_{\rm eig}}=\\{W(\bm{\theta}_{\rm
opt})\ket{v_{k}}\\}_{k=1}^{n_{\rm eig}}$.

Extracting the energy eigenvalues from $D$ is more subtle. Firstly, as
$WDW^{\dagger}$ and $U$, even in the limit of perfectly minimizing the cost
$C_{\rm fsVFF}$, may disagree by a global phase $\phi$, at best we can hope to
learn the difference between, rather than absolute values, of the energy
eigenvalues of $H$. For simple cases, where the diagonal ansatz $D$ is
composed of a polynomial number of terms, these energy value differences may
be extracted directly by rewriting $D$ in the computational basis. For
example, in our hardware implementation
$D(\gamma)=\exp\left(-i\frac{\gamma\Delta tZ_{1}}{2}\right)\otimes\openone$
and therefore the difference in energy between the two eigenvalues spanned by
$\ket{\psi_{0}}$ is given by $\gamma_{\rm opt}+\frac{k\pi}{\Delta t}$. Here
$k$ is an integer correcting for the arbitrary phase arising from taking the
$\log$ of $D$ that can be determined using the method described in Ref.
Commeau  _et al._ (2020). Using this approach, we obtain 1.9995 and 2.0019
from the training on IBM and Rigetti respectively, in good agreement with the
theoretically expected value of 2. For more complex cases, this simple post-
processing method will be become intractable and an algorithmic approach will
be necessary.

Figure 9: Quantum Phase Estimation. Using the 2-qubit diagonalization found
from training on ibmq_toronto, QPE was performed on ibmq_boeblingen on the
eigenvector $\ket{E_{1}}:=\frac{1}{\sqrt{2}}(\ket{01}+\ket{10})$. A phase of
$e^{\frac{2\pi i}{8}}$ is applied, so the measured output should be 001 with
probability 1. The variation distance from the target probability distribution
when using QPE with fsVFF was 0.578, compared to 0.917 using standard QPE.

Quantum Phase Estimation (QPE) Nielsen and Chuang (2000) and Quantum
Eigenvalue Estimation (QEE) Somma (2020) are fault tolerant quantum algorithms
for estimating the eigenvalues of a unitary operation. However, their
implementation on current quantum devices is limited by the reliance on the
execution of controlled unitaries from ancillary qubits. These controlled
unitaries require many entangling gates, and introduce too much noise to be
realized for large scale systems on current hardware. Once an evolution
operator has been diagonalized in the subspace of an initial state, fsVFF can
be used to significantly reduce the circuit depth of QPE and QEE, as shown in
Figure 8. In this manner, fsVFF provides a NISQ friendly means of estimating
the eigenvalues within a subspace of a Hamiltonian.

To demonstrate the power of fsVFF to reduce the depth of QPE, we perform QPE
using the diagonalization obtained from training on IBM’s quantum computer.
Specifically, we consider the input eigenvector
$\ket{E_{1}}:=\frac{1}{\sqrt{2}}(\ket{01}+\ket{10})$. This is one of the
eigenvectors spanned by the input state of our earlier hardware
implementation, $\ket{\psi_{0}}=\ket{10}$. We then consider evolving
$\ket{E_{1}}$ under $H_{\rm XY}$ for a time step of $\Delta t=1/8$. Since the
energy of the state $\ket{E_{1}}$ equals 1, we expect this to result in a
phase shift of $e^{2\pi i/8}$ being applied to $\ket{E_{1}}$. We implemented
QPE and fsVFF enhanced QPE to measure this phase using the circuits shown in
Fig 8. We chose to measure to 3 bits of precision and therefore the output
should be the measurement 001 with probability one. As Figure 9 shows, it
appears that the standard QPE implementation was unable to discern this phase.
In contrast, when fsVFF was used to reduce the circuit depth, the output
distribution was strongly peaked at the correct state.

Figure 10: Determining the eigenstates spanned by the initial state: The
3-qubit XY Hamiltonian was diagonalized on a quantum simulator in the subspace
of initial state $\ket{\psi_{0}}=\ket{110}$ to obtain $\bm{\theta}_{\rm opt}$
and $\bm{\gamma}_{\rm opt}$. Here we show the output of measuring
$W(\bm{\theta}_{\rm opt})^{\dagger}\ket{\psi_{0}}$ in the computational basis
on ibmq_boeblingen. The 4 non-zero states correspond to the 4 eigenvectors
spanned by $|\psi_{0}\rangle$. Figure 11: Eigenvalue Estimation using QEE.
Here we show the result of implementing QEE (using fsVFF as a pre-processing
step) on ibmq_santiago to calculate the eigenvalues of the eigenvectors in the
subspace spanned by $\ket{011}$. The solid yellow, red, blue and green lines
represent the eigenvalues obtained for the $\ket{000}$, $\ket{001}$,
$\ket{100}$ and $\ket{101}$ states, with exact corresponding energies of
{-2.828, 0, 0, 2.828}, indicated by the dotted lines. The eigenvalues are
plotted as phases since for $\Delta t=1$ there is a one to one correspondence.

QEE requires only one ancillary qubit, a single implementation of $e^{-iHt}$,
and no Quantum Fourier Transform and therefore is less resource intensive than
QPE. Nonetheless, we can again, as shown in Fig. 8, use fsVFF as a pre-
processing step to reduce the circuit depth.

We tested this on the 3-qubit XY Hamiltonian by first performing fsVFF on a
quantum simulator with the initial state $\ket{\psi_{0}}=\ket{110}$. Having
obtained an approximate diagonalization, we determined the eigenstates using
the sampling method described earlier. Figure 10 shows the results of the
measurement of $W(\bm{\theta}_{\rm opt})^{\dagger}\ket{\psi_{0}}$, with four
strong peaks corresponding to the four eigenvectors in this subspace.

Figure 11 shows the results of QEE implemented on ibmq_boeblingen. We use the
basis states found from the sampling method as our inputs to reduce the depth
of the circuit, and remove the need to use the time-series method originally
proposed for extracting the eigenvalues, as we could calculate the eigenvalues
individually by inputting their corresponding eigenvectors. A value of $\Delta
t=1$ was used so the phase calculated directly matched the eigenvalue. After
removing a global phase, QEE had accurately found the eigenvalues of the four
eigenvectors, with a mean-squared error from the true values of $4.37\times
10^{-3}$.

##  VII Discussion

In this work, we demonstrated that despite the modest size and noise levels of
the quantum hardware that is currently available, it is possible to perform
long time dynamical simulations with a high fidelity. Specifically, we have
introduced fsVFF, a new algorithm for NISQ simulations, which we used to
simulate a 2-qubit XY-model spin chain on the Rigetti and IBM quantum
computers. We achieved a fidelity of at least 0.9 for over 600 time steps.
This is a 150-fold improvement on the standard iterated Trotter approach,
which had a fidelity of less than 0.9 after only 4 time steps. Moreover, our
numerical simulations of the 4 qubit XY model and 8 qubit Fermi-Hubbard model
achieved fast-forwarding ratios of 87.5 and 174 respectively, indicating the
viability of larger implementations in the near future as hardware improves.

Central to the success of the fsVFF algorithm is the fact that it is tailored
to simulating a particular fixed initial state rather than an arbitrary
initial state. By focusing on this less demanding task, we showed that it is
possible to substantially reduce the width and depth of the previously
proposed VFF algorithm. In particular, since fsVFF only requires finding a
diagonalization of a short-time evolution unitary on the subspace spanned by
the initial state (compared to the entire Hilbert space in the case of VFF),
fsVFF can utilize much simpler ansätze. This is demonstrated in our hardware
implementation, where one CNOT and two parameterized single qubit rotations
proved sufficient for an effective ansatz for $W$ and one single qubit
rotation was sufficient for $D$.

The fsVFF algorithm, similarly to VFF, is fundamentally limited by the initial
Trotter error approximating the short time evolution of the system. The
Variational Diagonalization Hamiltonian (VHD) algorithm Commeau  _et al._
(2020) may be used to remove this error. However, like VFF, VHD is designed to
simulate any possible initial state. There are a number of different
approaches inspired by fsVFF that could be explored for reducing the resource
requirements of the VHD algorithm by focusing on simulating a particular
initial state. Such a “fixed state VHD” algorithm would allow for more
accurate long time simulations on NISQ hardware.

More generally, our work highlights the trade off between the universality of
an algorithm and the resources required to implement it. One can imagine a
number of alternative ways in which the universality of an algorithm can be
sacrificed, without significantly reducing its utility, in order to make it
more NISQ friendly. For example, one is often interested in studying the
evolution of a particular observable of interest, rather than all possible
observables. It would be interesting to investigate whether a fixed-observable
fsVFF could further reduce the resources required to implement long time high
fidelity simulations. More broadly, an awareness of this trade off may prove
useful beyond dynamical simulation for the ongoing challenge of adapting
quantum algorithms to the constraints of NISQ hardware.

##  VIII Acknowledgements

JG and KG acknowledge support from the U.S. Department of Energy (DOE) through
a quantum computing program sponsored by the Los Alamos National Laboratory
(LANL) Information Science & Technology Institute. ZH, BC and PJC acknowledge
support from the LANL ASC Beyond Moore’s Law project. ZH acknowledges
subsequent support from the Mark Kac Fellowship. We acknowledge the LANL
Laboratory Directed Research and Development (LDRD) program for support of AS
and initial support of BC under project number 20190065DR as well as LC under
project number 20200022DR. AA was supported by the U.S. Department of Energy
(DOE), Office of Science, Office of High Energy Physics QuantISED program
under Contract No. DE-AC52-06NA25396. LC and PJC were also supported by the
U.S. DOE, Office of Science, Basic Energy Sciences, Materials Sciences and
Engineering Division, Condensed Matter Theory Program. This research used
quantum computing resources provided by the LANL Institutional Computing
Program, which is supported by the U.S. Department of Energy National Nuclear
Security Administration under Contract No. 89233218CNA000001. This research
used additional quantum computational resources supported by the LANL ASC
Beyond Moore’s Law program and by the Oak Ridge Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under Contract DE-
AC05-00OR22725.

## References

  * Arute  _et al._ (2019) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, _et al._ , “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505–510 (2019).
  * Arute  _et al._ (2020) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Andreas Bengtsson, Sergio Boixo, Michael Broughton, Bob B Buckley, _et al._ , “Observation of separated dynamics of charge and spin in the fermi-hubbard model,” arXiv preprint arXiv:2010.07965 (2020).
  * Lloyd (1996) Seth Lloyd, “Universal quantum simulators,” Science , 1073–1078 (1996).
  * Sornborger and Stewart (1999) AT Sornborger and Ewan D Stewart, “Higher-order methods for simulations on quantum computers,” Physical Review A 60, 1956 (1999).
  * Low and Chuang (2019) Guang Hao Low and Isaac L Chuang, “Hamiltonian simulation by qubitization,” Quantum 3, 163 (2019).
  * Berry  _et al._ (2015) Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma, “Simulating hamiltonian dynamics with a truncated taylor series,” Physical Review Letters 114, 090502 (2015).
  * Cerezo  _et al._ (2020a) M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles, “Variational quantum algorithms,” arXiv preprint arXiv:2012.09265 (2020a).
  * Endo  _et al._ (2021) Suguru Endo, Zhenyu Cai, Simon C Benjamin, and Xiao Yuan, “Hybrid quantum-classical algorithms and quantum error mitigation,” Journal of the Physical Society of Japan 90, 032001 (2021).
  * Bharti  _et al._ (2021) Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik, “Noisy intermediate-scale quantum (nisq) algorithms,” arXiv preprint arXiv:2101.08448 (2021).
  * Peruzzo  _et al._ (2014) Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature communications 5, 1–7 (2014).
  * Farhi  _et al._ (2014) Edward Farhi, Jeffrey Goldstone, and Sam Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028 (2014).
  * McClean  _et al._ (2016) Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik, “The theory of variational hybrid quantum-classical algorithms,” New Journal of Physics 18, 023023 (2016).
  * Khatri  _et al._ (2019) Sumeet Khatri, Ryan LaRose, Alexander Poremba, Lukasz Cincio, Andrew T Sornborger, and Patrick J Coles, “Quantum-assisted quantum compiling,” Quantum 3, 140 (2019).
  * LaRose  _et al._ (2019) Ryan LaRose, Arkin Tikku, Étude O’Neel-Judy, Lukasz Cincio, and Patrick J Coles, “Variational quantum state diagonalization,” npj Quantum Information 5, 1–10 (2019).
  * Arrasmith  _et al._ (2019) Andrew Arrasmith, Lukasz Cincio, Andrew T Sornborger, Wojciech H Zurek, and Patrick J Coles, “Variational consistent histories as a hybrid algorithm for quantum foundations,” Nature communications 10, 1–7 (2019).
  * Cerezo  _et al._ (2020b) M. Cerezo, Alexander Poremba, Lukasz Cincio, and Patrick J Coles, “Variational quantum fidelity estimation,” Quantum 4, 248 (2020b).
  * Li and Benjamin (2017) Y. Li and S. C. Benjamin, “Efficient variational quantum simulator incorporating active error minimization,” Phys. Rev. X 7, 021050 (2017).
  * Endo  _et al._ (2020) Suguru Endo, Jinzhao Sun, Ying Li, Simon C Benjamin, and Xiao Yuan, “Variational quantum simulation of general processes,” Physical Review Letters 125, 010501 (2020).
  * Yao  _et al._ (2020) Yong-Xin Yao, Niladri Gomes, Feng Zhang, Thomas Iadecola, Cai-Zhuang Wang, Kai-Ming Ho, and Peter P Orth, “Adaptive variational quantum dynamics simulations,” arXiv preprint arXiv:2011.00622 (2020).
  * Heya  _et al._ (2019) Kentaro Heya, Ken M Nakanishi, Kosuke Mitarai, and Keisuke Fujii, “Subspace variational quantum simulator,” arXiv preprint arXiv:1904.08566 (2019).
  * Cirstoiu  _et al._ (2020) 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–10 (2020).
  * Commeau  _et al._ (2020) Benjamin Commeau, M. Cerezo, Zoë Holmes, Lukasz Cincio, Patrick J Coles, and Andrew Sornborger, “Variational hamiltonian diagonalization for dynamical quantum simulation,” arXiv preprint arXiv:2009.02559 (2020).
  * Trout  _et al._ (2018) Colin J Trout, Muyuan Li, Mauricio Gutiérrez, Yukai Wu, Sheng-Tao Wang, Luming Duan, and Kenneth R Brown, “Simulating the performance of a distance-3 surface code in a linear ion trap,” New Journal of Physics 20, 043038 (2018).
  * Benedetti  _et al._ (2020) Marcello Benedetti, Mattia Fiorentini, and Michael Lubasch, “Hardware-efficient variational quantum algorithms for time evolution,” arXiv preprint arXiv:2009.12361 (2020).
  * Bharti and Haug (2020) Kishor Bharti and Tobias Haug, “Quantum assisted simulator,” arXiv preprint arXiv:2011.06911 (2020).
  * Lau  _et al._ (2021) Jonathan Wei Zhong Lau, Kishor Bharti, Tobias Haug, and Leong Chuan Kwek, “Quantum assisted simulation of time dependent hamiltonians,” arXiv preprint arXiv:2101.07677 (2021).
  * Haug and Bharti (2020) Tobias Haug and Kishor Bharti, “Generalized quantum assisted simulator,” arXiv preprint arXiv:2011.14737 (2020).
  * Poland  _et al._ (2020) Kyle Poland, Kerstin Beer, and Tobias J Osborne, “No free lunch for quantum machine learning,” arXiv preprint arXiv:2003.14103 (2020).
  * Sharma  _et al._ (2020a) Kunal Sharma, M. Cerezo, Zoë Holmes, Lukasz Cincio, Andrew Sornborger, and Patrick J Coles, “Reformulation of the no-free-lunch theorem for entangled data sets,” arXiv preprint arXiv:2007.04900 (2020a).
  * Sharma  _et al._ (2020b) Kunal Sharma, Sumeet Khatri, M. Cerezo, and Patrick J Coles, “Noise resilience of variational quantum compiling,” New Journal of Physics 22, 043006 (2020b).
  * Sim  _et al._ (2019) Sukin Sim, Peter D. Johnson, and Alán Aspuru-Guzik, “Expressibility and entangling capability of parameterized quantum circuits for hybrid quantum-classical algorithms,” Advanced Quantum Technologies 2, 1900070 (2019).
  * Nakaji and Yamamoto (2020) Kouhei Nakaji and Naoki Yamamoto, “Expressibility of the alternating layered ansatz for quantum computation,” arXiv preprint arXiv:2005.12537 (2020).
  * Holmes  _et al._ (2021a) Zoë Holmes, Kunal Sharma, M. Cerezo, and Patrick J Coles, “Connecting ansatz expressibility to gradient magnitudes and barren plateaus,” arXiv preprint arXiv:2101.02138 (2021a).
  * Sahinoglu and Somma (2020) Burak Sahinoglu and Rolando D Somma, “Hamiltonian simulation in the low energy subspace,” arXiv preprint arXiv:2006.02660 (2020).
  * Suzuki (1976) Masuo Suzuki, “Generalized trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems,” Communications in Mathematical Physics 51, 183–190 (1976).
  * Cerezo  _et al._ (2021) M. Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J Coles, “Cost function dependent barren plateaus in shallow parametrized quantum circuits,” Nature Communications 12, 1791 (2021).
  * McClean  _et al._ (2018) Jarrod R McClean, Sergio Boixo, Vadim N Smelyanskiy, Ryan Babbush, and Hartmut Neven, “Barren plateaus in quantum neural network training landscapes,” Nature communications 9, 1–6 (2018).
  * Cerezo and Coles (2021) Marco Cerezo and Patrick J Coles, “Higher order derivatives of quantum neural networks with barren plateaus,” Quantum Science and Technology 6, 035006 (2021).
  * Arrasmith  _et al._ (2020a) Andrew Arrasmith, M. Cerezo, Piotr Czarnik, Lukasz Cincio, and Patrick J Coles, “Effect of barren plateaus on gradient-free optimization,” arXiv preprint arXiv:2011.12245 (2020a).
  * Holmes  _et al._ (2020) Zoë Holmes, Andrew Arrasmith, Bin Yan, Patrick J Coles, Andreas Albrecht, and Andrew T Sornborger, “Barren plateaus preclude learning scramblers,” arXiv preprint arXiv:2009.14808 (2020).
  * Volkoff and Coles (2021) Tyler Volkoff and Patrick J Coles, “Large gradients via correlation in random parameterized quantum circuits,” Quantum Science and Technology 6, 025008 (2021).
  * Sharma  _et al._ (2020c) Kunal Sharma, M. Cerezo, Lukasz Cincio, and Patrick J Coles, “Trainability of dissipative perceptron-based quantum neural networks,” arXiv preprint arXiv:2005.12458 (2020c).
  * Pesah  _et al._ (2020) Arthur Pesah, M. Cerezo, Samson Wang, Tyler Volkoff, Andrew T Sornborger, and Patrick J Coles, “Absence of barren plateaus in quantum convolutional neural networks,” arXiv preprint arXiv:2011.02966 (2020).
  * Uvarov and Biamonte (2020) Alexey Uvarov and Jacob Biamonte, “On barren plateaus and cost function locality in variational quantum algorithms,” arXiv preprint arXiv:2011.10530 (2020).
  * Marrero  _et al._ (2020) Carlos Ortiz Marrero, Mária Kieferová, and Nathan Wiebe, “Entanglement induced barren plateaus,” arXiv preprint arXiv:2010.15968 (2020).
  * Patti  _et al._ (2020) Taylor L Patti, Khadijeh Najafi, Xun Gao, and Susanne F Yelin, “Entanglement devised barren plateau mitigation,” arXiv preprint arXiv:2012.12658 (2020).
  * Krylov (1931) AN Krylov, “On the numerical solution of equation by which are determined in technical problems the frequencies of small vibrations of material systems,” News Acad. Sci. USSR 7, 491–539 (1931).
  * Kandala  _et al._ (2017) Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, and Jay M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature 549, 242–246 (2017).
  * Bilkis  _et al._ (2021) M Bilkis, M Cerezo, Guillaume Verdon, Patrick J Coles, and Lukasz Cincio, “A semi-agnostic ansatz with variable structure for quantum machine learning,” arXiv preprint arXiv:2103.06712 (2021).
  * Grant  _et al._ (2019) Edward Grant, Leonard Wossnig, Mateusz Ostaszewski, and Marcello Benedetti, “An initialization strategy for addressing barren plateaus in parametrized quantum circuits,” Quantum 3, 214 (2019).
  * Skolik  _et al._ (2020) Andrea Skolik, Jarrod R McClean, Masoud Mohseni, Patrick van der Smagt, and Martin Leib, “Layerwise learning for quantum neural networks,” arXiv preprint arXiv:2006.14904 (2020).
  * Kübler  _et al._ (2020) Jonas M Kübler, Andrew Arrasmith, Lukasz Cincio, and Patrick J Coles, “An adaptive optimizer for measurement-frugal variational algorithms,” Quantum 4, 263 (2020).
  * Arrasmith  _et al._ (2020b) Andrew Arrasmith, Lukasz Cincio, Rolando D Somma, and Patrick J Coles, “Operator sampling for shot-frugal optimization in variational algorithms,” arXiv preprint arXiv:2004.06252 (2020b).
  * Sweke  _et al._ (2020) Ryan Sweke, Frederik Wilde, Johannes Jakob Meyer, Maria Schuld, Paul K Fährmann, Barthélémy Meynard-Piganeau, and Jens Eisert, “Stochastic gradient descent for hybrid quantum-classical optimization,” Quantum 4, 314 (2020).
  * Defazio (2020) Aaron Defazio, “Understanding the role of momentum in non-convex optimization: Practical insights from a lyapunov analysis,” arXiv preprint arXiv:2010.00406 (2020).
  * Lieb  _et al._ (1961) Elliott Lieb, Theodore Schultz, and Daniel Mattis, “Two soluble models of an antiferromagnetic chain,” Annals of Physics 16, 407–466 (1961).
  * Nielsen and Chuang (2000) Michael A. Nielsen and Isaac L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2000).
  * Somma (2020) Rolando D. Somma, “Quantum eigenvalue estimation via time series analysis,” (2020), arXiv:1907.11748 [quant-ph] .
  * Holmes  _et al._ (2021b) Zoë Holmes, Andrew Arrasmith, Bin Yan, Patrick J Coles, Andreas Albrecht, and Andrew T Sornborger, “Barren plateaus preclude learning scramblers,” Physical Review Letters 126, 190501 (2021b).
  * Anselmetti  _et al._ (2021) Gian-Luca R Anselmetti, David Wierichs, Christian Gogolin, and Robert M Parrish, “Local, expressive, quantum-number-preserving vqe ansatze for fermionic systems,” arXiv preprint arXiv:2104.05695 (2021).

##  Appendix A Faithfulness of cost function

Here we demonstrate that the cost function is faithful in the limit that
leakage from the subspace spanned by the initial state can be disregarded.
That is, suppose we could compile the exact evolution of the system
$U=\exp(-iH\Delta t)$ for a short timestep $\Delta t$. Then the cost function
vanishes,

$C_{\rm fsVFF}(U,V,\psi_{0})=1-\frac{1}{n_{\rm eig}}\sum_{k=1}^{n_{\rm
eig}}|\bra{\psi_{0}}{V^{\dagger}}^{k}U^{k}\ket{\psi_{0}}|^{2}=0\,,$ (18)

if and only if the fidelity of the fast-forwarded simulation is perfect,

$F_{\tau}=|\bra{\psi_{0}}{V^{\dagger}}^{\tau}U^{\tau}\ket{\psi_{0}}|^{2}=1\,,$
(19)

for all times $\tau$. Note, the reverse direction is trivial. If $F_{\tau}=1$
for all $\tau$ then $C_{\rm fsVFF}(U,V,\psi_{0})=0$.

Before embarking on the core of the proof let us first emphasize that in the
definition of the cost (18) we average over $n_{\rm eig}$ terms, where $n_{\rm
eig}$ is the number of eigenstates spanned by the initial state
$\ket{\psi_{0}}$ corresponding to unique eigenvalues of the Hamiltonian $H$.
The restriction to eigenstates with unique eigenvalues is important since it
is this which determines the dimension of the subspace spanned by the future
evolution of $\ket{\psi_{0}}$ (which in turn determines the number of training
states required to learn $U=\exp(-iH\Delta t)$).

To see this consider a Hamiltonian $H$ with a spectrum
$\\{E_{k}\\}_{k=1}^{2^{n}}$ and corresponding eigenstates
$\\{\ket{E_{k}}\\}_{k=1}^{2^{n}}$. The initial state $\ket{\psi_{0}}$ can be
expanded in the energy eigenbasis as

$\ket{\psi_{0}}=\sum_{k=1}^{m}a_{k}\ket{E_{k}}\,.$ (20)

The future evolution of such a state, i.e. the set of states

$\mathcal{S}_{k}(\psi_{0},H):=\\{e^{-iHj\Delta
t}\ket{\psi_{0}}\\}_{j=0}^{k}\,$ (21)

with $k=\infty$, is solely contained within the subspace spanned by
$\\{\ket{E_{k}}\\}_{k=1}^{m}$ since

$e^{-iHj\Delta t}\ket{\psi_{0}}=\sum_{k=1}^{m}a_{k}e^{-iE_{k}j\Delta
t}\ket{E_{k}}\,.$ (22)

If the spectrum $\\{E_{k}\\}_{k=1}^{2^{n}}$ is non-degenerate then the
evolution generates relative phases between all the $m=n_{\rm eig}$
eigenstates spanned by $\ket{\psi_{0}}$. In that case, $\mathcal{S}_{\infty}$
will span the entirety of the $m$ dimensional subspace spanned by
$\\{\ket{E_{k}}\\}_{k=1}^{m}$, i.e., an $m=n_{\rm eig}$ dimensional space.
However, suppose the eigenstate expansion of $\ket{\psi_{0}}$ includes two
degenerate eigenstates, i.e. two eigenstates that share the same eigenvalue.
In that case the evolution generates relative phases between only $m-1$ of the
eigenstates $\\{\ket{E_{k}}\\}_{k=1}^{m}$ and therefore $\mathcal{S}_{\infty}$
is confined to an $m-1$ dimensional subspace. More generally, if the
eigenstate expansion of $\ket{\psi_{0}}$ includes $n_{\rm eig}$ states with
unique eigenvalues, then evolution under $U=\exp(-iH\Delta t)$ generates
relative phases between $n_{\rm eig}$ states and so $\mathcal{S}_{\infty}$
will span an $n_{\rm eig}$ dimensional subspace of the space. In this manner,
it is the number of eigenstates spanned by the initial state $\ket{\psi_{0}}$
corresponding to unique eigenvalues, $n_{\rm eig}$, that determines the
subspace spanned by $\mathcal{S}_{\infty}$.

We note that any initial state $\ket{\psi_{0}}$ can be written in the form Eq.
(20) where, crucially, the sum is over only eigenstates corresponding to
unique energies, i.e. $m=n_{\rm eig}$. For a Hamiltonian with a non-degenerate
spectrum this expansion is trivial. For a Hamiltonian with a degenerate
spectrum there is some freedom in how the eigenstates corresponding to
degenerate eigenvalues are defined, since the superposition of degenerate
eigenstates is also an eigenstate corresponding to the same energy. Therefore,
henceforth, for degenerate Hamiltonians, we suppose that the energy eigenbasis
$\\{\ket{E_{k}}\\}_{k=1}^{2^{n}}$ is defined such that the eigenstate
expansion of $\ket{\psi_{0}}$ only contains eigenstates with unique energies,
that is $m=n_{\rm eig}$ terms.

We remark that our approach here may also be framed in the language of Krylov
spaces. The Krylov subspace Krylov (1931) associated with the operator $U$ and
vector $\ket{\psi_{0}}$ is the linear subspace spanned by the vectors
generated by evolving $\ket{\psi_{0}}$ under $U$ up to $k$ times. That is

$\mathcal{K}_{k}(U,\psi_{0})=\text{span}\\{\mathcal{V}_{k}\\}\ \ \ \
\text{where}\ \ \ \ \mathcal{V}_{k}:=\\{\ket{\psi_{l}}\\}_{l=0}^{l=k}\,,$ (23)

with $\ket{\psi_{l}}:=U^{l}\ket{\psi_{0}}$. Now supposing $U=e^{-iH\Delta t}$
and $\ket{\psi_{0}}=\sum_{k=1}^{n_{\rm eig}}a_{k}\ket{E_{k}}$, we have that
$\mathcal{V}_{k}=\mathcal{S}_{k}$. Thus the future evolution of
$\ket{\psi_{0}}$ is confined to the Krylov space
$\mathcal{K}_{\infty}(U,\psi_{0})$. This is an $n_{\rm eig}$ dimensional
subspace, and therefore, as will become clear, we require $n_{\rm eig}$
training states in order to learn $U$ on this subspace.

To prove the forward direction, we first note that if $C_{\rm fsVFF}=0$ then
as
$0\leqslant|\bra{\psi_{0}}{V^{\dagger}}^{k}U^{k}\ket{\psi_{0}}|^{2}\leqslant
1$, we have that $|\bra{\psi_{0}}{V^{\dagger}}^{k}U^{k}\ket{\psi_{0}}|^{2}=1$
for all $k$. It thus follows that the action of $V^{k}$ on $\ket{\psi_{0}}$
agrees with the action of $U^{k}$ on $\ket{\psi_{0}}$ up to an unknown phase
$e^{i\phi_{k}}$, i.e. for $1\leqslant k\leqslant n_{\rm eig}$ we have that

$\displaystyle V^{k}\ket{\psi_{0}}=e^{i\phi_{k}}U^{k}\ket{\psi_{0}}\,.$ (24)

Or equivalently, the action of $V$ and $U$ agree on the training states
$\mathcal{S}_{\rm train}:=\\{\ket{\psi_{k}}\\}_{k=0}^{n_{\rm eig}-1}$ up to an
unknown phase $e^{i\tilde{\phi}_{k}}$, that is

$\displaystyle
V\ket{\psi_{k}}=e^{i\tilde{\phi}_{k}}U\ket{\psi_{k}}\,\forall\,\ket{\psi_{k}}\,\in\,\mathcal{S}_{\rm
train},$ (25)

where $\tilde{\phi_{k}}=\phi_{k+1}-\phi_{k}$.

Now by construction (see Section III.2) the $n_{\rm eig}$ training states
$\mathcal{S}_{\rm train}$ are linearly independent. Furthermore, since the
initial simulation time $\Delta t$ may be chosen freely, the unitary
$U=\exp(-iH\Delta t)$ can be chosen such that none of the states in
$\mathcal{S}_{\rm train}$ are orthogonal. In this case, the unknown phases all
agree and we have that $\tilde{\phi_{k}}=\phi$ for all $k$. To see this note
that given $\ket{\psi_{1}}$ and $\ket{\psi_{0}}$ are linearly independent but
non-orthogonal, the state $\ket{\psi_{1}}$ can be represented as follows:

$\displaystyle\ket{\psi_{1}}=c\ket{\psi_{0}}+c_{\perp}\ket{\psi_{0}^{\perp}},$
(26)

where $|c|^{2}+|c_{\perp}|^{2}=1$ and $|c|^{2}\neq 0$. Then from (25), we find
that

$\displaystyle e^{-i\tilde{\phi}_{1}}$
$\displaystyle=\langle\psi_{1}|W|\psi_{1}\rangle$ (27)
$\displaystyle=|c|^{2}e^{i\tilde{\phi_{0}}}+|c_{\perp}|^{2}\langle\psi_{0}^{\perp}|W|\psi_{0}^{\perp}\rangle$
(28)
$\displaystyle=|c|^{2}e^{i\tilde{\phi_{0}}}+(1-|c|^{2})\langle\psi_{0}^{\perp}|W|\psi_{0}^{\perp}\rangle\,,$
(29)

where $W:=V^{\dagger}U$ (note that this use of $W$ is distinct from
$W(\bm{\theta})$ used in the main text for a parameterized eigenvector
unitary). The above expression can be rearranged as

$e^{-i\tilde{\phi}_{1}}-\langle\psi_{0}^{\perp}|W|\psi_{0}^{\perp}\rangle=|c|^{2}(e^{-i\tilde{\phi}_{0}}-\langle\psi_{0}^{\perp}|W|\psi_{0}^{\perp}\rangle)\,.$
(30)

Since $|e^{i\theta_{1}}|=1$ and $|c|^{2}\neq 0$ the aforementioned equation is
satisfied if and only if

$\displaystyle\langle\psi_{0}^{\perp}|W|\psi_{0}^{\perp}\rangle=e^{-i\tilde{\phi}_{1}}=e^{-i\tilde{\phi}_{0}}\,,$
(31)

which implies that $\tilde{\phi}_{0}=\tilde{\phi}_{1}$. Then by recursively
applying this procedure to the rest of the states in the training set, we find
that $\tilde{\phi}_{k}=\tilde{\phi}_{j}:=\phi$ for $0\leqslant k\leqslant
n_{\rm eig}-1$ and $0\leqslant j\leqslant n_{\rm eig}-1$ .

To understand the constraints from the minimization of the cost function, it
is helpful to consider the form of the unitary matrix $W$. It follows from Eq.
(25), and the fact that since $W$ is unitary
($\sum_{j}|W_{ij}|^{2}=\sum_{i}|W_{ij}|^{2}=1$), that $W$ can be represented
as Poland  _et al._ (2020):

$W=\left(\begin{array}[]{@{}c|c@{}}\begin{matrix}e^{i\phi}&\dots&0\\\
\vdots&\ddots&\\\ 0&&e^{i\phi}\end{matrix}&\mbox{\Large 0}\\\
\hline\cr\mbox{\Large 0}&W_{\perp}\\\ \end{array}\right).$

Here the upper left hand block spans the $n_{\rm eig}$ dimensional subspace
spanned by the input training states $\mathcal{S}_{\rm train}$ and $W_{\perp}$
is an unknown unitary matrix acting on the $(d-n_{\rm eig})$ dimensional space
orthogonal to $\mathcal{S}_{\rm train}$. For later convenience let us also
note that the matrix $\tilde{W}:=UV^{\dagger}$ is also a unitary matrix of the
form

$\tilde{W}=\left(\begin{array}[]{@{}c|c@{}}\begin{matrix}e^{i\phi}&\dots&0\\\
\vdots&\ddots&\\\ 0&&e^{i\phi}\end{matrix}&\mbox{\Large 0}\\\
\hline\cr\mbox{\Large 0}&\tilde{W}_{\perp}\\\ \end{array}\right),$

where $\tilde{W}_{\perp}$ is again a $(d-n_{\rm eig})$ dimensional unitary
matrix.

We are now in a position to show that if $C_{\rm fsVFF}(U,V,\psi_{0})=0$, and
by construction the states in $\mathcal{S}_{\rm train}$ are linearly
independent and non-orthogonal, then

$F_{\tau}=|\bra{\psi_{0}}{V^{\dagger}}^{\tau}U^{\tau}\ket{\psi_{0}}|^{2}=1\,,$
(32)

for all times $\tau$. To see this first note that for all times $\tau\leqslant
n_{\rm eig}$ that $F_{\tau}=1$ follows directly from Eq. (24). Now by
construction, any state $\ket{\psi_{\tau}}$ for $\tau>n_{\rm eig}$ linearly
depends on the states in $\mathcal{S}_{\rm train}$ and so it can be written as

$\ket{\psi_{\tau}}=\sum_{j=0}^{n_{\rm eig}}b^{(\tau)}_{j}\ket{\psi_{j}}\,,$
(33)

where $b^{(\tau)}_{j}=\braket{\psi_{j}}{\psi_{\tau}}$. Thus we have that

$W\ket{\psi_{\tau}}=\tilde{W}\ket{\psi_{\tau}}=e^{i\phi}\ket{\psi_{\tau}}$
(34)

for any time $\tau$. It straightforwardly follows from Eq. (24) and Eq. (34)
that the simulation fidelity at time $\tau=n_{\rm eig}+1$ equals 1,

$\displaystyle F_{n_{\rm eig}+1}$
$\displaystyle=|\bra{\psi_{0}}{V^{\dagger}}^{n_{\rm eig}+1}U^{n_{\rm
eig}+1}\ket{\psi_{0}}|^{2}$ (35) $\displaystyle=|\bra{\psi_{n_{\rm
eig}}}W\ket{\psi_{n_{\rm eig}}}|^{2}=1\,.$ (36)

Now, let us consider the simulation fidelity at time $\tau=n_{\rm eig}+2$,

$\displaystyle F_{n_{\rm eig}+2}$
$\displaystyle=|\bra{\psi_{0}}{V^{\dagger}}^{n_{\rm eig}+2}U^{n_{\rm
eig}+2}\ket{\psi_{0}}|^{2}$ (37) $\displaystyle=|\bra{\psi_{n_{\rm
eig}}}V^{\dagger}WU\ket{\psi_{n_{\rm eig}}}|^{2}$ (38)
$\displaystyle=|\bra{\psi_{n_{\rm
eig}}}U^{\dagger}\tilde{W}WU\ket{\psi_{n_{\rm eig}}}|^{2}$ (39)
$\displaystyle=|\bra{\psi_{n_{\rm eig}+1}}\tilde{W}W\ket{\psi_{n_{\rm
eig}+1}}|^{2}=1\,,$ (40)

where we have again used Eq. (24) and Eq. (34). Finally, the simulation
fidelity at an arbitrary time $\tau>n_{\rm eig}$ can be evaluated as follows,

$F_{\tau}=|\bra{\psi_{0}}{V^{\dagger}}^{\tau}U^{\tau}\ket{\psi_{0}}|^{2}=|\bra{\psi_{0}}(U^{\dagger}\tilde{W})^{\tau}U^{\tau}\ket{\psi_{0}}|^{2}=1\,.$
(41)

Thus, as claimed, if the cost function $C_{\rm fsVFF}$ vanishes the simulation
fidelity is perfect for all times.

##  Appendix B Noise Resilience of Cost

To make a connection with prior results on noise resilience, we first note
that $C_{\rm fsVFF}$, Eq. (5), can be rewritten as

$C_{\rm fsVFF}(U,V,\psi_{0})=\frac{1}{n_{\rm eig}}\sum_{k=1}^{n_{\rm
eig}}C_{\rm fsVFF}^{(k)}(U,V,\psi_{0})\,,$ (42)

with

$C_{\rm fsVFF}^{(k)}(U,V,\psi_{0})=1-{\rm
Tr}\left[|\mathbf{0}\rangle\langle\mathbf{0}|Y_{k}^{\tiny\psi_{0}}|\mathbf{0}\rangle\langle\mathbf{0}|\left(Y_{k}^{\tiny\psi_{0}}\right)^{\dagger}\right]\,.$
(43)

Here
$Y^{\psi_{0}}_{k}:=V_{\psi_{0}}^{\dagger}WD^{k}W^{\dagger}U^{k}V_{\psi_{0}}$
with $V_{\psi_{0}}\ket{\mathbf{0}}=\ket{\psi_{0}}$. Let $QC_{k}$ denote the
circuit used to evaluate the cost term $C_{\rm fsVFF}^{(k)}$. Let
$\tilde{C}_{\mbox{\tiny fsVFF}}^{(k)}$ denote the noisy version of the cost
term $C_{\rm fsVFF}^{(k)}$, that is the cost evaluated when the circuit
$QC_{j}$ is run in the presence of noise. Let $\mathbb{V}_{k}^{\rm opt}$ and
$\mathbb{\tilde{V}}_{k}^{\rm opt}$ denote the sets of unitaries that optimize
$C_{\rm fsVFF}^{(k)}$ and $\tilde{C}_{\mbox{\tiny fsVFF}}^{(k)}$ respectively.
Now, it was shown in Sharma  _et al._ (2020b) that the costs $C_{\rm
fsVFF}^{(k)}$ are resilient to measurement noise, gate noise, and Pauli
channel noise in the sense that for circuits experiencing such noise we have
$\mathbb{\tilde{V}}_{k}^{\rm opt}\subseteq\mathbb{V}_{k}^{\rm opt}$. This
means that any set of parameters that minimize the noisy cost
$\tilde{C}_{\mbox{\tiny fsVFF}}^{(k)}$ are guaranteed to also minimize the
true exact cost $C_{\rm fsVFF}^{(k)}$.

We will now argue that this implies that $C_{\rm fsVFF}$ is also noise
resilient. To do so, we first note that the costs $C_{\rm fsVFF}^{(k)}$ can be
minimized simultaneously by any unitary $V^{\rm opt}$ that matches the target
unitary $U$ up to a global phase $\phi$, i.e. such that $V^{\rm
opt}=\exp(-i\phi)U$. Therefore, assuming that the ansatz for $V$ is
sufficiently expressive that the costs $C_{\rm fsVFF}^{(k)}$ can be
simultaneously minimized, the total cost $C_{\rm fsVFF}$ is minimized by the
set of unitaries that minimize each of the $C_{\rm fsVFF}^{(k)}$ costs
simultaneously, that is the set $\cap_{k}\mathbb{V}_{k}^{\rm
opt}:=\mathbb{V}^{\rm opt}$. (Note, if the ansatz is not sufficiently
expressive then it might not be possible to simultaneously minimize each of
the terms and therefore the intersection might be empty). Now, given that
$\mathbb{\tilde{V}}_{k}^{\rm opt}\subseteq\mathbb{V}_{k}^{\rm opt}$, it
follows that $\mathbb{\tilde{V}}^{\rm
opt}:=\cap_{k}\mathbb{\tilde{V}}_{k}^{\rm opt}\subseteq\mathbb{V}^{\rm opt}$.
Thus, as claimed the total cost $C_{\rm fsVFF}$ is also noise resilient in the
sense that any set of parameters that minimize the noisy cost
$\tilde{C}_{\mbox{\tiny fsVFF}}$ also minimize the true exact cost $C_{\rm
fsVFF}$.

##  Appendix C Local cost with trainability guarantee

While $C_{\rm fsVFF}$ was motivated in Section III.1 as a natural choice of
cost function to learn the evolution induced by a target unitary on a fixed
initial state, it is expected to encounter what is known as a barren plateau
for large simulation sizes McClean  _et al._ (2018); Cerezo  _et al._ (2021).
(See Refs. Cerezo and Coles (2021); Arrasmith  _et al._ (2020a); Holmes  _et
al._ (2021a); Volkoff and Coles (2021); Sharma  _et al._ (2020c); Pesah  _et
al._ (2020); Uvarov and Biamonte (2020); Marrero  _et al._ (2020); Patti  _et
al._ (2020); Holmes  _et al._ (2021b)) for further details about the barren
plateau phenomenon.) Namely, the gradient of the cost vanishes exponentially
with $n$. As a result, for large systems the cost landscape is prohibitively
flat and therefore an exponential precision is required to discern a
minimization direction. This precludes successful training.

In this appendix we introduce a local cost function to surmount this
difficulty. To motivate our local cost, we first recall that $C_{\rm fsVFF}$,
Eq. (5), can be rewritten as

$C_{\rm fsVFF}=1-\frac{1}{n_{\rm eig}}\sum_{k=1}^{n_{\rm eig}}{\rm
Tr}\left[|\mathbf{0}\rangle\langle\mathbf{0}|Y_{k}^{\tiny\psi_{0}}|\mathbf{0}\rangle\langle\mathbf{0}|\left(Y_{k}^{\tiny\psi_{0}}\right)^{\dagger}\right]\,,$
(44)

where
$Y^{\psi_{0}}_{k}:=V_{\psi_{0}}^{\dagger}WD^{k}W^{\dagger}U^{k}V_{\psi_{0}}$
with $V_{\psi_{0}}\ket{\mathbf{0}}=\ket{\psi_{0}}$. Analogously, we now define
the local fixed state VFF cost as

$C_{\rm fsVFF}^{\mbox{\tiny Local}}:=\frac{1}{n}\sum_{j=1}^{n}C_{\rm
fsVFF}^{\mbox{\tiny Local,j}}\,,$ (45)

with

$C_{\mbox{\tiny fsVFF}}^{\mbox{\tiny Local,j}}=1-\frac{1}{n_{\rm
eig}}\sum_{k=1}^{n_{\rm eig}}{\rm Tr}\left[\left(|0\rangle\langle
0|_{j}\otimes\openone_{\bar{j}}\right)Y_{k}^{\tiny\psi_{0}}|\mathbf{0}\rangle\langle\mathbf{0}|\left(Y_{k}^{\tiny\psi_{0}}\right)^{\dagger}\right]\,,$
(46)

and where $\bar{j}$ denotes all qubits except the $j_{\rm th}$ qubit. Each of
the $C_{\rm fsVFF}^{\mbox{\tiny Local, j}}$ terms can be measured using the
same Loschmidt echo circuit as $C_{\rm fsVFF}$ but the final measurement is
performed on just the $j_{\rm th}$ qubit (rather than all qubits).

Following the proof techniques of Khatri  _et al._ (2019); Sharma  _et al._
(2020b), it is possible to show that

$C_{\rm fsVFF}^{\mbox{\tiny Local}}\leqslant C_{\rm fsVFF}\leqslant nC_{\rm
fsVFF}^{\mbox{\tiny Local}}\,.$ (47)

It therefore follows that $C_{\rm fsVFF}^{\mbox{\tiny Local}}=0$ if and only
if $C_{\rm fsVFF}=0$. Thus, given that $C_{\rm fsVFF}$ is faithful (as shown
in Appendix A), its local variant $C_{\rm fsVFF}^{\mbox{\tiny Local}}$ is also
faithful.

Crucially, as $C_{\rm fsVFF}^{\mbox{\tiny Local}}$ is local, i.e. requires
only local rather than global measurements, as long as the ansatz is not too
deep, the cost landscape will be sufficiently featured for effective training.
We therefore advocate using $C_{\rm fsVFF}^{\mbox{\tiny Local}}$ for
simulations of larger systems.

##  Appendix D Cost function gradient derivation

Here we derive the analytic expressions for the gradient of the cost function
$C_{\rm fsVFF}(U,V,\psi_{0})$ for gradient descent optimisation. To emphasize
the independence of each of the terms in the cost function and its dependence
on $\bm{\theta}$ and $\bm{\gamma}$ we write

$C_{\rm fsVFF}(U,V,\psi_{0})=\frac{1}{n_{\rm eig}}\sum_{k=1}^{n_{\rm
eig}}C_{\rm fsVFF}^{(k)}(\bm{\theta},\bm{\gamma})\,,$ (48)

where we have

$\displaystyle C_{\rm fsVFF}^{(k)}(\bm{\theta},\bm{\gamma})$
$\displaystyle:=1-|\bra{\psi_{0}}W(\bm{\theta})D(k\bm{\gamma})W(\bm{\theta})^{\dagger}U^{k}\ket{\psi_{0}}|^{2}$
(49) $\displaystyle=1-{\rm
Tr}[XW(\bm{\theta})D(k\bm{\gamma})W(\bm{\theta})^{\dagger}Y_{k}W(\bm{\theta})D(k\bm{\gamma})^{\dagger}W(\bm{\theta})^{\dagger}]\,,$
(50)

with $X=\ket{\psi_{0}}\bra{\psi_{0}}$ and
$Y_{k}=U^{k}\ket{\psi_{0}}\bra{\psi_{0}}(U^{\dagger})^{k}$.

#### Expressions.

With the ansatz in (1), the partial derivative of $C_{\rm
fsVFF}(U,V,\psi_{0})$ with respect to $\theta_{l}$, a parameter of the
eigenvector operator $W(\bm{\theta})$, is

$\begin{split}\frac{\partial C_{\rm
fsVFF}(U,V,\psi_{0})}{\partial\theta_{l}}=\frac{1}{2}\Big{(}&C_{\rm
fsVFF}(U,W_{l+}DW^{\dagger})\,-\,C_{\rm fsVFF}(U,W_{l-}DW^{\dagger})\\\ +\
&C_{\rm fsVFF}(U,WD(W_{l+})^{\dagger})\ -\ C_{\rm
fsVFF}(U,WD(W_{l-})^{\dagger})\Big{)}\ .\end{split}$ (51)

The operator $W_{l+}$ ($W_{l-}$) is generated from the original eigenvector
operator $W(\bm{\theta})$ by the addition of an extra $\frac{\pi}{2}$
($-\frac{\pi}{2}$) rotation about a given parameter’s rotation axis:

$W_{l\pm}:=W\left(\bm{\theta}_{l\pm}\right)\ \ \text{with}\ \
(\theta_{l\pm})_{i}:=\theta_{l}\pm\frac{\pi}{2}\delta_{i,l}\;.$ (52)

Similarly, the partial derivative with respect to $\gamma_{l}$, a parameter of
the diagonal operator $D(\bm{\gamma})$, is

$\begin{split}\frac{\partial C_{\rm
fsVFF}}{\partial\gamma_{l}}=\frac{1}{n_{\rm eig}}\sum_{k=1}^{n_{\rm
eig}}\frac{k}{2}\Big{(}&C_{\rm
fsVFF}^{(k)}\left(U,WD_{l+}W^{\dagger}\right)-C_{\rm
fsVFF}^{(k)}\left(U,WD_{l-}W^{\dagger}\right)\Big{)}\,,\end{split}$ (53)

where

$C_{\rm
fsVFF}^{(k)}:=1-|\bra{\psi_{0}}WD^{k}W^{\dagger}U^{k}\ket{\psi_{0}}|^{2}\,$
(54)

and

$D_{l\pm}:=D\left(\bm{\gamma}_{l\pm}\right)\ \ \text{with}\ \
(\gamma_{l\pm})_{i}:=\gamma_{l}\pm\frac{\pi}{2}\delta_{i,l}\;.$ (55)

#### Derivative with respect to $\gamma_{l}$.

Here we provide the derivation of the partial derivative of $C_{\rm
fsVFF}(U,V,\psi_{0})$ with respect to $\gamma_{l}$ in (53). Taking the partial
derivative of $C_{\rm fsVFF}^{(k)}(\bm{\theta},\bm{\gamma})$ with respect to
an angle $\gamma_{l}$ gives

$\displaystyle\frac{\partial C_{\rm fsVFF}^{(k)}}{\partial\gamma_{l}}=$
$\displaystyle-{\rm Tr}\left[X\left(W\frac{\partial
D(k\bm{\gamma})}{\partial\gamma_{l}}W^{\dagger}\right)\ket{Y_{k}}\\!\bra{Y_{k}}\left(WD(k\bm{\gamma})^{\dagger}W^{\dagger}\right)\right]$
(56) $\displaystyle-{\rm
Tr}\left[X\left(WD(k\bm{\gamma})W^{\dagger}\right)\ket{Y_{k}}\\!\bra{Y_{k}}\left(W\frac{\partial
D(k\bm{\gamma})^{\dagger}}{\partial\gamma_{l}}W^{\dagger}\right)\right]\,.$

The eigenvector operator, $D$, consists of products of Pauli rotations and can
be decomposed as

$D=D_{L}\exp\left(-\frac{ik\gamma_{l}\sigma_{l}}{2}\right)D_{R^{\prime}}\equiv
D_{L}D_{R}\,,$ (57)

where the operators $D_{L}$ and $D_{R^{\prime}}$ consist of all Pauli
rotations to the left and right of the $\sigma_{l}$ rotation respectively and
we have defined $D_{R}=\exp(-i\theta_{l}\sigma_{l}/2)D_{R^{\prime}}$ for
convenience. It follows that the differential of $W$ with respect to
$\theta_{l}$ takes the form

$\frac{\partial D}{\partial\gamma_{l}}=-\frac{1}{2}ikD_{L}\sigma_{l}D_{R}\;,$
(58)

which on substituting into Eq. (56) gives

$\displaystyle\frac{\partial C_{\rm fsVFF}^{(k)}}{\partial\gamma_{l}}$
$\displaystyle=\frac{ik}{2}\bigg{(}{\rm
Tr}\left[XWD_{L}\sigma_{l}D_{R}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD^{\dagger}W^{\dagger}\right]-{\rm
Tr}\left[XWDW^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{R}^{\dagger}\sigma_{l}D_{L}^{\dagger}W^{\dagger}\right]\bigg{)}$
(59) $\displaystyle=\frac{ik}{2}{\rm
Tr}\left[XWD_{L}\left(\sigma_{l}D_{R}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{R}^{\dagger}-D_{R}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{R}^{\dagger}\sigma_{l}\right)D_{L}^{\dagger}W^{\dagger}\right]$
$\displaystyle=\frac{ik}{2}{\rm
Tr}\left[XWD_{L}[\sigma_{l},\rho_{1}^{(k)}]D_{L}^{\dagger}W^{\dagger}\right]\,,$

where we have defined

$\displaystyle\rho_{1}^{(k)}=D_{R}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{R}^{\dagger}\,.$
(60)

Eq. (53) is now obtained directly from Eq. (59) via the following identity,
which holds for any state $\rho$,

$\displaystyle i[\sigma_{l},\rho]$ $\displaystyle=e^{i\sigma_{l}\pi/4}\rho
e^{-i\sigma_{l}\pi/4}-e^{-i\sigma_{l}\pi/4}\rho e^{i\sigma_{l}\pi/4}\,.$ (61)

Specifically we find that

$\displaystyle\frac{\partial C_{\rm
fsVFF}^{(k)}}{\partial\gamma_{l}}=\frac{k}{2}$ $\displaystyle{\rm
Tr}\left[X\otimes
WD_{L}(e^{i\sigma_{l}\pi/4}D_{R}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{R}^{\dagger}e^{-i\sigma_{l}\pi/4})D_{L}W^{\dagger}\right]$
(62) $\displaystyle-\frac{k}{2}$ $\displaystyle{\rm
Tr}\left[XWD_{L}(e^{-i\sigma_{l}\pi/4}D_{R}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{R}^{\dagger}e^{i\sigma_{l}\pi/4})D_{L}W^{\dagger}\right]$
$\displaystyle=\frac{k}{2}$ $\displaystyle{\rm Tr}\left[X\otimes
WD_{l-}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{l-}W^{\dagger}\right]$
$\displaystyle-\frac{k}{2}$ $\displaystyle{\rm Tr}\left[X\otimes
WD_{l+}W^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD_{l+}W^{\dagger}\right]\,$
$\displaystyle=\frac{k}{2}\Big{(}$ $\displaystyle C_{\rm
fsVFF}^{(k)}\left(U,WD_{l+}W^{\dagger}\right)-C_{\rm
fsVFF}^{(k)}\left(U,WD_{l-}W^{\dagger}\right)\Big{)}\,.$

Thus we are left with

$\begin{split}\frac{\partial C_{\rm
fsVFF}}{\partial\gamma_{l}}=\frac{1}{n_{\rm eig}}\sum_{k=1}^{n_{\rm
eig}}\frac{\partial C_{\rm fsVFF}^{(k)}}{\partial\gamma_{l}}=\frac{1}{n_{\rm
eig}}\sum_{k=1}^{n_{\rm eig}}\frac{k}{2}\Big{(}&C_{\rm
fsVFF}^{(k)}\left(U,WD_{l+}W^{\dagger}\right)-C_{\rm
fsVFF}^{(k)}\left(U,WD_{l-}W^{\dagger}\right)\Big{)}\,.\end{split}$ (63)

#### Derivative with respect to $\theta_{l}$.

Here we provide the derivation of the partial derivative of
$C_{\mbox{fsVFF}}(U,V,\psi_{0})$ with respect to $\theta_{l}$ in (51). Taking
the partial derivative of $C_{\rm fsVFF}^{(k)}(\bm{\theta})$ with respect to
an angle $\theta_{l}$ gives

$\displaystyle\frac{\partial C_{\rm fsVFF}^{(k)}}{\partial\theta_{l}}=$
$\displaystyle-{\rm Tr}\left[X\left(\frac{\partial
W}{\partial\theta_{l}}DW^{\dagger}\right)\ket{Y_{k}}\\!\bra{Y_{k}}\left(WD^{\dagger}W^{\dagger}\right)\right]-{\rm
Tr}\left[X\left(WDW^{\dagger}\right)\ket{Y_{k}}\\!\bra{Y_{k}}\left(WD^{\dagger}\frac{\partial
W^{\dagger}}{\partial\theta_{l}}\right)\right]$ (64) $\displaystyle-{\rm
Tr}\left[X\left(WD\frac{\partial
W^{\dagger}}{\partial\theta_{l}}\right)\ket{Y_{k}}\\!\bra{Y_{k}}\left(WD^{\dagger}W^{\dagger}\right)\right]-{\rm
Tr}\left[X\left(WDW^{\dagger}\right)\ket{Y_{k}}\\!\bra{Y_{k}}\left(\frac{\partial
W}{\partial\theta_{l}}D^{\dagger}W^{\dagger}\right)\right].$

The eigenvector operator, $W$, consists of products of Pauli rotations and can
be decomposed as

$W=W_{L}\exp\left(-\frac{i\theta_{l}\sigma_{l}}{2}\right)W_{R^{\prime}}\equiv
W_{L}W_{R}\,,$ (65)

where the operators $W_{L}$ and $W_{R^{\prime}}$ consist of all Pauli
rotations to the left and right of the $\sigma_{l}$ rotation respectively and
we have defined $W_{R}=\exp(-i\theta_{l}\sigma_{l}/2)W_{R^{\prime}}$ for
convenience. It follows that the differential of $W$ with respect to
$\theta_{l}$ takes the form

$\frac{\partial W}{\partial\theta_{l}}=-\frac{1}{2}iW_{L}\sigma_{l}W_{R}\;,$
(66)

which on substituting into Eq. (64) gives

$\displaystyle\frac{\partial C_{\rm
fsVFF}^{(k)}}{\partial\theta_{l}}=\frac{i}{2}\bigg{(}{\rm
Tr}\left[XW_{L}[\sigma_{l},\rho_{2}^{(k)}]W_{L}^{\dagger}\right]-{\rm
Tr}\left[XWDW_{R}^{\dagger}[\sigma_{l},\rho_{3}^{(k)}]W_{R}D^{\dagger}W^{\dagger}\right]\bigg{)}\;,$
(67)

where we have defined

$\displaystyle\rho_{2}^{(k)}=W_{R}DW^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}WD^{\dagger}W_{R}^{\dagger}\
\ \ \text{and}\ \ \
\rho_{3}^{(k)}=W_{l}^{\dagger}\ket{Y_{k}}\\!\bra{Y_{k}}W_{L}\ .$ (68)

Eq. (51) is now obtained directly from Eq. (67) via Eq. (61).

##  Appendix E Fast forwarding error

Let us start by writing the state $\ket{\psi_{0}}$ as

$\ket{\psi_{0}}=\left(\begin{array}[]{@{}c@{}}\mathbf{v_{\psi_{0}}}\\\
\hline\cr\mathbf{0}\end{array}\right)$

to emphasise that $\ket{\psi_{0}}$ spans a subspace of the total Hilbert
space. We then note that in the limit in which leakage can be disregarded the
evolution unitary can be written in terms of its block decomposition

$U=\left(\begin{array}[]{@{}c|c@{}}U_{\parallel}&\mbox{\Large 0}\\\
\hline\cr\mbox{\Large 0}&U_{\perp}\\\ \end{array}\right).$

where $U_{\parallel}$ acts on the subspace spanned by $\psi_{0}$ and
$U_{\perp}$ acts on the rest of the Hilbert space. Let us also write the
learnt unitary in a block decomposition form as

$V=\left(\begin{array}[]{@{}c|c@{}}V_{\parallel}&\epsilon_{a}\\\
\hline\cr\epsilon_{b}&V_{\perp}\\\ \end{array}\right).$

In general, the off-diagonal blocks are expected to be close to the null
matrices, therefore we can take a perturbative approach to calculating
$V^{N}$. Expanding $V^{N}$ to first order in $\epsilon_{a}$ and $\epsilon_{b}$
gives

$V^{N}=\left(\begin{array}[]{@{}c|c@{}}V_{\parallel}&\epsilon_{a}\\\
\hline\cr\epsilon_{b}&V_{\perp}\\\
\end{array}\right)^{N}=\left(\begin{array}[]{@{}c|c@{}}V_{\parallel}^{N}&V_{\epsilon_{a}}(N)\\\
\hline\cr V_{\epsilon_{b}}(N)&V_{\perp}^{N}\\\ \end{array}\right)\,,$

where

$\displaystyle
V_{\epsilon_{a}}(N):=V_{\parallel}^{N-2}\epsilon_{a}V_{\parallel}+\sum_{k=1}^{N-1}V_{\parallel}^{N-1-k}\epsilon_{a}V_{\perp}^{k}$
(69) $\displaystyle
V_{\epsilon_{b}}(N):=V_{\perp}^{N-2}\epsilon_{b}V_{\perp}+\sum_{k=1}^{N-1}V_{\perp}^{N-1-k}\epsilon_{b}V_{\parallel}^{k}\,.$

Therefore we can now write

$V^{N}\ket{\psi_{0}}=\left(\begin{array}[]{@{}c|c@{}}V_{\parallel}^{N}&V_{\epsilon_{a}}(N)\\\
\hline\cr V_{\epsilon_{b}}(N)&V_{\perp}^{N}\\\
\end{array}\right)\left(\begin{array}[]{@{}c@{}}\mathbf{v_{\psi_{0}}}\\\
\hline\cr\mathbf{0}\end{array}\right)=\left(\begin{array}[]{@{}c@{}}V_{\parallel}^{N}\mathbf{v_{\psi_{0}}}\\\
\hline\cr V_{\epsilon_{b}}(N)\mathbf{v_{\psi_{0}}}\end{array}\right).$

We are interested in the fast forwarded simulation fidelity

$F_{N}=|\bra{\psi_{0}}{V^{\dagger}}^{N}U^{N}\ket{\psi_{0}}|^{2}\,,$ (70)

which to first order in $\epsilon_{a}$ and $\epsilon_{b}$ evaluates to

$F_{N}=|\bra{\psi_{0}}{V_{\parallel}^{\dagger}}^{N}U_{\parallel}^{N}\ket{\psi_{0}}|^{2}\,.$
(71)

The fidelity between any two quantum states $\ket{\psi_{a}}$ and
$\ket{\psi_{b}}$ can be related to the trace norm distance between the exact
and simulated fast forwarded states via the relation

$1-|\braket{\psi_{a}}{\psi_{b}}|^{2}=\left(\frac{1}{2}||\psi_{a}-\psi_{b}||_{1}\right)^{2}\,,$
(72)

where we use the shorthand $\psi_{a}=\ket{\psi_{a}}\bra{\psi_{a}}$ and
$\psi_{b}=\ket{\psi_{b}}\bra{\psi_{b}}$. Therefore we can write

$1-F_{N}=\left\|U_{\parallel}^{N}\psi_{0}{U_{\parallel}^{\dagger}}^{N}-V_{\parallel}^{N}\psi_{0}{V_{\parallel}^{\dagger}}^{N}\right\|_{1}\,.$
(73)

We further note that we can write

$\left\|U_{\parallel}^{N}\psi_{0}{U_{\parallel}^{\dagger}}^{N}-V_{\parallel}^{N}\psi_{0}{V_{\parallel}^{\dagger}}^{N}\right\|_{1}=\left\|(U_{\parallel}^{N}-V_{\parallel}^{N})\psi_{0}{U_{\parallel}^{\dagger}}^{N}+V_{\parallel}^{N}\psi_{0}\left({U_{\parallel}^{\dagger}}^{N}-{V_{\parallel}^{\dagger}}^{N}\right)\right\|_{1}\,,$
(74)

and therefore on applying Holder’s inequality,
$||XY||_{1}\leqslant||X||_{\infty}||Y||_{1}$, we have that

$\left\|U_{\parallel}^{N}\psi_{0}{U_{\parallel}^{\dagger}}^{N}-V_{\parallel}^{N}\psi_{0}{V_{\parallel}^{\dagger}}^{N}\right\|_{1}\leqslant\left\|U_{\parallel}^{N}-V_{\parallel}^{N}\right\|_{\infty}\left(\left\|\psi_{0}{U_{\parallel}^{\dagger}}^{N}\right\|_{1}+\left\|\psi_{0}{V_{\parallel}}^{N}\right\|_{1}\right)\,.$
(75)

Now from Lemma 1 of Cirstoiu  _et al._ (2020) we have that
$||U_{\parallel}^{N}-V_{\parallel}^{N}||_{\infty}\leqslant
N||U_{\parallel}-V_{\parallel}||_{\infty}$ and from the unitary invariance of
the Schatten norms we have that
$||\psi_{0}{U_{\parallel}^{\dagger}}^{N}||_{1}=||\psi_{0}{V_{\parallel}^{\dagger}}^{N}||_{1}=||\psi_{0}||_{1}=1$,
therefore we are left with

$\left\|U_{\parallel}^{N}\psi_{0}{U_{\parallel}^{\dagger}}^{N}-V_{\parallel}^{N}\psi_{0}{V_{\parallel}^{\dagger}}^{N}\right\|_{1}\leqslant
2N\left\|U_{\parallel}-V_{\parallel}\right\|_{\infty}\,.$ (76)

We conclude that the final simulation fidelity is bounded as

$1-F_{N}\leqslant
N^{2}\left(\left\|U_{\parallel}-V_{\parallel}\right\|_{\infty}\right)^{2}\,.$
(77)

Thus, the simulation fidelity scales sub-quadratically with the product of the
number of simulation steps and the infinity norm distance between the learnt
and target unitaries.

Figure 12: Quadratic Scaling of Infidelity: The 3-qubit diagonalization used
to perform the Eigenvalue Estimation from Figure 11 is fast-forwarded with its
infidelity, 1 - $F$, plotted against the number of timesteps, $N$. The close
alignment of infidelity against the dotted line with form $c\times N^{2}$
reproduces the scaling of Eq. 77.

##  Appendix F Symmetry preserving adaptive ansatz design

Here we provide further details on the adaptive ansatz used to optimise the
noisy simulations and for the batch sampling simulations. This method uses a
combination of discrete updates to the structure of the circuit and the
training of continuous rotation angles. To facilitate discrete structure
updates, at the start of the optimization’s inner loop a set of gates
initialised to the identity gate are randomly inserted into the ansatz, and
when the cost function has plateaued after optimizing the continuous
parameters, each gate is tested to see how much the cost function is increased
by removing that gate, with a change in cost $\Delta C_{1}$. In a manner
analogous to that used in simulated annealing, if
$\text{exp}(-\beta_{1}\frac{\Delta C_{1}}{C})$ < rand(0, 1) then the gate is
considered not to be significantly contributing to the diagonalization and is
deleted from the ansatz. To avoid getting trapped in a local minima, after the
cost has plateaued at the end of the inner loop the structure update is
accepted if 1-$\text{exp}(-\beta_{2}\frac{\Delta C_{2}}{C_{\text{best}}})$ <
rand(0, 1), where $\Delta C_{2}=C_{\text{plateau}}-C_{\text{best}}$. A
subfunction to compile together consecutive identical gates was used to
compress the circuit depth.

The Heisenberg class of Hamiltonians are known to observe the particle number-
conserving symmetry, so a dictionary of gates $\mathcal{D}$ that all maintain
this symmetry was chosen to restrict the size of the subspace the ansatz
explored. $\mathcal{D}=\\{R_{z},R_{zz},G\\}$, where $G$ is an entangling gate
conserving the total amplitude within the $\\{\ket{01},\ket{10}\\}$ subspace,
described in Ref. Anselmetti  _et al._ (2021). Linear nearest-neighbour
connectivity was enforced during the noisy simulations to avoid the insertion
of unwanted SWAP gates.

In our noiseless batched simulations, the same gate dictionary $\mathcal{D}$
was used as in the noisy simulations, although here all-to-all connectivity
was allowed. In a noisy setting, the accumulation of errors due to circuit
depth automatically adds a constraint on the circuit depth, however in an
unconstrained noise-free setting the insertion of gates may cause the ansatz
to arbitrarily grow in length as the optimisation progresses. To prevent this,
a regularised cost was used during the structure updates, with a term added to
the fsVFF cost to penalise deeper circuits,
$C_{\text{reg}}=C_{\text{fsVFF}}*(1+\frac{N_{gates}}{\lambda})$. The added 1
in the multiplicative factor serves to remove the effect of the regularisation
term for small circuits, avoiding early deletion of gates before
$C_{\text{fsVFF}}$ has begun to decrease.