Perfect Sampling for Quantum Gibbs States

TL;DR

This paper introduces a quantum algorithm for perfect sampling from quantum Gibbs states using a quantum adaptation of Coupling from the Past, with runtime depending on Hamiltonian degeneracy and stability under noise.

Contribution

It adapts classical perfect sampling algorithms to the quantum setting, enabling sampling without prior mixing time knowledge and analyzing its efficiency and noise stability.

Findings

Efficient for Hamiltonians with high degeneracy, polylogarithmic in dimension.

Runtime comparable to classical algorithms for non-degenerate spectra.

Algorithm is stable under implementation noise.

Abstract

We show how to obtain perfect samples from a quantum Gibbs state on a quantum computer. To do so, we adapt one of the `Coupling from the Past'-algorithms proposed by Propp and Wilson. The algorithm has a probabilistic run-time and produces perfect samples without any previous knowledge of the mixing time of a quantum Markov chain. To implement it, we assume we are able to perform the phase estimation algorithm for the underlying Hamiltonian and implement a quantum Markov chain such that the transition probabilities between eigenstates only depend on their energy. We provide some examples of quantum Markov chains that satisfy these conditions and analyze the expected run-time of the algorithm, which depends strongly on the degeneracy of the underlying Hamiltonian. For Hamiltonians with highly degenerate spectrum, it is efficient, as it is polylogarithmic in the dimension and linear in the mixing time. For non-degenerate spectra, its runtime is essentially the same as its classical counterpart, which is linear in the mixing time and quadratic in the dimension, up to a logarithmic factor in the dimension. We analyze the circuit depth necessary to implement it, which is proportional to the sum of the depth necessary to implement one step of the quantum Markov chain and one phase estimation. This algorithm is stable under noise in the implementation of different steps. We also briefly discuss how to adapt different `Coupling from the Past'-algorithms to the quantum setting.

Equations

$$[T(P_i),P_j]=0$$

$$\frac{\Vert \rho -\sigma {\Vert }_1}{2}=\underset{P\in {\mathcal{P}}_d}{\sup }\text{tr}[P(\rho -\sigma )].$$

$$\Vert T{\Vert }_{1\to 1}=\underset{X\in {\mathcal{M}}_d}{\sup }\frac{\Vert T(X){\Vert }_1}{\Vert X{\Vert }_1}.$$

$$t_{\text{mix}}=\min \{n\in \text{mathbbm}N:\underset{\rho \in D_d}{\sup }\Vert T^n(\rho )-\sigma {\Vert }_1\le 2e^{-1}\}.$$

$$T\left({\left(\frac{e^{-\beta H}}{{\mathcal{Z}}_{\beta }}\right)}^{\frac{1}{2}}X{\left(\frac{e^{-\beta H}}{{\mathcal{Z}}_{\beta }}\right)}^{\frac{1}{2}}\right)={\left(\frac{e^{-\beta H}}{{\mathcal{Z}}_{\beta }}\right)}^{\frac{1}{2}}{T}^{\ast }\left(X\right){\left(\frac{e^{-\beta H}}{{\mathcal{Z}}_{\beta }}\right)}^{\frac{1}{2}}$$

$$P(X_{n+1}=j|X_n=i)=\pi (i,j)$$

$$\pi (i,j)=\text{tr}\left(T(|{\psi }_i⟩⟨{\psi }_i|)|{\psi }_j⟩⟨{\psi }_j|\right).$$

$$\forall i,j\in S\exists n:{\pi }^n(i,j)>0.$$

$$\forall i\in S:\text{gcd}\{n\in \text{mathbbm}N\setminus \{0\}:{\pi }^n(i,i)>0\}=1.$$

$$\pi (i,j)\mu (i)=\mu (j)\pi (j,i).$$

$$\Vert \nu -\mu {\Vert }_1=\sum _{i\in S}|\mu (i)-\nu (i)|.$$

$$t_{\text{mix}}=\min \{n\in \text{mathbbm}N:\underset{\nu }{\sup }\Vert {\pi }^n\nu -\mu {\Vert }_1\le 2e^{-1}\}.$$

$$C=\min \{T|\forall i\in S\exists 1\le k\le T\text{ such that }{X}_k=i\}.$$

$$f(x)=i\iff x\in S_i.$$

$$\sum _{k\in S_k}\pi (l,k)=\sum _{k\in S_k}\pi (l^{\prime },k)$$

$$\text{tr}\left(T(|{\psi }_i⟩⟨{\psi }_i|)|{\psi }_j⟩⟨{\psi }_j|\right)=f(E_i,E_j)$$

$$\pi (i,j)=\text{tr}(T(|{\psi }_i⟩⟨{\psi }_i|)|{\psi }_j⟩⟨{\psi }_j|)$$

$$i\sim j\iff \text{tr}(H|{\psi }_i⟩⟨{\psi }_i|)=\text{tr}(H|{\psi }_j⟩⟨{\psi }_j|).$$

$$\tilde{\pi }(S_l,S_k)=|S_k|f(E_l,E_k).$$

$$\tilde{\mu }(S_i)=|S_i|\frac{e^{-\beta E_i}}{{\mathcal{Z}}_{\beta }},$$

$$\sum _{k\in S_k}\pi (l,k)=\sum _{k\in S_k}\pi (l^{\prime },k).$$

$$\tilde{\mu }(S_i)\tilde{\pi }(S_i,S_j)=|S_i|\frac{e^{-\beta E_i}}{{\mathcal{Z}}_{\beta }}|S_j|f(E_i,E_j).$$

$$\frac{e^{-\beta E_i}}{{\mathcal{Z}}_{\beta }}f(E_i,E_j)=f(E_j,E_i)\frac{e^{-\beta E_j}}{{\mathcal{Z}}_{\beta }}.$$

$$P(\text{algorithm terminates after at most }{N}_{ϵ}\text{ steps})\ge 1-ϵ.$$

$$P(X_0\ne Y)=P(X_0\ne Y|A_{ϵ})P(A_{ϵ})+P(X_0\ne Y|A_{ϵ}^C)P(A_{ϵ}^C).$$

