Gibbs Sampling of Continuous Potentials on a Quantum Computer

TL;DR

The paper tackles the challenge of Gibbs sampling for continuous, non-convex energy potentials by mapping the problem onto a torus and solving the toroidal Fokker–Planck equation with a quantum ODE solver to prepare a Gibbs-state quantum register. A key innovation is the semi-analyticity framework, which guarantees sub-exponential decay of Fourier coefficients and enables efficient Fourier-based interpolation and upsampling, yielding high-precision sampling with zeroeth-order oracle access to the energy function. The authors prove interpolation bounds for semi-analytic functions, derive complexity estimates that show quantum speedups under Morse-type structure and favorable temperatures, and extend the approach to non-periodic domains via Chebyshev methods. These results offer a pathway to efficient, high-precision sampling for periodic continuous potentials and inform PDE-based quantum machine learning methods, with implications for energy-based models and beyond.

Abstract

Gibbs sampling from continuous real-valued functions is a challenging problem of interest in machine learning. Here we leverage quantum Fourier transforms to build a quantum algorithm for this task when the function is periodic. We use the quantum algorithms for solving linear ordinary differential equations to solve the Fokker--Planck equation and prepare a quantum state encoding the Gibbs distribution. We show that the efficiency of interpolation and differentiation of these functions on a quantum computer depends on the rate of decay of the Fourier coefficients of the Fourier transform of the function. We view this property as a concentration of measure in the Fourier domain, and also provide functional analytic conditions for it. Our algorithm makes zeroeth order queries to a quantum oracle of the function. Despite suffering from an exponentially long mixing time, this algorithm allows for exponentially improved precision in sampling, and polynomial quantum speedups in mean estimation in the general case, and particularly under geometric conditions we identify for the critical points of the energy function.

Figures (4)

• Figure 1: Schematics of the circuits of quantum oracles. All registers receive float-point representations of real numbers. (a) The oracle for an energy function $E_\theta$. The first register receives the model parameters $\theta \in \mathbb R^m$, the second register receives a data sample $x \in \mathbb R^d$. And, the last register is used to evaluate the energy function. (b) A controlled variant of the same oracle, controlled on a single qubit represented by the top wire. (c) The oracle of figure (a) receiving an augmented sample in superposition as the state $\ket{\Psi_x}= \sum_{b \in \{0, 1\}^d}\ket{(-1)^{b_i} \arccos x_i}$.
• Figure 2: (a) and (b) show two families of functions considered respectively in \ref{['ex:cosine-function']} and \ref{['ex:inv-cos']}. The functions are normalized so that $\int_{x\in[0,1]} \mathrm dx\, \left(f(x)\right)^2 = 1$. That is, $f^2$ represents a distribution over one period. Note how in (a) the smoothness of the functions is controlled by the parameter $z$ and in (b) it is controlled by $(z-1)^{-1}$. (c) and (d) show the Fourier interpolation accuracy on the two respective families of functions considered in (a) and (b). We demonstrate the interpolation error given the state $\ket{f_N}$ for different $N$. Note that in both cases having $N$ larger than our upper bounds on $a$ results in a sampling error less than $0.1$. The sampling error is shown with respect to the smoothness parameters $z$ and $(z-1)^{-1}$, obtained by the application of the upsampling algorithm using $M=200$. Recall from \ref{['ex:cosine-function']} and \ref{['ex:inv-cos']} that we may think of $\max(1,\frac{z}{2})$ and $\max(8,\frac{8}{z-1})$ as upper bounds on the (average) inverse convergence radius of the respective functions in panels (a) and (b).
• Figure 3: Applying the Fourier interpolation of \ref{['thm:interpolation']} on the input function $u(x) = e^{\cos 2\pi x}$ of \ref{['ex:cosine-function']}. The plot shows the interpolation results with $N=3$ and $M=10$. Filled circles correspond to the initial samples, and the hollow circles represent the interpolation output. The solid blue line represents the graph of the underlying function $u$. And the dashed green lines show the Fourier derivative estimations.
• Figure 4: The circuit for the oracle of discrete generator $\mathbb L$ comprising $2d(2N+1)$ copies of the energy potential oracle, $O_E$. To query $\mathbb L[x_1,x_2] = \bra{x_2}\mathbb L \ket{x_1}$, first the controlled-$U$ gate checks for the difference between $x_1$ and $x_2$: the third register is set to $\ket{i}$ if $x_1$ and $x_2$ defer only on their $i$-th entry. The state remains unchanged, if $x_1=x_2$, and it is set to a null state $\ket \perp$ otherwise. Conditioned on this third register being at state $\ket{i}$, another register (the fourth register) computes the distance between $x_1$ and $x_2$ along the $i$-th axis on the lattice (the controlled-$D$ gate). Again, conditioned on the state of the third register, we query the energy function at specific lattice points to compute either $\partial_iE(x_1)$ (if the third register is in $\ket i$), or $\nabla^2 E(x_1)$ (if the third register is in $\ket{0}$) using the sequence of controlled-$V_j$ gates. The estimation of these derivatives exploits Fourier spectral method (see \ref{['sec:app']}) and is applied via a circuit performing simple arithmetic.

Theorems & Definitions (107)

• Definition 3.1
• Theorem 3.2: Main interpolation result
• Theorem 3.3
• Theorem 4.1: Main sampling result
• Corollary 4.2
• Corollary 4.3: Mean estimation
• Definition 1.1: \ref{['def:body-semi-analytic']} in the manuscript
• Example 1.1
• Example 1.2
• Proposition 1.2