# Block Lanczos method for excited states on a quantum computer

Thomas E. Baker  Department of Physics, University of York, Heslington,York
YO10 5DD, United Kingdom

###### Abstract

The method of quantum Lanczos recursion is extended to solve for multiple
excitations on the quantum computer. While quantum Lanczos recursion is in
principle capable of obtaining excitations, the extension to a block Lanczos
routine can resolve degeneracies with better precision and only costs
$\mathcal{O}(d^{2})$ for $d$ excitations on top of the previously introduced
quantum Lanczos recursion method. The formal complexity in applying all
operators to the system at once with oblivious amplitude amplification is
exponential, but this cost can be kept small to obtain the ground state by
incrementally adding operators. The error of the ground state energy based on
the accuracy of the Lanczos coefficients is investigated and the error of the
ground state energy. It is demonstrated to scale linearly with the uncertainty
of the Lanczos coefficients. Extension to non-Hermitian operators is also
discussed.

##  I Introduction

For quantum computation to lead to new discoveries Nielsen and Chuang (2010),
efficient means of solving for the ground state must be understood and
implemented. Some near term algorithms that have been used in the era of noisy
quantum devices have led to an increased interest in results obtained from
variational quantum eigensolvers O’Malley  _et al._ (2016); Kandala  _et al._
(2017); Gard  _et al._ (2020); Gokhale  _et al._ (2020); Ollitrault  _et al._
(2020); Klymko  _et al._ (2021); Kremenetski  _et al._ (2021) which ultimately
face noise and other limitations Bittel and Kliesch (2021). Looking to the
long term capabilities of a quantum computer when error correction is
available, other algorithms are useful to investigate.

By far the longest studied algorithms to obtain ground states is the
implementation of real time evolution. In this algorithm, an initial
Hamiltonian $\mathcal{H}_{0}$ is defined and an initial state $\Psi_{0}$ is
prepared on the quantum computer. The time dependent Hamiltonian is

$\mathcal{H}(t)=\mathcal{H}_{0}+\lambda(t)\mathcal{H}_{I}+\mathcal{C}$ (1)

with interaction term $\mathcal{H}_{I}$, constant $\mathcal{C}$, and time
dependent coupling constant (or similar form) $\lambda(t)$ at time $t$. By
adiabatically (slowly) increasing the interaction term in time, the
wavefunction will eventually arrive at the ground state for the fully
interacting Hamiltonian.

However, this solution strategy is known to be extremely slow for quantum
chemical systems Poulin  _et al._ (2015); Lemieux  _et al._ (2021). In order
to apply the Hamiltonian, a time evolution operator of the form
$\exp(-i\mathcal{H}(t)\delta t)$ must be applied to the wavefunction. The
Trotter-Suzuki decomposition of the time evolution operator must be decomposed
into many terms, $\mathcal{O}(N^{4})$ to capture the full electron-electron
interaction term, although this can be reduced as $N\rightarrow\infty$ to
$\mathcal{O}(N^{2})$ for the case of local basis functions Baker and Poulin
(2020). However, since the time step $\delta t$ must be very small depending
on the strength of correlation in the system, the resulting number of
operations makes the time necessary to solve for even small molecules
extremely long. This is true of other classical solution techniques such as
Hartree-Fock Gulania and Whitfield (2021). This is expected based on the
complexity of solving quantum chemistry systems Schuch and Verstraete (2009).

The question of how to obtain excited states is one that has been investigated
in recent papers Higgott  _et al._ (2019); Nakanishi  _et al._ (2019); Parrish
_et al._ (2019); Mizuta  _et al._ (2021); Kuroiwa and Nakagawa (2021); Greene-
Diniz and Muñoz Ramo (2021); Motta  _et al._ (2020); Jahangiri  _et al._
(2020); Ollitrault  _et al._ (2020); Roggero  _et al._ (2020); Sim  _et al._
(2018); Bauman  _et al._ (2019), and a direct solution would provide a means
to fully manipulate the wavefunction. If this can be accomplished without the
use of time evolutions, then avoiding the small time step necessary for those
methods could be possible.

In addition to ground state solvers, obtaining excited states is a highly
valuable quantity for quantum chemistry systems. Historically, excited states
have been more difficult to obtain. So, if a quantum algorithm could reliably
obtain the excitations, then this would represent a major improvement over
existing classical techniques Elliott  _et al._ (2011); Khoromskaia and
Khoromskij (2015); Helgaker  _et al._ (2014).

Recently, Lanczos methods have been investigated in the context of solving
quantum systems. One variety of Lanczos on the quantum computer, uses even
vectors of a Krylov subspace and imaginary time evolution techniques to obtain
the ground state and other quantities of interest Motta  _et al._ (2020); Sun
_et al._ (2021).

Another recently introduced variety of Lanczos algorithms Baker (2021) uses
oblivious amplitude amplification (OAA) Kothari (2014); Berry  _et al._
(2015); Childs  _et al._ (2017)–an algorithm similar to Grover’s search but
run on auxiliary qubits Grover (1996, 2001); Brassard  _et al._ (2002)–to
apply operators directly onto wavefunctions. Since OAA applies operators onto
wavefunctions, all the necessary operations for Lanczos can be implemented
directly. The full Lanczos recursion can be implemented to find the ground
state or the continued fraction representation of the Green’s function Baker
(2021). This second Lanczos technique called quantum Lanczos recursion (QLR)
avoids the use of a time evolution operator and therefore bypasses the
computational bottleneck in terms of the number of terms in the Trotter-Suzuki
decomposition.

One other advantages of the QLR is that the traditional limitations of Lanczos
on the classical computer are entirely circumvented Cahill  _et al._ (2000),
and only one number of wavefunctions. The only limitation is how accurately
the operations can be applied to the wavefunction in practice, and this has
recently been investigated on a quantum computer Jamet  _et al._ (2021).

This paper uses similar techniques as QLR to show that the traditional Lanczos
recursion can be replaced by a block Lanczos routine to resolve several
excitations, leading to quantum block Lanczos recursion (QBLR), although
“block” can be replaced by “banded” or some other name corresponding to a
trivial change of gauge of the unitary matrices involved. While it is
trivially demonstrated from QLR that excitations can be found, the performance
of block Lanczos will allow for the resolving of degeneracies with greater
ease. This was recently demonstrated in tensor network algorithms Baker  _et
al._ (2021a); Di Paolo  _et al._ (2021).

The algorithm here in the can also be performed with the preparation of a
single initialization of starting wavefunctions. A full collapse of the
eigenstates is avoided, meaning that a wavefunction can be preserved for the
next computation. This is accomplished by use of a state-preserving quantum
counting algorithm Baker (2021), referred to in some works as QMA-sampling
Marriott and Watrous (2005); Temme  _et al._ (2011).

Some additional discussion on how errors in the Lanczos coefficients will
affect the ground state energy and otherwise are also considered. The method
of slowly introducing terms into the Hamiltonian is also discussed in the
context of Lanczos to show that the OAA algorithm does not need to scale
exponentially for limited additional terms.

##  II Excitations from quantum Lanczos recursion

