# On analog quantum algorithms for the mixing of Markov chains

Shantanav Chakraborty1<EMAIL_ADDRESS>Kyle Luh2<EMAIL_ADDRESS>Jérémie Roland1<EMAIL_ADDRESS>1QuIC, Ecole Polytechnique de Bruxelles,
Université libre de Bruxelles  2Center for Mathematical Sciences and
Applications, Harvard University

###### Abstract

The problem of sampling from the stationary distribution of a Markov chain
finds widespread applications in a variety of fields. The time required for a
Markov chain to converge to its stationary distribution is known as the
classical mixing time. In this article, we deal with analog quantum algorithms
for mixing. First, we provide an analog quantum algorithm that given a Markov
chain, allows us to sample from its stationary distribution in a time that
scales as the sum of the square root of the classical mixing time and the
square root of the classical hitting time. Our algorithm makes use of the
framework of interpolated quantum walks and relies on Hamiltonian evolution in
conjunction with von Neumann measurements.

There also exists a different notion for quantum mixing: the problem of
sampling from the limiting distribution of quantum walks, defined in a time-
averaged sense. In this scenario, the quantum mixing time is defined as the
time required to sample from a distribution that is close to this limiting
distribution. Recently we provided an upper bound on the quantum mixing time
for Erdös-Renyi random graphs [Phys. Rev. Lett. 124, 050501 (2020)]. Here, we
also extend and expand upon our findings therein. Namely, we provide an
intuitive understanding of the state-of-the-art random matrix theory tools
used to derive our results. In particular, for our analysis we require
information about macroscopic, mesoscopic and microscopic statistics of
eigenvalues of random matrices which we highlight here. Furthermore, we
provide numerical simulations that corroborate our analytical findings and
extend this notion of mixing from simple graphs to any ergodic, reversible,
Markov chain.

##  I Introduction

Markov chain-based algorithms are applied in a plethora of fields ranging from
statistical physics binder1993monte , combinatorial optimization
aarts1985statistical  to network science page1999pagerank  and form the basis
of Markov chain Monte Carlo-based methods gilks1995markov . In many of these
applications, the underlying task is often to sample from the so-called steady
state (also known as the stationary distribution) of the associated Markov
chain.

One way to sample from a stationary distribution is by mixing. The Markov
chain, which is represented by a stochastic matrix $P$ is applied repeatedly
to some initial distribution. The resultant random walk reaches a final
distribution that is close to a stationary distribution of $P$, irrespective
of the initial distribution. For most applications, the Markov chain is
ergodic, implying that it has a unique stationary distribution and,
reversible, i.e. it satisfies detailed balance. (We refer the reader to Sec.
II for details on the definitions of these terms related to Markov chains).
Henceforth, unless stated otherwise, we shall restrict our attention to
ergodic, reversible Markov chains. For a given Markov chain $P$, the minimum
time after which the distribution is $\epsilon$-close to the stationary
distribution is known as the mixing time of the random walk on $P$. It is well
known that the mixing time is related to the spectral gap of $P$. For an
ergodic, Markov chain with spectral gap $\Delta$, the mixing time is in
$\widetilde{\mathcal{O}}\left(1/\Delta\right)$ 111Throughout the article, we
use $\tilde{\mathcal{O}}(.)$ and $\tilde{\Omega}(.)$ to suppress
polylogarithmic factors, i.e.
$\tilde{\mathcal{O}}(f(n))=\mathcal{O}(f(n)\mathrm{polylog}(f(n)))$ and
$\tilde{\Omega}(f(n))=\Omega(f(n)\mathrm{polylog}(f(n)))$..

The stationary distribution, by definition, is the limiting distribution of
the resultant random walk on $P$, i.e. once the stationary state is reached,
the random walk ceases to evolve. This implies that as $t\rightarrow\infty$,
$P^{t}$ applied to any initial distribution converges to the stationary
distribution. Thus the classical mixing time is also the time required to
sample from the limiting distribution of the underlying random walk.

In the context of quantum algorithms, there arises two notions of mixing and
hence of mixing time. First, it is natural to consider whether, given a Markov
chain $P$, a quantum algorithm can allow us to prepare a coherent encoding of
the stationary distribution of $P$. We shall refer to this problem as
$QSSamp$. Measuring the output state of such an algorithm would enable us to
sample from the (classical) stationary state of $P$. Preparing such a coherent
encoding also has other applications which we discuss later.

The other notion of mixing arises from considering the limiting distribution
of the underlying quantum walk itself. As quantum evolutions are unitary and
hence distance-preserving, there is no inherent limiting stationary
distribution for quantum walks. However, it turns out that one can define a
limiting distribution of the quantum walk in a time-averaged sense.

Starting from some initial state, one can obtain the probability that the
walker is in some final state after a time $T$ which is picked uniformly at
random in the interval $[0,T]$. This gives a time-averaged probability
distribution at any time $T$ and also a limiting probability distribution as
$T\rightarrow\infty$. The problem of sampling from this time-averaged limiting
distribution of a quantum walk gives rise to another notion of mixing and we
shall refer to this problem as $QLSamp$. The mixing time of a quantum walk is
then defined as the time after which the time-averaged probability
distribution is close to the limiting probability distribution, i.e. the time
required to solve $QLSamp$.

In this article, we deal with both $QSSamp$ and $QLSamp$ problems. We provide
the first purely analog quantum algorithm to solve the $QSSamp$ problem while
for the $QLSamp$ problem, we expand and extend upon the results of Ref.
chakraborty2020fast , where we prove an upper bound for the quantum mixing
time for almost all graphs.

Aharonov and Ta-Shma aharonov2003adiabatic  demonstrated that the existence of
an efficient quantum algorithm for $QSSamp$ would imply that problems in the
complexity class Statistical Zero Knowledge (SZK) such as Graph Isomorphism
would be solvable in polynomial-time using a quantum computer (BQP), i.e.
$\mathrm{SZK}\subseteq\mathrm{BQP}$. This would be a surprising result as such
a generic $QSSamp$ algorithm would be oblivious to the specific structure of
the underlying problem. For example, consider the problem of Graph Isomorphism
aho1974design  (deciding whether two graphs are isomorphic to each other).
Given graphs $G_{1}$ and $G_{2}$, a quantum algorithm for mixing could be used
to prepare states that are a uniform superposition of all graphs that are
isomorphic to them. If these states are equal, then $G_{1}$ and $G_{2}$ are
isomorphic. A simple SWAP Test could then be used in conjunction with a
quantum algorithm for $QSSamp$ to solve Graph Isomorphism. Thus, generic
quantum algorithms for $QSSamp$ are unlikely to be efficient.

Having said that, there do exist quantum algorithms that solve this problem
richter2007quantum ; wocjan2008speedup ; dunjko2015quantum , some of which
have even been instrumental in obtaining speedups for quantum machine learning
paparo2014quantum ; orsucci2018optimal ; dunjko2018machine . Richter
richter2007almost  conjectured that one could construct a quantum algorithm
for this problem that has a running time that is in
$\widetilde{\mathcal{O}}(1/\sqrt{\Delta})$, yielding a quadratic speedup over
its classical counterpart. Developing quantum algorithms that match this
conjectured bound have been challenging. Most of the existing quantum
algorithms are based on Szegedy’s framework for discrete-time quantum walks
szegedy2004quantum .

The key idea that encompasses all existing algorithms for $QSSamp$ is to make
use of the so called quantum spatial search algorithm magniez2011search .
Given an ergodic, reversible Markov chain $P$ with a set of marked nodes, a
spatial search algorithm finds an element from this marked set. Classically,
this task requires a time known as the hitting time of the corresponding
random walk on $P$. It has been shown that a discrete-time quantum walk-based
quantum algorithm for spatial search can accomplish this task quadratically
faster (up to logarithmic factors) krovi2016quantum ; ambainis2019quadratic .
Such quantum algorithms start from the coherent encoding of the stationary
state of $P$ (it inherently assumes that this state can be prepared
efficiently) and end up in a state that has a constant overlap with an element
from the marked set. Thus, intuitively, quantum spatial search algorithms can
be run in reverse to obtain quantum mixing algorithms. However, simply
obtaining a constant overlap with the stationary state is not enough and these
mixing algorithms require the use of quantum phase estimation cleve1998quantum
and quantum amplitude amplification brassard2002quantum  to solve the $QSSamp$
problem. Recently, Apers and Sarlette provided a quantum algorithm that can
quadratically fast-forward the dynamics of Markov chains which can also be
used to solve the $QSSamp$ problem simon2018quantum . The running time of
these algorithms scale as the square root of the hitting time of the
corresponding quantum walk on the underlying Markov chain.

To the best of our knowledge, there do not exist any analog quantum algorithm
for solving the $QSSamp$ problem. In this framework, key algorithmic
primitives such as quantum phase estimation and quantum amplitude
amplification are missing as they are inherently discrete-time. In order to
construct an analog quantum algorithm for $QSSamp$ we assume that, given an
ergodic, reversible Markov chain $P$, we have access to a time-independent
Hamiltonian that encodes the connectivity of $P$. This Hamiltonian, defined in
Sec. IV, corresponds to a quantum walk on the edges of $P$. Furthermore it has
been recently used to design continuous-time quantum walk-based quantum
algorithms for spatial search that can find a single marked node on any
ergodic, reversible Markov chain in square root of the hitting time
chakraborty2018finding . We use the time-evolution of this Hamiltonian as the
key primitive to our algorithm. The second key primitive is to use von Neumann
measurements von1955mathematical  for quantum state generation. Childs et al.
used a sequence of such von Neumann measurements as an alternative to
adiabatic quantum computation and for solving combinatorial search algorithms
childs2002quantum . In Sec. III we demonstrate that this scheme can be used to
prepare eigenstates of Hamiltonians.

We show (Sec. V) that these two primitives allow us to develop a continuous-
time quantum walk based algorithm for spatial search. This algorithm differs
from the one developed in Ref. chakraborty2018finding  which makes use of
quantum phase randomization somma2010quantum . It provides an alternative
scheme by which one can find an element in a marked set of states of any
ergodic, reversible Markov chain in square root of the extended hitting time.
Although this algorithm has the same running time as that of Ref.
chakraborty2018finding , it provides useful intuition about how to build an
analog quantum algorithm for $QSSamp$.

Our quantum algorithm for mixing, explained in detail in Sec. VI, avoids the
need for amplitude amplification by making use of the framework of
interpolated Markov chains and switching between two different values of the
interpolation parameter. The running time scales as the sum of the square root
of the classical mixing time and the square root of the hitting time.

We also discuss the problem of $QLSamp$ on generic graphs. The limiting
distribution of quantum walks can be quite different from that obtained from a
quantum algorithm for solving $QSSamp$. Unlike its classical counterpart, for
$QLSamp$, the limiting distribution is dependent on the initial state of the
quantum walk. Moreover, instead of being dependent on the spectral gap
$\Delta$, the quantum mixing time depends on all eigenvalue gaps of the
underlying Hamiltonian. Aharonov et al. aharonov2001quantum  were the first to
study this problem. They showed that a discrete-time quantum walk on the cycle
graph mixes faster than its classical counterpart. Since then several works
have considered the mixing time of both continuous and discrete-time quantum
walks on specific graphs ahmadi2003mixing ; kendon2003decoherence ;
fedichkin2005mixing ; marquezino2008mixing ; marquezino2010mixing ;
kieferova2012quantum  . The upper bound for the mixing time of quantum walks
have been proven to be slower than its classical counterpart for some graphs
while a quadratic speedup has been obtained for others.

Recently, we proved an upper bound for the mixing time of quantum walks for
almost all graphs chakraborty2020fast , i.e. the fraction of graphs of $n$
nodes for which our upper bound holds, goes to 1 as $n$ goes to infinity. We
proved this by obtaining the mixing time for quantum walks on Erdös-Renyi
random graphs: graphs of $n$ nodes such that the probability of an edge
existing between any two nodes is $p$, typically denoted as $G(n,p)$. Here, we
expand upon the results of the letter chakraborty2020fast . In particular, our
goal is to offer an intuitive explanation of our proof techniques with an
emphasis on the several recently developed random matrix theory tools that
were used to derive the aforementioned results. We also corroborate our
analytical findings numerically and also extend the notion of $QLSamp$ to any
ergodic, reversible Markov chain. In fact, our numerical findings confirm the
fact that the mixing time for quantum walks on $G(n,p)$ is in
$\widetilde{\mathcal{O}}\left(n^{3/2}\right)$ for dense random graphs
(constant $p$). Additionally they also show that the limiting probability
distribution is close to the uniform distribution.

This article is organized as follows: In Sec. II, we explain some basic
concepts and quantities related to Markov chains that we shall use in
subsequent sections. In Sec. III we show how von Neumann measurements can be
used for preparing eigenstates of Hamiltonians. In Sec. IV, we define a
Hamiltonian corresponding to a quantum walk on the edges of any ergodic,
reversible Markov chain. In Sec. V, we make use of von Neumann measurements
and Hamiltonian evolution to provide a quantum algorithm for spatial search.
This provides an intuitive understanding of our analog quantum algorithm for
solving $QSSamp$, which we describe in Sec. VI. Next, in Sec. VII, we deal
with solving the $QLSamp$ problem. Finally, we conclude with a brief
discussion and summary in Sec. VIII.

##  II Preliminaries

In this section we state some basic definitions about Markov chains which we
shall use subsequently.

###  II.1 Basics of Markov chains

A Markov chain on a discrete state space $X$, such that $|X|=n$, can be
described by a $n\times n$ stochastic matrix $P$ norris1998markov . Each entry
$p_{xy}$ of this matrix $P$ represents the probability of transitioning from
state $x$ to state $y$. Any distribution over the state space of the Markov
chain is represented by a stochastic row vector.

A Markov chain is irreducible if any state can be reached from any other state
in a finite number of steps. Any irreducible Markov chain is aperiodic if
there exists no integer greater than one that divides the length of every
directed cycle of the graph. A Markov chain is ergodic if it is both
irreducible and aperiodic. By the Perron-Frobenius Theorem, any ergodic Markov
chain $P$ has a unique stationary state $\pi$ such that $\pi P=\pi$. The
stationary state $\pi$ is a stochastic row vector and has support on all the
elements of $X$. Let us denote it as

$\pi=\left(\pi_{1}~{}~{}\pi_{2}~{}~{}\cdots~{}~{}\pi_{n}\right),$ (1)

such that $\sum_{j=1}^{n}\pi_{j}=1$. Starting from any initial probability
distribution $\mu$ over the state space $X$, the repeated application of $P$
leads to convergence to the stationary distribution $\pi$, i.e.
$\lim_{t\rightarrow\infty}\mu P^{t}=\pi$. This is known as the mixing of a
Markov chain. It follows from the Perron-Frobenius theorem that other than
$\pi$, all eigenvectors have eigenvalues of absolute value strictly less than
$1$. Thus, $\pi$ is the unique eigenvector with eigenvalue $1$ and all other
eigenvalues lie between $-1$ and $1$. Throughout the paper we shall be working
with the Markov chain corresponding to the lazy walk, i.e. we shall map
$P\mapsto(I+P)/2$. This transformation ensures that all the eigenvalues of $P$
lie between $0$ and $1$. This transformation will not affect our results other
than by a factor of two, which is irrelevant in the asymptotic limit.
Throughout the article, we shall denote the gap between the two highest
eigenvalues of $P$ (the spectral gap) by $\Delta$.

Let $p_{x,y}$ denotes the $(x,y)^{\text{th}}$-entry of the ergodic Markov
chain $P$ with stationary state $\pi$. Then the $(x,y)^{\text{th}}$ entry of
the time-reversed Markov chain of $P$, denoted by $P^{*}$, is

$p^{*}_{x,y}=p_{y,x}\dfrac{\pi_{y}}{\pi_{x}}.$

We shall concern ourselves with ergodic Markov chains that are also
reversible, i.e. Markov chains for which $P=P^{*}$. Any reversible $P$
satisfies the detailed balance condition

$\pi_{x}p_{xy}=\pi_{y}p_{yx},~{}\forall(x,y)\in X.$

This can also be rewritten as

$\text{diag}(\pi)P=P^{T}\text{diag}(\pi),$

where $\text{diag}(\pi)$ is a diagonal matrix with the $j^{\text{th}}$
diagonal entry being $\pi_{j}$. In other words, the reversibility criterion
implies that the matrix $\text{diag}(\pi)P$ is symmetric. Henceforth we shall
only deal with reversible (and hence ergodic) Markov chains.

Interpolated Markov chains: Let us assume that a subset of the elements of the
state space of the Markov chain $P$ is marked. Let $M\subset X$ denote the set
of marked elements. Given any $P$, we define $P^{\prime}$ as the absorbing
Markov chain obtained from $P$ by replacing all the outgoing edges from $M$ by
self-loops. If we re-arrange the elements of $X$ such that the unmarked
elements $U:=X\backslash M$ appear first, then we can write

$\displaystyle P=\begin{bmatrix}P_{UU}&P_{UM}\\\
P_{MU}&P_{MM}\end{bmatrix},~{}~{}~{}~{}~{}~{}~{}~{}~{}P^{\prime}=\begin{bmatrix}P_{UU}&P_{UM}\\\
0&I\end{bmatrix},$ (2)

where $P_{UU}$ and $P_{MM}$ are square matrices of size $(n-|M|)\times(n-|M|)$
and $|M|\times|M|$ respectively. On the other hand $P_{UM}$ and $P_{MU}$ are
matrices of size $(n-|M|)\times|M|$ and $|M|\times(n-|M|)$ respectively. Then
the interpolated Markov chain is defined as

$P(s)=(1-s)P+sP^{\prime},$ (3)

where $s\in[0,1]$. The interpolated Markov chain thus has a block structure

$\displaystyle P=\begin{bmatrix}P_{UU}&P_{UM}\\\
(1-s)P_{MU}&(1-s)P_{MM}+sI\end{bmatrix}.$ (4)

Clearly, $P(0)=P$ and $P(1)=P^{\prime}$. Notice that if $P$ is ergodic, so is
$P(s)$ for $s\in[0,1)$. This is because any edge in $P$ is also an edge of
$P(s)$ and so the properties of irreducibility and aperiodicity are preserved.
However when $s=1$, $P(s)$ has outgoing edges from $M$ replaced by self-loops
and as such the states in $U$ are not accessible from $M$, implying that
$P(1)$ is not ergodic. The spectral gap of $P(s)$ is denoted by $\Delta(s)$.

Now we shall see how the stationary state of $P$ is related to that of $P(s)$.
Since $X=U\cup M$, the stationary state $\pi$ can be written as

$\pi=(\pi_{U}~{}~{}\pi_{M}),$ (5)

where $\pi_{U}$ and $\pi_{M}$ are row-vectors of length $n-|M|$ and $|M|$
respectively. As mentioned previously, $P^{\prime}$ is not ergodic and does
not have a unique stationary state. In fact, any state having support over
only the marked set is a stationary state of $P^{\prime}$.

On the other hand $P(s)$ is ergodic for $s\in[0,1)$. Let $p_{M}=\sum_{x\in
M}\pi_{x}$ be the probability of obtaining a marked element in the stationary
state of $P$. Then it is easy to verify that the unique stationary state of
$P(s)$ is

$\pi(s)=\dfrac{1}{1-s(1-p_{M})}\left((1-s)\pi_{U}~{}~{}\pi_{M}\right).$ (6)

Discriminant matrix:  The discriminant matrix of $P(s)$ is defined as

$D(P(s))=\sqrt{P(s)\circ P(s)^{T}},$ (7)

where $\circ$ indicates the Hadamard product and the $(x,y)^{\mathrm{th}}$
entry of $D(P(s))$ is $D_{xy}(P(s))=\sqrt{p_{xy}(s)p_{yx}(s)}$. Thus $D(P(s))$
is a symmetric matrix.

###### Fact 1

For any ergodic, reversible Markov chain $P$, we have that for $s\in[0,1)$

$D(P(s))=\mathrm{diag}(\sqrt{\pi(s)})P(s)\mathrm{diag}(\sqrt{\pi(s)})^{-1},$

where $\sqrt{\pi(s)}$ is a row vector with its $j^{\text{th}}$-entry being
$\sqrt{\pi_{j}(s)}$.

For any $s\in[0,1)$ as $P(s)$ is reversible, the detailed-balance condition is
satisfied. So, each entry of $D(P(s))$ can be expressed as

$\displaystyle D_{xy}(P(s))$ $\displaystyle=\sqrt{p_{xy}(s)p_{yx}(s)}$ (8)
$\displaystyle=p_{xy}(s)\sqrt{\frac{\pi_{x}(s)}{\pi_{y}(s)}}.$ (9)

From this fact, we obtain that $D(P(s))$ is similar to $P(s)$, i.e. they have
the same set of eigenvalues 222Note that throughout this paper, we shall be
dealing with values of $s\in[0,1)$ and so the properties of $D(P(1))$ are not
relevant here..

Let the spectral decomposition of $D(P(s))$ be

$D(P(s))=\sum_{i=1}^{n}\lambda_{i}(s)\lvert v_{i}(s)\rangle\langle
v_{i}(s)\rvert,$ (10)

where $\lvert v_{i}(s)\rangle$ is an eigenvector of $D(P(s))$ with eigenvalue
$\lambda_{i}(s)$. Furthermore,
$\lambda_{n}(s)=1>\lambda_{n-1}(s)\geq\cdots\geq\lambda_{1}(s)$.

###### Fact 2

For $s\in[0,1)$, the eigenstate of $D(P(s))$ with eigenvalue $1$ is given by

$\lvert v_{n}(s)\rangle=\sqrt{\pi(s)^{T}},$

where $\sqrt{\pi(s)}$ is a row vector with its $j^{\text{th}}$-entry being
$\sqrt{\pi_{j}(s)}$.

This fact follows from the reversibility condition stated in Fact 1, i.e. for
$s\in[0,1)$ we have

$\displaystyle D(P(s))\sqrt{\pi(s)^{T}}$
$\displaystyle=\mathrm{diag}(\sqrt{\pi(s)})P(s)\mathrm{diag}(\sqrt{\pi(s)})^{-1}\sqrt{\pi(s)^{T}}$
(11) $\displaystyle=\sqrt{\pi(s)^{T}}.$ (12)

The $1$-eigenvector of $D(P(s))$ $\lvert v_{n}(s)\rangle$ can also be
expressed in a different form.

###### Proposition 3

The eigenstate of eigenvalue $1$ of $D(P(s))$ can be expressed as

$\lvert v_{n}(s)\rangle=\sqrt{\dfrac{(1-s)(1-p_{M})}{1-s(1-p_{M})}}\lvert
U\rangle+\sqrt{\dfrac{p_{M}}{1-s(1-p_{M})}}\lvert M\rangle,$ (13)

where $\lvert U\rangle$ and $\lvert M\rangle$ are defined as

$\displaystyle\lvert U\rangle$
$\displaystyle=\frac{1}{\sqrt{1-p_{M}}}\sum_{x\notin M}\sqrt{\pi_{x}}\lvert
x\rangle$ (14) $\displaystyle\lvert M\rangle$
$\displaystyle=\frac{1}{\sqrt{p_{M}}}\sum_{x\in M}\sqrt{\pi_{x}}\lvert
x\rangle.$ (15)

This follows directly from Fact 2.

###  II.2 Some quantities related to Markov chains: Hitting and mixing times

In this subsection, we define certain quantities related to Markov chains
which we shall use in subsequent sections for our analysis.

Spatial search problem and hitting time: Consider a graph $G(X,E)$ with
$|X|=n$ vertices and $|E|=e$ edges. Consider a subset $M\subset X$ of vertices
that are marked. Then the spatial search problem involves finding any of the
marked vertices in $M$. This problem can be solved by both classical random
walks and quantum walks.

Given an ergodic, reversible Markov chain $P$ with a stationary state $\pi$,
the random walk based algorithm to solve the spatial search problem is
described in Algorithm 1.

  * 1\.

Sample a vertex $x\in X$ from the stationary state $\pi$ of $P$.

  * 2\.

Check if $x\in M$.

  * 3\.

If $x$ is marked, output $x$.

  * 4\.

Otherwise update $x$ according to $P$ and go to step 2.

Algorithm 1 Spatial search by random walk

The hitting time of $P$ with respect to $M$ is the expected number of times
step $4$ of Algorithm 1 is executed. Let us denote this by $HT(P,M)$. Thus,
the random walk based algorithm finds a marked vertex in time
$\mathcal{O}\left(HT(P,M)\right)$. Note that the random walk algorithm stops
as soon as a marked element is reached. Thus, this is equivalent to applying
an absorbing Markov chain $P^{\prime}$ that is obtained by replacing all the
outgoing edges from the marked vertices of $P$ by self loops. From this we can
define $HT(P,M)$.

Hitting time of a Markov chain: The hitting time of any Markov chain $P$ with
respect to a set of marked elements $M$ can be expressed as

$HT(P,M)=\sum_{j=1}^{n-m}\dfrac{|\langle
v^{\prime}_{j}|U\rangle|^{2}}{1-\lambda^{\prime}_{j}},$ (16)

where $\lambda^{\prime}_{j}$ and $\lvert v^{\prime}_{j}\rangle$ are the
eigenvalues and eigenvectors of the matrix $D(P^{\prime})$ and

$\lvert U\rangle=\dfrac{1}{\sqrt{1-p_{M}}}\sum_{x\neq M}\sqrt{\pi_{x}}\lvert
x\rangle,$

where $p_{M}$ is the probability of sampling a marked vertex from the
stationary state of $P$.

Interpolated hitting time and Extended hitting time: For any interpolated
Markov chain $P(s)$, in refs. krovi2016quantum ; chakraborty2018finding , the
authors define a quantity known as the interpolated hitting time in the
context of spatial search which will also be useful here for subsequent
analysis. This is defined as

$HT(s)=\sum_{j=1}^{n-1}\dfrac{|\langle
v_{j}(s)|U\rangle|^{2}}{1-\lambda_{j}(s)}.$ (17)

There is a relationship between the spectral gap of the Markov chain and
$HT(s)$ since

$HT(s)\leq\dfrac{1}{\Delta(s)}\sum_{j=1}^{n-1}|\langle
v_{j}(s)|U\rangle|^{2}.$ (18)

For the spatial search algorithm, we shall find that the quantity that is of
interest is the extended hitting time. The extended hitting time of $P$ with
respect to a set $M$ of marked elements is given by

$HT^{+}(P,M)=\lim_{s\rightarrow 1}HT(s),$ (19)

Clearly for $|M|=1$, we have that $HT^{+}(P,M)=HT(P,M)$. Krovi et al. proved
an explicit relationship between $HT(s)$ and $HT^{+}(P,M)$ krovi2016quantum .
They showed that

$HT(s)=\dfrac{p_{M}^{2}}{\left(1-s(1-p_{M})\right)^{2}}HT^{+}(P,M).$ (20)

Combining Eqs. (18) and (20), we have

$HT^{+}(P,M)\leq\dfrac{1}{\Delta(s)}.\dfrac{\left(1-s(1-p_{M})\right)^{2}}{p_{M}^{2}}\sum_{j=1}^{n-1}|\langle
v_{j}(s)|U\rangle|^{2}.$ (21)

Mixing-time of a Markov chain: Given a reversible Markov chain $P$, any
initial probability distribution over the state space converges to the
stationary distribution $\pi$, i.e. $\lim_{t\rightarrow\infty}\mu=\pi$, for
any initial distribution $\mu$. Given $P$ and an initial state $\mu$, the
mixing-time of a classical random walk is defined as the minimum time
$T_{\mathrm{mix}}$ such that $\forall t\geq T_{\mathrm{mix}}$ we have that

$\dfrac{1}{2}\|\mu P^{t}-\pi\|_{1}\leq\epsilon,$

for some $\epsilon\in(0,1)$, where $\frac{1}{2}\|.\|_{1}$ is the total
variation distance.

That is, $T_{\mathrm{mix}}$ is the minimum time required for the Markov chain
to converge to a distribution that is $\epsilon$-close to the stationary
distribution which implies that aldous1982some :

$T_{\mathrm{mix}}\leq\dfrac{1}{\Delta}\log\left(\dfrac{1}{\epsilon\pi_{\min}}\right),$
(22)

where $\Delta$ is the spectral gap of $P$ and $\pi_{\min}=\min_{x}\pi_{x}$.
Thus, given an ergodic, reversible Markov chain $P$ with stationary state
$\pi$ and spectral gap $\Delta$, one can sample from a distribution that is
$\epsilon$-close to $\pi$ in time
$\widetilde{\mathcal{O}}\left(1/\Delta\right)$. Next we discuss how one can
use von Neumann measurements to prepare eigenstates of Hamiltonians, a tool
which will help us provide an analog quantum algorithm for solving $QSSamp$.

##  III Quantum state generation by von Neumann measurements

In this section, we make use of von Neumann measurements to prepare
eigenstates of a Hamiltonian. The goal would be to use this technique to
prepare the eigenstate of the quantum walk Hamiltonian (encoding an ergodic
reversible Markov chain $P$) that corresponds to a coherent encoding of the
stationary distribution of $P$.

In this framework, in order to measure any observable $\widehat{O}$, the
system of interest is coupled to a pointer, which is simply a free particle in
one dimension. If $H$ represents the Hamiltonian of the system and
$\widehat{p}$ the momentum operator corresponding to the pointer, then the
total Hamiltonian corresponding to the coupling between the system and the
pointer is given by

$\widetilde{H}=H+\dfrac{\widehat{p}^{2}}{2m}+g\
\widehat{O}\otimes\widehat{p},$ (23)

where $m$ is the mass of the free particle and $g$ is the interaction strength
between the observable and the pointer. Since we are interested in measuring
the energy of the system, we have $\widehat{O}=H$. We consider the particle as
“massive”, thereby enabling us to neglect the free Hamiltonian of the
particle. Furthermore, we assume that we are working with units such that the
interaction strength $g=1$. These imply that

$\widetilde{H}=H\otimes\widehat{p}.$ (24)

It is well known that the momentum operator, $\widehat{p}=-i\frac{d}{dx}$ is a
generator of translation in the position of the particle. In other words, the
operator $e^{-ix_{0}\widehat{p}}$ applied to a wavepacket whose wavefunction
is $\psi(x)$ results in

$\displaystyle e^{-ix_{0}\widehat{p}}\psi(x)$
$\displaystyle=e^{-x_{0}\frac{d}{dx}}\psi(x)$ (25)
$\displaystyle=\left(I-x_{0}\frac{d}{dx}+...\right)\psi(x)$ (26)
$\displaystyle=\psi(x-x_{0}).$ (27)

Thus the wavepacket is translated in position by $x_{0}$. Now consider that
the system Hamiltonian $H$ has eigenvalues

$\lambda_{n}=0<\lambda_{n-1}<\Delta\leq\lambda_{n-2}\leq\cdots\lambda_{1}\leq
1,$

such that $H\lvert v_{j}\rangle=\lambda_{j}\lvert v_{j}\rangle$. Furthermore,
suppose that we initialize the pointer to a state $\lvert x=0\rangle$, a
wavepacket centred around $0$. Then,

$e^{-i\widetilde{H}t}\lvert\lambda_{j}\rangle\lvert x=0\rangle=\lvert
v_{j}\rangle\lvert x=\lambda_{j}t\rangle.$ (28)

That is, the wavepacket is translated in position by $\lambda_{j}t$ and as
such, measuring the displacement of the pointer register can in principle
reveal information about the eigenstate of $H$ in the first register. By
linearity, for any initial state
$\lvert\psi_{0}\rangle=\sum_{j=1}^{n}\alpha_{j}\lvert v_{j}\rangle$, we have

$\displaystyle e^{-i\widetilde{H}t}\lvert\psi_{0}\rangle\lvert x=0\rangle$
$\displaystyle=e^{-iH\widehat{p}t}\sum_{j=1}^{n}\lvert\psi_{0}\rangle\lvert
x=0\rangle$ (29) $\displaystyle=\sum_{j=1}^{n}\alpha_{j}\lvert
v_{j}\rangle\lvert x=\lambda_{j}t\rangle.$ (30)

In order to implement this on a quantum computer, we assume that the pointer
register is of $l$-qubits. The choice of $l$ is crucial as it determines the
precision up to which the position of the pointer is obtained. In fact, if we
measure the position of the pointer with a high enough precision to resolve
all eigenvalue gaps, $(\lambda_{i}-\lambda_{j})t$, a measurement of the
position of the pointer results in a measurement of the system Hamiltonian
$H$.

For our purposes, we shall show how this formalism can be used to prepare the
$0$-eigenstate of $H$, i.e. $\lvert v_{n}\rangle$, in a purely analog fashion.
In other words, we demonstrate via Lemma 4 and Corollary 5 that by starting
from the state $\lvert\psi_{0}\rangle$ in the system register and $\lvert
x=0\rangle$ in the pointer register, and performing time evolutions of the
Hamiltonian $\widetilde{H}$ to prepare a quantum state that is
$\epsilon$-close to $\lvert v_{n}\rangle$ with probability $|\alpha_{n}|^{2}$.

###### Lemma 4

Let $H$ be a Hamiltonian with eigenvalues
$\lambda_{n}=0<\Delta<\lambda_{n-1}\leq\cdots\lambda_{1}\leq 1$ such that
$H\lvert v_{j}\rangle=\lambda_{j}\lvert v_{j}\rangle$. Let $\widehat{p}$
represent the momentum operator corresponding to a free particle in one
dimension with its mass large enough so that its free Hamiltonian can be
neglected and so that it can be represented in $l$-qubits as

$\widehat{p}=\sum_{q=0}^{2^{l}-1}\dfrac{q}{2^{l}}\lvert q\rangle\langle
q\rvert,$

where

$l=\lceil\log_{2}(1/\Delta)\rceil+1.$ (31)

Furthermore let

$\lvert\psi_{0}\rangle=\sum_{j=1}^{n}{\alpha_{j}}\lvert v_{j}\rangle.$

Then starting from the state $\lvert\psi_{0}\rangle\lvert x=0\rangle$ and
evolving for a time

$\tau=\dfrac{2\pi}{\Delta},$

according to the Hamiltonian $\widetilde{H}=H\otimes\widehat{p}$, results in a
state

$\lvert\widetilde{\psi}\rangle=\alpha_{n}\lvert v_{n}\rangle\lvert
0\rangle+\sum_{k=1}^{n-1}\alpha_{k}\lvert v_{k}\rangle\left(\gamma_{k}\lvert
0\rangle+\Gamma_{k}\lvert\Gamma_{k}\rangle\right),$

where $|\gamma_{k}|<1/2,~{}|\Gamma_{k}|>\sqrt{3}/2$ and
$\langle\Gamma_{k}|0\rangle=0$ for $1\leq k\leq n-1$.

Proof.

