Skip to content

Gibbs sampling

Rough overview (in words)

Gibbs sampling is the task of preparing a quantum state in thermal equilibrium. This task is interesting in its own right as a means of testing the thermodynamic properties of quantum systems in a controlled way, but it is also a subroutine that is surprisingly useful within other quantum algorithms. Formally, given a Hamiltonian and a temperature, the task is to prepare the Gibbs state (also known as the thermal state) of that Hamiltonian at the associated temperature, or equivalently, to sample eigenstates of the Hamiltonian with probability proportional to their Boltzmann weights (motivating the name Gibbs sampling).

Physically, Gibbs sampling is routinely achieved in experiments via cooling as a manifestation of open-system thermodynamics, although theoretical understanding of such processes has been largely heuristic. Computationally, quantum Gibbs sampling is the quantum analogue of the same classical task in the computational basis, often achieved by Markov chain Monte Carlo (MCMC) methods. As a representative example, the Metropolis–Hastings algorithm [1] uses rejection sampling to construct a Markov chain whose stationary state is the classical Gibbs distribution; the Gibbs distribution can be efficiently sampled if the Markov chain mixes rapidly. Nowadays, Monte Carlo methods have already surpassed their original intent (Ising model simulation) and found ubiquitous applications in optimization and machine learning due to their simplicity and robustness. It is natural to wonder if the same features will be present for quantum Gibbs sampling.

Surprisingly, quantum algorithms and theoretical understanding of Gibbs sampling are severely underdeveloped. The most direct quantum algorithms for Gibbs sampling suffer from an explicit cost exponential in the size of the system. Another approach is to quantize classical Monte Carlo algorithms [2], but this approach has faced serious technical challenges rooted in quantum mechanics: the energy-time uncertainty principle (for imposing the Boltzmann weights) and the no-cloning theorem (for "rejecting" a quantum state). Recently, a new wave [3, 4, 5, 6] of proposals revisits the issue from the angle of open-system thermodynamics and gives nature-inspired algorithms for Gibbs sampling. These more directly emulate the dynamical process of thermalization and have the potential to achieve better runtimes for specific systems where thermalization is expected to be fast.

Rough overview (in math)

Given a Hamiltonian \(H = \sum_i E_i\ket{\psi_i}\bra{\psi_i}\) over \(n\) qubits, a desired inverse temperature \(\beta\), and an error parameter \(\epsilon\), the Gibbs sampling task is to prepare an \(n\)-qubit quantum state \(\rho\) such that

\[\begin{align} \nrm{ \rho - \sigma_{\beta} }_{\mathrm{tr}} \le \epsilon \quad \text{where}\quad \sigma_\beta :=\frac{ \mathrm{e}^{-\beta H }}{\mathcal{Z}}\propto \sum_i \mathrm{e}^{-\beta E_i} \ket{\psi_i}\bra{\psi_i} \quad \text{and} \quad \mathcal{Z} := \text{tr}[\mathrm{e}^{-\beta H }]. \end{align}\]

The above uses the convenient error metric given by the trace norm \(\lVert \cdot \rVert_{\mathrm{tr}}\), which controls the error for arbitrary bounded (possibly nonlocal) observables. In some applications, it could be sufficient to give a state \(\rho\) that approximates all local observables up to high precision, even if the global distance between \(\rho\) and \(\sigma_\beta\) is large. Note that \(\sigma_\beta\) corresponds to an ensemble of eigenstates of \(H\), where an eigenstate with energy \(E_i\) occurs with probability proportional to the Boltzmann weight \(e^{-\beta E_i}\).

To solve this problem, the quantum algorithm requires access to \(H\), for example, through a block-encoding of \(H\). Block-encodings can often be efficiently constructed, for instance, when \(H\) is a sparse matrix or when \(H\) is given as a sum of \(\text{poly}(n)\) local interaction terms. Henceforth, assume that \(H\) is offset such that it is guaranteed to be a nonnegative operator (no negative energies).

An early approach [7] for Gibbs sampling relied on quantum phase estimation (QPE) and amplitude amplification. In particular, one starts with a \(2n\)-qubit maximally entangled state (for which the reduced density matrix on the first \(n\) qubits is the maximally mixed state) and applies QPE to the first \(n\) qubits, reading an estimate of the energy into an ancilla register. Under the simplification that QPE has perfect resolution, one now has the state

\[\begin{equation} \frac{1}{\sqrt{2^n}} \sum_{i} \ket{\psi_i} \ket{\phi_i} \ket{E_i} \end{equation}\]

where \(\ket{\psi_i}\) is the \(i\)th eigenstate of \(H\), \(E_i\) is the associated energy, and the states \(\ket{\phi_i}\) form an arbitrary (unimportant) orthonormal basis. Next, one coherently rotates an ancilla qubit to put the correct Boltzmann weight into the amplitude:

\[\begin{equation} \frac{1}{\sqrt{2^{n}}}\sum_{i} \ket{\psi_i} \ket{\phi_i} \ket{E_i}\left(e^{-\beta E_i/2}\ket{0} + \sqrt{1-e^{-\beta E_i}}\ket{1}\right)\,. \end{equation}\]

Note that the probability of measuring the final qubit in \(\ket{0}\) is precisely \(\mathcal{Z}/2^n\). Rather than measure and postselect, one now performs amplitude amplification on the ancilla being \(\ket{0}\) to produce

\[\begin{equation} \frac{1}{\sqrt{\mathcal{Z}}} \sum_i e^{-\beta E_i/2}\ket{\psi_i}\ket{\phi_i}\ket{E_i} \end{equation}\]

up to small error, which is a purification of the Gibbs state \(\sigma_\beta = \mathcal{Z}^{-1}\sum_i e^{-\beta E_i} \ket{\psi_i}\bra{\psi_i}\). While QPE does not exactly produce the operation described above, a more complete analysis in [7, 8] shows the idea still works. This approach is akin to classical rejection sampling (see also [9]), where a state is chosen at random and accepted with probability \(e^{-\beta E_i}\), such that repeating until acceptance yields a sample from the correct distribution. Due to amplitude amplification, the quantum algorithm enjoys a quadratic speedup.

More advanced methods that have exponentially better \(\epsilon\) dependence have since been developed. Reference [10] used a linear combination of unitaries approach to perform the imaginary time evolution operator \(e^{-\beta H}\), again followed by amplitude amplification. Technically, that work assumed access to an operator similar to \(\smash{\sqrt{H}}\), but this requirement was removed in Gibbs samplers appearing in [11, 12], which employ a method for implementing smooth Hamiltonian functions. Alternatively, one can use the quantum singular value transformation along with a polynomial approximation to the function \(e^{-\beta(1-x)/2}\) on the interval \(x \in [-1,1]\) [13, Section 5.3] and combine this with (fixed-point) amplitude amplification [14].

Another family of quantum algorithms is closer in spirit to classical Monte Carlo methods. They quantize the Metropolis–Hastings algorithm (quantum Metropolis sampling [2]) or simulate the dynamics arising from a system-bath interaction [3, 4, 5, 6]. These algorithms make fundamental usage of quantum phase estimation for probing the energy, but most importantly (and most nontrivially), they construct a detailed-balance "quantum Markov chain" via either discretely or continuously "rejecting" the quantum state. Care must be taken to perform the rejection step coherently and to handle the fact that the energies cannot be learned to infinite precision. Abstractly, Monte Carlo–style quantum algorithms emulate a discrete quantum channel (or a continuous Lindbladian) that converges to the Gibbs state after \(\ell\) iterations

\[\begin{align} \mathcal{N}[\sigma_\beta] \approx \sigma_\beta\quad \text{and} \quad \nrm{ \mathcal{N}^{\ell}[\rho_{0}] - \sigma_{\beta}}_{\mathrm{tr}} \le \epsilon \end{align}\]

for some initial state \(\rho_0\). Like the classical Metropolis–Hastings algorithm, for some systems, the number of iterations \(\ell\) for convergence can be exponentially large (or worse) in \(n\), while for other systems, the number of iterations needed can be much smaller. It is a generally difficult problem to determine \(\ell\), but it is expected that the size of \(\ell\) will be related to the natural thermalization rate of the system. Note that such a process can be further quantized to gain quadratic speedup [15, 6].

Dominant resource cost (gates/qubits)

Assuming one has access to a block-encoding of the Hamiltonian \(H\), that is, a unitary whose upper left block is the operator \(H/\alpha\), where \(\alpha\) is a normalization constant at least as large as the spectral norm of \(H\), one can accomplish the Gibbs sampling task using [11, Lemma 44] (see also [12, Corollary 16])

\[\begin{equation} \label{eq:Gibbs_complexity} \alpha\beta \sqrt{\frac{2^n}{\mathcal{Z}}} \cdot \text{poly}(\log(1/\epsilon), n) \end{equation}\]

calls to the block-encoding and a similar number of other gates. Note that since we have assumed \(H\) is non-negative, we have \(\mathcal{Z}\leq 2^n\). In the case that \(H\) is \(d\)-sparse, we can take \(\alpha = d\). In the case one has access to \(\sqrt{H}\), the \(\beta\) dependence can be reduced from \(\beta\) to \(\sqrt{\beta}\) [10]. This complexity statement might be regarded as a quadratic speedup compared to the classical method of rejection sampling, which requires \(2^n/\mathcal{Z}\) samples on average; however, note that this classical method only directly applies to diagonal (classical) Hamiltonians \(H\). Otherwise, a classical approach may need to resort to exact diagonalization of \(H\), which has \(\mathcal{O}\left( 2^n \right)\) space complexity and even worse time complexity.

