Mixing Time of Quantum Gibbs Sampling for Random Sparse Hamiltonians

The recently proposed CKG quantum MCMC algorithm addresses the problem of finding thermal states by imitating thermodynamic processes [11, 10]. In this process, a system of particles evolves in contact with a thermal bath at some fixed temperature ${\beta }^{-1}$. Due to interactions with the bath, the system is described by a probabilistic mixture of quantum states $𝝆$. This state evolves in time, by approximation, with Markovian dissipative dynamics, $\frac{\mathrm{d}𝝆}{\mathrm{d}t}={ℒ}_{\beta }[𝝆]$, given in terms of an operator ${ℒ}_{\beta }$ known as the Lindbladian. This operator involves a coherent term $𝑩$ that describes the interaction among the particles in the system. There is also a term in ${ℒ}_{\beta }$ that is specified by a series of Lindblad operators ${𝑳}^j$ that drive the dissipative transitions. The dynamics of the coherent term $𝑩$ are reversible, while the dissipative transitions drive all states toward some “stationary state”. These transitions can be understood as perturbations from the bath, and as all states converge to the Gibbs state, the information of the system is leaking via these perturbations to the bath. The expression, in terms of $𝑩$ and ${𝑳}^j$, is:

The summands ${𝑳}^j(\cdot ){𝑳}^{j†}$ are termed the transition part of the Lindbladian, and $-\frac{1}{2}\{{𝑳}^{j†}{𝑳}^j,\cdot \}$ are the decay part of the Lindbladian. The choice of the Lindbladian operator ${ℒ}_{\beta }$ can vary depending on the precise nature of interactions between the system and the bath. However, to prepare the Gibbs (thermal) state at temperature ${\beta }^{-1}$, the Lindbladian should satisfy

and moreover ${𝝆}_{\beta }$ should be the unique stationary state of the Lindbladian. The long-term evolution of the system under this Lindbladian, as a result, would converge to the Gibbs state of the Hamiltonian $𝑯$ at temperature $1/\beta$.

To ensure that the Lindbladian ${ℒ}_{\beta }$ converges to a state ${𝝆}_{\beta }$, [11] designs a Lindbladian that satisfies Kubo-Martin-Schwinger (KMS) detailed balance with respect to ${𝝆}_{\beta }$. KMS detailed balance is one of several ways of quantizing the notion of classical detailed balance for Markov chains. KMS detailed balance of ${ℒ}_{\beta }$ is self-adjointness with respect to the inner product

In particular, it is equivalent to the relation that

where ${ℒ}_{\beta }^{†}$ is the adjoint Lindbladian with respect to the Hilbert-Schmidt inner product $⟨{\sigma }_1,{\sigma }_2⟩=\mathrm{Tr}({\sigma }_1^{†}{\sigma }_2)$. The adjoint operator ${ℒ}_{\beta }^{†}$, in the Heisenberg picture, describes the dynamics of observables under evolution by ${ℒ}_{\beta }$. The Lindbladian evolution is described by some quantum channel and therefore the observable $I$ must always be fixed by $\exp ({ℒ}_{\beta }^{†})$. This implies that ${ℒ}_{\beta }^{†}[I]=0$. The detailed balance formula thereby implies that ${ℒ}_{\beta }[{𝝆}_{\beta }]=0$, as desired. Note that KMS detailed balance can be dually described as the self-adjointness of ${ℒ}_{\beta }^{†}$ with respect to the inner product ${⟨{\sigma }_1,{\sigma }_2⟩}_{{\rho }_{\beta }}=\mathrm{Tr}({\sigma }_1^{†}{𝝆}_{\beta }^{1/2}{\sigma }_2{𝝆}_{\beta }^{1/2})$.

The quantum Gibbs sampler in [11] constructs a Lindbladian that satisfies two properties:

1. 1.

   ${ℒ}_{\beta }$ satisfies detailed balance with respect to ${𝝆}_{\beta }$, and therefore ${ℒ}_{\beta }[{𝝆}_{\beta }]=0$.

2. 2.

   The dynamics of ${ℒ}_{\beta }$ can be efficiently implemented.

Their Lindbladian, which we term the CKG Lindbladian, can be simulated on a quantum computer with a cost per unit time $t=1$ roughly equal to that of simulating the Hamiltonian dynamics of $𝑯$. The CKG Lindbladian is closely related to the Davies generator, which is a physically motivated Lindbladian that satisfies detailed balance, but that is not efficiently implementable in general. A full description of both Lindbladians are given in the Appendix.

