Analog quantum algorithms for the mixing of Markov chains

TL;DR

This paper introduces an analog quantum algorithm for Markov chain mixing that improves sampling efficiency, and explores quantum mixing times, extending previous results with new theoretical insights and numerical validation.

Contribution

It presents a novel analog quantum algorithm for Markov chain mixing and extends understanding of quantum mixing times using advanced random matrix theory tools.

Findings

Quantum algorithm scales with the sum of square roots of classical mixing and hitting times.

Provided an intuitive understanding of eigenvalue statistics in random matrices.

Numerical simulations support the analytical results and extend to general Markov chains.

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\"os-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.

Equations

$$\pi =\left({\pi }_1{\pi }_2\cdots {\pi }_n\right),$$

$$p_{x,y}^{\ast }=p_{y,x}\frac{{\pi }_y}{{\pi }_x}.$$

$${\pi }_x{p}_{xy}={\pi }_y{p}_{yx},\forall (x,y)\in X.$$

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

$$P=\left[\begin{array}{cc}{P}_{UU}& P_{UM}\\ P_{MU}& P_{MM}\end{array}\right],P^{\prime }=\left[\begin{array}{cc}{P}_{UU}& P_{UM}\\ 0& I\end{array}\right],$$

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

$$P=\left[\begin{array}{cc}{P}_{UU}& P_{UM}\\ (1-s)P_{MU}& (1-s)P_{MM}+sI\end{array}\right].$$

$$\pi =({\pi }_U{\pi }_M),$$

$$\pi (s)=\frac{1}{1-s(1-p_M)}\left((1-s){\pi }_U{\pi }_M\right).$$

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

$$D_{xy}(P(s))$$

$$=p_{xy}(s)\sqrt{\frac{{\pi }_x(s)}{{\pi }_y(s)}}.$$

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

$$D(P(s))=\sum _{i=1}^n{\lambda }_i(s)|v_i(s)⟩⟨v_i(s)|,$$

$$|v_n(s)⟩=\sqrt{\pi (s{)}^T},$$

$$D(P(s))\sqrt{\pi (s{)}^T}$$

$$=\sqrt{\pi (s{)}^T}.$$

$$|v_n(s)⟩=\sqrt{\frac{(1-s)(1-p_M)}{1-s(1-p_M)}}|U⟩+\sqrt{\frac{p_M}{1-s(1-p_M)}}|M⟩,$$

$$|U⟩$$

$$|M⟩$$

$$HT(P,M)=\sum _{j=1}^{n-m}\frac{|⟨v_j^{\prime }|U⟩{|}^2}{1-{\lambda }_j^{\prime }},$$

$$|U⟩=\frac{1}{\sqrt{1-p_M}}\sum _{x\ne M}\sqrt{{\pi }_x}|x⟩,$$

$$HT(s)=\sum _{j=1}^{n-1}\frac{|⟨v_j(s)|U⟩{|}^2}{1-{\lambda }_j(s)}.$$

$$HT(s)\le \frac{1}{\Delta (s)}\sum _{j=1}^{n-1}|⟨v_j(s)|U⟩{|}^2.$$

$$H{T}^{+}(P,M)=\underset{s\to 1}{\lim }HT(s),$$

$$HT(s)=\frac{p_M^2}{{\left(1-s(1-p_M)\right)}^2}H{T}^{+}(P,M).$$

$$H{T}^{+}(P,M)\le \frac{1}{\Delta (s)}.\frac{{\left(1-s(1-p_M)\right)}^2}{p_M^2}\sum _{j=1}^{n-1}|⟨v_j(s)|U⟩{|}^2.$$

$$\frac{1}{2}\Vert \mu P^t-\pi {\Vert }_1\le ϵ,$$

$$T_{\mathrm{mix}}\le \frac{1}{\Delta }\log \left(\frac{1}{ϵ{\pi }_{\min }}\right),$$

$$\tilde{H}=H+\frac{{\hat{p}}^2}{2m}+g\hat{O}\otimes \hat{p},$$

$$\tilde{H}=H\otimes \hat{p}.$$

$$e^{-i{x}_0\hat{p}}\psi (x)$$

Full text

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 [1], combinatorial optimization [2] to network science [3] and form the basis of Markov chain Monte Carlo-based methods [4]. 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 arise 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. [6], where we prove an upper bound for the quantum mixing time for almost all graphs.

Aharonov and Ta-Shma [7] 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 [8] (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 [9, 10, 11], some of which have even been instrumental in obtaining speedups for quantum machine learning [12, 13, 14]. Richter [15] 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 [16].

The key idea that encompasses all existing algorithms for $QSSamp$ is to make use of the so called quantum spatial search algorithm [17]. 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) [18, 19]. 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 [20] and quantum amplitude amplification [21] 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 [22]. 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 [23]. 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 [24] 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 [25]. 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. [23] which makes use of quantum phase randomization [26]. 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. [23], 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. [27] 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 [28, 29, 30, 31, 32, 33]. 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 [6]. This implies that the fraction of graphs of $n$ nodes for which our upper bound holds, goes to 1 as $n$ goes to infinity or equivalently, if a graph is picked uniformly at random from the set of all graphs, our result provides an upper bound on the quantum mixing time almost surely, i.e. with probability $1-o(1)$. Throughout the article, we shall use the phrase almost all graphs to signify precisely this.

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 [6]. 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$ [34]. 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

[TABLE]

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 [math] 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}$ denote 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

[TABLE]

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

[TABLE]

This can also be rewritten as

[TABLE]

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$ to $X$ 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

**Discriminant matrix: ** We denote by

[TABLE]

the symmetric matrix whose $(x,y)^{\mathrm{th}}$ entry is $D_{xy}(P(s))=\sqrt{p_{xy}(s)p_{yx}(s)}$. Here $\circ$ indicates the Hadamard product.  

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

[TABLE]

This leads us to the following fact:

Fact 1

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

[TABLE]

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

From Fact 1, it follows 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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. [18, 23], 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

[TABLE]

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

[TABLE]

For the spatial search algorithm, we shall find that the quantity 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

[TABLE]

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)$ [18]. They showed that

[TABLE]

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

[TABLE]

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

[TABLE]

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 [36]:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math]. Then,

[TABLE]

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

[TABLE]

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 [math]-eigenstate of $H$, i.e. $\lvert v_{n}\rangle$, in a purely analog fashion. To that end, we formally state via Lemma 4 and Corollary 5.

Lemma 4

Let $H$ be a Hamiltonian with eigenvalues $\lambda_{n}=0<\Delta=\left|\lambda_{n-1}\right|\leq\cdots\left|\lambda_{1}\right|\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

[TABLE]

where

[TABLE]

for some $\tau>0$. Furthermore let

[TABLE]

Then starting from the state $\lvert\psi_{0}\rangle\lvert x=0\rangle$ and evolving for a time $\tau$ according to the Hamiltonian $\widetilde{H}=H\otimes\widehat{p}$, results in a state

[TABLE]

where $|\gamma_{k}|\leq\pi/(\left|\lambda_{k}\right|\tau),~{}|\Gamma_{k}|=\sqrt{1-|\gamma_{k}|^{2}}$ and $\langle\Gamma_{k}|0\rangle=0$ for $1\leq k\leq n-1$.

In particular, for

[TABLE]

we have $|\gamma_{k}|\leq 1/2$ and $|\Gamma_{k}|\geq\sqrt{3}/{2}$ for $1\leq k\leq n-1$.

**Proof. **

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

[TABLE]

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,

[TABLE]

Therefore,

[TABLE]

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

[TABLE]

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 eigenvalue corresponding 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 [math]-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.

[TABLE]

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

[TABLE]

Since $|1-e^{-iz}|\leq 2$ for any $z$ and $|1-e^{-iz}|\geq 2|z|/\pi$ for any $z\in[-\pi,\pi]$, we can bound this amplitude as

[TABLE]

where we have used the fact that $l\geq\log_{2}(\tau/\pi)$ so $\left|\lambda_{k}\right|\tau/2^{l}\leq\pi$.

Thus, after the time evolution for a time $\tau$, the state of the system and the pointer is given by

[TABLE]

