# Variational Quantum Time Evolution without the Quantum Geometric Tensor

Julien Gacon  IBM Quantum, IBM Research Europe – Zurich, CH-8803 Rüschlikon,
Switzerland  Institute of Physics, École Polytechnique Fédérale de Lausanne
(EPFL), CH-1015 Lausanne, Switzerland  Jannes Nys  Institute of Physics, École
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Riccardo Rossi  Sorbonne Université, CNRS, Laboratoire de Physique Théorique
de la Matière Condensée, LPTMC, F-75005 Paris, France  Stefan Woerner  IBM
Quantum, IBM Research Europe – Zurich, CH-8803 Rüschlikon, Switzerland
Giuseppe Carleo  Institute of Physics, École Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland

###### Abstract

The real- and imaginary-time evolution of quantum states are powerful tools in
physics, chemistry, and beyond, to investigate quantum dynamics, prepare
ground states or calculate thermodynamic observables. On near-term devices,
variational quantum time evolution is a promising candidate for these tasks,
as the required circuit model can be tailored to trade off available device
capabilities and approximation accuracy. However, even if the circuits can be
reliably executed, variational quantum time evolution algorithms quickly
become infeasible for relevant system sizes due to the calculation of the
Quantum Geometric Tensor (QGT). In this work, we propose a solution to this
scaling problem by leveraging a dual formulation that circumvents the explicit
evaluation of the QGT. We demonstrate our algorithm for the time evolution of
the Heisenberg Hamiltonian and show that it accurately reproduces the system
dynamics at a fraction of the cost of standard variational quantum time
evolution algorithms. As an application of quantum imaginary-time evolution,
we calculate a thermodynamic observable, the energy per site, of the
Heisenberg model.

††preprint: APS/123-QED

##  I Introduction

Quantum time evolution is a central task in physics. Real-time evolution
provides detailed insight into properties of quantum mechanical systems, such
as phase transitions [1, 2, 3] or thermalization [4, 5]. Imaginary-time
evolution is an important tool that enables the preparation of ground states
or thermal states [6, 7, 8]. These can, in turn, be used for the calculation
of thermodynamic observables [8, 9]. In particular, combining real- and
imaginary-time evolution would allow the direct calculation of dynamical
correlation functions at thermal equilibrium.

The range of applications of imaginary-time evolution extends beyond the field
of physics. Ground-state preparation with imaginary-time evolution for gapped,
non-degenerate Hamiltonians is guaranteed to converge in the generic case of
non-zero overlap between the ground state and the initial trial state. This
makes it a promising candidate in settings where a good initial state can be
constructed, e.g. in chemistry applications [10] or in classical optimization
problems [11]. In quantum machine learning, the preparation of Gibbs states
with imaginary-time evolution is a subroutine for quantum Boltzmann machines,
which can, for example, be used in distribution learning or classification
[12].

Figure 1: Estimated runtimes of variational imaginary-time evolution (VarQITE)
and our proposed dual method (DualQITE) as a function of the number of
parameters $d$ of the variational model, for an exemplary Heisenberg model and
200 timesteps. See Appendix A for more details.

Since performing quantum time evolution generally requires representing the
exponentially large wave function of a quantum system, quantum computers are a
promising platform for developing efficient algorithms [13]. In fact, in 1996,
the Trotter algorithm for real-time evolution was among the first proposed use
cases for a quantum computer [14]. However, the complexity of the quantum
circuits required for the Trotter algorithm depends on the Hamiltonian, and
the circuit depth scales with the simulation time and accuracy [15, 16]. This
renders the algorithm currently unsuitable for general time evolution on near-
term devices, which are characterized by limited qubit connectivity and
coherence times. The imaginary-time counterpart of Trotter suffers from the
same restriction [8].

Variational algorithms for quantum time evolution, on the other hand, allow to
choose a parameterized circuit as an ansatz to approximate the wave function,
that operates within the device’s capabilities. Using a variational principle,
variational quantum time evolution (VarQTE) maps the quantum state evolution
to the evolution of parameters in the model [17], both for _real_ -time
evolution (VarQRTE) and _imaginary_ -time evolution (VarQITE). The parameter
update rules depend on the evaluation of the Quantum Geometric Tensor (QGT)
and gradients of the current energy and state. For an ansatz with $d$
variational parameters, the number of circuits required to evaluate the QGT
and gradients scale as $\mathcal{O}(d^{2})$ and $\mathcal{O}(d)$,
respectively. While this does not pose a problem for small systems, the
evaluation of the QGT quickly becomes a bottleneck once the system size, and
therefore the number of variational parameters, increases.

Figure 1 shows a runtime estimation for VarQITE, assuming current
superconducting processor specifications, see Appendix A for details on the
derivation. For a few parameters the runtime is of the order of hours, but
already for only 200 parameters the computation time around 1 week, which
renders this algorithm currently impractical. With recent advances in
processor sizes exceeding 100 qubits, such as the IBM Quantum Eagle or Osprey
devices [18, 19], improving the resource requirements of quantum algorithms
becomes crucial for finding practically relevant applications of quantum
computers.

Recently, focus has shifted to optimization-based algorithms, which implement
partial steps or approximations of the full Suzuki-Trotter step [20, 21, 22,
23, 24]. In the case of real-time evolution, for example, the projected
variational quantum dynamics (p-VQD) algorithm [20] provides a scalable
alternative to VarQRTE on near term-devices, if a single Trotter step can be
efficiently implemented. However, the required quantum circuit gates in p-VQD
reflect the couplings of the Hamiltonian. This means that, for Hamiltonians
with long-distance interactions or numerous Pauli terms (e.g. in molecular
dynamics), even a single step could involve global connections or deep
circuits hindering the execution on near-term devices. Furthermore, the p-VQD
algorithm is not directly applicable to imaginary-time evolution.

Other approaches concerned with real-time evolution are Variational Fast
Forwarding (VFF) methods [25, 26, 27, 28] and classical pre-processing
approaches [29, 30]. VFF methods rely on diagonalizing the Hamiltonian or the
Trotterized time evolution operator with a variational ansatz. However,
finding the diagonalizing unitary remains challenging in practice, which
limits demonstrations to very few qubits. Classical pre-processing techniques,
on the other hand, impose additional restrictions on the simulated system,
such as translational invariance [30] or Hamiltonians with low entanglement
[29]. Within such systems, these techniques scale to large systems, but they
do not allow for general quantum time evolution.

Another line of work directly focuses on the preparation of thermal states by
minimizing the free energy of a variational ansatz [31]. This approach,
however, also does not implement general quantum time evolution.

In this paper, we propose a novel variational algorithm for quantum time
evolution based on a _dual_ optimization problem, which allows to replace the
QGT by evaluating the overlap of the variational ansatz for different
parameter values. This formulation applies equally real- and imaginary-time
evolution and does not require additional qubits or connections than already
present in the ansatz. We show that this new algorithm requires significantly
fewer measurements and thereby drastically reduces the expected runtime
compared to VarQTE. This is summarized in Fig. 1, where, under the same
assumptions, our proposed method can reduce the expected runtimes from several
weeks for VarQTE to only a few days. Following the naming conventions of
VarQTE, we name the algorithm DualQTE with specifiers DualQITE for imaginary-
time evolution and DualQRTE for real-time evolution.

The remainder of this paper is structured as follows. In Sec. II, we recap
VarQTE based on variational principles, derive the proposed dual formulation,
and discuss how to implement it on a quantum computer. Then, in Sec. III, we
demonstrate our proposed algorithm for the imaginary-time evolution of the
Heisenberg model and investigate the resource requirements. As a practical
application, we use the quantum minimally entangled typical thermal states
method (QMETTS) to calculate thermodynamic observables. Sec. IV demonstrates
the dual formulation for real-time evolution, including the calculation of
variational error bounds. Finally, Sec. V concludes the paper and gives an
outlook on possible applications and further research directions.

##  II Dual formulation of variational time evolution

For a time-independent Hamiltonian $H$ acting on $n$ qubits, an initial
quantum state $\ket{\Psi_{0}}$ and an evolution time $t$, the real-time
evolved quantum state is

$\ket{\Psi(t)}=e^{-itH}\ket{\Psi_{0}}.$ (1)

For an imaginary-time evolution, the time evolution operator is non-unitary,
and the normalized state reads

$\ket{\Psi(t)}=\frac{1}{\sqrt{\braket{\Psi_{0}}{e^{-2tH}}{\Psi_{0}}}}e^{-tH}\ket{\Psi_{0}}.$
(2)

Variational quantum time evolution maps the evolution of the quantum state
$\ket{\Psi(t)}$ to the evolution of parameters
$\bm{\theta}(t)\in\mathbb{R}^{d}$ of a parameterized quantum state
$\ket{\phi(\bm{\theta}(t))}$. The parameters’ dynamics can be derived with
variational principles such as the Dirac-Frenkel, McLachlan, or time-dependent
variational principle [17]. In McLachlan’s formulation, the derivatives of the
parameters are determined by the linear system of equations

$g(\bm{\theta}(t))\,\dot{\bm{\theta}}(t)=\bm{b}(\bm{\theta}(t)),$ (3)

where the matrix $g=\mathrm{Re}(G)\in\mathbb{R}^{d\times d}$ is the real part
of the QGT, and we call $b\in\mathbb{R}^{d}$ the evolution gradient.

The QGT is defined as

$G_{ij}(\bm{\theta})=\braket{\partial_{i}\phi(\bm{\theta})}{\partial_{j}\phi(\bm{\theta})}-\braket{\partial_{i}\phi(\bm{\theta})}{\phi(\bm{\theta})}\braket{\phi(\bm{\theta})}{\partial_{j}\phi(\bm{\theta})},$
(4)

where we use the notation $\partial_{i}:=\partial/(\partial\theta_{i})$ and do
not explicitly state the time dependence of the parameters. The evolution
gradient for VarQRTE is given by the expression

$b^{\text{R}}_{i}(\bm{\theta})=\mathrm{Im}\big{(}\braket{\partial_{i}\phi(\bm{\theta})}{H}{\phi(\bm{\theta})}-\braket{\partial_{i}\phi(\bm{\theta})}{\phi(\bm{\theta})}E(\bm{\theta})\big{)},$
(5)

whereas a VarQITE evolution yields the following

$b^{\text{I}}_{i}(\bm{\theta})=-\mathrm{Re}\big{(}\braket{\partial_{i}\phi(\bm{\theta})}{H}{\phi(\bm{\theta})}\big{)}=-\frac{\partial_{i}E(\bm{\theta})}{2},$
(6)

with the energy
$E(\bm{\theta})=\braket{\phi(\bm{\theta})}{H}{\phi(\bm{\theta})}$. From
hereon, we present general equations that apply to both real and imaginary-
time evolution; thus, unless specified, we simply use $b$ without a specific
superscript.

Note that these equations are introduced for a time-independent Hamiltonian,
but they can also be applied to the time-dependent case $H=H(t)$.

###  II.1 Dual formulation

Instead of solving the linear system defined in Eq. (3), we propose to solve
the dual formulation of the problem [32, 33] given by

$\dot{\bm{\theta}}=\underset{\dot{\bm{\theta}}}{\mathrm{argmin}}\frac{\dot{\bm{\theta}}^{T}g(\bm{\theta})\dot{\bm{\theta}}}{2}-\dot{\bm{\theta}}^{T}\bm{b}(\bm{\theta}).$
(7)

The term
$\|\dot{\bm{\theta}}\|^{2}_{g(\bm{\theta})}=\dot{\bm{\theta}}^{T}g(\bm{\theta})\dot{\bm{\theta}}$
is the squared norm of the parameter derivative in the metric of the QGT. This
quantity describes the magnitude of the derivative from an information
geometric point of view and is derived from the Fubini-Study metric. For
infinitesimal displacements $\bm{\delta\theta}$, we have

$\displaystyle||\bm{\delta\theta}||^{2}_{g(\bm{\theta})}$
$\displaystyle=\bm{\delta\theta}^{T}g(\bm{\theta})\bm{\delta\theta}$ (8)
$\displaystyle=1-|\braket{\phi(\bm{\theta})}{\phi(\bm{\theta}+\bm{\delta\theta})}|^{2}+\mathcal{O}(\|\bm{\delta\theta}\|_{2}^{3}),$