Monte Carlo–style quantum Gibbs sampling algorithms have complexity determined by

\[\begin{equation} \label{eq:cost_mixing_iteration} (\text{mixing time}) \cdot (\text{cost per iteration}). \end{equation}\]

The mixing time is expected to vary significantly for different systems of interest (based on classical Monte Carlo intuition), but for systems appearing in nature, one may be optimistic based on the observed fast thermalization of physical systems. The cost per iteration is dominated by the quantum phase estimation subroutine, which then scales with a certain energy resolution. An overall gate complexity can be roughly, e.g., \(\text{poly}(n, \beta,1/\epsilon)\). However, as new algorithms are still being proposed, we do not give more concrete estimates of the complexity. Indeed, to put together an end-to-end resource estimate, one needs to design better algorithms to reduce the cost per iteration as well as to estimate the mixing time (e.g., by exact diagonalization of the map for small system sizes). Of course, if Gibbs sampling is employed as a heuristic (as in many classical applications of Monte Carlo methods), the cost will be empirical.

Caveats

On the one hand, the superpolynomial \(\mathcal{O}(\sqrt{2^n})\) complexity for Gibbs sampling that appears explicitly in Eq. \(\eqref{eq:Gibbs_complexity}\) is necessary in general (for sufficiently large \(\beta\) it allows one to solve NP-hard or even QMA-hard problems in the general case). On the other hand, most physical Hamiltonians (if they appear to thermalize in nature) should be simulable without exponential hidden prefactors. The Monte Carlo–style approach to Gibbs sampling attempts to mimic nature more closely than the other algorithms with guaranteed complexities mentioned above; hence, it looks more promising for obtaining polynomial runtimes, but this must be verified through system-specific analysis or hardware demonstrations.

Finally, if the Hamiltonian comes from classical problems (such as solving semidefinite programs), loading the instance may have exponential cost (\(\mathrm{e}^{\Omega(n)}\)), which in the above presentation is hidden in the assumption of a block-encoding of classical data. Additionally, it is unclear whether systems arising from classical data, rather than underlying physical models, should be expected to "thermalize" quickly (i.e. whether Monte Carlo–style algorithms converge in a small number of iterations).

Example use cases

  • Multiplicative Weights Update (MWU) method and conic programming: Gibbs sampling is the main source of quantum speedup in the MWU method, which is used to solve semidefinite programs and other conic programs [16, 17, 11, 12, 18]. Existing analyses in this direction have employed Gibbs samplers with a guaranteed quadratic (but no larger) speedup, rather than the more heuristic and recent Monte Carlo–style algorithms.
  • Quantum chemistry: An important step of estimating the ground state energy of electronic structure Hamiltonians is generating an ansatz state that has a large overlap with the ground state. This might be done via Gibbs sampling at sufficiently low temperatures; the overlap with the ground state is \(e^{-\beta E_0}/\mathcal{Z}\).
  • Condensed matter physics: Similar to quantum chemistry, Gibbs sampling provides a method for producing ansatz states for ground state energy calculation. Furthermore, condensed matter physicists are often interested in material properties at finite temperatures so that the Gibbs state can be equally interesting as the ground state itself.
  • Computing partition functions: One of the early references to develop quantum Gibbs samplers [7] applied it to the problem of estimating the partition function \(\mathcal{Z}\) up to small relative error. The partition function contains all the relevant thermodynamic information of the system.

Further reading

Gibbs sampling has been studied in several specific cases. For example, [19] studied Gibbs sampling of local Hamiltonians in 1D. Moreover, [20] studied commuting spatially-local Hamiltonians and showed conditions under which they thermalize in polynomial time, suggesting efficient Gibbs sampling via Monte Carlo–style methods. These conditions hold for any 1D system at any temperature, and in any higher spatial dimension above a certain threshold temperature.