where $|\gamma_{k}|\leq\pi/(\left|\lambda_{k}\right|\tau),~{}|\Gamma_{k}|=\sqrt{1-|\gamma_{k}|^{2}}$ and $\langle\Gamma_{k}|0\rangle=0$ for $1\leq k\leq n-1$.

Since $\left|\lambda_{k}\right|\geq\Delta$, for $\tau=2\pi/\Delta$ we have $|\gamma_{k}|\leq 1/2$ and $|\Gamma_{k}|\geq\sqrt{3}/{2}$ for $1\leq k\leq n-1$. $\Box$

 

We shall use this lemma to derive the following corollary.

Corollary 5

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

[TABLE]

qubits initialized in the state $\lvert x=0\rangle^{\otimes m}$, where $l=\lceil\log_{2}(\tau/\pi)\rceil$.

Then repeating the Hamiltonian evolution of Lemma 4 with $\tau=\dfrac{2\pi}{\Delta}$ 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}}$, succeeds with probability at least $|\alpha_{n}|^{2}$ in constructing a quantum state $\lvert\phi\rangle$ such that

[TABLE]

in time

[TABLE]

**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

[TABLE]

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

[TABLE]

where $0\leq\epsilon_{k}\leq\epsilon^{\prime}$ and $\sqrt{1-\epsilon^{\prime 2}}\leq\delta_{k}\leq 1$. This takes time

[TABLE]

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

[TABLE]

where the state

[TABLE]

has norm $\delta$ with

[TABLE]

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

[TABLE]

with probability $|\alpha_{n}|^{2}(1+\delta^{2})\geq|\alpha_{n}|^{2}$, such that

[TABLE]

$\Box$

Thus Lemma 4 and corollary 5 can be used to prepare the eigenstate $\lvert v_{n}\rangle$ with success probability at least $|\alpha_{n}|^{2}$, and repeating this procedure $\Theta(1/|\alpha_{n}|^{2})$ times allows to boost this probability to $\Theta(1)$, which leads to a protocol using total time

[TABLE]

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 [26] and subsequently used in Refs. [37, 23]. 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$,

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

In Ref. [23], 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

[TABLE]

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. [37] 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$,

[TABLE]

Now, observe that

[TABLE]

This implies

[TABLE]

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

 

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

[TABLE]

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

[TABLE]

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 [math] 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)$.

The following connection, first observed in [37], between the eigenstates of $H(s)$ in $\mathcal{B}_{k}(s)\cup\mathcal{B}_{n}(s)$ and the interpolated hitting time $HT(s)$ introduced in Sec. II.2 will be useful to the complexity analysis of our algorithms

Proposition 6

Let $h(s)$ be defined as

[TABLE]

Then, we have

[TABLE]

**Proof. ** From Eqs. (59-60), we have

[TABLE]

so that

[TABLE]

The bounds on $h(s)$ then follow from the definition of $HT(s)$ and the fact that

[TABLE]

$\Box$

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 [38] introduced the first continuous-time quantum walk-based algorithm to tackle the spatial search problem for simple, unweighted graphs. They showed that this algorithm, 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 [39, 40, 41, 42, 43, 44, 45, 46, 47]. Recently in Ref. [48], 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. [23], 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$ [18, 19]. Given any $P$, their algorithm made use of the framework of interpolating quantum walks

[TABLE]

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 [math]-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. [23] in that both make use of the Somma-Ortiz Hamiltonian $H(s)$ defined in Eq. (50).

However, motivated by the problem of quantum state generation using von Neumann measurements, Algorithm 2 prepares the [math]-eigenstate of $H(s)$ by coupling this Hamiltonian to a free-particle in one dimension. Unlike the algorithm of Ref. [23], we evolve the Hamiltonian for a fixed time $\tau=\frac{\pi}{\delta}\sqrt{2HT(s^{*})}=\Theta(\sqrt{HT^{+}(P,M)})$, where

[TABLE]

Furthermore, the [math]-eigenstate $\lvert v_{n}(s^{*})\rangle$ can be written as