The Gibbs sampling algorithm evolves an initial state ${𝝆}_0$ according to the efficiently implemented Lindbladian ${ℒ}_{\beta }$, and produces the state ${𝝆}_t=e^{{ℒ}_{\beta }t}[{𝝆}_0]$ after time period $t$. The mixing time is roughly the time that it takes for the state ${𝝆}_t$ to approach the Gibbs state ${𝝆}_{\beta }$. That is, ${\mathrm{e}}^{{ℒ}_{\beta }{t}_{\text{mix}}}[{𝝆}_0]\approx {𝝆}_{\beta }.$ The efficiency of the algorithm therefore scales linearly with the unit time simulation cost and the mixing time. The algorithm has several parameters in the Lindbladian’s construction. In addition to the inverse temperature $\beta$, the algorithm specifies an energy resolution ${\sigma }_E$. A salient feature of [11]’s construction is that it can achieve detailed balance even though the algorithm only probes the energies of the Hamiltonian $𝑯$ with approximate precision. ${\sigma }_E$ quantifies this level of precision. The cost of the Lindbladian simulation depends linearly on ${\sigma }_E^{-1}$, but increasing the precision may also improve the mixing time. Taking ${\sigma }_E\to 0$ for absolute precision recovers the Davies generator – when distinguishing the energies of the system exactly is infeasible, this Lindbladian cannot be simulated efficiently.

A set of jump operators ${𝑨}^a$ must also be specified for the Lindbladian. These operators are decomposed by frequency and reassembled in a particular way to construct the Lindblad operators that help ${ℒ}_{\beta }$ satisfy detailed balance. They must appear in adjoint pairs: i.e., if $𝑨\in \{{𝑨}^a\}$, then ${𝑨}^{†}\in \{{𝑨}^a\}$. The cost of simulation scales with the cost of implementing the oracle $|a⟩\to |a⟩\otimes {𝑨}^a$. In particular, the jump operators must be normalized when implemented for the algorithm, satisfying ${\sum }_a\Vert {𝑨}^{a†}{𝑨}^a\Vert \le 1$. CKG Lindbladians are linear in their jump operators – if ${ℒ}_1$ has one jump operator ${𝑨}^1$ and ${ℒ}_2$ has one jump operator ${𝑨}^2$, then a Lindbladian $ℒ$ with jump operators ${𝑨}^1$ and ${𝑨}^2$ satisfies $ℒ={ℒ}_1+{ℒ}_2$. If ${ℒ}_{\beta }$ was constructed from jumps ${𝑨}^a$, then jump operators $\sqrt{s}{𝑨}^a$ produce the Lindbladian $s{ℒ}_{\beta }$, scaling the mixing time by $s$. So we may therefore assume that ${\sum }_a\Vert {𝑨}^{a†}{𝑨}^a\Vert =1$ exactly, since renormalizing can only improve the spectral gap. In its normalized form, the set of jump operators can be understood as a jumping distribution over ${𝑨}^a$ which we will notate $a\sim 𝒜$, where each is sampled with probability $\Vert {𝑨}^{a†}{𝑨}^a\Vert$.

In the description of the Davies generator for a given system $𝑯$, there is a coherent term and there are jump operators ${𝑨}^a$. The ${𝑨}^a$ terms must appear in adjoint pairs in the Davies generator. The dissipative part of the Lindbladian is expressed as follows:

where ${𝑨}_{\omega }^a$ is the Operator Fourier Transform (OFT) of jump operator ${𝑨}^a$:

The Davies’ generator chooses Lindblad operators $\sqrt{\gamma (\omega )}{𝑨}_{\omega }^a$, each of which selects the energy transitions, or Bohr frequencies, in ${𝑨}^a$ that are precisely $\omega$. Because it requires certainty in energy, by Heisenberg’s uncertainty principle of energy and time, in the general case simulating the evolution of the Davies generator efficiently is infeasible. In the above, $\gamma$ is some function satisfying $\gamma (\omega )=\gamma (\omega )=e^{-\beta \omega }\gamma (-\omega )$. The Lindblad operators are scaled by $\gamma (\omega )$ precisely to satisfy KMS detailed balance. Since $A_{\omega }^a$ represents jumps with Bohr frequency $\omega$, the functional equation of $\gamma$ ensures a desired ratio of jumps with Bohr frequency $\omega$ and $-\omega$. We choose the Metropolis filter, $\gamma (\omega )=\min (1,e^{-\beta \omega })$, though another common filter $\gamma (\omega )=\frac{1}{1+e^{-\beta \omega }}$ for “Glauber dynamics” could also be used for the same results.