If $\lvert q\rangle$ represents the momentum eigenstates, then the momentum
operator is represented by

$\widehat{p}=\sum_{m=0}^{2^{l}-1}\dfrac{q}{2^{l}}\lvert q\rangle\langle
q\rvert.$ (32)

Note that the position and momentum states are equivalent up to a Fourier
transform and so the localized wavepacket centred at $x=0$ is completely
delocalized in the momentum basis. That is,

$\lvert x=0\rangle=\dfrac{1}{\sqrt{2^{l}}}\sum_{q=0}^{2^{l}-1}\lvert
q\rangle.$ (33)

Therefore,

$\displaystyle e^{-i(H\otimes\widehat{p})\tau}\lvert\psi_{0}\rangle\lvert
x=0\rangle$
$\displaystyle=e^{-i(H\otimes\widehat{p})\tau}\lvert\psi_{0}\rangle\left(\dfrac{1}{\sqrt{2^{l}}}\sum_{q=0}^{2^{l}-1}\lvert
q\rangle\right)$ (34) $\displaystyle=\sum_{k=1}^{n}\alpha_{k}\lvert
v_{k}\rangle\left(\dfrac{1}{\sqrt{2^{l}}}\sum_{q=0}^{2^{l}-1}e^{\frac{-i\lambda_{k}\tau
q}{2^{l}}}\lvert q\rangle\right).$ (35)

Since we ultimately want to read off the position of the pointer variable, we
re-express the pointer register in the position-basis to obtain

$\displaystyle\begin{split}&e^{-i(H\otimes\widehat{p})\tau}\lvert\psi_{0}\rangle\lvert
x=0\rangle=\\\ &\sum_{k=1}^{n}\alpha_{k}\lvert
v_{k}\rangle\left(\dfrac{1}{2^{l}}\sum_{x=0}^{2^{l}-1}\sum_{q=0}^{2^{l}-1}e^{\frac{i(x-\lambda_{k}\tau)q}{2^{l}}}\lvert
x\rangle\right).\end{split}$ (36)

The pointer register has a measure of the displacement of the wavepacket which
was initially centered at $x=0$. In fact, as shown previously, the shift will
be proportional to the eigenstate in the first register (expressed in $l$
qubits). That is, we will have states of the form $\lvert
v_{j}\rangle\lvert\lambda_{j}\tau\rangle$. We are interested in preparing the
$0$-eigenstate $\lvert v_{n}\rangle$. We first observe that the amplitude of
obtaining $\lvert 0\rangle$ in the pointer register when the first register is
in the state $\lvert v_{n}\rangle\lvert 0\rangle$ is one, i.e.

$e^{-i\tau(H\otimes\widehat{p})}\lvert v_{n}\rangle\lvert
0\rangle\mapsto\lvert v_{n}\rangle\lvert 0\rangle.$

On the other hand, for any other eigenstate $\lvert v_{k}\rangle$, the
amplitude corresponding to measuring $\lvert 0\rangle$ in the second register
is

$\displaystyle\dfrac{1}{2^{l}}\left|\sum_{q=0}^{2^{l}-1}e^{\frac{i(x-\lambda_{k}\tau)q}{2^{l}}}\right|$
$\displaystyle=\dfrac{1}{2^{l}}\left|\sum_{q=0}^{2^{l}-1}e^{\frac{-i\lambda_{k}\tau
q}{2^{l}}}\right|$ (37)
$\displaystyle=\dfrac{1}{2^{l}}\left|\dfrac{1-e^{-i\lambda_{k}\tau}}{1-e^{-i\lambda_{k}\tau/2^{l}}}\right|$
(38)
$\displaystyle=\dfrac{1}{2^{l}}\left|\dfrac{1-e^{-i2\pi\lambda_{k}/\Delta}}{1-e^{-i2\pi\lambda_{k}/(2^{l}\Delta)}}\right|,$
(39)

where the last line follows from the fact that $\tau=2\pi/\Delta$. Let
$z=\frac{2\pi\lambda_{k}}{2^{l}\Delta}$. Then, we find that

$\displaystyle 0\leq z\leq\frac{2\pi}{2^{l}\Delta}\leq\pi,$ (40)

where we have used the fact that $l\geq\log_{2}(1/\Delta)+1$ and so
$2^{l}\Delta\geq 2$. Now we use the following facts: for $z\in[-\pi,\pi]$,
$|\sin(z/2)|\geq|z|/\pi$ which gives us that $|1-e^{-iz}|=2|\sin(z/2)|\geq
2|z|/\pi$. Also, observe that $|1-e^{-i2\pi\lambda_{k}/\Delta}|\leq 2$.

So combining these two, we obtain that

$\dfrac{1}{2^{l}}\left|\dfrac{1-e^{-i2\pi\lambda_{k}/\Delta}}{1-e^{-i2\pi\lambda_{k}/(2^{l}\Delta)}}\right|\leq\dfrac{\Delta}{2|\lambda_{k}|}<1/2,$
(41)

where the last expression follows from the fact that
$|\lambda_{k}|>\Delta,\forall 1\leq k\leq n-1$. This immediately implies that
the amplitude of the pointer register to be in a state different from $\lvert
0\rangle$ when the first register is in $\lvert v_{k}\rangle$ is at least
$\sqrt{3}/2$. In other words, in such a scenario the state of the pointer
register, denoted by $\lvert\Gamma_{k}\rangle$ will have at least one non-zero
qubit, ensuring that $\langle\Gamma_{k}|0\rangle=0$.

Thus, after the time evolution for a time $\tau=2\pi/\Delta$, the state of the
system and the pointer is given by

$\lvert\widetilde{\psi}\rangle=\alpha_{n}\lvert v_{n}\rangle\lvert
0\rangle+\sum_{k=1}^{n-1}\alpha_{k}\gamma_{k}\lvert v_{k}\rangle\lvert
0\rangle+\sum_{k=1}^{n-1}\alpha_{k}\Gamma_{k}\lvert
v_{k}\rangle\lvert\Gamma_{k}\rangle,$ (42)

where $\forall k\in[1,n-1]$, we have $|\gamma_{k}|<1/2$,
$|\Gamma_{k}|>\sqrt{3}/2$ and $\langle\Gamma_{k}|0\rangle=0$. $\Box$

We shall use this lemma to derive the following corollary.

###### Corollary 5

Let $\epsilon^{\prime}=\epsilon|\alpha_{n}|^{2}$ where $\epsilon\in(0,1)$ and
suppose that the pointer register contains

$m=l\ .\ \lceil\log_{2}(1/\epsilon^{\prime})\rceil$ (43)

qubits initialized in the state $\lvert x=0\rangle^{\otimes m}$.

Then repeating the Hamiltonian evolution of Lemma 4 a total of
$\lceil\log(1/\epsilon^{\prime})\rceil$-times using a fresh block of
$l$-pointer qubits each time, followed by post-selecting on the pointer
register to be in $\lvert 0\rangle^{\otimes{m}}$, results in a quantum state
$\lvert\phi\rangle$ such that

$\left\lVert\lvert
v_{n}\rangle-\lvert\phi\rangle\right\rVert_{2}\leq\Theta(\epsilon),$

in time

$T=\Theta\left(\dfrac{1}{\Delta|\alpha_{n}|^{2}}\log\left(\frac{1}{\epsilon|\alpha_{n}|^{2}}\right)\right).$
(44)

Proof.  After the application of $e^{-i(H\otimes\widehat{p})\tau}$ a total of
$\lceil\log(1/\epsilon^{\prime})\rceil$-times using $l$ blocks of qubits in
the pointer register each time, observe that for any $k\neq n$, the amplitude
for observing $\lvert 0\rangle^{\otimes m}$ in the pointer register when there
is $\lvert v_{k}\rangle$ in the first register is bounded by

$|\epsilon_{k}|=|\gamma_{k}|^{\lceil\log_{2}(1/\epsilon^{\prime})\rceil}\leq\left(\dfrac{1}{2}\right)^{\lceil\log_{2}(1/\epsilon^{\prime})\rceil}\leq\epsilon^{\prime}.$
(45)

This implies that the resulting state after this procedure is given by

$\displaystyle\begin{split}&\lvert\psi_{f}\rangle=\alpha_{n}\lvert
v_{n}\rangle\lvert 0\rangle^{\otimes
m}+\sum_{k=1}^{n-1}\epsilon_{k}\alpha_{k}\lvert v_{k}\rangle\lvert
0\rangle^{\otimes m}\\\ &~{}~{}~{}+\sum_{k=1}^{n-1}\alpha_{k}\delta_{k}\lvert
v_{k}\rangle\lvert\Gamma_{k}^{(m)}\rangle,\end{split}$ (46)

where $0\leq\epsilon_{k}\leq\epsilon^{\prime}$ and $\sqrt{1-\epsilon^{\prime
2}}\leq\delta_{k}\leq 1$. This takes time
$2\pi\lceil\log(1/\epsilon^{\prime})\rceil/\Delta$.

The state in Eq. (46) can be re-written as

$\lvert\psi_{f}\rangle=\alpha_{n}\left(\lvert
v_{n}\rangle+\lvert\text{err}\rangle\right)\lvert 0\rangle^{\otimes
m}+\sum_{k=1}^{n-1}\alpha_{k}\delta_{k}\lvert
v_{k}\rangle\lvert\Gamma_{k}^{(m)}\rangle,$ (47)

where the (unnormalized) state

$\lvert\text{err}\rangle=\sum_{k=1}^{n-1}\dfrac{\epsilon_{k}\alpha_{k}}{\alpha_{n}}\lvert
v_{k}\rangle.$

This implies that post-selecting on obtain $\lvert 0\rangle^{\otimes m}$ in
the pointer register we obtain the state

$\lvert\phi\rangle=\lvert v_{n}\rangle+\lvert\text{err}\rangle,$ (48)

with probability $|\alpha_{n}|^{2}$, such that

$\left\lVert\lvert
v_{n}\rangle-\lvert\phi\rangle\right\rVert^{2}=\|\lvert\text{err}\rangle\|^{2}=\sum_{k=1}^{n-1}\left|\dfrac{\epsilon_{k}\alpha_{k}}{\alpha_{n}}\right|^{2}\leq\dfrac{\epsilon^{\prime
2}}{|\alpha_{n}|^{2}}=\epsilon^{2}.$ (49)

The entire protocol takes time

$T=\Theta\left(\dfrac{1}{\Delta|\alpha_{n}|^{2}}\log\left(\frac{1}{|\alpha_{n}|^{2}\epsilon}\right)\right).$

$\Box$

Thus Lemma 4 and corollary 5 can be used to prepare the eigenstate $\lvert
v_{n}\rangle$. Note that it would have been possible to use quantum amplitude
amplification to reduce quadratically the dependency on $|\alpha_{n}|$.
However, we are interested in developing analog algorithms, assuming that we
have access to a time-independent Hamiltonian. Provided that the cost of
preparing the initial state $\lvert\psi_{0}\rangle$ is small, the cost of the
algorithm is the total time of Hamiltonian evolution. Moreover our protocol
(Sec. VI) to prepare the stationary state of any reversible Markov chain
ensures that $|\alpha_{n}|=\Theta(1)$, thereby resulting in at most a constant
slowdown with respect to amplitude amplification.

##  IV Hamiltonian for quantum walk on any reversible Markov chain

Given any ergodic, reversible Markov chain, we shall make use of the
Hamiltonian introduced by Somma and Ortiz somma2010quantum  and subsequently
used in Refs. krovi2010adiabatic ; chakraborty2018finding . We recall the
Hamiltonian and its spectral properties here for completeness and it will be
used in our quantum algorithm for $QSSamp$.

###  IV.1 Defining the Hamiltonian

Let $p_{xy}(s)$ denote the $(x,y)^{\mathrm{th}}$-entry of $P(s)$ and let $E$
be the set of edges of $P(s)$. Furthermore let
$\mathcal{H}=\mathrm{span}\\{\lvert x\rangle:x\in X\\}$ . Then one can define
a unitary $V(s)\in\mathcal{H}\times\mathcal{H}$ such that for all $x\in X$,

$V(s)\lvert x,0\rangle=\sum_{y\in X}\sqrt{p_{xy}(s)}\lvert x,y\rangle,$ (50)

where the state $\lvert 0\rangle$ represents a fixed reference state in
$\mathcal{H}$. Let us also define the swap operator

$S\lvert x,y\rangle=\begin{cases}\lvert y,x\rangle,&\text{if $(x,y)\in E$}\\\
\lvert x,y\rangle,&\text{otherwise}.\end{cases}$

Observe that $\langle
x,0|V(s)^{\dagger}SV(s)|y,0\rangle=\sqrt{p_{yx}(s)p_{xy}(s)}=D_{xy}(P(s))$.
Then, if $\Pi_{0}=I\otimes\lvert 0\rangle\langle 0\rvert$, we have,

$V^{\dagger}(s)SV(s)\Pi_{0}\lvert y,0\rangle=\sum_{x\in
X}\sqrt{p_{yx}(s)p_{xy}(s)}\lvert x,0\rangle+\lvert\Phi\rangle^{\perp},$ (51)

so that $\Pi_{0}\lvert\Phi\rangle^{\perp}=0$. We define the search Hamiltonian
as

$H(s)=i[V(s)^{\dagger}SV(s),\Pi_{0}].$ (52)

In Ref. chakraborty2018finding , we have shown that $H(s)$, in a rotated
basis, corresponds to a quantum walk on the edges of $P(s)$. That is, the
rotated Hamiltonian

$\displaystyle\overline{H}(s)$ $\displaystyle=V(s)H(s)V(s)^{\dagger}$ (53)
$\displaystyle=i[S,V\Pi_{0}V^{\dagger}],$ (54)

corresponds to a quantum walk on the edges of $P(s)$. If the walker is
localized in a directed edge from node $x$ to node $y$, i.e. $\lvert
x,y\rangle$, then the walker can move to a superposition of outgoing edges
from node $y$ of the form $\lvert y,.\rangle$. Note that our algorithms (See
Algorithm 2 and Algorithm 3) could be implemented using the Hamiltonian
$\overline{H}(s)$ instead of $H(s)$. In such a case, we need to apply the same
rotation to the initial state of the algorithm and the final state of the
algorithm. However, subsequently we shall be working with $H(s)$ as it
simplifies the analysis considerably. In the next subsection, we will
characterize the spectrum of $H(s)$.

###  IV.2 Spectrum of $H(s)$

As discussed in Sec. II.1, the spectrum of $H(s)$ is related to that of
$D(P(s))$ and in particular, the state $\lvert v_{n}(s),0\rangle$ is an
eigenstate of $H(s)$ with eigenvalue zero. The spectrum of $H(s)$ has been
explicitly described in Ref. krovi2010adiabatic  and we mention it here for
completeness. The total Hilbert space of $H(s)$ can be divided into the
following set of invariant subspaces:

For $1\leq k\leq n-1$,

$\displaystyle\mathcal{B}_{k}(s)$ $\displaystyle=\mathrm{span}\\{\lvert
v_{k}(s),0\rangle,V(s)^{\dagger}SV(s)\lvert v_{k}(s),0\rangle\\},$ (55)
$\displaystyle\mathcal{B}_{n}(s)$ $\displaystyle=\mathrm{span}\\{\lvert
v_{n}(s),0\rangle\\}$ (56) $\displaystyle\mathcal{B}^{\perp}(s)$
$\displaystyle=(\oplus_{k=1}^{n}\mathcal{B}_{k})^{\perp}.$ (57)

Now, observe that

$\displaystyle\Pi_{0}V(s)^{\dagger}SV(s)\lvert v_{n}(s),0\rangle=\lvert
v_{n}(s),0\rangle$ (58) $\displaystyle V(s)^{\dagger}SV(s)\Pi_{0}\lvert
v_{n}(s),0\rangle=\lvert v_{n}(s),0\rangle.$ (59)

This implies

$H(s)\lvert v_{n}(s),0\rangle=0,$ (60)

i.e. $\lvert v_{n}(s),0\rangle$ is an eigenstate with eigenvalue $0$.

On the other hand, note that for $1\leq k\leq n-1$, the eigenstates and
eigenvalues of $H(s)$ in $\mathcal{B}_{k}(s)$ are

$\displaystyle\lvert\Psi^{\pm}_{k}(s)\rangle$ $\displaystyle=\dfrac{\lvert
v_{k}(s),0\rangle\pm i\lvert v_{k}(s),0\rangle^{\perp}}{\sqrt{2}},$ (61)
$\displaystyle E^{\pm}_{k}(s)$ $\displaystyle=\pm\sqrt{1-\lambda_{k}(s)^{2}},$
(62)

where $\lvert v_{k}(s),0\rangle^{\perp}$ is a quantum state that is in
$\mathcal{B}_{k}(s)$ such that $\Pi_{0}\lvert v_{k}(s),0\rangle^{\perp}=0$.
Thus, if the underlying Markov chain has a spectral gap $\Delta(s)$, then in
this subspace $H(s)$ has a quadratically amplified spectral gap given by

$|E_{n}(s)-E^{\pm}_{n-1}(s)|=\sqrt{1-\lambda^{2}_{n-1}(s)}=\Theta\left(\sqrt{\Delta(s)}\right).$
(63)

Now, there are $n^{2}$ eigenvalues of $H(s)$ out of which $2n-1$ belong to
$\mathcal{B}_{k}(s)\cup\mathcal{B}_{n}(s)$. The remaining $(n-1)^{2}$
eigenvalues are $0$ and belong to $\mathcal{B}^{\perp}(s)$ which is the
orthogonal complement of the union of the invariant subspaces. We need not
worry about this subspace as we start from a state that has no support on
$B^{\perp}(s)$ which is an invariant subspace of $H(s)$. Thus, throughout the
evolution under $H(s)$, our dynamics will be restricted to
$\mathcal{B}_{k}(s)\cup\mathcal{B}_{n}(s)$.

##  V Spatial search by continuous-time quantum walk using von-Neumann
measurements

We first show how to make use of the state-generation scheme described in Sec.
III to provide a continuous-time quantum walk based algorithm to solve the
spatial search problem. This algorithm will provide an intuitive understanding
of our analog quantum algorithm for $QSSamp$.

Suppose we are given an ergodic, reversible Markov chain $P$ with the marked
set denoted by $M\subset X$. The spatial search algorithm on $P$ involves
finding a node within this marked set and is often tackled by the formalism of
random walks. We have seen previously in Sec. II.2 that the expected number of
steps taken by the walker to find a node within this marked set is known as
the hitting time of $P$ with respect to $M$. Quantum walks provide a natural
framework to tackle this problem. A natural question to ask is whether a
quantum walk can offer any speedup over its classical counterpart in order to
solve the spatial search problem. Here, we concentrate on the continuous-time
quantum walk framework to tackle this problem.

The spatial search algorithm by continuous-time quantum walk on $P$ involves
evolving a time-independent Hamiltonian (which encodes the connectivity of
$P$), starting from some initial state, for some time, and then measuring in
the basis spanned by the states of $P$.

Childs and Goldstone childs2004spatial  introduced the first continuous-time
quantum walk-based algorithm to tackle the spatial search problem for simple,
unweighted graphs. They showed that this algorirhtm, defined as a quantum walk
on the nodes of the underlying graph, could find a marked node in
$\mathcal{O}(\sqrt{n})$ time for certain graphs with $n$ nodes such as the
complete graph, hybercube and $d$-dimensional lattices with $d>4$. This
offered a quadratic speedup over classical random walks for the spatial search
problem on these graphs. When $d=4$, the running time of the Childs and
Goldstone algorithm is $\mathcal{O}(\sqrt{n}\log n)$, offering a less than
quadratic speedup whereas there is no substantial speedup for $d<4$. Since
then, a plethora of results have been published exhibiting a
$\mathcal{O}(\sqrt{n})$ running time on certain specific graphs each requiring
an ad-hoc analysis janmark2014global ; foulger2014quantum ; childs2014spatial-
crystal ; meyer2015connectivity ; novo2015systematic ; philipp2016continuous ;
wong2016quantum ; chakraborty2016spatial ; chakraborty2017optimal . Recently
in Ref. chakraborty2020optimality , the authors provided the necessary and
sufficient conditions for this algorithm to be optimal under very general
conditions on the spectrum of the quantum walk Hamiltonian. They showed that
attaining a generic quadratic speedup is impossible using this algorithm.

In Ref. chakraborty2018finding , the authors provided a different spatial
search algorithm by continuous-time quantum walk which finds a marked element
on any ergodic, reversible Markov chain in square-root of the extended hitting
time, thereby matching the running time of best known algorithms in the DTQW-
framework in the case of where a single node is marked, i.e. $|M|=1$
krovi2016quantum ; ambainis2019quadratic . Given any $P$, their algorithm made
use of the framework of interpolating quantum walks

$P(s)=(1-s)P+sP^{\prime},$

where $P^{\prime}$ is the absorbing Markov chain such that any state having
support only over the marked vertices is its stationary state. In fact, they
used the Somma-Ortiz Hamiltonian $H(s)$ (described in Sec. IV), to define
quantum walk on the edges of $P(s)$. The underlying technique is to use a
procedure called quantum phase randomization to (approximately) prepare a
(mixed) state that has a constant overlap with the $0$-eigenstate of $H(s)$.
For some specific value of $s$, this eigenstate has a constant overlap with
the marked subspace $M$. This required measurement at a time chosen uniformly
at random between $[0,~{}\sqrt{HT^{+}(P,M)}]$, where $HT^{+}(P,M)$ is the
extended hitting time, defined in Eq. (19).

In this section, we provide an alternative spatial search algorithm (Algorithm
2) by continuous-time quantum walk that finds an element from a marked set $M$
in time that scales as the square root of the extended hitting time. Algorithm
2 is similar in spirit to that of Ref. chakraborty2018finding  in that both
make use of the Somma-Ortiz Hamiltonian $H(s)$ defined in Eq. (52).

Let $\epsilon\in(0,1)$ and $\tau=2\pi/\sqrt{\Delta(s^{*})}$.

  * 1\.

Prepare the state $\lvert\psi_{0}\rangle=\lvert v_{n}(0),0\rangle\lvert
x=0\rangle^{\otimes\lceil\log_{2}(2/\epsilon)\rceil}$.

  * 2\.

Evolve according to $H(s^{*})\otimes\widehat{p}$ for time $\tau$ starting from
the state $\lvert\psi_{0}\rangle$.

  * 3\.

Repeat the Hamiltonian evolution in step 2, $m(s^{*})$ times, using a fresh
block of $l(s^{*})$-qubits in the pointer register each time.

  * 4\.

Measure in the basis of the state-space of the Markov chain in the first
register.

Algorithm 2 Spatial search by continuous-time quantum walk

However, motivated by the problem of quantum state generation using von
Neumann measurements, Algorithm 2 prepares the $0$-eigenstate of $H(s)$ by
coupling this Hamiltonian to a free-particle in one dimension. Unlike the
algorithm of Ref. chakraborty2018finding , we evolve the Hamiltonian for a
fixed time before making a measurement.

Note that the spectral gap of $H(s)$ is quadratically less than that of
$P(s)$. That is, if the discriminant matrix $D(P(s))$ has a spectral gap of
$\Delta(s)=|1-\lambda_{n-1}(s)|$, then the spectral gap of $H(s)$ is
$\sqrt{1-\lambda_{n-1}(s)^{2}}=\Theta(\sqrt{\Delta(s)})$.

Furthermore, for

$s=s^{*}=1-p_{M}/(1-p_{M}),$ (64)

the $0$-eigenstate $\lvert v_{n}(s^{*})\rangle$ can be written as

$\lvert v_{n}(s^{*})\rangle=\dfrac{\lvert U\rangle+\lvert
M\rangle}{\sqrt{2}},$ (65)

where $\lvert U\rangle$ and $\lvert M\rangle$ are as defined in Proposition 3.
Thus it has a constant overlap with both $\lvert U\rangle$ and $\lvert
M\rangle$. Consider the initial state

$\lvert v_{n}(0)\rangle=\sqrt{1-p_{M}}\lvert U\rangle+\sqrt{p_{M}}\lvert
M\rangle.$ (66)

Then from Lemma 4 and corollary 5, we choose

$l(s^{*})=\lceil\log_{2}(1/\Delta(s^{*}))\rceil+1,$ (67)

and

$m(s^{*})=l(s^{*})\ .\ \lceil\log_{2}(1/\epsilon^{\prime})\rceil,$ (68)

such that the time evolution of the Hamiltonian
$\widetilde{H}(s^{*})=H(s^{*})\otimes\widehat{p}$, starting from the state

$\lvert\psi_{0}\rangle=\lvert v_{n}(0)\rangle\lvert
x=0\rangle^{\otimes\lceil\log_{2}(2/\epsilon)\rceil},$

prepares a state that in $\epsilon$-close $\lvert v_{n}(s^{*})\rangle$ with
probability

$|\alpha_{n}|^{2}=|\langle
v_{n}(0)|v_{n}(s^{*})\rangle|^{2}\approx\dfrac{1}{2}.$ (69)

in time

$T=\Theta\left(\dfrac{1}{\sqrt{\Delta(s^{*})}}\log\dfrac{1}{\epsilon}\right).$

We prove in Lemma 6 that $\Delta(s^{*})$ is upper bounded by $HT^{+}(P,M)$.
Furthermore as $\lvert v_{n}(s^{*})\rangle$ has a constant overlap with
$\lvert M\rangle$, Algorithm 2 finds a marked node with a constant probability
in time

$T=\Theta\left(\sqrt{HT^{+}(P,M)}\log\dfrac{1}{\epsilon}\right).$

Now formally state Algorithm 2 and prove its correctness in Lemma 6.

###### Lemma 6

Algorithm 2 outputs a marked node with probability at least $1/4-\epsilon$ in
time

$T=\Omega\left(\sqrt{HT^{+}(P,M)}\log(1/\epsilon)\right).$

Proof.  We shall make use of Lemma 4 and Corollary 5. Observe that for
$s=s^{*}=1-p_{M}/(1-p_{M})$, the $0$-eigenstate of $H(s^{*})$ is simply

$\lvert v_{n}(s^{*})\rangle=\dfrac{\lvert U\rangle+\lvert
M\rangle}{\sqrt{2}},$ (70)

where $\lvert U\rangle$ and $\lvert M\rangle$ are as defined in Proposition 3.
Also, the initial state in the first register is

$\lvert v_{n}(0)\rangle=\sqrt{1-p_{M}}\lvert U\rangle+\sqrt{p_{M}}\lvert
M\rangle.$ (71)

Let $\alpha_{n}$ be the overlap of $\lvert v_{n}(s^{*})\rangle$ with $\lvert
v_{n}(0)\rangle$. Then,

$|\alpha_{n}|^{2}=\dfrac{1}{2}+\sqrt{p_{M}(1-p_{M})}.$ (72)

Clearly $|\alpha_{n}|^{2}\geq 1/2$. Also note that $|\langle
v_{n}(s^{*})|M\rangle|^{2}=1/2$. Consider the measurement operator

$\mathcal{M}=\Pi_{X}\otimes\lvert 0\rangle\langle 0\rvert\otimes\lvert
0\rangle\langle 0\rvert,$ (73)

where $\Pi_{X}$ is a projection on the states of the Markov chain. Thus, from
Lemma 4 and Corollary 5, we have that after executing steps $2$ and $3$,
conditioned on having $\lvert 0\rangle^{\otimes n}$ in the second register and
$\lvert 0\rangle^{\otimes m}$ in the pointer register, we end up in a state
such that when measured using the operator $\mathcal{M}$, we obtain a marked
element with probability at least $1/4-\epsilon$.

The total time required to execute steps $2$ and $3$ is

$T=\Theta\left(\dfrac{1}{\sqrt{\Delta(s^{*})}}\log(1/\epsilon)\right).$

Now we use the relationship between $\Delta(s^{*})$ and $HT^{+}(P,M)$ defined
in Eq. (21) for $s=s^{*}$.

Observe that for any $1\leq k\leq n-1$, $|\langle
v_{j}(s^{*})|v_{n}(s^{*})\rangle|=0$. Using this and the expression for
$\lvert v_{n}(s^{*})\rangle$ in Eq. (65) we have that $\langle
v_{k}(s^{*})|U\rangle=-\langle v_{k}(s^{*})|M\rangle$. Substituting this for
$s=s^{*}$, we can conclude that

$\displaystyle HT^{+}(P,M)$
$\displaystyle\leq\dfrac{4}{\Delta(s^{*})}\sum_{k=1}^{n-1}|\langle
v_{k}(s^{*})|M\rangle|^{2}$ (74)
$\displaystyle\leq\dfrac{4}{\Delta(s^{*})}\left(1-|\langle
v_{n}(s^{*})|M\rangle|^{2}\right)$ (75)
$\displaystyle\leq\dfrac{2}{\Delta(s^{*})}$ (76)
$\displaystyle\implies\dfrac{1}{\Delta(s^{*})}$
$\displaystyle\geq\dfrac{HT^{+}(P,M)}{2}.$ (77)

Thus, we obtain that

$T=\Omega\left(\sqrt{HT^{+}(P,M)}\log(1/\epsilon)\right).$

$\Box$
Algorithm 2 can be thought of as an analog version of the quantum spatial
search algorithm by Krovi et al krovi2016quantum . Although this algorithm
requires $m(s^{*})$ additonal ancillary qubits unlike the algorithm of Ref.
chakraborty2018finding , it will help create an intuitive understanding of our
algorithm for solving $QSSamp$, discussed in Sec. VI. In fact, the cost of
preparing the initial state of the spatial search algorithm corresponds to the
$QSSamp$ problem.

In Ref. chakraborty2018finding , it was shown that $H(s)$ can be simulated be
using only query access to the discrete-time quantum walk unitary $W(s)$
introduced in Ref. krovi2016quantum . This connection would allow us to
quantify the running time of our quantum algorithm in terms of basic Markov
chain operations.

To that end, given a Markov chain $P$, let us define the following oracular
operations:

  * •

Check ($M$): Cost of checking whether a given node is marked. We denote this
by $\mathcal{C}$.

  * •

Update ($P$): Cost of applying one step of the walk $P$, which we denote by
$\mathcal{U}$.

  * •

Setup ($P$): The cost of preparing the initial state $\lvert v_{n}(0)\rangle$,
denoted by $\mathcal{S}$.

As from Refs. krovi2016quantum ; ambainis2019quadratic , the cost of
implementing $W(s)$, and consequently, the cost of evolving $H(s)$ for
constant time, is in $\mathcal{O}(\mathcal{C}+\mathcal{U})$, the running time
of Algorithm 2 is

$T=\mathcal{O}\left(\mathcal{S}+\sqrt{HT^{+}(P,M)}\left(\mathcal{U}+\mathcal{C}\right)\right)$
(78)

The $QSSamp$ problem helps quantify the cost $\mathcal{S}$ and intuitively, a
quantum algorithm for this problem can be obtained by running the spatial
search algorithm in reverse.

##  VI Analog quantum algorithm to prepare coherent encoding of the stationary
state of a Markov chain

In this section we describe our algorithm which, given a reversible Markov
chain $P$ with stationary state $\pi=(\pi_{1},\cdots,\pi_{n})$, prepares a
state that is $\epsilon$-close to the state

$\lvert\pi\rangle=\sum_{x\in X}\sqrt{\pi_{x}}\lvert x\rangle.$ (79)

A measurement in the basis spanned by the states of the Markov chain will
allow us to sample from $\pi$, thereby solving the $QSSamp$ problem. From Fact
2 and Proposition 3, we have that

$\lvert\pi\rangle=\lvert v_{n}(0)\rangle.$ (80)

Thus, this is simply the highest eigenstate of the discriminant matrix $D(P)$
or equivalently, the $0$-eigenstate of $H(0)$. Therefore, given $P$, the
problem of preparing $\lvert\pi\rangle$ boils down to the state-generation
problem just as in the case of spatial search.

Following Lemma 4 and Corollary 5, one can think of an algorithm to prepare
$\lvert v_{n}(0)\rangle$ as follows.

Starting from some initial localized state $\lvert j,0\rangle$ where $(j\in
X)$, one can evolve according to the Hamiltonian $H(0)\otimes\widehat{p}$ for
a time that scales as $\widetilde{\mathcal{O}}\left(1/\sqrt{\Delta}\right)$ to
prepare $\lvert v_{n}(0)\rangle$ with probability $|\langle
v_{n}(0)|j\rangle|^{2}\geq\eta$. Then by using $\Theta(1/\sqrt{\eta})$-rounds
of (fixed-point) amplitude amplification yoder2014fixed , one can prepare
$\lvert v_{n}(0)\rangle$.

However, amplitude amplification is a discrete quantum algorithm and to the
best of our knowledge it has no analog counterpart. As such, while
constructing an analog quantum algorithm for this problem we cannot make use
of amplitude amplification. We shall switch the value of $s$ to get around the
need for amplitude amplification.

Consider the scenario where, given $P$, one marks a single state $j$, i.e. all
the outgoing edges from $j$ are replaced with self-loops. We denote the
absorbing Markov chain corresponding to this $P^{\prime}_{j}$. Then the
resulting interpolated Markov chain is

$P(s)=(1-s)P+sP_{j}^{\prime}.$ (81)

If the entry of the stationary state of $P$, corresponding to the marked
element is $\pi_{j}$, then we find that $p_{M}=\pi_{j}$ and so for

$s=s^{*}=1-\pi_{j}/(1-\pi_{j}),$ (82)

and from Eq. (13) we have that

$\lvert
v_{n}(s^{*})\rangle=\dfrac{1}{\sqrt{2}}\left(\dfrac{1}{\sqrt{1-\pi_{j}}}\sum_{x\neq
j}\sqrt{\pi_{x}}\lvert x\rangle+\lvert j\rangle\right).$ (83)

Thus, the state $\lvert j\rangle$ has a constant overlap with $\lvert
v_{n}(s^{*})\rangle$. Also observe that the initial state of Algorithm 2
contained $\lvert v_{n}(0)\rangle$ in the first register and our state-
generation scheme resulted in the preparation of a state that has a constant
overlap of $|\langle v_{n}(s^{*})|v_{n}(0)\rangle|=\Theta(1)$ with $\lvert
v_{n}(s^{*})\rangle$ .

Let $\epsilon\in(0,1)$ and $\tau(s)=2\pi/\Delta(s)$.

  * 1\.

Set $s=s^{*}=1-\pi_{j}/(1-\pi_{j})$:

    * (a)

Evolve according to $H(s^{*})\otimes\widehat{p}$ for time $\tau(s^{*})$
starting from the state $\lvert j,0\rangle\lvert
x=0\rangle^{\otimes\lceil\log_{2}(4/\epsilon)\rceil}$.

    * (b)

Repeat the Hamiltonian evolution in step (a)
$\lceil\log_{2}(4/\epsilon)\rceil$ times, using a fresh block of
$l(s^{*})$-qubits in the

pointer register each time.

* (c)

Post-select on obtaining $\lvert 0\rangle^{m(s^{*})}$ in the pointer register

Let the state obtained after step 1 be $\lvert\psi^{(1)}_{f}\rangle$. 2\.
Reinitialize the pointer register. 3\.  Set $s=0$: (a) Repeat steps 1(a)-1(c),
starting from the state $\lvert\psi^{(1)}_{f}\rangle\lvert
x=0\rangle^{\lceil\log_{2}(4/\epsilon)\rceil}$ 4\.  Output the state of the
first register.

Algorithm 3 Quantum algorithm to the prepare stationary state of any
reversible Markov chain

For our algorithm, we assume that for any $j\in X$, the state $\lvert
j,0\rangle$ is easy to prepare. The idea of the algorithm (See Algorithm 3) is
to invoke Lemma 4 and Corollary 5 twice. At the first stage, we set
$s=s^{*}=1-\pi_{j}/(1-\pi_{j})$ and then starting from the state $\lvert
j,0\rangle$, we prepare a state that is close to $\lvert v_{n}(s^{*})\rangle$.
At the second stage, we set $s=0$ and starting from the state obtained in
stage 1, we prepare the state $\lvert v_{n}(0)\rangle=\lvert\pi\rangle$. By
this two stage procedure, we can avoid the need to use amplitude
amplification. We formally state the algorithm in Algorithm 3 and prove its
correctness in Lemma 7.

###### Lemma 7

Algorithm 3 outputs a quantum state $\lvert\phi_{f}\rangle$ such that

$\|\lvert\phi_{f}\rangle-\lvert\pi\rangle\|\leq\mathcal{O}(\epsilon),$

in time

$T=\Theta\left(\dfrac{1}{\sqrt{\Delta(s^{*})}}+\dfrac{1}{\sqrt{\Delta}}\right).$

Proof.

First note that the 0-eigenstate of $H(s)$ is given by Eq. (13) and so for
$s^{*}=1-\pi_{j}/(1-\pi_{j})$, we have that

$\lvert
v_{n}(s^{*})\rangle=\dfrac{1}{\sqrt{2}}\left(\dfrac{1}{\sqrt{1-\pi_{j}}}\sum_{x\neq
j}\sqrt{\pi_{x}}\lvert x\rangle+\lvert j\rangle\right)$ (84)

Thus, on starting from the state $\lvert j,0\rangle$ we have that
$\alpha_{n}=1/\sqrt{2}$. Following Lemma 4 and Corollary 5, this implies that
at the end of step 1, we will prepare a state $\lvert\psi^{(1)}_{f}\rangle$
that is $\epsilon/2$-close to $\lvert v_{n}(s^{*})\rangle$ in a time
$\tau(s^{*}).\lceil\log_{2}(4/\epsilon)\rceil$.

Note that $\langle v_{n}(s^{*})|v_{n}(0)\rangle\geq 1/\sqrt{2}$. So for this
second stage, $s=0$ and $\alpha_{n}\geq 1/\sqrt{2}$. The total time taken in
the second stage is $\tau(0)\lceil\log_{2}(4/\epsilon)\rceil$ and we output a
state that is $\epsilon$-close to $\lvert v_{n}(0)\rangle=\lvert\pi\rangle$.

Clearly, the total time taken is

$\displaystyle T$
$\displaystyle=\Theta\left(\left(\tau(s^{*})+\tau(0)\right)\log(1/\epsilon)\right)$
(85)
$\displaystyle=\Theta\left(\left(\dfrac{1}{\sqrt{\Delta(s^{*})}}+\dfrac{1}{\sqrt{\Delta}}\right)\log(1/\epsilon)\right).$
(86)

For a single marked marked node, $HT^{+}(P,M)$ is the same as the average
hitting time $HT(P,\\{j\\})$. So from Eq. (77), we have that
$1/\Delta(s^{*})=\Theta(HT(P,\\{j\\}))$. From Eq. (22), we have that the
classical mixing time is in $\Theta(1/\sqrt{\Delta})$. These two facts imply
that the running time of our algorithm is actually the sum of the square root
of the classical hitting time and the square root of the classical mixing
time, i.e.

$T=\widetilde{\Theta}\left(\sqrt{HT(P,\\{j\\})}+\sqrt{T_{\mathrm{mix}}}\right).$

$\Box$

Note that in general, the hitting time is at least as large as the mixing time
of an ergodic, reversible Markov chain. Thus the running time is in fact,

$T=\Theta\left(\sqrt{HT(P,\\{j\\})}\log(1/\epsilon)\right).$ (87)

Furthermore, as mentioned in the previous section, the $QSSamp$ problem helps
quantify the Setup cost $\mathcal{S}$ of the spatial search problem (See Eq.
(78)). As such, Algorithm 3 implies that the setup cost of Algorithm 2 is
given by Eq. (87).

##  VII Time-averaged quantum mixing: Limiting distribution and mixing time

Now we shall deal with the $QLSamp$ problem and the notion of mixing time that
arises from this problem. Recently in Ref. chakraborty2020fast , we considered
this problem for simple, unweighted graphs and proved an upper bound on the
quantum mixing time for almost all graphs. In this section, we explain the
mathematical insights that are key to obtaining the aforementioned result and
provide numerical evidence to back up the analytical results. Furthermore, we
generalize the $QLSamp$ problem to any ergodic, reversible Markov chain.

For any ergodic, reversible Markov chain $P$, we have seen from Sec. II.2 that
it is possible to sample from its distribution at $T\rightarrow\infty$
(limiting distribution) after a time
$T_{\mathrm{mix}}=\widetilde{\mathcal{O}}\left(1/\Delta\right)$, known as the
mixing time of $P$, where $\Delta$ is the spectral gap of $P$. In fact, any
initial distribution converges to the stationary distribution $\pi$ after
$T_{\mathrm{mix}}$ applications of $P$. In a strict sense, such a limiting
distribution is absent for quantum walks as the underlying dynamics is unitary
and hence, distance preserving.

However, one can define the mixing of quantum walks on a graph is defined in a
time-averaged sense, i.e. one can obtain the probability that the walker is at
some node $f$ after a time $t$ which is picked uniformly at random in the
interval $[0,T]$ aharonov2001quantum . This gives a time-averaged probability
distribution at any time $t$ and also a limiting probability distribution as
$T\rightarrow\infty$. The mixing time of a quantum walk on any ergodic,
reversible Markov chain $P$ is the time after which the time-averaged
probability distribution is $\epsilon$-close to the limiting probability
distribution.

Consider any ergodic, reversible Markov chain $P$ with $|X|=n$. Given $P$,
suppose $H_{P}$ denotes the underlying Hamiltonian corresponding to a quantum
walk on $P$. We require that the eigenvalues of $H_{P}$ lie between $-1$ and
$1$, i.e. $\left\lVert H_{P}\right\rVert=1$. Let the spectral decomposition of
$H_{P}=\sum_{i}\lambda_{i}\lvert v_{i}\rangle\langle v_{i}\rvert$ where
$\lvert v_{i}\rangle$ is the eigenstate corresponding to the eigenvalue
$\lambda_{i}$, $i\in\\{1,2,...,n\\}$. Furthermore, suppose that the initial
state of the walker is $\lvert\psi_{0}\rangle$.

Consequently, the state of the walker after a time $t$ is governed by the
Schrödinger equation, i.e.

$\lvert\psi(t)\rangle=e^{-iH_{P}t}\lvert\psi_{0}\rangle.$ (88)

In order to define a limiting distribution for quantum walks, one obtains a
Césaro-average of the probability distribution, i.e. one evolves for a time
$t$ chosen uniformly at random between $0$ and $T$ followed by a measurement.
The average probability that the state of the walker is some localized node
$\lvert f\rangle$ is given by

$P_{f}(T)=\frac{1}{T}\int_{0}^{T}dt|\langle
f|e^{-iH_{P}t}|\psi_{0}\rangle|^{2}.$ (89)

Thus, as $T\rightarrow\infty$, this leads to a limiting probability
distribution, i.e.

$P_{f}(T\rightarrow\infty)=\lim_{T\rightarrow\infty}P_{f}(T)=\sum_{\lambda_{i}=\lambda_{l}}\langle
v_{l}|f\rangle\langle f|v_{i}\rangle\langle
v_{i}|\psi_{0}\rangle\langle\psi_{0}|v_{l}\rangle,$ (90)

where the sum is over all pairs of degenerate eigenvalues. So, if $H_{P}$ has
a simple spectrum, i.e. all its eigenvalues are distinct, then the sum is over
all its eigenvalues.

In order to calculate how fast the instantaneous time-averaged distribution of
the quantum walk converges to this limiting distribution, i.e. we need to
bound the quantity $\left\lVert
P_{f}(T\rightarrow\infty)-P_{f}(T)\right\rVert_{1}$.

In fact, it is easy to verify that they are $\epsilon$-close, i.e.

$\left\lVert P_{f}(T\rightarrow\infty)-P_{f}(T)\right\rVert_{1}\leq\epsilon,$

as long as

$\displaystyle
T=\Omega\left(\dfrac{1}{\epsilon}\sum_{\lambda_{i}\neq\lambda_{l}}\dfrac{\left|\langle
v_{i}|\psi_{0}\rangle\right|.\left|\langle\psi_{0}|v_{l}\rangle\right|}{\left|\lambda_{l}-\lambda_{i}\right|}\right).$
(91)

This naturally leads to the following upper bound on the quantum mixing time

$T_{\mathrm{mix}}=\mathcal{O}\left(\dfrac{1}{\epsilon}\sum_{\lambda_{i}\neq\lambda_{l}}\dfrac{\left|\langle
v_{i}|\psi_{0}\rangle\right|.\left|\langle\psi_{0}|v_{l}\rangle\right|}{\left|\lambda_{l}-\lambda_{i}\right|}\right).$
(92)

There do exist differences between the quantum and classical limiting
distributions. For example, in the quantum case, the limiting distribution is
dependent on the initial state of the quantum walk. Also, unlike classical
random walks, the quantum mixing time depends on all the eigenvalue gaps of
$H_{P}$ as opposed to only the spectral gap.

Ignoring the numerator in the right hand side of Eq. (92), we need to evaluate
the following quantity in order to upper-bound $T_{\mathrm{mix}}$:

$\Sigma=\sum_{\lambda_{i}\neq\lambda_{l}}\dfrac{1}{\left|\lambda_{l}-\lambda_{i}\right|}.$
(93)

As such we intend to obtain the best possible bounds for this quantity. To
this end, let us define $\Delta_{\min}$ as the minimum eigenvalue gap of
$H_{P}$, over all pairs of distinct eigenvalues, i.e.

$\Delta_{\min}=\min_{\\\
i,j,\lambda_{i}\neq\lambda_{j}}\left\\{|\lambda_{i}-\lambda_{j}|,s.t.~{}i\neq
j\right\\}.$ (94)

Note that this is different from the spectral gap $\Delta$, which is the
difference between the two highest eigenvalues of $H_{P}$. We prove the
following:

###### Lemma 8

If $\Sigma$ and $\Delta_{\min}$ are defined as in Eqs. (93) and (94) then

$\dfrac{1}{\Delta_{\min}}\leq\Sigma\leq\widetilde{\mathcal{O}}\left(\dfrac{n}{\Delta_{\min}}\right)$
(95)

Proof.  The lower bound is straightforward by noting that $\exists i,l$ such
that $|\lambda_{l}-\lambda_{i}|=\Delta_{\min}$.

For the upper bound, we have that

$\Sigma\leq\dfrac{1}{\Delta_{\min}}\sum_{l\neq
i}\dfrac{1}{|\lambda_{l}-\lambda_{i}|},$ (96)

where we have used the fact that for any $\lambda_{i}\neq\lambda_{l}$,
$|\lambda_{l}-\lambda_{i}|\geq|l-i|\Delta_{\min}$. This implies that if
$|l-i|=r$, we obtain

$\displaystyle\Sigma$
$\displaystyle\leq\dfrac{1}{\Delta_{\min}}\sum_{r=1}^{n-1}\sum_{l,i:|l-i|=r}\dfrac{1}{r}$
(97)
$\displaystyle\leq\dfrac{1}{\Delta_{\min}}(n-1)\left(1+\dfrac{1}{2}+\cdots+\dfrac{1}{n-1}\right)$
(98) $\displaystyle\leq\dfrac{n\log
n}{\Delta_{\min}}=\widetilde{\mathcal{O}}\left(\dfrac{n}{\Delta_{\min}}\right).$
(99)

$\Box$

The upper bound on $\Sigma$ obtained in Lemma 8 leads directly to an upper
bound on the quantum mixing time $T_{\mathrm{mix}}$. This can be seen from the
fact that the numerator in Eq. (92) is less than one and so

$T_{\mathrm{mix}}=\widetilde{\mathcal{O}}\left(n/\Delta_{\min}\right),$ (100)

for any $H_{P}$.

All prior works hitherto have analyzed the $QLSamp$ problem for simple,
unweighted graphs. Given a graph $G(V,E)$ of $|V|=n$ nodes and $|E|$ edges,
the underlying quantum walk is defined on the nodes of the graph with
corresponding Hamiltonian being the (normalized) adjacency of the graph. That
is, $H_{P}=A_{G}/\left\lVert A_{G}\right\rVert$, where $A_{G}$ is an $n\times
n$ matrix such that each entry

$a_{ij}=\begin{cases}&1,~{}~{}~{}(i,j)\in E\\\
&0,~{}~{}~{}\text{otherwise},\end{cases}$ (101)

and $\left\lVert A_{G}\right\rVert$ is the spectral norm of $A_{G}$.

In fact, as mentioned earlier, several prior works have been concerned with
upper bounding $T_{\mathrm{mix}}$ for quantum walks on specific graphs.
However, bounds on the mixing time of quantum walks on general graphs have
been absent.

Recently, we obtained an upper bound on the quantum mixing time for Erdös-
Renyi random graphs: graphs of $n$ nodes such that an edge between any two
nodes exist with probability $p$ irrespective of any other nodes, typically
denoted as $G(n,p)$ chakraborty2020fast . Our results implied a bound on the
quantum mixing time for almost all graphs, i.e. the fraction of graphs for
which it holds, goes to one as $n\rightarrow\infty$. Our proof techniques
involved using several state-of-the art results on random matrix theory, which
find independent applications in several areas of physics.

In this section, we elaborate on the results obtained therein and provide
numerical evidence to back up our analytical bounds. In particular, we focus
on the random matrix theory aspects of our proof, elaborating on the
underlying concepts. Finally, by defining a quantum walk on the edges as in
Sec. IV.1, we extend the notion of $QLSamp$ to any ergodic, reversible Markov
chain.

###  VII.1 Erdös-Renyi random graphs

Let us consider a graph $G$ with a set of vertices $V=\\{1,\dots,n\\}$. We
restrict ourselves to simple graphs, i.e. unweighted graphs which do not
contain self-loops or multiple edges connecting the same pair of vertices. The
maximum number of edges that a simple graph $G$ can have is $N={n\choose 2}$.
Thus, there are ${N\choose M}$ graphs of $M$ edges and the total number of
(labelled) graphs is $\sum_{M=0}^{N}{N\choose M}=2^{N}$ graph_enumeration . We
consider the random graph model $G(n,p)$, a graph with $n$ vertices where we
have an edge between any two vertices with probability $p$, independently of
all the other edges ER59 ; ER60 ; bollobas_book  (See Fig. 1). In this model,
a graph $G_{0}$ with $M$ edges appears with probability
$P\\{G(n,p)=G_{0}\\}=p^{M}(1-p)^{N-M}$. In particular, if we consider the case
$p=1/2$, each of the $2^{N}$ graphs appears with equal probability $P=2^{-N}$.
We shall refer to random graphs having a constant $p$ as a dense random graph.
On the other hand, random graphs for which $p=o(1)$, i.e. when $p$ decreases
with $n$ shall be referred to as sparse random graphs.

Figure 1: Erdös-Renyi random graph $G(50,0.2)$.

In their seminal papers, Erdös and Rényi introduced this model of random
graphs and studied the probability of a random graph to possess a certain
property $Q$ ER59 ; ER60 . For example, they investigated properties such as
the connectedness of the graph, the probability that a certain subgraph is
present, etc. They stated that _almost all graphs_ have a property $Q$ if the
probability that a random graph $G(n,p)$ has $Q$ goes to $1$ as
$n\rightarrow\infty$. Equivalently, it can be stated that $G(n,p)$ _almost
surely_ has property $Q$, i.e. property $Q$ holds with probability $1-o(1)$.

Interestingly, certain properties of random graphs arise suddenly for a
certain critical probability $p=p_{c}$, where this probability depends
typically on $n$. More precisely, if $p(n)$ grows faster than $p_{c}(n)$, the
probability that the random graph has property $Q$ goes to $1$ in the
asymptotic limit, whereas if it grows slower than $p_{c}(n)$ it goes to $0$.
For example above the percolation threshold, i.e. when $p>\log(n)/n$ the graph
is almost surely connected, whereas if $p<\log(n)/n$ the graph has almost
surely isolated nodes.

Here we shall concern ourselves with random graphs above the percolation
threshold and calculate an upper bound on the quantum mixing time for quantum
walks on such graphs. Observe that for a random graph, $G(n,p)$, its adjacency
matrix, which we denote as $A_{G(n,p)}$, is an $n\times n$ symmetric matrix
with each non-diagonal entry being $1$ with probability $p$ and $0$ with
probability $1-p$. All diagonal entries of $A_{G(n,p)}$ are $0$. Thus
$A_{G(n,p)}$ is a discrete random matrix and knowledge of its eigenvalues and
eigenvectors is crucial to obtaining the quantum mixing time.

Finally, from the aforementioned discussion, obtaining the quantum mixing time
for $G(n,p)$ can be interpreted as holding for almost all graphs.

###  VII.2 Random matrices: Spectral properties of $\mathbf{A_{G(n,p)}}$

Here we look at the eigenvalues and eigenvectors of the random matrix
$A_{G(n,p)}$.

As mentioned earlier, the Hamiltonian corresponding to the quantum walk on
$G(n,p)$ is simply the normalized Adjacency matrix of $G(n,p)$. The highest
eigenvalue of $A_{G(n,p)}$, converges to a Gaussian distribution with mean
$np$ and standard deviation $\sqrt{p(1-p)}$, as $n\rightarrow\infty$. This
fact was first shown in Ref. furedi1981eigenvalues  for constant $p$ and was
later improved for sparse random graphs ($p=o(1)$) in Ref. erdHos2013spectral
. In fact, as we shall show shortly it suffices to consider the matrix

$\bar{A}_{G(n,p)}=\dfrac{A_{G(n,p)}}{np},$ (102)

as the quantum walk Hamiltonian.

Let the eigenvalues of $\bar{A}_{G(n,p)}$ be
$\lambda_{n}>\lambda_{n-1}\geq\cdots\lambda_{1}$, such that $\lvert
v_{i}\rangle$ is the eigenvector corresponding to the eigenvalue
$\lambda_{i}$, $i\in\\{1,2,...,n\\}$, i.e. $\bar{A}_{G(n,p)}\lvert
v_{i}\rangle=\lambda_{i}\lvert v_{i}\rangle$. Then we have that for
$p\geq\log^{8}(n)/n$,

$\lambda_{n}=1+\sqrt{\dfrac{1-p}{np}}o(1)+o\left(\dfrac{1}{n\sqrt{p}}\right),$
(103)

with probability $1-o(1)$, which implies that
$\left\lVert\bar{A}_{G(n,p)}\right\rVert\approx 1$ erdHos2013spectral ;
chakraborty2020fast .

It can also be shown that for the same range of $p$, the second highest
eigenvalue $\lambda_{n-1}$ can be upper bounded as

$\lambda_{n-1}\leq\dfrac{6}{\sqrt{np}}+\mathcal{O}\left(\dfrac{\log(n)}{(np)^{3/4}}\right),$
(104)

with probability $1-o(1)$ furedi1981eigenvalues ; vu2007 ; chakraborty2020fast
. This immediately implies that the spectral gap of $\bar{A_{G(n,p)}}$,
$\Delta=\widetilde{\mathcal{O}}\left(1\right)$. Consequently, a classical
random walk on $G(n,p)$ mixes quite fast - in
$\widetilde{\mathcal{O}}\left(1\right)$ time.

However, it is clear from the expression for $T^{G(n,p)}_{\mathrm{mix}}$ in
Eq. (92) that the knowledge of all eigenvalue gaps are crucial in obtaining
the quantum mixing time. As such require the knowledge of the spacings between
all the eigenvalues of the random matrix $\bar{A}_{G(n,p)}$.

Semicircle law:  It is well known that as $np\rightarrow\infty$, the spectral
density of the bulk of the spectrum of $A_{G(n,p)}$ converges to the well
known semicircle distribution given by

$\rho_{sc}(\lambda)=\begin{cases}\dfrac{\sqrt{4np(1-p)-\lambda^{2}}}{2\pi
np(1-p)}&\mbox{if }|\lambda|<2\sqrt{np(1-p)}\\\ 0&\mbox{otherwise
}\end{cases}.$ (105)

This implies that $\Theta(n)$ eigenvalues of $A_{G(n,p)}$ lie within $[-R,R]$
where

$R=2\sqrt{np(1-p)},$ (106)

is the radius of the semicircle. On applying the appropriate normalization, we
find that the spectral density of $\bar{A}_{G(n,p)}$ converges to a semicircle
of radius

$\bar{R}=2\sqrt{\dfrac{1-p}{np}}.$

The fraction of eigenvalues of $\bar{A}_{G(n,p)}$ lying in some spectral
window $\mathcal{I}\in[-R,R]$ converges to the area of the semicircle within
$\mathcal{I}$ as $np\rightarrow\infty$. However, the semicircle law provides
only a macroscopic description of the eigenvalues of $\bar{A}_{G(n,p)}$, i.e.
the aforementioned result holds only when $|\mathcal{I}|\gg 1$. However, in
order to obtain the quantum mixing time, we need information about all
eigenvalue gaps including consecutive gaps where $|\mathcal{I}|\sim 1/n$,
which renders the semicircle law useless. As a result, for our purposes we
need to look at mesoscopic and microscopic statistics of eigenvalues of
$\bar{A}_{G(n,p)}$.

However, for subsequent analysis, we shall require two results that can be
obtained from the semicircle law itself which we state now. Note that there
are $\Theta(n)$ eigenvalues with a radius of $\bar{R}$. This directly gives
information about the average eigenvalue gap of $\bar{A}_{G(n,p)}$ given by

$\bar{\Delta}=\Theta\left(\dfrac{1}{n^{3/2}\sqrt{p}}\right).$ (107)

Also from the semicircle law itself, one can define the so-called classical
eigenvalue locations of $\bar{A}_{G(n,p)}$. For each $1\leq i\leq n-1$, we can
define the classical location $\gamma_{i}$ as the solution to the following
equation

$\int_{-\infty}^{\gamma_{i}}\rho_{sc}(x)\,dx=\frac{i}{n}.$ (108)

Thus, the position of the $i^{\mathrm{th}}$ classical location is obtained by
filling up $i/n$-area of the semicircle. From this condition one obtains that
for $i\leq n/2$, $r\leq n-2i$ and some universal constant $c>0$,

$\gamma_{i+r}-\gamma_{i}\geq c\dfrac{r}{n^{7/6}i^{1/3}\sqrt{p}}.$ (109)

An identical estimate holds for the other half of the spectrum by symmetry.

Eigenvalue rigidity criterion:  The semicircle law was shown to hold for
smaller spectral windows in Refs. erdHos2012rigidity ; erdHos2013spectral . An
immediate consequence of this fact is that the every eigenvalue (with the
exception of $\lambda_{n}$) of $\bar{A}_{G(n,p)}$ is located close to their
classical eigenvalue positions. Formally, they showed that for $n^{-1/3}\leq
p\leq 1-n^{-1/3}$ and any $\varepsilon\geq 0$ the eigenvalues of
$bar{A}_{G(n,p)}$ satisfy

$|\lambda_{i}-\gamma_{i}|\leq\frac{n^{\varepsilon}(n^{-2/3}\alpha_{i}^{-1/3}+n^{-1-\phi})}{(pn)^{1/2}}$
(110)

with probability $1-o(1)$, where

$\phi:=\frac{\log p}{\log n}\quad\text{and}\quad\alpha_{i}:=\max\\{i,n-i\\}.$

Figure 2: A pictorial representation of eigenvalue rigidity: The actual
eigenvalue locations of adjacency matrices of Erdös-Renyi random graphs,
$\lambda_{i}$ (denoted by black strokes) are close to their classical
eigenvalue positions, $\gamma_{i}$ (denoted by red strokes) as predicted by
the semicircle law. Although eigenvalue rigidity provides information about
eigenvalue value gaps that are far away, it does not provide any information
about the smallest eigenvalue gaps (such as consecutive eigenvalue gaps). As a
result the eigenvalue rigidity criterion provides information about
eigenvalues at a mesoscopic scale and information about the smallest gaps are
obtained from eigenvalue statistics at a microscopic scale.

Eigenvalue rigidity does reveal information about eigenvalue gaps of
$\bar{A}_{G(n,p)}$. Note that for any $r\geq 1$ one obtains

$\lambda_{i+r}-\lambda_{r}=\left(\lambda_{i+r}-\gamma_{i+r}\right)+\left(\gamma_{r}-\lambda_{r}\right)+\left(\gamma_{i+r}-\gamma_{r}\right).$
(111)

As a result, whenever $|\gamma_{i+r}-\gamma_{r}|$ scales larger than
$\left|\lambda_{i+r}-\gamma_{i+r}\right|+\left|\lambda_{r}-\gamma_{r}\right|$,
eigenvalue rigidity kicks in and an accurate estimate of
$\lambda_{i+r}-\lambda_{i}$ is given by the difference between their classical
eigenvalue locations, $\gamma_{i+r}-\gamma_{i}$. That is, there exists some
$r=r^{\star}(i)$, such that for all $r\geq r^{\star}(i)$,

$\lambda_{i+r}-\lambda_{r}\approx\left(\gamma_{i+r}-\gamma_{r}\right).$

From Eq. (109) and Eq. (110), we obtain that

$r^{\star}(i)=n^{\varepsilon}\max\\{1,n^{2/3}\alpha_{i}^{1/3}n^{-1-\phi}\\}\leq
n^{\varepsilon-\log p/\log n}.$ (112)

As such we cannot exploit eigenvalue rigidity to estimate gaps of the form
$|\lambda_{l}-\lambda_{i}|$ as long as $|l-i|\leq r^{\star}(i)$. Thus,
eigenvalue rigidity does not provide information about the smallest eigenvalue
gaps (for a pictorial representation of this fact, see Fig. 2) and reveals
eigenvalue statistics of $\bar{A}_{G(n,p)}$ at a mesoscopic scale. Thus in
order to obtain information about consecutive eigenvalue gaps of
$\bar{A}_{G(n,p)}$, we shall need to go to a microscopic scale.

Microscopic eigenvalue statistics of $\mathbf{\bar{A}_{G(n,p)}}$: At the
microscopic scale, results are notoriously difficult to obtain. Tao and Vu
tao2014random  showed that $\bar{A}_{G(n,p)}$ has a simple spectrum for dense
graphs, resolving a long-standing conjecture due to Babai babai1982isomorphism
. Recently, this was resolved also for sparse graphs luh2018sparse . We state
their results formally.

###### Fact 9

There exists a constant $C>0$ such that for $\frac{C\log^{6}(n)}{n}\leq p\leq
1-\frac{C\log^{6}(n)}{n}$, $\bar{A}_{G(n,p)}$ has a simple spectrum with
probability $1-o(1)$.

The fact that every eigenvalue gap is non-zero implies that the expression for
$\bar{A}_{G(n,p)}$, the double sum $\Sigma$ in Eq. (93) can now be re-written
as

$\displaystyle\Sigma$
$\displaystyle=\underbrace{\sum_{i=1}^{n-1}\dfrac{1}{\left|\lambda_{i+1}-\lambda_{i}\right|}}_{\Sigma_{1}}+\underbrace{\sum_{i=1}^{n-2}\dfrac{1}{\left|\lambda_{i+2}-\lambda_{i}\right|}}_{\Sigma_{2}}+\cdots$
(113)
$\displaystyle=\sum_{r=1}^{n-1}\Sigma_{r}=\sum_{r=1}^{n-1}\sum_{i=1}^{n-r}\dfrac{1}{\left|\lambda_{i+r}-\lambda_{i}\right|},$
(114)

while the limiting probability distribution is

$P_{f}(T\rightarrow\infty)=\lim_{T\rightarrow\infty}P_{f}(T)=\sum_{i=1}^{n}|\langle
f|v_{i}\rangle\langle v_{i}|\psi_{0}\rangle|^{2}.$ (115)

Furthermore, Nguyen, Tao and Vu nguyen2015random  proved that all the
eigenvalue gaps of $\bar{A}_{G(n,p)}$ for dense random graphs are not only
non-zero, but also separated. This was improved for the case of sparse random
graphs by Lopatto and Luh lopatto2019tail . In fact they asked the following
question: how likely is it for any eigenvalue gap
$\delta_{i}=\lambda_{i+1}-\lambda_{i}$ to be less than some $\delta$ times the
average gap $\bar{\Delta}$? They proved that there exists a constant $C>0$
such that for $n^{-1/3}\leq p\leq 1-n^{-1/3}$,

$\sup_{1\leq i\leq
n-1}\mathbb{P}\left(\delta_{i}\leq\frac{\delta}{n^{3/2}\sqrt{p}}\right)\leq
C\delta\log n,$ (116)

for all $\delta\geq n^{-C}$.
Applying a union bound to this gives a lower bound on the minimum eigenvalue
gap $\Delta_{\min}$ (defined in Eq. (94)) for $\bar{A}_{G(n,p)}$. We prove
that here.

###### Lemma 10

Lower bound on $\Delta_{\min}$: For $p\geq n^{-1/3}$,

$\Delta_{\min}\geq\dfrac{1}{n^{5/2+o(1)}\sqrt{p}}.$ (117)

with probability $1-o(1)$.

Proof.  Let $\mathcal{A}_{i}$ be the event that
$\delta_{i}\leq\frac{\delta}{n^{3/2}\sqrt{p}}$. Then using the union bound and
Eq. (116), we obtain

$\mathbb{P}\left(\bigcup_{i}\mathcal{A}_{i}\right)\leq\sum_{i}\mathbb{P}\left(\mathcal{A}_{i}\right)\leq
C\ n\ \delta\log n.$ (118)

This implies that the probability that at least one of the gaps is less than
$\frac{\delta}{n^{3/2}\sqrt{p}}$ is upper bounded by the right hand side of
Eq. (118). By choosing

$\delta=\dfrac{1}{n^{1+o(1)}},$

we have that

$\mathbb{P}\left(\bigcup_{i}\mathcal{A}_{i}\right)\leq o(1),$ (119)

i.e. with probability $1-o(1)$, no $\delta_{i}$ is less than
$\frac{1}{n^{5/2+o(1)}\sqrt{p}}$. This in turn implies that,

$\Delta_{\min}\geq\dfrac{1}{n^{5/2+o(1)}\sqrt{p}}.$ (120)

with probability $1-o(1)$. $\Box$

We are now equipped with the random matrix theory results, and in the next
subsection, we derive a tight upper bound on the double sum $\Sigma$, defined
in Eq. (114).

###  VII.3 Upper bound on $\mathbf{\Sigma}$

As mentioned previously, in order to obtain the quantum mixing time, we first
obtain bounds for the double sum $\Sigma$. Recall that we an obtain lower and
upper bounds for $\Sigma$ as

$\dfrac{1}{\Delta_{\min}}\leq\Sigma\leq\widetilde{\mathcal{O}}\left(\dfrac{n}{\Delta_{\min}}\right).$

In this subsection, our goal is to obtain an upper bound for $\Sigma$ that is
as close as possible to its lower bound. To that end, our strategy would be to
make use of the results on the eigenvalue statistics of $\bar{A}_{G(n,p)}$ at
macroscopic, mesoscopic and microscopic levels. In particular, in order to
evaluate $\lambda_{i+r}-\lambda_{i}$, for $r<r^{\star}(i)$, we shall make use
of the tail-bounds on consecutive eigenvalue gaps in Eq. (116). On the other
hand for $r>r^{\star}(i)$, the eigenvalue rigidity criterion (See Eq. (110))
kicks in and we can replace $\lambda_{i+r}-\lambda_{i}$ with
$\gamma_{i+r}-\gamma_{i}$.

Upper bound on $\mathbf{\Sigma_{1}}$:  We first obtain an upper bound on the
sum of the inverse of consecutive eigenvalue gaps, i.e.

$\Sigma_{1}=\sum_{i=1}^{n}\dfrac{1}{\lambda_{i+1}-\lambda_{i}}.$ (121)

In the Supplemental Material of Ref. chakraborty2020fast , we have explicitly
derived an upper bound for $\Sigma$. We restate the result here

######  Lemma 11 (Upper bound on $\mathbf{\Sigma_{1}}$ chakraborty2020fast )

$\Sigma_{1}=\sum_{i=1}^{n-1}\frac{1}{\lambda_{i+1}-\lambda_{i}}\leq
n^{5/2+o(1)}\sqrt{p},$ (122)

with probability $1-o(1)$.

The key idea is that we count the number of consecutive eigenvalue gaps
($\delta_{i}$) lying within an interval of $1/\log n$ times the average gap
and find that a high fraction of the $\delta_{i}$-s lie within this window
around the average gap. For a detailed derivation, we refer the readers to
Ref. chakraborty2020fast .
Now we can derive an upper bound on $\Sigma$ by combining mesoscopic and
microscopic eigenvalue statistics of $\bar{A}_{G(n,p)}$ at different scales of
$r$. In particular, we use the upper bound on $\Sigma_{1}$ along with the
eigenvalue rigidity condition. We state the upper bound on $\Sigma$ that we
obtained in Ref. chakraborty2020fast .

######  Lemma 12 (Upper bound on $\mathbf{\Sigma}$ chakraborty2020fast )

For $p\geq n^{-1/3}$, the eigenvalues of $\bar{A}_{G(n,p)}$ satisfy

$\Sigma=\sum_{i=1}^{n-1}\sum_{r=1}^{n-i}\dfrac{1}{\left|\lambda_{i+r}-\lambda_{i}\right|}\leq
n^{5/2-\frac{\log p}{\log n}+o(1)}\sqrt{p},$ (123)

with probability $1-o(1)$.

We provide an an intuition of the proof techniques and we refer the reader to
the Supplemental Material of Ref. chakraborty2020fast  for details. We first
split $\Sigma$ into two different parts:

$\displaystyle\Sigma=\sum_{r=1}^{r^{\star}(i)}\Sigma_{r}+\sum_{r=r^{\star}(i)+1}^{n-1}\Sigma_{r}$
(124)

For the first sum in the right hand side of Eq. (124), we are dealing with
small eigenvalue gaps and hence we make use of the microscopic eigenvalue
statistics, namely the upper bound on $\Sigma_{1}$, i.e. we replace this sum
with with the upper bound $r^{\star}(i).\Sigma_{1}$. On the other hand, for
the second double sum, eigenvalue rigidity provides kicks in and the gaps
between the classical eigenvalue locations $(\gamma_{i+r}-\gamma_{i})$, is a
better estimate of $\lambda_{i+r}-\lambda_{i}$ than the tail bounds. In fact,
an upper bound is obtained by replacing each eigenvalue gap
$\lambda_{i+r}-\lambda_{i}$ with the lower bound from Eq. (109).

Observe that for dense random graphs, the upper bound on $\Sigma$ is quite
close to its lower bound of $1/\Delta_{\min}$. Having obtained this bound, we
shall now upper bound the quantum mixing time for $G(n,p)$.

###  VII.4 Mixing of continuous-time quantum walks on $\mathbf{G(n,p)}$

Here, we shall obtain the (i) limiting distribution of the quantum walk and
the time after which the quantum walk converges (in a time-averaged sense) to
this distribution, the quantum mixing time.

Figure 3: The limiting probability distribution is close to the uniform
distribution for a quantum walk on $G(n,p)$. The figure shows that the
instantaneous time-averaged probability distribution (thick blue line) for a
quantum walk on $G(50,0.5)$ remains close to the uniform distribution
(horizontal dashed black) after a long enough time.

In order to obtain both these results, we make use of the fact that all the
eigenvectors of $\bar{A}_{G(n,p)}$ are completely delocalized. In fact, it was
conjectured in Ref. dekel2011eigenvectors , that for dense random graphs, the
eigenstates of $\bar{A}_{G(n,p)}$ are completely delocalized. This implies
that when any of its eigenvectors $\lvert v_{i}\rangle$ is expressed in the
basis of the nodes of the underlying graph, the absolute value of each entry
is at most $n^{-1/2}$ (up to logarithmic factors). Erdös et al.
erdHos2013spectral  answered this optimally even for sparse $p$ and the
results therein were subsequently extended for any $p$ above the percolation
threshold recently by He et al. he2018local . This implies that as long as
$p\geq\omega(\log(n)/n)$, for all $j\in\\{1,\cdots,n\\}$

$\|\lvert v_{j}\rangle\|_{\infty}\leq n^{-1/2+o(1)},$ (125)

with probability $1-o\left(\dfrac{1}{n}\right)$.

To the expression for the limiting probability distribution in Eq. (115), we
substitute the delocalization of eigenvectors from Eq. (125) to obtain

$\displaystyle P_{f}(T\rightarrow\infty)$
$\displaystyle=\sum_{i=1}^{n}|\langle f|v_{i}\rangle\langle
v_{i}|\psi_{0}\rangle|^{2}$ (126)
$\displaystyle\leq\widetilde{\mathcal{O}}\left(1/n\right)\sum_{i=1}^{n}|\langle
v_{i}|\psi_{0}\rangle|^{2}$ (127)
$\displaystyle\leq\widetilde{\mathcal{O}}\left(1/n\right),$ (128)

independent of $\lvert\psi_{0}\rangle$, i.e. the limiting distribution
converges to a (nearly) uniform distribution.

Observe that the upper bound on $\Sigma$ already provides an upper bound on
the quantum mixing time. However, we can improve the bound it further if we
assume that the quantum walk commences from an easy to prepare initial state.
By this we mean that the initial state $\lvert\psi_{0}\rangle$ is a
superposition over atmost a $\mathrm{polylog}(n)$ number of nodes. In fact,
generally it is assumed that the initial state is localized at some node of
the underlying graph, i.e. $\lvert\psi_{0}\rangle=\lvert l\rangle$, which is
standard.

Figure 4: Figure shows the time for the instantaneous time-averaged
probability distribution at any time $t$, denoted by $P_{f}(t)$ to be
$\epsilon$-close to the limiting probability distribution,
$P_{f}(t\rightarrow\infty)$ for Erdös-Renyi random graphs $G(n,p)$. The Y-axis
denotes the distance between these two distributions (as measured in one
norm), i.e. $D_{P}(t)=\left\lVert
P_{f}(t)-P_{f}(t\rightarrow\infty)\right\rVert_{1}$, while the X-axis denotes
time. We plot $D_{P}(t)$ for random graphs of $40$ nodes (solid thin green),
$60$ nodes (dot-dashed blue), $80$ nodes (thick, red) and $100$ nodes (dashed
pink), with $p=0.5$. The dotted horizontal line (dashed black) corresponds to
$\epsilon=0.1$ which helps indicate the time after which
$D_{P}(t)\leq\epsilon$ for the aforementioned instances. The inset plot shows
the exponent $c$ where $n^{c}$ corresponds to the minimum time after which
$D_{P}(t)\leq 0.1$ (quantum mixing time) for
$G(10,0.5),G(20,0.5),\cdots,G(100,0.5)$. The quantum mixing time is thus upper
bounded by $n^{3/2}$ which matches with our analytical predictions.

If the quantum walk commences from an easy to prepare state,

$\lvert\psi_{0}\rangle=\sum_{k=1}^{q}c_{k}\lvert k\rangle,$

where $q$ is in $\mathcal{O}(\mathrm{polylog}(n))$, we can use Eq. (125) to
obtain

$\displaystyle T^{G(n,p)}_{\mathrm{mix}}$
$\displaystyle=\mathcal{O}\left(\dfrac{1}{\epsilon}\sum_{i=1}^{n-1}\sum_{r=1}^{n-i}\dfrac{\left|\langle
v_{i}|\psi_{0}\rangle\right|.\left|\langle\psi_{0}|v_{i+r}\rangle\right|}{\left|\lambda_{i+r}-\lambda_{i}\right|}\right)$
(129)
$\displaystyle=\mathcal{O}\left(\dfrac{1}{n^{1-o(1)}\epsilon}\sum_{i=1}^{n-1}\sum_{r=1}^{n-i}\dfrac{\sum_{l=1}^{q}|c_{l}|.\sum_{m=1}^{q}|c^{*}_{m}|}{\left|\lambda_{i+r}-\lambda_{i}\right|}\right)$
(130) $\displaystyle=\widetilde{\mathcal{O}}\left(\frac{1}{\epsilon}\
\dfrac{\Sigma}{n}\right)=\widetilde{\mathcal{O}}\left(\frac{1}{\epsilon}\
n^{3/2-\frac{\log p}{\log n}+o(1)}\sqrt{p}\right),$ (131)

with probability $1-o(1)$.

Thus for $n^{-1/3}\leq p\leq 1-n^{-1/3}$,

$T^{G(n,p)}_{\mathrm{mix}}=\widetilde{\mathcal{O}}\left(n^{3/2-\log p/\log
n}\sqrt{p}/\epsilon\right),$ (132)

for $p\geq n^{-1/3}$.

Observe that for dense Erdös-Renyi random graphs,

$T^{G(n,p)}=\widetilde{\mathcal{O}}\left(\dfrac{n^{3/2}}{\epsilon}\right).$
(133)

Also, as $p$ decreases the upper bound on the mixing time increases.
Unfortunately for sparser random graphs, i.e. for $p=\log^{D}(n)/n$, such that
$D>8$, we cannot make use of eigenvalue rigidity. However simply using Lemma
11 along with the observation that

$\sum_{i=1}^{n-r-1}\dfrac{1}{\left|\lambda_{i+r}-\lambda_{i}\right|}\leq\sum_{i=1}^{n-1}\dfrac{1}{|\lambda_{i+1}-\lambda_{i}|},$

for $2\leq r\leq n-1$, gives us a weaker upper bound for the quantum mixing
time in such regimes of sparsity. We obtain that

$T^{G(n,p)}_{\mathrm{mix}}=\mathcal{O}\left(\dfrac{n^{5/2+o(1)}\sqrt{p}}{\epsilon}\right).$
(134)

In fact, the breakdown of rigidity estimates in erdHos2013spectral  is not an
artifact of the proof. For extremely sparse graphs, the optimal rigidity
estimates that hold in dense graphs are known to break down
huang2017transition .

Note that there exist weaker forms of rigidity of sparse graphs when $p\leq
n^{-1/3}$ which may lead to modest improvements of the exponent of $n$ in the
mixing time. However, we have not expended too much effort optimizing the
exponent as we are fundamentally limited by the smallest gap, $\Delta_{\min}$
(see lower bound of Eq. (120)) for which the bounds in lopatto2019tail  are
still quite far from the conjectured behaviour. Obtaining the conjectured
smallest gap behaviour appears to be a difficult problem in random matrix
theory.

Finally, we numerically verify the analytical results obtained. Fig. 3 shows
that for $G(50,0.5)$, the instantaneous time-averaged probability distribution
($P_{f}(t)$) converges to a distribution that is close to the uniform
distribution (horizontal dashed black line). While in Fig. 4 we plot
$D_{P}(t)=\left\lVert P_{f}(t)-P_{f}(t\rightarrow\infty)\right\rVert_{1}$ with
time and the inset plot depicts the exponent for the quantum mixing time
($D_{P}(t)\leq\epsilon$) for random graphs of various sizes and $p=0.5$. The
numerical results conform with the analytically obtained upper bound for the
quantum mixing time in Eq. (133).

###  VII.5 Mixing time for continuous-time quantum walks on any ergodic,
reversible Markov chain

Our results thus far have provided an upper bound on the quantum mixing time
for almost all simple unweighted graphs. Now we address the quantum mixing
time for any ergodic, reversible Markov chain. Given any such Markov chain
$P$, one can define the Hamiltonian $H=i[V^{\dagger}SV,\Pi_{0}]$ as stated in
Sec. IV (for $s=0$). In this section we consider the limiting distribution of
a continuous-time quantum walk under $H$, on the edges of $P$.

For any such $P$, one can also perform a quantum walk with the discriminant
matrix $D(P)$ as the Hamiltonian. For the quantum spatial search and the
$QSSamp$ problem, a quantum walk under $H$ offers a quadratic advantage over
$D(P)$. This is due to the fact that the spectral gap of $H$ is square root of
that of $P$ (or equivalently $D(P)$).

Here, we shall explore whether any generic speedup is obtained for the
$QLSamp$ problem. First we will compare the limiting distribution obtained
from performing a quantum walk on $H$ as compared to $D(P)$. The time
evolution of some initial state $\lvert\psi(0),0\rangle$, under the action of
$H$ is given by

$\displaystyle\lvert\psi(t)\rangle=\langle v_{n}|\psi_{0}\rangle\lvert
v_{n},0\rangle+\sum_{j=1,\sigma=\pm}^{n-1}e^{-itE^{\sigma}_{j}}\dfrac{\langle
v_{j}|\psi_{0}\rangle}{\sqrt{2}}\lvert\Psi^{\sigma}_{j}\rangle.$ (135)

The limiting probability distribution (note that now we are projecting on
obtaining $\lvert 0\rangle$ in the second register) is given by

$P_{f}(T)=\frac{1}{T}\int_{0}^{T}dt|\langle
f,0|e^{-iHt}|\psi_{0},0\rangle|^{2}.$

Figure 5: Comparison of the gaps between eigenvalues of a Markov chain $P$ and
the corresponding Hamiltonian $H$ defined in Sec. IV. The eigenvalues of $P$
lie between $0$ and $1$. Any such eigenvalue $\lambda$ of $P$ is mapped to the
eigenvalue pair $\pm\sqrt{1-\lambda^{2}}$ in the relevant subspace of $H$. As
a result the spectral gap, $\Delta$, of $P$ is mapped to
$\Theta(\sqrt{\Delta})$ for $H$. However, this is not the case for all
eigenvalue gaps. In fact, the minimum over all eigenvalue gaps of $P$,
$\Delta_{\min}$ is mapped to $\widetilde{\Delta}_{\min}$, such that
$\widetilde{\Delta}_{\min}>\Delta_{\min}$ if $\Delta_{\min}$ appears between
two eigenvalues that are close to $\lambda_{n-1}$. On the other hand,
$\widetilde{\Delta}_{\min}>\Delta_{\min}$ if it appears between two
eigenvalues that are close to $\lambda_{2}$. This has been elucidated in Sec.
VII.5.

This implies,

$P_{f}(T\rightarrow\infty)=\frac{1}{2}\sum_{\lambda_{i}=\lambda_{l}}\langle
v_{l}|f\rangle\langle f|v_{i}\rangle\langle
v_{i}|\psi_{0}\rangle\langle\psi_{0}|v_{l}\rangle,$ (136)

which is the same as the limiting distribution obtained for a quantum walk on
$D(P)$.

Also, the upper bound on the quantum mixing time is given by,

$\displaystyle T^{P}_{\mathrm{mix}}$
$\displaystyle=\mathcal{O}\left(\dfrac{1}{\epsilon}\sum_{i\neq
l}\dfrac{\left|\langle
v_{i}|\psi_{0}\rangle\right|.\left|\langle\psi_{0}|v_{l}\rangle\right|}{\left|E_{i}-E_{l}\right|}\right),$
(137)

where recall from Eq. (61) in Sec. IV.2 that $E_{j}=\sqrt{1-\lambda_{j}^{2}}$.

Now the generic upper bound on the quantum mixing time $T^{P}_{\mathrm{mix}}$
is upper bounded by the double sum $\Sigma$ and as such

$T^{P}_{\mathrm{mix}}\leq\Sigma=\sum_{i\neq
l}\dfrac{1}{\left|E_{i}-E_{l}\right|}\leq\widetilde{\mathcal{O}}\left(\dfrac{n}{\widetilde{\Delta}_{\min}}\right),$
(138)

where $\widetilde{\Delta}_{\min}$ is the minimum eigenvalue gap of the
Hamiltonian $H$. We now need bound $\widetilde{\Delta}_{\min}$ in terms of the
minimum eigenvalue gap of $D(P)$, $\Delta_{\min}$. To that end we have the
following lemma:

###### Lemma 13

Suppose $P$ is an ergodic, reversible Markov chain with eigenvalues
$\lambda_{n}=1>\lambda_{n-1}\geq\cdots\lambda_{1}\geq 0$. Suppose $\Delta$ is
the spectral gap of $P$ and the minimum of all gaps between distinct
eigenvalues of $P$ be $\Delta_{\min}$. Then the minimum eigenvalue gap of the
Hamiltonian $H=i[V^{\dagger}SV,\Pi_{0}]$, $\widetilde{\Delta}_{\min}$ is
bounded as

$\Theta\left(\lambda_{2}\Delta_{\min}\right)\leq\widetilde{\Delta}_{\min}\leq\Theta\left(\dfrac{\Delta_{\min}}{\sqrt{\Delta}}\right).$

Proof.  We know that for $H$, in the relevant subspace, each eigenvalue of
$P$, $\lambda_{j}$, maps to $\pm\sqrt{1-\lambda_{j}^{2}}$. Thus if
$\delta_{j}=|\lambda_{j+1}-\lambda_{j}|$, then we have

$\displaystyle\widetilde{\delta_{j}}$
$\displaystyle=\left|\sqrt{1-\lambda_{j+1}^{2}}-\sqrt{1-\lambda_{j}^{2}}\right|$
(139)
$\displaystyle=\left|\sqrt{1-\lambda_{j+1}^{2}}-\sqrt{1-(\lambda_{j+1}-\delta_{j})^{2}}\right|$
(140)
$\displaystyle=\sqrt{1-\lambda_{j+1}^{2}}\left|1-\sqrt{1+\dfrac{2\delta_{j}\lambda_{j+1}}{1-\lambda_{j+1}^{2}}-\dfrac{\delta_{j}^{2}}{1-\lambda_{j+1}^{2}}}\right|.$
(141)

We are concerned with the minimum eigenvalue gap $\widetilde{\Delta}_{\min}$.
Without loss of generality, we assume that $P$ has a simple spectrum
(consequently, so does $H$) and for some $1\leq j\leq n-1$, the eigenvalue gap
is minimum for two consecutive distinct eigenvalues $\lambda_{j}$ and
$\lambda_{j+1}$. That is, for some value of $j$, $\delta_{j}=\Delta_{\min}$
and henceforth we consider that value of $j$. Observe that in such a case, the
second term inside the square-root

$\dfrac{2\delta_{j}\lambda_{j+1}}{1-\lambda_{j+1}^{2}}=\dfrac{2\Delta_{\min}\lambda_{j+1}}{1-\lambda_{j+1}^{2}}<\dfrac{2\Delta_{\min}}{\Delta}<1.$
(142)

So expanding Eq. (141) according to Taylor series, we have

$\displaystyle\widetilde{\Delta}_{\min}$
$\displaystyle=\dfrac{2\Delta_{\min}\lambda_{j+1}}{\sqrt{1-\lambda^{2}_{j+1}}}+\Theta\left(\dfrac{\Delta^{2}_{\min}}{\sqrt{1-\lambda^{2}_{j+1}}}\right)$
(143)
$\displaystyle=\Theta\left(\dfrac{\Delta_{\min}\lambda_{j+1}}{\sqrt{1-\lambda^{2}_{j+1}}}\right).$
(144)

