# Quantum Krylov subspace algorithms
for ground and excited state energy estimation

Cristian L. Cortes  Center for Nanoscale Materials,
Argonne National Laboratory, Lemont, Illinois 60439, USA  Stephen K. Gray
Center for Nanoscale Materials,
Argonne National Laboratory, Lemont, Illinois 60439, USA

###### Abstract

Quantum Krylov subspace diagonalization (QKSD) algorithms provide a low-cost
alternative to the conventional quantum phase estimation algorithm for
estimating the ground and excited-state energies of a quantum many-body
system. While QKSD algorithms typically rely on using the Hadamard test for
estimating Krylov subspace matrix elements of the form,
$\langle\phi_{i}|e^{-i\hat{H}\tau}|\phi_{j}\rangle$, the associated quantum
circuits require an ancilla qubit with controlled multi-qubit gates that can
be quite costly for near-term quantum hardware. In this work, we show that a
wide class of Hamiltonians relevant to condensed matter physics and quantum
chemistry contain symmetries that can be exploited to avoid the use of the
Hadamard test. We propose a multi-fidelity estimation protocol that can be
used to compute such quantities showing that our approach, when combined with
efficient single-fidelity estimation protocols, provides a substantial
reduction in circuit depth. In addition, we develop a unified theory of
quantum Krylov subspace algorithms and present three new quantum-classical
algorithms for the ground and excited-state energy estimation problems, where
each new algorithm provides various advantages and disadvantages in terms of
total number of calls to the quantum computer, gate depth, classical
complexity, and stability of the generalized eigenvalue problem within the
Krylov subspace.

## Introduction

The eigenpair problem for large matrices, which often consists of finding the
smallest $k$ (or largest $k$) eigenvalues and eigenvectors of a matrix,
remains one of the most ubiquitous problems in science. Within physics and
chemistry, this problem is equivalent to finding the ground and low-lying
excited state energies of a quantum many-body system represented by a large
Hamiltonian matrix. While quantum computers provide a scalable route for
solving this problem through the multi-ancilla-based quantum phase estimation
algorithm, this approach will most likely require fault-tolerant quantum
computing hardware Lloyd (1996); Abrams and Lloyd (1997); Zalka (1998);
Aspuru-Guzik  _et al._ (2005); Nielsen and Chuang (2002); Kitaev (1997); Cleve
_et al._ (1998). As a result, variational quantum algorithms such as the
variational quantum eigensolver (VQE) Peruzzo  _et al._ (2014); McClean  _et
al._ (2016); Kandala  _et al._ (2017) and quantum approximate optimization
algorithm (QAOA) Farhi  _et al._ (2014) have emerged as possible candidate
algorithms capable of dealing with the constraints of the current hardware in
the noisy intermediate scale quantum (NISQ) era Preskill (2018).

Variational quantum algorithms aim to solve an optimization problem,
$\text{min}_{\boldsymbol{\theta}}\;C(\boldsymbol{\theta})$, encoded through a
cost function that is typically written in the form
$C(\boldsymbol{\theta})=\braket{\Psi(\boldsymbol{\theta})}{H}{\Psi(\boldsymbol{\theta})}$,
where $H$ represents a Hermitian operator that encodes the problem of interest
Cerezo  _et al._ (2021a). By defining a parameterized quantum circuit,
$\ket{\Psi(\boldsymbol{\theta})}=U(\boldsymbol{\theta})\ket{0}^{\otimes N}$,
with respect to a tunable set of parameters $\boldsymbol{\theta}$, e.g.
single-qubit Pauli rotation gates, the quantum computer provides estimates of
$C(\boldsymbol{\theta})$ while the classical computer performs an optimization
subroutine that provides an update rule for the parameters
$\boldsymbol{\theta}$ (e.g. using gradient descent). This methodology can be
used to solve the eigenpair problem to estimate the ground and excited-state
energies of quantum systems Santagati  _et al._ (2018); Higgott  _et al._
(2019, 2019); Sim  _et al._ (2018). While variational quantum algorithms have
a substantial advantage in terms of gate depth, they also have significant
drawbacks. For example, it has been shown that for a large class of quantum
circuits, the optimization landscapes are highly non-convex, making the
problem of finding the global minimum NP-hard Bittel and Kliesch (2021). It
has also been shown that barren plateaus, consisting of exponentially
vanishing gradients, can also arise in a wide range of conditions McClean  _et
al._ (2018); Wang  _et al._ (2020); Cerezo  _et al._ (2021b); Arrasmith  _et
al._ (2020). For such cases, the optimization problem becomes intractable due
to the inability to update the optimization parameters $\boldsymbol{\theta}$.

In recent years, quantum subspace diagonalization (QSD) methods have emerged
as an alternative way of solving the eigenvalue problem for large matrices,
capable of dealing with the aforementioned drawbacks Huggins  _et al._ (2020);
Parrish and McMahon (2019); Seki and Yunoki (2021); Stair  _et al._ (2020).
The basic idea consists of using a set of non-orthogonal quantum states,
easily preparable on quantum computers, which can be used to define a
generalized eigenvalue problem where the size of the corresponding matrices
are exponentially smaller. The hybrid quantum-classical algorithm consists of
using the quantum computer to compute the relevant matrix elements, while the
generalized eigenvalue problem is solved on the classical computer (see Figure
1). The solution of the generalized eigenvalue problem provides an estimate of
the relevant eigenpairs.

A number of interesting QSD methods have been proposed which can be classified
according to the numerous ways that the non-orthogonal states are defined. For
instance, McClean et al. showed that by using the set of non-orthogonal basis
states $a_{i}^{\dagger}a_{j}\ket{\Psi_{G}}$, it is possible to find low-lying
excited states based on the ground state
$\ket{\Psi_{G}}\approx\ket{\Psi(\boldsymbol{\theta})}$, which is found through
the standard variational quantum eigensolver McClean  _et al._ (2017); Colless
_et al._ (2018). $K$-moment states have also been proposed as an alternative
way of constructing the non-orthogonal basis states, which becomes scalable
when the $K$-moment unitaries are tensor products of Pauli operators Bharti
and Haug (2020a, b).

It is also possible to have provable guarantees for convergence if the set of
non-orthogonal states form a Krylov basis Golub and Van Loan (1996), which is
defined by the repeated application of the matrix of interest, $H$, acting on
the initial guess vector $\ket{\phi_{o}}$, resulting in the Krylov subspace
$\mathcal{K}_{M}=\text{span}\\{\ket{\phi_{o}},H\ket{\phi_{o}},H^{2}\ket{\phi_{o}},\cdots,H^{M-1}\ket{\phi_{o}}\\}$.
The Lanczos method Parlett (1980) is one of the most well-known algorithms
that uses this subspace to solve the eigenpair problem with convergence
properties that are dependent on the spectral properties of the matrix $H$ as
well as the overlap of the guess vector $\ket{\phi_{o}}$ with the true
solution. While this method is routinely executed on classical computers, its
implementation on a quantum computer is more challenging since $H$ is not
unitary. Nonetheless, interesting approaches that invoke sums of unitary
operators as approximations to the Hamiltonian matrix and its higher powers
have recently been suggested Seki and Yunoki (2021); Bespalova and Kyriienko
(2021) and remains an ongoing research direction. Motta et al. have also
proposed the QLanczos algorithm which defines the Krylov subspace by the
repeated application of the imaginary time evolution propagator,
$f(\hat{H})=e^{-\beta\hat{H}}$ Motta  _et al._ (2020). In this framework, the
non-unitary imaginary-time propagator is written as a unitary operator under
the condition that the Hamiltonian is $k$-local. A linear system of equations
must be solved classically for each imaginary time step, where the number of
measurements and size of such equations grows exponentially with the spreading
of entanglement.

Figure 1: Overview of quantum Krylov subspace algorithms.

In this manuscript, we focus on solving the eigenpair problem with sets of
Krylov basis states generated by real-time quantum dynamics. This idea, in the
context of quantum computing approaches to the eigenpair problem, was
pioneered by Parrish and McMahon Parrish and McMahon (2019). They referred to
their approach as the quantum filter diagonalization (QFD) algorithm, because
of similarities with the classical-computer-based filter diagonalization
methods developed in the 1990s Wall and Neuhauser (1995); Mandelshtam and
Taylor (1997); Mandelshtam (2001). Independently, Stair et al. proposed a
multi-reference selected quantum Krylov subspace (MRSQK) algorithm Stair  _et
al._ (2020) which can be viewed as a generalization of QFD. These methods
represent variants of the QLanczos algorithm where the real-time evolution
operator, $e^{-i\hat{H}\tau}$, is used to generate the Krylov basis, where
$\tau$ is equal to the time step size and we assume atomic units such that
$\hbar$ = 1 throughout this manuscript. While these methods have shown great
promise, they are not without practical issues with respect to NISQ-era
applications. First, quantum Krylov subspace algorithms based on real-time
dynamics require Hadamard test quantum circuits (see Appendix for more
details), which uses an ancilla qubit with controlled multi-qubit controlled
unitary operations Parrish and McMahon (2019); Stair  _et al._ (2020). This
approach substantially increases the circuit depth, making it more difficult
for NISQ-era hardware. Second, the number of calls to the quantum computer
that are required to construct the Krylov subspace matrices scales as
$\mathcal{O}(LM^{2})$ where $L$ is the number of terms in the Hamiltonian and
$M$ is the subspace matrix size Stair  _et al._ (2020). Third, single-
reference Krylov subspace algorithms also suffer from large condition numbers
(ratio of largest to smallest singular values) of the overlap matrix
$\mathbf{S}$ that become substantially worse as the number of time steps
increases. In principle, this could make the solution of the generalized
eigenvalue problem not possible for many problems of interest, such as
strongly correlated systems.

In this work, we provide several major contributions which address the
outstanding problems discussed above. First, we show that the Hadamard test is
not required to estimate the Krylov subspace elements of the form
$\braket{\phi_{i}}{e^{-i\tau\hat{H}}}{\phi_{j}}$ for a large class of
Hamiltonians relevant to nuclear physics, quantum chemistry, and condensed
matter physics. Our approach avoids the need for an ancilla qubit with
controlled unitary operations and, when combined with efficient fidelity
estimation protocols, provides a substantial reduction in circuit depth
compared to previous approaches. Our second major contribution includes the
proposal of three new generalized eigenvalue problems which can be used to
estimate both ground and excited-state energies. Each of these generalized
eigenvalue problems provides various advantages and disadvantages in terms of
the total number of calls to the quantum computer, gate depth, classical post-
processing complexity, as well as stability of the generalized eigenvalue
problem based on the condition number of the overlap matrix $\mathbf{S}$.

