Efficient Measurement-Driven Eigenenergy Estimation with Classical Shadows

Yizhi Shen,1, ∗ Alex Buzali,2, † Hong-Ye Hu,2 Katherine Klymko,3

Daan Camps,3 Susanne F. Yelin,2 and Roel Van Beeumen1
1
Applied Mathematics and Computational Research Division,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
2
Department of Physics, Harvard University, 17 Oxford Street, Cambridge, MA 02138, USA
3
NERSC, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
Quantum algorithms exploiting real-time evolution under a target Hamiltonian have demonstrated
remarkable efficiency in extracting key spectral information. However, the broader potential of these
methods, particularly beyond ground state calculations, is underexplored. In this work, we introduce
the framework of multi-observable dynamic mode decomposition (MODMD), which combines the
observable dynamic mode decomposition, a measurement-driven eigensolver tailored for near-term
arXiv:2409.13691v2 [quant-ph] 13 Nov 2024

implementation, with classical shadow tomography. MODMD leverages random scrambling in the
classical shadow technique to construct, with exponentially reduced resource requirements, a signal
subspace that encodes rich spectral information. Notably, we replace typical Hadamard-test circuits
with a protocol designed to predict low-rank observables, thus marking a new application of classical
shadow tomography for predicting many low-rank observables. We establish theoretical guarantees
on the spectral approximation from MODMD, taking into account distinct sources of error. In the
ideal case, we prove that the spectral error scales as exp(−∆Etmax ), where ∆E is the Hamiltonian
spectral gap and tmax is the maximal simulation time. This analysis provides a rigorous justification
of the rapid convergence observed across simulations. To demonstrate the utility of our framework,
we consider its application to fundamental tasks, such as determining the low-lying, i.e., ground or
excited, energies of representative many-body systems. Our work paves the path for efficient designs
of measurement-driven algorithms on near-term and early fault-tolerant quantum devices.

I. INTRODUCTION multiple eigenenergies may demand substantial quantum

or classical resources.
Quantum algorithms based on real-time evolution [1– In this work, we introduce a real-time evolution frame-
14] have gained increasing popularity due to the unitarity work that flexibly combines real-time quantum evolution
of real-time dynamics and their native implementation on with classical data post-processing to access both ground
quantum hardware. In particular, real-time eigensolvers and lower excited state energies of a quantum many-body
have demonstrated impressive efficacy in extracting the system. Specifically, our framework is efficient as it only
ground state information of physical many-body systems. requires the preparation of a single, simple initial state
This success raises a natural question: can the algorithms on the quantum computer, with overall resource demands
perform highly accurate excited state calculations? comparable to the leading ground state algorithms. Fur-
The real-time approaches were initially conceptualized thermore, our classical post-processing is ansatz-free and
from a subspace diagonalization perspective [1, 2]. While, completely circumvents the challenges of complex opti-
in theory, these subspace methods have the potential to mization landscapes commonly encountered in paramet-
uncover spectral information beyond the ground state, ric methods, particularly as the number of target eigenen-
they typically involve solving ill-conditioned generalized ergies increases [17–19].
eigenvalue problems [15] and can thus suffer from pertur- The main technical tools that we leverage are classi-
bative errors, including standard hardware noise. Alter- cal shadow tomography [20–26] and observable dynamic
natively, recent approaches [7, 8, 10, 11, 16] draw inspi- mode decomposition (ODMD) [9]. We develop a simple,
rations from a signal processing perspective, where real- single-ancilla shadow protocol to collect real-time signals
time evolution is primarily utilized to generate time se- associated with many observables, and utilize ODMD to
ries that resonate with the target eigenmodes. In general, unravel these signals into single-energy modes. Since os-
the signal processing approaches are resilient to pertur- cillating time signals with distinct power spectra are lin-
bative noise due to their robust regularization. On the early independent in the space of all sinusoidal functions,
other hand, they often rely on state preparation such that they form the basis of some signal subspace, which can
the initial state undergoing the Hamiltonian dynamics be systematically expanded to accommodate the desired
should overlap dominantly with the eigenstates of inter- frequencies. Our approach, termed the multi-observable
est. Without such a basic assumption, accurately finding dynamic mode decomposition (MODMD), can hence be
viewed as a unifying framework that enjoys the strength
of both subspace and signal processing algorithms.
Compared to alternative hybrid eigensolvers, our real-
∗ yizhis@lbl.gov; Equal contributions time framework has the following notable advantages: it
† Equal contributions (1) achieves near-exponential convergence of the eigenen-
 2

Hamiltonian Shallow
Simulation Shadows
⟨ϕ0 | O1 | ϕ0(Δt)⟩
⟨ϕ0 | O2 | ϕ0(Δt)⟩

Classical Post-processing
Quantum Measurements

⋮
⟨ϕ0 | OI | ϕ0(Δt)⟩
LS
OI =
⟨ϕ0 | O1 | ϕ0(2Δt)⟩
⟨ϕ0 | O2 | ϕ0(2Δt)⟩

⋮ O1
X′ ∈ ℝd×(K+1) A ∈ ℝd×d X ∈ ℝd×(K+1)
⟨ϕ0 | OI | ϕ0(2Δt)⟩

⋮ ⋮ Spectral
Energy Error Decomposition
⋮ ⋮
Excited State
⟨ϕ0 | O1 | ϕ0(7Δt)⟩
1
⟨ϕ0 | O2 | ϕ0(7Δt)⟩
− argλ˜n ≈ En
⋮ K Δt
⟨ϕ0 | OI | ϕ0(7Δt)⟩
Ground State

FIG. 1: MODMD for eigenenergy estimation. MODMD collects the expectations, Tr ρ(k∆t)ΓO = ℜ ⟨ϕ0 | Oe−iHk∆t |ϕ0 ⟩, with
 
respect to a simple reference state |ϕ0 ⟩. This data can be measured efficiently on a quantum processor through Hamiltonian
simulations combined with shadow tomographic techniques, where the shadow circuits can be shallow with a depth logarithmic
in the number of qubits. MODMD then constructs a pair of block Hankel matrices (X, X′ ), and computes the DMD system
matrix A that adopts a block companion structure. The eigenvalues λ̃n of A converge to the true eigenphases λn as the size of
X and X′ increases. The low-lying energies En are estimated as angles of λ̃n , Ẽn = − ∆t
1
arg(λ̃n ).

ergy estimates surpassing the conventional Fourier limit, In particular, classical shadows enable the simultaneous
(2) provides extensive knowledge of a many-body system prediction of expectations for many observables O of our
including eigenstate properties and dynamical responses, choice. These expectations generate a multivariate time
and (3) significantly saves quantum resources in terms of series, whose characteristic frequencies can be efficiently
the evolution time while showing stability against pertur- extracted by ODMD (see Fig. 1 for a summary).
bative noise. These features make our algorithm promis- From here on, we will use H ⊂ CN to denote the physi-
ing for near-term implementations on current quantum cal Hilbert space of many-body quantum states and L(H)
platforms, such as analog quantum simulators and early to denote the Liouville space, i.e., the space of linear op-
fault-tolerant quantum computers. erators acting on H (we work with finite-dimensional and
The manuscript is organized as follows. In Section II, thus bounded linear operators for simplicity). Additional
we first overview key concepts underlying recent advances essential notations are defined self-consistently in Table I
in real-time eigensolvers. We then establish in Section III and throughout the main text.
the building blocks of the real-time MODMD framework
and detail our core algorithm for eigenenergy and eigen-
state estimation. Theoretical guarantees on its conver- A. Real-time quantum eigensolvers
gence and preliminary error analysis are presented within
Section IV. Finally, we numerically demonstrate our algo-
rithm in Section V by focusing on many-body examples In this subsection, we first highlight the capabilities of
from condensed matter physics and quantum chemistry. real-time evolution in determining the eigenenergies of a
target Hamiltonian H ∈ CN ×N . We consider the spectral
PN −1
decomposition, H = n=0 En |ψn ⟩ ⟨ψn |, of the Hamilto-
nian with ordered energies E0 ≤ E1 ≤ . . . ≤ EN −1 . The
II. PRELIMINARIES
real-time approaches commonly require evaluation of the
expectation value [1–3, 5, 7–13, 16],
In this work, we develop an approach that gives highly
accurate eigenenergy estimates beyond the ground state. N
X −1
2
Specifically, we propose a real-time framework combining ⟨ϕ0 |e−iHt |ϕ0 ⟩ = |αn | e−iEn t , (1)
the observable dynamic mode decomposition (ODMD) n=0
and classical shadow tomography, where we fully leverage
the synergy between these two algorithmic components. where |ϕ0 ⟩ ∈ H is a reference state, αn = ⟨ψn |ϕ0 ⟩ is the
We will show that we can design a simple shadow protocol overlap between the reference state and eigenstate |ψn ⟩,
to predict the real-time expectation value ⟨ϕ0 |Oe−iHt |ϕ0 ⟩ and t = k∆t is typically an integer multiple of some time
for the problem Hamiltonian H and a reference state |ϕ0 ⟩. step ∆t. The expectation can be efficiently sampled from
 3

TABLE I: Glossary of key notations. B. Efficient measurement with classical shadows

Hilbert space H, dimH = N = 2L Classical shadow tomography [20, 21, 33–35] embodies
Liouville space L(H), dimL(H) = N 2 a powerful suite for efficiently measuring expectations of
Hamiltonian in Pauli basis H= M
P many observables simultaneously:
ν=1 κν Pν
{Oi = M I
{Γi }Ii=1 ,
P
ν=1 κi,ν Pi,ν }i=1
i
 
Hermitian observables Tr ρϕ Γi = ⟨ϕ|Γi |ϕ⟩ , (2)
N −1
Hamiltonian eigenstates {|ψn ⟩}n=0
where Γi ∈ L(H) has an efficient representation on a clas-
Pure quantum state ρϕ = |ϕ⟩ ⟨ϕ|
sical computer. Here ρϕ = |ϕ⟩ ⟨ϕ| ∈ L(H) is a pure state,
Multi-observable signal {⟨ϕ0 |Oi e−iHt |ϕ0 ⟩}Ii=1 though a similar trace evaluation applies to mixed states.
Classical shadow dataset {σ̂q = Uq† |bq ⟩⟨bq |Uq }Q
q=1 Classical shadow tomography consists of two key steps:
Systematic error (e.g., shot noise) ϵ1 (1) random quench evolution using U ∈ U from a unitary
Algorithmic error from eigensolver ϵ ensemble U, and (2) computational basis measurement.
Upon each measurement, the quantum state collapses to
a bitstring |bq ⟩. After repetitive experiments, one obtains
a classical shadow dataset, {Uq† |bq ⟩⟨bq |Uq }Q
q=1 , which can
repeated quantum measurements through the Hadamard be viewed as a classical sketch of the quantum state. It is
test [27] or mirror fidelity estimation [5, 28]. known [20] that Q = O(log(I) maxi ∥Γi ∥sh ϵ−2 1 ) shadows
There are several approaches for obtaining energy esti- can predict all the I expectation values given by Eq. (2)
mates from data of the type Eq. (1). Quantum subspace to uncertainty ϵ1 with high probability. Here, the shadow
diagonalization [1–3, 5, 12, 13] approximates the extremal norm, ∥Γ∥sh , depends on both the unitary ensemble U
eigenenergies by forming a projected eigenvalue problem. and the operator Γ. For example, when the ensemble U
In particular, a subspace can be built through successive is the L-qubit random Clifford unitaries, i.e., U = Cl(2L ),
real-time evolutions as new matrix elements with increas- and the operator Γ = Γ† is Hermitian, the shadow norm
ing k are added. With a single reference state |ϕ0 ⟩, these is ∥Γ∥sh = 3Tr[Γ2 ]. This is especially powerful if Γ is low-
methods often struggle to locate excited states, though rank, meaning that its operator rank stays independent
convergence can be potentially accelerated by the prepa- of the system size. Because Tr[Γ2 ] ≤ rank(Γ)∥Γ∥22 , one
ration of multiple reference states with a cost quadratic can use classical shadows to predict exponentially many
in their number. Moreover, we notice that the real-time low-rank expectations simultaneously even for large sys-
states {e−iHk∆t |ϕ0 ⟩}k≥0 lose orthogonality as a subspace tems. In Section III, we demonstrate how to transfer the
basis, which can thereby lead to ill-conditioning and sus- measurement of many expectations {ℜ⟨ϕ0 |Oi e−iHt |ϕ0 ⟩}i
ceptibility to noise. to the task of predicting low-rank Hermitian operators.
Overall, the framework presented within this work is
Alternatively, the signal processing methods [7, 8, 10, distinctive in various crucial aspects. First, our algorithm
11, 16] capitalize on Fourier or harmonic analysis to re- can directly estimate the individual eigenenergies and the
solve the eigenfrequencies or eigenenergies of interest. associated energy gaps. Conversely, resolving the energy
These methods mitigate the impact of noise by minimiz- levels from estimated excitation gaps is nearly infeasible.
ing a customized objective function, which sharpens into Second, instead of the standard usage of classical shadows
robust optima as the number of time steps increases. De- in predicting many local Pauli observables, we introduce
spite the resilience to noise, the optimization landscape a novel application of classical shadows by replacing the
depends critically on the choice of the reference state |ϕ0 ⟩. Hadamard-test-like circuits with those that predict many
For example, let us consider Eq. (1) as a spectral density low-rank observables, expanding their primary utility be-
over the unit circle, where the squared overlap |αn |2 in- yond entanglement witnesses in quantum information sci-
dicates the normalized spectral weight. Techniques such ence. Last, our post-processing scheme yields robust and
as peak finding on the spectral density may not yield a accurate eigenenergy estimates, achieving an exponential
unique solution when the reference state has nearly uni- error reduction in the low-noise regime. This surpasses
form eigenstate overlaps. To ensure accurate energy esti- an algebraic error decay in the conventional Fourier limit.
mations, more sophisticated reference state preparations Although recent work [29] has begun exploring the use
and post-processing designs need to be accounted for. of shadow techniques in spectroscopic calculations of en-
To access both ground and excited state properties, we ergy gaps, our approach distinguishes itself in the sense
aim for a simple, accurate, and robust real-time protocol. we have discussed above.
To this end, we will introduce a quantum signal subspace Table II presents a comprehensive comparison of state-
approach that utilizes signals of the form ⟨ϕ0 |Oe−iHt |ϕ0 ⟩ of-the-art real-time methods for extracting valuable spec-
for general operators O. Importantly, we seek to measure tral information, highlighting the efficiency and accuracy
many such operator expectations simultaneously with a of our eigensolver. In Section III we will outline the basic
cost logarithmic in the number of observables. This can construction of our measurement-driven approach, which
be precisely achieved via classical shadows, which we ex- deploys classical shadows to evaluate the real-time expec-
plore in the following subsection. tations.
 4

Real-time methods Measurement cost Algorithmic convergence Target overlaps System properties
Single-observable signal
O ϵ−2 Õ(ϵ−2


Subspace [1–3, 5, 12, 13] 1 1 ) in low-ϵ1 limit non-vanishing ✓ eigenstates
unstable in high-ϵ1 limit
O ϵ−2 Õ ϵ−1



Signal processing [7, 8, 10, 11] 1 dominant ✗

Multi-observable signal
O log(I)3loc[O] ϵ−2


Shadow spectroscopy [29] 1 uncertain non-vanishing ✓ dynamics
−2

Õ ϵ−1 in some cases


Signal subspace (this work) O log(I)ϵ1 non-vanishing ✓ eigenstates
stable with respect to ϵ1 [30–32] dynamics

TABLE II: Landscape of real-time hybrid eigensolver leveraging quantum measurements of ⟨ϕ0 |Oe−iHt |ϕ0 ⟩ and classical post-
processing schemes. We consider estimating a target set of Hamiltonian eigenergies, {En : n ∈ Neig ⊂ [N ]}, where |Neig | = Neig
for Neig ≪ N . On quantum hardware, real-time evolution of a reference state |ϕ0 ⟩ is performed and different expectation values
are measured with an uncertainty of ϵ1 . The classical post-processing then extracts from these measurements our eigenenergy
estimates with an algorithmic error of ϵ. Here Õ(·) denotes asymptotic upper bound with multiplicative polylogarithmic factors
neglected.

III. MODMD FRAMEWORK We will show that one can evaluate all the I expectation
values in Eq. (4) up to a small error simultaneously, with
Here we present a novel perspective to the problem of each expressed as the expectation value of a low-rank ob-
eigenenergy estimation, pushing the limits of convergence servable. Before elaborating on the favorable exponential
and robustness via a quantum signal subspace approach. cost reduction, we first introduce our basic measurement-
Our quantum signal space is composed of time correlation driven framework.
functions of the form,

⟨ϕ0 |Oe−iHt |ϕ0 ⟩ , O = O† ∈ CN ×N , (3) A. Basic construction

which captures the system quantum dynamics. Since the 1. Primer: ODMD
expectation value oscillates over time, it can be uniquely
expressed in the natural basis e−iEn t . Representing real- In recent work [9] we examined the observable dynamic
time data in this eigenfrequency basis establishes a clear mode decomposition (ODMD) as a powerful extension of
notion of linear independence in the space of signals. By the classical DMD formalism. ODMD exploits quantum
evaluating the expectation values in Eq. (3) for multiple resources to efficiently measure the expectations of time
independent operators {O1 , O2 , . . . , OI }, we thus manage evolution operators, s(t) = ℜ ⟨ϕ0 |e−iHt |ϕ0 ⟩, rather than
to construct a signal subspace from which we can extract directly tracking the dynamics of the full quantum state
spectral information significantly better. We will refer to e−iHt |ϕ0 ⟩ ∈ H. The extremal (both maximum and mini-
the vector of expectation values, mum) eigenphases e−iEn t and, thereby, the eigenenergies
  En can be inferred via the least-squares (LS) solution to
⟨ϕ0 |O1 e−iHt |ϕ0 ⟩ the following system of linear homogeneous equations,
⟨ϕ0 |O2 e−iHt |ϕ0 ⟩
     
⃗s(t) = 
 .. ,
 (4) s1 s2 · · · sK+1 s0 s1 · · · sK
 .  s2 s3 · · · sK+2  LS  s1 s2 · · · sK+1 
   
⟨ϕ0 |OI e−iHt |ϕ0 ⟩ .
. .. .. ..  =A  ..
  .. . . .. ,
. .

. . .   . . . 
as a multi-observable signal associated with {Oi }Ii=1 . The sd sd+1 · · · sK+d sd−1 sd · · · sK+d−1
central ingredient of our real-time framework is hence the
| {z } | {z }
X′ ∈Rd×(K+1) X∈Rd×(K+1)
efficient collection of a multi-observable signal ⃗s(t), whose (5)
dimensionality I reflects independence or richness of the
underlying spectral information. In particular, the state where sk = s(k∆t) is sampled at a regular time step ∆t,
overlap ⟨ϕ0 |e−iHt |ϕ0 ⟩ in Eq. (1) can be viewed as a simple and (X, X′ ) are time-shifted data matrices containing the
‘one-dimensional’ signal which oscillates over time. overlap evaluations. We note that the Hankel structure
To efficiently measure the multi-observable signal ⃗s(t), of X and X′ immediately implies a compact companion
we leverage classical shadow tomography, specifically the structure of the system matrix, A ∈ Rd×d . The extremal-
shallow shadows recently demonstrated on hardware [35]. phase eigenvalues λ̃ ∈ C of the system matrix A converge
 5

rapidly to the extremal eigenphases λn = e−iEn ∆t of the which establishes a fundamental connection between the
target Hamiltonian as we increase the dimensions d and expectation in Eq. (1) and the trace in Eq. (2). Notably,
K of the data matrices. Γ = Γ† is Hermitian, and its classical simulability entirely
The ODMD algorithm excels in extracting the ground depends on that of |ϕ0 ⟩ and |ϕ⊥ ⟩. Following Eq. (8), it is
state information. As we move towards the interior of the straightforward to show that the density operator ρ(t) in
spectrum, its convergence often slows down progressively fact encodes the time correlation function ⟨ϕ0 |O|ϕ0 (t)⟩ of
or even stagnates. To achieve compact and accurate ex- any Hermitian operator O = O† ∈ L(H), since
cited state estimation – a task more intricate than ground
state estimation – we take advantage of classical shadows N
X −1
cn e−iEn t ,
 
to generate extensive real-time data with minimal quan- ℜ ⟨ϕ0 |O|ϕ0 (t)⟩ = Tr ρ(t)ΓO = ℜ (9)
tum resource overhead. n=0

PN −1 ∗
where cn = m=0 αm αn ⟨m|O|n⟩ and
2. MODMD with classical shadows
ΓO = (Ida ⊗ O) |1, ϕ0 ⟩ ⟨0, ϕ⊥ | + h.c., (10)
Now we generalize the ODMD approach from the time with the identity Ida acting on the ancilla Hilbert space.
evolution operator e−iHt to an operator pool of arbitrary The classical simulability of ΓO now depends on that of
size. For the signal collection process, we invoke the idea |ϕ0 ⟩, |ϕ⊥ ⟩, and O. In particular, ΓO has an efficient clas-
of classical shadow tomography originally introduced [20] sical representation if the two states and the operator are
to extract arbitrary many-body properties from random sparse in the computational and Pauli basis, respectively.
projective measurements. Similar to the Hadamard test, For instance, if |ϕ0 ⟩ and |ϕ⊥ ⟩ are simple computational
we introduce a single ancilla qubit to control elementary basis states, as in the context of quantum chemistry, ΓO
operations on the system registers, demonstrating for the can be classically represented for any Pauli string O.
first time that this setup can in fact efficiently predict the If one measures each term from {ℜ⟨ϕ0 |Oi |ϕ0 (t)⟩}Ii=1
multi-observable signal ⃗s(t) in Eq. (4). individually with a Hadamard-test circuit, it will take a
With the ancilla initialized in |0a ⟩ ∈ Ha , we first create total of O(Iϵ−2 1 ) measurements to achieve a uniform error
a superposition state |Φ(t)⟩ ∈ Ha ⊗ H given by, of ϵ1 . With classical shadow tomography, we significantly
reduce this measurement sample complexity. As shown
|0, ϕ⊥ ⟩ + |1, ϕ0 (t)⟩
|Φ(t)⟩ = √ , (6) in Eq. (9), we can rewrite each expectation value as the
2 trace of a Hermitian, low-rank observable ΓO over a state
ρ(t) ∈ L(Ha ⊗H). Moreover, we emphasize that the state
where |1, ϕ0 (t)⟩ = |1a ⟩ ⊗ |ϕ0 (t)⟩ denotes the product of
preparation circuit of ρ(t) only involves a single ancillary
the ancillary state and the time-evolved reference state
qubit, matching the overhead of a Hadamard-test circuit.
|ϕ0 (t)⟩ = e−iHt |ϕ0 ⟩, and |ϕ⊥ ⟩ ∈ H is any residual state.
As shown in Fig. 1, after preparing ρ(t) with the quantum
In many-body systems, we can conveniently initialize |ϕ0 ⟩
time evolution circuit, we apply a global Clifford shadow
and |ϕ⊥ ⟩ according to particular symmetry sectors of the
protocol on the quantum state and collect shadows {σ̂q =
Hamiltonian H. For the case of molecular Hamiltonians,
it suffices to prepare |ϕ0 ⟩ and |ϕ⊥ ⟩ with distinct fermionic Uq† |bq ⟩⟨bq |Uq }Q
q=1 . On a classical computer, we employ
occupations which, after second quantization, correspond the stored shadows to construct an empirical estimator
to two computational basis states with distinct Hamming for the trace,
weights. In particular, by preparing |ϕ⊥ ⟩ in the vacuum Q
state with zero particle occupation, one can implement 1 X
⟨ΓO ⟩est = Tr[M−1 (σ̂q )ΓO ], (11)
time evolutions e−iHt directly on the system registries in Q q=1
actual experimental settings. For general Hamiltonians,
we note that the requirement for controlled evolution can where M ∈ L(Ha ⊗ H) is linear. For classical shadows
be formally relaxed by jointly time-evolving the system obtained from global random circuits or shallow random
and ancilla qubits, circuits, M can be calculated efficiently [20, 23, 34, 36]. It
can be proven [20] that when Q = O(log(I) ϵ−2 1 ), we can
|Φ(t)⟩ = e−iH̃t |Φ(0)⟩ , (7) estimate all the I expectation values to ϵ1 error with high
probability using the global Clifford classical shadow, i.e.,
under the total Hamiltonian H̃ = 12 (Ida −Za )⊗H, where |⟨ΓOi ⟩est − Tr(ρ(t)ΓOi )| < ϵ1 for all ΓOi . In addition, one
Ida and Za are the identity and Pauli Z operators acting can also read out the imaginary part, ℑ⟨ϕ0 |O|ϕ0 (t)⟩, by
on the single ancilla, respectively. setting ΓO = i (Ida ⊗ O) |1, ϕ0 ⟩ ⟨0, ϕ⊥ | + h.c..
For the composite state ρ(t) := |Φ(t)⟩ ⟨Φ(t)| ∈ L(Ha ⊗ To realize global Clifford random unitaries, one needs
H) and the operator Γ = |1, ϕ0 ⟩ ⟨0, ϕ⊥ | + |0, ϕ⊥ ⟩ ⟨1, ϕ0 | ∈ linear-depth quantum circuits. This could pose a serious
L(Ha ⊗ H), we recognize the pivotal relation that challenge on the near-term quantum platforms due to the
  severe two-qubit gate errors. Fortunately, recent finding
ℜ ⟨ϕ0 |ϕ0 (t)⟩ = Tr ρ(t)Γ , (8) shows that shallow quantum circuits with depth O(log L)
 6

can form global random unitaries on L qubits [37]. More- Algorithm 1: MODMD eigenenergy estimation
over, experiments have demonstrated that classical shad-
ows using log-depth quantum circuits can be made robust Input: Reference state |ϕ0 ⟩, operator pool {Oi }Ii=1 ,
against various quantum errors through new theoretical time step ∆t, noise threshold δ̃.
advancements [35]. These emerging developments enable Output: Estimated energies Ẽn and eigenstates |ψ̃n ⟩.
one to fully leverage the robust shallow shadow technique 1 k←0
in experiments to achieve a low measurement overhead while not converged do
for MODMD. Without loss of generality, we will focus on 2 ⃗sk ← e−iHk∆t |ϕ0 ⟩ ; /* classical shadows */
the global Clifford classical shadow tomography for our 3 X, X′ ← Hankel (⃗s0 , ⃗s1 , . . . , ⃗sk ) ; /* Eq. (13) */
our analysis. 4 Xδ̃ ← X ; /* least-squares regularization */
5 A ← X′ X+δ̃
; /* update system matrix */
6 Ẽn ← −ℑ log(λ̃n )/∆t ; /* energies */
B. Main algorithm 7 |ψ̃n ⟩ ← (a,i) z(a,i),n Oi e−iHa∆t |ϕ0 ⟩ ; /* states
P
*/
8 k ←k+1
With our efficient shadow implementation, we can esti- end
mate the density operator ρ(t) and the associated expec-
tation values for an arbitrary pool of operators {Oi }Ii=1 .
By doing so, we facilitate an exponential expansion of the
where λ̃ℓ give the DMD eigenvalues while (ΨR,ℓ , ΨL,ℓ )
signal subspace relative to the measurement cost for the
denote the corresponding right and left eigenvectors. We
shadow reconstruction, as estimating I observables only
note that λ̃ℓ ∈ C since A, being block companion, is not
requires Q = O(log(I)) samples {σ̂q }Q q=1 . The collection Hermitian.
of expectation values takes the form,
We can read off our eigenpair approximations from the
ordering of the phases arg(λ̃ℓ ). Without loss of generality,
 
⟨ϕ0 |O1 |ϕ0 (t)⟩
−1 we assume that the phases are arranged in a descending
⟨ϕ0 |O2 |ϕ0 (t)⟩ NX
 
⃗s(t) =  = ⃗cn e−iEn t , (12) order, arg(λ̃0 ) ≥ arg(λ̃1 ) ≥ . . ., such that the eigenvalue
.. 
with the maximal phase, λ̃0 := |λ̃0 |e−iẼ0 ∆t ≈ e−iE0 ∆t ,
 
 .  n=0
⟨ϕ0 |OI |ϕ0 (t)⟩ encodes the DMD approximation Ẽ0 to the exact ground
PN −1 ∗ state energy E0 . Likewise, the eigenvalue λ̃1 provides an
for ⃗cn ∈ CI with cn,i = m=0 αm αn ⟨ψm |Oi |ψn ⟩. approximation to the first excited state energy E1 . The
Given our access to real-time expectations ⃗sk = ⃗s(k∆t) eigenstates, on the other hand, can be accessed from the
sampled at time step ∆t, we formulate a LS problem as DMD eigenvectors ΨL,n = [z0,n , z1,n , . . . , zdI−1,n ]. The
the multi-dimensional variant of Eq. (5), left eigenvectors satisfy the eigenvalue equation,
   
⃗s1 ⃗s2 · · · ⃗sK+1 ⃗s0 ⃗s1 · · · ⃗sK ΨL,n X′ = λ̃n ΨL,n X, (15)
⃗s2 ⃗s3 · · · ⃗sK+2  LS  ⃗s1 ⃗s2 · · · ⃗sK+1 
   
. .. ..  =A , where Eq. (15) can be viewed as an equality relating the
..  .. .. . . ..
 
.
. . matrix elements of X and X′ . Such an equality restricted

. . .   . . . 
⃗sd ⃗sd+1 · · · ⃗sK+d ⃗sd−1 ⃗sd · · · ⃗sK+d−1 to the first columns of the data matrices, for example,
| {z } | {z } implies
X′ ∈RdI×(K+1) X∈RdI×(K+1)    
(13) ⃗s1 ⃗s0
⃗s2   ⃗s1 
   
where (X, X′ ) are time-shifted data matrices containing ΨL,n  = λ̃
 ..  Ψ
n L,n  .  , (16)
  
the observable evaluations through shadows. The system .  .. 
matrix A ∈ CdI×dI is now block companion with dI 2 free
⃗sd ⃗sd−1
parameters, i.e., its last I rows. In the special case where
the observable pool contains a single operator O1 , namely which can be expressed in terms of the eigenvector coor-
the identity Id on the system Hilbert space, our approach dinates and real-time observables,
reduces to the original ODMD setting.  
The system matrix A captures the evolution of multi- d−1 X
X I
dimensional expectations propagated by unitary dynam- ⟨ϕ0 |e−iH∆t  z(a,i),n Oi e−iHa∆t |ϕ0 ⟩
ics e−iH∆t . Hence, the eigenenergies and corresponding a=0 i=1
eigenstates of the Hamiltonian H, as the generator of the   (17)
d−1 X
I
dynamics, can be simply estimated via the eigenvalue de- X
composition of A, = λ̃n ⟨ϕ0 |  z(a,i),n Oi e−iHa∆t |ϕ0 ⟩ ,
a=0 i=1
dI−1
X
A= λ̃ℓ ΨR,ℓ ΨL,ℓ , (14) where z(a,i),n := zaI+i−1,n are the vectorized coefficients
ℓ=0 of ΨL,n . Since λ̃n ≈ e−iEn ∆t , the dynamic mode above
 7

closely follows the eigenstate oscillation driven at desired Additionally, ∆t must satisfy a further compatibility con-
frequency, ⟨ϕ0 |e−iH∆t |ψn ⟩ = e−iEn ∆t ⟨ϕ0 |ψn ⟩. Thus, we dition [9],
can approximate |ψn ⟩ by Eq. (17),
2π
∆t ≲ , (20)
d−1 X
I ∆0,N −1 + max (∆n,n+1 − ∆0,n )
X 0≤n<Neig
|ψ̃n ⟩ = z(a,i),n Oi e−iHa∆t |ϕ0 ⟩ , (18)
a=0 i=1 where Neig counts the eigenenergies {E0 , . . . , ENeig −1 } of
interest. For unambiguous and appropriately ordered (so
where z(a,i),n are now scaled to give a normalized state. that Ẽn approximates En ) estimation, we suggest bound-
Any eigenstate properties can in turn be derived in terms ing the spectral range of the Hamiltonian and then lin-
of the pool of operators {Oi }Ii=1 and time-evolved states early shifting the range to be in [−Cπ, Cπ] for some pos-
{e−iHa∆t |ϕ0 ⟩}a . itive constant C < 1. In this case, the time step can be
The formal solution to Eq. (13) entails computing the set to 1, uniquely restricting eigenangles En ∆t in the 2π-
Moore–Penrose pseudo-inverse X+ = (X† X)−1 X† of the window (−π, π). Note that Eq. (20) holds for the relevant
data matrix X. To ensure stability and filter out pertur- energy lower and upper bounds so the exact eigenenergies
bative noise, we employ the following truncated singular do not need to be known in advance.
value decomposition (SVD) of the data matrix, The choice of observables {Oi }Ii=1 also determines the
convergence. Drawing upon the signal subspace intuition
σℓ uℓ v †ℓ ,
X
Xδ̃ = (19) from Eq. (12), the observables should be ‘independent’ in
ℓ:σℓ >δ̃σmax
the sense that the matrix of oscillation amplitudes,
 
c0,1 c0,2 . . . c0,I
where σℓ and (uℓ , v ℓ ) are the singular values and vectors
 c1,1 c1,2 . . . c1,I 
 
respectively. Here δ̃ > 0 is a truncation threshold defined c=  .. .. .. .. , (21)
relative to the largest singular value σmax = maxℓ σℓ of  . . . . 
X. This thresholding procedure, which removes smaller cN −1,1 cN −1,2 . . . cN −1,I
singular values associated with noise, serves to regularize
the LS problem of Eq. (13). maintains a full column rank of I. Otherwise, the multi-
In summary, our shadow-based algorithm requires as dimensional signals contain redundant information. As a
input the selected observables {Oi }Ii=1 , time step ∆t, and convention, we always fix O1 = Id corresponding to the
singular value threshold δ̃. The algorithm is described in ODMD algorithm.
Algorithm 1, which we call the multi-observable dynamic It is worth noting that the shadow reconstruction in-
mode decomposition (MODMD). volves strictly classical computation, thus requiring each
observable Oi to have a sparse representation in the Pauli
basis,
C. Selection of hyperparameters Mi
X
Oi = κi,ν Pi,ν , (22)
ν=1
The performance of our algorithm clearly relies on the
choice of the input parameters. We first remark that the where {Pi,ν }M ν=1 is a set of Mi = O(poly(L)) distinct
i