where $\|\cdot\|_{2}$ is the $\ell_{2}$ norm [33]. By writing
$\dot{\bm{\theta}}=\bm{\delta\theta}/\delta\tau$, for some time perturbation
$\delta\tau>0$, we can now reformulate the optimization in terms of the
fidelity
$F(\bm{\theta},\bm{\theta^{\prime}})=|\braket{\phi(\bm{\theta})}{\phi(\bm{\theta}^{\prime})}|^{2}$
as

$\displaystyle\bm{\delta\theta}$
$\displaystyle\approx\underset{\bm{\delta\theta}}{\mathrm{argmin}}\frac{1-F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})}{2(\delta\tau)^{2}}-\frac{\bm{\delta\theta}^{T}\bm{b}(\bm{\theta})}{\delta\tau}$
(9)
$\displaystyle=\underset{\bm{\delta\theta}}{\mathrm{argmin}}\frac{\mathcal{L}(\bm{\delta\theta})}{(\delta\tau)^{2}},$

where we directly optimize for the parameter update $\bm{\delta\theta}$ and we
introduced the loss function

$\mathcal{L}(\bm{\delta\theta})=\frac{1-F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})}{2}-\delta\tau\cdot\bm{\delta\theta}^{T}\bm{b}(\bm{\theta}).$
(10)

In practice, the optimization problem can be solved without the factor
$(\delta\tau)^{-2}$, which decouples the shape of the locally quadratic
infidelity term from the time perturbation and improves the numerical
stability of the optimization.

Note that this dual formulation can alternatively be obtained from the
derivation of quantum natural gradients [33, 31], which is detailed in
Appendix B. For an intuitive understanding of the relationship of the
infidelity and QGT the effect of approximating
$\|\bm{\delta\theta}\|_{g(\bm{\theta})}\approx
1-F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})$ in an illustrative example is
demonstrated in Appendix C.

Instead of computing the QGT at each timestep, which requires
$\mathcal{O}(d^{2})$ circuit evaluations, we now have to solve an optimization
problem where the loss function requires only one fidelity evaluation. The
required resources of DualQTE per timestep are therefore $\mathcal{O}(d)$ for
the computation of the evolution gradient $b$, times the number of iterations
in the optimization. Thus, we improve upon the direct QGT approach if the
number of iterations scales better than $\mathcal{O}(d)$, which, as we show in
the following sections, is the case for the examples we investigate in this
work.

###  II.2 Evaluating the loss function

The evaluation of the loss function $\mathcal{L}$, defined in Eq. (10),
requires the calculation of the evolution gradient $b$ and the fidelity of the
ansatz $\ket{\phi(\bm{\theta})}$ for two different parameter sets. For
imaginary-time evolution, the evolution gradient can, for example, be
evaluated with analytic gradient rules, such as the parameter-shift rule or a
linear combination of unitaries (LCU), or with finite difference methods [34].
In the case of real-time evolution, however, we are restricted to an LCU
approach, as this is the only method that allows the calculation of the
imaginary part of gradients [12].

The fidelity $F$ can, for example, be estimated using the swap test [35] and
its variants [36], where the states are prepared in separate qubit registers
followed by entangling gates across these registers, or with the Hadamard
test, which adds only a single auxiliary qubit, but requires controlling the
state-preparing unitary [37]. A more near-term-friendly option is the compute-
uncompute method [38], which does not introduce additional global operations.
If the states are given by $\ket{\phi(\bm{\theta})}=U(\bm{\theta})\ket{0}$ for
a parameterized unitary $U$ and two different parameter values $\bm{\theta}$
and $\bm{\theta}^{\prime}$, the fidelity can be calculated by preparing
$U^{\dagger}(\bm{\theta})U(\bm{\theta}^{\prime})\ket{0}$ and measuring the
probability of obtaining $\ket{0}$ on all qubits.

If the state $\ket{\phi(\bm{\theta})}$ has $n$ qubits and the preparing
unitary $U$ has depth $m$, the swap test variants require a circuit width of
$2n$ with depth of $m+\mathcal{O}(1)$, whereas the evaluated circuits for the
compute-uncompute method are of only width $n$, but of depth $2m$. The
Hadamard test for fidelities between the same circuit with different
parameters can be evaluated by controlling the parameterized gates, resulting
in depth of $m$ and depth of $n+1$, plus the overhead of controlling the
gates. For sparse device connectivities, this can be a challenge. To avoid
increasing the circuit complexity, the overlap can also be estimated via
randomized measurements of two independent state preparations [39]. However,
this technique requires an exponential number of measurements.

Evaluating the QGT for VarQTE, however, suffers from similar issues. The QGT
can be evaluated as the Hessian of the infidelity [40] using a parameter-shift
or finite difference technique, which comes with the restrictions for fidelity
evaluations as described above. Alternatively, Eq. (4) can be directly
computed with an LCU approach, which adds two auxiliary qubits and two
entangling gates [12]. This method is less demanding than e.g. a Hadamard
test, but still comes with additional connectivity requirements. In practice,
for both VarQTE and DualQTE a suitable combination of parameterized quantum
state $\ket{\phi(\bm{\theta})}$ and gradient method must be selected, such
that the resulting circuits can be executed reliably.

Depending on the topology and coherence times of the available hardware and
the structure and size of the unitary, either method for gradient and fidelity
calculations can be advantageous. In this work, we focus on near-term friendly
methods and use the parameter-shift rule for gradients (if possible) and the
compute-uncompute method for the fidelity, as these do not require additional
gate connections or an exponential number of measurements. Note that, for
systems with a large number of qubits, this method might become unsuitable as
it measures the global zero projector. Then, approaches using only local
measurements, such as the Hadamard test, could be the better choice.

###  II.3 Solving for the update step

The infidelity-based loss function $\mathcal{L}$ is locally convex around
$\bm{\delta\theta}=\bm{0}$, as its Hessian at this point is
$\bm{\nabla}\bm{\nabla}^{T}\mathcal{L}(0)=g/2$, and $g$ is positive semi-
definite. To leverage this property, we use gradient descent as a local
optimization routine, which also allows the use of analytic gradient formulas
that have proven more stable in presence of shot noise. The gradient of
$\mathcal{L}$ with respect to the parameter update $\bm{\delta\theta}$ is

$\bm{\nabla}_{\bm{\delta\theta}}\mathcal{L}(\bm{\delta\theta})=-\frac{\bm{\nabla}_{\bm{\delta\theta}}F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})}{2}-\delta\tau\cdot\bm{b}(\bm{\theta}).$

The gradient of the fidelity can be evaluated with a parameter-shift rule

$\frac{\partial
F}{\partial(\delta\theta)_{i}}=\frac{F(\bm{\theta},\bm{\theta}+\bm{\delta\theta}+\bm{e}_{i}s)-F(\bm{\theta},\bm{\theta}+\bm{\delta\theta}-\bm{e}_{i}s)}{2\sin(s)},$

where $\bm{e}_{i}$ is the $i$-th unit vector and $s$ is the parameter shift,
which can be chosen as, e.g., $\pi/2$ for single-qubit Pauli rotations [34].

At each timestep, the gradient descent update for the update step
$\bm{\delta\theta}$ is

$\bm{\delta\theta}^{(k+1)}=\bm{\delta\theta}^{(k)}-\eta_{k}\bm{\nabla}\mathcal{L}\left(\bm{\delta\theta}^{(k)}\right),$

where $\eta_{k}>0$ is the learning rate at step $k$. This iteration is
continued until a maximum number of iterations or a convergence criterion is
met. An example of the latter is a minimum tolerance in the change of the cost
function or the norm of the gradient.

An intuitive choice for the initial guess $\bm{\delta\theta}^{(0)}$ is the
zero vector, which corresponds to no change in the parameters. However, a more
efficient choice can be to introduce momentum by warm starting the
optimization with the update step from the previous timestep. This heuristic
is motivated by the idea that, especially for small timesteps, we do not
expect the parameter derivatives $\dot{\bm{\theta}}$ to change significantly.

Methods that approximate the gradient, such as finite difference or SPSA [41],
may face challenges in the optimization. For small timesteps, the fidelity is
close to $1$ and the noise in the readout, e.g. from finite sampling
statistics or device noise, can easily mask changes in the cost function.
Parameter-shift gradients suffer less from this problem, as they allow to
evaluate the cost function over larger perturbations, and do not amplify the
noise by dividing by a small constant.

The ideal choice of the time perturbation $\delta\tau$ is a trade-off: the
error in approximating the QGT scales with $(\delta\tau)^{3}$, but a smaller
perturbation amplifies any measurement noise in the loss function as the
update step is obtained as $\bm{\delta\theta}/\delta\tau$. Appendix C displays
this trade-off for an illustrative example.

###  II.4 Trainability

Recently, there has been a lot of research showing that, in certain settings,
the loss function gradients of variational algorithms decay to zero
exponentially and cannot be evaluated efficiently, as they would require an
exponential number of measurements. These so-called barren plateaus can be
induced, for example, if the loss function requires measuring a global
observable [42], if the quantum circuit preparing the parameterized state is
too deep or generates too much entanglement [43, 42, 44], or if the
measurements are too noisy [45].

Since variational quantum dynamics is driven by the evolution gradient defined
Eqs. (5) and (6), it can be affected by a barren plateau and fail to track the
true evolution of the quantum state. However, it is important to note that the
gradients only vanish on average for a random initialization, whereas in time-
evolution the initial quantum state is typically specifically chosen.
Furthermore, Hamiltonians of physical systems are usually local, as they
reflect the interactions of the quantum mechanical system, and exponentially
vanishing gradients can be avoided by choosing a circuit depth scaling
logarithmically in system size [42]. Alternatively, an application-specific
ansatz with few variational parameters can help mitigate barren plateaus, such
as circuits based on Hamiltonian evolutions [46, 47].

In addition to the evolution gradient, the DualQTE loss function gradient
$\bm{\nabla}_{\bm{\delta\theta}}\mathcal{L}$ depends on the gradient of the
fidelity, which relies on measuring a global observable. This can be seen by
writing the fidelity of two $n$-qubit states prepared by unitaries
$U(\bm{\theta})$ and $U(\bm{\theta}^{\prime})$ as
$|\braket{0}{U^{\dagger}(\bm{\theta})U(\bm{\theta}^{\prime})}{0}|^{2}=\braket{\lambda}{P_{0}}{\lambda}$,
where $\ket{\lambda}=U^{\dagger}(\bm{\theta}^{\prime})U(\bm{\theta})\ket{0}$
and $P_{0}=\ket{0}\bra{0}^{\otimes n}$ is the global projector on the all-zero
state. Thus, evaluating the fidelity gradient for two randomly selected
parameter sets $\bm{\theta}$ and $\bm{\theta}^{\prime}$ would exhibit barren
plateaus at any circuit depth [42]. However, the optimization in DualQTE
starts at zero perturbations, $\bm{\theta}=\bm{\theta}^{\prime}$ where the
total state preparing unitary is the identity,
$\ket{\lambda}=U^{\dagger}(\bm{\theta})U(\bm{\theta})\ket{0}=\mathbb{I}\ket{0}$,
which is an initialization that is proven to not exhibit barren plateaus even
for global cost functions [48]. Together with the fact that the DualQTE loss
function is locally convex, the non-vanishing gradients at the initial point
of the optimization is a strong motivation for the efficient trainability of
DualQTE.

In Appendix E.4, we provide numerical evidence that for a local Hamiltonian
and a logarithmic-depth circuit, neither the evolution gradient or the
fidelity gradients decay exponentially with system size.

###  II.5 Sample complexity

The implementation of VarQTE on quantum hardware has several sources of
errors: the model $\ket{\phi(\bm{\theta})}$ could lack expressitivity to
capture the dynamics, the time integration scheme introduces errors, the QGT
and evolution gradient are subject to sampling error from a finite number of
measurement, and each operation is affected by hardware noise. If we denote
the ideal VarQTE parameters without sampling or hardware noise by
$\bm{\theta}(t)$ and the noisy parameters by $\tilde{\bm{\theta}}(t)$, the
error contributions can be split as