The Lanczos algorithm from Ref. Baker, 2021 can be used to find excitations of
a given model by a simple alternative of coefficients. In this section, QLR
will be reviewed and an extension to finding excited states will be shown.

A Lanczos recursion relation to find subsequent elements of the Krylov
subspace, $\\{\psi_{0},\psi_{1},\ldots,\psi_{N}\\}$ is

$|\psi_{n+1}\rangle=\mathcal{H}|\psi_{n}\rangle-\alpha_{n}|\psi_{n}\rangle-\beta_{n}|\psi_{n-1}\rangle$
(2)

where the resulting Hamiltonian in the basis of the Krylov subspace forms a
tridiagonal matrix. When diagonalized, the ground state is found to a high
accuracy even if only a few $n$ are determined.

##### Applying operators to wavefunctions.–

The only ingredient that is necessary for this form of the ground state
solution is to ensure that the operators could be applied to the wavefunction.
The method to do this was introduced in Refs. berry2014exponential; Kothari,
2014 where an oblivious amplitude amplification (OAA) algorithm Berry  _et
al._ (2015); Childs  _et al._ (2017) was used to improve on repeat-until-
success methods which have the same goal Guerreschi (2019); Paetznick and
Svore (2013). Essentially, when an operator is applied to a wavefunction as in
a repeat-until-success strategy, the operation is controlled on some auxiliary
qubits that give 0 if the application of an operator is successful and 1 if it
is not. OAA searches for a state on some auxiliary qubits corresponding to all
0s. By searching over the auxiliary qubit states, the algorithm remains
oblivious to the wavefunction itself. This essentially guarantees the outcome
from the related method of repeat-until-success to apply the operator in one
step Guerreschi (2019); Paetznick and Svore (2013). For how to construct the
operators, see the discussion in Ref. Low and Chuang, 2019.

##### Minimal measurements of the wavefunctions.–

In order to sample the coefficients $\alpha_{n}$ and $\beta_{n}$ without
completely measuring (and therefore destroying) the wavefunction, a state-
preserving quantum counting algorithm can be used. The algorithm applies a
generic operator $\mathcal{A}$ as in the form

$\mathcal{A}|\Psi\rangle=p|\Psi\rangle+p^{\perp}|\Psi^{\perp}\rangle$ (3)

which signifies a superposition of eigenstates. The operator $\mathcal{A}$ is
normalized and represented as a unitary such that it can be represented on the
qubits Baker and Poulin (2020); Baker (2021). The resulting probability $p$ is
the expectation value of $\mathcal{A}$. All states orthogonal to the original
state are marked with a $\perp$ symbol. The basis of the superposition chosen
here is that of the eigenbasis because this will be the natural basis to pick
for quantum phase estimation (QPE) Nielsen and Chuang (2010); Low and Chuang
(2019); Poulin  _et al._ (2018).

The key to finding $p$ is to count the number of transitions from $\Psi$ to
$\Psi$ (wavefunction to same wavefunction) after the application of
$\mathcal{A}$ by OAA or some other algorithm. To verify that the same state is
obtained after this procedure, the energy of the state can be computed with
QPE at the start of the algorithm and stored on a register. After the operator
is applied, the energy is then computed on a separate register. The two
registers are compared and represented on a single qubit, called a pointer
qubit. The pointer qubit is then measured. If the energies match, one value is
returned and the original wavefunction is recovered. This is counted as a
“success” in the algorithm. The ratio of successes to the total times the
algorithm is run is $p$. Once $p$ is found, the expectation value of
$\mathcal{A}$ can be determined.

If the pointer qubit demonstrates that the wrong state was recovered, then a
recovery procedure is used to find the original state Baker (2021); Baker and
Poulin (2020); Temme  _et al._ (2011); Paetznick and Svore (2013). In essence,
the unitary of all operations applied onto the wrong wavefunction and the
process above is repeated until the correct wavefunction is found Temme  _et
al._ (2011).

##### Operators for Lanczos recursion.–

To find the Lanczos coefficients, the following relations hold Baker (2021)

$\alpha_{n}=\langle\psi_{n}|\mathcal{H}|\psi_{n}\rangle\quad\mathrm{and}\quad\beta_{n}=\langle\psi_{n-1}|\mathcal{H}|\psi_{n}\rangle$
(4)

and therefore provide a means to use state-preserving quantum counting to
obtain the coefficients. The map between the Krylov states and the original
ground state is crucial but can be summarized as

$|\psi_{n}\rangle=\hat{G}_{n}|\Psi\rangle$ (5)

where some examples of these operators are shown in Ref. Baker, 2021. The
coefficient $\alpha_{n}$ can be regarded as applying quantum counting on the
state $\Psi$ with the operator
$(\hat{G}_{n}^{\dagger}\mathcal{H}\hat{G}_{n})$. Each operator $\hat{G}_{n}$
depends on coefficients from $\\{0,1,\ldots,n-1\\}$ iterations and therefore
the algorithm discovers the coefficients iteratively. To find $\beta_{n}$, the
operator $(\hat{G}_{n-1}^{\dagger}\mathcal{H}\hat{G}_{n})$ is used instead.

##### Energies from QLR.–

Once the Lanczos coefficients are obtained, the Hamiltonain matrix in its
tridiagonal form is then known. Diagonalizing this matrix and retaining both
the energies and eigenvectors makes the new ground state in terms of the
Krylov basis chosen.

Defining an operator $\hat{Y}^{(g)}$ for an excitation $g\in\mathbb{Z}^{+}$,
the operator for the $g$th excitation would be defined as

$\hat{Y}^{(g)}=\sum_{n}\gamma^{(g)}_{n}\hat{G}_{n}$ (6)

in terms of coefficients $\gamma_{n}$ found from the diagonalization. The
summation over $n$ is over as many elements in the Krylov basis as kept.

Thus, the QLR algorithm can access ground states ($g=0$) or excited states
($g\geq 1$) by simply selecting different coefficients from the
diagonalization of $\mathcal{H}$. At the end of this procedure, the ground
state is recovered and can be used without fully re-preparing it.

###  II.1 Reduced computational complexity of applying operators

##### First strategy: apply interaction incrementally.–

When applying an operator to a wavefunction, several auxiliary qubits are
used. OAA is a Grover search over auxuiliary qubits, so it will scale
exponentially, $\mathcal{O}(\sqrt{2^{q}})$ for $q$ qubits, in terms of the
number of auxiliary qubits. While Grover’s algorithm promises a square root
speedup over classical search methods, the size of the database on the quantum
computer is exponentially large in terms of the number of qubits. This calls
into question whether operators can actually be applied to wavefunctions in
the manner prescribed earlier in total, but the particular details of the
algorithm will show that this can be done efficiently. This clearly will not
present an issue for finding quantities necessary for density functional
theory Baker and Poulin (2020), since the required operators are small ($q=2$
for coefficients of the one-body reduced density matrix).

By applying operators from the full Hamiltonian incrementally in groups of
size $D$, the cost can be reduced to $\mathcal{O}(D2^{\frac{\lceil
q/D\rceil}{2}})$. Thus, the OAA is feasible here since the Hamiltonian can be
expressed as a linear combination of unitaries, this strategy applies to any
system of interest Berry  _et al._ (2015). The minimal size of a group of
operators is the same as a single term in the Hamiltonian. For many-body
problems, this will be $q=$ 2 or 4 terms for the quadratic single-body terms
or the quartic electron-interaction terms Baker and Poulin (2020).

##### Second strategy: decrease strength of interaction.–

Note that the operator applied onto the original wavefunction $\Psi$ does not
need to be an operator whose eigenvalues include $\Psi$. This means that
starting from some known state $\Psi$ with associated operator
$\mathcal{H}_{0}$,

$\mathcal{H}=\mathcal{H}_{0}+\lambda\mathcal{H}_{I}$ (7)

can be defined as a new Hamiltonian. This form of the Hamiltonian has a
similar form as the the adiabatic evolution in Eq. (1), except that the
$\lambda$ coefficient can be much larger than the step size in time evolution.

The interaction term, $\mathcal{H}_{I}$, can be the full interaction term or
one term of that interaction. By adding the terms incrementally, the final
state is closer to the initial state provided and will allow for the algorithm
to find the new ground state. It is crucial that each step of this procedure
begins with an eigenvector so that the QPE can be applied correctly in the
quantum counting procedure.

$\mathbf{J_{xy}=J_{z}}$

Figure 1:  Convergence in energy ($E$) while adding one interaction term
$\hat{S}^{z}_{i}\cdot\hat{S}^{z}_{i+1}$ to an XY model on a 10 site lattice (9
interaction terms to add). Two solvers are used at each step. DMRG is shown as
a solid blue line and a Lanczos solver is shown as a red dashed line. The
initial wavefunction for the Lanczos algorithm is taken from the previous
iteration. A black dashed line shows the energy of the full Heisenberg model.
An inset shows the difference between the DMRG and Lanczos solutions. Only one
Lanczos iteration was used for each added term. Inset: energy differences from
the exact value, $\Delta E$.

###  II.2 Demonstration on model systems

A full demonstration on real systems by implementing a linear combination of
unitaries will be delayed for a future work, but it will be shown here taht
Lanczos can be applied iteratively according to the suggestions in the
previous paragraphs.

In order to demonstrate that only a limited number of Lanczos steps can be
used to solve a model that includes a small number of extra terms from the
starting wavefunction’s Hamiltonian, a numerical study on a 10-site XXZ model
of the form

$\mathcal{H}=\sum_{i}J_{xy}\left(\hat{S}^{x}_{i}\hat{S}^{x}_{i+1}+\hat{S}^{y}_{i}\hat{S}^{y}_{i+1}\right)+J_{z}\hat{S}^{z}_{i}\hat{S}^{z}_{i+1}$
(8)

where $J_{xy}=J_{z}$ is the Heisenberg Hamiltonian and spin matrices
$\mathbf{S}=(\hat{S}^{x},\hat{S}^{y},\hat{S}^{z})=\frac{1}{2}\boldsymbol{\sigma}$
(with $\hbar=1$) for the vector of Pauli matrices $\boldsymbol{\sigma}$
Townsend (2000). To start, an XY model,

$\mathcal{H}_{\mathrm{XY}}=J_{xy}\sum_{i}\hat{S}^{x}_{i}\hat{S}^{x}_{i+1}+\hat{S}^{y}_{i}\hat{S}^{y}_{i+1},$
(9)

will be solved for the initial wavefunction and a single interaction term of
the form $\hat{S}^{z}_{i}\hat{S}^{z}_{i+1}$ will be added. In all, 9 terms
will be added to the Hamiltonian. Computations were made with the DMRjulia
library using the density matrix renormalization group (DMRG), Lanczos, and
exact diagonalization (ED) routines Baker  _et al._ (2021b); *baker2019m;
Baker and Thompson (2021); Baker .

There are three cases of study here: $J_{xy}=J_{z}$ (small perturbations),
$J_{xy}\ll J_{z}$ (large perturbations), and when the initial $\psi_{0}$ is
far from the starting state.

$\mathbf{J_{xy}\ll J_{z}}$

Figure 2:  Convergence in energy ($E$) while adding a large interaction term
$\hat{J}_{z}S^{z}_{i}\cdot\hat{S}^{z}_{i+1}$ with $J_{xy}\ll J_{z}$. The
parameters are the same as FIg. 1. By using 4 Lanczos steps instead of 2, the
accuracy is greatly improved in this case. Decreasing the incremental
interaction strength, $d\lambda$ (here, 0.5 with two rounds of
$N_{\mathrm{Lanczos}}=2$), would require more than double-precision (see
text).

####  II.2.1 Small perturbations ($J_{xy}=J_{z}$)

Figure 1 shows how the energy converges to the energy of the full Heisenberg
model with increasing numbers of terms. In each case, using the initial
wavefunction provided from the previous iteration, only one Lanczos iteration
must be run to obtain a highly accurate energy with the new interaction term
included. The small difference in energies between the DMRG solution and the
Lanczos solution is shown in the inset of Fig. 1. There is also the
possibility to add in a partial term by using a small term $d\lambda$ for each
of $N$ times such that $Nd\lambda=\lambda$. This would require more
applications of the Lanczos algorithm with new Hamiltonians and may be useful
for long-range interactions or other cases.

In this one-dimensional example, the convergence is aided by only adding one
term. In realistic systems, the same strategy of adding interaction terms
piecemeal will both allow the Lanczos algorithm to converge quickly and also
to reduce the time spent applying the operator with OAA. Still, this should be
expected in general due to the rapid convergence of Lanczos techniques.

####  II.2.2 Large perturbations ($J_{xy}\ll J_{z}$)

In this case, $100J_{xy}=J_{z}$ to simulate a large perturbation. For the
example here, increasing the number of Lanczos iterations to 4 for each new
interaction term added allows for the energy to be obtained to a good accuracy
here. This is a large improvement over using only 2 Lanczos iterations as
shown by comparing the curves in Fig. 2.

It is possible to use only 2 Lanczos updates if one adds portions of the
interaction term. In this case, $J_{z}/2$ is added to the Hamiltonian twice
with 2 Lanczos steps each (total 4 Lanczos steps). The results could be
improved if more Lanczos iterations were used here, but this would require
greater precision in the coefficients than is available here with a double-
precision classical implementation.

Note that tests with coefficients found on the classical computer will be
subject to numerical instabilities that will eventually degrade the accuracy
of the resulting ground state, so a full solution with this method will
generate imprecise answers. For example, when increasing the number of Lanczos
iterations to 10 for each step, this can occur in this example.

random $\psi_{0}$

Figure 3:  Convergence in energy ($E$) similar to Fig. 1 except that a random
starting state was used. Tests on larger systems or with more Lanczos
recursion steps quickly suffer from precision errors that will not be present
on the quantum computer.

####  II.2.3 Arbitrary starting wavefunction

So far, eigenstates for the starting wavefunction were used. The natural
question is whether the starting states could be replaced by some arbitrary
state and how the convergence of the algorithm is affected. When implementing
this type of initialization on the quantum computer, the only key is to
determine the energy of the wavefunction. If using QPE to do this, then the
starting state must have an associated Hamiltonian.

For demonstration purposes, the initial state of alternating spins here is

$|\Psi\rangle=|\psi_{0}\rangle=|\uparrow\uparrow\downarrow\downarrow\downarrow\uparrow\downarrow\downarrow\uparrow\uparrow\rangle$
(10)

is taken as the initial state with $J_{xy}=J_{z}$ in the Hamiltonian. While
the Hamiltonian which has an eigenstate of Eq. (10) is not known and would
therefore prevent the use of this state in the quantum algorithm, it is
instructive to observe the convergence in this case where the initial state is
not close to the final problem to solve.

Figure 3 shows the convergence for this starting wavefunction. The first
Lanczos iteration will be the most difficult since it must change the
wavefunction the most. Subsequent iterations converge more easily to the true
ground-state just as in the previous cases.

These examples demonstrate that applying operators incrementally is possible
and will yield a large cost reduction on the quantum computer. The examples
here were on the ground-state, but higher excited states can be found similar
to the discussion around Eq. (6) once the Lanczos coefficients are known. In
this example, the precision that can be found from the classical simulation is
restricted Cahill  _et al._ (2000). The quantum computer is only limited is
not a limitation on the quantum computer assuming perfect application of the
operators.

##  III Quantum block Lanczos recursion

In principle should be able to discover all excitations from QLR as described
above Chepiga and Mila (2017), but the extension of this method to a block or
banded Lanczos algorithm is known to resolve degeneracies to a higher degree
and aid convergence generally Baker  _et al._ (2021a). Block Lanczos is used
in quantum chemistry Bai  _et al._ (2000) but also in physics, particularly
for dynamical mean field theory computations Sénéchal (2008). The scalar
coefficients of the Lanczos recursion are extended to a matrix of
coefficients.

To formulate the problem, consider a set of wavefunctions grouped as a
supervectors $\mathbf{\Psi}$ is a vector of $d$ excitations,

$\mathbf{\Psi}_{n}=\Big{(}|\psi_{n}\rangle_{1},|\psi_{n}\rangle_{2},\ldots,|\psi_{n}\rangle_{d}\Big{)}$
(11)

and therefore the task is to find matrices $\mathbf{A}$ and $\mathbf{B}$ which
are of dimension $d\times d$ block such that Lanczos can be performed. This
means that $d$ registers, each with a wavefunction are also available.

Extending the 3-term recursion from Eq. 2 to a block or banded (or other)
Lanczos scheme with more terms in the recursion would be possible to prepare
more than one excitation at a time Sénéchal (2008). The expanded Lanczos
recursion relation appears as Weikert  _et al._ (1996)

$\mathbf{\Psi}_{n+1}\mathbf{B}_{n+1}=\mathcal{H}\mathbf{\Psi}_{n}-\mathbf{\Psi}_{n}\mathbf{A}_{n}-\mathbf{\Psi}_{n-1}\mathbf{B}_{n}^{\dagger}$
(12)

for a vector of wavefunctions $\mathbf{\Psi}$ with matrices

$\mathbf{A}_{n}=\mathbf{\Psi}^{\dagger}_{n}\mathcal{H}\mathbf{\Psi}_{n}$ (13)

and $\mathbf{B}_{n}$ defined recursively from

$\mathbf{B}_{n}=\mathbf{\Psi}_{n-1}^{\dagger}\mathcal{H}\mathbf{\Psi}_{n}=\mathbf{\Psi}_{n-1}^{\dagger}\mathcal{H}\mathbf{\Psi}_{n-1}=\mathbf{B}_{n}^{\dagger}$
(14)

as can be most immediately seen from the tensor network diagrams in Ref. Baker
_et al._ , 2021a. It should be understood that when applying $\mathcal{H}$
onto the wavefunction supervector that the Hamiltonian is applied onto each
wavefunction in the vector. The implication of Eqs. 13 and 14 is that there
are simply $d^{2}$ applications of QLR to obtain QBLR.

Because the coefficients can be obtained by an operator acting on a ground
state of an initial Hamiltonian, this implies that the coefficients of both
$\mathbf{A}$ and $\mathbf{B}$ can be recovered by performing a state-
preserving quantum counting operation on each term of the matrix. This was
also true in the scalar case, and this establishes that the block Lanczos case
is in the RWMP category of algorithms thus avoiding repeated wavefunction
preparation Baker and Poulin (2020). Each coefficient must be obtained
individually: one quantum counting process for each operator. However, If
operations can be performed in parallel on the quantum computer, then this
could improve performance to find all coefficients of a block simultaneously.

On the classical computer, the blocks can be represented as a block diagonal
super-matrix representing the Hamiltonian

$\boldsymbol{\mathcal{H}}=\left(\begin{array}[]{cccccc}\mathbf{A}_{0}&\mathbf{B}_{1}^{\dagger}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{B}_{1}&\mathbf{A}_{1}&\mathbf{B}_{2}^{\dagger}&\cdots&\mathbf{0}\\\
\mathbf{0}&\mathbf{B}_{2}&\mathbf{A}_{2}&\ddots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\ddots&\mathbf{B}_{n}^{\dagger}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{B}_{n}&\mathbf{A}_{n}\\\
\end{array}\right)$ (15)

which can be diagonalized to find the energies of the excitations. From the
coefficients of the eigenvectors of $\boldsymbol{\mathcal{H}}$ (denoted as
$\gamma_{n}$ still here too), the excitation wavefunctions can also be
obtained.

In order to understand how best to apply the operators on the quantum
computer, the block Lanczos equation can be written as

$\displaystyle\sum_{j}B^{(n+1)}_{ij}|\psi_{n+1}\rangle_{j}=$
$\displaystyle\mathcal{H}|\psi_{n}\rangle_{i}-\sum_{j}A_{ij}^{(n)}|\psi_{n}\rangle_{j}$
(16) $\displaystyle-\sum_{j}(B^{*})_{ij}^{(n)}|\psi_{n-1}\rangle_{j}$

for a single row $i$ of the left hand side. This form makes clear that the $d$
registers containing the $d$ excitations can be acted on with the appropriate
operators (sum over $j$) and added with the wavefunction on other registers
with the appropriate operators. Note that the elements of $\mathbf{A}$ and
$\mathbf{B}$ are stored classically and therefore the inverse of each matrix
can be found. If this is done, then apply $\mathbf{B}^{-1}_{n+1}$ once the
coefficients are found from Eq. (14).

The operators from Eq. (6) will be denoted as $\mathbf{G}_{n}$ for a given
level and are extended from the scalar Lanczos case as
$\hat{G}_{n}\rightarrow\mathbf{G}_{n}$. Similarly,
$\psi\rightarrow\mathbf{\Psi}$ for the wavefunction. One additional operation
is included $\mathbf{G}$, compared with the operators listed in Ref. Baker,
2021, that of the inverse $\mathbf{B}^{-1}$ operator. Since the operators can
be determined on the classical computer, the operators can simply be prepared
in a slightly different way. Again, the inverse is converted to a unitary as
is standard for applying operators Baker (2021); Baker and Poulin (2020);
Berry  _et al._ (2015); Kothari (2014); Guerreschi (2019).

In all, the algorithm costs a practical amount of $d^{2}$ over the scalar
Lanczos case. Each element of the matrix equations is simply one application
of a set of equations slightly modified from QLR. Note that the effort
expended to apply operators onto wavefunctions is not wasted at the end with a
single measurement. The state-preserving quantum counting algorithm can be
used to obtain the correct expectation values without completely measuring the
wavefunction.

###  III.1 Algorithm summary

The following summarizes the previous discussion to demonstrate how to
implement block Lanczos on the quantum computer.

Quantum block Lanczos recursion: Algorithm for multiple excited-states

  1. 1.

A set of starting eigenfunctions $\\{|\Xi\rangle\\}$ is prepared on $\nu$
registers (one for each wavefunction).

  2. 2.

Start a counter at $n=0$.

  3. 3.

Construct the operator $\mathbf{G}_{n}$ representing the appropriate step from
Eq. (12). This operator spans several excitations.

  4. 4.

$\mathbf{G}_{n}$ is applied onto the current state and quantum counting is
used to determine the Lanczos coefficients, given by Eq. (13) and (13).

  5. 5.

With the new coefficient(s), the algorithm returns to step 2, increments $n$,
and finds the next step of coefficients.

  6. 6.

Storing the coefficients classically, the block diagonal Hamiltonian from Eq.
(15) can be formed and diagonalized. A set of coefficients $\gamma_{n}$ can
then be used to determine the energies.

The scaling of the method is no worse than the application of the operators at
each step, albeit the number of times this must be run for $d$ excitations is
$\mathcal{O}(d^{2})$. As pointed out in Sec. II.1, the incautious application
of OAA will scale exponentially, but to aid convergence and avoid large
computations of operations the interaction term can be applied in smaller
pieces. This drastically reduces the number of Lanzcos steps required to
obtain the next wavefunction and limits the time necessary to run the OAA
algorithm. It is also possible to use a repeat-until-success strategy here
Baker (2021); Guerreschi (2019); Paetznick and Svore (2013), but the
guaranteed result from application of OAA probably will save time in the long
run.

It is not known beforehand the coefficients of the next step, so each step
must be performed iteratively. Once the coefficients are obtained with the
quantum counting process, then they can be stored classically and used without
obtaining them again, although the $d^{2}$ coefficients of a given block can
be parallelized. Note also that the algorithm can be restarted at any time
since the coefficients are stored classically.

The ground state energy can be checked when diagonalizing Eq. (15) on the
classical computer. If the ground state energy is converged within the
accuracy needed for QPE, then there are sufficient Lanczos steps run to find
an accurate ground state energy.

###  III.2 Error analysis

In this section, a careful analysis of the relevant parameters of the system
will show how uncertainty in the Lanczos coefficients will influence the
resulting energies.

Figure 4: Mean average error for a system of block size $d=20$ for various
numbers of iterations $\\{4,5,6,\ldots,20\\}$. All lines are nearly the same
in magnitude. The general linear trend is mostly unchanging with the input
parameters.

To investigate the effects of noise on the resulting energies, a block
Hamiltonian of a particular size $a$ is created for $b$ such blocks. When
$a=1$, the QLR algorithm is used. When $a>1$, QBLR is represented. A set of
matrices are generated with values between $[0,1]$ to mimic the values on the
operators for the quantum computation. A noise term selected from a Gaussian
distribution of width $\eta$ is selected for each element.

The energies with noise $E$ are computed with respect to the true energies
without noise $E_{0}$ to find the mean absolute error (MAE) $\bar{E}$ defined
as

$\bar{E}=\frac{1}{M}\sum_{i=1}^{M}\left|E^{(i)}-E^{(i)}_{0}\right|$ (17)

for a number of eigenvalues $M$.

The MAE is shown in Fig. 4. As expected, as the noise decreases in each of the
parameters, the accuracy of the energy eigenstates increases. In both cases,
the error decreases linearly with the noise applied to each term. This is a
straightforward but useful result. If the quantum counting obtains
coefficients to a precision $\delta$, then the resulting energy is also
obtained to an error of roughly $\delta$. The change in this trend for other
parameters that could have been picked (larger block sizes, more blocks, etc.)
is very stable.

What this implies is that for a QPE to obtain the correct eigenvalue, the
Lanczos coefficients must be sufficiently accurate (i.e., $\delta$ is a
threshold for how accurate the coefficients will be). This is a crucial aspect
to know how much precision must be obtained from the quantum counting
algorithm.

###  III.3 Discussion

The advantages of many excitations being discovered is chiefly one of
stability. The ability to resolve degeneracies and generate orthogonal
wavefunctions that is valuable in a variety of contexts and now available for
the quantum computer. That this method is rapidly convergent Ye (1991); Cullum
(1996), resolves degeneracies, and avoids repeated wavefunction preparation
means that this method could be a valuable alternative to time evolution
methods.

Note that unlike in Ref. Baker  _et al._ , 2021a, a singular value, QR, or LQ
decomposition Press  _et al._ (1992) was not used. This is primarily because
it was not necessary here, but it has several advantages. First, many of these
algorithms on the quantum computer require a quantum memory, which may be far
off in terms of development Chia  _et al._ (2020); Tang (2018), particularly
with the growing consensus that probably approximately correct (PAC) learning
on the quantum computer is not likely to be the best way forward for quantum
machine learning Baker and Poulin (2020); Huang  _et al._ (2021a). Although,
expanding the definition has been suggested to lead to a logarithmic
dependence on the number of measurements for tomography in terms of the total
Hilbert space size Aaronson (2019); Huang  _et al._ (2021b), finding the
initial state for the tomography could be accomplished with QLR or QBLR.

The goal of QBLR is to retain the accuracy of the wavefunction, which is
controlled directly through the gate fidelities and other quantities on the
quantum computer. The quantum advantage sought here is for an exponential
reduction in memory, which is a property that quantum computers exhibit under
the extended Church-Turing thesis which asks if an algorithm can be computed
with significantly less memory resources over classical computing Aaronson and
Arkhipov (2011). Lanczos techniques do have this property when applied on a
quantum computer theoretically if perfect application of gates is assumed and
with sufficient time to sample the coefficients. The typical errors on a
classical computer Cahill  _et al._ (2000); Paige (1976); Bai (1994); Chepiga
and Mila (2017); Aaronson and Arkhipov (2011) do not appear since values are
retained to quantum precision.

The opportunity to obtain ground states from QLR places added emphasis on the
development of state-preserving quantum counting and methods of applying
operators to wavefunctions such as the linear combination of unitaries (and
included subroutines such as OAA) to obtain both the Green’s function and
ground-state.

The possibility to discover ground states or excited states deserve more
consideration since these methods are free of time steps and Trotter-Suzuki
decompositions, a common feature in many other proposals. The Lanczos
algorithms also do not need any notion of locality to aid convergence as would
be expected for a matrix product state’s relationship to local Hamiltonians
Baker and Thompson (2021); Verstraete and Cirac (2006); Hastings (2004),
particularly for time-evolving block decimation (TEBD) Baker  _et al._
(2021b). The only instance where locality plays any role is in the writing of
the application of the operators to the wavefunction, which could play a role
in near-term implementation. It is not clear how many Lanczos coefficients
would need to be obtained to find accurate ground states from this method for
phase-estimation to give accurate results. Further discussion would require
specialization to a specific problem and is deferred to a future study.

Note that extensions to a fitting function for the remaining coefficients is
possible in principle Viswanath and Müller (2008).

##  IV Non-Hermitian operators

The previous discussion was for the block Lanczos algorithm that is formulated
for the problem of a Hermitian operator. The extension of the above methods to
non-Hermitian operators is possible without introducing too much additional
computational cost. Instead, only one additional set of matrices (or in the
case of scalar Lanczos, one additional set of coefficients) should be found.
The procedure will be define for the block-matrix case, but it can be reduced
to the scalar case when the block size is of size 1.

Given a non-Hermitian operator $\tilde{\mathcal{H}}$ with eigenvalue $E$, a
set of two eigenvalues can be defined. One is known as the left-eigenvectors
(a transpose is used even if there are complex valued entries) Bai (1994)

$\mathbf{\Psi}_{L}^{T}\tilde{\mathcal{H}}=\mathbf{\Psi}_{L}^{T}E$ (18)

and the other is for the right-eigenvectors,

$\tilde{\mathcal{H}}\mathbf{\Psi}_{R}=E\mathbf{\Psi}_{R}$ (19)

where block matrices $\mathbf{\Psi}_{L}$ and $\mathbf{\Psi}_{R}$ have the
orthogonality relation Bai (1994)

$\mathbf{\Psi}_{L}^{T}\mathbf{\Psi}_{R}=\mathbb{I}_{d\times d}$ (20)

with an identity matrix $\mathbb{I}_{d\times d}$ of dimension $d\times d$.

The block operator represented analogous to Eq. (15) appears in this case as

$\mathbf{\check{T}}_{n}=\left(\begin{array}[]{cccccc}\mathbf{A}_{0}&\mathbf{C}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{B}_{1}&\mathbf{A}_{1}&\mathbf{C}_{2}&\cdots&\mathbf{0}\\\
\mathbf{0}&\mathbf{B}_{2}&\mathbf{A}_{2}&\ddots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\ddots&\mathbf{C}_{n}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{B}_{n}&\mathbf{A}_{n}\\\
\end{array}\right).$ (21)

There are two recursion relations for each of the left and right eigenvectors.
They take the form of Bai  _et al._ (1999); Grüning  _et al._ (2011)

$\displaystyle\vec{\mathbf{\Psi}}_{L,(n)}^{T}\tilde{\mathcal{H}}$
$\displaystyle=\mathbf{\check{T}}_{n}\vec{\mathbf{\Psi}}_{L,(n)}^{T}+\mathcal{I}_{n}\mathbf{C}_{n+1}\mathbf{\Psi}^{T}_{L,(n+1)}$
(22) $\displaystyle\tilde{\mathcal{H}}\vec{\mathbf{\Psi}}_{R}$
$\displaystyle=\vec{\mathbf{\Psi}}_{R}\mathbf{\check{T}}_{n}+\mathbf{\Psi}_{R,(n+1)}\mathbf{B}_{n+1}\mathcal{I}_{n}^{T}$
(23)

where

$\displaystyle\vec{\mathbf{\Psi}}_{R,(n)}$
$\displaystyle=\left(\begin{array}[]{ccccc}\mathbf{\Psi}_{R}^{(0)},&\mathbf{\Psi}_{R}^{(1)},&\ldots,&\mathbf{\Psi}_{R}^{(n)}\end{array}\right),$
(25) $\displaystyle\vec{\mathbf{\Psi}}_{L,(n)}$
$\displaystyle=\left(\begin{array}[]{ccccc}\mathbf{\Psi}_{L}^{(0)},&\mathbf{\Psi}_{L}^{(1)},&\ldots,&\mathbf{\Psi}_{L}^{(n)}\end{array}\right),$
(27) $\displaystyle\mathrm{and}\quad\mathcal{I}_{n}$
$\displaystyle=\left(\begin{array}[]{cccc}\mathbf{0}_{d\times
d},&\mathbf{0}_{d\times d},&\ldots,&\mathbb{I}_{d\times d}\end{array}\right).$
(29)

The final term is an identity matrix to ensure that the dimension of the
resulting matrix is consistent with all terms in the recursion. Only the last
block of the matrix is non-zero where blocks of zero matrices,
$\mathbf{0}_{d\times d}$.

Writing out the explicit terms of the recursion relation in the style of the
$\mathbf{G}$ operators from earlier is lengthy, but a single operator can be
written to be applied on the starting wavefunction just as for QLR Baker
(2021). The only additional cost is the determination of an extra block
matrix, $\mathbf{C}$, and double the number of registers to represent both the
left- and right-eigenvectors.

The extension of the block Lanczos algorithm to the non-Hermitian case has
therefore only introduced an extra matrix $\mathbf{C}_{n}$ which must also be
found with the other two sets of matrices. However, this does not increase the
overall computational cost of the algorithm presented here. The same forms for
the operators connecting the original wavefunction provided to the algorithm
can be derived similarly as for the Hermitian case.

##  V Conclusion

Lanczos recursion methods on the quantum computer were extended to solve for
many excitations simultaneously. This uses a block Lanczos technique that is
good at resolving degeneracies in quantum states. This comes at only a cost of
the number of excitations squared sought on the quantum computer. The use of
the quantum counting algorithm here allows for the wavefunction to not be
collapsed at each step, cutting out a major cost of many other algorithms.
Further, interaction terms can be applied in small groups to aid convergence
and keep the process of applying those operators to the wavefunction less than
exponentially long. The error of diagonalizing the Hamiltonian in the Krylov
basis was demonstrated to scale linearly with the noise of the coefficients.
The method can also be applied onto non-Hermitian operators with a moderate
additional cost. Due to the feasible cost of this algorithm and the rapid
convergence of Lanczos techniques, it is expected that quantum block Lanczos
recursion could be an alternative to existing methods to find the ground and
excited states.

##  VI Acknowledgements

The author thanks Alexandre Foley, David Sénéchal, Anirban Chowdhury, and
David Poulin for enlightening conversations and to Rex Godby for hosting.

The author is grateful to the US-UK Fulbright Commission for financial support
under the Fulbright U.S. Scholarship programme as hosted by the University of
York. This research was undertaken in part thanks to funding from the Bureau
of Education and Cultural Affairs from the United States Department of State.

This project was undertaken on the Viking Cluster, which is a high performance
compute facility provided by the University of York. T.E.B. is grateful for
computational support from the University of York High Performance Computing
service, Viking and the Research Computing team.

## References

  * Nielsen and Chuang (2010) Michael A Nielsen and Isaac L Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2010).
  * O’Malley  _et al._ (2016) Peter JJ O’Malley, Ryan Babbush, Ian D Kivlichan, Jonathan Romero, Jarrod R McClean, Rami Barends, Julian Kelly, Pedram Roushan, Andrew Tranter, Nan Ding, _et al._ , “Scalable quantum simulation of molecular energies,” Physical Review X 6, 031007 (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).
  * 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).
  * Gokhale  _et al._ (2020) Pranav Gokhale, Olivia Angiuli, Yongshan Ding, Kaiwen Gui, Teague Tomesh, Martin Suchara, Margaret Martonosi, and Frederic T Chong, “$o(n^{3})$ measurement cost for variational quantum eigensolver on molecular hamiltonians,” IEEE Transactions on Quantum Engineering 1, 1–24 (2020).
  * Ollitrault  _et al._ (2020) Pauline J Ollitrault, Abhinav Kandala, Chun-Fu Chen, Panagiotis Kl Barkoutsos, Antonio Mezzacapo, Marco Pistoia, Sarah Sheldon, Stefan Woerner, Jay M Gambetta, and Ivano Tavernelli, “Quantum equation of motion for computing molecular excitation energies on a noisy quantum processor,” Physical Review Research 2, 043140 (2020).
  * 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).
  * Kremenetski  _et al._ (2021) Vladimir Kremenetski, Carlos Mejuto-Zaera, Stephen J Cotton, and Norm M Tubman, “Simulation of adiabatic quantum computing for molecular ground states,” arXiv preprint arXiv:2103.12059 (2021).
  * Bittel and Kliesch (2021) Lennart Bittel and Martin Kliesch, “Training variational quantum algorithms is np-hard,” Physical Review Letters 127, 120502 (2021).
  * Poulin  _et al._ (2015) David Poulin, Matthew B Hastings, Dave Wecker, Nathan Wiebe, Andrew C Doherty, and Matthias Troyer, “The Trotter step size required for accurate quantum simulation of quantum chemistry,” Quantum Information and Computation 15, 0361–0384 (2015).
  * Lemieux  _et al._ (2021) Jessica Lemieux, Guillaume Duclos-Cianci, David Sénéchal, and David Poulin, “Resource estimate for quantum many-body ground-state preparation on a quantum computer,” Physical Review A 103, 052408 (2021).
  * Baker and Poulin (2020) Thomas E Baker and David Poulin, “Density functionals and Kohn-Sham potentials with minimal wavefunction preparations on a quantum computer,” Physical Review Research 2, 043238 (2020).
  * Gulania and Whitfield (2021) Sahil Gulania and James Daniel Whitfield, “Limitations of hartree–fock with quantum resources,” The Journal of Chemical Physics 154, 044112 (2021).
  * Schuch and Verstraete (2009) Norbert Schuch and Frank Verstraete, “Computational complexity of interacting electrons and fundamental limitations of density functional theory,” Nature Physics 5, 732–735 (2009).
  * Higgott  _et al._ (2019) Oscar Higgott, Daochen Wang, and Stephen Brierley, “Variational quantum computation of excited states,” Quantum 3, 156 (2019).
  * Nakanishi  _et al._ (2019) Ken M Nakanishi, Kosuke Mitarai, and Keisuke Fujii, “Subspace-search variational quantum eigensolver for excited states,” Physical Review Research 1, 033062 (2019).
  * Parrish  _et al._ (2019) Robert M Parrish, Edward G Hohenstein, Peter L McMahon, and Todd J Martínez, “Quantum computation of electronic transitions using a variational quantum eigensolver,” Physical review letters 122, 230401 (2019).
  * Mizuta  _et al._ (2021) Kaoru Mizuta, Mikiya Fujii, Shigeki Fujii, Kazuhide Ichikawa, Yutaka Imamura, Yukihiro Okuno, and Yuya O Nakagawa, “Deep variational quantum eigensolver for excited states and its application to quantum chemistry calculation of periodic materials,” arXiv preprint arXiv:2104.00855 (2021).
  * Kuroiwa and Nakagawa (2021) Kohdai Kuroiwa and Yuya O Nakagawa, “Penalty methods for a variational quantum eigensolver,” Physical Review Research 3, 013197 (2021).
  * Greene-Diniz and Muñoz Ramo (2021) Gabriel Greene-Diniz and David Muñoz Ramo, “Generalized unitary coupled cluster excitations for multireference molecular states optimized by the variational quantum eigensolver,” International Journal of Quantum Chemistry 121, e26352 (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).
  * Jahangiri  _et al._ (2020) Soran Jahangiri, Juan Miguel Arrazola, Nicolás Quesada, and Alain Delgado, “Quantum algorithm for simulating molecular vibrational excitations,” Physical Chemistry Chemical Physics 22, 25528–25537 (2020).
  * Roggero  _et al._ (2020) Alessandro Roggero, Chenyi Gu, Alessandro Baroni, and Thomas Papenbrock, “Preparation of excited states for nuclear dynamics on a quantum computer,” Physical Review C 102, 064624 (2020).
  * 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).
  * Bauman  _et al._ (2019) Nicholas P Bauman, Guang Hao Low, and Karol Kowalski, “Quantum simulations of excited states with active-space downfolded hamiltonians,” The Journal of chemical physics 151, 234114 (2019).
  * Elliott  _et al._ (2011) Peter Elliott, Sharma Goldson, Chris Canahui, and Neepa T Maitra, “Perspectives on double-excitations in tddft,” Chemical Physics 391, 110–119 (2011).
  * Khoromskaia and Khoromskij (2015) Venera Khoromskaia and Boris N Khoromskij, “Tensor numerical methods in quantum chemistry: from hartree–fock to excitation energies,” Physical Chemistry Chemical Physics 17, 31491–31509 (2015).
  * Helgaker  _et al._ (2014) Trygve Helgaker, Poul Jorgensen, and Jeppe Olsen,  _Molecular electronic-structure theory_ (John Wiley & Sons, 2014).
  * Sun  _et al._ (2021) Shi-Ning Sun, Mario Motta, Ruslan N Tazhigulov, Adrian TK Tan, Garnet Kin-Lic Chan, and Austin J Minnich, “Quantum computation of finite-temperature static and dynamical properties of spin systems using quantum imaginary time evolution,” PRX Quantum 2, 010317 (2021).
  * Baker (2021) Thomas E Baker, “Lanczos recursion on a quantum computer for the Green’s function and ground state,” Physical Review A 103, 032404 (2021).
  * Kothari (2014) Robin Kothari,  _Efficient algorithms in quantum query complexity_ , Ph.D. thesis, University of Waterloo (2014).
  * Berry  _et al._ (2015) Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma, “Simulating hamiltonian dynamics with a truncated taylor series,” Phys. Rev. Lett. 114, 090502 (2015).
  * Childs  _et al._ (2017) Andrew M Childs, Robin Kothari, and Rolando D Somma, “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision,” SIAM Journal on Computing 46, 1920–1950 (2017).
  * Grover (1996) Lov K Grover, “A fast quantum mechanical algorithm for database search,” in  _Proceedings of the twenty-eighth annual ACM symposium on Theory of computing_ (1996) pp. 212–219.
  * Grover (2001) Lov K Grover, “From schrödinger’s equation to the quantum search algorithm,” Pramana 56, 333–348 (2001).
  * Brassard  _et al._ (2002) Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp, “Quantum amplitude amplification and estimation,” Contemporary Mathematics 305, 53–74 (2002).
  * Cahill  _et al._ (2000) Eamonn Cahill, Alan Irving, Christopher Johnston, James Sexton, Ukqcd Collaboration, _et al._ , “Numerical stability of lanczos methods,” Nuclear Physics B-Proceedings Supplements 83, 825–827 (2000).
  * Jamet  _et al._ (2021) Francois Jamet, Abhishek Agarwal, Carla Lupo, Dan E Browne, Cedric Weber, and Ivan Rungger, “Krylov variational quantum algorithm for first principles materials simulations,” arXiv preprint arXiv:2105.13298 (2021).
  * Baker  _et al._ (2021a) Thomas E Baker, Alexandre Foley, and David Sénéchal, “Direct solution of multiple excitations in a matrix product state with block Lanczos,” arXiv preprint arXiv:2109.08181 (2021a).
  * Di Paolo  _et al._ (2021) Agustin Di Paolo, Thomas E Baker, Alexandre Foley, David Sénéchal, and Alexandre Blais, “Efficient modeling of superconducting quantum circuits with tensor networks,” npj Quantum Information 7, 1–11 (2021).
  * Marriott and Watrous (2005) Chris Marriott and John Watrous, “Quantum Arthur–Merlin games,” Comp. Complexity 14, 122–152 (2005).
  * Temme  _et al._ (2011) Kristan Temme, Tobias J Osborne, Karl G Vollbrecht, David Poulin, and Frank Verstraete, “Quantum metropolis sampling,” Nature 471, 87 (2011).
  * Guerreschi (2019) Gian Giacomo Guerreschi, “Repeat-until-success circuits with fixed-point oblivious amplitude amplification,” Physical Review A 99, 022306 (2019).
  * Paetznick and Svore (2013) Adam Paetznick and Krysta M Svore, “Repeat-until-success: Non-deterministic decomposition of single-qubit unitaries,” arXiv preprint arXiv:1311.1074 (2013).
  * 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).
  * Townsend (2000) John S Townsend,  _A modern approach to quantum mechanics_ (University Science Books, 2000).
  * Baker  _et al._ (2021b) Thomas E Baker, Samuel Desrosiers, Maxime Tremblay, and Martin P Thompson, “Méthodes de calcul avec réseaux de tenseurs en physique,” Canadian Journal of Physics 99, 4 (2021b).
  * Baker  _et al._ (2019) Thomas E Baker, Samuel Desrosiers, Maxime Tremblay, and Martin P Thompson, “Basic tensor network computations in physics,” arXiv preprint arXiv:1911.11566, p. 19 (2019).
  * Baker and Thompson (2021) Thomas E Baker and Martin P Thompson, “DMRjulia I. Basic construction of a tensor network library for the density matrix renormalization group,” arXiv preprint arXiv:2109.03120 (2021).
  * (51) Thomas E. Baker, “DMRjulia,” https://github.com/bakerte/DMRjulia.jl.
  * Chepiga and Mila (2017) Natalia Chepiga and Frédéric Mila, “Excitation spectrum and density matrix renormalization group iterations,” Physical Review B 96, 054425 (2017).
  * Bai  _et al._ (2000) Zhaojun Bai, James Demmel, Jack Dongarra, Axel Ruhe, and Henk van der Vorst,  _Templates for the solution of algebraic eigenvalue problems: a practical guide_ (SIAM, 2000).
  * Sénéchal (2008) David Sénéchal, “An introduction to quantum cluster methods,” arXiv preprint arXiv:0806.2690 (2008).
  * Weikert  _et al._ (1996) H-G Weikert, H-D Meyer, LS Cederbaum, and F Tarantelli, “Block lanczos and many-body theory: Application to the one-particle green’s function,” The Journal of chemical physics 104, 7122–7138 (1996).
  * Ye (1991) Qiang Ye, “A convergence analysis for nonsymmetric lanczos algorithms,” Mathematics of Computation 56, 677–691 (1991).
  * Cullum (1996) Jane Cullum, “Arnoldi versus nonsymmetric lanczos algorithms for solving matrix eigenvalue problems,” BIT Numerical Mathematics 36, 470–493 (1996).
  * Press  _et al._ (1992) William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery,  _Numerical recipes in C++_ (1992).
  * Chia  _et al._ (2020) Nai-Hui Chia, András Gilyén, Tongyang Li, Han-Hsuan Lin, Ewin Tang, and Chunhao Wang, “Sampling-based sublinear low-rank matrix arithmetic framework for dequantizing quantum machine learning,” in  _Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing_ (2020) pp. 387–400.
  * Tang (2018) Ewin Tang, “Quantum-inspired classical algorithms for principal component analysis and supervised clustering,” arXiv preprint arXiv:1811.00414 (2018).
  * Huang  _et al._ (2021a) Hsin-Yuan Huang, Richard Kueng, and John Preskill, “Information-theoretic bounds on quantum advantage in machine learning,” Physical Review Letters 126, 190505 (2021a).
  * Aaronson (2019) Scott Aaronson, “Shadow tomography of quantum states,” SIAM Journal on Computing 49, STOC18–368 (2019).
  * Huang  _et al._ (2021b) Hsin-Yuan Huang, Richard Kueng, Giacomo Torlai, Victor V Albert, and John Preskill, “Provably efficient machine learning for quantum many-body problems,” arXiv preprint arXiv:2106.12627 (2021b).
  * Aaronson and Arkhipov (2011) Scott Aaronson and Alex Arkhipov, “The computational complexity of linear optics,” in  _Proceedings of the forty-third annual ACM symposium on Theory of computing_ (2011) pp. 333–342.
  * Paige (1976) Christopher C Paige, “Error analysis of the lanczos algorithm for tridiagonalizing a symmetric matrix,” IMA Journal of Applied Mathematics 18, 341–349 (1976).
  * Bai (1994) Zhaojun Bai, “Error analysis of the lanczos algorithm for the nonsymmetric eigenvalue problem,” Mathematics of Computation 62, 209–226 (1994).
  * Verstraete and Cirac (2006) Frank Verstraete and J Ignacio Cirac, “Matrix product states represent ground states faithfully,” Physical Review B 73, 094423 (2006).
  * Hastings (2004) Matthew B Hastings, “Locality in quantum and markov dynamics on lattices and networks,” Physical review letters 93, 140402 (2004).
  * Viswanath and Müller (2008) VS Viswanath and Gerhard Müller,  _The Recursion Method: Application to Many-Body Dynamics_ , Vol. 23 (Springer Science & Business Media, 2008).
  * Bai  _et al._ (1999) Zhaojun Bai, David Day, and Qiang Ye, “Able: an adaptive block lanczos method for non-hermitian eigenvalue problems,” SIAM Journal on Matrix Analysis and Applications 20, 1060–1082 (1999).
  * Grüning  _et al._ (2011) Myrta Grüning, Andrea Marini, and Xavier Gonze, “Implementation and testing of lanczos-based algorithms for random-phase approximation eigenproblems,” Computational materials science 50, 2148–2156 (2011).