# Multi-reference Quantum Davidson Algorithm for Quantum Dynamics

Noah Berthusen  Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, Maryland 20742, USA  Faisal Alam
Department of Physics, University of Illinois at Urbana-Champaign, IL 61801,
USA  Yu Zhang<EMAIL_ADDRESS>Theoretical Division, Los Alamos National
Laboratory, Los Alamos, NM 87545, USA

###### Abstract

Simulating quantum systems is one of the most promising tasks where quantum
computing can potentially outperform classical computing. However, the
robustness needed for reliable simulations of medium to large systems is
beyond the reach of existing quantum devices. To address this, Quantum Krylov
Subspace (QKS) methods have been developed, enhancing the ability to perform
accelerated simulations on noisy intermediate-scale quantum computers. In this
study, we introduce and evaluate two QKS methods derived from the QDavidson
algorithm, a novel approach for determining the ground and excited states of
many-body systems. Unlike other QKS methods that pre-generate the Krylov
subspace through real- or imaginary-time evolution, QDavidson iteratively adds
basis vectors into the Krylov subspace. This iterative process enables faster
convergence with fewer iterations and necessitates shallower circuit depths,
marking a significant advancement in the field of quantum simulation.

##  I Introduction

The time evolution of quantum systems is a central problem in condensed matter
physics and quantum chemistry [1, 2, 3, 4]. In particular, the exploration of
strongly correlated quantum systems driven out of equilibrium is essential for
uncovering new physics [5, 6, 7], which has garnered significant research
interest in recent days. This surge in interest requires efficient numerical
modeling capabilities to understand and harness the rich physical phenomena
inherent in these systems, including non-equilibrium topological phases [8,
9], time crystalline order [10], and high harmonic generations [11]. However,
the intrinsic complexity of these problems, characterized by strong
correlations and non-equilibrium, makes the simulations intractable on
classical computers because representations of quantum states and operators
grow exponentially with system size.

Quantum computers were conceived to efficiently and precisely solve such
problems. It is known that the problem of quantum dynamics belongs in the
complexity class, Bounded-error Quantum Polynomial time (BQP) [12], meaning it
can be solved on a quantum computer with a circuit depth that scales
polynomially with system size and evolution time, promising a significant
quantum advantage in efficiency over classical computing methods. However,
near-term quantum computers have limited coherence times and suffer from gate
infidelities, severely limiting the depth of the circuit they can execute
reliably. This restricts the application of quantum algorithms for dynamics to
short evolution times [13]. While foundational in the realm of quantum
simulation, the conventional Trotter methods encounter significant limitations
when applied to the complex dynamics of driven quantum systems. These methods,
hinging on the Trotter-Suzuki decomposition [14], discretize the continuous
evolution of a quantum system into a series of small, manageable steps.
However, this approach inherently assumes a slowly varying Hamiltonian, which
does not align with the rapidly changing dynamics of driven systems. In cases
where the Hamiltonian varies significantly over time, as is common in driven
systems, the Trotter approximation loses its fidelity, leading to large
cumulative errors.

In recent years, a suite of algorithms has emerged under the umbrella of fast
forwarding [15, 16, 17, 18] that attempt to overcome this limitation and
perform evolution to arbitrary times with constant depth circuits. While these
algorithms are not universally applicable to generic Hamiltonians [19], fast
forwarding allows us to effectively exploit currently available quantum
resources for many relevant systems of interest. In general, fast forwarding
refers to approximating the evolution operator for a short time $t$ as [15]:

$e^{-iHt}\approx WDW^{\dagger},$ (1)

where $D$ and $W$ are quantum circuits of polynomial depth, with $D$
consisting of mutually commuting gates. This condition ensures that $D^{r}$
has the same circuit depth as $D$ and allows us to approximate the evolution
operator for long times, $T=rt$ as $e^{-iHT}\approx WD^{r}W^{\dagger}$. The
challenge lies in identifying the circuits $W$ and $D$ for a given
Hamiltonian. Early work on fast-forwarding proposed building the circuits from
parametrized gates, which would then be optimized with a quantum-classical
feedback loop to satisfy Eq. (1), akin to the Variational Quantum Eigensolver
(VQE) technique [20, 21, 22]. However, achieving good optimization requires
executing a substantial number of quantum circuits. More recent efforts have
explored using Krylov subspace (QKS) methods [23, 24, 25, 26, 27] to determine
the $W$ and $D$ circuits that approximate the evolution operator well within a
small subspace of the larger Hilbert space.

In this paper, we propose enhancements to Krylov subspace methods by utilizing
the QDavidson algorithm. Initially developed to identify the ground and
excited states of quantum chemistry Hamiltonians, QDavidson enables the
intelligent selection of additions to our Krylov subspace, thereby reducing
quantum resource demands for fast-forwarding. A significant benefit of the
Krylov subspace approach is its adaptability to open quantum dynamics, wherein
the quantum system evolves under the influence of a bath. Our study further
investigates this generalization, highlighting the practicality and
versatility of the QDavidson-enhanced Krylov subspace methods in addressing
complex quantum dynamics.

This paper is structured as follows: We start with detailed descriptions of
both the closed and open quantum versions of the algorithms we aim to
benchmark in Section II. Section III presents our numerical experiments. The
paper concludes with discussing our findings and suggestions for future
research in Section IV.

##  II Quantum Krylov Subspace Methods