$$G(k,j)=\{\begin{array}{ll}j& \text{mbox}ifk=0\\ 0& \text{mbox}otherwise\end{array}.$$

$$|S_j|\frac{e^{-\beta E_j}}{{\mathcal{Z}}_{\beta }}.$$

$${\rho }_j=\frac{1}{|S_j|}{P}_j.$$

$$\sigma (H)\subset \underset{i=1}{\overset{d^{\prime }}{⨆}}(e_i-ϵ,e_i+ϵ)$$

Full text

\frefformat

plainlemlemma\fancyrefdefaultspacing#1 \FrefformatplainlemLemma\fancyrefdefaultspacing#1 \frefformatplainthmtheorem\fancyrefdefaultspacing#1 \FrefformatplainthmTheorem\fancyrefdefaultspacing#1 \frefformatplaincorcorollary\fancyrefdefaultspacing#1 \Frefformatplaincorcorollary\fancyrefdefaultspacing#1 \frefformatplaindefidefinition\fancyrefdefaultspacing#1 \FrefformatplaindefiDefinition\fancyrefdefaultspacing#1 \Frefformatplaindefialgorithm\fancyrefdefaultspacing#1

PERFECT SAMPLING FOR QUANTUM GIBBS STATES

Daniel Stilck França

Department of Mathematics, Technische Universität München, Boltzmannstrasse 3

Garching, 85748, Germany

We show how to obtain perfect samples from a quantum Gibbs state on a quantum computer. To do so, we adapt one of the “Coupling from the Past”-algorithms proposed by Propp and Wilson. The algorithm has a probabilistic run-time and produces perfect samples without any previous knowledge of the mixing time of a quantum Markov chain. To implement it, we assume we are able to perform the phase estimation algorithm for the underlying Hamiltonian and implement a quantum Markov chain such that the transition probabilities between eigenstates only depend on their energy. We provide some examples of quantum Markov chains that satisfy these conditions and analyze the expected run-time of the algorithm, which depends strongly on the degeneracy of the underlying Hamiltonian. For Hamiltonians with highly degenerate spectrum, it is efficient, as it is polylogarithmic in the dimension and linear in the mixing time. For non-degenerate spectra, its runtime is essentially the same as its classical counterpart, which is linear in the mixing time and quadratic in the dimension, up to a logarithmic factor in the dimension. We analyze the circuit depth necessary to implement it, which is proportional to the sum of the depth necessary to implement one step of the quantum Markov chain and one phase estimation. This algorithm is stable under noise in the implementation of different steps. We also briefly discuss how to adapt different “Coupling from the Past”-algorithms to the quantum setting.

Keywords: quantum Gibbs states, perfect sampling, quantum algorithms

1 Introduction

Markov chain Monte Carlo methods are ubiquitous in science. They have a similar structure: the solution to a problem is encoded in the stationary distribution of a Markov chain that can be simulated. The chain is then simulated for a “long enough” time until the current state of the chain is “close enough” to a sample of the stationary distribution of interest.

It is expected that with the advent of quantum computers one could use similar methods to develop algorithms to simulate quantum many-body systems that do not suffer from the sign problem [1], and many quantum algorithms with this property were proposed  [2, 3, 4, 5, 6, 7, 8, 9].

However, as for classical Monte Carlo methods, it is in general difficult to obtain rigorous bounds on how long is “long enough”, as the huge literature dedicated to Markov chain mixing attests [10]. This prompted research on an algorithm that would “decide for itself” when the current state of the Markov chain is close to or is even a perfect sample of the stationary distribution, without any prior knowledge of the mixing properties of the chain.

One of the first algorithms to do so was developed in [11]. Later Propp and Wilson proposed the “Coupling from the Past” (CFTP)-algorithm [12] and showed how it can be applied to efficiently obtain perfect samples for the Ising model. They also showed how to sample perfectly from the stationary distribution of an unknown Markov chain for which we can only observe transitions in a subsequent paper [13] and many perfect sampling algorithms were developedaaaThe website http://dimacs.rutgers.edu/~dbwilson/exact/, maintained by Wilson, contains a comprehensive list of references concerning the topic and other related material.  since. There are also proposals of quantum speedups to these algorithms [14].

In this article, we will generalize some of these algorithms to get perfect samples from a Gibbs state of a Hamiltonian on a quantum computer. By this, we mean that we are able to perform any measurement and observe the same statistics for the outcomes as if we were measuring the actual Gibbs state.

To implement them, we will need to be able to perform the phase estimation algorithm [15] for the underlying Hamiltonian. Furthermore, we assume we can implement a quantum Markov chain that drives the system to the desired Gibbs state and fulfills certain assumptions, such as reversibility of the chain, which we elaborate on below. We comment on which of the current proposals to prepare Gibbs states on quantum computers may be adapted for our purposes.

Like it is the case for the classical algorithms, our quantum algorithms do not require any previous knowledge about the mixing properties of the quantum Markov chain. They “decide for themselves” when the current state of the system corresponds to a perfect sample of the target Gibbs state and their run-time is probabilistic.

We will focus on adapting the “voter CFTP” algorithm. For the classical voter CFTP the expected run-time is quadratic in the dimension of the system and linear in the mixing time. The run-time of our version will depend highly on the number of distinct eigenvalues of the Hamiltonian and the dimension of the eigenspaces. In the worst case, which corresponds to Hamiltonians with a non-degenerate spectrum, our version of this algorithm will turn out to have the same expected run-time as its classical counterpart, up to a logarithmic factor in the dimension. However, for Hamiltonians with degenerate spectra, our algorithm can be more efficient than the classical CFTP. In the case of Hamiltonians with an extremely degenerate spectrum, our algorithm can even have a run-time which is proportional to the time necessary to obtain approximate samples. We discuss how to explore this fact to obtain certifiably good samples efficiently for Hamiltonians whose spectrum can be “lumped together” into a small number of intervals. We also briefly discuss how to generalize other variations of CFTP.

These algorithms are stable under noise and we give bounds on their stability with respect to the implementation of the quantum Markov chain and the phase estimation steps. A potential advantage of the “voter CFTP” in comparison to other methods proposed in the literature is that it only requires the implementation of a quantum circuit of low depth a (potentially prohibitive) number of times and significant classical post-processing to obtain a perfect sample without any prior knowledge of the mixing time. Other methods require the implementation of a circuit of (potentially prohibitive) length one time to obtain an approximate sample, but only under the previous knowledge of or assumptions on the mixing time. Therefore, these algorithms are qualitatively different from ours, as our algorithm provides a certificate that we are obtaining good samples.

Another motivation for this work is to explore how coupling techniques for classical Markov chains may be applied and generalized to the quantum setting. These are one of the most useful tools to derive mixing times [10, Section 4.2] and lie at the heart of many perfect sampling algorithms, as their name already suggests. The fact that we still lack a notion of a quantum analog of a coupling is, therefore, one of the main technical hurdles to generalize many results in the theory of classical Markov chains and is by itself an interesting open problem which hopefully this work can shed some light on.

2 Preliminaries

2.1 Notation and basic concepts

We begin by introducing some basic concepts we will need and fixing the notation. Throughout this paper, $\mathcal{M}_{d}$ will denote the space of $d\times d$ complex matrices and $[d]=\{1,\ldots,d\}$. We will denote by $\mathcal{D}_{d}$ the set of $d$-dimensional quantum states, i.e. positive semidefinite matrices $\rho\in\mathcal{M}_{d}$ with trace $1$. We will call a Hermitian operator $H\in\mathcal{M}_{d}$ a Hamiltonian. We will always denote its spectral decomposition by $H=\sum_{i=1}^{d^{\prime}}E_{i}P_{i}$, with $P_{i}$ orthogonal projections. Here $d^{\prime}$ denotes the number of distinct eigenvalues of $H$. As we will see later, the expected run-time of the algorithm will depend more on this number than the dimension. The eigenspace corresponding to the energy level $E_{i}$ of $H$ will be denoted by $S_{i}$ and we will denote its dimension by $|S_{i}|$. When we write $H=\sum_{i=1}^{d}E_{i}|\psi_{i}\rangle\langle\psi_{i}|$ we mean that $\{\ket{\psi_{i}}\}_{i=1}^{d}$ is an orthonormal eigenbasis of $H$. For an inverse temperature $\beta>0$, we define $\mathcal{Z}_{\beta}=\text{tr}[e^{-\beta H}]$ to be its partition function and $e^{-\beta H}/\mathcal{Z}_{\beta}$ its Gibbs state. A quantum channel $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ is a trace preserving completely positive map. We will also refer to such a map as quantum Markov chain. A state $\sigma\in\mathcal{D}_{d}$ is a stationary state of $T$ if we have $T(\sigma)=\sigma$. The channel is called primitive if we have $\forall\rho\in\mathcal{D}_{d}:\lim\limits_{n\to\infty}T^{n}(\rho)=\sigma$ and $\sigma>0$. There is an equivalent spectral characterization of primitive quantum channels. A quantum channel is primitive if $\sigma>0$ is the only eigenvector corresponding to eigenvalues of modulus $1$ of the channel [16]. In particular, this implies that the property of being primitive is stable under small perturbations of the channel.

A collection of self-adjoint operators $\{F_{i}\}_{i\in I}$ is called a POVM (positive-operator valued measure) if the $F_{i}\in\mathcal{M}_{d}$ are all positive semidefinite and $\sum_{i\in I}F_{i}=\mathds{1}$. Here $\mathds{1}\in\mathcal{M}_{d}$ is the identity matrix. A state $\rho\in\mathcal{D}_{d}$ and a POVM $\{F_{i}\}_{i\in I}$ induce a probability distribution $p$ through $p(i)=\text{tr}[F_{i}\rho]$. All the algorithms we will discuss have as their goal to produce exact samples of the distribution $p$ generated by an arbitrary POVM in the case that $\rho$ is a Gibbs state.

The following class of quantum channels will be one of the backbones of the algorithms we will present later.

Definition 2.1** (Eigenbasis preserving quantum channels).**

A quantum channel $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ is called eigenbasis preserving for a Hamiltonian $H=\sum\limits_{i=1}^{d^{\prime}}E_{i}P_{i}$ and inverse temperature $\beta>0$ if we have that for all $i,j\in[d^{\prime}]$ the commutator

[TABLE]

vanishes and $T\left(\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}\right)=\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}$.

By the commutator property, we can model the dynamics under $T$ on states that commute with $e^{-\beta H}/\mathcal{Z}_{\beta}$ as a classical Markov chain. One should not take this condition to imply that the dynamics under $T$ are classical, as will become clear in subsection 3.2, where we present some examples of eigenbasis preserving channels. We will first suppose that we can implement these channels exactly, but will later relax this condition and discuss the influence of noise in section 4.1.

We will need some distinguishability measures for quantum states and channels and convergence speed measures for primitive quantum channels. One of the main ones is through the Schatten $1-$Norm $\|X\|_{1}=\text{tr}[|X|]$ for $X\in\mathcal{M}_{d}$. This is justified by the clear operational interpretation given by its variational expression [17, p. 404]. If we denote by $\mathcal{P}_{d}$ the set of orthogonal projections in $\mathcal{M}_{d}$, we have for $\rho,\sigma\in\mathcal{D}_{d}$

[TABLE]

That is, $\|\rho-\sigma\|_{1}/2$ expresses the maximal probability of correctly distinguishing two states $\sigma,\rho$ by a projective measurement. This norm induces the $1\to 1$ norm on operators $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$:

[TABLE]

As a measure of the convergence speed of a quantum channel, we define the $l_{1}$-mixing time threshold of a primitive quantum channel $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ with unique stationary state $\sigma$, which is given by