[TABLE]

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$. The intuition behind the algorithm is as follows. In step 1, it either detects a marked node, in which case it succeeds with probability 1, or it prepares state $\lvert U,0\rangle$. Since this state has overlap $1/\sqrt{2}$ with $\lvert v_{n}(s^{*}),0\rangle$, step 2 allows to prepare a state such that, if we were to project the third register on $\lvert 0\rangle$, we would obtain with probability close to $1/2$ a state close to $\lvert v_{n}(s^{*}),0\rangle$. This state therefore has overlap close to $1/\sqrt{2}$ over $\lvert M,0\rangle$, so that step 3 will yield a marked node with probability close to $1/2$.

We now formally state Algorithm 2 and prove its correctness in Lemma 7.

Lemma 7

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

[TABLE]

**Proof. ** We shall make use of Lemma 4. Observe that for $s=s^{*}=1-p_{M}/(1-p_{M})$, the [math]-eigenstate of $H(s^{*})$ is simply $\lvert v_{n}(s^{*}),0\rangle$ with

[TABLE]

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

[TABLE]

so at the end of step 1, we project on $\lvert M\rangle$ with probability $p_{M}$ (in which case we will immediately obtain a marked node in step 3), or we project on $\lvert U\rangle$.

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

[TABLE]

where $\alpha_{n}=\langle v_{n}(s^{*})|U\rangle=\frac{1}{\sqrt{2}}$ and $\alpha_{k}^{\sigma}=\langle\Psi_{k}^{\sigma}(s^{*})|U,0\rangle$, prepares a state

[TABLE]

where $|\gamma_{k}^{\sigma}|\leq\pi/(\left|E_{k}^{\sigma}(s^{*})\right|\tau)$ and $\langle\Gamma_{k}^{\sigma}|0\rangle=0$ for $1\leq k\leq n-1$.

Consider the measurement operator

[TABLE]

where $\Pi_{X}$ is a projection on the states of the Markov chain. Since $\langle\Gamma_{k}^{\pm}|0\rangle=0$, the probability to measure [math] in the third register is

[TABLE]

From $|\alpha_{n}|^{2}=\frac{1}{2}$ and $|\gamma_{k}^{\sigma}|\leq\pi/(\left|E_{k}^{\sigma}(s^{*})\right|\tau)$, we can bound $p_{0}$ as

[TABLE]

where we have used Prop. 6 and the fact that $\tau=\frac{\pi}{\delta}\sqrt{2HT(s^{*})}$.

In summary, we prepare with probability at least $1/2$ the state

[TABLE]

where, from the calculation above, the second term has norm at most $\delta/\sqrt{2}$. This state is therefore close to $\lvert v_{n}(s^{*}),0\rangle$, more precisely, we have

[TABLE]

Let us finally consider the measurement in the state-space of the Markov chain. If the state we had prepared was $\lvert v_{n}(s^{*}),0\rangle$ we would obtain a marked node with probability $1/2$ since $|\langle v_{n}(s^{*})|M\rangle|^{2}=1/2$. Equivalently, if we define the projector $\Pi_{M}$ to be the projector over marked nodes of the Markov chain, we have $\left\lVert\Pi_{M}\lvert v_{n}(s^{*}),0\rangle\right\rVert=\frac{1}{\sqrt{2}}$. The triangle inequality then implies

[TABLE]

This implies that measuring $\lvert\psi_{2}\rangle$ in the state-space of the Markov chain yields a marked node with probability at least $\left\lVert\Pi_{M}\lvert\psi_{2}\rangle\right\rVert^{2}\geq\frac{1}{2}-2\delta$. Since, in the previous step we had prepared this state with probability at least $1/2$, overall the algorithm outputs a marked node with probability at least $1/4-\delta$.

Finally, Eq. (20) implies that for $s=s^{*}=1-p_{M}/(1-p_{M})$, we have $HT(s^{*})=\frac{HT^{+}(P,M)}{4}$, so that the time to execute Step 2 of the algorithm is $\tau=\frac{\pi}{\delta}\sqrt{2HT(s^{*})}=\frac{\pi}{\delta}\sqrt{\frac{HT^{+}(P,M)}{2}}$, which concludes the proof. $\Box$

Algorithm 2 can be thought of as an analog version of the quantum spatial search algorithm by Krovi et al [18]. Although this algorithm requires $\Theta(\log_{2}HT^{+}(P,M))$ additional ancillary qubits unlike the algorithm of Ref. [23], 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. [23], 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. [18]. 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. [18, 19], 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