In this section, we delve into the intricacies of the QKS methods,
particularly their multi-reference counterparts, outlining their fundamental
principles. Multi-reference methods have played a pivotal role in advancing
the accuracy and applicability of quantum chemistry simulations [28]. Unlike
single-reference methods that are limited by the description of a single
electronic configuration, multi-reference methods consider multiple
configurations simultaneously, thereby capturing complex electronic
correlation effects essential for describing a wide range of chemical
phenomena. These methods are particularly advantageous in handling systems
with strong electron correlation, open-shell species, and systems undergoing
bond-breaking and forming processes. We briefly introduce various strategies
for constructing the Krylov subspace, each with distinct circuit complexities,
providing a comprehensive understanding of their respective advantages,
limitations, and implications for computational efficiency and accuracy.
Following this, we discuss the application of the multi-reference Quantum
Krylov subspace in fast-forwarding the quantum dynamics of both closed and
open systems.

###  II.1 Multi-reference quantum Krylov subspace methods

Within the QKS framework, the low-lying states of a system are represented as
a linear combination within the Krylov subspace. This is mathematically
expressed as,

$\ket{\psi(t)}\approx\ket{\psi_{K}(t)}=\sum^{M}_{m=1}\sum^{R}_{r=1}\bm{c}_{mr}(t)\ket{\chi_{mr}}$
(2)

where $\bm{c}(t)\in\mathbb{C}^{RM}$ are the corresponding coefficients
associated with each Krylov subspace $\\{\ket{\chi_{mr}}\\}$. The order-$M$
Krylov subspace for a single reference state is generated using real-time
evolution,

$\mathcal{K}_{M}=\text{span}\\{\ket{r},e^{-iH\tau}\ket{r},...,e^{-iH(M-1)\tau}\ket{r}\\}.$
(3)

In practice, we must approximate $e^{-iH\tau}$ using a Trotter decomposition
of the exact time evolution. Additional reference states are iteratively added
to the Krylov subspace until some convergence criteria are met; this could be
when some physical observable or dynamical property is estimated to be within
some tolerance, $\epsilon$. In Ref. [16], the authors propose a simple
selection process that samples the time-evolved reference state
$\ket{\chi_{M-1,R-1}}=e^{-iH(M-1)\tau}\ket{r_{R-1}}$ by performing
measurements in the computational basis on the quantum computer. The
measurement outcomes will follow the probability distribution
$p(x)=|\bra{x}e^{-iH(M-1)\tau}\ket{r_{R-1}}|^{2}$. The bitstring with the
largest observed probability has its order-$M$ Krylov subspace, Eq. (3), added
to the full Krylov subspace.

###  II.2 Quantum Davidson Algorithm

Unlike the previous method, the Krylov subspace is not pre-generated using
time evolution. Instead, the QDavidson algorithm iteratively grows the Krylov
subspace in such a way that it stays close to the exact eigenspace. Doing this
has the advantage of faster convergence and fewer reference states required at
the cost of more circuit evaluations. Previously, the QDavidson algorithm has
been used to express the ground and excited states of the Hamiltonian $H$ as a
linear combination of reference states in the Krylov subspace [27],

$\ket{\psi_{i}}=\sum_{k=1}^{M}\bm{v}_{ki}\ket{\chi_{k}}.$ (4)

in much the same way as the original single-reference Quantum Davidson
algorithm [27]. The intuition is that if these states are well expressed, then
time evolutions (within the eigenspace found by the algorithm) will be
similarly precise. An iteration of the QDavidson algorithm begins by computing
on the quantum computer the subspace matrices $D$ and $E$, whose matrix
elements are defined as:

$D_{k\ell}=\bra{\chi_{k}}H\ket{\chi_{\ell}}\ \text{and}\
E_{k\ell}=\braket{\chi_{k}}{\chi_{\ell}}.$ (5)

This can be done using repeated SWAP tests [29], Hadamard tests [30], or the
MFE protocol [16]. We now solve the generalized eigenvalue problem
$D\bm{v}=\lambda E\bm{v}$ to get the approximate eigenstates and associated
energies. Afterward, the residue of each of the computed eigenstates can be
evaluated,

$\ket{R_{i}}=H\ket{\psi_{i}}-\lambda_{i}\ket{\psi_{i}}=(H-\lambda_{i})\sum_{k=1}^{M}\bm{v}_{ik}\ket{\chi_{k}}.$
(6)

The norm of the residue is then measured on the quantum computer in a similar
manner to the $D$ matrix,

$\displaystyle|\ket{R_{i}}|$
$\displaystyle=\braket{\bm{v_{i}}}{(H-\lambda_{i})^{2}}{\bm{v_{i}}}$ (7)
$\displaystyle=\sum_{k,\ell=1}^{M}\bm{v}_{ki}^{*}\bm{v}_{\ell
i}\braket{\chi_{k}}{(H-\lambda_{i})^{2}}{\chi_{\ell}}$ (8)

Consequently, the correction vector ($\ket{\delta_{i}}$) and the associated
pre-conditioned correction vector ($\ket{\delta^{\prime}_{i}}$) can be derived
as

$\displaystyle\ket{\delta_{i}}$
$\displaystyle=e^{-\Delta\tau(\hat{H}-E_{I})}\ket{\Psi_{I}},$ (9)
$\displaystyle\ket{\delta_{i}^{\prime}}$
$\displaystyle=\ket{\delta_{I}}-\sum_{KJ}\ket{\Psi_{K}}(S^{-1})_{KJ}\bra{\Psi_{J}}\delta_{I}\rangle.$
(10)