[TABLE]

We will say that a channel $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ satisfies detailed balance or is reversible with respect to $\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}$ if we have that

[TABLE]

holds for all $X\in\mathcal{M}_{d}$. Here $T^{*}$ is the adjoint of the channel with respect to the Hilbert-Schmidt scalar product. Satisfying detailed balance with respect to $\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}$ implies in particular that $\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}$ is a stationary state of the channel.

A crucial ingredient for our sampling algorithm is the phase estimation algorithm, discovered originally in [15]. There are now many variations of it [18, 19, 20] and it is still the subject of active research. We will neither discuss in detail how to implement it nor its complexity and refer to [20] for that. For our purposes, we will just suppose that for a given Hamiltonian $H$ acting on $\mathbbm{C}^{d}$ we may implement a unitary $U$ on $\mathbbm{C}^{d}\otimes(\mathbbm{C}^{2})^{\otimes t}$ for some $t\in\mathbbm{N}$ that acts as follows:

For $\ket{\psi_{i}}$ an eigenstate of a Hamiltonian $H$ with $H\ket{\psi_{i}}=E_{i}\ket{\psi_{i}}$ we have $U\ket{\psi_{i}}\otimes\ket{0}=\ket{\psi_{i}}\otimes\ket{E_{i}}$, where $\ket{E_{i}}$ is the binary expansion of $E_{i}$ in the computational basis of $(\mathbbm{C}^{2})^{\otimes t}$. We will first assume that we may implement $U$ exactly, but later discuss how imperfections in the implementation of the phase estimation algorithm influence the output of the sampling algorithm in section 4.1.

We will now fix some notation and terminology for classical Markov chains. A sequence $X_{0},X_{1},X_{2},\ldots$ of random variables taking values in a (finite) set $S$, referred to as the state space, is called a Markov chain if we have

[TABLE]

for a $|S|\times|S|$ matrix $\pi$. $\pi$ is called the transition matrix of the chain. We will always denote by $\pi$ the transition matrix of a Markov chain that should be clear from context. Most of the times it will be the one induced by an eigenbasis-preserving channel through

[TABLE]

A probability distribution $\mu$ on $S$ is called stationary if we have that $\pi\mu=\mu$. A Markov chain is said to be irreducible if

[TABLE]

It is aperiodic if

[TABLE]

Analogously to the quantum case, we say that the transition matrix $\pi$ satisfies detailed balance with respect to $\mu$ if

[TABLE]

Satisfying detailed balance again implies that $\mu$ is stationary. It is a well-known fact that if a Markov chain is aperiodic and irreducible there exists a unique stationary distribution $\mu$ such that for any other distribution $\nu$ on $S$ we have that $\lim\limits_{n\to\infty}\pi^{n}\nu=\mu$. We define the variational distance between probability distributions $\nu,\mu$ as

[TABLE]

With a slight abuse of notation, we will also denote the $l_{1}$-mixing time threshold in variation distance for a Markov chain by

[TABLE]

Let

[TABLE]

We will denote by $E_{i}(C)$ the expected time it takes to observe all states starting from $X_{0}=i$ and by $T_{C}=\max_{i\in S}E_{i}(C)$ the cover time of the chain. We refer to e.g. [10, Chapter 1] for a review of these basic concepts of Markov chains.

2.2 Lumpable Channels and Chains

Given an equivalence relation or equivalently a partition of the state space $S=\bigsqcup_{i=1}^{d^{\prime}}S_{i}$ of a Markov chain $X_{0},X_{1},X_{2},\ldots$, it is possible to define a new stochastic process in which the state space is given by the equivalence classes in the following way. Define the function $f:S\to[d^{\prime}]$ which maps a state to its equivalence class, that is

[TABLE]

This new stochastic process is then given by the random variables $f(X_{0}),(X_{1}),f(X_{2}),\ldots$. If this stochastic process is again a Markov chain for all possible initial probability distributions on $S$, the chain is said to be lumpable with respect to this equivalence relation. We refer to [21, Section 6] for more on lumpable chains. These are sometimes also called projective chains. The next Theorem gives necessary and sufficient conditions for lumpability.

Theorem 2.1** (Lumpable Chain).**

A necessary and sufficient condition for a Markov chain to be lumpable with respect to a partition $S=\bigsqcup_{i=1}^{d^{\prime}}S_{i}$ is that for every pair $S_{l},S_{k}$ we have for all $l,l^{\prime}\in S_{l}$

[TABLE]

Moreover, the transition probability between $S_{l}$ and $S_{k}$ in the lumpable chain is given by Eq. (5).

Proof: We refer to [21, Theorem 6.3.2] for a proof. ∎

In order to perfectly sample from Gibbs states with degenerate spectra, we will need the concept of a lumpable channel, which we introduce here.

Definition 2.2** (Lumpable channel).**

An eigenbasis preserving quantum channel $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ is called lumpable for a Hamiltonian $H=\sum_{i=1}^{d^{\prime}}E_{i}P_{i}$ at inverse temperature $\beta>0$ if it is reversible and there is a function $f:\mathbbm{R}\times\mathbbm{R}\to[0,1]$ such that

[TABLE]

for any unit vectors $\ket{\psi_{i}}\in P_{i}(\mathbbm{C}^{d})$ and $\ket{\psi_{j}}\in P_{j}(\mathbbm{C}^{d})$.

Here $P_{j}(\mathbbm{C}^{d})$ denotes the image of the projection. That is, the transition probabilities for eigenstates of $H$ depend only on their respective energies. Notice that as we demand that the quantum channel satisfies detailed balance, the classical transition matrix induced by the channel will also satisfy detailed balance. The definition of lumpable channels is motivated by the following lemma.

Lemma 2.1** (Lumping energy levels).**

Let $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a lumpable quantum channel for a Hamiltonian $H=\sum_{i=1}^{d}E_{i}|\psi_{i}\rangle\langle\psi_{i}|$ and inverse temperature $\beta>0$. Here $\{\ket{\psi_{i}}\}_{i=1}^{d}$ is an orthonormal eigenbasis of $H$. Define the classical Markov chain on $[d]$ with transition matrix given by

[TABLE]

and partition the state space according to their energy into $\{S_{1},\ldots,S_{d^{\prime}}\}$, that is, the equivalence relation on the state space is given by

[TABLE]

Then the Markov chain is lumpable with respect to this partition and the lumped chain has transition matrix

[TABLE]

Moreover, the stationary distribution of the lumped chain is given by

[TABLE]

where $|S_{i}|$ is the degeneracy of the energy level, and the chain satisfies detailed balance with respect to $\tilde{\mu}$.

Proof: It follows from Theorem 2.1 that it is sufficient and necessary for the chain to be lumpable that for $S_{l},S_{k}$ we have for all $l,l^{\prime}\in S_{l}$

[TABLE]

Eq. (9) holds for lumpable quantum channels, as we have $\sum\limits_{k\in S_{k}}\pi(l,k)=|S_{k}|f(E_{l},E_{k})$, which clearly only depends on the equivalence class of $l$. Therefore, we may define a Markov chain with respect to this partition and from Theorem 2.1 it follows that the transition matrix of the lumpable chain is given by (7). We will now show that it satisfies detailed balance with respect to $\tilde{\mu}$ and therefore $\tilde{\mu}$ is the stationary distribution of the chain. We have

[TABLE]

Now, as the original chain satisfied detailed balance, it holds that

[TABLE]

Plugging Eq. (11) into (10) we see that the lumpable chain satisfies detailed balance with respect to $\tilde{\mu}$. This implies that $\tilde{\mu}$ is the stationary distribution of the lumped Markov chain. ∎

When working with lumpable channels, $t_{\textrm{mix}}$ will always refer to the mixing time of the lumped chain. It is in general not clear how the mixing time of the lumpable chain relates to the mixing time of the original chain and this is a topic of current research. Surprisingly, the mixing time may even increase under lumping, as was shown recently in [22]. However, as is shown e.g. in [10, Lemma 12.8] and remarked in [22], important parameters that describe the convergence of a Markov chain, such as the spectral gap or Cheeger constant [10, Chapter 12], can only increase when lumping chains. This implies that the mixing time cannot increase significantly by lumping. For example, in the counterexample found in [22], lumping increases the mixing time by a factor of $\Theta(\log(d))$.

2.3 (Classical) Voter CFTP

We will briefly describe a perfect sampling algorithms based on CFTP for Markov chains introduced in [13], called “voter CFTP”. We mostly stick to their terminology and notation. The goal of this algorithm is to produce perfect samples of the stationary distribution $\mu$ of some Markov chain. One of the main advantages of this algorithm is that we only need to be able to observe valid transitions of this Markov chain to obtain perfect samples of the target distribution.