$\displaystyle\varepsilon(t)$
$\displaystyle=D_{B}(\phi(\tilde{\bm{\theta}}(t)),\Psi(t))$ (11)
$\displaystyle\leq\varepsilon_{M}(t)+\varepsilon_{S}(t),$

where we measured the error in Bures distance

$D_{B}(\psi,\phi)=\sqrt{2(1-|\braket{\psi}{\phi}|)},$ (12)

and distinguish in error due to lack of model expressitivity plus integration
error $\varepsilon_{M}(t)=D_{B}(\phi(\bm{\theta}(t)),\Psi(t))$ and error due
to a noisy implementation of VarQTE
$\varepsilon_{S}(t)=D_{B}(\phi(\tilde{\bm{\theta}}(t)),\phi(\bm{\theta}(t))$
[49].

Since the proposed DualQTE algorithm promises a reduction in measurement cost
of VarQTE, but is not concerned with the ansatz selection or hardware noise,
we here focus on investigating the scaling of the sampling error. The model
error $\varepsilon_{M}$ can be bounded with a-posteriori errorbounds [50],
which we also investigate for the real-time evolution case in Sec. IV.

Deriving a concrete bound in terms of system quantities such as the energy or
the number of parameters requires an assumption on the circuit structure. Here
we assume a circuit with only Pauli rotations $R_{P}(\theta_{i})$ where each
parameter $\theta_{i}$ is unique and does not have coefficients. Note that the
bounds can be adjusted for different circuit structures and parameterizations.
In addition, we assume a cutoff $\delta_{c}>0$ on the smallest eigenvalue of
$g$ due to a regularization of the linear system.

Then, we can state the following upper bound on the number of samples required
to achieve a sampling error of $\varepsilon_{S}$,

$N\leq\mathcal{O}\left(\frac{d^{3}E_{\text{max}}^{2}\Delta_{t}^{2}}{\delta_{c}^{4}\varepsilon_{S}^{2}}\right),$
(13)

where $E_{\text{max}}$ is the maximal eigenvalue of the Hamiltonian.

In contrast to a similar approach described in Ref. [49], we state the upper
bound in terms of the number of parameters in the model or the system’s
Hamiltonian, instead of the QGT and evolution gradients. Further, we are able
to derive a tighter result by leveraging Latala’s theorem from random matrix
theory to upper bound the sampling error in $g$ [51].

Since the DualQTE algorithm does not construct the QGT directly but only the
evolution gradient, we expect a reduction of a factor $d$ in the complexity,
and an additional factor for the number of optimization steps $K$ in each
timestep. Indeed we can show that the upper bound for the number of samples is

$N\leq\mathcal{O}\left(\frac{d^{2}K^{2}\Delta_{t}^{2}}{\delta\tau^{2}\varepsilon_{S}^{2}}\left(\frac{1}{\delta\tau}+E_{\text{max}}\right)^{2}\right).$
(14)

The detailed derivation of both bounds is described in Appendix D.

While it is possible to construct circuits where each component of these
bounds are tight Sec. III shows that in practice the actual number of required
samples scales less than this upper bound, which is further discussed in the
Appendix.

Figure 2:  (a) The mean and standard deviation of DualQITE and VarQITE, each
averaged over 5 independent experiments for a varying number of shots. (b) The
accuracy measured in integrated Bures distance $I_{B}$ (see Eq. (16)) for
DualQITE and VarQITE, with mean and standard deviation of 5 experiments. The
resources are measured in total number of measurements and are shown for
evaluation of gradients (and the QGT, in the case of VarQITE) using parameter-
shift rules (dashed) or a LCU approach (solid lines).

##  III Imaginary-time evolution

In this section, we show the results of DualQITE and investigate the circuit
costs compared to VarQITE. As an application, we use our algorithm as a
subroutine to prepare typical thermal states of the Heisenberg model, which
are then used to calculate the energy per site as a thermodynamic observable.
All circuits are constructed and simulated using Qiskit [52].

###  III.1 Heisenberg model

We simulate the imaginary-time evolution of the Heisenberg model with nearest-
neighbor interaction on a 12-qubit circle in a transverse field:

$H=J\sum_{\braket{ij}}\left(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\right)+g\sum_{i}Z_{i},$
(15)

with interaction strength $J=1/4$ and transversal field strength $g=-1$. As a
variational ansatz, we use a circuit with Pauli-$Y$ and Pauli-$Z$ single qubit
rotation layers that alternate with pairwise CNOT entangling gates. The
circuit structure is shown in Appendix E.1 and we use $r=3$ repetitions. The
initial state for the evolution is the equal superposition of all qubits,
$\ket{+}^{\otimes n}$, which we prepare by setting the rotation angles of the
last Pauli-$Y$ layer to $\pi/2$ and the remaining angles to $0$.

The optimization problems in DualQITE are solved using gradient descent with a
fixed learning rate of $\eta=0.1$ and time perturbation $\delta\tau=0.01$. The
initial iteration performs $100$ update steps and the subsequent, warmstarted
iterations, only $10$. These values are motivated by the simulations shown in
the Appendix E.2 and are partially heuristic, as a termination criterion is
challenging to define with access only to noisy loss function and gradient
evaluations. The parameters are integrated with an explicit Euler scheme with
timestep $\Delta_{t}=0.01$, i.e.

$\bm{\theta}(t+\Delta_{t})=\bm{\theta}(t)+\Delta_{t}\dot{\bm{\theta}}(t)=\bm{\theta}(t)+\Delta_{t}\frac{\bm{\delta\theta}}{\delta\tau}.$

Note that the integration timestep $\Delta_{t}$, which determines the accuracy
of the time integration, can be chosen differently from the time perturbation
$\delta\tau$, which affects the approximation error of the QGT metric with the
infidelity.

We compare the performance to VarQITE with the same integration scheme and use
an L-curve regularization [53] for a stable solution of the linear system.
Among all regularization techniques we attempted, such as adding a diagonal
shift, truncating small or negative singular values or solving on a stable
subsystem, the L-curve regularization provided the most accurate and stable
results.

In Fig. 2(a), we present the results for a varying number of shots along with
the exact time evolution based on exact diagonalization. Already with as
little as $100$ measurements per circuit evaluation (shots) on the $12$-qubit
model, the dual time evolution is able to qualitatively follow the imaginary-
time evolution and, up to time $t\approx 1$, even outperform VarQITE with
$1024$ shots. Increasing the number of measurements of DualQITE to $1024$
shots allows the dual method to closely track the exact solution towards the
ground state, with a higher accuracy than VarQITE with $8192$ shots.

###  III.2 Resource requirements

In the above experiment, DualQITE requires fewer circuit evaluations to
achieve the same accuracy as VarQITE. To investigate the total resource
requirements, we perform both DualQITE and VarQITE with different resources
and compute the achieved error. Since we are interested in following the
imaginary-time dynamics as closely as possible at each timestep, we define the
error as the average integrated Bures distance to the exact solution over the
time evolution,

$I_{B}(T)=\frac{1}{T}\int_{0}^{T}D_{B}(\phi(\bm{\theta}(t)),\psi(t))\mathrm{d}t.$
(16)

The state fidelity is computed exactly, i.e., we compute the state vector of
the model $\ket{\phi}$ at variational parameters $\bm{\theta}(t)$, and take
the inner product with the exact time-evolved state $\ket{\Psi(t)}$.

The results for an integration time of $T=2$ are shown in Fig. 2(b). We show
the integrated Bures distance with respect to the total number of measurements
recorded during the time evolution. In DualQITE, the resources can be split
between using more optimization steps in each timestep or more shots to
evaluate the gradients. The algorithm settings are detailed in Appendix E.3.
The figure shows the resource counts for gradient calculations via the
parameter-shift rule (PSR) and linear combination of unitaries (LCU). The LCU
technique requires additional auxiliary qubits and additional non-local
operations, but less overall circuits than PSR. For $P$ Pauli terms in the
Hamiltonian, the total number of required circuits $C$ per timestep

$\displaystyle C^{\text{VarQITE}}_{\text{LCU}}$
$\displaystyle=\frac{d(d+5)}{2}+Pd$ (17) $\displaystyle
C^{\text{VarQITE}}_{\text{PSR}}$ $\displaystyle=2d(d+P+1).$

For DualQTE the number of circuits is

$C^{\text{DualQITE}}_{\text{LCU}}=Pd+Kd,$ (18)

and $C^{\text{DualQITE}}_{\text{PSR}}=2C^{\text{DualQITE}}_{\text{LCU}}$,
where $K$ is the number of optimization steps per timestep. The total number
of measurements $N$ is obtained by multiplying the number of circuits with the
number of shots.

We see that, on average, DualQITE requires about one order of magnitude fewer
measurements to achieve the same accuracy as VarQITE. With an increasing
number of parameters, we expect this difference to grow, since VarQITE scales
as $\mathcal{O}(d^{2})$ whereas our algorithm, with warm starting, only
computes small corrections at each time step.

###  III.3 Sample complexity

In addition to the fixed-size model with 12 qubits, we investigate how the
resource requirements scale with system size. We compare VarQITE and DualQITE
for the Heisenberg model from Eq. (15) with varying number of spins $n$ and
the same circuit structure as before, but with an adjusted number of
repetitions of $r=\lceil\log_{2}(n)\rceil$ times, plus a final rotation layer.
We then tune the settings of VarQITE and DualQITE to achieve a mean accuracy
of $I_{B}\leq 0.1$ over 5 experiments and count the total number of required
measurements $N$. This threshold corresponds to a per-timestep fidelity of
0.995.

The results are presented in Fig. 3, which show the improved scaling of
DualQITE compared to VarQITE. For small system sizes and few parameters, the
overhead of solving the optimization problem in DualQITE is larger than
evaluating the QGT. But, as we increase the problem size, the quadratic
scaling of VarQITE takes over and our algorithm becomes more efficient.

This experiment allows to validate the upper bound on the number of
measurements of Sec. II.5. As shown in Fig. 3, the model error
$\varepsilon_{M}$ is negligible in comparison to the sampling error
$\varepsilon_{S}$ and we approximately have $\varepsilon_{S}\approx
I_{B}\approx 0.1$. The maximal energy of the Heisenberg model on a periodic
chain scales with the number of spins, and can be bounded by
$E_{\text{max}}=\mathcal{O}(n)\leq\mathcal{O}(n\log(n))=\mathcal{O}(d)$.
Inserting these values in Eqs. (13) and (14) we expect the scaling to be upper
bounded by $\mathcal{O}(d^{5})$ for VarQTE and $\mathcal{O}(d^{4}K)$ for
DualQTE. In practice, we observe approximately a scaling of $d^{3.56}$ for
VarQTE and $d^{2.29}$ for DualQTE, which shows the expected improved scaling
for our algorithm. The measured scaling also suggests that the bounds are not
yet tight, which we discuss further in Appendix D.

Figure 3: (a) Total number of measurements required to achieve a mean accuracy
of $I_{B}\leq 0.1$ over an average of 5 experiments. See Table 2 for the exact
algorithm settings. Dotted and dashed lines show fits for the number of
measurements. The bumps in the fits are due to the discontinuity of the
circuit depth, which depends on $\lceil\log_{2}(n)\rceil$. VarQITE is not
evaluated for $n=14$ qubits as it requires too many measurements. (b) Mean
accuracy and standard deviation of each point of the top panel. The grey line
indicates the infidelity threshold of $I_{B}=0.1$.

###  III.4 Calculating thermodynamic observables

As an application of imaginary-time evolution, we calculate thermodynamic
observables using the quantum minimally entangled thermal states algorithm
(QMETTS) [54, 8]. While the classical METTS algorithm has been specifically
developed for Matrix Product State simulations, the thermal state preparation
is still costly, and classical simulations fail if the system produces
macroscopic entanglement during the imaginary time evolution (e.g. for low-
temperature, 2D systems). Due to these restrictions, QMETTS is a promising
application for quantum imaginary-time evolution algorithms.

For an observable $A$ and inverse temperature $\beta$, the QMETTS algorithm
generates samples $\\{A_{m}\\}_{m}$ using a Markov chain whose average
approximate the ensemble average:

$\braket{A}_{\text{ens}}=\frac{\mathrm{Tr}(e^{-\beta
H}A)}{\mathrm{Tr}(e^{-\beta H})}\approx\frac{1}{M}\sum_{m=1}^{M}A_{m}.$

The sampling process to obtain the sample $A_{m}$ is

  1. 1.

Start from a product state $\ket{\phi_{m}(t=0)}$.

  2. 2.

Evolve up to imaginary time $t=\beta/2$

$\ket{\phi_{m}(\beta/2)}\propto e^{-\beta H/2}\ket{\phi_{m}(0)}.$

  3. 3.

Evaluate the observable to obtain the sample

$A_{m}=\braket{\phi_{m}(\beta/2)}{A}{\phi_{m}(\beta/2)}.$

  4. 4.

Measure $\ket{\phi_{m}(\beta/2)}$ in some basis to obtain the next random
product state $\ket{\phi_{m+1}(0)}$.

We investigate the Heisenberg model from Eq. (15) on a chain with $n=6$ spins
with parameters $J=1/4$ and $g=-1$. As a thermodynamic observable we compute
the energy per site, $\braket{H}/n$. To reduce the auto-correlation length in
the QMETTS Markov chain, and for faster convergence to the ensemble average,
it is favorable to measure in different bases in each step. Since the
Heisenberg Hamiltonian conserves the number of qubits in the $\ket{1}$ state,
avoiding the $Z$ basis greatly reduces the standard deviation of the Markov
chain. Thus, we here alternate between the $X$ and $Y$ basis for each sample.

As ansatz for DualQITE, we use problem-inspired circuits with pairwise CNOT
couplings and $r=2$ repetitions of rotation and entanglement layers, plus a
final rotation layer, see Appendix E.1 for a circuit diagram. For evolutions
of product states $\ket{\pm}$ in the $X$ basis, the rotation layers are single
qubit $R_{\mathrm{Y}}R_{\mathrm{Z}}$ gates, and for the states $\ket{\pm i}$
in the $Y$ basis, the layers implement $R_{\mathrm{X}}R_{\mathrm{Z}}$ gates.
The initial product states $\ket{\phi_{m}(0)}$ are prepared by setting the
parameters in the final layer rotation layer of the ansatz as follows,

$\displaystyle\ket{\pm}$ $\displaystyle\rightarrow
R_{\mathrm{Y}}\left(\frac{\pm\pi}{2}\right)R_{\mathrm{Z}}(0),$
$\displaystyle\ket{+i}$ $\displaystyle\rightarrow
R_{\mathrm{X}}\left(\frac{\pi}{2}\right)R_{\mathrm{Z}}(\pi),$
$\displaystyle\ket{-i}$ $\displaystyle\rightarrow
R_{\mathrm{X}}\left(\frac{\pi}{2}\right)R_{\mathrm{Z}}(0).$

Each energy sample is evaluated with $1024$ measurements per basis. The
optimization problem in DualQITE is solved with a time perturbation
$\delta\tau=0.01$ and gradient descent with a learning rate of $\eta=0.1$ and
$100$ iterations in the first timestep, followed by $10$ iterations in the
following, warmstarted timesteps. We integrate with a fixed timestep of
$\Delta_{t}=0.01$.

Figure 4: Energy per site for the Heisenberg model on a 6-spin chain,
comparing mean and standard deviation of QMETTS with DualQITE (blue circle and
errorbars) with a reference METTS implementation (black line and grey shade).

Figure 4 shows the estimated energy per site, along with the standard
deviation of the samples, for different inverse temperatures $\beta$. For the
alternating $X-Y$ basis, the Markov chain converges quickly and $M=25$ samples
suffice for an accurate estimate of the observable. For the imaginary-time
evolution, we compare DualQITE with the same settings as in the previous
sections to an exact evolution performed with matrix exponentials. It shows
that using the dual method allows to reliably reproduce the mean and standard
deviation of the Markov chain samples compared to the exact reference METTS.

##  IV Real-time evolution

The focus of this paper is on imaginary-time evolution as, to date, no other
QGT-free time evolution algorithms exist in this setting. For real-time
evolution, p-VQD has a similar structure as our algorithm and solves an
optimization problem rather than evaluating the QGT. However, there are key
differences to the dual algorithm applied to real-time evolution.

The p-VQD algorithm [20] projects a single Suzuki-Trotter step onto the
circuit model by solving the following optimization problem:

$\bm{\theta}(t+\Delta_{t})=\underset{\bm{\theta}^{\prime}}{\mathrm{argmax}}\big{|}\braket{\phi(\bm{\theta}^{\prime})}{e^{-iH\Delta_{t}}}{\phi(\bm{\theta}(t))}\big{|}^{2}.$

For Hamiltonians with many Pauli terms or long-range interactions, such as
those arising in molecular dynamics, the single step might already lead to
large circuits with non-local gates. While DualQRTE requires an LCU method to
evaluate the imaginary part of the energy and state gradients, see Eq. (5),
this only adds a single entangling gate, compared to a full Suzuki-Trotter
step is required. Furthermore, our dual time evolution allows the evaluation
of error bounds at almost no additional cost, which is not possible in p-VQD.
Due to these differences, this section presents DualQRTE: the dual time
evolution for real-time evolution.

###  IV.1 Heisenberg model

We present the real-time evolution under the Heisenberg Hamiltonian of Eq.
(15) on a linear chain with $n=4$ spins with parameters $J=1/4,\;g=-1$. As
variational model, we use a circuit with alternating Pauli-$X$ and Pauli-$Y$
rotation layers, and Pauli-$ZZ$ entangling gates that reflect the connectivity
of the spins. The circuit structure is visualized in Appendix F and, in this
experiment, all algorithms use $r=3$ repetitions of the rotation as well as
entangling gates. To prepare the initial state, $\ket{+}^{\otimes 4}$, we set
the parameters of the final Pauli-$Y$ rotations to $\pi/2$ and the rest to 0.

During the evolution, we track the average magnetization in the $X$ and $Z$
direction,

$\braket{X}=\frac{1}{n}\sum_{i=1}^{n}\braket{X_{i}},~{}~{}\braket{Z}=\frac{1}{n}\sum_{i=1}^{n}\braket{Z_{i}}.$

Since this Heisenberg Hamiltonian preserves the qubit excitations, and the
initial state is the equal superposition, the $\braket{Z}$ expectation value
should remain $0$ throughout the evolution.

The results of the different time evolution algorithms for an integration time
of $T=2$ and timestep $\Delta_{t}=T/100$ are presented in Fig. 5. Both
DualQRTE and p-VQD accurately track the observables using only 200 shots per
circuit. With the same resources, VarQRTE, on the other hand, has lower
accuracy and we need to use 1024 shots per circuit to match the result of the
optimization-based algorithms.

Figure 5: Average magnetization in $X$ and $Z$ direction as tracked by
different variational algorithms.

###  IV.2 Error bounds

In variational real-time evolution, the model error in terms of Bures distance
$D_{B}$ due to restriction to the variational manifold can be expressed as
[50]

$\displaystyle\dot{\varepsilon}_{M}$
$\displaystyle:=\left\|\sum_{k=1}^{d}\dot{\theta}_{k}\ket{\partial_{k}\phi(\bm{\theta})}+iH\ket{\phi(\bm{\theta})}\right\|_{2}^{2}$
(19)
$\displaystyle=\mathrm{Var}(H|\phi(\bm{\theta}))+\dot{\bm{\theta}}^{T}g(\bm{\theta})\dot{\bm{\theta}}-2\dot{\bm{\theta}}^{T}\bm{b}(\bm{\theta}),$

where we set $\hbar\equiv 1$. Integrating this error rate provides an upper
bound on the Bures distance, that is

$D_{B}(\phi(\bm{\theta}(T)),\Psi(T))\leq\int_{0}^{T}\dot{\varepsilon}_{M}(t)\mathrm{d}t,$

where $\ket{\Psi(t)}$ is the exact time-evolved state and the time-dependence
of $\varepsilon_{M}$ is due to the time-dependence of the parameters
$\bm{\theta}=\bm{\theta}(t)$.

Up to the variance
$\mathrm{Var}(H|\phi(\bm{\theta}))=\braket{\phi(\bm{\theta})}{H^{2}}{\phi(\bm{\theta})}-(\braket{\phi(\bm{\theta})}{H}{\phi(\bm{\theta})})^{2}$
, this error is proportional to the loss function used in DualQRTE. By using
the same expansion $\dot{\bm{\theta}}=\bm{\delta\theta}/\delta\tau$ and using
the infidelity to approximate the inner product with respect to the geometric
tensor, we can rewrite the error as

$\displaystyle\dot{\varepsilon}_{M}$
$\displaystyle=\mathrm{Var}(H|\phi(\bm{\theta}))+\frac{1-F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})}{(\delta\tau)^{2}}-\frac{2\bm{\delta\theta}^{T}\bm{b}(\bm{\theta})}{\delta\tau}+\mathcal{O}(\delta\tau)$
$\displaystyle=\mathrm{Var}(H|\phi(\bm{\theta}))+\frac{2\mathcal{L}(\bm{\delta\theta})}{(\delta\tau)^{2}}+\mathcal{O}(\delta\tau).$