convergence of MODMD, just as any subspace method, is L-qubit Pauli strings with associated weights {κi,ν }M ν=1 .
influenced by the choice of reference state |ϕ0 ⟩, where the Although it is rather convenient to select the observables
estimation error scales inversely with the squared overlap Oi randomly from the 4L Pauli strings, the resulting real-
with the eigenstates of interest. Although having a larger time signals may suffer from diminished utility because of
overlap is ideal, it suffices to prepare |ϕ0 ⟩ using simplified probable suppression of the target oscillation amplitudes
single-particle calculations, which generally translates to |cn,i | = |⟨ψn |ϕ0 ⟩||⟨ϕ0 |Oi |ψn ⟩|. For ground state estimate,
a sparse sum of product states. Next, the SVD threshold this occurs when |⟨ϕ0 |Oi |ψ0 ⟩| ≈ 0, which deteriorates the
δ̃ is entirely subject to the noise level, which can be con- quality of the signals. As an example, a 1-local Pauli X
trolled as we collect real-time data via low-rank shadow operator changes the Hamming weight of reference state,
techniques. For practical purposes, we set δ̃ ≈ 10εnoise to and can hence lead to zero amplitude if the Hamiltonian
be roughly an order of magnitude above the uncertainty preserves the total Z-spin.
εnoise due to statistical/shot noise. Thus, our focus shifts Alternatively, we propose the systematic generation of
to optimizing the time step and choice of operator pool. observable pools starting from the problem Hamiltonian,
The time step ∆t impacts the algorithmic convergence
M
as it sets the separation of eigenphases e−iEn ∆t over the X
H= κν Pν , (23)
unit circle. A larger time step is advantageous for better
ν=1
distinguishing the eigenphases, until an ambiguity arises
when ∆0,N −1 ∆t ≥ 2π, where the energy gaps are defined where the terms are ordered by the magnitude of their
as ∆m,n = En − Em (so ∆0,N −1 is the spectral range). coefficients, |κ1 | ≥ ... ≥ |κM |. Such sorting induces a
 8
P
family of partial sums, ν∈V κν Pν , with V ⊆ {1, . . . , M } dependent expectations ⟨ϕ0 (t)|Oi |ϕ0 (t)⟩ (note that they
labeling a subset of Pauli strings. Our observables can be differ from ⟨ϕ0 |Oi |ϕ0 (t)⟩). Specifically, we recognize that
selected from these partial sums based on importance of
the M Pauli weights {κν }M ν=1 . P
That is, we may consider Q
2 X
M −1
Tr[M−1 (σ̂q )(|1a ⟩ ⟨1a | ⊗ Oi )],
PM
O1 = ν=1 κν Pν = H, O2 = ν=1 κν Pν , etc. Let us ⟨Oi (t)⟩est = (25)
Q q=1
assume that the target Hamiltonian is linearly shifted,
as discussed for selection of time step, such that the low-
lying energies are large in magnitude. In contrast to ran- gives an unbiased estimate of ⟨Oi (t)⟩ = ⟨ϕ0 (t)|Oi |ϕ0 (t)⟩.
domly selecting a Pauli string, the low energy amplitudes These additional data can be taken as an input to a sep-
of interest, for example |c0,1 | = |⟨ψn |ϕ0 ⟩|2 |E0 |, are effec- arate set of MODMD calculations.
tively ‘magnified’ relative to amplitudes associated with Augmenting the capability of real-time subspace meth-
energies interior in the spectrum. ods to compactly represent the eigenstates, the gen-
We remark that integer powers of the partial sums and eral framework of dynamic mode decomposition (DMD)
their linear combinations can also be desired additions to moreover enables the prediction of system dynamics over
the observable pool for generating high-quality real-time longer timescales [38–40]. A multi-observable signal ⃗s(t),
signals. This imposes no computational bottleneck since composed of time correlation functions, contains essential
the observable predictions can be performed in a parallel dynamical fingerprint that characterizes, for instance,
and distributed manner on classical computers, and the how a many-body system reacts to an external pertur-
variance only depends on the 1-norm, max1≤i≤I ∥⃗κi ∥1 , of bation in the linear-response regime [41, 42]. Here, the
the Pauli weights (see Section IV C). response represents a dynamical property distinct from
the stationary properties governed by a single eigenstate.
To further study these dynamical properties, we analyze
D. Hamiltonian properties beyond energies how time correlation functions are predicted within the
MODMD framework.
The MODMD framework extends beyond estimating First, the one-point correlators ⟨ϕ0 |Oi e−iHt |ϕ0 ⟩ can be
the eigenenergies, providing access to useful Hamiltonian propagated forward in time in increments of ∆t: Eq. (13)
properties including eigenstate properties and dynamical suggests that integer powers of the system matrix A can
responses. To illustrate, we first recall from Section III B be used to (approximately) fast-forward the observables
that the MODMD eigenpairs (λ̃ℓ , ΨL,ℓ ) of system matrix beyond the measurement window. That is, for any inte-
A can also be leveraged to construct compact approx- ger k > K + d, MODMD predicts the dynamics at a later
imations |ψ̃n ⟩ to the low-lying eigenstates |ψn ⟩. Such time k∆t via
eigenstate information can be explicitly translated and
implemented as linear combination of time evolutions on MODMD
⟨ϕ0 |Oi e−iHk∆t |ϕ0 ⟩ ≈ e†i Ak−K−d+1 xK , (26)
quantum hardware, and is unavailable from typical signal
processing methods. As a consequence, arbitrary eigen-
state properties can be predicted as where xK ∈ RdI is the last column of the data matrix X,
and ei = [0, . . . , 0, 1, 0, . . . , 0]⊤ ∈ RdI labels the standard
f (|ψn ⟩) ≈ f (|ψ̃n ⟩), (24) basis vector with the [(d−1)I +i]th entry equal to 1. The
system matrix A functions as the forward ∆t-propagator
for any scalar-valued function f : H → C. This predictive for evolving observables in time, establishing A−1 as the
capability straightforwardly applies to any state property respective backward propagator. Such predictive abilities
within the low-lying energy subspace. of MODMD are illustrated in Appendix F.
In addition, the state shadows {σ̂q }Q
q=1 stored on the Next, the two-point correlators can be approximated,
classical computer can be utilized to calculate the time- with MODMD eigenenergy and eigenstate estimates, as

dI−1
MODMD X
⟨ϕ0 |Oi (k∆t)Oj (l∆t)|ϕ0 ⟩ ≈ e−iẼn (k−l)∆t ⟨ϕ0 |eiHk∆t Oi |ψ̃n ⟩ ⟨ψ̃n |Oj e−iHl∆t |ϕ0 ⟩ , (27)
n=0

where O(t) = eiHt Oe−iHt labels a time-evolved operator the approximate eigenstates |ψ̃n ⟩ alongside |ϕ0 ⟩. Thus by
in the Heisenberg picture. The behavior of the one-point Eq. (18), measuring these snapshots in general incurs an
correlators ⟨ψ̃n |Oi e−iHk∆t |ϕ0 ⟩ over longer times can then additional cost of O(I 2 (K+d)2 ), which makes predictions
be predicted via short-time snapshots ⟨ψ̃n |Oi e−iHa∆t |ϕ0 ⟩ of two-point correlators more demanding. However in the
for 0 ≤ a ≤ K +d. Notice that the inner products involve case where |ϕ0 ⟩ = |ψ̃m ⟩, we can fully leverage information
 9

in the signal subspace with no extra measurements as, for for λn = e−iEn ∆t . While for the complementary case of
instance, e−iHa∆t |ψ̃0 ⟩ ≈ e−iẼ0 a∆t |ψ̃0 ⟩. This reduces the (d, I) = (1, N ), CA (z) = det(z − A) where A ∈ CN ×N
cost back to O(log(I)(K +d)). Importantly, the favorable satisfies the linear homogeneous equation,
log(I) scaling is absent in conventional subspace methods    
where matrix elements of the forms ⟨ϕ0 |Oi e−iHk∆t Oj |ϕ0 ⟩ ⟨ϕ0 |O1 |ϕ0 (t + ∆t)⟩ ⟨ϕ0 |O1 |ϕ0 (t)⟩
or ⟨ϕ0 |eiHk∆t Oi e−iH∆t Oj e−iHl∆t |ϕ0 ⟩ can be measured at  ⟨ϕ0 |O2 |ϕ0 (t + ∆t)⟩ 
 
 ⟨ϕ0 |O2 |ϕ0 (t)⟩ 
 
the costs of at least O(I log(I)(K +d)) or O(I 2 (K +d)2 )), 
 .. 
 = A 
 ..  , (30)

respectively. While MODMD offers versatile capabilities  .   . 
discussed in this section, our work focuses on estimation ⟨ϕ0 |ON |ϕ0 (t + ∆t)⟩ ⟨ϕ0 |ON |ϕ0 (t)⟩
of eigenenergies, leaving the specific explorations of other
applications for future studies. for any time t ∈ R. Such a matrix A indeed exists because
the real-time expectations ⟨ϕ0 |Oi |ϕ0 (t)⟩ reside within the
N -dimensional space spanned by single-frequency signals
IV. THEORETICAL GUARANTEES e−iEn t driven at the individual eigenfrequencies En . To
investigate the general case 1 ≪ dI ≪ N , we start with a
We establish in this section fundamental connections simpler version of Eq. (13), assuming the d matrix blocks
between MODMD and modern spectral approaches, fur- to be diagonal:
nishing a theoretical framework that guarantees its con-  
vergence. We first explicitly show a speedup of MODMD Aℓ,11
over ODMD by having an expansive pool of observables. 
 Aℓ,22


Next, we exploit the linearity of quantum dynamical evo- Aℓ = Diag  ..
, (31)
.
 
lution and consider MODMD as a multi-reference scheme  
within a suitably defined function space through Koop- Aℓ,II
man operator analysis [43–48]. These two analytical
viewpoints reinforce each other, underpinning the reliable where Aℓ,ij ≡ 0 for i ̸= j. In this case, the individual
performance of our algorithm for excited state problems. LS residuals associated with the observables Oi are inde-
pendent, and we can show that the resulting MODMD
estimates are bounded by the I single-observable ODMD
A. Multi-observable dynamic mode decomposition estimates, e.g.,

For the least-squares (LS) problem of Eq. (13) with min |Ẽi,0 − E0 | ≤ |Ẽ0 − E0 | ≤ max |Ẽi,0 − E0 |, (32)
1≤i≤I 1≤i≤I
observables {Oi }Ii=1 , the multivariate solution is,
where Ẽ0 and Ẽi,0 indicate, respectively, the ground state
 
0
energy estimate using the entire observable pool (the full
 ..
 
 . Idd−1 ⊗ IdI 
system matrix A) and one single observable (only ith row
A= , (28)

 0  Aℓ,ii of the submatrices Aℓ ).
More interesting convergence arises when the I single-
 
−A0 −A1 · · · −Ad−1
observable residuals are coupled to one another, allowing
reductions in the individual residuals due to the flex-
where each Aℓ ∈ RI×I for ℓ = 0, 1, · · · , d − 1 represents ibility of off-diagonal elements in the submatrices Aℓ .
a submatrix, and Idd−1 and IdI are identity operators of Specific scenarios in which MODMD substantially im-
respective dimensions. The block companion structure of proves upon the ODMD residuals are considered within
the system matrix A ∈ CdI×dI immediately follows from Appendix C. Intuitively, we expect a reduced total resid-
the block Hankel structure of the matrices X and X′ . ual to, in turn, improve the eigenenergy estimates, where
The multi-observable system matrix has a characteristic Eq. (32) holds accordingly with tighter lower and upper
polynomial, bound – ideally, both approaching zero. Moreover, the

d
 reduction in total LS residual from MODMD indicates a
X more expressive system matrix as a proxy for the under-
CA (z) = det(z − A) = det  z ℓ Aℓ  , (29)
lying dynamics, which is crucial for accurately predicting
ℓ=0
a multi-observable signal over longer times as discussed
where Ad ≡ IdI . Hence, the roots of CA (z) correspond in Section III D.
precisely to the dI eigenvalues of the system matrix A.
To understand the convergence of the A-spectrum to the
actual eigenphases, we first consider two limiting cases. B. Koopman operator analysis
For the case (d, I) = (N, 1), we essentially recover the
celebrated Prony’s method [49] as the single-observable We now establish convergence properties of MODMD
QN −1
ODMD with Aℓ ∈ C, such that CA (z) = n=0 (z − λn ) from a functional-theoretic perspective. Our analysis will
 10

revolve around the study of the Koopman operator, a piv- spanned by these eigenfunctions. By choosing |ψ⟩ and O,
otal mathematical object in understanding the complexi- we recast the task of identifying Hamiltonian eigenmodes
ties of a dynamical system. The Koopman operator K∆t as an equivalent, finite-dimensional Koopman eigenvalue
probes the underlying dynamics of a system by acting on problem. We use ⃗g = [g1 , g2 , . . . , gD ]⊤ to denote a vector
scalar-valued functions. For any function f : H → C, its of D distinct scalar functions gi in the invariant function
action is subspace, all taking the form of fψ,O .
We first examine the case where gi = (K∆t )i−1 [fψ,Id ],
K∆t [f ](|ϕ⟩) = f (e−iH∆t |ϕ⟩), (33) which gives rise to the ODMD approach when D = d. We
seek to determine the closest approximation K∆t ∈ Cd×d
of the Koopman operator when restricted to the invariant
which gives a push-forward of the dynamics via the time subspace. The closest-fitting problem has a least-squares
evolution operator e−iH∆t . For a quantum-dynamical formulation,
system, we first observe that fn (|ϕ⟩) = ⟨ψn |ϕ⟩ constitutes
an eigenfunction of K∆t since K
X 2
K∆t = arg min K∆t [⃗g ](|ϕk ⟩) − K⃗g (|ϕk ⟩) 2
, (35)
−iEn ∆t K∈Cd×d
K∆t [fn ] = e fn , (34) k=0

where K∆t acts component-wise on ⃗g , and the ℓ2 -residual

−iEn ∆t
where e is the corresponding eigenvalue. There- is being minimized with respect to some states {|ϕk ⟩}Kk=0
fore scalar functions of the form fψ,O (|ϕ⟩) := ⟨ψ|O|ϕ⟩ = sampled from the Hilbert space. Formally, we can express
PN −1
n=0 ⟨ψ|O|ψn ⟩ ⟨ψn |ϕ⟩ lie in a K∆t -invariant subspace Eq. (35) as a matrix equation

   
   
⃗g (e−iH∆t |ϕ0 ⟩) · · · ⃗g (e−iH∆t |ϕK ⟩) = K∆t ⃗g (|ϕ0 ⟩) · · · ⃗g (|ϕK ⟩), (36)
   

| {z } | {z }
G′ ∈Cd×(K+1) G∈Cd×(K+1)

whose solution K∆t = G′ (G)+ involves computing the eigenvalue equation,

pseudo-inverse G+ . Utilizing G+ = (G† G)−1 G† , we can
directly show that K∆t yields an equivalent factorization W⃗vK = λK Z⃗vK , (41)
K∆t = W(Z)+ , where Z, W ∈ Cd×d are constructed by
where Z and W can now be reinterpreted, respectively,
K
1 X as the matrix representation of the functional overlap and
Z= ⃗g (|ϕk ⟩)⃗g † (|ϕk ⟩), (37) d−1
Koopman operator in the finite basis {g1 , · · · , K∆t [g1 ]}.
K +1
k=0 This special nonorthogonal basis is known as the order-d
K Krylov basis, and the associated subspace is the Krylov
1 X
W= ⃗g (e−iH∆t |ϕk ⟩)⃗g † (|ϕk ⟩). (38) subspace Kd (g1 ) ⊂ span{fn }n .
K +1
k=0 Projecting the full Koopman eigenvalue problem onto
this subspace of size d ≪ N allows an efficient retrieval of
This alternative factorization is intimately related to the
spectral information, including the extremal eigenvalues.
powerful Krylov approaches for operator diagonalization
In the ODMD algorithm, we choose |ϕk ⟩ = e−iHk∆t |ϕ0 ⟩
[50–52]. Suppose we sample the states {|ϕk ⟩}k according
for a fixed state |ϕ0 ⟩ and take the sampling measure µ to
to a probability measure µ defined over the Hilbert space.
be the empirical measure (a sum of Dirac-delta meaures),
With sufficiently many samples, we have
Z K
1 X
Zij → dµgi gj∗ = ⟨gj , gi ⟩µ , (39) µ= δ|ϕk ⟩ , (42)
|ϕ⟩∈H K +1
Z k=0

Wij → dµK∆t [gi ]gj∗ = ⟨gj , K∆t [gi ]⟩µ , (40)

|ϕ⟩∈H along the orbit (|ϕ0 ⟩ , e−iH∆t |ϕ0 ⟩ , · · · , e−iHK∆t |ϕ0 ⟩). In
this case, limK→∞ ∥K∆t [f ]∥L2 (µ) = limK→∞ ∥f ∥L2 (µ) for
where the L2 (µ) inner product gives the continuum limit a continuous function f (we expect the limit to exist [53]
of Monte-Carlo averages in Eqs. (37) and (38). Therefore for a generic many-body Hamiltonian H and time step
eigenpairs (λK , ⃗vK ) of K∆t should satisfy the generalized ∆t), since
 11