The Davies generator satisfies detailed balance with respect to ${𝝆}_{\beta }$, the thermal state. In some presentations of the Davies generator, it contains a coherent term $-i[𝑯,\cdot ]$. If this term is included, the generator does not satisfy detailed balance, so we do not follow this convention. However, the term does not affect the fixed point of the generator, since ${𝝆}_{\beta }$ commutes with $𝑯$ and therefore $-i[𝑯,{𝝆}_{\beta }]=0$.

The CKG Lindbladian is defined almost identically to the Davies generator, but is altered slightly so that it still obeys detailed balance, but is efficiently implementable.

where ${\widehat{𝑨}}^a(\omega )$ is now the Gaussian-supported OFT of jump operator ${𝑨}^a$:

To ensure that the jump operators do not have infinite precision in energy, a Gaussian supported OFT is performed instead to obtain ${\widehat{𝑨}}^a(\omega )$, which selects a Gaussian band energies of around $\omega$.

Here, $f(t)=e^{-{\sigma }_E^2{t}^2}\sqrt{{\sigma }_E\sqrt{2/\pi }}$, with Fourier transform $\widehat{f}(\omega )=\frac{1}{\sqrt{{\sigma }_E\sqrt{2\pi }}}\exp (-\frac{{\omega }^2}{4{\sigma }_E^2})$. As a result, the operator ${\widehat{𝑨}}^a(\omega )$ can be shown to be equal to ${\sum }_{\nu }\widehat{f}(\omega -\nu ){𝑨}_{\nu }^a$. The function $f(t)$ was chosen so that its squared Fourier transform ${\widehat{f}}^2(\omega )$ is a Gaussian with standard deviation ${\sigma }_E$, which features prominently in the Lindbladian (since it consists of quadratic terms in ${\widehat{𝑨}}^a(\omega )$). Taking ${\sigma }_E=\mathrm{\Theta }({\beta }^{-1})$ yields an efficient simulation algorithm with the assumption of a block-encoding of $𝑯$ and a block-encoding for the jump operators ${\sum }_{a\in [M]}|a⟩\otimes {𝑨}^a$, so ${\sigma }_E$ is taken to be on the order of ${\beta }^{-1}$ in this paper.

Since ${𝑨}^a(\omega )$ is a noisy decomposition of ${𝑨}^a$ into frequencies, it is not immediately clear whether there is a choice of function $\gamma (\omega )$ for which they can be recombined to achieve detailed balance. Indeed, as shown in [11], there is! The choice of $\gamma (\omega )$ is such that the transition part of ${ℒ}_{\beta }$, the summand ${\sum }_{a\in [M]}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\gamma (\omega ){\widehat{𝑨}}^a(\omega )(\cdot ){\widehat{𝑨}}^a{(\omega )}^{†}𝑑\omega$, still satisfies KMS detailed balance. [11] proved that there is a unique choice of $𝑩$, up to translation by a scalar, such that $-i[𝑩,\cdot ]-\frac{1}{2}{\sum }_{a\in [M]}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\gamma (\omega )\{{\widehat{𝑨}}^a{(\omega )}^{†}{\widehat{𝑨}}^a(\omega ),\cdot \}𝑑\omega$ also satisfies detailed balance. For the Davies generator, this coherent term $𝑩$ is simply $0$ (or corresponds to a Lamb shift that commutes with the Hamiltonian), and the decay term by itself already satisfies detailed balance. $𝑩$ can be expressed in general as:

The choice of $\gamma$ for this algorithm, for which the filter is efficiently implementable, is $\gamma (\omega )=\exp \left(-\beta \max (\omega +\frac{\beta {\sigma }_E^2}{2},0)\right)$. As ${\sigma }_E\to 0$ it converges to the Metropolis filter of the Davies generator. In particular, this $\gamma$ is precisely Metropolis filter for the Davies generator shifted by $\beta {\sigma }_E^2$. The CKG Lindbladian requires a choice of $\gamma$ for which $\alpha =\gamma \ast g$ satisfies the functional equation that was originally satisfied by $\gamma$ in the Davies generator.

We also note that ${ℒ}_{\beta }$ is bounded in operator norm.