If $|\ket{\delta^{\prime}}|<\epsilon$, then the correction vector
$\ket{\delta}$ cannot be expressed as a linear combination of the reference
states in the current Krylov subspace. As such, $\ket{\delta}$ is then added
as a new reference state after mapping the non-unitary operator in Eq. (9)
into unitaries [24]. Doing this for each of the current reference states
completes a full iteration of the algorithm. Iterations are continued until
the resulting simulation is sufficiently accurate for some desired physical
observable and target final time.

With the subspace now constructed, we can estimate the time dynamics of the
system by expressing it as a linear combination of reference states, as we did
in Eq. (2). However, in this case, the coefficients are the expansion of the
Wavefunction in the eigenstate computed from the Quantum Davidson algorithm,

$\ket{\psi(t)}\approx\ket{\psi_{K}(t)}=\sum^{M}_{i=1}\bm{c}_{i}(t)\ket{\chi_{i}}.$
(11)

###  II.3 Multi-reference Quantum Davidson method

A natural extension of the QDavidson algorithm involves integrating with the
multi-reference Krylov method, incorporating time-evolved reference states
into the subspace. After a certain number of iterations, QDavidson yields a
progressively expanded set of reference states $\\{\ket{\chi_{i}}\\}$. Each
reference state then forms the basis for constructing the order-$M$ Krylov
subspace, following the process described in Eq. (3). The accuracy of the
fast-forwarding process is evaluated to determine if it meets the desired
level of precision. If the accuracy is lacking, an additional QDavidson
iteration is undertaken, utilizing the existing reference states
$\\{\ket{\chi_{i}}\\}$. Compared to the conventional multi-reference Krylov
approach, the distinct aspect of this method lies in the selection process of
the reference states. By adopting this strategy, we aim to leverage the
advantageous features of both methodologies, thus enhancing the efficacy and
precision of the fast-forwarding process in quantum dynamics simulations.

###  II.4 Fast Forwarding for closed-system dynamics

Once the Krylov subspace is constructed using any of the mentioned algorithms,
the subsequent steps in the fast-forwarding process remain the same. With the
subspace matrices $D$ and $E$, one can address a generalized eigenvalue
problem to represent the matrix exponential $e^{-iHt}$ in the eigenbasis of
the Krylov subspace [31]. Alternatively, the quantum subspace Schrödinger
equation can be solved,

$iS\partial_{t}\bm{c}(t)=D\bm{c}(t),$ (12)

which has the following analytical solution:

$\bm{c}(t)=e^{-iS^{-1}Dt}\bm{c}(0).$ (13)

where $\bm{c}(0)$ are the initial expansion coefficients. When the initial
state is one of the Krylov basis states $\ket{\chi_{i}}$, then $\bm{c}(0)$ is
the all zeros vector with a one in the $i$th position. Generally, the subspace
overlap matrix $S$ can be poorly conditioned and might not have a well-defined
inverse. To remedy this, we use the singular value decomposition of $E$ and
zero out small singular values before taking the inverse (see [32, 33]). Then,
to fast-forward the system to arbitrary times $t$, we apply the solution of
the subspace Schrödinger equation Eq. (13) to the expansion coefficients
themselves. This yields a time-evolution of the coefficients $\bm{c}(t)$ which
can be used to express an estimate of $\ket{\psi(t)}=e^{-iHt}\ket{\psi(0)}$ as
a linear combination of the reference states in the Krylov subspace as noted
in Eq. (11).

One of the main motivations for performing time dynamics of quantum systems is
to calculate time-dependent observables. In other fast-forwarding methods
where the time-evolved state is available [15, 17], the observables of
interest can be computed directly. However, this is not true for QKS methods
since generating the corresponding quantum state from the time-evolved
expansion coefficients and reference states is not straightforward.
Nevertheless, a general time-evolved observable $O(t)$ can be expressed as
follows:

$\displaystyle O(t)$ $\displaystyle=\braket{\psi_{K}(t)}{O}{\psi_{K}(t)}$
$\displaystyle\approx\sum_{ij}\bm{c}_{i}^{*}(t)\bm{c}_{j}(t)\braket{\chi_{i}}{O}{\chi_{j}}.$
(14)

We have used Eq. (11) for simplicity, but this technique works for any of the
QKS methods, multi-reference or not. The values of
$\braket{\chi_{i}}{O}{\chi_{j}}$ can be calculated on the quantum computer in
the same way as the subspace matrix $D$, and then the results are post-
processed using the time-evolved coefficients to give the observable $O(t)$.

A benefit of using QKS methods for fast-forwarding is that once the $D$ and
$E$ matrices and Eq. (14) are calculated, only classical post-processing is
required to determine $O(t)$ for any time $t$ less than the target final time.
Computing the time-evolved coefficients, Eq. (13), is done entirely on the
classical computer. If a different observable is of interest, it is only
necessary to recompute Eq. (14) on the quantum computer. This is an advantage
over other fast-forwarding methods, where a new time-evolved state must be
prepared for each observable and each time.

###  II.5 Quantum dynamics of open quantum system

This subsection will briefly discuss how the QKS method (particularly Quantum
Davidson) for quantum dynamics developed above can be applied to open quantum
systems within Lindblad and generalized quantum master equation (GQME)
schemes. The main challenge in simulating open quantum system dynamics on
quantum computers is managing the non-Hermitian time propagation [34, 35, 36,
37, 38], as the fundamental operations (gates) of quantum computers are
inherently unitary. To address this challenge, we can leverage the techniques
developed for mapping non-unitary operators into unitary ones. Specifically,
here, we extend the mapping of non-unitary operators into unitaries, as
described in Eq (9), to the non-Hermitian time propagation to simulate the
dynamics of open quantum systems.

The standard approach for deriving the quantum dynamics of an open quantum
system interacting with its environment starts from the von Neumann equation
of the combined system,

$i\hbar\partial\rho_{\text{tot}}(t)=[\hat{H}_{\text{tot}},\rho_{\text{tot}}]\equiv\mathcal{L}\rho_{\text{tot}}(t),$
(15)

where $\hat{H}_{\text{tot}}$ is the total Hamiltonian of the combined system,
$\hat{H}_{\text{tot}}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{SB}$, that includes the
Hamiltonians of the system $\hat{H}_{S}$, bath $\hat{H}_{B}$, and system-bath
interaction $\hat{H}_{SB}$. $\rho_{\text{tot}}$ is the density matrix of the
entire system ($S+B$). In general, the system-bath interaction can be written
as $\hat{H}_{SB}=\sum_{k\alpha}g_{k\alpha}\hat{A}_{k}\hat{B}_{\alpha}$, where
$\hat{A}_{k}$ are the system operators which the bath operators
$\hat{B}_{\alpha}$ couple to, and $g_{k\alpha}$ is the corresponding coupling
strength. Since we are only interested in the dynamics of the system, we can
perform a partial trace over the bath degrees of freedom in Eq. (15) and
thereby obtain a master equation for the motion of the original system density
matrix.

The non-Markovian dynamics of an open quantum system can be described by the
generalized quantum master equation (GQME) (Nakajima-Zwanzig) [39, 40, 41]. We
define a projector operator $\mathcal{P}$ that projects the total density
matrix into the system subspace. The complementary projection superoperator is
denoted as $\mathcal{Q}=1-\mathcal{P}$. $\mathcal{P}$ and $\mathcal{Q}$
satisfy $\mathcal{P}^{2}=\mathcal{P},\mathcal{Q}^{2}=\mathcal{Q}$, and
$[\mathcal{P},\mathcal{Q}]=0$. Hence, projecting Eq. (15) into the subsystem
and following the standard derivation within the Nakajima-Zwanzig formalism
[39, 40, 41], we can get the dynamics of the reduced density matrix (for the
system only) as

$i\hbar\partial\rho_{\text{S}}(t)=\mathcal{L}\rho_{\text{S}}(t)=\mathcal{L}_{0}\rho_{\text{S}}(t)+\mathcal{L}_{m}[\rho_{\text{S}}(t)],$
(16)

where $\mathcal{L}$ is the Liouvillian superoperator, including the overall
system Liouvillian $\mathcal{L}_{0}=\text{Tr}_{B}[\hat{H},\cdot]$ (time-
independent) and memory kernel (time-dependent) $\mathcal{L}_{m}$ parts. Here,
$\text{Tr}_{B}$ denotes the partial trace over the bath degrees of freedom.
The memory kernel, in general, can be written as

$\mathcal{L}_{m}[\rho_{\text{S}}(t)]=-i\hbar\int^{t}_{0}\mathcal{K}(\tau)\rho_{\text{S}}(t-\tau)d\tau.$
(17)

The memory kernel is given by

$\mathcal{K}(\tau)=\frac{1}{\hbar^{2}}\text{Tr}_{B}\left\\{\mathcal{L}e^{-i\mathcal{Q}\mathcal{L}\tau/\hbar}\mathcal{Q}\mathcal{L}\rho_{B}(0)\right\\},$
(18)

where $\mathcal{Q}$ is the projection operator on the bath degrees of freedom,
and $\mathcal{L}$ is the overall system-bath Liouvillian (note: the overall
system-bath Liouvillian is Hermitian).

When memory effect is not essential, we can further simpilfy the quantum
dynamics with Markovian approximation. The most general trace-preserving and
completely positive form of Markovian evolution is the Lindblad master
equation for the reduced density matrix
$\rho_{\text{S}}=\text{Tr}_{B}[\rho_{\text{tot}}]$,

$\displaystyle i\hbar\dot{\rho}_{\text{S}}(t)=$
$\displaystyle\mathcal{L}_{0}\rho_{\text{S}}(t)$
$\displaystyle+i\hbar\sum_{k}\left(C_{k}\rho_{\text{S}}(t)C^{\dagger}_{k}-\frac{1}{2}\\{C^{\dagger}_{k}C_{k},\rho_{\text{S}}(t)\\}\right).$
(19)

where $C_{k}=\sqrt{\gamma_{k}}\hat{A}_{k}$ are the collapse operators and
$\hat{A}_{k}$ are the operators through which the bath couples to the system
in $\hat{H}_{int}$.

Nevertheless, no matter whether it is in Lindblad or GQME formalism, the
resulting time-dependent density matrix can be written in the following
generalized form with a non-unitary propagator, i.e.,

$\rho_{\text{S}}(t)=e^{-i/\hbar\mathcal{L}t}\rho_{\text{S}}(0)=e^{-i/\hbar(\mathcal{L}_{1}+\mathcal{L}_{2})t}\rho_{\text{S}}(0).$
(20)

where $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ describe the unitary and non-
unitary evolution operators, respectively. It should be noted that the
Lindblad master equation is identically the vectorization mapping in the
Liouville space, resulting in

$\displaystyle\frac{d\ket{\rho_{\text{S}}}}{dt}=$
$\displaystyle\Big{[}-iI\otimes\hat{H}+i\hat{H}^{\intercal}\otimes
I+\sum_{k}\Big{(}\overline{L_{k}}\otimes L_{k}$
$\displaystyle-\frac{1}{2}I\otimes(L^{\dagger}_{k}L_{k})-\frac{1}{2}(L^{\intercal}_{k}\overline{L_{k}})\otimes
I\Big{)}\Big{]}\ket{\rho_{\text{S}}}.$ (21)

where the bar indicates entrywise complex conjugation.

However, this requires a doubling of the qubits and an overhead of an ancilla
and controlled operators for evaluating observables. Alternatively, we can map
the density matrix operator as (isomorphism mapping)

$\rho_{\text{S}}=\sum_{p}C_{p}U\ket{p}\bra{p}U^{\dagger}\rightarrow\quad\ket{\rho_{\text{S}}}=\sum_{p}C_{p}U\ket{p}\otimes\bar{U}\ket{p},$
(22)

where $\ket{p}$ are the $n$-qubit computational basis states in the $2^{n}$
possible bit strings. The corresponding Lindblad equation then becomes Eq.
(II.5).

After the isomorphism mapping, the propagator can be again Trotterized, and
each term can be applied term by term.

$\displaystyle e^{-i(I\otimes\hat{H}-i\hat{H}^{\intercal}\otimes
I)\tau}\sum_{p}C_{p}U\ket{p}\otimes\bar{U}\ket{p}$
$\displaystyle=\sum_{p}C_{p}e^{i\hat{H}^{\intercal}\tau}U\ket{p}\otimes
e^{-i\hat{H}\tau}\bar{U}\ket{p}.$ (23)

Where $e^{i\hat{H}^{\intercal}\tau}=\overline{{e^{-i\hat{H}\tau}}}$ for
Hermitian $\hat{H}$. The remaining Lindblad terms in the Trotterized
propagator are of the form $e^{-L^{\intercal}_{k}L_{k}\tau/2}\otimes
e^{-L^{\dagger}_{k}L_{k}\tau/2}$ and $e^{\bar{L}_{k}\otimes L_{k}\tau}$. The
first term preserves the ansatz but is non-unitary, while the second term does
not preserve the ansatz and is non-unitary. The QITE algorithm can then be
used to map the non-unitary terms into unitaries. Suppose the non-unitary term
is denoted as $V_{k}$. For example,
$V_{k}=e^{-L^{\intercal}_{k}L_{k}\tau/2}\otimes
e^{-L^{\dagger}_{k}L_{k}\tau/2}$, we can use the QITE algorithm to find a set
of numbers $D_{p}$ and a Hermitian operator $A$ such that

$\displaystyle V_{k}\sum_{p}C_{p}U\ket{p}\otimes\bar{U}\ket{p}=$
$\displaystyle\sum_{p}(D_{p}+C_{p})e^{iA}U\ket{p}$ $\displaystyle\otimes
e^{-i\bar{A}}\overline{U}\ket{p}+\mathcal{O}(\tau^{2}).$ (24)

where

$D_{q}=-\tau
C_{p}\text{Re}\left[\bra{p}U^{\dagger}L^{\intercal}_{k}\overline{L_{k}}U\ket{p}\right].$
(25)

Like the QITE method, the $A$ operator can be solved using a linear algebra
[42].

##  III Results

Figure 1: Results of performing fast-forwarding of the XYZ model on $N=8$
qubits, setting $J_{X}=J_{Y}=J_{Z}=h=1$, and using a Krylov subspace generated
by the QDavidson algorithm after running for several iterations. We compare
subspaces of dimension 10 (blue), 25 (yellow), and 40 (red) against an exact
time-evolved state (black). The top three panels show the expectation value
$\langle Z_{1}\rangle$ for each state. The bottom panel indicates the state
overlap between the exact time-evolved state with the fast-forwarded states,
$\mathcal{F}=|\braket{\psi(t)}{\psi_{K}(t)}|^{2}$. The gray dashed line
indicates a state overlap with a time evolution using a first-order Trotter
approximation and a comparable number of Trotter steps. The inset is a semilog
plot of infidelity at time $t_{f}=10$ as a function of the Krylov subspace
dimension.

We first benchmark the QDavidson algorithm by simulating the time evolution of
the Heisenberg model,

$H=\sum_{i=1}^{N-1}J_{x}X_{i}X_{i+1}+J_{y}Y_{i}Y_{i+1}+J_{z}Z_{i}Z_{i+1}+\sum_{i=1}^{N}hZ_{i},$
(26)

where $X_{i},Y_{i},Z_{i}$ are Pauli matrices acting on qubit $i$,
$J_{X},J_{Y},J_{Z},h$ are coefficients, and $N$ is the number of qubits. We
initialize the algorithm in the Néel state $\ket{\psi_{i}}=\ket{0101\cdots}$
and simulate to a final time $t_{f}=10$. In this section, all computations
were done on the classical computer, including building the $D$ and $E$
matrices, Eq. (5), and computing the residues. To verify the accuracy of the
fast-forwarded state obtained by Eq. (11), we compute an exact time evolution
of the system, $\ket{\psi(t)}=e^{-iHt}\ket{\psi_{i}}$.

In Fig. 1, we present the results of fast-forwarding the Heisenberg model on
$N=8$ qubits while varying the number of iterations (the dimension of the
Krylov subspace) of the QDavidson algorithm. After performing the fast-
forwarding, one can calculate time-dependent properties of interest by using
Eq. (14). The upper three panels show the computed expectation value $\langle
Z_{1}\rangle$ for each fast-forward state and a black line indicating the true
value. The bottom panel reports the fidelity of each fast-forwarded state with
respect to the exact time-evolved state,
$\mathcal{F}=|\braket{\psi(t)}{\psi_{K}(t)}|^{2}$. Increasing the number of
iterations predictably increases the accuracy of the fast-forwarding, which is
reflected in a more precise match with the true expectation value and higher
fidelity with the true state. The gray dashed line indicates the overlap with
the exact Wavefunction of an evolution done using a first-order Trotter
approximation and 40 Trotter steps. We see that the evolution done with
QDavidson vastly outperforms Trotter, even when using a significantly shorter
circuit. The inset in the bottom panel shows infidelity at time $t_{f}=10$ as
a function of the dimension of the Krylov subspace. We observe a slow
convergence rate before a sharp drop-off. At short times, the fidelity remains
high even when using a less expressive (smaller) subspace; if one is only
interested in the dynamics during this regime, then using fewer iterations and
a smaller Krylov subspace is acceptable.

Figure 2: (a)-(b) Scaling of QKS methods with system size on the Heisenberg
XXX model with $J_{x}=J_{y}=J_{z}=h=1$. (a) The semilog plot required Krylov
dimension to fast-forward to time $t_{f}=10$ with fidelity 0.9 as a function
of system size. Here, we compare the QDavidson algorithm (black markers)
against multi-reference QDavidson (blue markers) and multi-reference Krylov
(yellow markers). Both multi-reference methods use $M=10$ reference states.
(b) The semilog plot of a corresponding number of algorithm iterations to
reach the desired fidelity. (c)-(d) Similar scaling of QKS methods with system
size on the Heisenberg XYZ model with $J_{x}=h=1,J_{y}=2,J_{z}=3$.

To compare the scaling of (multi-reference) QDavidson and multi-reference
Krylov for time dynamics, we simulate systems of up to $N=13$ qubits. We focus
here on the Heisenberg model where the Hamiltonian coefficients are set to be
equal, $J_{x}=J_{y}=J_{z}=h=1$. In Fig. 2(a), we plot the minimum Krylov
dimension required for the fast-forwarded state to have 90% overlap with the
exact Wavefunction at the final time $t_{f}=10$. We compare the standard
QDavidson algorithm (black markers) against multi-reference Krylov with $M=10$
reference states (yellow markers) as well as multi-reference QDavidson again
with $M=10$ reference states (blue markers). For the multi-reference methods,
the order-$M$ Krylov subspace was generated using an exact time evolution; we
later investigated the effect of Trotter error on the methods’ performance.
Regarding the number of states needed to recreate the exact time-evolved state
accurately, QDavidson outperforms the other methods by almost an order of
magnitude. The semilog plot clearly shows sub-linear behavior from the
QDavidson algorithm, indicating sub-exponential scaling of the Krylov
dimension for this Hamiltonian. The same cannot be said about the multi-
reference methods, where the scaling is unclear. Fig. 2(b) displays the
corresponding number of algorithm iterations required to reach the desired 90%
overlap. Here, we see that the multi-reference methods need fewer iterations
than the QDavidson algorithm, although, for large enough system sizes, this
will no longer hold.

Figs. 2(c)-(d) show similar trends for the Heisenberg XYZ model with
$J_{x}=h=1,J_{y}=2,J_{z}=3$. The scaling behavior of the required Krylov
dimension and the number of iterations is consistent with Figs. 2(a)-(b). The
QDavidson algorithm continues to demonstrate sub-exponential scaling in the
Krylov dimension, outperforming the multi-reference methods regarding the
number of states needed to represent the time-evolved state accurately. The
number of iterations required for the multi-reference methods remains lower,
but this difference diminishes with increasing system size. These results
indicate that the QDavidson algorithm is robust and efficient across different
models and parameter settings, providing a significant advantage for
simulating quantum dynamics in larger systems.

Figure 3: Semilog plot comparing the effect of Trotter error on the required
Krylov dimension to fast-forward to time $t_{f}=10$ with a fidelity of 0.9.
Here we again use the Heisenberg XXX model with $J_{x}=J_{y}=J_{z}=h=1$. We
compare the multi-reference Krylov method with $M=10$ reference states where
the order-$M$ Krylov subspace was generated with exact time evolution (solid
yellow) against one where the Krylov subspace was generated with a first-order
Trotter product formula and Trotter step size $\tau=0.1$ (dashed yellow). We
also compare an exact time evolution-generated multi-reference QDavidson
(solid blue) with a Trotter-generated multi-reference QDavidson (dashed blue).

###  III.1 Effect of Trotter error

In the previous section, we focused on the performance of the algorithms by
neglecting the effects of noise and error. Here, we discuss the impact of the
Trotter error on the Krylov subspace dimensions. Fig. 3 compares the required
Krylov dimension to fast-forward to time $t_{f}=10$ with a fidelity of 0.9 for
the Heisenberg XXX model with $J_{x}=J_{y}=J_{z}=h=1$. We compare the multi-
reference Krylov method with $M=10$ reference states, where the order-$M$
Krylov subspace was generated with exact time evolution (solid yellow) against
one where the Krylov subspace was generated with a first-order Trotter product
formula and Trotter step size $\tau=0.1$ (dashed yellow). We also compare an
exact time evolution-generated multi-reference QDavidson (solid blue) with a
Trotter-generated multi-reference QDavidson (dashed blue).

The results indicate that the multi-reference QDavidson algorithm is robust
against Trotter errors, maintaining its efficiency and accuracy even when
approximations are introduced. For this method, the required Krylov dimension
is practically unchanged in the presence of Trotter errors. While still
effective, the multi-reference methods show a more pronounced increase in the
required Krylov dimension when Trotter errors are present.

##  IV Discussion

In this work, we have developed and benchmarked the QDavidson algorithm for
fast-forwarding quantum simulations. Compared with the recently introduced
selected quantum Krylov fast-forwarding (sQKFF) algorithm, we find that the
QDavidson algorithm can achieve comparable accuracies while using a smaller
Krylov subspace. However, this comes at the cost of more iterations and
circuit evaluations on the quantum computer. To address the trade-off between
the Krylov subspace dimension and the number of iterations, we introduced a
multi-reference approach that combines properties of both QDavidson and multi-
reference Krylov methods. We found it to be a compromise, balancing the Krylov
subspace dimension and the number of iterations required. Depending on the
limitations of the available device, the number of reference states can be
tuned to optimize performance.

We also examined the effect of Trotter error on the performance of these
algorithms. Interestingly, we found that Trotter error can benefit the
algorithm by making the resulting states less dependent. The multi-reference
QDavidson algorithm, in particular, proved robust against Trotter errors,
maintaining its efficiency and accuracy even when approximations are
introduced. This robustness highlights the potential of QDavidson for
practical implementations on near-term quantum devices.

Future work on using Krylov subspace methods would benefit from investigating
the performance of these algorithms on real quantum hardware. Understanding
the impact of hardware-specific noise and errors and developing error
mitigation techniques will be crucial for successfully applying these methods
in practical quantum simulations. Additionally, exploring the scalability of
the QDavidson algorithm for larger systems and its integration with other
quantum algorithms could further enhance its utility in solving complex
quantum problems.

## Acknowledgements

We acknowledge support from the US DOE, Office of Science, Basic Energy
Sciences, Chemical Sciences, Geosciences, and Biosciences Division under Triad
National Security, LLC (“Triad”) contract Grant 89233218CNA000001 (FWP:
LANLECF7) and the Laboratory Directed Research and Development (LDRD) program
of Los Alamos National Laboratory (LANL). This research used computational
resources provided by the Institutional Computing (IC) Program and the Darwin
testbed at Los Alamos National Laboratory (LANL), which is funded by the
Computational Systems and Software Environments subprogram of LANL’s Advanced
Simulation and Computing program. LANL is operated by Triad National Security,
LLC, for the National Nuclear Security Administration of the U.S. Department
of Energy (Contract No. 89233218CNA000001).

## Data Availability

The source code and data to generate the figures in the paper are provided
freely at https://github.com/noahberthusen/qdavidson-dynamics.