K+1 2 K 2
1 X 1 X
lim ⟨K∆t [f ], K∆t [f ]⟩µ = lim f (e−iHk∆t |ϕ0 ⟩) = lim f (e−iHk∆t |ϕ0 ⟩) = lim ⟨f, f ⟩µ ,
K→∞ K→∞ K + 1 K→∞ K + 1 K→∞
k=1 k=0
(43)

where the middle equality holds due to a vanishing differ- can now simultaneously compute the Koopman action on
ence, limK→∞ K+1 1
[|f (e−iH(K+1)∆t |ϕ0 ⟩)|2 −|f (|ϕ0 ⟩)|2 ] = a vector ⃗g of scalar functions. Each function corresponds
0. In other words, the Koopman operator K∆t , when re- to a unique choice of observable. This key algorithmic
stricted to the invariant subspace, is isometric and hence improvement enables a block Krylov scheme, significantly
normal. By the spectral theorem, the K∆t -eigenfunctions accelerating the energy convergence. The rate of conver-
−1
{fn }N
n=0 are orthogonal to each other. Under this condi- gence in the MODMD setting is described by the theorem
tion, exponentially rapid convergence of standard Krylov below.
approach has been thoroughly analyzed [54], which aligns
consistently with our ODMD observations. Theorem 2. Let Ẽn (d) be the approximate nth eigenen-
Now we explicitly assert an exponential convergence of ergy extracted in the dI-dimensional function subspace
d−1
the functional Krylov approach when the eigenfunctions span{⃗g , K∆t [⃗g ], · · · , K∆t [⃗g ]}, and δEn (d) = Ẽn (d) − En
are orthogonal, i.e., be the approximation error. Consider the diagonal error
Z matrix
δmn
⟨fm , fn ⟩µ = dµ ⟨ψn |ψ⟩ ⟨ψ|ψm ⟩ = , (44)  
|ψ⟩∈H N δE0 (d)
∆I (d) = 
 .. 
, (47)
where Eq. (44) holds if we randomly sample the states |ψ⟩  . 
i.i.d. via a quantum 1-design. For the single observable δEI−1 (d)
case, we refer to the following theorem for the exponential
decay of the estimation error in the ground state energy. which contains approximations to the lowest I energies.
For d ≥ 1, there exists a time step ∆t so that the spectral
approximation ∆I (d) is bounded by
Theorem 1. Let Ẽ0 (d) be the approximate ground state
energy extracted in the d-dimensional function subspace
d−1 sin[(EN −1 − E0 )∆t]
spanned by {g1 , K∆t [g1 ], · · · , K∆t [g1 ]}. For d ≥ 1, there ∥∆I (d)∥2 ≤ 2(d−1)
∥tan2 Θ∥2 , (48)
exists time step ∆t so that the error δE0 (d) = Ẽ0 (d)−E0 ϵ̃I−17→I ∆t
is bounded by
for the operator norm ∥·∥2 . Here Θ denotes the canonical
sin[(EN −1 − E0 )∆t] angle between the two subspaces span{fn : 0 ≤ n ≤ I−1}
δE0 (d) ≤ 2(d−1)
tan2 Ξ, (45) and span{gi : 1 ≤ i ≤ I}, which generalizes the squared
ϵ̃0→1 ∆t
overlap in Theorem 1 (see Appendices D 2 and D 3). In
where cos2 Ξ = |⟨N f0 , g1 ⟩µ |2 = |⟨ψ0 |ϕ0 ⟩|2 is the squared the denominator, ϵ̃I−17→I ∈ [1, 2] depends on the Ith
overlap between the reference function g1 and the true spectral gap (EI − EI−1 ) of the Hamiltonian H.
ground state K∆t -eigenfunction, while ϵ̃0→1 = 1 + 3(E1 − Proof. The proof is provided in Appendix D 3.
E0 )∆t/2π ∈ [1, 2] characterizes the normalized spectral
gap of the Hamiltonian H. Despite the formal similarity between the bounds from
Proof. The proof is provided in Appendix D 1. Theorem 1 and Theorem 2, we highlight that convergence
in a multi-observable setting can offer distinct advantages
Building on the single-observable Krylov idea (ODMD) for excited state calculations. The ODMD bound for the
above, the MODMD approach can thus be viewed as an Ith lowest eigenenergy, as in standard subspace methods
enriched extension, where we allow for an extra degree of employing the reference state |ϕ0 ⟩, includes an additional
2I−2


freedom in selecting multiple functions multiplicative factor of Ω 2I−1 ϵ̃I−17 →I [4]. Notably, this
  prefactor grows exponentially as we approach the higher
fϕ0 ,Id excited states, counteracting the single-observable error
fϕ0 ,O2 
  decay from Theorem 1 unless d > I. While it is nat-
−iHt ural to extend the standard subspace methods using I
⃗g =  . 

 =⇒ ⃗g (e |ϕ0 ⟩) = ⃗s(t). (46)
 ..  reference states |ϕi ⟩ ∝ Oi |ϕ0 ⟩, the quantum cost of mea-
fϕ0 ,OI suring the relevant expectations ⟨ϕi |e−iHt |ϕj ⟩ is at least
O(I log(I)). This is because each state |ϕi ⟩ must be time-
By leveraging the classical shadows, each quantum circuit evolved, making it exponentially more expensive than
e−iHk∆t originally capable of computing a single overlap MODMD.
 12

C. Error analysis specifically the Trotter scheme for illustration and com-
ment that similar analysis should hold for other schemes.
PM PM
To account for noisy quantum hardware, we present in For H = ν=1 Hν = ν=1 κH,ν Pν , a first-order Trotter
this section a preliminary error analysis. For our basic formula gives,
considerations, we examine error components of two dis- !r !
M 2
tinct kinds, the statistical error arising from the random- −iH∆t
Y
−iHν ∆t/r M 2 ∥⃗κH ∥∞ ∆t2
izing shadow protocol and the deterministic error due to e − e =O ,
ν=1
r
the imperfect compilation of the time evolution. We show 2
that both errors remain independent of the problem size (51)
across practical range of observable selections. Moreover,
the latter grows linearly with the maximal evolution time where the Trotterized Hamiltonian simulation can be
when we implement Trotterized evolution as viable proxy performed efficiently on a quantum computer if M =
to the actual evolution. O(polylog(N )). Specifically for
2
r = O(M 2 ∥⃗κH ∥∞ ∆t2 ϵ−1
1 ) (52)
1. Statistical noise
Trotterized blocks with a time discretization ∆tTrotter =
∆t/r, we make an compilation error of at most ϵ1 in the
First, we note that the prediction error associated with operator norm. We further assume that the time step ∆t
randomized measurements in the classical shadow tech- is suitably chosen such that,
niques admits, in our case, a variance independent of the
system size L. The classical shadows can predict observ-
E0 − E0Trotter ∆t ≤ O(ϵ1 ), (53)
ables, γi (t) = Tr[ρ(t)ΓOi ], of the form Eq. (9), where we
recall that rank(ΓOi ) ≤ 2 even when rank(Oi ) = 2L = N .
The low-rank property implies [20], where E0Trotter represents the approximate ground state
energy associated with the Trotterized time evolution op-
Trotter
Tr[ΓOi ]
2
U =Cl(N ) erator. As Eq. (51) implies Id − ei(H−H )∆t
2
≤ ϵ1
V[γi ] ≤ ΓOi − Id = 3Tr[Γ2Oi ], (49) and thus (H − H Trotter
)∆t 2 ≤ 2 arcsin [ϵ1 /2] = O(ϵ1 ),
N sh
our assumption holds as a straightforward consequence
where ∥·∥sh is the shadow norm conditional on the mea- of Weyl’s theorem; a necessary condition on the time step
surement primitive and the trace Tr[Γ2Oi ] on the RHS of ∆t to avoid any eigenphase ambiguity is ∥H∥2 ∆t < π.
Eq. (49) can be unfolded by the defining relation ΓOi = Defining γi (t) = ℜ ⟨ϕ0 |Oi e−iHt |ϕ0 ⟩, we have for a fixed
Ida ⊗Oi |1, ϕ0 ⟩ ⟨0, ϕ⊥ |+|0, ϕ⊥ ⟩ ⟨1, ϕ0 | Ida ⊗Oi . Hence, the Trotter time interval ∆tTrotter ,
NL
prediction variance for a Pauli string, P = ℓ=1 σℓ with
σℓ ∈ {Idℓ , Xℓ , Yℓ , Zℓ }, can be uniformly bounded through γi (k∆t) − γiTrotter (k∆t) ≤ kϵ1 ∥Oi ∥2 ≤ kϵ1 ∥⃗κi ∥1 , (54)
2
the Cauchy–Schwarz inequality, Tr[Γ2P ] ≤ 4∥P ∥2 ≤ 4.
This immediately implies which corresponds to a systematic deviation of O(ϵ1 ) in
the data matrix element.
Mi
X 2
Oi = κi,ν Pν =⇒ V[γi ] ≤ O(∥⃗κi ∥1 ), (50)
ν=1 V. APPLICATIONS

for a general Hermitian operator Oi = Oi† . Observe that In this section, we detail numerical studies conducted
the variance remains rather insensitive to the operator on representative many-body systems from condensed-
locality or the system size. Accordingly, each data matrix matter physics and quantum chemistry to demonstrate
element in the MODMD setting incurs an additive error the efficacy of the MODMD framework. Our numerical
2
of at most ϵ1 if we take Q = O(log(I) max1≤i≤I ∥⃗κi ∥1 ϵ−2
1 ) experiments precisely follow prescriptions in Section III.
quantum measurements.

A. Condensed-matter physics
2. Trotter error
We examine the convergence of the ground and excited
A second error source pertains to inexact implementa- state eigenenergies of the 1D transverse field Ising model
tion of the unitary evolution e−iHt , which perturbs both (TFIM) with a total of L = 15 spins and open boundary
the eigenfrequencies e−iEn t and the time-evolved states conditions. The system Hamiltonian is given by
ρ(t). For near-term implementation, we assume query ac-
cess to an approximate compilation of the evolution, for L−1
X L
X
example, through Trotter–Suzuki factorization [55, 56] or HTFIM = −J Zi Zi+1 − h Xi , (55)
linear combination of unitaries (LCU) [57]. We consider i=1 i=1
 13

101 | E0| 101

| E1|
100 | E2| 100
| E3|
Absolute Error

Absolute Error
10 1
noise 10 1

10 2 10 2

10 3 10 3

10 4 10 4

0 100 200 300 400 500 0 100 200 300 400 500
K ( Simulation Time) K ( Simulation Time)
(a) MODMD algorithm (b) ODMD algorithm

FIG. 2: Convergence of eigenenergies for the transverse-field Ising model (TFIM). To obtain eigenenergy estimates Ẽn , we fix
the (M)ODMD parameters K d
= 25 and δ̃ = 10−2 for constructing and thresholding the pair of data matrices X, X′ ∈ RdI×(K+1) .
Gaussian N (0, εnoise ) noise with εnoise = 10−3 is added independently to the real or/and imaginary parts of the matrix elements.
2

The absolute errors, |Ẽn −En |, in the first four eigenenergies of the TFIM Hamiltonian are shown with respect to K proportional
to the non-dimensional maximal simulation time. The reference state |ϕ0 ⟩ is an equal superposition of six computational basis
states and we use a time step of ∆t ≈ 0.08. The shading above the solid/dashed lines shows the standard deviation across trials
for each quantity. Left Absolute errors from the multi-observable (MODMD) algorithm with I = 6, where the convergence
results are averaged over 20 trials, each involving a Gaussian noise realization and a selection of I random observables. Right.
Absolute errors from the single-observable (ODMD) algorithm where the convergence results are averaged over 20 trials, each
involving a Gaussian noise realization.

for coupling constant J and external field strength h. 100

We first demonstrate the performance of MODMD for
TFIM parameters fixed at J = h = 1, examining how the
algorithm behaves with varying number of time steps, or 10 1
Absolute Error in E1

equivalently the maximal simulation time. As discussed

in Sections II and IV, increasing the dimensions d and K
of the data matrix X ∈ RdI×(K+1) facilitates convergence
10 2
of the MODMD eigenvalues λ̃n = |λ̃n |e−iẼn ∆t (note that
these eigenvalues are not necessarily confined to the unit
circle in complex plane) to the eigenphases λn = e−iEn ∆t ODMD
10 3
of the Hamiltonian. Fig. 2 illustrates this convergence for MODMD
the first Neig = 4 eigenenergy estimates. Specifically, we
noise
report the absolute error |δEn | = |Ẽn − En | as a function
of the dimension K, with the ratio K 5
d = 2 held constant
10 4 10 3 10 2 10 1
throughout subsequent calculations as this ratio provides E1 E0
near-optimal performance for ODMD [9].
The left panel illustrates convergence in the multi- FIG. 3: Convergence of the first excited state energy of the
TFIM. Absolute error in the first excited state energy is shown
observable setting, where we select I = 6 operators ran-
as a function of the spectral gap between the ground and the
domly from the set of 1-local Pauli gates and initialize a first excited state for fixed K = 500. The vertical dotted line
reference state |ϕ0 ⟩ composed equally of 6 computational marks the noise level εnoise = 10−3 . Convergence results are
basis states. The use of simple random Pauli observables averaged over 20 trials each involving a Gaussian noise real-
proves effective here since (1) the TFIM does not preserve ization and, in the MODMD case, also a selection of I random
the Hamming weight (local bit-flip X does not annihilate observables. The shading above the solid/dashed lines shows
⟨ϕ0 |Oi |ϕ0 (t)⟩) and (2) our reference is a superposition (so the standard deviation across trials for each quantity.
local phase-flip Z does not add a trivial ±1 overall phase).
The reference |ϕ0 ⟩ contains relatively small but sufficient
overlap with the first few Hamiltonian eigenstates, where nal without requiring substantial similarity between the
PNeig −1
the squared overlap sums to n=0 |⟨ψn |ϕ0 ⟩|2 ≈ 10−1 . reference state and target eigenstates. To stimulate the
This allows for the algorithm to generate an adequate sig- shadow-induced errors, we additionally introduce Gaus-
 14

101 | E0| 101

| E1|
100 | E2| 100
10 1 | E3| 10 1
Absolute Error

Absolute Error
noise
10 2 10 2

10 3 10 3

10 4 10 4

10 5 10 5
0 50 100 150 200 250 300 0 50 100 150 200 250 300
K ( Simulation Time) K ( Simulation Time)
(a) MODMD algorithm (b) ODMD algorithm

FIG. 4: Convergence of eigenenergies for the lithium hydride (LiH) molecule. To obtain eigenenergy estimates Ẽn , we fix the
(M)ODMD parameters K d
= 52 and δ̃ = 10−2 for constructing and thresholding the pair of data matrices X, X′ ∈ RdI×(K+1) .
Gaussian N (0, εnoise ) noise with εnoise = 10−3 is added independently to the real or/and imaginary parts of the matrix elements.
2

The reference state |ϕ0 ⟩ is an equal superposition of the lowest four Hartree-Fock states with identical particle occupation number
and we use a time step of ∆t ≈ 0.33. The absolute errors, |Ẽn − En |, in the first four eigenenergies of the LiH Hamiltonian are
shown with respect to K proportional to the non-dimensional maximal simulation time. Convergence results are averaged over
20 trials of Gaussian noise realizations with shading above the solid/dashed lines showing the standard deviation across trials
for each quantity. Left. Absolute errors from the multi-observable (MODMD) algorithm with I = 6. Right. Absolute errors
from the single-observable (ODMD) algorithm.

sian N (0, ε2noise ) noise with εnoise = 10−3 to the multi- modynamic limit L → ∞, the gap closes at Jh = 1 and
observable signal. The absolute energy errors shown are increases monotonically with Jh . To investigate the gap
averaged over a total of 20 realizations of both the Gaus- dependence near a phase transition, we demonstrate the
sian noise and 1-local Pauli observables. For comparison, difference in performance between the single-observable
the right panel displays the single-observable convergence (ODMD) and multi-observable (MODMD) algorithms in
from the standard ODMD algorithm as our benchmark. the presence of near-degenerate target energies. In ??,
We observe considerably faster convergence in the excited we focus on comparing the error |δE1 | in the first excited
state energies in MODMD for cases where ODMD nearly energy against the gap E1 − E0 . The convergence results
stagnates. We highlight that the quantum cost, or total indicate that, in a multi-observable approach, the first
number of shots required, is at most 2(K +d) log(I)ε−2noise , excited state energy E1 can be accurately estimated if
only a factor of 2 log(I) < 4 more than that of ODMD in the noise level is slightly smaller than the spectral gap.
this case. In contrast, ODMD requires a visibly higher gap-to-noise
The convergence of MODMD naturally divides into ratio to achieve a comparable accuracy. This highlights
two regimes, each defined by a distinct error scaling. The a significant improvement of MODMD in distinguishing
noise level εnoise essentially determines the crossover be- near-degenerate eigenstates.
tween the two convergence regimes. When |δEn | > εnoise ,
we observe an exponentially decaying error typical of the
classical subspace methods [50, 52, 54]. Conversely, in the B. Quantum chemistry
regime where |δEn | < εnoise , increasing simulation time
leads to slower, algebraic error decay with precision ulti- Electronic structure calculation is a fundamental prob-
mately limited by Heisenberg scaling. This crossover be- lem in quantum chemistry. Here we evaluate performance
tween the exponential and algebraic error decay is shown of the MODMD algorithm for molecular Hamiltonians,
explicitly within Appendix E 2. In practice, the onset of whose second quantization can be efficiently implemented
the algebraic error behavior can be considerably delayed, on the quantum computer. Unlike our condensed-matter
for example, via increased and noise-mitigated sampling example where the first few Hamiltonian eigenvalues map
of classical shadows. directly to the energy levels of interest, we only consider
Despite its simplicity, the TFIM is an instructive toy eigenvalues corresponding to eigenstates with the correct
model as it undergoes a quantum phase transition, where particle number. This symmetry is preserved under time
the spectral gap between the first two eigenstates can be evolution.
systematically tuned by varying the Jh ratio. In the ther- Moreover, we construct observables based on the Pauli
 15