Note that the error scales linearly in time perturbation $\delta\tau$ as the
infidelity approximation has a cubic error term [33], which is divided by the
square of the perturbation. If we, for example, use a forward Euler rule with
timestep $\Delta_{t}$ to integrate the variational error, the integration
error scales as $\mathcal{O}(\Delta_{t}+T\delta\tau)$. This highlights the
importance of differentiating between the timestep $\Delta_{t}$ for the
integration, and the time-perturbation $\delta\tau$ to approximate the
derivative.

In Fig. 6, we show the error bounds along with the true error for the time
evolution of the Heisenberg model. The bounds are computed for different
timesteps $\Delta_{t}$ for VarQRTE, and for DualQRTE for a fixed time
perturbations $\delta\tau=10^{-3}$ in exact simulations. Firstly, we can
verify that the error bounds hold. Secondly, the larger the timestep relative
to the time-perturbation, the more accurate the approximation of the dual time
evolution, as the error $\mathcal{O}(\Delta_{t}+T\delta\tau)$ is dominated by
the integration error.

Figure 6: Error of the real-time evolution in Bures distance, plus error
bounds obtained with VarQRTE and DualQRTE.

##  V Conclusion

In this paper, we present a novel algorithm for variational quantum time
evolution that does not require the evaluation of the QGT, but instead solves
a dual optimization problem in each timestep. The proposed dual time evolution
algorithm, DualQTE, is particularly interesting for imaginary-time evolution,
as there is currently no alternative variational algorithm able to circumvent
the $\mathcal{O}(d^{2})$ cost of VarQITE. For real-time evolution, p-VQD also
offers an optimization-based approach by projecting a single Trotter step onto
the variational form. In comparison, the dual time evolution has the advantage
that no Suzuki-Trotter step has to be implemented, which could require deep
circuits or non-local operations, depending on the Hamiltonian. Furthermore,
our algorithm allows to evaluate variational error bounds [50], although how
accurately they can be evaluated in the presence of shot noise remains an open
question.

We demonstrated DualQTE for the imaginary-time evolution of a Heisenberg
Hamiltonian on 12 qubits, and found that, in this setting, it requires about
one order of magnitude less measurements to achieve the same accuracy as
VarQITE. As a practical application of imaginary-time evolution, we calculated
thermodynamic observables with the QMETTS algorithm and showed that the
DualQITE is suitable to reproduce the sampling distributions. Finally, we
applied our algorithm to an illustrative example for real-time evolution,
where it produced comparable results to p-VQD for the same amount of
resources, while both algorithms outperformed VarQRTE.