## References

  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Quantum chemistry in the age of quantum computing, Chem. Rev. 119, 10856 (2019).
  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Quantum computational chemistry, Rev. Mod. Phys. 92 (2020).
  * Ollitrault  _et al._ [2021] P. J. Ollitrault, A. Miessen, and I. Tavernelli, Molecular quantum dynamics: A quantum computing perspective, Acc. Chem. Res. 54, 4229 (2021).
  * Nelson  _et al._ [2020] T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A. E. Roitberg, and S. Tretiak, Non-adiabatic excited-state molecular dynamics: Theory and applications for modeling photophysics in extended molecular materials, Chem. Rev. 120, 2215 (2020).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The hubbard model: A computational perspective, Annu. Rev. Condens. Matter Phys. 13, 275 (2022).
  * Rudner and Lindner [2020] M. S. Rudner and N. H. Lindner, Band structure engineering and non-equilibrium dynamics in floquet topological insulators, Nat. Rev. Phys. 2, 229 (2020).
  * Liu  _et al._ [2018] H. Liu, J.-T. Sun, C. Cheng, F. Liu, and S. Meng, Photoinduced nonequilibrium topological states in strained black phosphorus, Phys. Rev. Lett. 120, 237403 (2018).
  * Zaletel  _et al._ [2023] M. P. Zaletel, M. Lukin, C. Monroe, C. Nayak, F. Wilczek, and N. Y. Yao, Colloquium: Quantum and classical discrete time crystals, Rev. Mod. Phys. 95, 031001 (2023).
  * Gorlach  _et al._ [2020] A. Gorlach, O. Neufeld, N. Rivera, O. Cohen, and I. Kaminer, The quantum-optical nature of high harmonic generation, Nat. Commun. 11, 4598 (2020).
  * [12] A. Childs, Lecture notes on quantum algorithms.
  * Miessen  _et al._ [2023] A. Miessen, P. J. Ollitrault, F. Tacchino, and I. Tavernelli, Quantum algorithms for quantum dynamics, Nature Computational Science 3, 25 (2023).
  * Suzuki [1976] M. 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 (1976).
  * Cîrstoiu  _et al._ [2020] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, Variational fast forwarding for quantum simulation beyond the coherence time, npj Quantum Information 6 (2020).
  * Cortes and Gray [2022] C. L. Cortes and S. K. Gray, Quantum krylov subspace algorithms for ground- and excited-state energy estimation, Phys. Rev. A 105, 022417 (2022).
  * 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, arXiv preprint arXiv:2102.04313 (2021).
  * Atia and Aharonov [2017] Y. Atia and D. Aharonov, Fast-forwarding of hamiltonians and exponentially precise measurements, Nat. Commun. 8, 1572 (2017).
  * Berry  _et al._ [2006] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Efficient quantum algorithms for simulating sparse hamiltonians, Communications in Mathematical Physics 270, 359 (2006).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5, 4213 (2014).
  * Tkachenko  _et al._ [2021] N. V. Tkachenko, J. Sud, Y. Zhang, S. Tretiak, P. M. Anisimov, A. T. Arrasmith, P. J. Coles, L. Cincio, and P. A. Dub, Correlation-informed permutation of qubits for reducing ansatz depth in the variational quantum eigensolver, PRX Quantum 2, 020337 (2021).
  * Zhang  _et al._ [2022] Y. Zhang, L. Cincio, C. F. A. Negre, P. Czarnik, P. J. Coles, P. M. Anisimov, S. M. Mniszewski, S. Tretiak, and P. A. Dub, Variational quantum eigensolver with reduced circuit complexity, npj Quantum Information 8, 96 (2022).
  * Stair  _et al._ [2020] N. H. Stair, R. Huang, and F. A. Evangelista, A multireference quantum krylov algorithm for strongly correlated electrons, J. Chem. Theory Comput. 16, 2236 (2020).
  * 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. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nat. Phys. 16, 205 (2020).
  * Yeter-Aydeniz  _et al._ [2021] K. Yeter-Aydeniz, B. T. Gard, J. Jakowski, S. Majumder, G. S. Barron, G. Siopsis, T. S. Humble, and R. C. Pooser, Benchmarking quantum chemistry computations with variational, imaginary time evolution, and krylov space solver algorithms, Adv. Quantum Technol. 4, 2100012 (2021).
  * Epperly  _et al._ [2022] E. N. Epperly, L. Lin, and Y. Nakatsukasa, A theory of quantum subspace diagonalization, SIAM J. Matrix Anal. Appl. 43, 1263 (2022).
  * Tkachenko  _et al._ [2024] N. V. Tkachenko, L. Cincio, A. I. Boldyrev, S. Tretiak, P. A. Dub, and Y. Zhang, Quantum davidson algorithm for excited states, Quantum Sci. Technol. 9, 035012 (2024).
  * Lischka  _et al._ [2018] H. Lischka, D. Nachtigallová, A. J. A. Aquino, P. G. Szalay, F. Plasser, F. B. C. Machado, and M. Barbatti, Multireference approaches for excited states of molecules, Chem. Rev. 118, 7293 (2018).
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. de Wolf, Quantum fingerprinting, Phys. Rev. Lett. 87, 167902 (2001).
  * Aharonov  _et al._ [2006] D. Aharonov, V. Jones, and Z. Landau, A polynomial quantum algorithm for approximating the jones polynomial, in  _Proceedings of the Thirty-Eighth Annual ACM Symposium on Theory of Computing_, STOC ’06 (Association for Computing Machinery, New York, NY, USA, 2006) pp. 427–436.
  * Lim  _et al._ [2021] K. H. Lim, T. Haug, L. C. Kwek, and K. Bharti, Fast-forwarding with NISQ processors without feedback loop, Quantum Sci. Technol. 7, 015001 (2021).
  * Klymko  _et al._ [2022] K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, J. Moussa, N. Wiebe, W. A. de Jong, and N. M. Tubman, Real-time evolution for ultracompact hamiltonian eigenstates on quantum hardware, PRX Quantum 3, 020323 (2022).
  * Cortes  _et al._ [2022] C. L. Cortes, A. E. DePrince, and S. K. Gray, Fast-forwarding quantum simulation with real-time quantum krylov subspace algorithms, Phys. Rev. A 106 (2022).
  * Wang  _et al._ [2023] Y. Wang, E. Mulvihill, Z. Hu, N. Lyu, S. Shivpuje, Y. Liu, M. B. Soley, E. Geva, V. S. Batista, and S. Kais, Simulating open quantum system dynamics on nisq computers with generalized quantum master equations, J. Chem. Theory Comput. 19, 4851 (2023).
  * Hu  _et al._ [2020] Z. Hu, R. Xia, and S. Kais, A quantum algorithm for evolving open quantum dynamics on quantum computing devices, Sci. Rep. 10, 3301 (2020).
  * Schlimgen  _et al._ [2021] A. W. Schlimgen, K. Head-Marsden, L. M. Sager, P. Narang, and D. A. Mazziotti, Quantum simulation of open quantum systems using a unitary decomposition of operators, Phys. Rev. Lett. 127, 270503 (2021).
  * Cattaneo  _et al._ [2023] M. Cattaneo, M. A. C. Rossi, G. García-Pérez, R. Zambrini, and S. Maniscalco, Quantum simulation of dissipative collective effects on noisy quantum computers, PRX Quantum 4, 010324 (2023).
  * Del Re  _et al._ [2024] L. Del Re, B. Rost, M. Foss-Feig, A. F. Kemper, and J. K. Freericks, Robust measurements of n-point correlation functions of driven-dissipative quantum systems on a digital quantum computer, Phys. Rev. Lett. 132, 100601 (2024).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Nakajima [1958] S. Nakajima, On quantum theory of transport phenomena: Steady diffusion, Prog. Theor. Phys. 20, 948 (1958).
  * Zwanzig [1960] R. Zwanzig, Ensemble method in the theory of irreversibility, J. Chem. Phys. 33, 1338 (1960).
  * Kamakari  _et al._ [2022] H. Kamakari, S.-N. Sun, M. Motta, and A. J. Minnich, Digital quantum simulation of open quantum systems using quantum imaginary–time evolution, PRX Quantum 3, 010320 (2022).