In particular, we show that two of the newly proposed generalized eigenvalue
problems involving the unitary function $e^{-i\hat{H}\tau}$ only require
$\mathcal{O}(M)$ calls to the quantum computer compared to the
$\mathcal{O}(LM^{2})$ calls that are required for QKSD algorithms which use
Hamiltonian-based generalized eigenvalue problems. We also show that two of
the newly proposed generalized eigenvalue problems, which also use real-time
quantum dynamics to generate the Krylov space as in the QFD approach Parrish
and McMahon (2019) more closely resemble the original classical filter
diagonalization method (FDM) originally proposed by Wall and Neuhauser Wall
and Neuhauser (1995) and subsequently further developed and elaborated upon by
Mandelshtam and Taylor Mandelshtam and Taylor (1997) in that they target
specific energy ranges by energy filtering. Compared to typical Krylov
subspace algorithms which have a classical computational complexity scaling
that is polynomial in the total number of time steps, the FDM method provides
a constant time scaling $\mathcal{O}(1)$. We also find empirical evidence that
the corresponding generalized eigenvalue problems have condition numbers that
are orders of magnitude smaller than the single-reference Krylov subspace
counterparts as the number of time steps increases. To test the efficacy of
the proposed algorithms, we numerically compare these four methods for the
problem of finding the ground-state and excited-state energies of various
quantum chemistry Hamiltonians, showing fast convergence with a small number
of discrete time steps.

It is worth noting that while we were finalizing this manuscript, we became
aware of the paper by Klymko _et al._ Klymko  _et al._ (2021), currently in
preprint, which is similar in spirit to our work. Klymko _et al._ provide a
theoretical basis for non-orthogonal states generated by real-time dynamics,
and independently proposed a hybrid quantum-classical algorithm based on the
generalized eigenvalue problem with the unitary function, $e^{-i\hat{H}\tau}$,
which they refer to as variational QPE or VQPE. This method is equivalent to
one of the three new methods that we present in this manuscript, which we
refer to as KDM U. In addition, their major contributions include a
comprehensive study of the effects of noise, Trotter-Suzuki error, and a
comparison between VQPE and the conventional QPE algorithm in terms of the
total simulation time and total number of time steps required to reach
chemical accuracy. Our major contributions include the multi-fidelity
estimation protocol which may be used to avoid the Hadamard test, as well as
the two hybrid quantum-classical algorithms based on Fourier filter energies,
which more closely resemble the classical filter diagonalization method of
Wall and Neuhauser, and Mandelshtam and Taylor Wall and Neuhauser (1995);
Mandelshtam and Taylor (1997), which we refer to as the FDM H and FDM U
methods.

##  I Quantum Krylov subspace diagonalization method

We aim to find the ground and excited-state energies of a general many-body
Hamiltonian written as a sum of $N$-qubit Pauli terms,
$\hat{H}=\sum_{i}^{L}h_{i}\hat{P}_{i}$, where $h_{i}$ is a weighting
coefficient and $\hat{P}_{i}$ is a general tensor product of $N$ Pauli
operators,
${\hat{P}}_{i}=\otimes_{k=1}^{N_{i}}\hat{\sigma}_{i_{k}}^{(\mu_{k})}$, with
$\mu_{k}$ denoting the qubit number and $i_{k}$ acts as a label for the type
of Pauli operator
$\\{\hat{I},\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}\\}$. We do not
impose any type of restrictions on the locality of the Hamiltonian, thereby
allowing for the implementation of the proposed algorithms for many problems
relevant to nuclear physics, condensed matter physics, and quantum chemistry.
For illustration purposes, we will focus on the canonical quantum chemistry
Hamiltonian in second quantization which is able represent the electronic
structure problem of a wide variety of molecular systems (see appendix for
details).

The Hamiltonian $\hat{H}$ obeys the standard eigenvalue equation,
$\hat{H}\ket{\psi_{k}}=E_{k}\ket{\psi_{k}}$, with the energy eigenvalue
$E_{k}$ and corresponding eigenvector $\ket{\psi_{k}}$, assumed to satisfy the
orthonormality condition,
$\braket{\psi_{k^{\prime}}}{\psi_{k}}=\delta_{kk^{\prime}}$. General functions
of the Hamiltonian $f(\hat{H})$ will also obey the eigenvalue equation,

$f(\hat{H})\ket{\psi_{k}}=f(E_{k})\ket{\psi_{k}}.$ (1)

The matrix size for this eigenvalue problem scales exponentially with the
total number of qubits. However, as we show below, this equation can be used
to define a wide variety of generalized eigenvalue problems in a subspace that
is exponentially smaller. The basic idea requires expanding the eigenvector
$\ket{\psi_{k}}$ as a linear combination of non-orthogonal states
$\ket{\phi_{n}}$,

$\ket{\psi_{k}}\approx\sum_{n=0}^{M-1}c_{n}\ket{\phi_{n}}.$ (2)

Substituting this result into (1) and multiplying from the left by
$\bra{\phi_{n^{\prime}}}$, we find the generalized eigenvalue problem,

$\mathbf{F}(\hat{H})\mathbf{c}=f(E_{k})\mathbf{S}\mathbf{c},$ (3)

where $\mathbf{c}=(c_{0},c_{1},\cdots,c_{M-1})^{T}$ is a column vector of
expansion coefficients and the subspace matrices $\mathbf{F}(\hat{H})$ and
$\mathbf{S}$ are defined by the matrix elements,

$\displaystyle[\mathbf{F}(\hat{H})]_{nn^{\prime}}=\braket{\phi_{n}}{f(\hat{H})}{\phi_{n^{\prime}}}\;\text{and}\;[\mathbf{S}]_{nn^{\prime}}=\braket{\phi_{n}}{\phi_{n^{\prime}}}.$
(4)

Naturally, the subspace matrix size is much smaller in the non-orthogonal
basis when $M\ll 2^{N}$, and is also more general than Hamiltonian-based
generalized eigenvalue problems derived in previous works. The choice of
Hamiltonian function $f(\hat{H})$ and non-orthogonal basis, $\ket{\phi_{n}}$,
will ultimately lead to a wide variety of different hybrid quantum-classical
algorithms with trade-offs in terms of convergence, number of calls to the
quantum computer, circuit depth, and classical complexity required for post-
processing. Although it might not be implementable in the near-term, it is
worth mentioning that a non-orthogonal basis defined by the function,
$f(\hat{H})=\exp\left(i\arccos{(\hat{H}\tau)}\right)$, as used in qubitization
Low and Chuang (2019), would provide an effective way to estimate both ground
and excited-state energies with equivalent circuit depths as a single Trotter
time step but avoiding Trotter error Poulin  _et al._ (2018). However, this
comes at a cost adding a register of ancilla qubits with controlled multi-
qubit unitary gates. For our purposes, we will focus on the standard
Hamiltonian and real-time evolution operators, $f(\hat{H})=\hat{H}$ and
$f(\hat{H})=e^{-i\hat{H}\tau}$. We will also consider two different sets of
non-orthogonal states, thereby obtaining four different quantum-classical
algorithms which provide various advantages and disadvantages as discussed
below.

##### Krylov subspace diagonalization method (KDM).

The Krylov subspace diagonalization method assumes the eigenvector
$\ket{\psi_{k}}$ may be written as a linear combination of real-time evolved
Krylov basis states,

$\ket{\psi_{k}}\approx\ket{\psi_{K}}=\sum_{n=0}^{M-1}c_{n}e^{-in\hat{H}\tau}\ket{\phi_{o}}=\sum_{n=0}^{M-1}c_{n}\ket{\phi_{n}}$
(5)

where $\ket{\phi_{o}}$ is the initial single-reference state. Using the steps
outlined above, the corresponding generalized eigenvalue problem may be
written as,
$\mathbf{F}_{K}(\hat{H})\mathbf{c}_{K}=f(E_{k})\mathbf{S}_{K}\mathbf{c}_{K}$,
where the subscript $K$ denotes the real-time Krylov basis with the subspace
matrix elements defined using equation (4). We emphasize that this basis,
specifically when $f(\hat{H})=\hat{H}$, which we refer to as the KDM H method,
recovers the QFD and MRSQK methods in the limit of a single reference state
Parrish and McMahon (2019); Stair  _et al._ (2020). Our new contribution
therefore corresponds to the KDM method with $f(\hat{H})=e^{-i\hat{H}\tau}$,
which we refer to as the KDM U method, which also coincides with the recent
proposal by Klymko et al. Klymko  _et al._ (2021).

##### Filter Diagonalization Method (FDM).

The filter diagonalization method, on the other hand, approximates
$\ket{\psi_{k}}$ as a linear combination of time-evolved wavefunctions that
are Fourier transformed with respect to a set of filter energies $E_{j}$,

$\displaystyle\ket{\psi_{k}}\approx\ket{\psi_{F}}=\sum_{j=1}^{J}\sum_{n=0}^{M-1}c_{j}e^{-in(\hat{H}-E_{j})\tau}\ket{\phi_{o}}=\sum_{j}^{J}c_{j}\ket{\phi_{j}}$
(6)

resulting in the generalized eigenvalue problem,
$\mathbf{F}_{J}(\hat{H})\mathbf{c}_{J}=f(E_{k})\mathbf{S}_{J}\mathbf{c}_{J}$
with matrix elements defined by the non-orthogonal basis of filter energies.
This type of basis has an interesting property in that $\ket{\phi_{j}}$ will
be dominated by eigenvectors whose eigenvalues are close to the filter
energies $E_{j}$. If we expand the starting state, $\ket{\phi_{o}}$, in terms
of the true eigenvectors $\ket{\psi_{k^{\prime}}}$ of the Hamiltonian, such
that $\ket{\phi_{o}}=\sum_{k^{\prime}}c_{k^{\prime}}\ket{\psi_{k^{\prime}}}$,
then it is possible to show that,

$\ket{\phi_{j}}=\sum_{k^{\prime}}c_{k^{\prime}}\frac{1-e^{-iM(E_{k^{\prime}}-E_{j})\tau}}{1-e^{-i(E_{k^{\prime}}-E_{j})\tau}}\ket{\psi_{k^{\prime}}},$
(7)

highlighting the fact that eigenvectors with eigenvalues $E_{k^{\prime}}$ that
are close to the filter energies $E_{j}$ ($E_{k^{\prime}}\sim E_{j}$) will
have the largest contribution. In total, the FDM method will give rise to two
new quantum-classical algorithms which we refer to as the FDM H method if
$f(\hat{H})=\hat{H}$ and FDM U method if $f(\hat{H})=e^{-i\hat{H}\tau}$.

##### Relationship between both non-orthogonal bases.

Analyzing the two generalized eigenvalue problems, it is possible to show that
the two approaches are related by the $M\times J$ transformation matrix,

$\mathbf{W}=\begin{pmatrix}1&1&\cdots&1\\\
e^{-iE_{1}\tau}&e^{-iE_{2}\tau}&\cdots&e^{-iE_{J}\tau}\\\
e^{-i2E_{1}\tau}&e^{-i2E_{2}\tau}&\cdots&e^{-i2E_{J}\tau}\\\
\vdots&\vdots&\ddots&\vdots\\\
e^{-i(M-1)E_{1}\tau}&e^{-i(M-1)E_{2}\tau}&\cdots&e^{-i(M-1)E_{J}\tau}\\\
\end{pmatrix}$