In the presented experiments, we used standard gradient descent algorithms
with a fixed learning rate. We expect that the performance could be further
improved by using more advanced optimization schemes, or methods that also
take into account information from previous iterations. Another possible
improvement would be a suitable termination criterion for noisy evaluations of
the loss function. As for other optimization-based time evolution algorithms,
such as p-VQD, it remains challenging to accurately measure the fidelity in
the presence of hardware noise.

In conclusion, the proposed DualQTE is an efficient variational algorithm for
quantum time evolution that does not suffer from the quadratic complexity of
evaluating the QGT. This cost reduction enables scaling imaginary-time
evolution to larger, practically relevant system sizes and allows the
simulation and demonstration of a wide variety of important tasks such as
Gibbs state preparation, mixed time evolution, or the evaluation of
thermodynamic observables. Improving the resource requirements for near-term
algorithms is an important step for scaling demonstrations to the full size of
today’s quantum computers and work towards practical applications.

##  VI Acknowledgements

We thank Christa Zoufal, Stefano Barison, Almudena Carrera Vazquez, David
Sutter, Caroline Tornow, Laurin Fischer and Daniel Egger for insightful
conversations on this project.

We acknowledge the use of IBM Quantum services for this work. The views
expressed are those of the authors, and do not reflect the official policy or
position of IBM or the IBM Quantum team.

IBM, the IBM logo, and ibm.com are trademarks of International Business
Machines Corp., registered in many jurisdictions worldwide. Other product and
service names might be trade- marks of IBM or other companies. The current
list of IBM trademarks is available at https://www.ibm.com/legal/copytrade.

## References

  * Zhang  _et al._ [2017] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Nature 551, 601 (2017).
  * Dborin  _et al._ [2022] J. Dborin, V. Wimalaweera, F. Barratt, E. Ostby, T. E. O’Brien, and A. G. Green, Nature Communications 13, 5977 (2022).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 595, 227 (2021).
  * Altman [2018] E. Altman, Nature Physics 14, 979 (2018).
  * de Jong  _et al._ [2022] W. A. de Jong, K. Lee, J. Mulligan, M. Płoskoń, F. Ringer, and X. Yao, Phys. Rev. D 106, 054508 (2022).
  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, npj Quantum Information 5, 75 (2019).
  * Jones  _et al._ [2019] T. Jones, S. Endo, S. McArdle, X. Yuan, and S. C. Benjamin, Phys. Rev. A 99, 062304 (2019).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O. Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandao, and G. K.-L. Chan, Nature Physics 16, 205 (2020).
  * Getelina  _et al._ [2023] J. C. Getelina, N. Gomes, T. Iadecola, P. P. Orth, and Y.-X. Yao, arXiv:2301.02592 (2023).
  * Barkoutsos  _et al._ [2018] P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis, A. Fuhrer, M. Ganzhorn, D. J. Egger, M. Troyer, A. Mezzacapo, S. Filipp, and I. Tavernelli, Phys. Rev. A 98, 022322 (2018).
  * Zoufal  _et al._ [2023] C. Zoufal, R. V. Mishmash, N. Sharma, N. Kumar, A. Sheshadri, A. Deshmukh, N. Ibrahim, J. Gacon, and S. Woerner, Quantum 7, 909 (2023).
  * Zoufal  _et al._ [2021a] C. Zoufal, A. Lucchi, and S. Woerner, Quantum Machine Intelligence 3, 7 (2021a).
  * Miessen  _et al._ [2023] A. Miessen, P. J. Ollitrault, F. Tacchino, and I. Tavernelli, Nature Computational Science 3, 25 (2023).
  * Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
  * Zhao  _et al._ [2022] H. Zhao, M. Bukov, M. Heyl, and R. Moessner, arXiv:2209.12653 (2022).
  * Mukhopadhyay  _et al._ [2023] P. Mukhopadhyay, N. Wiebe, and H. T. Zhang, npj Quantum Information 9, 31 (2023).
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, Quantum 3, 191 (2019).
  * IBM Quantum [2021] IBM Quantum, IBM Unveils Breakthrough 127 Qubit Quantum Processor (2021).
  * IBM Quantum [2022] IBM Quantum, IBM Unveils 400 Qubit Plus Quantum Processor and Next Generation IBM Quantum System Two (2022).
  * Barison  _et al._ [2021] S. Barison, F. Vicentini, and G. Carleo, Quantum 5, 512 (2021).
  * Barratt  _et al._ [2021] F. Barratt, J. Dborin, M. Bal, V. Stojevic, F. Pollmann, and A. G. Green, npj Quantum Information 7, 10.1038/s41534-021-00420-3 (2021).
  * Lin  _et al._ [2021] S.-H. Lin, R. Dilip, A. G. Green, A. Smith, and F. Pollmann, PRX Quantum 2, 010342 (2021).
  * Benedetti  _et al._ [2021] M. Benedetti, M. Fiorentini, and M. Lubasch, Phys. Rev. Res. 3, 033083 (2021).
  * Slattery  _et al._ [2022] L. Slattery, B. Villalonga, and B. K. Clark, Phys. Rev. Res. 4, 023072 (2022).
  * Commeau  _et al._ [2020] B. Commeau, M. Cerezo, Z. Holmes, L. Cincio, P. J. Coles, and A. Sornborger, arXiv:2009.02559 (2020).
  * Cîrstoiu  _et al._ [2020] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, npj Quantum Information 6, 82 (2020).
  * Gibbs  _et al._ [2021] J. Gibbs, K. Gili, Z. Holmes, B. Commeau, A. Arrasmith, L. Cincio, P. J. Coles, and A. Sornborger, Long-time simulations with high fidelity on quantum hardware (2021).
  * Lim  _et al._ [2021] K. H. Lim, T. Haug, L. C. Kwek, and K. Bharti, Quantum Science and Technology 7, 015001 (2021).
  * Mc Keever and Lubasch [2023] C. Mc Keever and M. Lubasch, Phys. Rev. Res. 5, 023146 (2023).
  * Mansuroglu  _et al._ [2023] R. Mansuroglu, T. Eckstein, L. Nützel, S. A. Wilkinson, and M. J. Hartmann, Quantum Science and Technology 8, 025006 (2023).
  * Sbahi  _et al._ [2022] F. M. Sbahi, A. J. Martinez, S. Patel, D. Saberi, J. H. Yoo, G. Roeder, and G. Verdon, arXiv:2206.04663 (2022).
  * Amari [1998] S.-i. Amari, Neural Computation 10, 251 (1998).
  * Stokes  _et al._ [2020] J. Stokes, J. Izaac, N. Killoran, and G. Carleo, Quantum 4, 269 (2020).
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Phys. Rev. A 99, 032331 (2019).
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. de Wolf, Phys. Rev. Lett. 87, 167902 (2001).
  * Cincio  _et al._ [2018] L. Cincio, Y. Subaşı, A. T. Sornborger, and P. J. Coles, New Journal of Physics 20, 113022 (2018).
  * Cleve  _et al._ [1998] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 339 (1998).
  * Havlíček  _et al._ [2019] V. Havlíček _et al._ , Nature 567, 209 (2019).
  * Elben  _et al._ [2019] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Phys. Rev. A 99, 10.1103/PhysRevA.99.052323 (2019).
  * Gacon  _et al._ [2021] J. Gacon, C. Zoufal, G. Carleo, and S. Woerner, Quantum 5, 567 (2021).
  * Spall [1998] J. C. Spall, JOHNS HOPKINS APL TECHNICAL DIGEST 19, 11 (1998).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Nature Communications 12, 1791 (2021), arXiv: 2001.00550.
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nature Communications 9, 4812 (2018), arXiv: 1803.11173.
  * Ortiz Marrero  _et al._ [2021] C. Ortiz Marrero, M. Kieferová, and N. Wiebe, PRX Quantum 2, 040316 (2021).
  * Wang  _et al._ [2021] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Nature Communications 12, 6961 (2021).
  * Ollitrault  _et al._ [2020] P. J. Ollitrault, A. Baiardi, M. Reiher, and I. Tavernelli, Chemical Science 11, 6842 (2020).
  * Park and Killoran [2023] C.-Y. Park and N. Killoran, arXiv:2302.08529 (2023).
  * Grant  _et al._ [2019] E. Grant, L. Wossnig, M. Ostaszewski, and M. Benedetti, Quantum 3, 214 (2019), arXiv:1903.05076 [quant-ph].
  * Endo  _et al._ [2020] S. Endo, J. Sun, Y. Li, S. C. Benjamin, and X. Yuan, Phys. Rev. Lett. 125, 010501 (2020).
  * Zoufal  _et al._ [2021b] C. Zoufal, D. Sutter, and S. Woerner, arXiv:2108.00022 (2021b).
  * Latała [2005] R. Latała, Proceedings of the American Mathematical Society 133, 1273 (2005).
  * Qiskit contributors [2023] Qiskit contributors, Qiskit: An open-source framework for quantum computing, 10.5281/zenodo.2562110 (2023).
  * Cultrera and Callegaro [2020] A. Cultrera and L. Callegaro, IOP SciNotes 1, 025004 (2020).
  * Stoudenmire and White [2010] E. M. Stoudenmire and S. R. White, New Journal of Physics 12, 055026 (2010).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Annual Review of Condensed Matter Physics 11, 369 (2020).
  * IBM Quantum [2023] IBM Quantum, https://quantum-computing.ibm.com/ (2023), accessed on 07/26/2023.
  * Tornow  _et al._ [2022] C. Tornow, N. Kanazawa, W. E. Shanks, and D. J. Egger, Phys. Rev. Appl. 17, 064061 (2022).
  * Egger  _et al._ [2018] D. J. Egger, M. Werninghaus, M. Ganzhorn, G. Salis, A. Fuhrer, P. Müller, and S. Filipp, Phys. Rev. Appl. 10, 044030 (2018).
  * Magnard  _et al._ [2018] P. Magnard, P. Kurpiers, B. Royer, T. Walter, J.-C. Besse, S. Gasparinetti, M. Pechal, J. Heinsoo, S. Storz, A. Blais, and A. Wallraff, Phys. Rev. Lett. 121, 060502 (2018).
  * Gentinetta  _et al._ [2022] G. Gentinetta, A. Thomsen, D. Sutter, and S. Woerner, arXiv:2203.00031 (2022).
  * Wang [1993] Y. H. Wang, Statistica Sinica 3, 295 (1993), publisher: Institute of Statistical Science, Academia Sinica.

##  Appendix A Runtime estimates of variational time evolution

The benchmark in Sec. III.3 provides a scaling for the total number of
measurements $N$ required by VarQITE and DualQITE, which allows a runtime
estimation on the algorithms on quantum hardware. In this estimation we
neglect the overhead of classical processors and assume a superconducting
quantum computer with a basis gate set including $\sqrt{X}$, $R_{\mathrm{Z}}$
and CX gates, as reported by several IBM Quantum backends, for example. This
gate set allows to compile any sequence of single qubit gates into two
$\sqrt{X}$ gates and three virtual $R_{\mathrm{Z}}$ gates. For an $n$-qubit
simulation of the Heisenberg model and the considered circuit model (see Fig.
10) with $r$ repetitions, the time for a single measurement can then be
approximated as

$t_{\text{shot}}=2rt_{\text{CX}}+2(r+1)t_{\sqrt{X}}+t_{\text{meas}}+t_{\text{reset}},$
(20)

where $t_{\text{CX}}$ is the duration of a CX gate, $t_{\sqrt{X}}$ the
duration of a $\sqrt{X}$ gate, $t_{\text{meas}}$ the time of a measurement and
$t_{\text{reset}}$ the time to reset the qubits for the next execution. Since
the $R_{\mathrm{Z}}$ gates are virtual they do not contribute to the runtime.
The total runtime is then estimated by $Nt_{\text{shot}}$.