One should note that this is in general not the most efficient algorithm for perfect sampling [13], but arguably the simplest to understand. Besides the pedagogical motivation to present it, it turns out that this version is of interest in the quantum case, as we will see later. For this algorithm we suppose we have access to a randomized procedure $\mathbf{RandomSuccessor}:S\to S$ such that $P(\mathbf{RandomSuccessor}(i)=j)=\pi(i,j)$, where $\pi$ is a transition matrix having $\mu$ as a stationary measure. Let $G$ be a vertex-labeled graph with vertices $-\mathbbm{N}_{0}\times S$ and labels $S$. We will define the labels and edges as the algorithm runs and denote by $G(k,i)$ the label of the vertex $(k,i)$. Pseudocode for the algorithm is provided below in algorithm LABEL:algo:cftppropp .

One does not need to add the edge in step LABEL:step:addedgeclassical. This only helps to visualize the process. The expected run-time of this algorithm and its complexity, of course, depend on properties of $\mathbf{RandomSuccessor}:S\to S$. We will discuss these when we analyze the same questions for our algorithm in the quantum case. We now provide a proof that algorithm LABEL:algo:cftppropp indeed produces a perfect sample if it terminates almost surely.

Theorem 2.2**.**

Suppose algorithm LABEL:algo:cftppropp terminates with probability $1$ and denote the output by $Y$. Then $P(Y=i)=\mu(i)$.

Proof: Let $\epsilon>0$. As the algorithm terminates with probability $1$, there is a $N_{\epsilon}$ such that

[TABLE]

Denote by $A_{\epsilon}$ the event that the algorithm terminates after at most $N_{\epsilon}$ steps. Define a Markov chain $X_{-M},X_{-M+1},X_{-M+2},\ldots,X_{0}$ for some $M\in\mathbbm{N}$ and choose $X_{-M}$ according to $\mu$, i.e. $P(X_{-M}=i)=\mu(i)$. The transitions are defined by the graph, which we suppose was labeled for all $(k,i)$ with $k>-M$. Given $X_{k}=j$ we set $X_{k+1}=i$, where $\{(k,j),(k+1,i)\}$ is an edge of the graph $G$. As we chose $X_{-M}$ according to $\mu$ and $\mathbf{RandomSuccessor}$ has $\mu$ as a stationary distribution, $P(X_{k}=i)=\mu(i)$ for all $-M\leq k\leq 0$. We have

[TABLE]

One can check that the label on the graph at $(-M,i)$ is nothing but the value of $X_{0}$. Thus, if we assume that the algorithm has terminated, the value of $X_{0}$ does not depend on the initial value and will always be equal to $Y$. Therefore $P(X_{0}\not=Y|A_{\epsilon})=0$ if $-M\leq N_{\epsilon}$. Also, by construction, $P(X_{0}\not=Y|A_{\epsilon}^{C})P(A_{\epsilon}^{C})\leq\epsilon$. We then conclude $P(X_{0}\not=Y)\leq\epsilon$ and so the value of $Y$ and $X_{0}$ coincide, as $\epsilon$ was arbitrary. As $X_{0}$ is distributed according to $\mu$, so is $Y$.∎

3 CFTP for quantum Gibbs states

3.1 Voter CFTP

We will now show how to adapt the voter CFTP algorithm to quantum Gibbs state. We start by focusing on Hamiltonians that have a non-degenerate spectrum, as we need less assumptions in this case and the proof is simpler. We later generalize to arbitrary Hamiltonians. Given a Hamiltonian $H\in\mathcal{M}_{d}$ with non-degenerate spectrum, an eigenbasis preserving quantum channel $T$ for some inverse temperature $\beta>0$ and a POVM $\mathcal{F}=\{F_{i}\}_{i\in I}$, the following algorithm allows us to obtain perfect samples from the distribution $p(i)=\text{tr}\left[F_{i}\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}\right]$. The algorithm uses three registers corresponding to the tensor factors $\mathbbm{C}^{d}\otimes\left(\mathbbm{C}^{2}\right)^{\otimes t}\otimes\left(\mathbbm{C}^{2}\right)^{\otimes t}$, where $t$ is large enough to perform phase estimation for $H$ and tell apart the different eigenvalues of $H$. We will discuss how to choose $t$ in section 4.1. The first one will encode the current state of our system, while the other two will be used to record the output of two phase estimation steps. Define a labeled graph $G$ with vertices $V=-\mathbbm{N}_{0}\times\{1,\ldots,d\}$ and labels given by $\{0,\ldots,d\}$. We assume that $G$ has no edges at the beginning of the algorithm and the vertices are labeled as

[TABLE]

We assume we can prepare the maximally mixed state $\frac{\mathds{1}}{d}$. This can be done by picking a uniformly distributed integer between $1$ and $d$ and preparing the corresponding state of the computational basis, for example. We will assume that the Hamiltonian has a spectral decomposition given by $H=\sum_{i=1}^{d}E_{i}|\psi_{i}\rangle\langle\psi_{i}|$. The number $n$ denotes how many samples we wish to obtain in total and $c$ will denote a counter for the number of samples we still wish to obtain. The pseudocode for the algorithm is below in algorithm 3.1.

We now prove it indeed outputs perfect samples.

Theorem 3.1**.**

Let $T$ be a primitive, eigenbasis preserving quantum channel for a Hamiltonian $H$ and inverse temperature $\beta>0$. Then algorithm 3.1 terminates with probability 1 and generates $n$ perfect samples of the distribution $p$ defined above.

Proof: We will first show that with probability $1$ there is a $k\in-\mathbbm{N}$ and $l\in[d]$ such that $\forall i\in[d]$ $G(k,i)=l$. The probability that we observe an eigenstate $\ket{\psi_{i}}$ at step 6 is $\frac{1}{d}$, so with probability 1 we will observe it if we run the loop at step 3 often enough. This implies that we will assign a label different to [math] to arbitrary vertices of the graph $G$ if we run the while-loop at step 3 for long enough. Observe that as $T$ is an eigenbasis preserving quantum channel, the dynamics on the eigenbasis of $H$ under $T$ is just a classical Markov chain. As $T$ is primitive and the stationary state has full rank, this Markov chain is aperiodic and irreducible [16]. Because of that, the probability that we will obtain a $k$ such that $G(k,i)=l$ $\forall i\in[d]$ is $1$, using the same argument as the one given in [13] for the classical case. By the same argument, the probability that this label is $l$ is given by $\frac{e^{-\beta E_{l}}}{\mathcal{Z}_{\beta}}$, as this is the stationary distribution of the underlying classical Markov chain. As before, we will observe $\ket{\psi_{l}}\ket{E_{l}}\ket{0}$ at step 6 with probability 1 if we run the while-loop at step 3 often enough and this will then be a perfect sample by the previous discussion. ∎

Note that we could check if the measurement outcome we observe at step 13 is one of the desired outcomes to increase our probability of observing it.

Given how many distinct eigenvalues the Hamiltonian has and that we are able to implement a lumpable channel, we may run a modified version of algorithm 3.1 and obtain perfect samples. The steps of the algorithm are exactly the same and we do not write them out in detail. The only difference is the graph we feed the transitions to and what we feed. Let $d^{\prime}$ be again the number of distinct eigenvalues of $H$. In the case of degenerate Hamiltonians, we define a labeled graph $G$ with vertices $V=-\mathbbm{N}_{0}\times\{1,\ldots,d^{\prime}\}$ and labels given by $\{0,\ldots,d^{\prime}\}$. We assume that $G$ has no edges at the beginning of the algorithm and the vertices are labeled as

[TABLE]

That is, the graph is essentially the same as in the non-degenerate case but with $d^{\prime}$ instead of $d$ labels and vertices. At step 14 we then label the graph according to the energy levels we measured before at steps 6 and 13, as we can only tell apart states with different energies using phase estimation. We then have:

Theorem 3.2**.**

Let $T$ be a primitive, lumpable quantum channel for a Hamiltonian $H$ and inverse temperature $\beta>0$. Suppose further that $H$ has $d^{\prime}$ distinct eigenvalues. Then, if we run algorithm 3.1 with a graph modified as explained above, it terminates with probability 1 and generates $n$ perfect samples of the distribution $p(i)=\text{tr}\left(F_{i}\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}\right)$.

Proof: It should be clear that in this case the classical CFTP algorithm we are running based on the measurement outcomes will generate perfect samples from the stationary distribution of the lumped chain defined in Lemma 2.1. The convergence is guaranteed by the same argument as in the proof of Theorem 3.1. From Lemma 2.1 it follows that we will obtain the sample $S_{j}$ with probability

[TABLE]

Now, given that we have observed the label associated to $S_{j}$ after the first phase estimation step, we know that the state of the first register is given by

[TABLE]

Measuring $F_{i}$ on the outputs of the algorithm, therefore, gives perfect samples from the distribution $p$. ∎

3.2 Examples of eigenbasis preserving and lumpable channels