truncation strategy can be used to protect the actual sig-

10 1 | E0| nal from noise across a wide range of noise magnitudes.
| E1| By comparison, the ODMD algorithm yields larger errors
| E2| for any given εnoise . Furthermore, reducing noise only im-
10 3
| E3| prove the ODMD performance when εnoise reaches small
Absolute Error

values. In other words, a more aggressive truncation δ̃

10 5 must be employed, which may cause a serious loss of the
signal information and thereby slower convergence.
10 7

10 9
VI. CONCLUSION
102 103 104 105 106 107 108
1/ noise
In this work, we developed a hybrid quantum-classical
FIG. 5: Convergence of eigenergies for LiH under varying measurement-driven framework for effectively extracting
noise level. The absolute errors in the first four eigenener- information about the low-energy eigenspaces of quan-
gies of the LiH Hamiltonian are shown as a function of the tum many-body systems. Our novel MODMD approach
noise level εnoise for fixed K = 300 and δ̃ = 10εnoise . Con- leverages real-time evolution on quantum hardware and
vergence results are averaged over 20 trials of Gaussian noise classically unravels multi-dimensional signals, composed
realizations. The solid lines depict the results for the multi- of real-time observables, from a limited number of ran-
observable (MODMD) algorithm with I = 6 and the dashed domized measurements. The simultaneous prediction of
lines depict the results for the single-observable (ODMD) al- many observables leads to accurate estimates of eigenen-
gorithm.
ergies and shallower circuits with shorter evolution time.
We explored the theoretical underpinnings of MODMD,
which exponentially suppresses spectral error in the low-
representation of the Hamiltonian as per Eq. (23), and noise regime. We numerically demonstrated its rapid
deterministically select a subset of I = 6 Pauli operators convergence in the presence of perturbative noise using
Pν with medium-magnitude weights κν to maximize the examples from condensed matter physics and quantum
signal independence among the observables. The refer- chemistry.
ence state |ϕ0 ⟩ is prepared as an equal superposition of
the lowest four Hartree-Fock states with identical particle Compared to state-of-the-art real-time approaches, we
occupation number. These Hartree-Fock states are com- highlight the unique strengths of our method in addition
putational basis states readily derived from a mean-field to its reliable convergence and noise resilience. To our
calculation. This choice ensures both a nonzero overlap best knowledge, MODMD is among the most resource-
with the target eigenstates and viability for experimental efficient for generating real-time signals. This is be-
preparation. cause (1) we evolve a single reference state for a duration
We illustrate MODMD for the lithium hydride (LiH) shorter than required by single-observable approaches,
molecule. Fig. 4 shows the convergence for the first where the reduction in the simulation time becomes more
Neig = 4 eigenenergy estimates, constrained to fixed par- substantial as the number of observables included in the
ticle occupation and bond length of 1.59Å at equilibrium signal subspace increases, and (2) the reference state does
geometry. Similar to the TFIM results in Fig. 2, Fig. 4 not have to possess large overlaps with the low energy
exhibits, roughly, an exponential-to-algebraic crossover eigenstates of interest. Furthermore, our classical post-
in the error decay, with transition occurring around the processing consists of solving a simple least-squares prob-
noise level εnoise (see Appendix E 2). The plateaus for lem followed by a standard eigenvalue problem, which is
the higher energy levels at intermediate K values can be ansatz-free and thus circumvents an exponential growth
attributed to the anomalous convergence to different but in optimization costs, whether quantum or classical, as-
higher-lying eigenenergies, which we verified numerically. sociated with the number of desired eigenenergies.
As we approach the interior of the Hamiltonian spectrum, In fact, our MODMD framework is capable of retriev-
MODMD can resolve the higher excited state energies for ing ground and excited state properties beyond eigenen-
sufficiently large K, whereas ODMD stagnates. ergy levels, demonstrating an extensive and timely appli-
In Fig. 5 we assess the noise robustness of the MODMD cation of the low-rank shadow. Building upon recent the-
algorithm in our LiH example. Adopting the same hy- oretical progress in adaptive time scheduling, we finally
perparameters as used in Fig. 4, we now fix K = 300 and comment that our algorithm may in principle saturate
vary the noise level εnoise . Importantly, we maintain the the optimal Heisenberg-limited scaling for phase estima-
SVD threshold δ̃ = 10εnoise at a consistently larger value. tion (here the multi-eigenvalue estimation) under some
We observe a power law scaling of the absolute error with additional spectral gap assumptions [16]. This promising
respect to the noise level, suggesting that a conservative prospect warrants further analysis in the future work.
 16

DATA AVAILABILITY

Our numerical experiments are available through the

GitHub repository: ‡MODMD.

ACKNOWLEDGEMENTS

This work was funded by the U.S. Department of En-

ergy (DOE) under Contract No. DE-AC02-05CH11231,
through the Office of Science, Office of Advanced Scien-
tific Computing Research (ASCR) Exploratory Research
for Extreme-Scale Science (YS, KK, DC, RVB). This re-
search used resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Depart-
ment of Energy Office of Science User Facility located
at Lawrence Berkeley National Laboratory, operated un-
der Contract No. DE-AC02-05CH11231. HYH would
like to acknowledge the support from the Harvard Quan-
tum Initiative. SFY and HYH would like to acknowledge
funding from the DOE through the QUACQ center (DE-
SC0025572).
 17

APPENDIX

Appendix A: Shadow under locally scrambling dynamics

Here we adopt a dynamical perspective on the efficient generation of the classical shadow, which introduces random
scrambling via real-time evolution under a disordered local Hamiltonian Hs (t). For concreteness, we focus on quantum
spin glass, in particular the disordered transverse field Ising model (TFIM) with Hamiltonian,
L
X L
X
Hs (τ ) = − Jij (τ )Xi Xj − h(τ ) Zi , (A1)
⟨i,j⟩ i=1

where Jij and h set the spin-spin P coupling and external field strength respectively. We assume for simplicity that Jij =
χ χ χ
P
k Jij,k [τk ,τk+1 ) (τ ) and h = k hk [τk ,τk+1 ) (τ ) remain piecewise constant in time, with [τk ,τk+1 ] being an indicator
function of the kth-step interval. To incorporate scrambling disorder, we employ Q0 quenched random interactions, e.g.,
Jij,k ∼ N (0, σJ2 ) and hk ∼ N (0, σh2 ). Thus under the time evolution Us = k=Ks e−iHs (τk )(τk+1 −τk ) , local quantum
information gets scrambled and the resulting entanglement produced by the disorder dynamics facilitates the recovery
of ρΦ (t). We remark that the number of scrambling time steps Ks serves the role analogous to the circuit depth when
considering Clifford gates [58]. For any integer T ≥ 0, limKs →∞ U({τk }0≤k≤Ks ) is a T -design, whereby its first T
statistical moments align with those of the Haar measure. To simulate long-time dynamics, we draw upon our recent
tool of algebraic circuit compression [59, 60], particularly suited to a disordered TFIM Hamiltonian Hs (τ ) in Eq. (A1),
to keep the effective depth of the time evolution Us shallow. Remarkably, the depth post compression should exhibit
no dependence on the maximal runtime τKs , therefore allowing a significantly more efficient exploration of the Haar
limit.

Appendix B: Observable dynamic mode decomposition

The standard dynamic mode decomposition (DMD), originally developed in the field of numerical fluid dynamics, is
a measurement-driven approximation for the temporal progression of a classical dynamical system [61–65]. Specifically,
DMD samples the system snapshots at regular time intervals ∆t and uses them to construct an efficient representation
of the full dynamical trajectory. For simplicity, we consider a system whose N -dimensional state manifold is CN . The
optimal linear approximation for the discretized time step k 7→ k + 1 is expressed as the least-squares (LS) relation,
LS
ϕk+1 = Ak ϕk , (B1)

where ϕk ∈ CN specifies the system state at time k∆t and Ak ∈ CN ×N is the system matrix, i.e., the linear operator
that minimizes the residual ∥ϕk+1 −Ak ϕk ∥2 to yield the LS relation above. Similarly, the optimal linear approximation
for a sequence of successive snapshots k = 0, 1, . . . , K + 1 can be determined by the solution,
   
  LS  
K+1  = A ϕ0 ϕ1 · · · ϕK  , (B2)
ϕ ϕ · · · ϕ   
 1 2

where the system matrix A minimizes the sum of squared residuals over the length-K sequence. The linear flow
t
described by Eq. (B2) naturally generates approximate dynamics ϕDMD (t) = A ∆t ϕ0 governed by eigenmodes of A.
DMD-based approaches can be remarkably effective despite their formal simplicity, since they are rooted in the general
Koopman operator theory developed to describe the behavior of general (non)linear dynamical systems [43–45, 47, 48].
The standard DMD approach described above for classical dynamics cannot be immediately translated to quantum
dynamics. The DMD approximation of the system evolution would require complete knowledge of the system state,
as specified by an N -dimensional complex vector at each time step. However, we do not have direct access to the full
many-body quantum state. Instead, we can only access the state of a quantum system via measurement sampling of
observables [66]. To address this challenge, we employ a technique motivated by Takens’ embedding theorem [67–69]
to obtain an effective state vector consisting of an operator measured at a sequence of successive times. We reformulate
the linear model underpinning DMD in terms of these observable-vectors to approximate the system dynamics.
Takens’ embedding theorem [67, 68] establishes a connection between the manifold of states, which an observer
cannot directly access, and time-delayed measurements of an observable. In particular, the theorem asserts, under
 18

generous conditions, that a state on an N -dimensional (sub)manifold can be completely determined using a sequence
of at most d⋆ ≤ 2N + 1 time-delayed observables. This correspondence reads
 
o(t)
o(t + ∆τ )
 
 
ϕ(t) ↔ ot,d⋆ = 
 .. ,
 (B3)
 . 
o(t + (d∗ − 1)∆τ )

where ∆τ > 0 is the time delay, ϕ(t) is the system state, o(t) = o[ϕ(t)] is the measured observable, and ot is the
d⋆ -dimensional “observable trajectory” containing the dynamical information. The RHS of Eq. (B3) is known as a
d⋆ -dimensional delayed embedding of the observable. Takens’ theorem relates the evolution of microscopic degrees of
freedom to the evolution history of macroscopic observables, providing a concrete probe into the dynamical properties
of the system without direct access to the full states. Here we adopt the term Takens’ embedding technique to refer to
the method of applying time delays on the system observables, motivated by the rigorous results of Takens’ embedding
theorem.
In anticipation of efficiently leveraging near-term quantum resources, we choose the time delay in Takens’ embedding
technique to equal the DMD time interval, i.e., ∆τ = ∆t. Given this choice, we then measure the system along time
steps {tk = k∆t}K+1 K+1
k=0 and acquire the sequence of observable trajectories {otk ,d }k=0 , each of some length d ≤ d⋆ ,
 
o(tk )
o(tk+1 )
 
 
otk ,d =
 .. ,
 0 ≤ k ≤ K + 1. (B4)
 . 
o(tk+d−1 )

By construction, the first (d − 1) entries of otk ,d are identical to the last (d − 1) entries of otk−1 ,d . Consequently, the
matrix assembled by arranging successive trajectories otk as columns
h i
Xk1 :k2 = otk1 ,d otk1 +1 ,d · · · otk2 ,d , (B5)

has a Hankel form, i.e., the matrix elements on each anti-diagonal are equal. In the embedding space, we can identify
the closest linear flow,
LS
X1:K+1 = AX0:K =⇒ A = X1:K+1 (X0:K )+ , (B6)

where + denotes the Moore–Penrose pseudo-inverse. Here the system matrix A assumes a companion structure with
just d free parameters. The approximation to the system dynamics is then stored in the d parameters inferred from
measurements of K + d + 1 delayed observables. We hence name our least-squares embedding in the observable space
the observable dynamic mode decomposition (ODMD).

Appendix C: Multi-observable dynamic mode decomposition

Recall that the multi-observable system matrix has a characteristic polynomial,

 
Xd
CA (z) = det(z − A) = det  z ℓ Aℓ  , (C1)
ℓ=0

where Ad ≡ IdI . We examine the simpler version of Eq. (13) in which we assume the d matrix blocks,
 
Aℓ,11
Aℓ,22
 
 
Aℓ = Diag 
..
, (C2)
.

 
Aℓ,II
 19

to be diagonal with Aℓ,ij ≡ 0 for i ̸= j. In this case, the multivariate residual is

 
d d N −1 N −1 d N −1 I d
−−→ X X X X X X X X
Resk (A) = Aℓ⃗sℓ+k = Aℓ ⃗cn λℓ+k
n = λkn λℓn Aℓ⃗cn = λkn cn,i  λℓn Aℓ,ii  ⃗ei , (C3)
ℓ=0 ℓ=0 n=0 n=0 ℓ=0 n=0 i=1 ℓ=0
 
N
X −1 I
X I
X N
X −1
= λkn cn,i CA|i (λn )⃗ei =  cn,i λkn CA|i (λn ) ⃗ei , (C4)
n=0 i=1 i=1 n=0

where Eq. (C3) directly follows from the submatrices Aℓ being diagonal and ⃗ei denotes the canonical basis vector of
CI . Moreover, we define in Eq. (C4) the I single-observable characteristic polynomials,
d
X d−1
Y 
CA|i (z) = z ℓ Aℓ,ii = z − λ̃ℓ|i , (C5)
ℓ=0 ℓ=0

each with d roots λ̃ℓ|i (A0,ii , A1,ii , · · · , Ad−1,ii ). Clearly, Eq. (C4) indicates that the multivariate residual admits I
−−→
decoupled components corresponding to the different observables. Minimizing the total residual, ∥Resk (A)∥2 , is hence
PN −1 k
equivalent to minimizing each of the I observable residuals, n=0 cn,i λn CA|i (λn ) . This is consistent with the fact
that the matrix determinant det(z − A) factorizes into independent contributions,
 
YI Xd YI
CA (z) =  z ℓ Aℓ,ii  = CA|i (z), (C6)
i=1 ℓ=0 i=1

such that the eigenvalues of A also factorize into clusters based on single-observable residuals. The resulting MODMD
estimates are bounded by the I single-observable ODMD estimates, e.g.,
min |Ẽi,0 − E0 | ≤ |Ẽ0 − E0 | ≤ max |Ẽi,0 − E0 |, (C7)
1≤i≤I 1≤i≤I

where Ẽ0 and Ẽi,0 designate, respectively, the ground state energy estimate using the entire observable pool (the full
system matrix A) and a single observable (one row of the system matrix Aℓ,ii ).
More interesting convergence arises when the I observable residuals are coupled and give a total residual,
 
N −1 I d
−−→ X X X
Resk (A) = λkn cn,j  λℓn Aℓ,ij  ⃗ei , (C8)
n=0 i,j=1 ℓ=0
 
I
X N
X −1 N
X −1 d−1
X I
X
=  cn,i λkn CA|i (λn ) + λkn λℓn Aℓ,ij cn,j  ⃗ei , (C9)
i=1 n=0 n=0 ℓ=0 j̸=i

where we can lower the residual by utilizing the flexibility of off-diagonal elements in the submatrices Aℓ . For example,
the total residual may vanish completely when the multi-observable signal is d(I −1)-sparse in the eigenfrequency basis,
where, for 1 ≤ i ≤ I, the coefficients cn,i are supported on at most d(I−1) eigenindices Ni = {0 ≤ n ≤ N −1 : cn,i ̸= 0}.
That is, max1≤i≤I |Ni | ≤ d(I − 1). In this case, a vanishing residual is possible, provided that (a) |∪Ii=1 Ni | ≤ d(I − 1)
and (b) the matrix elements {Aℓ,ij }ℓ,j̸=i can be set appropriately such that ∀1 ≤ i ≤ I,
   

β0,i1 (1) . . . βℓ,ij̸=i (1) . . . βd−1,iI (1) A0,i1 β0,ii (1)CA|i λn1 [i]
   

 β0,i1 (2) . . . βℓ,ij̸=i (2) . . . βd−1,iI (2)   A0,i2   β0,ii (2)CA|i λn2 [i] 