Depending on the architecture and the gate decomposition the duration and
fidelity of single- and two-qubit operation, as well as measurements, varies
on superconducting qubit chips [55]. Here we use gate times of
$t_{\text{CX}}=451ns$, $t_{\sqrt{X}}=36ns$ and $t_{\text{meas}}=860ns$ as
reported by ibm_peekskill (v2.6.5), which is an IBM Quantum Falcon processors
[56]. For shallow circuits in particular, the time to reset qubits for the
following execution is a crucial bottleneck. The reset operation can, for
example, be implemented by waiting 5-10 T1 times and let the qubits decay to
the computational ground-state, but with T1 times of the order of $400\mu s$
the reset via relaxation is orders of magnitude slower than the other circuit
operations. Active resets instead measure the qubit state and apply an
$X$-operation conditionally if the state $\ket{1}$ is measured. This technique
allows to reduce the reset times to typically $50$ to $250\mu s$ on IBM
hardware [57], which, however, still dominates the overall runtime for the
considered circuits. By using a second excited state it is possible to
implement reset schemes with $500ns$ to $2\mu s$ [58, 59] and we therefore use
$t_{\text{reset}}=2\mu s$ in our estimation.

##  Appendix B Derivation via quantum natural gradient descent

VarQITE is inherently connected to the quantum natural gradient (QNG)
algorithm [33]. In fact, this connection is a motivation for the convergence
of the QNG as imaginary-time evolution is guaranteed to converge to the ground
state, if there is sufficient initial overlap with it.

With a forward Euler integration the VarQITE update rule is

$\bm{\theta}^{(t+1)}=\bm{\theta}^{(t)}+\Delta_{t}g^{-1}(\bm{\theta}^{(t)})\left(-\frac{\bm{\nabla}E(\bm{\theta}^{(t)})}{2}\right).$

This coincides with the QNG update step for the loss function
$\ell(\bm{\theta})=E(\bm{\theta})/2$ and a learning rate of $\eta=\Delta_{t}$,

$\bm{\theta}^{(t+1)}=\bm{\theta}^{(t)}-\eta
g^{-1}(\bm{\theta}^{(t)})\bm{\nabla}\ell(\bm{\theta}^{(t)}).$ (21)

The natural gradients step can be expressed in a dual formulation as

$\bm{\theta}^{(t+1)}=\underset{\bm{\theta}}{\mathrm{argmin}}\braket{\bm{\nabla}\ell(\bm{\theta}^{(t)},\bm{\theta}-\bm{\theta}^{(t)}}+\frac{1}{2\eta}d^{2}(\bm{\theta},\bm{\theta}^{(t)}),$

with a distance metric $d$. In this equation we see that the update step is
going into the opposite direction of the gradient $\bm{\nabla}\ell$, while the
magnitude is limited by the distance metric and the learning rate.

Standard gradient descent uses the model-agnostic $\ell_{2}$ norm as distance
metric. Natural gradients on the other hand limit the update step by the
amount of change it induces in the model. To measure the induced change the
metric $d$ is chosen to be the Fubini-Study metric, which, as shown in Ref.
[33], if locally approximated, yields the QGT:

$\displaystyle d^{2}(\phi(\bm{\theta}),\phi(\bm{\theta}+\bm{\delta\theta}))$
$\displaystyle=\arccos^{2}|\braket{\phi(\bm{\theta})}{\phi(\bm{\theta}+\bm{\delta\theta})}|$
$\displaystyle=1-|\braket{\phi(\bm{\theta})}{\phi(\bm{\theta}+\bm{\delta\theta})}|^{2}+\mathcal{O}(||\bm{\delta\theta}||_{2}^{4})$
$\displaystyle=\braket{\bm{\delta\theta},g(\bm{\theta})\bm{\delta\theta})}+\mathcal{O}(||\bm{\delta\theta}||_{2}^{3}).$

The formulation in Eq. (21) is then obtained by solving the minimization
problem.

To circumvent the explicit evaluation of the QGT the natural gradient update
can instead be calculated without the quadratic local approximation, and
instead solve the optimization problem directly. If we use the infidelity as
distance metric and replace the loss function gradient by the evolution
gradient
$\bm{\nabla}\ell(\bm{\theta})=\bm{\nabla}E(\bm{\theta})/2=-\bm{b}(\bm{\theta})$,
we obtain the same update rule as the main text

$\bm{\theta}^{(t+1)}=\underset{\bm{\theta}}{\mathrm{argmin}}-\braket{\bm{b}(\bm{\theta}^{(t)}),\bm{\theta}-\bm{\theta}^{(t)}}+\frac{1-F(\bm{\theta},\bm{\theta}^{(t)})}{2\Delta_{t}}.$

##  Appendix C Illustrative example

For an intuitive understanding of the approximations of the QGT norm, we
investigate an illustrative example with the variational model
$\ket{\phi(\theta)}=R_{\mathrm{Z}}(\theta)R_{\mathrm{Y}}(\theta)\ket{0}$, the
Hamiltonian $H=Z$ and a timestep of $\delta\tau=1/2$. In Fig. 7(a) we compare
the exact values of the loss function $\mathcal{L}$ for imaginary-time
evolution around $\theta=\pi/4$ obtained by using the metric
$\braket{\delta\theta,g(\theta)\delta\theta}$ or the infidelity
$1-F(\theta,\theta+\delta\theta)$ as norm.

Figure 7: (a) Values of the loss function $\mathcal{L}$ for evaluation with
the QGT metric, and with introduced infidelity approximation. (b) Difference
of the QGT metric and infidelity as function of the perturbation
$\delta\theta$.

At $\delta\theta=0$ the approximation is exact and in the vicinity the
difference scales as $(\delta\theta)^{3}$, see also Fig. 7(b). Note, that the
infidelity is periodic and bounded in $[0,1]$ but the linear term
$b^{T}\delta\theta$ is unbounded, which leads to the fact that the minimum of
the infidelity-based loss function close to $\delta\theta=0$ is not the global
minimum. This is well visible in Fig. 7(a) where the infidelity-based loss
function achieves lower values for large $\delta\theta$ than the minimum of
the QGT close to $\delta\theta=0$. Since we aim to find the same minimum as
the QGT-based loss function using a local optimization routine, such as
gradient descent, is crucial for the dual time evolution.

### Impact of the time perturbation

The approximation error scales with the norm of $\bm{\delta\theta}$ and,
therefore, solving for the update step
$\dot{\bm{\theta}}=\bm{\delta\theta}/\delta\tau$ with a smaller time
perturbation $\delta\tau$ should result in a smaller error in the update step.
Remembering the definition of the loss function

$\mathcal{L}(\bm{\delta\theta})=\frac{1-F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})}{2}-\delta\tau\cdot\bm{b}^{T}(\bm{\theta})\bm{\delta\theta},$

we see that a smaller $\delta\tau$ moves the minimum closer to the minimum of
the infidelity at $\bm{\delta\theta}=\bm{0}$, leading to a smaller
approximation error. Since the fidelity is bounded but the linear part
$\bm{b}^{T}(\bm{\theta})\bm{\delta\theta}$ is not, there is a maximum feasible
range for the value of $\delta\tau$. A necessary condition for the existence
of the minimum is that the gradient of the loss function vanishes,
$\bm{\nabla}\mathcal{L}=0$, which requires

$\forall
i\in\\{1,\dots,d\\}:\frac{1}{2}\frac{\partial}{\partial(\delta\theta)_{i}}F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})=-\delta\tau\cdot
b_{i}(\bm{\theta}).$ (22)

For a circuit with unique parameters and only Pauli rotations gates, the
gradient of the fidelity can be bounded via the parameter-shift rule to be in
$[-1/2,1/2]$ (see also Appendix D). Thus, a necessary condition for the
timestep perturbation is

$\forall
i\in\\{1,\dots,d\\}:\delta\tau\in\left[\frac{-1}{4|b_{i}(\bm{\theta})|},\frac{1}{4|b_{i}(\bm{\theta})|}\right],$
(23)

which can be generalized to circuits with repeated parameters or other than
Pauli gates. Note that this is only a necessary and not a sufficient condition
for the existence of a minimum since, depending on the circuit structure, the
fidelity gradient may not support the full range $[-1/2,1/2]$.

In Fig. 8(a) we visualize the impact of $\delta\tau$ on the loss landscape.
For small time perturbations the QGT-based and dual loss landscapes almost
coincide, but if $\delta\tau$ is chosen too large the dual loss function has
no minimum. If the loss function can be evaluated exactly, choosing
$\delta\tau$ as small as possible therefore minimizes the approximation error.
In Fig. 8(b), we find the the error in the parameter derivative
$\dot{\bm{\theta}}=\bm{\delta\theta}/\delta\tau$ scales approximately as
$\mathcal{O}(\delta\tau)$. In practice, however, the loss function is subject
to measurement noise and errors in the solution $\bm{\delta\theta}$ are
amplified by $1/\delta\tau$. Hence, for a finite number of measurements there
is a trade-off between QGT approximation error and controlling the noise
amplification.

Figure 8: (a) Loss landscapes and optimal solutions of the original, QGT-based
loss function and the dual loss function for different $\delta\tau$. (b) Error
in calculating the parameter derivative $\dot{\theta}$ depending on
$\delta\tau$ and the number of measurements $N$.

##  Appendix D Bound the sample complexity of VarQTE

In this section, we present the derivation on the upper bound of the sample
complexity of VarQTE and DualQTE. The target error is measured in integrated
Bures distance,