[TABLE]

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

[TABLE]

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

[TABLE]

Thus, this is simply the highest eigenstate of the discriminant matrix $D(P)$ or equivalently, the [math]-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 [49], 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

[TABLE]

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

[TABLE]

and from Eq. (13) we have that

[TABLE]

Thus, the state $\lvert j\rangle$ has a constant overlap with $\lvert v_{n}(s^{*})\rangle$. We can use this fact to establish the following connection between the average hitting time and the eigenstates of $H(s^{*})$, following Proposition 6.

Proposition 8

Let $h(s)$ be as defined in Proposition 6 and $HT(P,\{j\})$ be the average hitting time of $P$ with respect to the vertex $j$. Then, for $s^{*}=1-\pi_{j}/(1-\pi_{j})$,

[TABLE]

Furthermore,

[TABLE]

**Proof. ** For $s=s^{*}$, we have that for any $1\leq k\leq n-1$, $\langle v_{k}(s^{*})|v_{n}(s^{*})\rangle=0$ and substituting the expression of $\lvert v_{n}(s^{*})\rangle$, we obtain that for any $k$ such that $1\leq k\leq n-1$,

[TABLE]

So, we have:

[TABLE]

Following proposition 6, we have

[TABLE]

Since we have a single marked node $j$, the extended hitting time with respect to $j$ is equal to $HT(P,\{j\})$. The final expression of the proposition can be obtained by observing that from Eq. (20), we have $HT(s^{*})=HT(P,\{j\})/4$. $\Box$

Recall 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$. 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 first invoke Lemma 4 to prepare an intermediate state that has a constant overlap with $v_{n}(s^{*})$. We will prove that this can be achieved by setting $s=s^{*}=1-\pi_{j}/(1-\pi_{j})$, $\tau=\Theta(\sqrt{HT(P,\{j\})})$ and choosing the initial state to be $\lvert j,0\rangle$. For the second stage, we set $s=0$ and start from the state obtained in stage 1. By invoking Corollary 5, we prepare a state that is $\epsilon$-close to $\lvert v_{n}(0)\rangle=\lvert\pi\rangle$ with a constant probability. 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 9.

Lemma 9

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

[TABLE]

in time

[TABLE]

*with success probability at least $1/4-\delta$. *

**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

[TABLE]

Thus, the state $\lvert j,0\rangle$ can be written in the eigenbasis of $H(s^{*})$ as

[TABLE]

where $\alpha_{n}=1/\sqrt{2}$ and $\alpha^{\sigma}_{k}=\langle\Psi^{\sigma}_{k}(s^{*})|j,0\rangle$.

Now, we shall evolve the state $\lvert j,0\rangle\lvert x=0\rangle$ according to the Hamiltonian $\tilde{H}(s^{*})=H(s^{*})\otimes\hat{p}$ for a time $\tau$. From Lemma 4, we have that the resulting state

[TABLE]

where $|\gamma_{k}^{\sigma}|\leq\pi/(\left|E_{k}^{\sigma}(s^{*})\right|\tau)$ and $\langle\Gamma_{k}^{\sigma}|0\rangle=0$ for $1\leq k\leq n-1$.

The probability to measure $\lvert 0\rangle$ in the third register is

[TABLE]

From $|\alpha_{n}|^{2}=\frac{1}{2}$, $|\gamma_{k}^{\sigma}|\leq\pi/(\left|E_{k}^{\sigma}(s^{*})\right|\tau)$ and the analysis in Lemma 7, we can bound $p_{0}$ as

[TABLE]

where we have used Prop. 8 and the fact that $\tau=\frac{\pi}{\delta}\sqrt{HT(P,\{j\})/2}$. Thus, after stage one, with probability $p_{0}\geq 1/2$, we have prepared the state

[TABLE]

in the first register, where $\left\lVert\lvert\text{err}\rangle\right\rVert\leq\sqrt{2}\delta$. The running time of stage one of the algorithm is

[TABLE]