resulting in the following relationship between the real-time evolution Krylov
diagonalization method and the filter diagonalization method,

$\displaystyle\mathbf{W}^{\dagger}\mathbf{F}_{K}(\hat{H})\mathbf{W}$
$\displaystyle=\mathbf{F}_{J}(\hat{H})$ (8)
$\displaystyle\mathbf{W}^{\dagger}\mathbf{S}_{K}\mathbf{W}$
$\displaystyle=\mathbf{S}_{J}.$ (9)

In the case that the filter frequencies are chosen with an equidistant grid,
such that $E_{j}=\frac{2\pi}{M\tau}j$ where $j=0,\cdots,M-1$, then the
transformation $\mathbf{W}$ becomes a unitary matrix up to a normalization
factor equivalent to the discrete Fourier transform matrix. It is important to
emphasize, however, that the total number of discrete energies $J$ can be much
smaller than the total number of time steps, $M$, resulting in a constant-time
$\mathcal{O}(1)$ computational complexity for solving the FDM-based
generalized eigenvalue problem on the classical computer, compared to the
polynomial scaling $\mathcal{O}(\text{poly}(M))$ for the KDM method. Moreover,
the choice of filter window with a suitable number of filter energies can also
stabilize the generalized eigenvalue problem resulting in smaller condition
numbers for the overlap matrix $\mathbf{S}$. We provide evidence of this
result in the numerical experiments section of the manuscript.

$f(\hat{H})$ | Non-orthogonal Basis | $N_{K}$ | Classical Complexity
---|---|---|---
$\hat{H}$ | real-time dynamics | $\mathcal{O}(LM^{p}/\epsilon^{2})$ | $\mathcal{O}(\text{poly}(M))$
$e^{-i\hat{H}\tau}$ | real-time dynamics | $\mathcal{O}(M/\epsilon^{2})$ | $\mathcal{O}(\text{poly}(M))$
$\hat{H}$ | Fourier filter energies | $\mathcal{O}(LM^{p}/\epsilon^{2})$ | $\mathcal{O}(1)$
$e^{-i\hat{H}\tau}$ | Fourier filter energies | $\mathcal{O}(M/\epsilon^{2})$ | $\mathcal{O}(1)$
Table 1: Summary of results. $N_{K}$ is equal to the number of
measurements/calls to the quantum computer required to construct the Krylov
subspace matrix elements. $L$ is equal to the total number of terms in the
Hamiltonian $\hat{H}$. $M$ is equal to the order of real-time evolved Krylov
subspace. The exponent $p$ is equal to one if the Hamiltonian, $\hat{H}$,
perfectly commutes with the quantum circuit unitary which approximates the
time-evolution operator, $e^{-i\tau\hat{H}}$, otherwise it is equal to two.

##### Number of calls to the quantum computer.

In the following, we provide an estimate of the number of calls to the quantum
computer, $N_{K}$, required to construct the Krylov subspace matrix elements
defined by Equation (4). We assume that the subspace matrix elements are
computed with Hadamard-test quantum circuits, or equivalently, with the multi-
fidelity estimation protocol which we describe in the following section. For
the latter, our estimates are based on single fidelity estimation circuits,
such as a SWAP-test circuit or a mirror-like quantum circuit, as outlined in
the Discussion section below. In general, the estimation of these quantities
to precision $\epsilon$ will require $1/\epsilon^{2}$ samples, which will
result in a $1/\epsilon^{2}$ multiplicative factor for all of the cases we
consider below. Furthermore, we will restrict ourselves to the single-
reference Krylov subspace algorithm, though more general estimates of the
multi-reference case may be done with the same arguments.

We first consider the estimation of the overlap matrix elements in
$\mathbf{S}$, noting that they will be the same for all four methods. These
matrix elements are equivalent to correlation functions of the form
$C_{n}(\tau)=\braket{\phi_{o}}{e^{-in\hat{H}\tau}}{\phi_{o}}$. By assuming
that a single call to the quantum computer provides an estimate of both the
real and imaginary parts of the correlation function $C_{n}(\tau)$, then $M-1$
calls to the quantum computer are required to construct the overlap matrix
$\mathbf{S}$.

For $f(\hat{H})=\hat{H}$, the matrix elements of the subspace matrix
$\mathbf{F}(\hat{H})$ may be written as
$\braket{\phi_{n}}{\hat{H}}{\phi_{n^{\prime}}}=\braket{\phi_{o}}{e^{in\hat{H}\tau}\hat{H}e^{-in^{\prime}\hat{H}\tau}}{\phi_{o}}$.
Here, the number of calls to the quantum computer will depend on whether the
Hamiltonian $\hat{H}$ and the quantum circuit implementation of the time-
evolution operator, $e^{-in\hat{H}\tau}$, commute. If we assume that they
commute, these elements may be written as
$\braket{\phi_{o}}{\hat{H}e^{-i(n^{\prime}-n)\hat{H}\tau}}{\phi_{o}}=\sum_{i}^{L}h_{i}\braket{\phi_{o}}{\hat{P}_{i}e^{-i(n^{\prime}-n)\hat{H}\tau}}{\phi_{o}}$,
resulting in a Toeplitz matrix structure that requires $\mathcal{O}(LM)$ calls
to the quantum computer. However, in the case of Trotterized quantum circuits
where the commutation relation does not hold exactly, the total number of
calls would scale as $\mathcal{O}(LM^{2})$.

Methods using the real-time evolution function,
$f(\hat{H})=e^{-i\hat{H}\tau}$, will have a complexity that is substantially
less. In this case, the matrix elements will also correspond to correlation
functions, $C_{n}(\tau)$. The matrix elements from the overlap matrix,
$\mathbf{S}$, can therefore be used to construct the majority of the matrix
elements in $\mathbf{F}(\hat{H})$. The off-diagonal elements in the top-right
and bottom-left corners, however, will require an additional call to the
quantum computer for the estimation of the $C_{M}(\tau)$ correlation function.
In total, the $f(\hat{H})=e^{-i\hat{H}\tau}$ method will require $M$ calls to
the quantum computer. This, however, comes at the cost of requiring a single
quantum circuit evaluation with an increased circuit depth equivalent to a
single time step (assuming a Trotterized time-evolution circuit).

Finally, it is worth noting that for fixed $f(\hat{H})$, the KDM and FDM
methods will have an equivalent number of calls to the quantum computer
because they are related by equations (8) and (9). This highlights the fact
that both methods only differ in the post-processing methodology used to
estimate the ground and excited-state energies and, as a result, both methods
can be carried out in parallel on a classical computer. A summary of these
results is shown in Table I, underlining how each of these four methods
carries different complexities due to quantum and classical computational
resources. Here, we wrote the scaling of the number of calls $N_{K}$ for
$f(\hat{H})=\hat{H}$ as $\mathcal{O}(LM^{p})$ where $p=1$ if the Hamiltonian
and time-evolution unitary circuit commutes and $p=2$ otherwise.

##  II Multi-fidelity estimation protocol

We now consider the evaluation of the complex-valued matrix elements (4),
equivalent to a single call the quantum computer as defined above. The
Hadamard test is the standard method used for estimating these types of non-
Hermitian quantities (see Appendix D for more details), which originates from
the single-ancilla-based quantum phase estimation algorithm from Kitaev. This
approach requires an ancilla qubit with controlled unitary operations that
substantially increases the total number of multi-qubit gates in the overall
circuit. In the near term, multi-qubit gates (e.g. CNOT gates) represent an
expensive resource. In the following, we propose a method that avoids the
Hadamard test thereby making a wide variety of quantum Krylov subspace
diagonalization methods more amenable to near-term devices.

The motivation for our proposed method stems from the field of interferometry
which aims to measure a target phase $\theta_{t}$ that encodes a physical
parameter of interest. Interference pattern measurements can only provide
information about the phase difference, $\Delta\theta=\theta_{r}-\theta_{t}$.
A reference laser is typically used to provide a controllable reference phase
$\theta_{r}$, allowing for the proper estimation of $\theta_{t}$. While the
Hadamard test provides a reference phase through use of the ancilla qubit, the
multi-fidelity estimation (MFE) protocol generates the reference phase through
the superposition state $\tfrac{1}{\sqrt{2}}(\ket{R}+\ket{\phi_{k}})$, where
the reference state $\ket{R}$ and the target state $\ket{\phi_{k}}$ originate
from different symmetry sectors of the Hamiltonian. If the time evolution of
the reference state is classically simulatable, it will be possible to have a
reference phase without an ancilla qubit. Below, we provide a more detailed
mathematical description.

The proposed approach assumes that the Hamiltonian contains a symmetry
$\hat{S}$ such that, $[\hat{H},\hat{S}]=0$, with quantum states
$\ket{\phi_{k}}$ that have a definite symmetry, such that
$\braket{\phi_{k}}{\hat{S}}{\phi_{k}}=s_{k}$, where $s_{k}$ corresponds to an
eigenvalue of the symmetry operator $\hat{S}$. We also assume that there
exists a reference state $\ket{R}$ where $\braket{R}{e^{-in\tau\hat{H}}}{R}$
is efficient to calculate on the classical computer. We emphasize that the
reference and target states, $\ket{R}$ and $\ket{\phi_{k}}$, are not required
to be eigenstates of the Hamiltonian, but they do need to originate from
different symmetry sectors such that $\braket{R}{\phi_{k}}=0$. Note that this
condition relaxes the requirements from conventional quantum phase estimation
algorithms that normally require the initial state to be an eigenstate of the
Hamiltonian. If all of these conditions hold, then it will be possible to
implement the multi-fidelity estimation protocol as shown below. For many
nuclear physics, quantum chemistry and condensed matter physics applications,
the particle number, total spin and spin projection symmetries may be
applicable and can be used in this approach due to the fact that they contain
symmetry sectors that are classically tractable.

As a concrete example, we consider the quantum chemistry Hamiltonian which
conserves the electron number
$\hat{N}=\sum_{i}\hat{a}_{i}^{\dagger}\hat{a}_{i}$. We assume that the quantum
states $\ket{\phi_{k}}$ have a definite electron number,
$\braket{\phi_{k}}{\hat{N}}{\phi_{k}}=n_{k}$, that is not equal to zero. To
estimate the most general off-diagonal element,
$\braket{\phi_{i}}{e^{-in\hat{H}\tau}}{\phi_{j}}$, as required by multi-
reference Krylov-based methods, the MFE protocol requires measuring the
following state fidelities on the quantum computer,

$\displaystyle F_{1}$
$\displaystyle=|\braket{\phi_{i}}{e^{-in\hat{H}\tau}}{\phi_{j}}|^{2},$ (10)
$\displaystyle F_{2}$
$\displaystyle=\tfrac{1}{4}|(\bra{\phi_{i}}+\bra{R})|e^{-in\hat{H}\tau}|(\ket{R}+\ket{\phi_{j}})|^{2}.$
(11)