$\varepsilon_{S}=\frac{1}{T}\int_{0}^{T}\sqrt{2(1-|\braket{\phi(\bm{\theta})}{\phi(\tilde{\bm{\theta}})}|}\mathrm{d}t.$
(24)

Assuming a forward Euler integration, the Bures distance can be formulated in
terms of the QGT as

$\displaystyle\varepsilon_{S}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\sqrt{2\left(1-\sqrt{1-\Delta_{t}^{2}\Delta\dot{\bm{\theta}}^{T}g(\bm{\theta})\Delta\dot{\bm{\theta}}}\right)}\mathrm{d}t$
(25)
$\displaystyle=\frac{1}{T}\int_{0}^{T}\Delta_{t}\sqrt{\Delta\dot{\bm{\theta}}g(\bm{\theta})\Delta\dot{\bm{\theta}}}\mathrm{d}t$
$\displaystyle\leq\frac{1}{T}\int_{0}^{T}\Delta_{t}\|g(\bm{\theta})\|_{2}\|\Delta\dot{\bm{\theta}}\|_{2}\mathrm{d}t$
$\displaystyle\leq\Delta_{t}\sqrt{\lambda_{\text{max}}}\|\Delta\dot{\bm{\theta}}_{\mathrm{max}}\|_{2},$

where we introduced
$\Delta\dot{\bm{\theta}}=\tilde{\dot{\bm{\theta}}}-\dot{\bm{\theta}}$,
$\lambda_{\text{max}}\geq 0$ is a bound on the largest eigenvalue of $g$ for
any parameter value $\bm{\theta}$ and, similarly,
$\|\Delta\dot{\bm{\theta}}_{\mathrm{max}}\|_{2}$ an upper bound on the norm
$\Delta\dot{\bm{\theta}}$. In the first line, we dropped
$\mathcal{O}(\Delta_{t}^{3})$ error terms, in the second line we used a first
order Taylor-expansion and in the last line we use the definition of the
operator norm to bound the inner product of $\bm{\delta\theta}$ in the metric
of $g$.

###  D.1 VarQTE

In each VarQTE step we solve a linear system for the update step, where the
measurements of the QGT and evolution gradient are subject to sampling error.
We define the noisy quantities as
$\tilde{g}(\bm{\theta})=g(\bm{\theta})+\Delta g(\bm{\theta})$ and
$\tilde{\bm{b}}(\bm{\theta})=\bm{b}(\bm{\theta})+\Delta\bm{b}(\bm{\theta})$
and, then, solve the noisy linear system

$\tilde{g}(\bm{\theta})\tilde{\dot{\bm{\theta}}}=\tilde{\bm{b}}(\bm{\theta}),$
(26)

with the noisy update
$\tilde{\dot{\bm{\theta}}}=\dot{\bm{\theta}}+\Delta\dot{\bm{\theta}}$. To
stabilize the linear system and ensure the QGT and it’s estimate are
invertible, we assume a regularization of $g$ and $\tilde{g}$ in form of a
diagonal shift $\delta_{c}$. This shift is a trade-off of stability and bias,
which is also discussed in Ref. [40], Appendix D.

We can write the error in the update step using the difference of the noisy
and exact linear system solutions, as

$\displaystyle\|\Delta\dot{\bm{\theta}}\|_{2}$ $\displaystyle=\|(g+\Delta
g)^{-1}(\bm{b}+\Delta\bm{b})-g^{-1}\bm{b}\|_{2}$ (27)
$\displaystyle\approx\|(g^{-1}-g^{-1}\Delta
gg^{-1})(\bm{b}+\Delta\bm{b})-g^{-1}\bm{b}\|_{2}$
$\displaystyle=\|g^{-1}\Delta\bm{b}-g^{-1}\Delta gg^{-1}\bm{b}-g^{-1}\Delta
g^{-1}g^{-1}\Delta\bm{b}\|_{2}$
$\displaystyle\approx\|g^{-1}\Delta\bm{b}-g^{-1}\Delta
g\dot{\bm{\theta}}\|_{2}$
$\displaystyle\leq\|g^{-1}\|_{2}\left(\|\Delta\bm{b}\|_{2}+\|\Delta
g\|_{2}\|\dot{\bm{\theta}}\|_{2}\right),$

where we dropped the explicit parameter dependence for legibility. In the
second line we used the Neumann series to approximate $(g+\Delta
g)^{-1}=g^{-1}-g^{-1}\Delta gg^{-1}+\mathcal{O}(\|\Delta
g\|_{2}^{2}\|g^{-1}\|_{2}^{3})$ and dropped quadratic error terms on the
fourth line. In the following we derive upper bounds on the maximal value of
the individual contributions in the error bound, such that we finally obtain a
bound $\|\Delta\dot{\bm{\theta}}_{\mathrm{max}}\|$ on the error in the update
step.

#### Spectrum of $g$

Each QGT entry can be computed as [40]

$g_{ij}(\bm{\theta})=-\frac{1}{2}\partial_{i}\partial_{j}F(\bm{\theta}^{\prime},\bm{\theta})\Bigg{|}_{\bm{\theta}^{\prime}=\bm{\theta}}.$
(28)

For a circuit with unique, non-interacting parameter and only plain Pauli
rotation gates $R_{P}(\theta)$, we can use the parameter-shift rule [34] to
write the entry as

$g_{ij}(\bm{\theta})=-\frac{1}{2}\frac{F^{(++)}_{ij}-F^{(+-)}_{ij}-F^{(-+)}_{ij}+F^{(--)}_{ij}}{4},$
(29)

where we abbreviated
$F^{(\pm\pm)}_{ij}=F(\bm{\theta},\bm{\theta}\pm\bm{e}_{i}\pi/2\pm\bm{e}_{j}\pi/2)$
for the $i$th unit vector $e_{i}$ (and $j$th unit vector $e_{j}$). Since the
fidelity is in $[0,1]$ we can bound each entry by

$-\frac{1}{4}\leq g_{ij}(\bm{\theta})\leq\frac{1}{4},$ (30)

for any value of $\bm{\theta}$. Gershgorin’s circle theorem tells us that the
maximal eigenvalue of $g$ is bounded from above by the maximal sum over the
columns or rows, which in this case is achieved by setting all elements of a
column/row to $1/4$. This gives the bound

$\lambda_{\text{max}}\leq\sum_{i=1}^{d}\frac{1}{4}=\frac{d}{4}.$ (31)

This bound can be generalized to circuits with coefficients or repeated
parameters by applying the chain and product rules. For example, for a
coefficient-free circuit where parameters can be repeated up to $m$ times, the
bound becomes $md/4$.

#### Norm of the update step

The update step can be bounded as
$\|\dot{\bm{\theta}}\|_{2}\leq\|g^{-1}\|_{2}\|b\|_{2}$, where
$\|g^{-1}\|_{2}\leq\delta_{c}^{-1}$. The evolution gradient can be bounded
using the parameter-shift rule, under the circuit structure assumptions as the
previous section. Each element in the gradient is bounded by

$\displaystyle|b_{i}|=\frac{|E^{(+)}_{i}-E^{(-)}_{i}|}{2}\leq\frac{|E^{(+)}_{i}|+|E^{(-)}_{i}|}{2}\leq
E_{\text{max}},$ (32)

where $E^{(\pm)}_{i}=E(\bm{\theta}\pm\bm{e}_{i}\pi/2)$ and $E_{\text{max}}$ is
the absolute maximum system energy. The norm over all elements is then
$\|b\|_{2}\leq\sqrt{d}E_{\text{max}}$, leading to an overall bound of

$\|\dot{\bm{\theta}}\|_{2}\leq\frac{\sqrt{d}E_{\text{max}}}{\delta_{c}},$ (33)

for any parameter value $\bm{\theta}$.

#### Sampling errors

Since the measurement noise is unbiased, the random variable $\Delta
g=\tilde{g}-g$ has zero mean with i.i.d. entries. This allows to apply
Latala’s theorem [51], which states that

$\mathbb{E}[\|\Delta g\|_{2}]\leq
C\left(\max_{i}\sqrt{\sum_{j=1}^{d}\mathbb{E}[(\Delta
g)_{ij}^{2}]}+\max_{j}\sqrt{\sum_{i=1}^{d}\mathbb{E}[(\Delta
g)_{ij}^{2}]}+\sqrt[4]{\sum_{i,j=1}^{d}\mathbb{E}[(\Delta
g)_{ij}^{4}]}\right),$ (34)

for some constant $C\in\mathbb{R}$.

Ref. [60] is concerned with the similar case of sampling the matrix
$[F(x_{i},x_{j})]_{i,j=1}^{d}$ for a set of parameters
$\\{x_{i}\\}_{i=1}^{d}$. There, the matrix entries are Bernoulli distributed
with probability $F(x_{i},x_{j})$. Using QGT representation as Hessian and
applying the parameter-shift rule, we can see that the entries of $g_{ij}$ are
Poisson binomial distributed [61] with probabilities
$[F^{(++)}_{ij},1-F^{(+-)}_{ij},1-F^{(-+)}_{ij},F^{(--)}_{ij}]$ over a shifted
support $[0,1,2,3,4]\rightarrow[-2,-1,0,1,2]$. Since the of this distribution
are independent of the number of circuit parameters, it can be shown analogous
to Ref. [60] that

$\mathbb{E}[|(\Delta g)_{ij}|^{2}]=\mathcal{O}\left(\frac{1}{N}\right)\text{
and }\mathbb{E}[|(\Delta
g)_{ij}|^{4}]=\mathcal{O}\left(\frac{1}{N^{2}}\right),$ (35)

which leads to a total bound of

$\mathbb{E}[\|\Delta g\|_{2}]=\mathcal{O}\left(\sqrt{\frac{d}{N}}\right).$
(36)

The bound on $\|\Delta\bm{b}\|_{2}$ does not need to be tighter than $\|\Delta
g\|_{2}\|\dot{\bm{\theta}}\|_{2}$, which is straightforward to achieve via the
sampling error. Using the product rule we have

$\displaystyle|\Delta b_{i}|$
$\displaystyle=|\tilde{b}_{i}-b_{i}|=\frac{|\tilde{E}^{(+)}_{i}-\tilde{E}^{(-)}_{i}-E^{(+)}_{i}+E^{(-)}_{i}|}{2}$
(37)
$\displaystyle\leq\frac{|\tilde{E}^{(+)}_{i}-E^{(+)}_{i}|+|\tilde{E}^{(-)}_{i}-E^{(-)}_{i}|}{2}$
$\displaystyle=\mathcal{O}\left(\frac{\sqrt{\mathrm{Var}(E^{(+)}_{i})}+\sqrt{\mathrm{Var}(E^{(-)}_{i})}}{2\sqrt{N}}\right).$

The variance of any state $\ket{\psi}$ can be upper bounded by

$\mathrm{Var}(E)=\braket{\psi}{H^{2}}{\psi}-E^{2}\leq\braket{\psi}{H^{2}}{\psi}\leq
E_{\text{max}}^{2}.$ (38)

Summing over all gradient elements we obtain

$\|\Delta\bm{b}_{\mathrm{max}}\|_{2}=\mathcal{O}\left(\frac{\sqrt{d}E_{\text{max}}}{\sqrt{N}}\right).$
(39)

#### Final bound

Plugging the bounds in the previous paragraphs into Eq. (27) and then into Eq.
(25), we obtain the final bound of

$\varepsilon_{S}\leq\mathcal{O}\left(\frac{d^{3/2}E_{\text{max}}\Delta_{t}}{\delta_{c}^{2}\sqrt{N}}\right).$
(40)

The same asymptotic bound can be derived by performing a moment expansion on
the expectation $\mathbb{E}[\dot{\bm{\theta}}-\tilde{\dot{\bm{\theta}}}]$.

As an example we investigate a simple product-state model, which allows to
show the tightness of several of the above bounds. We look a the first
timestep of the $n$-qubit Hamiltonian $H=\sum_{i=1}^{n}Z_{i}$ with an ansatz
that consists of a single layer of Pauli-$Y$ rotations, each with an
individual parameter. The initial state is $\ket{+}^{\otimes n}$ which is
prepared by setting each of the parameters to $\pi/2$. Each expectation value
is computed with $N=1000$ measurements and we use a regularization of
$\delta_{c}=10^{-2}$. We then vary the number of qubits from $n=2$ to $10$ and
measure the error term contributions over 10 averages, since $\Delta g$ and
$\Delta\bm{b}$ are random variables.

The QGT is measure of the correlation between the parameter derivatives and as
there is not light-cone connecting any two parameterized gates in the product
state ansatz, the QGT is a diagonal matrix. Its norm is therefore
$\|g\|_{2}=1/4$ for any system size. With this restriction, we observe in Fig.
9(a) that the bound on the Bures metric in Eq. (25) is tight as
$\varepsilon_{S}\propto\|\Delta\dot{\bm{\theta}}\|_{2}$. While all bounds are
obeyed, we observe that in particular the bound on $\|\dot{\bm{\theta}}\|_{2}$
is loose, since the bound in Eq. (33) scales with
$\sqrt{d}E_{\text{max}}\propto d^{1.5}$, but we only observe a $d^{0.5}$
scaling. This bound could potentially be further improved by taking into
account that the magnitude of the update step is bounded by the change induced
of the evolution operator $\exp(-\Delta_{t}H)$, which is independent of the
number of parameters $d$.

Figure 9: Scaling of different error bound contributions for a product state
setting. The labels include the scaling with number of parameters, i.e.
$\propto\alpha$ indicates a scaling with $d^{\alpha}$, where $d$ is the number
of parameters.

###  D.2 DualQTE

Assume we require $K$ steps to converge. Then the error in the update step
$\bm{\delta\theta}$ is

$\displaystyle\|\Delta(\bm{\delta\theta})\|_{2}=\|\Delta(\bm{\delta\theta}^{(K)})\|_{2}$
$\displaystyle=\|\widetilde{\bm{\delta\theta}}^{(K-1)}-\eta\widetilde{\bm{\nabla}\mathcal{L}}(\widetilde{\bm{\delta\theta}}^{(K-1)})-\bm{\delta\theta}^{(K-1)}+\eta\bm{\nabla}\mathcal{L}(\bm{\delta\theta}^{(K-1)})\|_{2}$
(41)
$\displaystyle\leq\|\Delta(\bm{\delta\theta}^{(K-1)})\|_{2}+\eta\|\widetilde{\bm{\nabla}\mathcal{L}}(\widetilde{\bm{\delta\theta}}^{(K-1)})-\bm{\nabla}\mathcal{L}(\bm{\delta\theta}^{(K-1)})\|_{2}$
$\displaystyle\leq\|\Delta(\bm{\delta\theta}^{(K-1)})\|_{2}+\eta\|\Delta(\bm{\nabla}\mathcal{L})_{\mathrm{max}}\|_{2}$
$\displaystyle\leq\eta
K\|\Delta(\bm{\nabla}\mathcal{L})_{\mathrm{max}}\|_{2},$

where we used that the error at the initial point is zero,
$\|\Delta(\bm{\delta\theta}^{(0)}\|_{2}=0$. The error in the loss function
gradient can then be written as

$\displaystyle\|\Delta(\bm{\nabla}\mathcal{L}(\bm{\delta\theta}))\|_{2}$
$\displaystyle=\left\|\frac{\Delta(\bm{\nabla}F(\bm{\theta},\bm{\theta}+\bm{\delta\theta}))}{2}+\delta\tau\Delta\bm{b}(\bm{\theta})\right\|_{2}$
(42)
$\displaystyle\leq\frac{\|\Delta(\bm{\nabla}F(\bm{\theta},\bm{\theta}+\bm{\delta\theta}))\|}{2}+\delta\tau\|\Delta\bm{b}(\bm{\theta})\|_{2}$
$\displaystyle\leq\frac{\|\Delta(\bm{\nabla}F)_{\mathrm{max}}\|}{2}+\delta\tau\|\Delta\bm{b}_{\mathrm{max}}\|_{2},$

where
$\Delta(\bm{\nabla}F(\bm{\theta},\bm{\theta}+\bm{\delta\theta}))=\widetilde{\bm{\nabla}F}(\bm{\theta},\bm{\theta}+\bm{\delta\theta})-\bm{\nabla}F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})$
and $\|\Delta(\bm{\nabla}F)_{\mathrm{max}}\|_{2}$ is an upper bound on the
maximum fidelity gradient error for any parameter.