This expression implies that the minimum eigenvalue gap of $P$ is mapped to
the minimum eigenvalue gap of $H$ multiplied by the ratio of the corresponding
eigenvalues of $P$ and $H$. The upper and lower bounds follow from observing
that for all $1\leq j\leq n-1$,
$\sqrt{1-\lambda^{2}_{j+1}}\leq\Theta(\sqrt{\Delta})$ and
$\lambda_{j+1}=\Omega\left(\lambda_{2}\right)$, respectively. $\Box$

Assuming that the eigenvalues of $P$ (and hence of $D(P)$) are ordered, for a
continuous-time quantum walk on $H$, from Eq. (138) and Lemma 13, we observe
that the quantum mixing time may be faster or slower than a quantum walk
performed on $D(P)$ depending on where the minimum eigenvalue gap appears (See
Fig. 5 for a pictorial representation). If $\Delta_{\min}$ happens to be
between two eigenvalues that are close to $\lambda_{n-1}$, then the quantum
mixing time on $H$ is faster than that of $D(P)$. However, if $\Delta_{\min}$
is in the vicinity of $\lambda_{2}$, the quantum mixing time is slower. So
depending on the underlying Markov chain, it may happen that performing a
quantum walk on $H$ can be disadvantageous with respect to performing the walk
on $D(P)$ itself.

This is in contrast to the $QSSamp$ problem, where using $H$ offers a generic
quadratic speedup over using $D(P)$ as the Hamiltonian in Algorithm 3. This
shows a fundamental difference between the two different notions of mixing for
quantum algorithms as elucidated by $QSSamp$ and $QLSamp$ problems.

##  VIII Discussion

In this article we have discussed the two notions of quantum mixing and
designed analog quantum algorithms to tackle these problems. First, using
Hamiltonian evolution and von Neumann measurements, we have presented an
analog quantum algorithm that, given an ergodic, reversible Markov chain
outputs a coherent encoding of its stationary state. The running time of our
algorithm matches that of its discrete-time counterparts. Secondly, we have
also discussed the problem of sampling from the limiting distribution of a
(time-averaged) continuous-time quantum walk. We have offered an intuitive
explanation of the tools used in Ref. chakraborty2020fast  to derive upper
bounds on the mixing time for random graphs. We have also backed up the
analytical results therein with numerical simulations and extended the time-
averaged notion of mixing to any ergodic, reversible Markov chain.

Our results could pave the way for further research. For example, quantum
state-generation using von Neumann measurements can be used to develop novel
analog quantum algorithms. Note that our methods could be used to obtain other
analog quantum algorithms for solving the $QSSamp$ problem. One could reverse
the spatial search algorithm by Childs and Goldstone childs2004spatial ;
chakraborty2020optimality  and use von Neumann measurements to prepare a
coherent encoding of the highest eigenstate of underlying Hamiltonian. In the
case of state-transitive graphs, this will allow for uniform sampling.

It would be interesting to explore whether using our framework, one can
construct an analog quantum algorithm to fast-forward the dynamics of any
ergodic, reversible Markov chain much like the results of Apers and Sarlette
in discrete-time simon2018quantum . The challenge is that most of the
underlying techniques that enable this, such as the recently developed
techniques in the context of quantum simulation low2016hamiltonian ;
chakraborty2018power ; gilyen2018quantum , are absent in continuous-time.
However, the fact that the Hamiltonian defining our continuous-time quantum
walk can be efficiently simulated using query access to the unitary defining
the discrete-time quantum walk of Ref. krovi2016quantum  might offer useful
insights towards designing such algorithms.

Our algorithm can also be used to prepare stationary states of slowly-evolving
Markov chains, i.e. given a sequence of Markov chains
$\\{P_{1},\cdots,P_{n}\\}$, such that there is a significant overlap between
the stationary distributions of any two consecutive Markov chains, meaning
$|\langle\pi_{j+1}|\pi_{j}\rangle|$ is large aharonov2003adiabatic ;
wocjan2008speedup ; orsucci2018optimal . Given that one can prepare
$\lvert\pi_{1}\rangle$ efficiently, the task is to prepare
$\lvert\pi_{n}\rangle$. Such situations arise in a host of approximation
algorithms for counting as has been pointed out in Ref. aharonov2003adiabatic
. Our algorithm will provide a quadratic speedup over that of Ref.
aharonov2003adiabatic  as given any $P_{j}$, the spectral gap of the
Hamiltonian defined in Sec. IV is amplified quadratically over the
corresponding discriminant matrix, which acts as the Hamiltonian for the
approach in aharonov2003adiabatic .

For the problem of time-averaged mixing, it would be interesting to explore
the possibility of obtaining better bounds on the quantum mixing time for any
ergodic, reversible Markov chain. Furthermore, this notion of quantum mixing
is closely related to the problem of equilibration of isolated quantum
systems, a widely studied problem in quantum statistical mechanics
wilming2018equilibration . As a result, our results can help obtain better
upper bounds for the equilibration times of isolated quantum systems defined
by random Hamiltonians.

###### Acknowledgements.

S.C. acknowledges funding from F.R.S.-FNRS. S.C. and J.R. are supported by the
Belgian Fonds de la Recherche Scientifique - FNRS under grants no F.4515.16
(QUICTIME) and R.50.05.18.F (QuantAlgo). K.L. has been partially supported by
NSF postdoctoral fellowship DMS-1702533.

## References

  * [1] Kurt Binder, Dieter Heermann, Lyle Roelofs, A John Mallinckrodt, and Susan McKay.  Monte carlo simulation in statistical physics.  Computers in Physics, 7(2):156–157, 1993.
  * [2] Emile HL Aarts and Peter JM Van Laarhoven.  Statistical cooling: A general approach to combinatorial optimization problems.  Philips J. Res., 40(4):193–226, 1985.
  * [3] Lawrence Page, Sergey Brin, Rajeev Motwani, and Terry Winograd.  The pagerank citation ranking: Bringing order to the web.  Technical report, Stanford InfoLab, 1999.
  * [4] Walter R Gilks, Sylvia Richardson, and David Spiegelhalter.  Markov chain Monte Carlo in practice.  Chapman and Hall/CRC, 1995.
  * [5] Throughout the article, we use $\mathaccentV{tilde}07E{\mathcal{O}}(.)$ and $\mathaccentV{tilde}07E{\Omega}(.)$ to suppress polylogarithmic factors, i.e. $\mathaccentV{tilde}07E{\mathcal{O}}(f(n))=\mathcal{O}(f(n)\mathrm{polylog}(f(n)))$ and $\mathaccentV{tilde}07E{\Omega}(f(n))=\Omega(f(n)\mathrm{polylog}(f(n)))$.
  * [6] Shantanav Chakraborty, Kyle Luh, and Jérémie Roland.  How fast do quantum walks mix?  Physical Review Letters, 124(5):050501, 2020.
  * [7] Dorit Aharonov and Amnon Ta-Shma.  Adiabatic quantum state generation and statistical zero knowledge.  In Proceedings of the thirty-fifth annual ACM symposium on Theory of computing, pages 20–29. ACM, 2003.
  * [8] Alfred V Aho and John E Hopcroft.  The design and analysis of computer algorithms.  Pearson Education India, 1974.
  * [9] Peter C Richter.  Quantum speedup of classical mixing processes.  Physical Review A, 76(4):042306, 2007.
  * [10] Pawel Wocjan and Anura Abeyesinghe.  Speedup via quantum sampling.  Physical Review A, 78(4):042336, 2008.
  * [11] Vedran Dunjko and Hans J Briegel.  Quantum mixing of markov chains for special distributions.  New Journal of Physics, 17(7):073004, 2015.
  * [12] Giuseppe Davide Paparo, Vedran Dunjko, Adi Makmal, Miguel Angel Martin-Delgado, and Hans J Briegel.  Quantum speedup for active learning agents.  Physical Review X, 4(3):031002, 2014.
  * [13] Davide Orsucci, Hans J Briegel, and Vedran Dunjko.  Faster quantum mixing for slowly evolving sequences of markov chains.  Quantum, 2:105, 2018.
  * [14] Vedran Dunjko and Hans J Briegel.  Machine learning & artificial intelligence in the quantum domain: a review of recent progress.  Reports on Progress in Physics, 81(7):074001, 2018.
  * [15] Peter C Richter.  Almost uniform sampling via quantum walks.  New Journal of Physics, 9(3):72, 2007.
  * [16] Mario Szegedy.  Quantum speed-up of markov chain based algorithms.  In Proceedings. 45th Annual IEEE Symposium on Foundations of Computer Science, 2004., pages 32–41. IEEE, 2004.
  * [17] Frédéric Magniez, Ashwin Nayak, Jérémie Roland, and Miklos Santha.  Search via quantum walk.  SIAM Journal on Computing, 40(1):142–164, 2011.
  * [18] Hari Krovi, Frédéric Magniez, Maris Ozols, and Jérémie Roland.  Quantum walks can find a marked element on any graph.  Algorithmica, 74(2):851–907, 2016.
  * [19] Andris Ambainis, András Gilyén, Stacey Jeffery, and Martins Kokainis.  Quadratic speedup for finding marked vertices by quantum walks.  arXiv preprint arXiv:1903.07493, 2019.
  * [20] 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(1969):339–354, 1998.
  * [21] Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  Contemporary Mathematics, 305:53–74, 2002.
  * [22] Simon Apers and Alain Sarlette.  Quantum fast-forwarding markov chains.  arXiv preprint arXiv:1804.02321, 2018.
  * [23] Shantanav Chakraborty, Leonardo Novo, and Jérémie Roland.  Finding a marked node on any graph by continuous time quantum walk.  arXiv preprint arXiv:1807.05957, 2018.
  * [24] John Von Neumann.  Mathematical foundations of quantum mechanics.  Princeton university press, 1955.
  * [25] Andrew M Childs, Enrico Deotto, Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Andrew J Landahl.  Quantum search by measurement.  Physical Review A, 66(3):032314, 2002.
  * [26] RD Somma and G Ortiz.  Quantum approach to classical thermodynamics and optimization.  In Quantum Quenching, Annealing and Computation, pages 1–20. Springer, 2010.
  * [27] Dorit Aharonov, Andris Ambainis, Julia Kempe, and Umesh Vazirani.  Quantum walks on graphs.  In Proceedings of the thirty-third annual ACM symposium on Theory of computing, pages 50–59. ACM, 2001.
  * [28] Amir Ahmadi, Ryan Belk, Christino Tamon, and Carolyn Wendler.  On mixing in continuous-time quantum walks on some circulant graphs.  Quantum Information & Computation, 3(6):611–618, 2003.
  * [29] Viv Kendon and Ben Tregenna.  Decoherence can be useful in quantum walks.  Physical Review A, 67(4):042315, 2003.
  * [30] Leonid Fedichkin, Dmitry Solenov, and Christino Tamon.  Mixing and decoherence in continuous-time quantum walks on cycles.  arXiv preprint quant-ph/0509163, 2005.
  * [31] FL Marquezino, R Portugal, G Abal, and R Donangelo.  Mixing times in quantum walks on the hypercube.  Physical Review A, 77(4):042312, 2008.
  * [32] Franklin L Marquezino, Renato Portugal, and Gonzalo Abal.  Mixing times in quantum walks on two-dimensional grids.  Physical Review A, 82(4):042341, 2010.
  * [33] Maria Kieferova and Daniel Nagaj.  Quantum walks on necklaces and mixing.  International Journal of Quantum Information, 10(02):1250025, 2012\.
  * [34] James R Norris.  Markov chains.  Number 2. Cambridge university press, 1998.
  * [35] Note that throughout this paper, we shall be dealing with values of $s\in[0,1)$ and so the properties of $D(P(1))$ are not relevant here.
  * [36] David J Aldous.  Some inequalities for reversible markov chains.  Journal of the London Mathematical Society, 2(3):564–576, 1982\.
  * [37] Hari Krovi, Maris Ozols, and Jérémie Roland.  Adiabatic condition and the quantum hitting time of markov chains.  Physical Review A, 82(2):022333, 2010.
  * [38] Andrew M Childs and Jeffrey Goldstone.  Spatial search by quantum walk.  Physical Review A, 70(2):022314, 2004.
  * [39] Jonatan Janmark, David A Meyer, and Thomas G Wong.  Global symmetry is unnecessary for fast quantum search.  Physical Review Letters, 112(21):210502, 2014.
  * [40] Iain Foulger, Sven Gnutzmann, and Gregor Tanner.  Quantum search on graphene lattices.  Physical review letters, 112(7):070504, 2014.
  * [41] Andrew M. Childs and Yimin Ge.  Spatial search by continuous-time quantum walks on crystal lattices.  Phys. Rev. A, 89:052337, May 2014.
  * [42] David A Meyer and Thomas G Wong.  Connectivity is a poor indicator of fast quantum search.  Physical review letters, 114(11):110503, 2015.
  * [43] Leonardo Novo, Shantanav Chakraborty, Masoud Mohseni, Hartmut Neven, and Yasser Omar.  Systematic dimensionality reduction for quantum walks: Optimal spatial search and transport on non-regular graphs.  Scientific reports, 5:13304, 2015.
  * [44] Pascal Philipp, Luís Tarrataca, and Stefan Boettcher.  Continuous-time quantum search on balanced trees.  Physical Review A, 93(3):032305, 2016.
  * [45] Thomas G Wong.  Quantum walk search on johnson graphs.  Journal of Physics A: Mathematical and Theoretical, 49(19):195303, 2016.
  * [46] Shantanav Chakraborty, Leonardo Novo, Andris Ambainis, and Yasser Omar.  Spatial search by quantum walk is optimal for almost all graphs.  Physical review letters, 116(10):100501, 2016.
  * [47] Shantanav Chakraborty, Leonardo Novo, Serena Di Giorgio, and Yasser Omar.  Optimal quantum spatial search on random temporal networks.  Physical review letters, 119(22):220503, 2017.
  * [48] Shantanav Chakraborty, Leonardo Novo, and Jérémie Roland.  On the optimality of spatial search by continuous-time quantum walk.  arXiv preprint arXiv:2004.12686, 2020.
  * [49] Theodore J Yoder, Guang Hao Low, and Isaac L Chuang.  Fixed-point quantum search with an optimal number of queries.  Physical review letters, 113(21):210501, 2014.
  * [50] Frank Harary and Edgar M Palmer.  Graphical enumeration.  Elsevier, 2014.
  * [51] P. Erdős and A. Rényi.  On random graphs. I.  Publ. Math. Debrecen, 6:290–297, 1959.
  * [52] P. Erdős and A. Rényi.  On the evolution of random graphs.  Publications of the Mathematical Institute of the Hungarian Academy of Sciences, 5:17–61, 1960.
  * [53] Béla Bollobás.  Random graphs.  Springer, 1998.
  * [54] Zoltán Füredi and János Komlós.  The eigenvalues of random symmetric matrices.  Combinatorica, 1(3):233–241, 1981.
  * [55] László Erdös, Antti Knowles, Horng-Tzer Yau, Jun Yin, et al.  Spectral statistics of erdős–rényi graphs i: local semicircle law.  The Annals of Probability, 41(3B):2279–2375, 2013.
  * [56] Van H. Vu.  Spectral norm of random matrices.  Combinatorica, 27(6):721–736, 2007.
  * [57] László Erdős, Horng-Tzer Yau, and Jun Yin.  Rigidity of eigenvalues of generalized wigner matrices.  Advances in Mathematics, 229(3):1435–1515, 2012.
  * [58] Terence Tao and Van Vu.  Random matrices have simple spectrum.  arXiv preprint arXiv:1412.1438, 2014.
  * [59] László Babai, D Yu Grigoryev, and David M Mount.  Isomorphism of graphs with bounded eigenvalue multiplicity.  In Proceedings of the fourteenth annual ACM symposium on Theory of computing, pages 310–324. ACM, 1982.
  * [60] Kyle Luh and Van Vu.  Sparse random matrices have simple spectrum.  arXiv preprint arXiv:1802.03662, 2018.
  * [61] Hoi Nguyen, Terence Tao, and Van Vu.  Random matrices: tail bounds for gaps between eigenvalues.  Probability Theory and Related Fields, pages 1–40, 2015.
  * [62] Patrick Lopatto and Kyle Luh.  Tail bounds for gaps between eigenvalues of sparse random matrices.  arXiv preprint arXiv:1901.05948, 2019.
  * [63] Yael Dekel, James R Lee, and Nathan Linial.  Eigenvectors of random graphs: Nodal domains.  Random Structures & Algorithms, 39(1):39–58, 2011.
  * [64] Yukun He, Antti Knowles, and Matteo Marcozzi.  Local law and complete eigenvector delocalization for supercritical erdös-rényi graphs.  arXiv preprint arXiv:1808.09437, 2018.
  * [65] Jiaoyang Huang, Benjamin Landon, and Horng-Tzer Yau.  Transition from tracy-widom to gaussian fluctuations of extremal eigenvalues of sparse erd$\backslash$h $\\{$o$\\}$ sr$\backslash$’enyi graphs.  arXiv preprint arXiv:1712.03936, 2017.
  * [66] Guang Hao Low and Isaac L Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, 2019.
  * [67] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The power of block-encoded matrix powers: improved regression techniques via faster hamiltonian simulation.  Proceedings of the 46th International Colloquium on Automata, Languages and Programming (ICALP), pages 33:1–33:14, 2019.
  * [68] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.
  * [69] Henrik Wilming, Thiago R de Oliveira, Anthony J Short, and Jens Eisert.  Equilibration times in closed quantum many-body systems.  In Thermodynamics in the Quantum Regime, pages 435–455. Springer, 2018.