.. ..   ..  ..
    
 ..  
 . . .  .   . 
   

 β0,i1 (|Ni |) . . . βℓ,ij̸=i (|Ni |) . . . βd−1,iI (|Ni |)   Aℓ,ij  = − β0,ii (|Ni |)CA|i λn [i]  , (C10)
    |Ni | 
 β (1) . . . β
ℓ,1j̸=i (1) . . . βd−1,1I (1)   0
    
 0,11   
 .. .. ..  . 
.
 .. 

 . . .  . 
   
 .


β0,I1 (|NI |) . . . βℓ,Ij̸=i (|NI |) . . . βd−1,II (|NI |) Ad−1,iI 0

where nξ [i] numerates support eigenindices in the set Ni for observable Oi , and βℓ,ij (ξ) = λℓnξ [i] cnξ [i],j . Given {Aℓ,ii }ℓ ,
a formal solution {Aℓ,ij }ℓ,j̸=i to Eq. (C10) exists if the matrix on the LHS has a full rank. This directly implies that the
 20

total residual may vanish completely for a dI-sparse multi-observable signal in the eigenfrequency basis (as expected),
because {Aℓ,ii }ℓ can always be chosen so that the single-observable characteristic polynomial CA|i (z) encompasses the
roots {λnξ [i] }ξ corresponding to at most d members of Ni . Without additional degrees of freedom from the off-diagonal
matrix elements of A, the residual only vanishes for a d-sparse signal.
We now relax the sparsity assumption on the signal to explore the general case. Observe that
2
K K X
I N −1 I N −1
X −−→ 2 X X X X
min Resk (A) = min cn,i λkn CA|i (λn ) + cn,j λkn CA|ij (λn ) , (C11)
A 2 A
k=0 k=0 i=1 n=0 j̸=i n=0
2
K X
X I X I
X X
≤ min cn,i λkn CA|i (λn ) + cn,j λkn C¯A|ij (λn ) , (C12)
{Aℓ,ii :0≤ℓ≤d−1,1≤i≤I } / i⋆
k=0 i=1 n∈N / i⋆
j̸=i n∈N

Pd Pd−1
where CA|ij (z) = ℓ=0 z ℓ Aℓ,ij ≡ ℓ=0 z ℓ Aℓ,ij for j ̸= i defines a degree-(d − 1) polynomial via the off-diagonal matrix
Pd−1
elements {Aℓ,ij }ℓ , Ni⋆ ⊂ [N ] is an eigenindex subset of size |Ni⋆ | = d(I − 1), and C¯A|ij (z) = ℓ=0 z ℓ Aℓ,ij (Ni⋆ ; {Aℓ,ii }ℓ )
fixes CA|ij (z) by evaluating specific coefficients {Aℓ,ij }ℓ such that,
I
X
cn,j C¯A|ij (λn ) = −cn,i CA|i (λn ), ∀n ∈ Ni⋆ . (C13)
j̸=i

Resembling Eq. (C10), we reserve the notation n⋆ξ [i] for the support eigenindices in Ni⋆ . Eq. (C13) then holds only if
the following coefficient matrix has full rank, i.e.,
 
⋆ ⋆ ⋆
β0,i1 (1) . . . βℓ,ij̸ = i (1) . . . βd−1,iI (1)
 
 β ⋆ (2) . . . β ⋆ ⋆
0,i1 ℓ,ij̸=i (2) . . . βd−1,iI (2) 
rank[βNi⋆ ] = rank   = d(I − 1), (C14)
 
.. .. ..

 . . . 

⋆
β0,i1 (|Ni⋆ |) . . . βℓ,ij̸
⋆ ⋆ ⋆ ⋆
=i (|Ni |) . . . βd−1,iI (|Ni |)

⋆
for βℓ,ij (ξ) = λℓn⋆ [i] cn⋆ξ [i],j . An arbitrary index set Ni⋆ can be assigned as long as the full-rank condition is met; ideally,
ξ
we aim for Ni⋆ = Ni⋆ ({Aℓ,ii }ℓ ) to satisfy the optimal property,

I
X X I X
X
cn,j λkn C¯A|ij (λn ) = min cn,j λkn C¯A|ij (λn ) , (C15)
N ⊂[N ]:|N |=d(I−1)
j=1 / i⋆
n∈N j=1 n∈N
/

which minimizes the effective residual parametrized by the d diagonal matrix elements {Aℓ,ii }ℓ , with C¯A|ii (z) = CA|i (z).
By the triangle inequality, we can bound the RHS of Eq. (C12) as,

X I
X X I
X X I
X
cn,i λkn CA|i (λn ) + cn,j λkn C¯A|ij (λn ) ≤ cn,j C¯A|ij (λn ) ≤ c̄i C̄ij , (C16)
/ i⋆
n∈N / i⋆
j̸=i n∈N / i⋆
j=1 n∈N j=1

P
where c̄i = max1≤j≤I n∈N ¯
/ i⋆ |CA|ij (λn )|. With suitable choices of the reference state |ϕ0 ⟩
/ i⋆ cn,j and C̄ij = maxn∈N
I
and operator pool {Oi }i=1 , we can strategically adjust the coefficients {cn,i }n,i , therefore controlling the constant c̄i
−1 ⋆
and the conditioning of the matrix βNi⋆ (or equivalently ∥βN ⋆ ∥2 ). Accordingly, for given index sets {Ni }i , we have
i

 2
K X
X I I
X
min c̄i C̄ij 
{Aℓ,ii :0≤ℓ≤d−1,1≤i≤I } k=0 i=1 j=1
 2 (C17)
I
X d−1 X
X I
≤ (K + 1) c̄2i min  max CA|i (λn ) + Aℓ,ij (Ni⋆ ; {Aℓ,ii }ℓ )  ,
/ i⋆
{Aℓ,ii :0≤ℓ≤d−1} n∈N
i=1 ℓ=0 j̸=i
 21

where the residual from the off-diagonal matrix elements can be bounded by,
d−1 X
X I
−1
Aℓ,ij (Ni⋆ ; {Aℓ,ii }ℓ ) ≤ d(I − 1) max⋆ |cn,i |∥βN ⋆ ∥2 max CA|i (λn ) ,
⋆
(C18)
n∈Ni i n∈Ni
ℓ=0 j̸=i

using an identity analogous to Eq. (C10). To further bound the RHS of Eq. (C17), we attempt to tightly approximate
the minimizer over the d diagonal matrix elements {Aℓ,ii }ℓ , or alternatively over the set of degree-d monic polynomials
defined along the circular arc Aθ = {z ∈ C : |z| = 1, −θ ≤ arg(z) ≤ θ}. In particular let us consider the complex-valued
Chebyshev polynomials,
ˆ =
C(z) arg min sup C(z) , (C19)
C monic:deg(C)=d z∈Aθ

whose parametric representations are explicitly constructible by tools such as Jacobi’s elliptic and theta functions [70].
For notational convenience, we denote ∥C∥θ = supz∈Aθ C(z) . Similar to the real-valued Chebyshev polynomials over
ˆ retains the minimal-norm property on Aθ , where ∥C∥
the interval [−1, 1], C(z) ˆ θ ∼ 2 sind (θ/2) cos2 (θ/4) asymptotically
for large d [71]. Setting θH (∆t) = (EN −1 − E0 )∆t/2, the total residual decays at least exponentially:
K I
" #2
X −−→ 2 X
min Resk (A) ≤ (K + 1)∥C∥ ˆ 2 c̄2i 1 + d(I − 1) max |cn,i |∥β −1⋆ ∥2 ,
θH (∆t) Ni (C20)
A 2 n∈Ni⋆
k=0 i=1

where the prefactor improves significantly, compared to the single-observable case derived in [9], when
X 1
max cn,j ≪ −1 , ∀1 ≤ i ≤ I. (C21)
1≤j≤I
/ ⋆
n∈N
1 + d(I − 1) max⋆ |cn,i |∥βN ⋆ ∥2
i n∈Ni i

Appendix D: Proof of theorems

1. Theorem 1

Theorem 1. Let Ẽ0 (d) be the approximate ground state energy extracted from the d-dimensional function subspace
d−1
spanned by {g1 , K∆t [g1 ], · · · , K∆t [g1 ]}. For d ≥ 1, there exists time step ∆t such that the error δEg (d) = Ẽ0 (d) − E0
is bounded by
sin[(EN −1 − E0 )∆t]
δE0 (d) ≤ 2(d−1)
tan2 Ξ, (D1)
ϵ̃0→1 ∆t
where cos2 Ξ = |⟨N f0 , g1 ⟩µ |2 = |⟨ψ0 |ϕ0 ⟩|2 is the squared overlap between the reference function g1 and the true ground
state K∆t -eigenfunction while ϵ̃0→1 = 1 + 3(E1 − E0 )∆t/2π ∈ [1, 2] characterizes the normalized spectral gap of the
Hamiltonian H.
[Proof.] Let us define,
⟨f, K∆t f ⟩µ
θ(f ) = arg ∈ [−π, π], (D2)
⟨f, f ⟩µ
which returns the expected dynamical phase factor associated with a function f under Koopman evolution, where f
N −1
belongs to the invariance subspace span{fn }n=0 . Given a suitable symmetrizing spectral shift such that ∥H∥ ∆t ≤ π,
we note that θ(f0 ) = maxf θ(f ) = −E0 ∆t. This variational principle implies,
1 1
Ẽ0 (d) = − max θ(f ) = − max θ(p(K∆t )g1 ), (D3)
∆t f ∈Krd ∆t p∈Pd−1
d−1
for which we use Krd to denote our d-dimensional function Krylov subspace span{g1 , · · · , K∆t [g1 ]} and Pd−1 the set
of degree-(d − 1) polynomials over C. By Eq. (44), we recall that g1 can be expanded in N orthogonal eigenfunctions
−1
{fn }N
n=0 , i.e.,

N
X −1
g1 = zn fn , (D4)
n=0
 22
PN −1 2
with n=0 |zn | = 1. Thus Eq. (D3) in the eigenbasis reads,

1 ⟨p(K∆t )g1 , K∆t p(K∆t )g1 ⟩µ

Ẽ0 (d) = − max arg , (D5)
∆t p∈Pd−1 ⟨p(K∆t )g1 , p(K∆t )g1 ⟩µ
PN −1 2
1 n=0 λn zn p(λn )
=− max arg P N −1 2 , (D6)
∆t p∈Pd−1 zn p(λn )
n=0

where λn = e−iE∆t label the eigenvalues of the Koopman operator. We note that
PN −1 2 2 PN −1 2
n=0 λn zn p(λn ) λ0 z0 p(λ0 ) + λN −1 n=1 zn p(λn )
max arg PN −1 2 ≥ max arg PN −1 2 , (D7)
p∈Pd−1 p∈Pd−1
n=0 zn p(λn ) n=0 zn p(λn )

by convexity of the circular sector with angle [−EN −1 ∆t, −E0 ∆t] in the complex plane. Then
 
PN −1 2
1 zn p(λ n )
Ẽ0 (d) ≤ − max arg λ0 + (λN −1 − λ0 ) Pn=1
N −1 2 ,
 (D8)
∆t p∈Pd−1 zn p(λ n )
n=0
 
PN −1 2
1 z n p(λ n )
=− arg λ0 + (λN −1 − λ0 ) min Pn=1
N −1 2 ,
 (D9)
∆t p∈Pd zn p(λn )
n=0

where we observe arg(λN −1 − λ0 ) = 21 [arg λN −1 + arg(−λ0 )] ≤ arg λ0 assuming E0 ≤ 0 WLOG. To proceed, we seek a
family of polynomials defined over the unit circle S1 = {z ∈ C : |z| = 1} to bound the fraction on the RHS of Eq. (D9).
Now let us fix some q ∈ (0, 1] and consider the handy choice of complex-valued Rogers-Szegő polynomials [72, 73],



 1 d=0

z + 1 d=1
Wd (z|q) = , (D10)
(1 + z)Wd−1 (z|q)


−(1 − q j−1 )zWd−2 (z|q) d ≥ 2



over the circle S1,q = {z : |z| = q −1/2 }. For simplicity, we rewrite z = −q −1/2 exp (−iϑ) where ϑ ∈ [−π, π]) denotes an
angular phase. Here a prefactor of −1 is included to periodically translate the polynomials so that Wd (ϑ|q) adapts the
symmetry Wd (−ϑ) = Wd (ϑ)∗ (we also omit a conditional dependence of Wd on q for notational clarity). Such family
of polynomials shares the key properties that (i) |Wd (ϑ)| remains bounded below unity over some angular window
W = [−Ω, Ω] ⊂ [−π, π] and (ii) |Wd (ϑ)| grows rapidly outside W.
Note that the constant q controls the width of our truncated angular window W. In the limit of q → 1, one can
verify that these polynomials converge to,

d  
X d  
Wd (ϑ) → exp −ik(ϑ + π) , (D11)
k
k=0

which simply gives the sum of evenly spaced points along S1 weighted by binomial coefficients. As a consequence,
supϑ |Wd (ϑ)| ≈ 2d for q ≈ 1. To bound the fraction in Eq. (D9) tightly, we look for a suitable linear transformation
L : S1 → S1 acting on the N individual eigenphases, λn = e−iϑn (∆t) , such that L nudges excited state angles ϑn≥1 all
inside the truncated window W while keeping the ground state angle ϑ0 outside. It is safe to assume ϑN −1 − ϑ1 ≤ 2Ω
and 2Ω ≤ ϑN −1 − ϑ0 ≤ π + Ω by adapting a suitable time step size ∆t, e.g.,
n o
∆t = sup τ ∈ R+ : ϑN −1 (τ ) − ϑ1 (τ ) ≤ 2Ω, ϑ1 (τ ) − ϑ0 (τ ) ≤ Ωc , (D12)
τ

with Ωc = π −Ω. Therefore a natural choice of L is the phase multiplication or angular translation, L : ϑ 7→ ϑ−Ω+ϑ1 ,
−1
which circularly moves {λn }N n=0 so that L(ϑ0 ) ≤ −Ω = L(ϑ1 ) ≤ −L(ϑN −1 ) as desired. With L chosen above, we
establish a variational upper bound on the ground state energy error by substituting the trial polynomials p = Wd−1 ◦L
 23

into Eq. (D9),

 
PQ−1 2
1 n=1 z n W d−1 (L(λ n ))
Ẽ0 (d) ≤ − arg λ0 + (λN −1 − λ0 ) PN −1 2 ,
 (D13)
∆t z n W d−1 (L(λ n ))
n=0
 
PQ−1 2
1 zn W d−1 (Ω)
≤− arg λ0 + (λN −1 − λ0 ) n=1 2
, (D14)
∆t z0 Wd−1 (L(λ0 ))
 
2
1 sin Ξ
=− arg λ0 + (λN −1 − λ0 ) 2
, (D15)
∆t Wd−1 (L(λ0 )) cos2 Ξ

where in arriving at Eq. (D15) we have utilized the property (i) of Wd−1 (ϑ) and defined the angle Ξ by cos2 Ξ =
|z0 |2 = |⟨N f0 , g1 ⟩µ |2 which specifies L2 -projection of our reference function g1 onto the ground state eigenfunction.
For the limiting case q = 1, it is rather straightforward to show that Ω = π/3 and,
1/d p
Wd (L(λ0 )) = 2 − 2 cos (ϑ1 − ϑ0 + Ω) ≥ 1 + Cϵ0→1 , (D16)

where ϵ0→1 = (ϑ1 − ϑ0 )/Ωc = (E1 − E0 )∆t/Ωc gives the normalized spectral gap times the time step and C is a
constant for which Eq. (D16) holds with ϵ0→1 ∈ [0, 1]. For example, C = 1 can be justified by concavity of the LHS
of the inequality above with respect to (ϑ1 − ϑ0 ). Hence we can further bound Eq. (D15) using Eq. (D16),
  
1 −2(d−1)

−2(d−1)
Ẽ0 (d) ≤ − arg λ0 1 − ϵ̃0→1 tan2 Ξ + λN −1 ϵ̃0→1 tan2 Ξ , (D17)
∆t
 
1 (1 − ζ) sin ϑ0 + ζ sin ϑN −1
≤ arctan , (D18)
∆t (1 − ζ) cos ϑ0 + ζ cos ϑN −1
sin(ϑN −1 − ϑ0 )
= E0 + ζ + O(ζ 2 ), (D19)
∆t
−2(d−1)
where we have defined ϵ̃0→1 = 1 + Cϵ0→1 > 1 and ζ = ϵ̃0→1 tan2 Ξ. The last equality can be derived from a Taylor
expansion up to leading order in ζ. So we have proved the statement as claimed. □

2. Aside: canonical angles

To set up the proof of Theorem 2, we first introduce the notion of subspace overlap for our subsequent discus-
sions. [74] Suppose that we have two subspaces X and Y of CN with dimX ≤ dimY. Let BX ∈ CN ×dimX and
BY ∈ CN ×dimY be the orthornormal basis matrices of X and Y (BX and BY are thus not unique). The canonical
angles between the subspaces are defined as
π
θℓ [X , Y] = arccos σℓ ∈ [0, ], 1 ≤ ℓ ≤ dimX , (D20)
2
†
where σ1 ≤ σ2 ≤ . . . ≤ σdimX denote the singular values of BX BY . For convenience, we use a diagonal matrix
 
θ1
θ2
 
 
Θ[X , Y] =  ..
 (D21)
.
 
 
θdimX

to record the set of dimX canonical angles. Notice when dimX = 1, the canonical angle is determined by the familiar
ℓ2 -inner product on the Hilbert space. Here we present elementary and known results [75] that are helpful for deriving
the block convergence bound.
Result C.2.1. For Y0 ⊆ Y with dimY0 = dimX , θℓ [X , Y] ≤ θℓ [X , Y0 ] and the equality can be saturated.
Result C.2.2. Suppose dimX = dimY and θℓ [X , Y] ̸= π2 .
 24

(a) For Y0 ⊆ Y, there exists unique X0 ⊆ X with dimX0 = dimY0 such that ΠY X0 = Y0 for orthogonal projection ΠY
onto Y. Moreover,

θℓ [X0 , Y0 ] ≤ θℓ [X , Y], (D22)

for 1 ≤ ℓ ≤ dimY0 and hence

∥sin Θ[X0 , Y0 ]∥ ≤ ∥sin Θ[X , Y]∥ (D23)

for any unitarily invariant norm ∥·∥.

(b) For any orthonormal vectors [y1 , . . . , yQ ] in Y, there exists linearly independent vectors [x1 , . . . , xQ ] in X so that
ΠY xℓ = yℓ and (a) holds for X0 = span{x1 , . . . , xQ } and Y0 = span{y1 , . . . , yQ } with dimX0 = dimY0 = Q.

3. Theorem 2

Theorem 2. Let Ẽn (d) be the approximate nth eigenenergy extracted from the dI-dimensional function subspace
d−1
span{⃗g , K∆t [⃗g ], · · · , K∆t [⃗g ]}, and δEn (d) = Ẽn (d) − En be the approximation error. Consider the error matrix,
 
δE0 (d)

∆I (d) = diag  .. 
, (D24)
 . 
δEI−1 (d)

which contains approximations to the lowest I energies. For d ≥ 1, there exists time step ∆t such that the spectral
approximation ∆I (d) is bounded by

sin[(EN −1 − E0 )∆t]
∥∆I (d)∥2 ≤ 2(d−1)
∥tan2 Θ∥2 , (D25)
ϵ̃I−17→I ∆t

for the operator norm ∥·∥2 . Here Θ denotes the canonical angle between the two subspaces span{fn : 0 ≤ n ≤ I − 1}
and span{gi : 1 ≤ i ≤ I}, which generalizes the squared overlap in Theorem 1 (see Appendices D 2 and D 3). In the
denominator, ϵ̃I−17→I ∈ [1, 2] depends on the Ith spectral gap (EI − EI−1 ) of the Hamiltonian H.
[Proof.] It suffices to establish the more general result for the error matrix ∆ = diag[δEn1 (d), δEn1 +1 (d), . . . , δEn2 (d)]
for n2 − n1 + 1 = I. From here on, we will use the operator norm as a specific example of a unitarily invariant norm by
setting ∥·∥ = ∥·∥2 . Our key contribution lies in extending the proof of the theorem, building on existing results for block
Krylov methods [74], to the real-time setting. Applying result C.2.2(b), we know there exists BI = [x1 , . . . , xI ] with
spanBI = span{g1 , . . . , gI } so that Πinvariant subspace BI = [fn1 , . . . , fn2 ]. WLOG we assume that BI 7→ BI (BI† BI )−1/2
has orthonormal columns, i.e., ⟨xi , xj ⟩µ = δij for 1 ≤ i, j ≤ I. For a degree-(d − 1) polynomial p, we define

N
X −1
RI = p(K∆t )BI = fn p(λn )fn† BI , (D26)
n=0

where the dual fn† are defined with respect to the functional inner product ⟨·, ·⟩µ . We will further write the equality
above as
h i
RI = F p(Λ)F † BI = Fa p(Λa )Fa† + Fb p(Λb )Fb† + Fc p(Λc )Fc† BI , (D27)

where F = [Fa Fb Fc ] and Λ = Λa ⊕ Λb ⊕ Λc contain the N Koopman eigenfunctions and eigenvalues respectively in a
block form. Here we use three subscripts a, b, and c to label the partition of Hamiltonian spectrum into three disjoint
blocks with energies below En1 , from En1 to En2 , and above En2 . In particular, the middle block b is represented by
Fb = [fn1 , . . . , fn2 ] and Λb = diag[λn1 , . . . , λn2 ].
By our construction, Fb† BI is nonsingular so
h i−1 h i−1 h i−1
RI (Fb† RI )−1 = RI p(Λb )Fb† BI = Fa p(Λa )Fa† BI p(Λb )Fb† BI + Fb + Fc p(Λc )Fc† BI p(Λb )Fb† BI , (D28)
 25

whenever p(Λb ) is nonsingular by suitable choice of the polynomial p. Since spanRI is clearly a subspace of the Krylov
d−1
subspace Krd,I := span{⃗g , · · · , K∆t [⃗g ]}, result C.2.1 implies
 
tan Θ spanFb , Krd,I ≤ tan Θ[spanFb , spanRI ] , (D29)

−1
= sin Θ[spanFb , spanRI ] cos Θ[spanFb , spanRI ] , (D30)
" #" #
† † −1 p(Λa ) 0 Fa† BI (Fb† BI )−1
= [Fa Fc ] RI (Fb RI ) = p(Λb )−1 , (D31)
0 p(Λc ) Fc† BI (Fb† BI )−1
" # " #
p(Λa ) 0 Fa† BI (Fb† BI )−1
≤ p(Λb )−1 , (D32)
0 p(Λc ) Fc† BI (Fb† BI )−1
max |p(λn )|
n<n1 ∨n>n2
≤ tan Θ[spanFb , spanBI ] , (D33)
min |p(λn )|
n1 ≤n≤n2

where Eq. (D31) results from the substitution of Eq. (D28).

The inequality above characterizes the convergence of the Krylov subspace (towards to the target eigenstates), and
our task is to choose a polynomial p that can tightly bound the fraction in Eq. (D33). As in Appendix D 1, we exploit
the Rogers-Szegő polynomials Wd−1 (z|q) and work in the limit q = 1 to derive relevant results. For n2 = N − 1, let
us consider the polynomial
Wd−1 (e−iφ z)
p(z) = , (D34)
Wd−1 (e−iφ λn1 )
where e−iφ is a constant phase offset shifting the angles ϑn<n1 (recall that ϑn = En ∆t) inside the window W = [−Ω, Ω]
while keeping the angles ϑn1 ≤n≤n2 outside W. In this case, observe that
min |p(λn )| = |p(λn1 )| = 1, (D35)
n1 ≤n≤n2
1
max |p(λn )| = max |p(λn )| ≤ . (D36)
n<n1 ∨n>n2 n<n1 Wd−1 (e−iφ λn1 )
By employing suitable time step size ∆t, e.g., that from Eq. (D12), we assume ϑn1 −1 −ϑ0 ≤ 2Ω and 2Ω ≤ ϑN −1 −ϑ0 ≤
−1
π +Ω. Specifically, φ = Ω−ϑn1 −1 circularly shifts the eigenphases {λn }N
n=0 . This ensures that |φ + ϑ0 | ≤ φ+ϑn1 −1 =
Ω ≤ φ + ϑn1 as intended. From Theorem 1, we recall Ω = π/3 and |Wd−1 (e−iφ λn1 )|1/(d−1) ≥ 1 + Cϵn1 −17→n1 where
ϵn1 −17→n1 = 3(En1 − En1 −1 )∆t/2π encodes the normalized n1 th energy gap and C is some appropriate constant that
can be set to 1 (c.f. Eq. (D16)). This immediately implies an exponential convergence to the eigenspaces,
  1  
tan Θ spanFb , Krd,I ≤ (d−1) tan Θ spanFb , span{g1 , . . . , gI } , (D37)
ϵ̃n1 −17→n1
for the base case n2 = N − 1. When n2 < N − 1, we consider the polynomial p = p1 · p2 with
Wd−N +n2 (e−iφ z)
p1 (z) = , (D38)
Wd−N +n2 (e−iφ λn1 )
and
N −1
Y z − λn
p2 (z) = , (D39)
n=n2
λ
+1 n2
− λn

where p is factorized into lower order polynomials p1 of degree (d − N + n2 ) and p2 of degree (N − n2 − 1). Observe
that by design p(λn ) = 0 for n ≥ n2 + 1 and thus
N −1
Y λn − λm
min |p(λn )| ≥ min |p1 (λn )| min p2 (λn ) = min , (D40)
n1 ≤n≤n2 n1 ≤n≤n2 n1 ≤n≤n2 n1 ≤n≤n2
m=n2 +1
λ n2 − λ m
N −1
1 Y λn − λm
max |p(λn )| ≤ max |p1 (λn )| max |p2 (λn )| ≤ max , (D41)
n<n1 ∨n>n2 n<n1 n<n1 Wd−N +n2 (e−iφ λn1 ) n<n1 m=n +1 λn2 − λm
2
 26

which directly follows from applying our base case result to the eigensector {λn }nn=0
2
. Note that our argument requires
d ≥ N − n2 : a symmetric bound can be readily established using the exact same argument with the roles of n1 and
n2 exchanged, i.e., for d ≥ n1 + 1, upon a spectral flip H 7→ −H.
Now we proceed to analyze the convergence of the Krylov energy approximation. For simplicity, we focus again on
approximating from below Hamiltonian eigenvalues En ≲ EN −1 near the right edge of the spectrum. The argument
for approximating from above eigenvalues near the left edge of the spectrum can be easily adapted with a spectral flip.
We observe that up to the trivial spectral rotation K∆t 7→ λ−1n2 K∆t leaving the functional Krylov subspace invariant,
we can assume WLOG En2 −1 ≤ 0 < En2 +1 for the remainder of the proof. For V0 := RI (Fb† RI )−1 , we first construct
V− = V0 (V0† V0 )−1/2 , which contains orthonormal columns. This allows us to introduce eigenvalues, λ̂n = |λ̂n |e−iÊn ∆t
for n1 ≤ n ≤ n2 , of the spanV0 -projected Koopman operator V−† K∆t V− , with the ordering arg(λ̂n1 ) ≥ . . . ≥ arg(λ̂n2 ).
These eigenvalues depend explicitly on choices of the degree-(d − 1) polynomial p (as RI does). Since spanV0 ⊆ Krd,I
irrespective of p, the variational principle gives

Ên (d) ≤ Ẽn (d) ≤ En =⇒ 0 ≤ δEn (d) ≤ En − Ên (d), n1 ≤ n ≤ n2 . (D42)

The variational characterization, as we shall see, leads to an upper bound on the spectral error
 
δEn1
∆I,n1 ,n2 = diag 
 .. 
, (D43)
 . 
δEn2

which generalizes ∆I in Theorem 2 (∆I,n1 ,n2 = ∆I when n1 = N −I and n2 = N −1). With our previous assumption
En2 = 0, we have Ên ≤ 0 for n1 ≤ n ≤ n2 . Consequently, ∀x ∈ CI

1 1 ⟨V0 x, K∆t V0 x⟩µ

− θ(V0 x) = − arg ≤ 0, (D44)
∆t ∆t ⟨V0 x, V0 x⟩µ

where we adopt the same notation as in Theorem 1 (c.f. Eq. (D2)). By unfolding V0 via Eq. (D28), we arrive at

V0† K∆t V0 = G†a Λa Ga + Λb + G†c Λc Gc , (D45)

and

V0† V0 = G†a Ga + IdI + G†c Gc , (D46)

for Ga = p(Λa )Fa† BI [p(Λb )Fb† BI ]−1 and Gc = p(Λc )Fc† BI [p(Λb )Fb† BI ]−1 . Substituting Eqs. (D45) and (D46) into the
variational inequality, we can derive

1 x† (G†a Λa Ga + ηb λn2 IdI )x 1 x† Λb x 1 1 x† (G†a Λa Ga + Λb )x

− arg †
− arg † + argλn2 ≤ − arg , (D47)
∆t x x ∆t x x ∆t ∆t x† x
1 x† (G†a Λa Ga + Λb + G†c Λc Gc )x
≤ − arg , (D48)
∆t x† (G†a Ga + IdI + G†c Gc )x
1 ⟨V0 x, K∆t V0 x⟩µ
=− arg ≤ 0, (D49)
∆t ⟨V0 x, V0 x⟩µ

where ηb = inf x∈CI :x† x=1 ∥x† Λb x∥ is a constant of O(1) when ϑn2 −ϑn1 < π/2, and we recall that En2 = 0 =⇒ λn2 = 1
in our current setting. In particular, the inequality Eq. (D47) above holds since for any unit vector x 7→ √ x† ,
x x

λn2 λn2 λ n2 λn2

† + +
λn x (G†a Λa Ga
+ Λb )x † †
x Λb x x Ga Λa Ga x
† † †
x Λb x x Ga Λa Ga x
†
arg † † 2 = arg ≤ arg ≤ 0, (D50)
x (Ga Λa Ga + ηb λn2 IdI )x · x† Λb x 1+
ηb λn2 ∥x† Λb x∥λn2
1+
x† G†a Λa Ga x x† G†a Λa Ga x

where the numerator has a smaller phase than the denominator. The inequality Eq. (D49), on the other hand, follows
from our basic assumption that En ≤ 0 for n ≤ n2 . The inequalities indicate that the G†a Λa Ga -eigenvalues, denoted as
 27

λ̄n for n1 ≤ n ≤ n2 (with the eigenvalues ordered by decreasing arguments), are related to the V−† K∆t V− -eigenvalues
through
1
− arg(λ̄n + ηb λn2 ) + En − En2 ≤ Ên , n1 ≤ n ≤ n2 . (D51)
∆t
Consequently, we may relax the error bound in Eq. (D42),
1
δEn ≤ En − Ên ≤ arg(λ̄n + ηb λn2 ) + En2 , (D52)
∆t " #
1 e†n G†a Ga en
≤ arg λ0 + λn2 + En2 , (D53)
∆t ηb e†n en
" #
1 yn† G†a Ga yn
≤ arg λ0 + λn2 + En2 , (D54)
∆t ηb yn† yn

where en denotes the eigenvector of G†a Λa Ga corresponding to λ̄n and yn the eigenvector of G†a Ga corresponding to
the (n − n1 + 1)th largest eigenvalue. Importantly, we remark that Eqs. (D47) and (D54) can be established without
phase ambiguity of 2πZ by choosing suitable time step size ∆t.
We therefore seek a tight bound on the eigenvalues of G†a Ga (or equivalently the singular values of Ga ). Notice that
the singular values of p(Λa )−1 · Ga · p(Λb ) = p(Λa )−1 · p(Λa )Fa† BI [p(Λb )Fb† BI ]−1 · p(Λb ) = Fa† BI [Fb† BI ]−1 are bounded
above by the singular values of the canonical tangents (c.f. Eq. (D33))
" #
Fa† BI (Fb† BI )−1
tan Θ[spanFb , spanBI ] = . (D55)
Fc† BI (Fb† BI )−1

This allows us to effectively control the spectral radius of G†a Ga ,

∥G†a Ga ∥ ≤ ∥p(Λa )∥2 ∥p(Λb )−1 ∥2 tan2 Θ[spanFb , spanBI ] , (D56)

max |p(λn )|2
n<n1
= tan2 Θ[spanFb , spanBI ] , (D57)
min |p(λn )|2
n1 ≤n≤n2

where it suffices to identify a degree-(d−1) polynomial p that tightly bounds the fraction on the RHS of the expression
above. In addition, here p must satisfy the key orthogonality constraints,

f˜n† RI = [0 0 . . . 0], n2 + 1 ≤ n ≤ N − 1, (D58)

where f˜n ∈ Krd,I denotes the approximate eigenfunction corresponding to the Krylov eigenvalue λ̃n = |λ̃n |e−iẼn ∆t .
The set of (N − n2 ) constraints ensure that spanRI = p(K∆t )spanBI = p(K∆t )span{g1 , . . . , gI } belongs to the correct
eigensector of the Krylov subspace. Let us consider the familiar polynomial p = p1 · p2 with p1 of degree-(d − N + n2 )
and p2 of degree (N − n2 − 1) given by

Wd−N +n2 (e−iφ z)

p1 (z) = , (D59)
Wd−N +n2 (e−iφ λn1 )
and
N −1
Y z − λ̃n
p2 (z) = , (D60)
n=n2 +1 n2 − λ̃n
λ

as introduced in Eq. (D38) and Eq. (D39) respectively. Accordingly, we can generalize Eq. (D17) from Theorem 1,
" #
1 ζn1 ,n2 sin ϑ0 + sin ϑn2
δEn ≤ − arctan + En2 , (D61)
∆t ζn1 ,n2 cos ϑ0 + cos ϑn2
sin(ϑn2 − ϑ0 )
= ζn1 ,n2 + O(ζn21 ,n2 ), (D62)
∆t
 28

−2(d−N +n2 )
where ζn1 ,n2 = ηb−1 ι2n1 ,n2 ϵ̃n1 −1→n1 ∥tan2 Θ∥ for



 1 when n2 = N − 1,
N −1




 max
 Y λn − λ̃m
ιn1 ,n2 = n<n1
m=n2 +1 λn2 − λ̃m (D63)

 N −1
otherwise.

 Y λn − λ̃m
min