The error in the gradient of $F$ can be derived via the parameter-shift rule,
as

$\displaystyle|\Delta\partial_{i}F|$ $\displaystyle=\frac{\Delta
F^{(+)}_{i}-\Delta F^{(-)}_{i}}{2}$ (43)
$\displaystyle=\mathcal{O}\left(\sqrt{\frac{\mathrm{Var}(F^{(+)}_{i})}{N}}+\sqrt{\frac{\mathrm{Var}(F^{(-)}_{i})}{N}}\right)$
$\displaystyle=\mathcal{O}\left(\frac{1}{\sqrt{N}}\right),$

where we used that the variance of the fidelity can be bounded for any state
$\ket{\psi}$ as

$\displaystyle\mathrm{Var}(F)=\braket{\psi}{P_{0}^{2}}{\psi}-\braket{\psi}{P_{0}}{\psi}^{2}=\braket{\psi}{P_{0}}{\psi}-\braket{\psi}{P_{0}}{\psi}^{2}=F(1-F)\leq\frac{1}{4}.$
(44)

Hence the total error of the fidelity gradient in $\ell_{2}$ norm is

$\|\Delta(\bm{\nabla}F)_{\mathrm{max}}\|_{2}=\mathcal{O}\left(\sqrt{\frac{d}{N}}\right).$
(45)

The bound on $\|\Delta\bm{b}_{\mathrm{max}}\|_{2}$ is already derived in the
previous subsection, which gives then a total of

$\|\Delta(\bm{\delta\theta})\|_{2}=\mathcal{O}\left(\frac{\sqrt{d}K(1+\delta\tau
E_{\mathrm{max}})}{\sqrt{N}}\right).$ (46)

Using the definition $\dot{\bm{\theta}}=\bm{\delta\theta}/\delta\tau$ we then
obtain

$\varepsilon_{S}\leq\Delta_{t}\sqrt{\lambda_{\mathrm{max}}}\frac{\|\Delta(\bm{\delta\theta})\|_{2}}{\delta\tau}=\mathcal{O}\left(\sqrt{\frac{\lambda_{\mathrm{max}}d}{N}}\frac{\Delta_{t}K(1+\delta\tau
E_{\mathrm{max}})}{\delta\tau}\right).$ (47)

##  Appendix E Imaginary-time evolution of the Heisenberg model

###  E.1 Circuit diagram

The circuit used as variational model is schematically presented in Fig. 10.
Each Pauli rotation gate has an independent parameter and the dotted box is
repeated several times. For $r$ repetitions and $n$ qubits the total number of
tunable parameters is thus $2n(r+1)$. The CNOT entangling gates are arranged
in a pairwise manner to minimize the total depth to 2 per entangling layer.

Figure 10: The hardware efficient ansatz for the imaginary-time evolution
experiments. In the QMETTS experiments the Pauli-$Y$ rotations is replaced by
Pauli-$X$ rotations, if the evolution starts in the $Y$ basis states $\ket{\pm
i}$.

###  E.2 Termination and warmstarting

Termination criteria for gradient descent algorithms are typically defined as
achieving a minimal threshold in the difference in the loss function between
update steps or in the gradient norm. However, if only noisy readout of the
loss function is available these criteria become unreliable as the noise in
the evaluation might prevent the termination criterion to be fulfilled even
though the algorithm converged.

One possible resolution would be to consider a moving average over a past
batch of iterations. However, depending on the level of noise, this could
require a large batchsize and therefore many iterations until the termination
can be checked. Since the dual time evolution only has to compute small
corrections, if small timesteps are performed and the optimization are
warmstarted, we only expect a few iterations and a moving average is not a
resource-efficient solution. Therefore we use a heuristic where the first
optimization uses a large number of steps and the subsequent ones perform a
fixed number of few iterations.

To demonstrate the effectiveness of warmstarting and to calibrate the number
of required steps for noisy evaluations we investigate the dual time evolution
in an ideal setting with exact statevector simulations and no finite-sampling
statistics. First, we perform the time evolution for a Heisenberg Hamiltonian
with periodic boundary conditions with $n=12$ sites, $J=1/4$, $g=-1$ and the
initial state $\ket{+}^{\otimes n}$. As circuit model we use the hardware
efficient circuit from Fig. 10 with $r=6$ repetitions and optimize the update
step with a gradient descent routine with a fixed learning rate of $\eta=0.1$.
In each timestep we iterate until the change in loss function
$\Delta\mathcal{L}$ is below the threshold of $10^{-4}\Delta_{t}=10^{-6}$. The
results are presented in Fig. 11(a), and we observe that warmstarting
drastically reduces the number of required iterations until the convergence
criterion is reached.

In a second experiment we analyze how the required number of optimization
steps scales with the system size. In Fig. 11(b) we repeat the above
experiment for $n=3$ to 12 spins and track the number of steps in the first
iteration and the mean and standard deviation of the warmstarted iterations.
We see that the number of steps scales sublinearly in the number of parameters
$d$ and is almost constant for the warmstarted iterations.

Figure 11: (a) The number of iterations required per timestep until
convergence is reached with different initialization techniques. (b) The
number of iterations for different numbers of qubits and the first iteration
and warmstarted iterations. The warmstarted points show mean and standard
deviation of the number of iteration of all steps after the first.

###  E.3 Resource requirements for the dual time evolution

This section shows the detailed VarQITE and DualQITE settings for the resource
estimations in Sec. III. For VarQITE we only varied the number of shots and in
the dual method we additionally allowed to vary the number of iterations in
the optimization in each time step. Especially DualQITE has a lot of
additional degrees of freedom that could be optimized, such as the kind of
optimizer, in addition to settings shared with VarQITE, such as timestep size.

Table 1 shows the settings for VarQITE and DualQITE for the resource
estimation in Fig. 2(b) and Table 2 the settings for the scaling with system
size in Fig. 3.

$I_{B}$ | shots | $N$
---|---|---
1.601 | 100 | $\sim 10^{8}$
0.558 | 1024 | $\sim 10^{9}$
0.149 | 8192 | $\sim 8\cdot 10^{9}$
(a) Settings for VarQITE.

$I_{B}$ | shots | $K_{0}$ | $K_{>0}$ | $N$
---|---|---|---|---
0.937 | 100 | 100 | 10 | $\sim 2.5\cdot 10^{7}$
0.735 | 100 | 200 | 20 | $\sim 5\cdot 10^{7}$
0.305 | 1024 | 100 | 10 | $\sim 2.5\cdot 10^{8}$
0.236 | 1024 | 200 | 20 | $\sim 5\cdot 10^{8}$
0.153 | 2048 | 250 | 25 | $\sim 10^{9}$
(b) Settings for DualQITE.

Table 1: Detailed settings for the resource comparison of VarQITE and DualQITE
at fixed number of qubits $n=12$: the achieved Bures distance $D_{B}$, the
number of shots per circuit and the total number of measurements $N$. The dual
method additionally shows the number of iterations $K_{0}$ in the first
optimization and $K_{>0}$ in the subsequent, warmstarted optimizations. Each
optimization used gradient descent with a learning rate of $\eta=0.1$.

$n$ | shots | $N$
---|---|---
4 | 500 | $4.2\cdot 10^{7}$
6 | 1500 | $4.2\cdot 10^{8}$
8 | 2500 | $1.2\cdot 10^{9}$
10 | 6000 | $6.7\cdot 10^{9}$
12 | 8000 | $1.3\cdot 10^{10}$
(a) Settings for VarQITE.

n | shots | $K_{0}$ | $K_{>0}$ | $\eta$ | $N$
---|---|---|---|---|---
4 | 500 | 100 | 15 | 0.07 | $8.8\cdot 10^{7}$
6 | 600 | 200 | 25 | 0.07 | $3.3\cdot 10^{8}$
8 | 1000 | 100 | 20 | 0.1 | $6\cdot 10^{8}$
10 | 1500 | 200 | 25 | 0.12 | $1.7\cdot 10^{9}$
12 | 2500 | 200 | 25 | 0.1 | $3.5\cdot 10^{9}$
14 | 3000 | 250 | 25 | 0.12 | $4.9\cdot 10^{9}$
(b) Settings for DualQITE.

Table 2: Algorithm settings for the size scaling experiments of VarQITE and
DualQITE.

###  E.4 Gradient benchmark

In this section, we measure how the norm of the loss function gradient,
defined as

$\bm{\nabla}_{\bm{\delta\theta}}\mathcal{L}(\bm{\theta})=-\frac{\bm{\nabla}_{\bm{\delta\theta}}F(\bm{\theta},\bm{\theta}+\bm{\delta\theta})}{2}-\delta\tau\cdot\bm{b}(\bm{\theta})$

scales with the number of qubits $n$ in the Heisenberg Hamiltonian. This
Hamiltonian is 2-local and since the ansatz depth grows logarithmic with the
number of qubits we do not expect exponentially vanishing gradients for the
evolution gradient $b$ [42]. Further, as discussed in the main text, the
initial point of each timestep $\bm{\delta\theta}$ is close to $\bm{0}$, which
ensures the circuit required to measure to fidelity gradient is close to the
identity and we do not expect to encounter a barren plateau [48].

Since both parts of the loss functions are not in a barren plateau setting, we
expect the loss function to be measurable efficiently. In Fig. 12, we measure
the $\ell_{2}$ norm of both the evolution gradient and fidelity gradients and
we find that neither gradient vanishes exponentially. Instead, the evolution
gradient increases with system size, which reflects the extensiveness of
energy in the Heisenberg model. Since we perform imaginary time evolution the
system converges towards the ground state and we expect the energy gradients
to vanish, once converged. In this case, we do not infer a barren plateau as
the gradient norm does not systematically decrease faster for larger systems.
Similarly, the fidelity gradients are expected to decay as the optimal
parameter update $\bm{\delta\theta}$ is found. Note that for system sizes of
12 qubits the exponential decay of the gradients is typically clearly visible
[43].

Figure 12: The $\ell_{2}$ norms for increasing number of qubits $n$ of the
imaginary evolution gradient, $\|\bm{b}\|_{2}$, during the entire evolution
and the norm of the fidelity gradients $\|\bm{\nabla}F^{(t)}_{k}\|_{2}$ at
selected times $t$, where $k$ indicates the iteration in the optimization
within the timestep. Lengths for the fidelity gradients differ as the
optimization were performed using a different number of steps, see Table 2.

##  Appendix F Real-time evolution of the Heisenberg model

###  F.1 Circuit diagram

In the real-time evolution of the Heisenberg model we use the circuit sketched
in Fig. 13, which is the same model used in Ref. [20]. The dotted box is
repeated three times and the rotation layer alternates between Pauli-$X$
rotations (starting from the first layer) and Pauli-$Y$ rotations.

Figure 13: The circuit model used for the real-time evolution experiments.