In order to run algorithm 3.1, we need to be able to implement a primitive eigenbasis preserving quantum channel for the Gibbs state we want to sample from in the case of non-degenerate spectrum and further that it is lumpable for the general case. In recent years many algorithms have been proposed to approximately prepare quantum Gibbs states on a quantum computer [2, 3, 4, 5, 6, 7, 8, 9]. We will here briefly discuss how some of them provide us with eigenbasis preserving or lumpable quantum channels for Gibbs states.

One class of eigenbasis preserving channels in the non-degenerate case are quantum dynamical semigroups with Davies generators. These are Markovian approximations for a quantum system weakly coupled to a thermal reservoir. A detailed description of the derivation and structure of Davies generators is beyond the scope of this article and can be found in [23, 24]. Under some conditions on the Hamiltonian and the coupling of the system to the bath, the Davies semigroup is primitive. The exact speed of this convergence is the subject of current research. We refer to [25] for a discussion of the conditions under which the Davies generators are primitive and some bounds on the convergence speed. In [9] Davies generators are proposed as a way of preparing thermal states on a quantum computer.

For our purposes, their main relevant property is that if the underlying Hamiltonian has a non-degenerate spectrum, the dynamics in the eigenbasis of the Hamiltonian does not couple diagonal terms to off-diagonal terms. They are therefore eigenbasis preserving. This was observed by many authors since the beginning of their study [24, 26, 27] and we refer to those for a proof of this claim.

More generally, it can be shown that dynamical semigroups that satisfy a quantum version of the detailed balance condition and whose stationary state has a non-degenerate spectrum are always eigenbasis preserving [28]. Our two previous examples fall into that category. This gives us a simple sufficient criterion to check whether a given implementation is eigenbasis preserving.

Note, however, that it is not a priori clear that a quantum dynamical semigroup can be implemented efficiently or by only using local operations. We refer to [9, 29] for a discussion of these topics.

An example of a lumpable channel is given by the implementation of the quantum Metropolis algorithm proposed in [7], as the quantum channel implemented at each step maps eigenstates of $H$ to eigenstates of $H$ and the transition probabilities are a function of their energy difference. However, it can be simplified for our purposes. As in the usual Metropolis algorithm, at each step we have to accept or reject a move that was made. One of the main difficulties to implement the quantum algorithm is reversing the evolution of the system if we reject the move. This is because, by the No-Cloning Theorem [30], we can’t make a copy of the previous state of the system. But the information that we rejected the move is enough for our algorithm, as we may simply copy the previous label when labeling the vertices. Therefore, we may skip the procedure of reversing the move.

3.3 Lumping Eigenstates together to obtain good samples

Until now we assumed we are able to implement phase estimation exactly and know the number of distinct eigenvalues of $H$. We may loosen this assumption and lump different eigenvalues together.

Definition 3.1** ($\epsilon$-Spectral Covering).**

Let $H\in\mathcal{M}_{d}$ be a Hamiltonian and $\epsilon>0$ be given. We call $\{e_{1},\ldots,e_{d^{\prime}}\}\subset\mathbbm{R}$ a $\epsilon$-Spectral Partition for $H$ if

[TABLE]

and for all $i\in[d^{\prime}]:\sigma(H)\cap(e_{i}-\epsilon,e_{i}+\epsilon)\not=\varnothing$. Here $\sigma(H)$ denotes the spectrum of $H$. We will refer to $d^{\prime}$ as the size of the covering.

It should be clear that $\epsilon$-spectral coverings are not unique and may have different sizes for fixed $\epsilon$. Although an $\epsilon$-spectral covering will not be readily available in most cases, there are some methods to obtain them. One can use e.g. the Gershgorin circle Theorem [31, Section VIII] to obtain a covering. If we can decompose $H$ into local commuting terms it is also possible to obtain an $\epsilon$-spectral covering by considering that the spectrum of $H$ must consist of sums of the eigenvalues of the local terms. Spectral coverings will be useful later to quantify the stability of algorithm 3.1 with respect to measuring the wrong energy with phase estimation. Here we will focus on showing how we may still obtain good samples based on an $\epsilon$-spectral covering and that the algorithm is stable w.r.t. introducing degeneracies into the spectrum because we can only obtain an estimate of it to a finite precision. Given an $\epsilon$-spectral covering for the Hamiltonian, we may run algorithm 3.1 with the number of labels being given by the size of the covering. If we use the Metropolis algorithm from [7] with the $e_{i}$ as the possible energies to define the transition probabilities we obtain:

Theorem 3.3**.**

Let $H\in\mathcal{M}_{d}$ be a Hamiltonian and $\{e_{1},\ldots,e_{d^{\prime}}\}\subset\mathbbm{R}$ be an $\epsilon$-spectral covering. Suppose we run algorithm 3.1 with this $\epsilon$-spectral covering as described above. Then the probability distribution $\tilde{p}$ of samples obtained from outputs of algorithm 3.1 satisfies

[TABLE]

Proof: The stationary distribution of the lumped chain will be

[TABLE]

with $\tilde{\mathcal{Z}}_{\beta}=\sum_{i=1}^{d^{\prime}}|\sigma(H)\cap(e_{i}-\epsilon,e_{i}+\epsilon)|e^{-\beta e_{i}}$. Let $\tilde{P}_{i}$ be the projection onto the subspace spanned by the eigenvectors of $H$ corresponding to eigenvalues in $\sigma(H)\cap(e_{i}-\epsilon,e_{i}+\epsilon)$. From the proof of Theorem 3.2, it follows that algorithm 3.1 will output the state

[TABLE]

We will now show $\left\|\tilde{\rho}-\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}\right\|_{1}\leq\sqrt{4\epsilon\beta}$, from which the claim again follows from the variational definition of the trace norm. From Pinsker’s inequality [32, Theorem 3.1], it follows that

[TABLE]

where $D\left(\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}||\tilde{\rho}\right)=\text{tr}\left(\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}(\log(\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}})-\log(\tilde{\rho})\right)$ is the relative entropy. As $\tilde{\rho}$ and

[TABLE]

commute, we have

[TABLE]

As we have an $\epsilon$-spectral covering, we have $E_{j}-e_{i}\leq\epsilon$ for $E_{j}\in\sigma(H)\cap(e_{i}-\epsilon,e_{i}+\epsilon)$ and $\frac{\tilde{\mathcal{Z}}_{\beta}}{\mathcal{Z}_{\beta}}\leq e^{\beta\epsilon}$. From this we obtain

[TABLE]

Plugging Eq. (20) into (18) we obtain the claim. ∎

This result may be interpreted as a first stability result. This shows that if we lump eigenvalues that are very close together, the Gibbs state does not change a lot. That is, if we introduce artificial degeneracies by not being able to tell apart eigenvalues that are very close through phase estimation this will not change the output of the algorithm significantly. As observed in [7], one could argue that a similar effect could in principle also affect classical Markov chain methods, as we are only able to compute the transition probabilities up to a finite precision. This does not seem to affect them in practice. Moreover, if we want samples that are certifiably at most $\delta$ apart in total variation distance at inverse temperature $\beta>0$, we may lump together eigenvalues that are at most $\delta^{2}/4\beta$ apart. As we will see later, high levels of degeneracy can reduce the run-time of the algorithm and this can be used to obtain good samples more efficiently.

4 Stability of the Algorithm

We will now address two possible sources of noise for algorithm 3.1 and show it is stable under these two. First, in the implementation of the channel and second in the phase estimation steps.

4.1 Stability in the implementation of the Channel

As shown in [33], one may quantify the stability of primitive quantum Markov chains with the following constant:

Definition 4.1**.**

Let $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a primitive quantum channel with stationary state $\sigma\in\mathcal{D}_{d}$. We define

[TABLE]

with $T_{\infty}(X)=\text{tr}(X)\sigma$.

We refer to [33] for bounds on it and how it can be used to quantify the stability of a quantum Markov chain with respect to different perturbations. Note that due to the spectral characterization of primitive quantum channels [34], the set of primitive quantum channels is relatively open in the convex set of quantum channels.

Theorem 4.1**.**

Let $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a primitive eigenbasis preserving channel for a Hamiltonian $H$ and inverse temperature $\beta>0$ and $T^{\prime}:\mathcal{M}_{d}\to\mathcal{M}_{d}$ a quantum channel satisfying

[TABLE]

for some $\epsilon>0$ small enough for $T^{\prime}$ to be primitive too. For a POVM $\{F_{i}\}_{i\in I}$, let $p$ and $p^{\prime}$ be probability distributions we obtain by measuring $\{F_{i}\}_{i\in I}$ on the output of algorithm 3.1 using $T$ and $T^{\prime}$ respectively. Then

[TABLE]

Proof: Let $\{P_{i}\}_{1\leq i\leq d^{\prime}}$ be the eigenprojections of $H$ and define $Q:\mathcal{M}_{d}\to\mathcal{M}_{d}$ to be the quantum channel given by

[TABLE]

Note that as $T$ is an eigenbasis preserving channel, $QTQ$ is an eigenbasis preserving channel with stationary state $\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}$. As $T^{\prime}$ is assumed to be primitive, $QT^{\prime}Q$ is primitive too, as $\|QT^{\prime}Q-QTQ\|_{1\to 1}\leq\|T^{\prime}-T\|_{1\to 1}$. Denote by $\rho$ the stationary state of the channel $QT^{\prime}Q$. By the variational expression for the trace distance, we have that

[TABLE]

From theorem 1 in [33] it follows that

[TABLE]

As $Q$ is a quantum channel, it follows that $\|Q\|_{1\to 1}\leq 1$ and so

[TABLE]

Eq. (22) would then follow from $\kappa(QTQ)\leq 2+\kappa(T)$. Note that as $T$ is primitive, we have that $\|T-T_{\infty}\|<1$, where we use the operator norm. Also, $QT_{\infty}Q=T_{\infty}$ and $Q$ is a projection. Thus

[TABLE]

As $T$ is an eigenbasis preserving channel, we have that

[TABLE]

and so

[TABLE]

We therefore have

[TABLE]

As $\|Q\|_{1\to 1},\|\text{id}\|_{1\to 1}\leq 1$ and from (27) we obtain

[TABLE]

[TABLE]

which completes the proof. ∎

Theorem 4.1 shows that the algorithm is stable under perturbations of the eigenbasis preserving channel. The stability for lumpable channels follows by observing that every lumpable channel is in particular eigenbasis preserving.

4.2 Faulty Phase Estimation

We will now analyze the errors stemming from faulty phase estimation. It is important to differentiate two different types of error that are caused by the phase estimation procedure. The first type of error comes from the fact that we are only able to obtain an estimate of the energy up to $t$ bits from the phase estimation procedure. This leads to round-off errors and may introduce degeneracies. As discussed in section 3.3 in theorem 3.3, algorithm 3.1 is stable against this sort of error. Moreover, as we will see later, this can even lead to the algorithm being more efficient.

The second kind of error comes from the fact that the phase estimation procedure only gives the correct energy of the state with high probability. This can cause some transitions we record to be corrupted. We will now show that algorithm 3.1 is stable against this kind of error.

We note that the exact distribution of the outcomes of the phase estimation procedure depends on which version is being used and this is still a topic of active research [20].

However, it is reasonable to assume that for any phase estimation routine the distribution of the outcomes will concentrate on the correct output. We will give bounds in terms of how large this peak is and explicit bounds for the implementation discussed in [17, Section 5.2]. We now assume we have some rule to assign the labels based on the outcome of the phase estimation step. In case we have an $\epsilon$-spectral covering, this might just be a function which assigns the label based on which interval of the covering the outcome of the measurement belongs to. Let $X_{1}\in[d^{\prime}]$ be the random variable which describes which label we assign to the graph after the measurement and $Y_{1}\in\{E_{1},\ldots,E_{d^{\prime}}\}$ the random variable which describes in which eigenspace the system finds itself after the first measurement at step 6. Let analogously $Y_{2}$ be the random variable which is distributed according to the probability of each eigenspace at step 13 and $X_{2}$ the second label which we assign. We will now assume that the errors stemming from the phase estimation steps are independent and have the same distribution. That is, given that the system is in a given eigenstate, the probability distribution of the measurement outcomes is the same in the two steps. Let the stochastic matrix $\Xi\in\mathcal{M}_{d^{\prime}}$ be given by

[TABLE]

Then, given that we have assigned the label $i$ to the graph after step 16 in algorithm 3.1, the state of the system is described by

[TABLE]

where $P_{j}$ is the projection onto the eigenspace corresponding to $E_{j}$. After we apply an eigenbasis preserving channel $T$, the state of the system is described by the state

[TABLE]

Furthermore, denote by $\Xi^{\prime}\in\mathcal{M}_{d^{\prime}}$ the stochastic matrix

[TABLE]

From Eq. (30) it then follows that the probability that the second label is $l$ given that the first label was $i$ is

[TABLE]

From Eq. (31) it is clear that the transition matrix for the labels is given by

[TABLE]

when we have faulty phase estimation.

As mentioned before, we expect $P(X_{1}=i|Y_{1}=E_{j})\simeq\delta_{i,j}$, that is, that the distribution peaks around the right outcome. To quantify this we define

[TABLE]

and $\xi^{\prime}$ analogously. We then have

Lemma 4.1**.**

Let $\xi$ and $\xi^{\prime}$ be defined as above and $\{F_{i}\}_{i\in I}$ a POVM. For a primitive lumpable channel $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ for a Hamiltonian $H$ and inverse temperature $\beta>0$, let $p(i)=\text{tr}\left(F_{i}\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}\right)$ and $p^{\prime}(i)$ be the probability of observing $F_{i}$ at the output of algorithm 3.1 with faulty phase estimation. Then

[TABLE]

Proof: As discussed in Eq. (32), the transition matrix for the observed energy labels is given by $\pi^{\prime}=\Xi^{\prime}\pi\Xi$. It easily follows from the definition of $\xi$ and $\xi^{\prime}$ that

[TABLE]

where $\tilde{\Xi}$ and $\tilde{\Xi}^{\prime}$ are again stochastic matrices. We may therefore write

[TABLE]

This transition matrix will still be primitive for $\xi$ and $\xi^{\prime}$ sufficiently large. Let $\mu$ be the stationary distribution of $\pi$ and $\mu^{\prime}$ the one of $\pi^{\prime}$. Observe that, as $\kappa(\pi)\leq\kappa(QTQ)$, we may use the bound $\kappa(QTQ)\leq 2+\kappa(T)$ from the proof of theorem 4.1 and obtain

[TABLE]

Here we have used that $\|\pi\|_{1\to 1}\leq 1$ for a stochastic matrix $\pi$. At the output of the algorithm, we would be measuring the POVM on the state $\rho^{\prime}=\sum_{i=1}^{d^{\prime}}\mu^{\prime}(i)\frac{P_{i}}{|S_{i}|}$ if no error occurs at step 6 of algorithm 3.1. But as an error might occur when we try to identify a given eigenstate, we will be measuring the state $\rho_{EM}=\sum_{i=1}^{d^{\prime}}(\Xi^{\prime}\mu^{\prime})(i)\frac{P_{i}}{|S_{i}|}$. By the definition of $\xi^{\prime}$, we have $\Xi^{\prime}\mu^{\prime}=\xi^{\prime}\mu^{\prime}+(1-\mu^{\prime})\tilde{\Xi^{\prime}}\mu^{\prime}$. We will measure the POVM on the state

[TABLE]

Here $\rho^{\prime\prime}$ is some density matrix. It then follows that

[TABLE]

The claim then follows from the variational expression for the trace distance as in the proof of theorem 4.1. ∎

Using Bayes’ rule it is possible to express the entries of the matrix $\Xi$ in terms of those of $\Xi^{\prime}$, which are more readily accessible. We have

[TABLE]

As the initial state is the maximally mixed one, we have that $P(Y_{1}=E_{j})=|S_{j}|d^{-1}$. From this discussion it follows that:

Theorem 4.2**.**

Let $\xi$ be defined as in Eq. (33). Then

[TABLE]

Proof: See the discussion above. ∎

This shows that the algorithm is stable if we do not have eigenvalues that we can misidentify with considerable probability and s.t. the degeneracy levels are of different order.

We now give estimates of $\xi$ and $\xi^{\prime}$ for the implementation of phase estimation considered in [17, Section 5.2] in case we have an $\epsilon$-spectral covering of the Hamiltonian or know that different eigenvalues are $\epsilon$ apart. In [17] it is shown that if we use

[TABLE]

qubits to perform phase estimation, then we obtain $E_{i}$ accurate to $n$ bits with probability at least $1-\delta$. This implies that $\xi^{\prime}\geq 1-\delta$. To estimate $\xi$ we need to control the terms of the form $P(X_{1}=i|Y_{1}=E_{j})$ for $i\not=j$. To this end we define

[TABLE]

with

[TABLE]

Lemma 4.2**.**

Let $H\in\mathcal{M}_{d}$ be a Hamiltonian and $\{e_{1},\ldots,e_{d^{\prime}}\}$ be an $\epsilon$-spectral covering of it. Suppose we implement phase estimation for $H$ using $t$ qubits. Then

[TABLE]

for $j\not=i$ and $\Delta_{i,j}$ defined as in Eq. (36).

Proof: In  [17, Section 5.2] it is shown that given that the eigenstate of the system is $E_{i}$, we have that the probability that the observed outcome is $E$ is bounded by

[TABLE]

For any point of the spectrum of $H$ in $(e_{i}-\epsilon,e_{i}+\epsilon)$ and for a point in $E\in(e_{j}-\epsilon,e_{j}+\epsilon)$ we have that $|2^{t}(e_{i}-E)\mod 2^{t}|^{-2}\geq\Delta(i,j)$. There are at most $2\epsilon 2^{t}+1$ possible outcomes that lie in the interval $(e_{j}-\epsilon,e_{j}+\epsilon)$. We therefore have

[TABLE]

Note that we have $2^{t}\leq\Delta(i,j)$, so the probability of misidentifying the labels goes to zero exponentially fast with the number of qubits for fixed $\epsilon$. We then obtain for $\xi$ and $\xi^{\prime}$:

Corollary 4.1**.**

Let $H\in\mathcal{M}_{d}$ be a Hamiltonian and $\{e_{1},\ldots,e_{d^{\prime}}\}$ be an $\epsilon$-spectral covering of it with $\epsilon\geq 2^{-n}$. Suppose we implement phase estimation for $H$ using $t\geq n+1+\log\left(2+(2\delta)^{-1}\right)$ qubits. Then $\xi^{\prime}\geq 1-\delta$ and

[TABLE]

Proof: As $\epsilon\geq 2^{-n}$ and with probability at least $1-\delta$ we will obtain an output which is accurate up to $n+1$ bits, with probability at least $1-\delta$ we will correctly identify in which element of the covering we are from the output. From this, it follows that $\xi^{\prime}\geq 1-\delta$. As the function $(x,y)\mapsto\frac{x}{x+y}$ is monotone increasing in $x$ and decreasing in $y$ for $x,y>0$, we obtain Eq. (37) by inserting the bound on $\xi^{\prime}$ and the result of Lemma 4.2 into the expression we derived for $\xi$ in theorem 4.2. ∎

Corollary 4.1 clarifies that the algorithm requires a larger number of qubits to be reasonably stable if we have close $E_{j}$ and $E_{j^{\prime}}$ s.t. $|S_{j}|\ll|S_{j^{\prime}}|$. The converse is also true; if $|S_{j}|\backsimeq|S_{j^{\prime}}|$ for all $j,j^{\prime}$ the algorithm is already stable with a small precision.

5 Expected run-time, Memory Requirements and Circuit Depth

We will now address the expected run-time of algorithm 3.1. To this end, we will only consider the number of calls of the phase estimation and eigenbasis preserving or lumpable channel and not the necessary classical post-processing, as we consider the quantum routines the more expensive resources.

In [13, Theorem 5], it was shown that the expected time to obtain a sample using algorithm LABEL:algo:cftppropp is $\mathcal{O}(t_{\textrm{mix}}|S|^{2})$ steps, where again $t_{\textrm{mix}}$ is the time such that the chain is $e^{-1}$ close to stationarity and $|S|$ the size of our state space. In the case of Hamiltonians with degenerate spectrum $t_{\textrm{mix}}$ will denote the mixing time of the classical lumped chain induced by the lumpable channel (see Eq. (7)).

Recall that $d^{\prime}$ denotes the number of distinct eigenvalues of the Hamiltonian or the size of the $\epsilon$-spectral covering being used. That is, $d^{\prime}$ is just the number of different labels of the graph. We will say that a column indexed by $k\in-\mathbbm{N}$ of $G$ is complete if $\forall i\in[d^{\prime}]$ $G(k,i)\not=0$. In [13], it is shown that we need to complete on average $\mathcal{O}(t_{\textrm{mix}}d^{\prime})$ columns of the graph $G$ before the labels on a column become constant. As each step to complete a column needs $\mathcal{O}(d^{\prime})$ calls of $\mathbf{RandomSuccessor}$, this leads to a total of $\mathcal{O}(t_{\textrm{mix}}{d^{\prime}}^{2})$ calls of $\mathbf{RandomSuccessor}$. The dynamics in the eigenbasis of $H$ is classical, so we may use the exact same reasoning to conclude that we will need an expected $\mathcal{O}(t_{\textrm{mix}}d^{\prime})$ number of complete columns until we obtain one perfect sample.

But in our case, we may need more uses of the channel and phase estimation, as we may not prepare an arbitrary eigenstate of $H\in\mathcal{M}_{d}$ which might be necessary to complete a column deterministically. We will denote the expected number of measurements necessary to complete a column by $\phi(H)$ and in theorem A.1 in the Appendix A we give an explicit expression for this quantity.

In Appendix A we prove bounds on $\phi(H)$ for various cases of interest and remark that in the worst case, namely Hamiltonians with a non-degenerate spectrum, $\phi(H)=\mathcal{O}\left(d\log(d)\right)$. Preparing the initial states probabilistically does not significantly change the overall efficiency of the algorithm, as illustrated by the next theorem.

Theorem 5.1**.**

Let $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a lumpable quantum channel for a Hamiltonian $H$ at inverse temperature $\beta>0$ with mixing time $t_{\textrm{mix}}$. Then the expected number of steps until algorithm 3.1 returns a perfect sample is $\mathcal{O}(t_{\textrm{mix}}d^{\prime}\phi(H))$.

Proof: We will need an average of $\phi(H)$ measurements to complete a column. From the result [13, Theorem 5] we know that we will need an expected number of $\mathcal{O}(t_{\textrm{mix}}d^{\prime})$ number of complete columns to obtain a sample. As the number of measurements needed to complete a column and complete columns to obtain a sample are independent, we have an expected $\mathcal{O}(t_{\textrm{mix}}d^{\prime}\phi(H))$ number of steps to obtain a sample.∎

It should be clear from theorem 5.1 that algorithm 3.1 is considerably less efficient than other algorithms such as quantum Metropolis [7] if we are willing to settle for an approximate sample for Hamiltonians with a non-degenerate spectrum. In this case we have $d^{\prime}=d$ and $\phi(H)=\mathcal{O}\left(d\log(d)\right)$, giving a total complexity of $\mathcal{O}(t_{\textrm{mix}}d^{2}\log(d))$ in the worst case. After all, to obtain a sample that is $e^{-1}$ close in trace distance to the Gibbs state, one only needs $t_{\textrm{mix}}$ steps of the Metropolis algorithm instead of the $\mathcal{O}(t_{\textrm{mix}}d^{2}\log(d))$ needed for CFTP. Therefore, it is important to stress again that these algorithms are very different in nature. Algorithm 3.1 provides us with perfect, not approximate samples, and it is the first algorithm of this form for quantum Gibbs states to the best of our knowledge. It provides a certificate that we are indeed sampling from the right distribution when it terminates, while most other algorithms require some mixing time bounds to obtain a sample that can be considered close to the target distribution. Moreover, it only requires us to be able to implement one step of the chain.

However, for the case of Hamiltonians with a highly degenerate spectrum, our algorithm is efficient, as is illustrated by the next theorem:

Theorem 5.2**.**

Let $T:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a lumpable quantum channel for a Hamiltonian $H\in\mathcal{M}_{d}$ at inverse temperature $\beta>0$ with mixing time $t_{\textrm{mix}}$. Moreover, assume that we have $d\leq|S_{i}|r(d)$ for some function $r:\mathbbm{R}\to\mathbbm{R}$ and all eigenspaces $S_{i}$. Then the expected number of steps until algorithm 3.1 returns a perfect sample is $\mathcal{O}(t_{\textrm{mix}}r(d)^{2}\log(r(d)))$.

Proof: From $d\leq|S_{i}|r(d)$ it follows that $d^{\prime}\leq r(d)$. It follows from theorem A.2 that we will need an average of $\mathcal{O}\left(r(d)\log(r(d))\right)$ measurements to complete a column. From the result [13, Theorem 5] we know that we will need an expected number of $\mathcal{O}(t_{\textrm{mix}}d^{\prime})$ number of complete columns to obtain a sample. As the number of measurements needed to complete a column and complete columns to obtain a sample are independent, we have an expected $\mathcal{O}(t_{\textrm{mix}}r(d)^{2}\log(r(d)))$ number of steps to obtain a sample. ∎

In particular, for the cases $r(d)=c$ for some $c\in\mathbbm{R}$, which corresponds to having eigenspaces with a degeneracy proportional to the dimension, we have that we only need $\mathcal{O}\left(t_{\textrm{mix}}\right)$ steps to obtain a perfect sample. That is, the time necessary to obtain perfect samples with our algorithm and approximate ones are the same up to a constant factor. Slightly more generally, for $r(d)=c\log(d)^{m}$ our algorithm still has a polylogarithmic runtime and is efficient. Admittedly such level of degeneracy is not usual for Hamiltonians of physical relevance. One could use the strategy discussed in section 3.3 and still obtain certifiably good samples by lumping together eigenvalues that are close. Moreover, as we will only need to run a $r(d)$ dimensional version of classical CFTP, the classical part of the algorithm will be efficient.

Although the worst case $\mathcal{O}(t_{\textrm{mix}}d^{2}\log(d))$ scaling is prohibitive for applications, this is still more efficient than explicitly diagonalizing $H$ as long as $t_{\textrm{mix}}\log(d)=\mathcal{O}(d^{\omega-2})$. Here $2<\omega<2.373$ is the optimal exponent of matrix multiplication, which has the same complexity as diagonalization [35]. That is, as long as approximate sampling is efficient, obtaining perfect samples is faster than diagonalizing even in the worst case.

We now analyze the circuit depth and memory requirements to obtain a sample.

Theorem 5.3**.**

Let $C_{PT}$ and $C_{T}$ be the circuit depth needed to implement the phase estimation for $H$ and the eigenbasis preserving channel, respectively. Then one needs to implement a quantum circuit of depth $\mathcal{O}(C_{PT}+C_{T})$ to obtain a sample and moreover an expected $\mathcal{O}(\phi(H))$ classical memory.

Proof: The circuit length part follows easily from just going through the steps of algorithm 3.1, as to label the new vertex we need to implement two phase estimation steps and apply the eigenbasis preserving channel once.

To see the that we only need $\mathcal{O}(\phi(H))$ classical memory, notice that we only need to store the information contained in the last complete column to perform the later steps. This is because it contains all possible labels for future columns. By corollary A.2, we have that the expected number of labels we obtain before completing a column is $\mathcal{O}(\phi(H))$, and so we need a total classical memory of size $\mathcal{O}(\phi(H))$. ∎

The quantum part of algorithm 3.1 can be easily parallelized, as we could use different quantum computers feeding a classical computer with valid transitions. Note that the classical resources necessary to run the algorithm are also not very large in the cases in which we have a highly degenerate spectrum, as discussed before.

6 Adapting other Variations of CFTP

In [13] the authors discuss other variations of CFTP that can be more efficient, such as the cover time CFTP algorithm. We will not discuss in detail how to adapt these other proposals, but it should be straightforward to do so from the results in the last sections. In this section, we will just mention the main ideas. Note that the only thing necessary to implement all these variations is a valid $\mathbf{RandomSuccessor}$ function and the outputs of the measurements in steps 6 and 13 of algorithm 3.1 do exactly that. This information could then be fed to a classical computer running a variation of CFTP. The only difference to the classical case is that we may not choose arbitrary initial states, but do so probabilistically. However, by waiting until each initial state is observed, we may circumvent this and do not have a significant overhead by the result of corollary A.2.

For some variations of CFTP, like again the cover time CFTP, one needs to iterate $\mathbf{RandomSuccessor}$. This is also straightforward. If we want to obtain a given number of iterations of $\mathbf{RandomSuccessor}$, we just apply an eigenbasis preserving or lumpable channel $T$ to the first register, repeated by a phase estimation step and a measurement in the computational basis. We then repeat this procedure to obtain the iterations.

One could then repeat the analysis done in this section and see that the run-time is again of the same order of magnitude as the classical version of the CFTP algorithm and obtain a perfect sampling algorithm with a run-time proportional to the cover time of the lumped chain.

7 Conclusion and Open Problems

We have shown how to adapt perfect sampling algorithms for classical Markov chains to obtain perfect samples of quantum Gibbs states on a quantum computer. These algorithms have an average run-time which in the worst case is similar to their classical counterparts. For highly degenerate Hamiltonians this algorithm gives an efficient sampling scheme and in the extreme case of having degeneracies proportional to the dimension, the time required to sample perfectly is even proportional to the time necessary to obtain an approximate sample. In these cases, the classical post-processing required can be done efficiently. We showed how to increase the efficiency of the sampling scheme and still obtain certifiably good samples by lumping close eigenstates together. We argue that one of its main advantages is its short circuit depth. We show that the algorithm is stable under noise in different steps of the implementation. It would be interesting to find sampling applications or models that satisfy the conditions under which our algorithms are efficient. Moreover, it would be worthwhile to investigate if there is a class of models to which we can tailor the perfect sampling algorithms to be efficient, as was done with success for attractive spin systems [12].

Acknowledgments

D.S.F. would like to thank Michael M. Wolf for raising the question of how to perfectly sample from a quantum Gibbs state and the constant support and Andreas Bluhm for helpful discussions. D.S.F acknowledges support from the graduate program TopMath of the Elite Network of Bavaria, the TopMath Graduate Center of TUM Graduate School at Technische Universität München, the Deutscher Akademischer Austauschdienst(DAAD) and by the Technische Universität München – Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement no. 291763.

References

Appendix A Expected number of observations to sample all possible outputs

To estimate the expected run-time of algorithm 3.1, we need to determine how often, on average, we must measure projections $\{P_{i}\}_{1\leq i\leq d^{\prime}}$ on the state $\frac{\mathds{1}}{d}$ until we observe all possible outcomes $i$. Here the $P_{i}$ correspond to projections onto different eigenspaces of the Hamiltonian Gibbs state $\frac{e^{-\beta H}}{\mathcal{Z}_{\beta}}$ we are trying to sample from. The number of measurements corresponds to the time necessary to complete a column in algorithm 3.1 and will be denoted by $\phi(H)$.

One can see that this corresponds to the classical problem of determining how many coupons one must collect to obtain at least one of each, the coupon collector problem [36]. But in our case, we have unequal probabilities for the coupons or outcomes, as the probability $q(i)$ of observing $i$ is given by

[TABLE]

Theorem A.1** (Coupon collector with unequal probabilities).**

Let $q\in\mathbbm{R}^{d^{\prime}}$ be the probability distribution of the measurement outcomes as in Eq. (38). Let $Y$ be the random variable given by the number of measurements necessary to observe all possible outputs. Then

[TABLE]

where $Q_{J}$ is defined as $Q_{J}=\sum_{i\in J}q(i)$ for $J\subset[d^{\prime}]$ and the second sum is over all subsets of size $j$.

Proof: We refer to [37, Corollary 4.2] for a proof. ∎

Although the expression in Eq. (39) is exact, its asymptotic scaling is not clear from it. Therefore, we show the following bound which is more directly accessible.

Theorem A.2**.**

Let $q\in\mathbbm{R}^{d^{\prime}}$ be the probability distribution of the measurement outcomes as in Eq. (38). Let $Y$ be the random variable given by the number of measurements necessary to observe all possible outputs. Moreover, assume that $d\leq|S_{i}|r(d)$ for some $r(d):\mathbbm{R}\to\mathbbm{R}$. Then

[TABLE]

where $\varphi_{d^{\prime}}=\sum\limits_{l=1}^{d^{\prime}}\frac{1}{l}$.

Proof: We mimic the proof of the classical result for the coupon collector problem with uniform probability distribution. Denote by $t_{l}$ the expected time to collect a new coupon after $l-1$ have been collected. We then have

[TABLE]

We clearly have $E(t_{1})=1$. Define $k(d)=\frac{d}{r(d)}$. For $l\geq 2$, note that as we have that we get each coupon with probability at least $\frac{k(d)}{d}$, the probability of getting a new coupon after having collected $l-1$ is at least $(d^{\prime}-l+1)\frac{k(d)}{d}$. From this, it follows that $E(t_{l})\leq\frac{d}{k(d)}\left(d^{\prime}-l+1\right)^{-1}$ and so

[TABLE]

The claim follows from observing that $\varphi_{d}=\mathcal{O}(\log(d))$ and that $d^{\prime}\leq r(d)$. ∎

From this, it is easier to get estimates for cases that might be of interest. Here we collect the bounds for the extreme cases of highly degenerate spectra, that is, with each eigenspace having dimension $\Omega\left(d\log(d)^{-m}\right)$ for $m\in\mathbbm{N}$ and the non-degenerate case.

Corollary A.1**.**

Let $q\in\mathbbm{R}^{d^{\prime}}$ be the probability distribution of the measurement outcomes as in Eq. (38). Let $Y$ be the random variable given by the number of measurements necessary to observe all possible outputs. Moreover, assume that $|S_{i}|\geq c\frac{d}{\log(d)^{m}}$ for some $c\in\mathbbm{R}$ and all $i\in[d^{\prime}]$. Then

[TABLE]

Proof: Just take $r(d)=c^{-1}\log(d)^{m}$ in theorem A.2. ∎

Corollary A.2**.**

Let $q\in\mathbbm{R}^{d^{\prime}}$ be the probability distribution of the measurement outcomes as in Eq. (38). Let $Y$ be the random variable given by the number of measurements necessary to observe all possible outputs. Moreover, assume that $q(i)=\frac{1}{d}$. Then

[TABLE]

Proof: Just take $r(d)=d$ in theorem A.2. ∎

That is, we might go from a constant number of samples necessary to complete a column in the case of degeneracies proportional to the dimension to a scaling like $d\log(d)$. One should note that applying the bound in corollary A.2 to analyze the runtime of algorithm 3.1 probably leads to bounds that are too pessimistic for spectra that are not very degenerate. To see why this is the case, note that in algorithm 3.1 we do not discard measurements outcome we have already observed, but rather use them to complete other columns, which we do not take into account in this analysis.