Bibliography

  1. W Keith Hastings. Monte carlo sampling methods using markov chains and their applications. Biometrika, 57:97–109, 4 1970. doi:10.1093/biomet/57.1.97.

  2. K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete. Quantum metropolis sampling. Nature, 471(7336):87–90, 3 2011. arXiv: https://arxiv.org/abs/0911.3635. doi:10.1038/nature09770.

  3. Chi-Fang Chen and Fernando G. S. L. Brandão. Fast thermalization from the eigenstate thermalization hypothesis. arXiv: https://arxiv.org/abs/2112.07646, 2021.

  4. Oles Shtanko and Ramis Movassagh. Algorithms for gibbs state preparation on noiseless and noisy random quantum circuits. arXiv: https://arxiv.org/abs/2112.14688, 2021.

  5. Patrick Rall, Chunhao Wang, and Pawel Wocjan. Thermal state preparation via rounding promises. arXiv: https://arxiv.org/abs/2210.01670, 2022.

  6. Chi-Fang Chen, Michael J. Kastoryano, Fernando G. S. L. Brandão, and András Gilyén. Quantum thermal state preparation. arXiv: https://arxiv.org/abs/2303.18224, 2023.

  7. David Poulin and Pawel Wocjan. Sampling from the thermal quantum gibbs state and evaluating partition functions with a quantum computer. Physical Review Letters, 103(22):220502, 2009. arXiv: https://arxiv.org/abs/0905.2199. doi:10.1103/PhysRevLett.103.220502.

  8. Chen-Fu Chiang and Pawel Wocjan. Quantum algorithm for preparing thermal gibbs states-detailed analysis. In Quantum Cryptography and Computing, volume 26, 138–147. 2010. arXiv: https://arxiv.org/abs/1001.1130. doi:10.3233/978-1-60750-547-1-138.

  9. Maris Ozols, Martin Roetteler, and Jérémie Roland. Quantum rejection sampling. ACM Trans. Comput. Theory, 8 2013. arXiv: https://arxiv.org/abs/1103.2774. URL: https://doi.org/10.1145/2493252.2493256, doi:10.1145/2493252.2493256.

  10. Anirban Narayan Chowdhury and Rolando D. Somma. Quantum algorithms for gibbs sampling and hitting-time estimation. Quantum Information and Computation, 17(1&2):41–64, 2017. arXiv: https://arxiv.org/abs/1603.02940. doi:10.26421/QIC17.1-2.

  11. Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf. Quantum sdp-solvers: better upper and lower bounds. Quantum, 4:230, 2020. Earlier version in FOCS'17. arXiv: https://arxiv.org/abs/1705.01843. doi:10.22331/q-2020-02-14-230.

  12. Joran van Apeldoorn and András Gilyén. Improvements in quantum sdp-solving with applications. In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP), 99:1–99:15. 2019. arXiv: https://arxiv.org/abs/1804.05058. doi:10.4230/LIPIcs.ICALP.2019.99.

  13. András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st ACM Symposium on the Theory of Computing (STOC), 193–204. 2019. arXiv: https://arxiv.org/abs/1806.01838. doi:10.1145/3313276.3316366.

  14. Theodore J. Yoder, Guang Hao Low, and Isaac L. Chuang. Fixed-point quantum search with an optimal number of queries. Physical Review Letters, 113(21):210501, 2014. arXiv: https://arxiv.org/abs/1409.3305. doi:10.1103/PhysRevLett.113.210501.

  15. Pawel Wocjan and Kristan Temme. Szegedy walk unitaries for quantum maps. Communications in Mathematical Physics, 2023. arXiv: https://arxiv.org/abs/2107.07365. URL: https://doi.org/10.1007/s00220-023-04797-4, doi:10.1007/s00220-023-04797-4.

  16. Fernando G. S. L. Brandão and Krysta M. Svore. Quantum speed-ups for solving semidefinite programs. In Proceedings of the 58th IEEE Symposium on Foundations of Computer Science (FOCS), 415–426. 2017. arXiv: https://arxiv.org/abs/1609.05537. URL: http://ieee-focs.org/FOCS-2017-Papers/3464a415.pdf, doi:10.1109/FOCS.2017.45.

  17. Fernando G. S. L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M. Svore, and Xiaodi Wu. Quantum sdp solvers: large speed-ups, optimality, and applications to quantum learning. In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP), 27:1–27:14. 2019. arXiv: https://arxiv.org/abs/1710.02581. doi:10.4230/LIPIcs.ICALP.2019.27.

  18. Joran van Apeldoorn and András Gilyén. Quantum algorithms for zero-sum games. arXiv: https://arxiv.org/abs/1904.03180, 2019.

  19. Ersen Bilgin and Sergio Boixo. Preparing thermal states of quantum systems by dimension reduction. Physical Review Letters, 105:170405, 10 2010. arXiv: https://arxiv.org/abs/1008.4162. URL: https://link.aps.org/doi/10.1103/PhysRevLett.105.170405, doi:10.1103/PhysRevLett.105.170405.

  20. Michael J. Kastoryano and Fernando G. S. L. Brandão. Quantum gibbs samplers: the commuting case. Communications in Mathematical Physics, 344(3):915–957, 2016. arXiv: https://arxiv.org/abs/1409.3435. URL: https://doi.org/10.1007/s00220-016-2641-8, doi:10.1007/s00220-016-2641-8.