Combining both results yields the magnitude and phase of the off-diagonal
matrix element $\braket{\phi_{i}}{e^{-in\hat{H}\tau}}{\phi_{j}}=re^{i\theta}$,
written in complex polar coordinates,

$\displaystyle r$ $\displaystyle=\sqrt{F_{1}}$ (12) $\displaystyle\theta$
$\displaystyle=\cos^{-1}\left(\frac{4F_{2}-F_{1}-r_{R}^{2}}{2r_{R}\sqrt{F_{1}}}\right)+\theta_{R}$
(13)

where $r_{R}$ and $\theta_{R}$ represent the reference amplitude and phase
defined as, $\braket{R}{e^{-in\hat{H}\tau}}{R}=r_{R}e^{i\theta_{R}}$. If we
take the reference state as the zero particle number (vacuum) state,
$\ket{R}\equiv\ket{0}^{\otimes N}$, then the reference amplitude $r_{R}$ will
be equal to one while the reference phase will be equal to
$\theta_{R}=-n\tau\braket{0}{\hat{H}}{0}$, where $\braket{0}{\hat{H}}{0}$
denotes the expectation value of the Hamiltonian with respect to the vacuum
state which can be evaluated efficiently on a classical computer.

###  II.1 Preparation of $\tfrac{1}{\sqrt{2}}(\ket{R}+\ket{\phi_{k}})$

One of the key requirements of this protocol is the preparation of the state
$\tfrac{1}{\sqrt{2}}(\ket{R}+\ket{\phi_{k}})$ on the quantum computer. Note
that because we have imposed the requirement that $\ket{R}$ and
$\ket{\phi_{k}}$ belong to different symmetry sectors (i.e. contain different
particle numbers), it is possible to prepare such states using GHZ-state-
preparation circuits. For instance, we consider the preparation of the state
$\tfrac{1}{\sqrt{2}}(\ket{R}+\ket{\phi_{\text{HF}}})$ where $\ket{R}$ is the
vacuum state and $\ket{\phi_{\text{HF}}}$ is the Hartree-Fock state for a
system of $N$ spin-orbitals (represented by $N$ qubits) and $\eta$ electrons.
In the Jordan-Wigner basis, the Hartree-Fock state takes the simple product-
state form, $\ket{\phi_{\text{HF}}}=\ket{0}^{\otimes
N-\eta}\otimes\ket{1}^{\otimes\eta}$, where the first $\eta$ qubits are
prepared in the one state and the rest of the qubits remain in the zero state.
To prepare the target superposition state,
$\tfrac{1}{\sqrt{2}}(\ket{0}^{\otimes N}+\ket{\phi_{\text{HF}}})$, a Hadamard
gate is applied to the first qubit, followed by a ladder of CNOT gates applied
up to the $\eta$th qubit, resulting in a total of $\eta-1$ CNOT gate
operations. More general states can also be prepared by subsequently applying
a symmetry-conserving quantum circuit $U_{S}$ Gard  _et al._ (2020), which
could in principle represent a parameterized quantum circuit originating from
a VQE pre-processing step. As an example, we provide the quantum circuit that
prepares the state, $\tfrac{1}{\sqrt{2}}(\ket{000000}+\ket{000111})$, where
the Hartree-Fock state represents a system of six spin-orbitals with three
electrons.

$\Qcircuit@C=2em@R=1em@!R{\lstick{\ket{0}}&\gate{H}\ctrl{1}\qw\multigate{5}{U_{S}}\qw\\\
\lstick{\ket{0}}\qw\targ\ctrl{1}\ghost{U_{S}}\qw\\\
\lstick{\ket{0}}\qw\qw\targ\ghost{U_{S}}\qw\\\
\lstick{\ket{0}}\qw\qw\qw\ghost{U_{S}}\qw\\\
\lstick{\ket{0}}\qw\qw\qw\ghost{U_{S}}\qw\\\
\lstick{\ket{0}}\qw\qw\qw\ghost{U_{S}}\qw\\\
\mbox{$\tfrac{1}{\sqrt{2}}(\ket{00000}+\ket{000111})$}\gategroup{1}{2}{6}{4}{1.5em}{--}}$

Figure 2: Comparison of KDM and FDM methods for predicting the ground-state
potential energy curves of a (a) 6-site linear hydrogen chain, (b) H2O
molecule, and (c) BeH2 molecule. The top row plots the absolute energy scale
measured in Hartrees. The second and third rows plot the energy error $\Delta
E=|E_{\text{FCI}}-E_{\text{approx}}|$ between the full configuration
interaction calculation and the proposed approximate methods, where the second
row uses a narrow frequency window,
$[E_{\text{HF}}-0.3E_{h},E_{\text{HF}}+0.2E_{h}]$ for the FDM method, while
the third row uses a large frequency window,
$[E_{\text{HF}}-20E_{h},E_{\text{HF}}+20E_{h}]$.

## Numerical Experiments

In Figure 2, we compare these four methods for the estimation of the ground-
state energy curves for three different molecular systems: (a) a one-
dimensional $H_{6}$ chain, (b) an H2O molecule with fixed bond angle
$\phi=104.45^{\circ}$, and (c) a BeH2 molecule. The energy curves are plotted
as a function of interatomic distance $R$ labelled in the insets. See Appendix
A for a discussion of the various levels of electronic structure theory used
in these examples. In all cases, we numerically simulate the hybrid quantum-
classical algorithm as described in Figure 1 using the Hartree-Fock (HF) state
as the initial single-reference state $\ket{\phi_{o}}$ and an ideal time-
evolution quantum circuit $U(n\tau)=e^{-in\hat{H}\tau}$ where the time step
size, $\tau=0.1$ a.u., is used for all of the simulations. The choice of time
step is very important but can be a subtle issue that is discussed in Appendix
B. In the top row, we plot the absolute energy scale measured in Hartrees
while the bottom two rows show the energy error $\Delta
E=|E_{\text{FCI}}-E_{\text{approx}}|$. We also compare the FDM method with a
small frequency window $[E_{\text{HF}}-0.3E_{h},E_{\text{HF}}+0.2E_{h}]$ in
the second row and a large frequency window
$[E_{\text{HF}}-20E_{h},E_{\text{HF}}+20E_{h}]$ in the third row. The
conventional implementation of the FDM method uses a narrow frequency window
to select the filter energies $E_{j}$ Mandelshtam and Taylor (1997),
corresponding to the methodology used in the second row. Using this approach,
we find that the FDM method converges much more slowly than the equivalent KDM
method, thereby requiring more time steps to achieve similar accuracy. Instead
we found that by using large frequency windows to choose the filter energies
$E_{j}$, it is possible to improve the convergence rate of the FDM method so
that it is comparable to the KDM method. This is shown by the nearly
overlapping energy error curves in the third row. In the latter case, we find
that we are able to achieve chemical accuracy within 3 to 6 time steps for all
three molecular systems. While we do not consider Trotter error in this work,
it is worth noting that previous work Parrish and McMahon (2019); Stair  _et
al._ (2020) has shown that Trotterized quantum circuits provide an additional
error, $\Delta_{t}$, that increases the energy error $\Delta E$ of the ground-
state energy estimate. We expect similar behavior for the quantum-classical
algorithms that we have proposed in this manuscript, but will leave such
studies for future work.

Figure 2 demonstrates the dependence of the energy error of the FDM method as
a function of the hyperparameters: $E_{\text{min}}$, $E_{\text{max}}$, and $J$
which define the energy window $[E_{\text{min}},E_{\text{max}}]$ as well as
the total number of filter energies, $J$. The choice of energy window and
number of filter energies affects the convergence rate of the FDM method quite
dramatically, which is also illustrated clearly in Figure 3. In practice, a
hyperparameter optimization loop can be added to the conventional hybrid
quantum-classical algorithm depicted in Figure 1, which would be performed on
the classical computer exclusively without any additional calls to the quantum
computer. For instance, this loop may consist of an exhaustive grid search of
the parameters $E_{\text{min}}$, $E_{\text{max}}$, and $J$ while monitoring
either the ground-state energy estimate or the variance of the function
$f(\hat{H})$ as a function of the FDM wavefunction,
$\ket{\psi_{F}}=\sum_{j}^{J}c_{j}\ket{\phi_{j}}$. For the ground-state energy
estimation problem, $c_{j}$ would correspond to the eigenvector coefficients
from the generalized eigenvalue equation (3) with the smallest eigenvalue.
Explicitly, the hyperparameter optimization loop for the FDM U method would
monitor the variance of the function, $f(\hat{H})=e^{-i\hat{H}\tau}$, written
as

$\text{Var}[e^{-i\hat{H}\tau}]=1-|\\!\\!\braket{\psi_{F}}{e^{-i\hat{H}\tau}}{\psi_{F}}\\!\\!|^{2},$
(14)

where we have assumed that the function $\ket{\psi_{F}}$ has been normalized,
such that $\braket{\psi_{F}}{\psi_{F}}=1$. This definition of the variance has
a maximum value of 1 and minimum value of zero when $\ket{\psi_{F}}$ is an
eigenstate of the Hamiltonian, $\hat{H}$. It is important to note that the
second term on the right hand side of equation (14) may be written as a sum of
the matrix elements defined by equation (4),
$\braket{\phi_{n}}{e^{-i\hat{H}\tau}}{\phi_{n^{\prime}}}$, highlighting the
fact that the variance is easily calculated on the classical computer with no
additional calls to the quantum computer.

Figure 3: (a) Energy error $\Delta E$, (b) condition number of the overlap
matrix $\mathbf{S}$, and (c) variance of the approximate wavefunction as a
function of time steps, $n$. The results are shown for the water molecule for
the extended bond length, $R=2.8$ a.u., from Figure 1-b. The operator,
$f(\hat{H})=e^{-i\hat{H}\tau}$, was used for the KDM method (blue) and FDM
method (red). The FDM results are decomposed in terms of various total number
of filter energies ($J=3,4,5$) for the small frequency window and $J=6,7,8$
for the large frequency window. The transparency of the KDM curve is
determined by a threshold value, $\kappa_{th}=1\times 10^{17}$, for
illustration purposes, where the dark blue lines correspond to
$\kappa<\kappa_{th}$.