For any constant $\delta\in(0,1)$, we have $T_{1}=\Theta(\sqrt{HT(P,\{j\}})$.

We shall use this state as the initial state for stage two of our algorithm. In this stage, we set $s=0$ and prepare a state that is $\epsilon$-close to $\lvert v_{n}(0)\rangle$. Observe that the overlap $\langle v_{n}(0)|v_{n}(s^{*})\rangle=1/\sqrt{2}$, so

[TABLE]

We begin by expressing $\lvert\psi^{(1)}_{f}\rangle$ in the eigenbasis of $H$. We have

[TABLE]

where $\left|\delta_{n}\right|\geq(1-2\delta)/\sqrt{2}$ and for any $1\leq k\leq n-1$, $\delta_{k}^{\sigma}=\langle\Psi_{k}^{\sigma}(0)|\psi^{(1)}_{f}\rangle$.

Let $\Delta(0)=\sqrt{1-\lambda^{2}_{n-1}(0)}$, be the eigenvalue gap between the [math]-eigenstate of $H(0)$ and the rest of the spectrum. So, if $\Delta=1-\lambda_{n-1}(0)$ is the spectral gap of $P$, $\sqrt{\Delta}\leq\Delta(0)\leq\sqrt{2\Delta}$.

From Lemma 4, by choosing $\tau=2\pi/\Delta(0)$, we obtain a state

[TABLE]

where now, for any $k$ such that $1\leq k\leq n-1$, $|\gamma_{k}^{\sigma}|\leq 1/2$, $|\Gamma_{k}|\geq\sqrt{3}/2$ and $\langle\Gamma_{k}|0\rangle=0$.

We can now use Corollary 5, with $l=1+\lceil\log_{2}(1/\Delta(0))\rceil$, $\epsilon^{\prime}=\epsilon/4\leq|\delta_{n}|\epsilon/\sqrt{2}$ and $m=l\cdot\lceil\log(1/\epsilon^{\prime})\rceil$. By choosing $\tau=2\pi/\Delta(0)$, we prepare a quantum state $\lvert\phi_{f}\rangle$ that is $\epsilon$-close to $\lvert v_{n}(0)\rangle$ with probability at least $\left|\delta_{n}\right|^{2}\geq 1/2-2\delta$. The running time of the second stage of the algorithm is in

[TABLE]

Now putting things together, Algorithm 3 prepares a state $\lvert\phi_{f}\rangle$ which is $\epsilon$-close to $\lvert v_{n}(0)\rangle$. The first stage succeeded with probability $p_{0}\geq 1/2$ and so the overall success probability is $\geq 1/4-\delta$. The overall running time

[TABLE]

Interestingly, by Eq. (22) this also implies 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.

[TABLE]

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,

[TABLE]

$\Box$

 

 

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

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. 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: the probability that the walker is at some node $f$ after some time $t$, picked uniformly at random in the interval $[0,T]$ [27]. 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.

[TABLE]

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 [math] 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

[TABLE]

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

[TABLE]

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.

[TABLE]

as long as

[TABLE]

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

[TABLE]

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. (86), we need to evaluate the following quantity in order to upper-bound $T_{\mathrm{mix}}$:

[TABLE]

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.

[TABLE]

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 10

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

[TABLE]

**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

[TABLE]

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

[TABLE]

$\Box$

 

The upper bound on $\Sigma$ obtained in Lemma 10 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. (86) is less than one and so

[TABLE]

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

[TABLE]

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

In this section, we elaborate on the results of Ref. [6] 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}$ [50]. 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 [51, 52, 53] (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.

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$ [51, 52]. 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 [math]. For example 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 [math] with probability $1-p$. All diagonal entries of $A_{G(n,p)}$ are [math]. 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. [54] for constant $p$ and was later improved for sparse random graphs ($p=o(1)$) in Ref. [55]. In fact, as we shall show shortly it suffices to consider the matrix

[TABLE]

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$,

[TABLE]

with probability $1-o(1)$, which implies that $\left\lVert\bar{A}_{G(n,p)}\right\rVert\approx 1$ [55, 6].

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

[TABLE]

with probability $1-o(1)$ [54, 56, 6]. 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. (86) 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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$,

[TABLE]

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. [57, 55]. 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

[TABLE]

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

[TABLE]

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

[TABLE]

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)$,

[TABLE]