m=n2 +1 λn2 − λ̃m
 n
 1 ≤n≤n 2

This immediately implies

ι2n1 ,n2 sin[(En2 − E0 )∆t] ∥tan2 Θ∥

∥∆I,n1 ,n2 ∥ ≤ 2(d−N +n )
, (D64)
ηb ϵ̃n1 −1→n12 ∆t

which simplifies to the desired result when n2 = N − 1 as claimed. □

Appendix E: Additional simulation details

1. Reference states

Below we provide the exact reference states |ϕ0 ⟩ employed for the TFIM and LiH calculations. The TFIM reference
is an equal superposition of 6 bitstring states,
1
|ϕ0 ⟩TFIM = √ (|000000000000000⟩ + |111111111111111⟩ + |100000000000000⟩
6
+ |000000001111111⟩ + |000000011111111⟩ + |000000111111111⟩),

while the LiH reference is an equal superposition of 4 Hartree-Fock states,

1
|ϕ0 ⟩LiH = (|0000100001⟩ + |0000000011⟩ + |0001000001⟩ + |0010100000⟩),
2
where we work with the STO-3G basis set (2 core and 2 valence electrons) to represent the second-quantized molecular
Hamiltonian. Note that both reference states are superpositions of a small number of computational basis states and
thus have efficient classical representation.

2. Convergence of eigenenergies

In order to better illustrate the regime division in error convergence, as explained in Section V, we provide alternative
visualizations of Fig. 2a and Fig. 4a. Fig. 6 shows the absolute energy error as a function of time steps on the log-log
scale for a larger range of K values, highlighting the transition from the exponential to algebraic convergence. As the
number of time steps increases, the log-log plot demonstrates the linear behavior expected for algebraic error decay.

Appendix F: Prediction of system dynamics

To demonstrate the power of MODMD algorithm beyond eigenenergy estimation, we provide the results of additional
experiments where we construct the system matrix A from a set of observable measurements and then use it to predict
future system dynamics. As per Eq. (12) and Eq. (13), we measure the multi-observable expectations ⃗sk = ⃗s(k∆t) for
k = 0, 1, . . . , k ∗ and arrange these real-time snapshots to construct X and X′ . Here k ∗ is a hyperparameter controlling
how many snapshots are collected and therefore how much past information is used to predict future dynamics. After
constructing the system matrix A, we follow the MODMD prescription ∗
⃗sk+1 ≈ A⃗sk as outlined in Section III D. This
allows us to predict ⃗sk for k = k ∗ + 1, k ∗ + 2, . . . by computing Ak−k X′ . In Fig. 7, we show the real part of the ith
element of ⃗sk , i.e., ℜ⟨ϕ0 |Oi |ϕ0 (k∆t)⟩, up to 200 time steps beyond k ∗ . We observe that for both the transverse-field
Ising and lithium hydride Hamiltonians, the predicted system dynamics become more accurate as k ∗ increases. These
 29

101 | E0| 101 | E0|

100 | E1| 100 | E1|
| E2| 10 1 | E2|
10 1 | E3| | E3|
Absolute Error

Absolute Error
noise 10 2 noise
10 2 1/K 1/K
10 3
10 3
10 4
10 4 10 5
10 5 10 6
10 6 10 7
101 102 103 101 102 103
K ( Simulation Time) K ( Simulation Time)
(a) TFIM MODMD convergence (b) LiH MODMD convergence

FIG. 6: Convergence of eigenenergies for TFIM and LiH. To obtain eigenenergy estimates Ẽn , we fix the MODMD parameters
K
d
= 52 and δ̃ = 10−2 for constructing and thresholding the data matrix pair X, X′ ∈ RdI×(K+1) . Gaussian N (0, ε2noise ) noise
with εnoise = 10−3 is added independently to distinct matrix elements. The absolute errors, |Ẽn − En |, in the first four lowest
target eigenenergies of the TFIM and LiH Hamiltonians are shown with respect to K proportional to the non-dimensional
maximal simulation time. The reference states |ϕ0 ⟩ are sparse in the computational basis as prescribed in Appendix E 1. The
shading above the solid lines indicates standard deviation across repeated trials for each quantity. Left. TFIM absolute errors
from the MODMD algorithm with ∆t ≈ 0.08 and I = 6 (plotted on a log-log scale), where the convergence results are averaged
over 5 trials, each involving a Gaussian noise realization and a random selection of I observables. Right. LiH absolute errors
from the MODMD algorithm with ∆t ≈ 0.33 and I = 6 (on a log-log scale), where the results are averaged over 20 trials, each
involving a Gaussian noise realization.

results align with our basic intuition that collecting measurements over longer simulation time should grant the ability
to accurately predict dynamics further into the future. Though we only predict the dynamics for the particular sets of
observables used in Section V, one could perform similar signal extrapolation for other quantities of interest depending
on the system and dynamics in question.

0.05 True Signal

MODMD Predicted Signal

k*=500
0.00
0.05
0 25 50 75 100 125 150 175 200
0.05
Im[ 0|Oi| 0(k t) ]

0.00 k*=1500
0.05
0 25 50 75 100 125 150 175 200
0.05
k*=2500

0.00
0.05
0 25 50 75 100 125 150 175 200
k k*
(a) Dynamics of TFIM
 30

0.2 True Signal

MODMD Predicted Signal

k*=100
0.0

0.2
0 25 50 75 100 125 150 175 200
0.2
Re[ 0|Oi| 0(k t) ]

k*=500
0.0

0.2
0 25 50 75 100 125 150 175 200
0.2

k*=1000
0.0

0.2
0 25 50 75 100 125 150 175 200
k k*
(b) Dynamics of LiH

FIG. 7: Predicted observable dynamics from MODMD. For different values of k∗ (as specified on the right side of each panel),
we construct the system matrix A from the multi-observable snapshots ⃗sk = ⃗s(k∆t) for k = 0, 1, . . . , k∗ . We then predict ⃗sk for
∗
k = k∗ + 1, k∗ + 2, . . . , k∗ + 200 by computing ⃗sk = Ak−k X′ and plot the ith element of predicted ⃗sk , i.e., ℜ⟨ϕ0 |Oi |ϕ0 (k∆t)⟩.
The figure shows the results for i = 5, though similar convergence behavior is observed for all observables. The MODMD
parameters K d
= 52 and δ̃ = 10−2 are fixed for constructing and thresholding the pair of data matrices X, X′ ∈ RdI×(K+1) .
Gaussian N (0, ε2noise ) noise with εnoise = 10−3 is added independently to the matrix elements. The reference states, time steps,
and observable sets are the same as those used in Section V. The results are averaged over 20 trials, each involving a Gaussian
noise realization. Top. Dynamics predictions for the transverse-field Ising model (TFIM) Hamiltonian. Bottom. Dynamics
predictions for the lithium hydride (LiH) Hamiltonian.
 31

[1] R. M. Parrish and P. L. McMahon, arXiv preprint [27] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proc.
10.48550/arXiv.1909.08925 (2019). Roy. Soc. Lond. A 454, 339 (1998).
[2] N. H. Stair, R. Huang, and F. A. Evangelista, J. Chem. [28] V. Havlı́ček, A. D. Córcoles, K. Temme, A. W. Harrow,
Theory Comput. 16, 2236 (2020). A. Kandala, J. M. Chow, and J. M. Gambetta, Nature
[3] K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, 567, 209 (2019).
M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, [29] H. H. S. Chan, R. Meister, M. L. Goh, and B. Koczor,
J. Moussa, N. Wiebe, et al., PRX Quantum 3, 020323 arXiv preprint 10.48550/arXiv.2212.11036 (2023).
(2022). [30] A. Moitra, in Proceedings of the Forty-Seventh Annual
[4] Y. Shen, K. Klymko, J. Sud, D. B. Williams-Young, ACM Symposium on Theory of Computing, STOC ’15
W. A. d. Jong, and N. M. Tubman, Quantum 7, 1066 (Association for Computing Machinery, New York, NY,
(2023). USA, 2015) p. 821–830.
[5] C. L. Cortes and S. K. Gray, Phys. Rev. A 105, 022417 [31] W. Li, W. Liao, and A. Fannjiang, IEEE Trans. Inform.
(2022). Theory 66, 4593 (2020).
[6] N. H. Stair, C. L. Cortes, R. M. Parrish, J. Cohn, and [32] Z. Ding, E. N. Epperly, L. Lin, and R. Zhang, arXiv
M. Motta, Phys. Rev. A 107, 032414 (2023). preprint 10.48550/arXiv.2404.03885 (2024).
[7] Z. Ding and L. Lin, PRX Quantum 4, 020331 (2023). [33] H.-Y. Hu, S. Choi, and Y.-Z. You, Phys. Rev. Res. 5,
[8] Z. Ding and L. Lin, Quantum 7, 1136 (2023). 023027 (2023).
[9] Y. Shen, D. Camps, A. Szasz, S. Darbha, K. Klymko, [34] A. A. Akhtar, H.-Y. Hu, and Y.-Z. You, Quantum 7,
D. B. Williams-Young, N. M. Tubman, and R. Van Beeu- 1026 (2023).
men, arXiv preprint 10.48550/arXiv.2306.01858 (2023). [35] H.-Y. Hu, A. Gu, S. Majumder, H. Ren, Y. Zhang, D. S.
[10] G. Wang, D. S. França, R. Zhang, S. Zhu, and P. D. Wang, Y.-Z. You, Z. Minev, S. F. Yelin, and A. Seif,
Johnson, Quantum 7, 1167 (2023). arXiv preprint 10.48550/arXiv.2402.17911 (2024).
[11] Z. Ding, H. Li, L. Lin, H. Ni, L. Ying, and R. Zhang, [36] M. Ippoliti, Y. Li, T. Rakovszky, and V. Khemani, Phys.
Quantum 8, 1487 (2024). Rev. Lett. 130, 230403 (2023).
[12] Z. Zhang, A. Wang, X. Xu, and Y. Li, Quantum 8, 1438 [37] T. Schuster, J. Haferkamp, and H.-Y. Huang, arXiv
(2024). preprint 10.48550/arXiv.2407.07754 (2024).
[13] M. Motta, W. Kirby, I. Liepuoniute, K. J. Sung, J. Cohn, [38] N. Gomes, J. Yin, S. Niu, C. Yang, and W. A. de Jong,
A. Mezzacapo, K. Klymko, N. Nguyen, N. Yoshioka, and arXiv preprint (2023).
J. E. Rice, Electronic Structure 6, 013001 (2024). [39] A. F. Kemper, C. Yang, and E. Gull, Phys. Rev. Lett.
[14] N. Maskara, S. Ostermann, J. Shee, M. Kalinowski, A. M. 132, 160403 (2024).
Gomez, R. A. Bravo, D. S. Wang, A. I. Krylov, N. Y. Yao, [40] R. Kaneko, M. Imada, Y. Kabashima, and T. Ohtsuki,
M. Head-Gordon, M. D. Lukin, and S. F. Yelin, arXiv arXiv preprint (2024).
preprint https://doi.org/10.48550/arXiv.2312.02265 [41] D. Chandler, Introduction to Modern Statistical Mechan-
(2023). ics (Oxford University Press, 1987).
[15] E. N. Epperly, L. Lin, and Y. Nakatsukasa, SIAM J. Ma- [42] D. Limmer, Statistical Mechanics and Stochastic Ther-
trix Anal. Appl. 43, 1263 (2022). modynamics: A Textbook on Modern Approaches in and
[16] H. Li, H. Ni, and L. Ying, Phys. Rev. A 108, 062408 Out of Equilibrium, Oxford Graduate Texts (Oxford Uni-
(2023). versity Press, 2024).
[17] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, [43] B. O. Koopman, Proc. Natl. Acad. Sci. USA 17 5, 315
and H. Neven, Nat. Commun. 9, 4812 (2018). (1931).
[18] A. Arrasmith, Z. Holmes, M. Cerezo, and P. J. Coles, [44] G. D. Birkhoff and B. O. Koopman, Proceedings of the
Quantum Sci. Technol. 7, 045015 (2022). National Academy of Sciences 18, 279 (1932).
[19] H. Qi, L. Wang, H. Zhu, A. Gani, and C. Gong, Quantum [45] S. Brunton, B. Brunton, J. Proctor, and J. Kutz, PloS
Inf. Process. 22, 435 (2023). one 11 (2015).
[20] H.-Y. Huang, R. Kueng, and J. Preskill, Nat. Phys. 16, [46] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, J.
1050 (2020). Nonlinear Sci. 25, 1307 (2015).
[21] H.-Y. Hu and Y.-Z. You, Phys. Rev. Res. 4, 013054 [47] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Proceedings
(2022). of the National Academy of Sciences 113, 3932 (2016).
[22] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, [48] H. Arbabi and I. Mezić, SIAM J. Appl. Dyn. Syst. 16,
J. Preskill, B. Vermersch, and P. Zoller, Nat. Rev. Phys. 2096 (2017).
5, 9 (2023). [49] G. Prony, Journal de l’Ecole Polytechnique (Parı́s) 1, 24
[23] C. Bertoni, J. Haferkamp, M. Hinsche, M. Ioannou, (1795).
J. Eisert, and H. Pashayan, Phys. Rev. Lett. 133, 020602 [50] C. C. Paige (1971).
(2024). [51] Y. Saad, SIAM J. Numer. Anal. 17, 687 (1980).
[24] K. Van Kirk, J. Cotler, H.-Y. Huang, and M. D. Lukin, [52] R. A. Horn and C. R. Johnson, Matrix Analysis (Cam-
arXiv preprint 10.48550/arXiv.2212.06084 (2022). bridge University Press, 1985).
1
PK −iHk∆t
[25] A. Zhao, N. C. Rubin, and A. Miyake, Phys. Rev. Lett. [53] limK→∞ K+1 k=0 f (e |ϕ0 ⟩) converges uniformly
127, 110504 (2021). R
to H dµ̄f for all |ϕ0 ⟩ ∈ H, provided that the limit mea-
[26] Z. Liu, Z. Hao, and H.-Y. Hu, arXiv preprint sure µ̄ is uniquely invariant. For quantum systems, the
10.48550/arXiv.2311.00695 (2023). action of dynamical evolution can be regarded as a col-
 32

lection of independent rotations up to a change of basis. [64] P. J. Schmid, J. Fluid Mech. 656, 5–28 (2010).
If the Hamiltonian spectrum and evolution time step re- [65] I. Mezić, Annu. Rev. Fluid Mech. 45, 357 (2013).
main general, µ̄ is then the uniform measure on H as the [66] While “observable” is typically used in quantum mechan-
unique invariant measure. ics to refer specifically to Hermitian operators (with real
[54] Y. Saad, Numerical Methods for Large Eigenvalue Prob- expectation value), here we use a broader definition, en-
lems (Society for Industrial and Applied Mathematics, compassing also complex scalar quantities that can be
2011). computed from measurements on a quantum computer.
[55] M. Suzuki, J. Math. Phys. 32, 400 (1991). [67] D. Ruelle and F. Takens, Comm. Math. Phys. 20, 167
[56] N. Hatano and M. Suzuki, Finding exponential prod- (1971).
uct formulas of higher orders, in Quantum Annealing [68] F. Takens, in Lecture Notes in Mathematics (Springer
and Other Optimization Methods, edited by A. Das and Berlin Heidelberg, 1981) pp. 366–381.
B. K. Chakrabarti (Springer Berlin Heidelberg, Berlin, [69] I. L. Gutiérrez, F. Dietrich, and C. B. Mendl, Quant. Inf.
Heidelberg, 2005) pp. 37–68. Proc. 22, 251 (2023).
[57] A. M. Childs and N. Wiebe, Quantum Info. Comput. 12, [70] N. Akhiezer, Izv. Kaz. fiz. mat. ob-va 3 (1928).
901–924 (2012). [71] K. Schiefermayr, Acta Scientiarum Mathematicarum 85,
[58] Y. Nakata, C. Hirche, M. Koashi, and A. Winter, Phys. 629 (2019).
Rev. X 7, 021006 (2017). [72] C. Schmid and K. J. DeMars, Mathematics 8,
[59] E. Kökcü, D. Camps, L. Bassman Oftelie, J. K. Freericks, 10.3390/math8020171 (2020).
W. A. de Jong, R. Van Beeumen, and A. F. Kemper, [73] M. Foupouagnigni and W. Koepf, Orthogonal Polynomi-
Phys. Rev. A 105, 032420 (2022). als: 2nd AIMS-Volkswagen Stiftung Workshop, Douala,
[60] D. Camps, E. Kökcü, L. Bassman Oftelie, W. A. de Jong, Cameroon, 5-12 October, 2018 (2020).
A. F. Kemper, and R. Van Beeumen, SIAM J. Matrix [74] R.-C. Li and L.-H. Zhang, Numer. Math. 131, 83 (2015).
Anal. Appl. 43, 1084 (2022). [75] G. W. Stewart and J.-g. Sun, Matrix perturbation theory,
[61] I. Mezić and A. Banaszuk, Physica D: Nonlinear Phe- Computer science and scientific computing (Academic
nomena 197, 101 (2004). Press, Boston, 1990).
[62] I. Mezić, Nonlinear Dynamics 41, 309 (2005).
[63] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and
D. S. Henningson, J. Fluid Mech. 641, 115 (2009).