To better compare the KDM and FDM-based algorithms, we plot the energy error,
condition number, and variance in Figure 3 as a function of total number of
time steps. We only plot the results for the water molecule case from Figure
1-b at the long interatomic distance of $R=2.8$ a.u. because it represents the
toughest point with the largest energy error for both methods. In general, we
found similar behavior for the other two molecules for different interatomic
distances. In the Figure, we only plot the $f(\hat{H})=e^{-i\hat{H}\tau}$
operator case, though similar behavior is seen for the Hamiltonian case,
$f(\hat{H})=\hat{H}$. In the top row, we present the result for the narrow
filter frequency window, $[E_{\text{HF}}-0.3E_{h},E_{\text{HF}}+0.2E_{h}]$,
while the bottom row presents the results for the large frequency window,
$[E_{\text{HF}}-20E_{h},E_{\text{HF}}+20E_{h}]$. We emphasize that we have not
performed the hyperparameter optimization loop as discussed above, but chose
the appropriate frequency window and total number of filter energies $J$ for
illustration purposes. We first analyze the KDM U method results which is the
same for the top and bottom rows. We find that the KDM U method reaches
chemical accuracy within four time steps while also reaching the minimum
variance of $2\times 10^{-6}$ within eight time steps. The condition number of
the KDM method, however, displays a troubling increasing trend as a function
of time steps. We highlight this point by controlling the transparency of the
KDM curve based on the threshold value, $\kappa_{th}=1\times 10^{17}$, where
the dark blue lines correspond to $\kappa<\kappa_{th}$. In general, we find
that the KDM method becomes unstable as the number of time steps increases, as
shown by the fluctuating energy error (light blue lines) and relatively large
variance values $\sim 10^{-2}$. Comparatively, the FDM method displays
drastically different behavior based on the chosen frequency window. While the
small frequency window stabilizes the FDM method quite well, as shown by the
decreasing condition numbers in the top row of Figure 3-b, the corresponding
energy error (Figure 3-a) and variance (Figure 3-c) decreases much more
slowly, never quite reaching the minimum energy error and variance values of
the KDM method. It is also worth noting the convergence behavior for the
different numbers of filter energies ($J=3,4,5$). The overall condition number
for the $J=3$ FDM method is much smaller than the $J=5$ case, though for a
large number of time steps, the energy error for the latter is over two orders
of magnitude smaller. On the other hand, the FDM method, based on the large
frequency window, is able to achieve fast convergence comparable to the KDM
method. As the total number of time steps increases, the FDM method displays a
much smaller energy error, condition number, and variance compared to the KDM
method. These results highlight the fact that for certain problems of
interest, where a large number of time steps are required to estimate the
ground (or excited-state) energies, the FDM method might be preferable when
compared to the KDM method, resulting in much smaller energy errors even
within 10 to 20 time steps.

Figure 4: Comparison of KDM and FDM methods for the first four singlet excited
state energies of H2O as a function of interatomic distance $R$ between the
oxygen and hydrogen atoms. (a) Absolute energy scale (measured in Hartrees) as
a function of interatomic distance. (b) Energy error of the $n$th excited
state, $\Delta E_{n}=|E^{(n)}_{\text{FCI}}-E^{(n)}_{\text{approx}}|$, as a
function of interatomic distance.

In Figure 4, we numerically simulate the estimation of the first four singlet
excited-state energies for the water molecule as a function of interatomic
distance $R$. It is important to emphasize that the excited-state energy
estimation problem has a wide variety of important applications across physics
and chemistry, where classical algorithms are often much less developed than
the corresponding ground-state energy estimation algorithms. For illustration
purposes, we only consider the singlet ($S=0$) excited-state energy levels,
and use singlet states with zero total angular momentum as the initial
starting states (see Appendix E for more details) where a time step size of
$\tau=0.1$ a.u. is used for all of the simulations. We emphasize that our
approach does not require knowledge of the ground-state wavefunction (or any
other eigenfunction), making this approach distinct from conventional excited-
state energy estimation algorithms, specifically within the scope of
variational quantum algorithms, which typically use the ground-state
wavefunction as a prerequisite. Here, we present the results for the FDM
method with a large frequency window. We found that a small frequency window
did not give very good results for the FDM method. We find that both the KDM
(blue) and FDM (red) methods perform comparably well for the first three
excited states. For the first two excited-state energy levels, both methods
achieve chemical accuracy within 4 time steps across all interatomic
distances. The third and fourth excited state energy curves are quite
challenging due to an avoided crossing at approximately $R=2.6$ a.u. which
results in nearly degenerate energy curves for the last four data points. We
find that both the KDM and FDM methods achieve chemical accuracy for the third
excited-state energy curve for nearly all of the interatomic distances, all
within three time steps. On the other hand, the fourth excited state energy
curve estimation was particularly challenging for both methods. Although it is
not shown here, we found that the KDM method only remained stable up to the
fifth time step. However, even within five time steps it was unable to reach
chemical accuracy for any of the interatomic distances. For this reason, we
present the excited-state energy predictions for the tenth time step, where
the KDM U method reaches chemical accuracy for nine out of the twenty data
points. Furthermore, as evidenced by Figure 4-a and Figure 4-b, the KDM H
method, denoted by the blue crosses, displays a large amount instability where
it is only able to reach chemical accuracy for four out of the twenty data
points. For the fourth excited-state energy curve, the FDM U method
outperforms all of the other methods, achieving chemical accuracy for
seventeen out of the twenty data points within 10 time steps. The FDM H method
is only able to reach chemical accuracy for nine out of the twenty data points
within ten time steps.

Finally, it is worth emphasizing that increasing the total number of time
steps does not help improve the convergence of the KDM method due the large
condition numbers, however, using step sizes that are much larger (e.g.
$\tau=0.5$ a.u.) does improve the performance of the KDM method quite
dramatically. This, however, comes at a cost of requiring larger Trotter
circuits for a single time step to maintain an equivalent Trotter error. We
should also point out that within all of the numerical experiments, the
presented results correspond to a fixed number of time steps across all
interatomic distances. In practice, the quantum-classical algorithm could have
a predetermined stopping criterion, as illustrated in Figure 1, which would
imply that every distance point could terminate with a different number of
time steps. This would also hold for the hyperparameter optimization loop,
which may be implemented independently for every time step and distance point.
In summary, while both methods provide adequate excited-state energy
estimation, the FDM method might be the preferable choice and perhaps even the
only option for reaching a high level of accuracy based on these numerical
simulations.

## Discussion

##### Single Fidelity Estimation Protocols.

In the following, we discuss various single fidelity estimation protocols that
are relevant to the MFE protocol presented above. Over the past few years,
there have been a wide range of resource-efficient fidelity estimation
protocols that have been developed for the problem of measuring the state
fidelity between two pure states $\ket{\phi_{i}}$ and $\ket{\phi_{j}}$. Here,
we highlight a few of these methods which could be used to provide a
substantial reduction in the overall circuit depth in conjunction with the
multi-fidelity estimation protocol, specifically when compared to the standard
Hadamard test. Before discussing modern approaches, it is worth noting that in
the context of fault-tolerant quantum computing, it is sufficient to use the
SWAP test to estimate the state fidelity
$|\\!\braket{\phi_{i}}{\phi_{j}}\\!|^{2}$, to precision $\epsilon$ with
$\mathcal{O}(1/\epsilon^{2})$ repeated measurements. However, the original
implementation of the SWAP test requires a $2N+1$ qubit circuit, where two
$N$-qubit registers are used for storing $\ket{\phi_{i}}$ and $\ket{\phi_{j}}$
with a controlled-SWAP gate that uses an additional ancilla qubit.

####  II.1.1 Destructive SWAP Test

As shown in Garcia-Escartin and Chamorro-Posada (2013); Cincio  _et al._
(2018), the destructive SWAP test achieves the same outcome as the original
SWAP test but without the use of an ancilla qubit. This approach uses set of
parallel Bell measurements and classical post-processing to achieve the same
result as the original SWAP test with constant depth, $\mathcal{O}(1)$. As a
concrete example, we provide the quantum circuit that implements the
destructive SWAP test for two 3-qubit registers:

$\Qcircuit@C=2em@R=1em@!R{\lstick{\ket{0}}&\multigate{2}{U_{j}}\ctrl{4}\gate{H}\qw\qw\meter\\\
\lstick{\ket{0}}\ghost{U_{j}}\qw\ctrl{4}\gate{H}\qw\meter\\\
\lstick{\ket{0}}\ghost{U_{j}}\qw\qw\ctrl{4}\gate{H}\meter\\\ \\\
\lstick{\ket{0}}\multigate{2}{U_{i}}\targ\qw\qw\qw\meter\\\
\lstick{\ket{0}}\ghost{U_{i}}\qw\targ\qw\qw\meter\\\
\lstick{\ket{0}}\ghost{U_{i}}\qw\qw\targ\qw\meter}$

These Bell measurements only require one-to-one connectivity between each of
the qubits in each of the two registers (as shown above). In principle, this
might be a more hardware-friendly connectivity compared to the Hadamard test
quantum circuit which requires one-to-all connectivity between the ancilla
qubit and the rest of the $N$-qubit register. In addition, the destructive
SWAP test also allows for an additional reduction in the gate depth by
allowing the implementation of Trotterized quantum circuits where only $n/2$
time steps are implemented on each register. Overall, the destructive SWAP
test avoids the use of controlled multi-qubit unitaries, and reduces the depth
of the time-stepping circuits by a half, resulting in a substantial reduction
in overall circuit depth compared to the original Hadamard test.

####  II.1.2 Mirror-type Quantum Circuits

As discussed in Havlíček  _et al._ (2019), it is also possible to estimate
pure-state fidelities by using the following mirror-type quantum circuits
(shown for a 3-qubit register):

$\begin{array}[]{c}\Qcircuit<EMAIL_ADDRESS>\lstick{\ket{0}}&\ghost{U_{j}}&\ghost{U^{\dagger}}&\qw&\meter\\\
\lstick{\ket{0}}&\ghost{U_{j}}&\ghost{U^{\dagger}}&\qw&\meter}\end{array}$

By assuming that the all-zero state $\ket{0}^{\otimes N}$ is the initial state
of the quantum circuit, it is then possible to estimate $F_{1}$ by measuring
the transition probability of returning to the same state at the output, i.e.
measuring $F_{1}=|\\!\braket{0}{U_{i}^{\dagger}U_{j}}{0}\\!|^{2}$. The same
type of circuit can be used to measure the second fidelity $F_{2}$, however,
$U_{i}$ and $U_{j}$ must be replaced with the appropriate circuits that
prepare $\tfrac{1}{\sqrt{2}}(\ket{0}^{\otimes N}+\ket{\phi_{i}})$ and
$\tfrac{1}{\sqrt{2}}(\ket{0}^{\otimes N}+\ket{\phi_{j}})$ respectively. While
this approach requires higher depth quantum circuits, it does not require a
second qubit register as in the destructive SWAP test.

####  II.1.3 Randomized Measurements

Fidelity estimation protocols involving randomized measurements Flammia and
Liu (2011); Elben  _et al._ (2019, 2020); Zhu  _et al._ (2021); Huang  _et
al._ (2020) may also be used for the proposed quantum Krylov subspace
algorithms. These types of protocols simply require preparing the quantum
state $\ket{\phi_{k}}=U_{k}\ket{0}^{\otimes N}$ ($k=i,j$), and applying a
randomized measurement quantum circuit at the output, labelled
$U_{\text{meas}}$, as shown below for a three-qubit register.

$\begin{array}[]{c}\Qcircuit<EMAIL_ADDRESS>\lstick{\ket{0}}&\ghost{U_{k}}&\ghost{U_{\text{meas}}}&\qw&\meter\\\
\lstick{\ket{0}}&\ghost{U_{k}}&\ghost{U_{\text{meas}}}&\qw&\meter}\end{array}$

The choice of measurement quantum circuit will result in different approximate
classical descriptions of the quantum state $\ket{\phi_{k}}$, a so-called
_classical shadow_ , which can be used to estimate an exponential number of
local observables as well as state fidelities with a small number of
measurements Huang  _et al._ (2020). For instance, it was shown in Huang  _et
al._ (2020) that $N$-qubit Clifford circuits achieve a precision $\epsilon$
with $1/\epsilon^{2}$ measurements, while a measurement circuit consisting of
randomized single-qubit Pauli gates achieves a precision $\epsilon$ with
$4^{k}/\epsilon^{2}$ measurements, i.e. scaling exponentially with the
locality $k$ of the observable which is relevant for fidelity estimation.
While the Clifford circuit approach is the most powerful, it comes with
increased depth requirements, making it less desirable for near-term quantum
computing applications. In principle, this randomized measurement approach can
reduce the number of calls/measurements $N_{K}$ even further from the estimate
that we provided in Table 1, though a more careful analysis of the trade-off
between the circuit depth (required for Clifford circuits) and exponential
number of samples (as required by the Pauli circuits) would need to be
performed in order to understand the benefits of this approach in the near
term.

## Concluding Remarks

Building on important previous work on the eigenpair problem on quantum
computers Parrish and McMahon (2019); Stair  _et al._ (2020), and the
classical filter diagonalization method Wall and Neuhauser (1995); Mandelshtam
and Taylor (1997); Mandelshtam (2001), we have presented a unified view of
quantum subspace diagonalization methods including the introduction of three
new generalized eigenvalue problems that can be used to estimate the ground
and excited-state energies of quantum many-body Hamiltonians. Numerical
illustrations of the approaches were carried out for three quantum chemistry
problems. Each of these methods provide various advantages and disadvantages
in terms of number of calls to the quantum computer, gate depth, classical
complexity, and numerical stability on the classical computer portion of the
algorithms. A key new aspect of our approach is a multi-fidelity estimation
protocol that avoids the use of the Hadamard test to estimate the off-diagonal
subspace matrix elements, which allows for a substantial reduction in gate
depth.

Overall, we observed that all of the methods, when tuned properly, converged
to the correct eigenstate and provided excellent estimates of the energy
eigenvalues with a small number time steps for both the ground and excited
state energy estimation problems. While the KDM H and FDM H methods require
less circuit depth (equivalent to a single Trotter time step), the KDM U and
FDM U methods require a lot less calls to the quantum computer, scaling as
$\mathcal{O}(M/\epsilon^{2})$ compared to the
$\mathcal{O}(LM^{p}/\epsilon^{2})$ scaling. This difference is particularly
striking for quantum chemistry Hamiltonians in a localized basis where $L$
scales as $\mathcal{O}(N^{4})$ where $N$ is equal to the total number of spin-
orbitals. For instance, a 50 spin-orbital problem may have approximately
$10^{6}$ terms in the Hamiltonian, which is substantial. This not only
increases the total time to solution, but also increases the overall cost of
the algorithm when the quantum computer requires a per-shot fee.

We also found that by properly tuning the FDM hyperparameters, both FDM
methods converged as quickly as the KDM methods while also remaining much more
stable over a large number of time steps. As a result, we found that the FDM
algorithms achieved a much smaller energy error and variance magnitude
compared to the KDM methods as the number of time steps increased. The FDM
methods also provide a clear advantage in their classical complexity which
might make them more appealing for much larger problems. For example, it is
quite possible that a multi-reference method with tens of thousands (or over a
million) initial reference states could be used to ensure that the initial
wavefunction has a strong overlap with the true ground state wavefunction. In
such a case, the classical cost of solving the generalized eigenvalue problem
would become a severe limitation for the KDM methods, while the FDM methods
would not suffer from such problems. The FDM method should also provide an
advantage in finding interior excited-state energies (eigenvalues) where a
large number of time steps would be required and the filtering technique
becomes more important.

It should be noted that the QKSD algorithms, and in particular the ones we
introduced based on diagonalization of the time evolution propagator (the “U”
forms) share some commonalities with the quantum phase estimation (QPE)
algorithm, as also pointed by Klymko _et al._ Klymko  _et al._ (2021). QPE is
similar in spirit to Fourier analysis of a correlation function and as a
consequence one often sees arguments that the maximum propagation times
(equivalent to maximum circuit depth) must be on the order of of
$\mathcal{O}(1/\epsilon)$ to achieve an eigenvalue estimation accurate to
$\epsilon$ McArdle  _et al._ (2020); Bauer  _et al._ (2020). This would
suggest propagation times several orders of magnitude greater would be
required for QPE algorithms compared to the short propagation times shown here
for QKSD algorithms, which was also confirmed in Klymko  _et al._ (2021).
Variations on QPE that use more sophisticated signal processing ideas O’Brien
_et al._ (2019); Somma (2019) could potentially overcome the
$\mathcal{O}(1/\epsilon)$ limit. Such methods effectively use prior knowledge
about the form of the correlation function or related quantities and, as a
result, the QKSD approaches discussed here could also be put in this class.

Lastly, it is worth mentioning that the proposed algorithms can also be used
in conjunction with variational quantum algorithms that use parameterized
quantum circuits. This approach would use the final optimized quantum circuit,
obtained from a specific variational quantum algorithm, as the initial state
$\ket{\phi_{o}}$ and would subsequently build a Krylov basis to improve the
energy or cost function estimate. A more detailed analysis of the effect of
Trotter error, noise, as well as the development of algorithmic-specific error
mitigation strategies might also represent important next steps for future
work.

##  III ACKNOWLEDGMENTS

This material is based upon work supported by Laboratory Directed Research and
Development (LDRD) funding from Argonne National Laboratory, provided by the
Director, Office of Science, of the U.S. Department of Energy under Contract
No. DE-AC02-06CH11357. Use of the Center for Nanoscale Materials, an Office of
Science user facility, was supported by the U.S. Department of Energy, Office
of Science, Office of Basic Energy Sciences, under Contract No. DE-
AC02-06CH11357. We thank Dmitry A. Fedorov and Yuri Alexeev for helpful
discussions.

## References

  * Lloyd (1996) Seth Lloyd, “Universal quantum simulators,” Science , 1073–1078 (1996).
  * Abrams and Lloyd (1997) Daniel S Abrams and Seth Lloyd, “Simulation of many-body fermi systems on a universal quantum computer,” Physical Review Letters 79, 2586 (1997).
  * Zalka (1998) Christof Zalka, “Efficient simulation of quantum systems by quantum computers,” Fortschritte der Physik: Progress of Physics 46, 877–879 (1998).
  * Aspuru-Guzik  _et al._ (2005) Alán Aspuru-Guzik, Anthony D Dutoi, Peter J Love, and Martin Head-Gordon, “Simulated quantum computation of molecular energies,” Science 309, 1704–1707 (2005).
  * Nielsen and Chuang (2002) Michael A Nielsen and Isaac Chuang, “Quantum computation and quantum information,” (2002).
  * Kitaev (1997) Aleksei Yur’evich Kitaev, “Quantum computations: algorithms and error correction,” Uspekhi Matematicheskikh Nauk 52, 53–112 (1997).
  * Cleve  _et al._ (1998) Richard Cleve, Artur Ekert, Chiara Macchiavello, and Michele Mosca, “Quantum algorithms revisited,” Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 339–354 (1998).
  * Peruzzo  _et al._ (2014) Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications 5, 1–7 (2014).
  * McClean  _et al._ (2016) Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik, “The theory of variational hybrid quantum-classical algorithms,” New Journal of Physics 18, 023023 (2016).
  * Kandala  _et al._ (2017) Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M Chow, and Jay M Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature 549, 242–246 (2017).
  * Farhi  _et al._ (2014) Edward Farhi, Jeffrey Goldstone, and Sam Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028 (2014).
  * Preskill (2018) John Preskill, “Quantum computing in the nisq era and beyond,” Quantum 2, 79 (2018).
  * Cerezo  _et al._ (2021a) Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, _et al._ , “Variational quantum algorithms,” Nature Reviews Physics , 1–20 (2021a).
  * Santagati  _et al._ (2018) Raffaele Santagati, Jianwei Wang, Antonio A Gentile, Stefano Paesani, Nathan Wiebe, Jarrod R McClean, Sam Morley-Short, Peter J Shadbolt, Damien Bonneau, Joshua W Silverstone, _et al._ , “Witnessing eigenstates for quantum simulation of hamiltonian spectra,” Science advances 4, eaap9646 (2018).
  * Higgott  _et al._ (2019) Oscar Higgott, Daochen Wang, and Stephen Brierley, “Variational quantum computation of excited states,” Quantum 3, 156 (2019).
  * Sim  _et al._ (2018) Sukin Sim, Jonathan Romero, Peter D Johnson, and Alán Aspuru-Guzik, “Quantum computer simulates excited states of molecule,” Physics 11, 14 (2018).
  * Bittel and Kliesch (2021) Lennart Bittel and Martin Kliesch, “Training variational quantum algorithms is np-hard–even for logarithmically many qubits and free fermionic systems,” arXiv preprint arXiv:2101.07267 (2021).
  * McClean  _et al._ (2018) Jarrod R McClean, Sergio Boixo, Vadim N Smelyanskiy, Ryan Babbush, and Hartmut Neven, “Barren plateaus in quantum neural network training landscapes,” Nature Communications 9, 1–6 (2018).
  * Wang  _et al._ (2020) Samson Wang, Enrico Fontana, Marco Cerezo, Kunal Sharma, Akira Sone, Lukasz Cincio, and Patrick J Coles, “Noise-induced barren plateaus in variational quantum algorithms,” arXiv preprint arXiv:2007.14384 (2020).
  * Cerezo  _et al._ (2021b) Marco Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J Coles, “Cost function dependent barren plateaus in shallow parametrized quantum circuits,” Nature Communications 12, 1–12 (2021b).
  * Arrasmith  _et al._ (2020) Andrew Arrasmith, M Cerezo, Piotr Czarnik, Lukasz Cincio, and Patrick J Coles, “Effect of barren plateaus on gradient-free optimization,” arXiv preprint arXiv:2011.12245 (2020).
  * Huggins  _et al._ (2020) William J Huggins, Joonho Lee, Unpil Baek, Bryan O’Gorman, and K Birgitta Whaley, “A non-orthogonal variational quantum eigensolver,” New Journal of Physics 22, 073009 (2020).
  * Parrish and McMahon (2019) Robert M Parrish and Peter L McMahon, “Quantum filter diagonalization: Quantum eigendecomposition without full quantum phase estimation,” arXiv preprint arXiv:1909.08925 (2019).
  * Seki and Yunoki (2021) Kazuhiro Seki and Seiji Yunoki, “Quantum power method by a superposition of time-evolved states,” PRX Quantum 2, 010333 (2021).
  * Stair  _et al._ (2020) Nicholas H Stair, Renke Huang, and Francesco A Evangelista, “A multireference quantum krylov algorithm for strongly correlated electrons,” Journal of Chemical Theory and Computation 16, 2236–2245 (2020).
  * McClean  _et al._ (2017) Jarrod R McClean, Mollie E Kimchi-Schwartz, Jonathan Carter, and Wibe A De Jong, “Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states,” Physical Review A 95, 042308 (2017).
  * Colless  _et al._ (2018) James I Colless, Vinay V Ramasesh, Dar Dahlen, Machiel S Blok, Mollie E Kimchi-Schwartz, Jarrod R McClean, Jonathan Carter, Wibe A de Jong, and Irfan Siddiqi, “Computation of molecular spectra on a quantum processor with an error-resilient algorithm,” Physical Review X 8, 011021 (2018).
  * Bharti and Haug (2020a) Kishor Bharti and Tobias Haug, “Iterative quantum assisted eigensolver,” arXiv preprint arXiv:2010.05638 (2020a).
  * Bharti and Haug (2020b) Kishor Bharti and Tobias Haug, “Quantum assisted simulator,” arXiv preprint arXiv:2011.06911 (2020b).
  * Golub and Van Loan (1996) Gene H. Golub and Charles F. Van Loan,  _Matrix Computations_ , 3rd ed. (The Johns Hopkins University Press, 1996).
  * Parlett (1980) Beresford N. Parlett,  _The Symmetric Eigenvalue Problem_ (Prentice-Hall, Englewood Cliffs, N.J., 1980).
  * Bespalova and Kyriienko (2021) Tatiana A Bespalova and Oleksandr Kyriienko, “Hamiltonian operator approximation for energy measurement and ground-state preparation,” PRX Quantum 2, 030318 (2021).
  * Motta  _et al._ (2020) Mario Motta, Chong Sun, Adrian TK Tan, Matthew J O’Rourke, Erika Ye, Austin J Minnich, Fernando GSL Brandão, and Garnet Kin-Lic Chan, “Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution,” Nature Physics 16, 205–210 (2020).
  * Wall and Neuhauser (1995) Michael R Wall and Daniel Neuhauser, “Extraction, through filter-diagonalization, of general quantum eigenvalues or classical normal mode frequencies from a small number of residues or a short-time segment of a signal. i. theory and application to a quantum-dynamics model,” The Journal of Chemical Physics 102, 8011–8022 (1995).
  * Mandelshtam and Taylor (1997) Vladimir A Mandelshtam and Howard S Taylor, “A low-storage filter diagonalization method for quantum eigenenergy calculation or for spectral analysis of time signals,” The Journal of Chemical Physics 106, 5085–5090 (1997).
  * Mandelshtam (2001) Vladimir A Mandelshtam, “FDM: the filter diagonalization method for data processing in nmr experiments,” Progress in Nuclear Magnetic Resonance Spectroscopy 38, 159–196 (2001).
  * Klymko  _et al._ (2021) Katherine Klymko, Carlos Mejuto-Zaera, Stephen J Cotton, Filip Wudarski, Miroslav Urbanek, Diptarka Hait, Martin Head-Gordon, K Birgitta Whaley, Jonathan Moussa, Nathan Wiebe, _et al._ , “Real time evolution for ultracompact hamiltonian eigenstates on quantum hardware,” arXiv preprint arXiv:2103.08563 (2021).
  * Low and Chuang (2019) Guang Hao Low and Isaac L Chuang, “Hamiltonian simulation by qubitization,” Quantum 3, 163 (2019).
  * Poulin  _et al._ (2018) David Poulin, Alexei Kitaev, Damian S Steiger, Matthew B Hastings, and Matthias Troyer, “Quantum algorithm for spectral measurement with a lower gate count,” Physical Review Letters 121, 010501 (2018).
  * Gard  _et al._ (2020) Bryan T Gard, Linghua Zhu, George S Barron, Nicholas J Mayhall, Sophia E Economou, and Edwin Barnes, “Efficient symmetry-preserving state preparation circuits for the variational quantum eigensolver algorithm,” npj Quantum Information 6, 1–9 (2020).
  * Garcia-Escartin and Chamorro-Posada (2013) Juan Carlos Garcia-Escartin and Pedro Chamorro-Posada, “Swap test and hong-ou-mandel effect are equivalent,” Physical Review A 87, 052330 (2013).
  * Cincio  _et al._ (2018) Lukasz Cincio, Yiğit Subaşı, Andrew T Sornborger, and Patrick J Coles, “Learning the quantum algorithm for state overlap,” New Journal of Physics 20, 113022 (2018).
  * Havlíček  _et al._ (2019) Vojtěch Havlíček, Antonio D Córcoles, Kristan Temme, Aram W Harrow, Abhinav Kandala, Jerry M Chow, and Jay M Gambetta, “Supervised learning with quantum-enhanced feature spaces,” Nature 567, 209–212 (2019).
  * Flammia and Liu (2011) Steven T Flammia and Yi-Kai Liu, “Direct fidelity estimation from few pauli measurements,” Physical Review Letters 106, 230501 (2011).
  * Elben  _et al._ (2019) Andreas Elben, Benoît Vermersch, Christian F Roos, and Peter Zoller, “Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states,” Physical Review A 99, 052323 (2019).
  * Elben  _et al._ (2020) Andreas Elben, Benoît Vermersch, Rick van Bijnen, Christian Kokail, Tiff Brydges, Christine Maier, Manoj K Joshi, Rainer Blatt, Christian F Roos, and Peter Zoller, “Cross-platform verification of intermediate scale quantum devices,” Physical Review Letters 124, 010504 (2020).
  * Zhu  _et al._ (2021) Daiwei Zhu, Ze-Pei Cian, Crystal Noel, Andrew Risinger, Debopriyo Biswas, Laird Egan, Yingyue Zhu, Alaina M Green, Andrii Maksymov, Yunseong Nam, _et al._ , “Cross-platform comparison of arbitrary quantum computations,” arXiv preprint arXiv:2107.11387 (2021).
  * Huang  _et al._ (2020) Hsin-Yuan Huang, Richard Kueng, and John Preskill, “Predicting many properties of a quantum system from very few measurements,” Nature Physics 16, 1050–1057 (2020).
  * McArdle  _et al._ (2020) Sam McArdle, Suguru Endo, Alán Aspuru-Guzik, Simon C. Benjamin, and Xiao Yuan, “Quantum computational chemistry,” Rev. Mod. Phys. 92, 015003 (2020).
  * Bauer  _et al._ (2020) Bela Bauer, Sergey Bravyi, Mario Motta, and Garnet Kin-Lic Chan, “Quantum algorithms for quantum chemistry and quantum materials science,” Chemical Reviews 120, 12685–12717 (2020).
  * O’Brien  _et al._ (2019) Thomas E O’Brien, Brian Tarasinski, and Barbara M Terhal, “Quantum phase estimation of multiple eigenvalues for small-scale (noisy) experiments,” New Journal of Physics 21, 023022 (2019).
  * Somma (2019) Rolando D Somma, “Quantum eigenvalue estimation via time series analysis,” New Journal of Physics 21, 123025 (2019).
  * Sun  _et al._ (2018) Qiming Sun, Timothy C Berkelbach, Nick S Blunt, George H Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D McClain, Elvira R Sayfutyarova, Sandeep Sharma, _et al._ , “PySCF: the python-based simulations of chemistry framework,” Wiley Interdisciplinary Reviews: Computational Molecular Science 8, e1340 (2018).
  * McClean  _et al._ (2020) Jarrod R McClean, Nicholas C Rubin, Kevin J Sung, Ian D Kivlichan, Xavier Bonet-Monroig, Yudong Cao, Chengyu Dai, E Schuyler Fried, Craig Gidney, Brendan Gimby, _et al._ , “OpenFermion: the electronic structure package for quantum computers,” Quantum Science and Technology 5, 034014 (2020).
  * Bergholm  _et al._ (2018) Ville Bergholm, Josh Izaac, Maria Schuld, Christian Gogolin, M Sohaib Alam, Shahnawaz Ahmed, Juan Miguel Arrazola, Carsten Blank, Alain Delgado, Soran Jahangiri, _et al._ , “Pennylane: Automatic differentiation of hybrid quantum-classical computations,” arXiv preprint arXiv:1811.04968 (2018).
  * Press  _et al._ (2007) William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery,  _Numerical Recipes 3rd Edition: The Art of Scientific Computing_ , 3rd ed. (Cambridge University Press, Cambridge, 2007).

## Appendix A: Quantum chemistry and molecular systems studied

Within the Born-Oppenheimer approximation, a molecule is comprised of $\eta$
electrons interacting within a potential produced by nuclei at fixed
positions. Using the second quantization formalism, the problem may be cast in
terms of $N$ single-particle spin orbitals that can be occupied or empty. In
the absence of external fields, the non-relativistic molecular Hamiltonian is
written as:

$H=h_{\text{nuc}}+\sum_{pq}h_{pq}a_{p}^{\dagger}a_{q}+\tfrac{1}{2}\sum_{pqrs}h_{pqrs}a^{\dagger}_{p}a^{\dagger}_{q}a_{r}a_{s}$
(15)

where $a_{p}$ and $a_{p}^{\dagger}$ correspond to fermionic annihilation and
creation operators that obey the anti-commutation relations,
$\\{a_{j},a_{k}\\}=0$, $\\{a_{j}^{\dagger},a_{k}^{\dagger}\\}=0$, and
$\\{a_{j},a_{k}^{\dagger}\\}=\delta_{jk}$. $h_{\text{nuc}}$ corresponds to the
classical electrostatic repulsion between nuclei, while $h_{pq}$ and
$h_{pqrs}$ correspond to one- and two-electron integrals and written
explicitly as:

$\displaystyle h_{\text{nuc}}$ $\displaystyle=\frac{1}{2}\sum_{i\neq
j}\frac{Z_{i}Z_{j}}{|\mathbf{R}_{i}-\mathbf{R}_{j}|},$ (16) $\displaystyle
h_{pq}$
$\displaystyle=\int\\!d\sigma\;\phi_{p}^{*}(\sigma)\left(-\frac{\nabla_{r}^{2}}{2}-\sum_{i}\frac{Z_{i}}{|\mathbf{R}_{i}-\mathbf{r}|}\right)\phi_{q}(\sigma),$
(17) $\displaystyle h_{pqrs}$
$\displaystyle=\int\\!d\sigma_{1}d\sigma_{2}\;\;\frac{\phi_{p}^{*}(\sigma_{1})\phi_{q}^{*}(\sigma_{2})\phi_{r}(\sigma_{1})\phi_{s}(\sigma_{2})}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}.$
(18)

where the summations run over all nuclei. $Z_{i}$ represents the nuclear
charge, $\mathbf{r}$ and $\mathbf{R}$ denote electronic and nuclear spatial
coordinates, and $\sigma$ is a generalized coordinate consisting of the
spatial and spin degrees of freedom, $\sigma_{i}=(\mathbf{r}_{i},s_{i})$. The
function $\phi(\sigma)$ represents a one-electron spin-orbital. In the main
text, we calculate the potential energy surface $E(\mathbf{R})$, which can be
used to describe a wide range of processes such as reaction dynamics, bond-
breaking, and chemical dynamics. In order to predict thermochemical properties
accurately at room temperature and atmospheric pressure, we require an
estimate of the potential energy surface within a chemical accuracy of $1$
kcal/mol or equivalently $1.59\times 10^{-3}$ Hartrees or 43.4 meV,
represented by the horizontal lines in the energy error plots in the
manuscript.

For our illustrative calculations, we consider (i) a 6-atom linear hydrogen
chain, H6, (ii) the bent H2O molecule, and (iii) the linear BeH2 molecule. We
employ minimal STO-3G bases such that there are six spatial molecular orbitals
for H6, and seven molecular orbitals each for H2O and BeH2. The number of
qubits, after performing the Jordan-Wigner transformation is equal to the
number of (active) spin-orbitals being treated. The hydrogen chain can thus be
treated with twelve qubits corresponding to the total number of spin-orbitals
and within full configuration interaction (FCI). For H2O and BeH2 we employ a
complete active space approach wherein there are $n$ active electrons
distributed in $m$ spatial orbitals with a full configuration interaction
being carried out only among the active electrons and orbitals, sometimes
denoted as CASCI$(n,m)$. Generally the $m$ orbitals are the highest energy
ones and the remaining electrons are effectively frozen in the lower, inactive
orbitals. In the case of H2O we use CASCI(8,6) and for BeH2 we use CASCI(6,6),
both corresponding to twelve qubit Hamiltonians in both cases. Of course the
levels of theory here, in particular the use of STO-3G basis sets, is quite
crude by computational chemistry standards, but it is sufficient to illustrate
and compare the various diagonalization procedures of interest to us.

To obtain the coefficients of the second quantized Hamiltonian, we defined the
structure of the molecule by its constituent elements and nuclear coordinates
of the atoms. The initial basis of single-particle states was obtained via the
Hartree-Fock method where each electron is treated as an independent particle
that moves under the influence of the Coulomb potential due to the nuclei as
well as a mean field generated by all of the other electrons. We used the
quantum chemistry package PySCF Sun  _et al._ (2018) to obtain the optimized
coefficients of the linear combination of the atomic orbitals, obtained via
the Hartree-Fock method. To map the fermionic second quantization Hamiltonian
to a qubit-basis, it is possible to use the Jordan-Wigner, Parity, or Bravyi-
Kitaev transformations. We used the OpenFermion Python package McClean  _et
al._ (2020) to obtain the qubit representation of the Hamiltonian, where the
Jordan-Wigner basis was chosen for all of the numerical simulations. All of
these packages can also be found within the Pennylane cross-platform Python
library Bergholm  _et al._ (2018).

## Appendix B: Choice of time step

First we consider the case that the generalized eigenvalue problem involves
$f(\hat{H})$ = $\hat{H}$. If an absolute energy scale is used for the
corresponding Hamiltonian matrix such that the energy eigenvalues of interest
could be large in magnitude, the time step $\tau$ should ideally be equal to
or smaller than $\pi/\Delta E$ where $\Delta E$ is an upper bound to the
spectral range (difference between maximum and minimum eigenvalues) and atomic
units ($\hbar$ = 1) are assumed. This parallels considerations of discrete
Fourier transforms wherein the possible range of eigenvalues is determined
from the Nyquist interval (in terms of angular frequencies or energies),
$-\pi/\tau$ to $+\pi/\tau$ Press  _et al._ (2007). Note that it is generally
not difficult to make such an upper limit estimate to the spectral range given
some knowledge of the problem, e.g., the Hartree-Fock and molecular orbital
energies in the case of the electronic structure problem.

If sufficient accuracy can still be achieved with the corresponding Trotter
time steps, however, one can actually use larger $\tau$ values than
$\pi/\Delta E$ if the energy zero of the Hamiltonian is simply shifted to be
in the center of the desired energy range, i.e. the Hamiltonian in all
equations is taken to be the energy-shifted one. Aliasing, i.e., mapping of
eigenvalues outside the corresponding Nyquist interval into it can be easily
and cheaply monitored for: a slightly different $\tau$ will yield the same
“true” eigenvalues but the aliased ones will shift. Actually, in the three
chemical systems considered in this paper there were no issues with aliasing
at all because even with the energy shifting permitting the larger $\tau$
values, the Nyquist interval was still tens of Hartrees wide and sufficient to
capture all relevant eigenvalues present.

Similar considerations to the above apply when the generalized eigenvalue
problem involves $f(\hat{H})$ = $\exp(-i\hat{H}\tau)$. That is, ideally $\tau$
$\leq$ $\pi/\Delta E$ where $\Delta E$ is an upper bound to the Hamiltonian’s
spectral range but, with appropriate recognition of aliasing issues larger
$\tau$ could be used. We note that in this case the eigenvalues $f(E_{k})$ are
complex and lie on the unit circle. Let $\theta_{k}$ =
$\tan^{-1}(-\rm{Im}~{}f(E_{k})/\rm{Re}~{}f(E_{k}))$ be the angle between
$-\pi$ and $\pi$ that results from typical numerical arctangent function
(“$\rm{atan2(y,x)}$”) calls. The possible energy eigenvalues are $E_{k}$ =
$(\theta_{k}+j2\pi)/\tau$, where $j$ = 0, $\pm 1$, $\pm 2$, etc. If some
absolute energy scale that leads to large magnitude physical energy
eigenvalues is used, then $|j|$ $\geq$ 1 value may have to be used to obtain
energy eigenvalues in the desired physical range. The Hamiltonian energy
shifting noted above is particularly convenient for this case because then,
typically, it is not necessary to shift the eigenvalues and $j$ = 0 suffices.

## Appendix C: Numerical solution of the generalized eigenvalue problem

The numerical simulations shown in Figures 2 to 4 were performed with an in-
house Python code using standard NumPy and SciPy numerical linear algebra
packages. The classical solution of the complex generalized eigenvalue
problem, $\bf{F}\bf{c}$. = $f\bf{S}\bf{c},$ was performed using the QZ
algorithm or generalized Schur decomposition Golub and Van Loan (1996) of the
complex matrices $\bf{F}$ and $\bf{S}$, which we found to be much more stable
compared to conventional generalized eigenvalue solvers.

An alternative, more “hands on” approach to solving the generalized eigenvalue
problem is to carry out singular value decomposition Golub and Van Loan (1996)
of the overlap matrix, $\bf{S}$, thereby allowing construction of its inverse
and turning the problem into an eigenvalue problem for a complex matrix of the
form, e.g., $\bf{S}^{-1}\bf{F}\bf{c}$ = $f\bf{c}$. The latter problem can be
numerically solved with standard numerical software but may require
appropriate zeroing of some of the singular values when the inverse is
constructed.

## Appendix D: Overview of Hadamard Test

In the following, we describe the Hadamard test used for estimating the matrix
element, $\braket{\phi_{i}}{U}{\phi_{i}}$. This method uses an ancilla qubit
initially prepared in the $\ket{0}$ state, with the rest of the quantum
circuit shown below:

$\Qcircuit@C=1em@R=1em{\lstick{\ket{0}}&\gate{H}\ctrl{1}\meter\qw\\\
\lstick{\ket{\phi_{i}}}\qw{/}\gate{U}\qw}$

This type of approach will effectively require an ancilla with one-to-all
connectivity in order to keep the circuit gate depth as short as possible. The
real and imaginary parts of the matrix element are obtained by measuring the
ancilla qubit in the X and Y Pauli basis, where it can be shown that
$\braket{\sigma_{x}+i\sigma_{y}}=\braket{\phi_{i}}{U}{\phi_{i}}$. This type of
Hadamard quantum circuit is sufficient for single-reference-based Krylov
subspace algorithms. For more general off-diagonal matrix elements of the
form, $\braket{\phi_{i}}{\phi_{j}}=\braket{0}{U_{i}^{\dagger}U_{j}}{0}$, as
required by multi-reference Krylov subspace algorithms Stair  _et al._ (2020),
the following quantum circuit would be used:

$\Qcircuit@C=1em@R=1em@!R{\lstick{\ket{0}\;\;\;}&\gate{H}\gate{X}\ctrl{1}\gate{X}\ctrl{1}\meter\qw\\\
\lstick{\ket{0}^{\small{\otimes
N}}\\!\\!\\!}\qw{/}\qw\gate{U_{j}}\qw\gate{U_{i}}\qw}$

where $U_{i}$ and $U_{j}$ correspond to the quantum circuits that prepare
$\ket{\phi_{i}}$ and $\ket{\phi_{j}}$, defined through the relations
$\ket{\phi_{i}}=U_{i}\ket{0}$ and $\ket{\phi_{j}}=U_{j}\ket{0}$ respectively.
The real and imaginary parts are obtained through the measurement of the
ancilla qubit in the $X$ and $Y$ Pauli basis as before.

## Appendix E: Details of excited-state energy calculations

In the main manuscript, we performed excited-state energy calculations for a
water molecule as a function of the bond length $R$ between the oxygen and
hydrogen atoms. We only considered the singlet excited-state energies where
the total angular momentum is zero. In principle, it is possible to use the
Hartree-Fock state as the initial starting state $\ket{\phi_{o}}$ since it
corresponds to a singlet state and, by symmetry, it will be connected to the
excited singlet states. However, the convergence towards the excited-state
energies will be slow because the Hartree-Fock state is nearly orthogonal to
the desired excited-state wavefunctions. Here, we propose the use of
physically-motivated ansatz states to find the first four singlet excited-
state energies of water. The ansatze consists of the following four states:

$\displaystyle\ket{\Phi_{0}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\left(\ket{000110111111}-\ket{001001111111}\right)$
(19) $\displaystyle\ket{\Phi_{1}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\left(\ket{000111101111}-\ket{001011011111}\right)$
(20) $\displaystyle\ket{\Phi_{2}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\left(\ket{010010111111}-\ket{100001111111}\right)$
(21) $\displaystyle\ket{\Phi_{3}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\left(\ket{010011101111}-\ket{100011011111}\right)$
(22)

where we use an alternating alpha (spin up) beta (spin down) ordering for the
spin-orbitals with increasing energies from right to left. Note that each of
these states has zero total angular momentum, therefore, they will be
decoupled from any triplet states that might otherwise slow down the
convergence of the Krylov-based algorithm. For the calculations, we used
states (18)-(21) as the initial single-reference starting state
$\ket{\phi_{o}}$. Each of the calculations are independent from one another.
To perform the multi-fidelity estimation protocol from the main manuscript,
each of these four states would require a state preparation step that
initializes the superposition state,
$\tfrac{1}{\sqrt{2}}(\ket{0}^{\otimes{N}}+\ket{\Phi_{k}})$, where $k=0,1,2,3$,
on the quantum computer.