From Eq. (103) and Eq. (104), we obtain that

[TABLE]

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 [58] showed that $\bar{A}_{G(n,p)}$ has a simple spectrum for dense graphs, resolving a long-standing conjecture due to Babai [59]. Recently, this was resolved also for sparse graphs [60]. We state their results formally.

Fact 11

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. (87) can now be re-written as

[TABLE]

while the limiting probability distribution is

[TABLE]

Furthermore, Nguyen, Tao and Vu [61] 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 [62]. 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}$,

[TABLE]

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. (88)) for $\bar{A}_{G(n,p)}$. We prove that here.

Lemma 12

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

[TABLE]

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. (110), we obtain

[TABLE]

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. (112). By choosing

[TABLE]

we have that

[TABLE]

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,

[TABLE]

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. (108).

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

[TABLE]

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. (110). On the other hand for $r>r^{\star}(i)$, the eigenvalue rigidity criterion (See Eq. (104)) 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.

[TABLE]

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

Lemma 13** (Upper bound on $\mathbf{\Sigma_{1}}$ [6])**

[TABLE]

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. [6].  

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. [6].

Lemma 14** (Upper bound on $\mathbf{\Sigma}$ [6])**

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

[TABLE]

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. [6] for details. We first split $\Sigma$ into two different parts:

[TABLE]

For the first sum in the right hand side of Eq. (118), 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. (103).

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.

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. [63], 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. [55] 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. [64]. This implies that as long as $p\geq\omega(\log(n)/n)$, for all $j\in\{1,\cdots,n\}$

[TABLE]

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

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

with probability $1-o(1)$.  

 

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

[TABLE]

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

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

[TABLE]

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 13 along with the observation that

[TABLE]

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

[TABLE]

In fact, the breakdown of rigidity estimates in [55] 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 [65].

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. (114)) for which the bounds in [62] 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. (127).

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 $P$. Any symmetric matrix that captures the local connectivity of $P$ can be a used as a Hamiltonian for performing a CTQW on $P$. As $P$ need not be symmetric in general, one cannot consider a quantum walk on $P$ directly. 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$.

Here, we shall explore whether any generic speedup is obtained for the $QLSamp$ problem. Note that the time evolution of some initial state $\lvert\psi(0),0\rangle$, under the action of $H$ is given by

[TABLE]

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

[TABLE]

This implies,

[TABLE]

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

[TABLE]

where recall from Eq. (59) 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

[TABLE]

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 $P$, $\Delta_{\min}$. To that end we have the following lemma:

Lemma 15

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

[TABLE]

**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

[TABLE]

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

[TABLE]

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

[TABLE]

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$

 

So from Lemma 15, we have that for any ergodic, reversible Markov chain $P$

[TABLE]

Let us now consider that $P$ is a symmetric, i.e. $P=P^{T}$. Then the underlying quantum walk can also be performed on $P$ itself. Assuming that the eigenvalues of $P$ are ordered, for a continuous-time quantum walk on $H$, from Eq. (132) and Lemma 15, we observe that the upper bound for the quantum mixing time may be faster or slower than a quantum walk performed on $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}$, $\tilde{\Delta}_{\min}\approx\Delta_{\min}/\sqrt{\Delta}$ and hence the upper bound on the quantum mixing time is in $\widetilde{\mathcal{O}}\left(n\sqrt{\Delta}/\Delta_{\min}\right)$ which is faster than the bound in Eq. (94). On the other hand, if $\Delta_{\min}$ is in the vicinity of $\lambda_{2}$, the upper bound on the quantum mixing time is given by Eq. (139).

For generic ergodic, reversible Markov chains however, this comparison is inapplicable as $P$ may not be symmetric and cannot be used as a CTQW Hamiltonian.

This is in contrast to the $QSSamp$ problem, where using $H$ offers a generic quadratic speedup over using $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. [6] 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 [38, 48] 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 [22]. The challenge is that most of the underlying techniques that enable this, such as the recently developed techniques in the context of quantum simulation [66, 67, 68], 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. [18] 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 [7, 10, 13]. 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. [7]. Our algorithm will provide a quadratic speedup over that of Ref. [7] 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 [7].

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 [69]. 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.