With high probability at constant temperature, a randomly selected $d=\log (n)$-regular graph, with $\mathrm{poly}(d)$ random 1-design jumps, has a Lindbladian spectral gap of $\mathrm{\Omega }(d^{-3/4})$. This gives a polynomial algorithm to prepare the Gibbs state for most such graphs at constant temperature.

The gap of $\mathrm{\Omega }(d^{-3/4})$ arises because a random $d$-regular graph, for $d\to \mathrm{\infty }$, has one eigenvalue at $d$ and the rest distributed from $-2\sqrt{d-1}$ to $2\sqrt{d-1}$ in a distribution that converges to a (normalized) semicircle. This semicircular distribution frequently appears in random matrix theory, for instance in the Gaussian unitary ensemble (GUE), which models the spectrum of many chaotic quantum systems. When the spectrum of a quantum system indeed follows this distribution, it implies that $\delta (n)=\mathrm{\Omega }(d^{-3/4})$ of the eigenvalues lie within a constant of the minimum eigenvalue.

As an immediate result of Theorem 4.2 and Theorem 2.6, we obtain an algorithm polynomial in $d$ to prepare the Gibbs state of a $d$-regular graph. To prove the corollary, we use Theorem 2 in [17]. For this context, the following statement suffices.

In the above, ${\rho }_{sc}$ is the asymptotic distribution as $d\to \mathrm{\infty }$ of a random $d$-regular graph is the semicircular distribution with radius $2\sqrt{d-1}$, ${\rho }_{sc}(x)=\frac{1}{2\pi (d-1)}\sqrt{4(d-1)-x^2}.$ From this result, Theorem 4.2 follows. As described more explicitly in [32], the estimated density ${\rho }_{sc}$ near the bottom of the semicircle can be estimated. Then, using the theorem, this can be related to the fraction of eigenvalues near the minimum eigenvalue as well, proving the result.

We now mention another ensemble of Hamiltonians studied by [9] in the context of low-energy state preparation. In [9], efficient low energy state preparation with phase estimation is demonstrated under the same conditions as our efficient Gibbs sampling in Theorem 2.6. Indeed, if many eigenvectors are close to the ground-state energy, as we require, then performing phase estimation on the maximally mixed state has a high probability of measuring a low-energy state, so low-energy state preparation is possible as well. They study the following ensemble of Hamiltonians on $n_0$ qubits, ${𝑯}_{PS}={\sum }_{j=1}^m\frac{r_j}{\sqrt{m}}{𝝈}_j$, where $𝝈$ is a random Pauli string on $n_0$ qubits, each $r_j$ is sampled randomly from $\{-1,1\}$, and $m=\lfloor c_2\frac{n_0^5}{{ϵ}^4}\rfloor$. The parameter $ϵ$ satisfies $ϵ\ge 2^{-n_0/c_1}$, and $c_1,c_2$ are absolute constants. The resulting spectrum is again close enough to a semicircular distribution to obtain an efficient Gibbs sampler for certain temperatures that depend on $ϵ$. As $ϵ$ decreases, Gibbs sampling becomes efficient for even larger values of $\beta$ (lower temperatures), since the ensemble’s spectrum converges closer to a perfect semicircular distribution at the edge of the spectrum.

Using the results in their paper, we establish that Gibbs sampling is efficient in $n_0$ for certain values of $ϵ$ and corresponding temperatures ${\beta }^{-1}$.

For hypercube with varying dimension at a constant temperature, using unitary 1-design jumps would yield an exponentially large runtime. The spectrum of a hypercube with dimension $d$ and $2^d$ vertices consists of the integers $-d,-d+2\mathrm{\dots },d-2,d$. The eigenvalue $j$ has multiplicity $\left(\genfrac{}{}{0pt}{}{d}{\frac{d+j}{2}}\right)$. In particular, for any constant $C$, only an exponentially small fraction of the eigenvalues $\delta (d)$ lie below $-d+C$. This leads to a naive algorithm with at worst exponential complexity in $d$.

However, a better result can be obtained considering the hypercube as a system of $d$ qubits. The graph with dimension $d$ has $2^d$ vertices, which can be considered length $d$ bitstrings. Then, the adjacency matrix is the sum of Pauli $X$ operators on each qubit, ${\sum }_{i=1}^d{X}_i$, since the hypercube has an edge between any two bitstrings of Hamming distance 1. Choosing $d$ jump operators as $\frac{1}{\sqrt{d}}{Z}_a$, the mixing time can be improved to $\mathrm{poly}(d,\log ({ϵ}^{-1}))$: