Connection between single-layer Quantum Approximate Optimization Algorithm interferometry and thermal distributions sampling

The Quantum Approximate Optimization Algorithm (QAOA) is an algorithm originally proposed to find approximate solutions to Combinatorial Optimization problems on quantum computers. However, the algorithm has also attracted interest for sampling purposes since it was theoretically demonstrated under reasonable complexity assumptions that one layer of the algorithm already engineers a probability distribution beyond what can be simulated by classical computers. In this regard, a recent study has shown as well that, in universal Ising models, this global probability distribution resembles pure but thermal-like distributions at a temperature that depends on internal correlations of the spin model. In this work, through an interferometric interpretation of the algorithm, we extend the theoretical derivation of the amplitudes of the eigenstates, and the Boltzmann distributions generated by single-layer QAOA. We also review the implications that this behavior has from both a practical and fundamental perspective.


INTRODUCTION
Variational quantum algorithms (VQAs) are a promising framework for solving computationally hard tasks (Cerezo et al., 2021;Bharti et al., 2022).These algorithms are based on a quantum circuit with tunable parameters acting as ansatz.Finding the minimum energy state among a discrete set of possible solutions (combinatorial optimization problems) (Nemhauser and Wolsey, 1988;Moll et al., 2018), or drawing samples from a classical probability distribution (sampling problems) (Lund et al., 2017;Wild et al., 2021) are two examples of worthwhile challenges addressed in this framework.
Among the plethora of proposed VQAs, the Quantum Approximate Optimization Algorithm (QAOA) (Farhi et al., 2014) has received special attention from the scientific community, with remarkable empirical and theoretical results on the algorithm's performance (Harrigan et al., 2021;Farhi et al., 2022;Blekos et al., 2023).The ansatz of QAOA is inspired by a trotterized adiabatic evolution capable of approximating the minimum energy state of a given cost function.An advantage of the practical implementation of QAOA is that the number of variational parameters in the quantum circuit to be optimized is independent of the problem size, on the other hand, the depth of the circuit increases with the quality of the approximation of the fundamental state.Beyond QAOA's original use in approximating ground states, QAOA has been recently shown to have interest from a global sampling perspective (Farhi and Harrow, 2019).It has been shown that the shallowest version of the algorithm produces a probability distribution that cannot be simulated by classical computation since, otherwise, the Polynomial Hierarchy would collapse.Therefore, sampling the QAOA circuit in some variational parameters range could exhibit quantum supremacy (Preskill, 2018), understood as a task that can be more efficiently done in a quantum computer than in a classical one.This result is connected to previous demonstrations of the hardness of simulating quantum circuits, for instance Boson Sampling (Brod, 2015) or Instantaneous Quantum Polynomial time (IQP) circuits (Bremner et al., 2010(Bremner et al., , 2016)).
Despite the importance of this theoretical result, it gives no clue about the nature of the probability distribution generated by QAOA, nor whether sampling from it has any practical interest beyond combinatorial optimization or a quantum supremacy demonstration.Following this direction, in (Díez-Valle et al., 2023), we studied the global structure and probability distribution in the energy space of the singlelayer QAOA ansatz.In that letter, we provide evidence for a clear connection between the QAOA and thermal distribution sampling with probability amplitude distributions resembling Boltzmann distributions at low temperatures.We call such distributions QAOA pseudo-Boltzmann states.The excellent performance of QAOA revealed at temperatures below the state-of-the-art theoretical limit to fast mixing of Markov Chain Monte Carlo methods (Eldan et al., 2021).
In the current manuscript, we expand the theoretical derivation of pure thermal-like QAOA states introduced in (Díez-Valle et al., 2023).The paper is structured as follows.First, in Section 2, we review the main ingredients of the Quantum Approximate Optimization Algorithm and introduce an intuitive picture of QAOA as an interferometer in energy space.In Section 3.1 we extend the interferometric picture to the multilayer scenario and derive the analytical amplitudes of the final wavefunction.In Section 3.2, we specify the results of the previous section to the case where an Ising Model Hamiltonian defines the cost function.We detail the internal correlations between the eigenstates of the Hamiltonian in degenerate and nondegenerate scenarios, leading to a closed expression for the probability amplitudes of the QAOA wavefunction.We connect such probability amplitudes to the sampling of Boltzmann or thermal-like distributions in Section 4, also analyzing the features that these results reveal about the distribution behavior as we change the variational parameters of the algorithm.We conclude with a perspective on the importance of the pseudo-Boltzmann QAOA states and approximate thermal sampling in Section 5.

The Quantum Approximate Optimization Algorithm
The Quantum Approximate Optimization Algorithm (QAOA) aims to find good approximate solutions to combinatorial optimization problems (Farhi et al., 2014), which are defined by a classical objective function E(x) mapping N -binary strings to real values: The algorithm's goal is to find a binary string x * that minimizes the function E min ≡ min x E(x), or at least that achieves a good approximation ratio r ∈ [0, 1] This optimization problem is equivalent to finding the ground state of a spin Hamiltonian Ê, where each binary variable x i describes the state of one spin A well-known approach to solving these problems is the Quantum Adiabatic Algorithm (QAA), which ensures the achievement of the global minimum given a sufficiently long run time T (Farhi et al., 2000(Farhi et al., , 2001)).The adiabatic protocol is guided by a time-dependent Hamiltonian such as where Ĥx = N i σ x i with σ x i the Pauli X operator acting on the ith qubit, and f (x) is the schedule function with f (0) = 0 and f (1) = 1.The adiabatic theorem guarantees that a sufficiently slow time evolution will keep the system in its ground state.Thus, an evolution that starts with the ground state of Ĥx , will approximately bring the system to the ground state of Ê, To ensure success, the runtime T must scale as T = O ∆ −2 min , where ∆ min is the minimum spectral gap (Albash and Lidar, 2018).QAOA proposes a trotterized approximation to adiabatic evolution consisting of a quantum circuit built by the alternation of the following two operators with H denoting single-qubit Hadamard gates, and where γ = (γ 1 , γ 2 , ..., γ p ) and θ = (θ 1 , θ 2 , ..., θ p ) are set of variational angles that can be tuned to approximate the ground state of the cost Hamiltonian Ê, and p is the number of layers that determines the depth of the quantum circuit and the accuracy of the algorithm.The infinite depth limit p → ∞ returns the adiabatic evolution showing that we can get a good enough approximation to the minimum of the optimization problem for sufficiently large p.

QAOA as an interferometer in energy space
Let us highlight QAOA's potential in small-depth regimes by introducing an interpretation QAOA's circuit as an interferometer operating globally in energy space.We will derive analytical results for the interference amplitude in later sections, focusing on the p = 1 single-layer ansatz (see Eq. ( 7)), As sketched in Fig. 1, the Hadamard gate splits the quantum state |0⟩ into a superposition of all states in the computational basis, acting like a multidimensional "mirror" Then the evolution with the diagonal cost Hamiltonian Ê |x⟩ = E x |x⟩ , given by the operator U E ( Ê, γ), imparts phases on all states, acting like the branches of the interferometer Finally, the mixing operator R x (θ) ⊗N recombines the energy states so that the interference transforms the energy-dependent relative phases in measurable probability amplitudes This intuition is much clearer when looking at the algorithm operating on a single qubit, with rotation R x (θ), Consider the following two-state subspace where x i = {0, 1} (so x 0 and x 1 differ by only one bit), with energies The probability amplitudes before the local interference are respectively.After applying the local gate in Eq. ( 14) on the qubit which corresponds to the different bit, the interference shifts the amplitudes to N x 0 and N x 1 , where Thus, the probability amplitudes become Note that the relative phase between the states γF x 1 ,x 0 = γ (E x 1 − E x 0 ) controls the amplification.Therefore, assuming sin (2θ) < 0 and γ > 0 or vice versa and provided the angle γ is sufficiently small, the interference term always increases the population of the lowest energy state and reduces that of the other.Furthermore, Eq. ( 17) also reflects how the relative sign between sin (2θ) and γ controls the direction of optimization.The probability amplitude of the lowest/highest energy state is enhanced when the signs are opposite/equal.
This interferometric behavior expanded to the N -level scenario imposes nontrivial, structure-dependent upper bounds on the |θ| and |γ| angles to avoid symmetries and random scrambling of the energy states.In particular, Eq. ( 17) shows that for the single-qubit interference with n ∈ Z.The N -qubit interference makes this picture richer and more complicated as derived in the next sections.Analytical and numerical results for the single layer QAOA on general Ising spin models show that the optimal |γ| for an N qubit system actually scales as O(N (N |e|) −1/2 ) (Ozaeta et al., 2022;Díez-Valle et al., 2023), where |e| is the number of edges or non-null elements in the coupling matrix J (see Eq. ( 27)).For instance, in the two-level system Ê = 1 2 ∆σ z the probability of measuring the ground state |σ z = −1⟩ and the highest energy state |σ z = +1⟩ is Figure 2. Quantum Approximate Optimization Algorithm interferometry on the two-level spin system defined by the Hamitonian Ê = 1 2 ∆σ z with ∆ = 1.We plot the probability amplitude of the ground state σ z = −1 (dashed lines), and the highest energy state σ z = +1 (solid lines) for different values of the QAOA angles θ and γ.We highlight the points of maximum amplification θ = ±π/4, γ = π/2 (dotted brown lines).

General scenario
Let us extend the two-level picture, to a general framework dealing with the complete interference spectrum of a single-layer QAOA circuit.We focus on how the Quantum Approximate Optimization Algorithm transforms the wavefunction amplitudes As previously shown, a local R x rotation acting on the i-th qubit mixes the amplitude between each state x with the corresponding state x ′ that differs only on the value of the i-th bit (see Eq. ( 17)).Thus, the action of the whole mixing operator R ⊗N x is that of a complete mixing between all states in the computational basis The weight of each mixing process depends on the number of bits that must be flipped to go from one state x to the other x ′ -i.e., the Hamming distance between both states-, and also on the oscillating terms induced by the angles θ and γ.The amplitude of the state generated by the algorithm after the mixing operator U ( Ĥx , θ) is given by an interference formula, where H x,x ′ is the Hamming distance between two bit configurations x and x ′ that represent eigenstates of a spin model with energy E x and E x ′ .We can unify the weight sum terms into a single exponential by the following change of variables in which the rotation angle θ is reparameterized in terms of an exponent r and a normalization R, where due to the symmetries of the interference we consider θ ∈ (0, π 2 ) without loss of generality.The sign of γ allows us to define if the interference increases the population of lower or higher energy states so that we cover all the possibilities.By this transformation, the interference amplitudes in Eq. ( 22) become a sum of exponentials over the entire configuration space, Since Eq. ( 24) encompasses all eigenstates of the system encoded in its energy E x ′ and Hamming distance H x,x ′ with the reference state x, it is convenient to introduce a probability distribution relating distances in the computational base to the energy spectrum, This expression represents the relative number of eigenstates that, given a reference state x, have a Hamming distance H to that state and energy equal to E. In other words, this distribution captures the internal correlations inherent to the specific spin model.With this definition at hand, we can express the previous sum in Eq. ( 24) as the average over this probability distribution Hence, the interference amplitude in Eq. ( 20) depends on the structure of the energy levels of the spin system manifested in the probability distribution p(H, E; x).In the following, we study such a structure to derive a unified expression of the QAOA probability amplitudes for universal Ising spin models.

Ising Models
Next, we focus on a universal set of nondeterministic polynomial time hard (NP-Hard) Ising Models (Barahona, 1982) and derive this scenario's interference probability amplitude formula.

Frontiers
A broad spectrum of combinatorial optimization problems can be represented as a spin model described by an Ising Hamiltonian (Lucas, 2014), where s = {−1, +1} N , N is the number of variables, J is an N -by-N square coupling matrix, and the magnetic field h is a vector of N coefficients.The quantum version of this Hamiltonian ÊI (σ z ) is simply obtained by replacing the spin variables s with the corresponding Pauli-Z matrices σ z .The coupling matrix defines a structure in the problem that can be represented by a weighted graph with N vertices, connected by undirected edges i ↔ j that have associated weights J ij , J ji .In this work, we study families of models in which the J ij coefficients are randomly drawn from a normal distribution N (µ = 0, σ 2 ).
This structure together with the magnetic field values defines families of optimization problems with different inner correlations and degeneracies.In the following subsections, we derive the interference amplitude for two scenarios that involve distinct energy level structure due to intrinsic global symmetries in the model.These differences necessitate a separate study of both cases.Nevertheless, as explained below the slightly different behaviors of these models converge to a common expression of the QAOA probability amplitude distribution |F (x)| 2 in Eq. ( 20).The two scenarios are represented by two well-known combinatorial optimization problems: Quadratic Unconstrained Binary Optimization (QUBO) (Kochenberger et al., 2004(Kochenberger et al., , 2014)), with a non-degenerate energy spectrum, and the maximum cut (MaxCut) (Nannicini, 2019;Sung et al., 2020;Harrigan et al., 2021), which exhibits a global Z 2 symmetry.

Non-degenerate Ising Models
The family of QUBO problems is composed of NP-Hard binary combinatorial optimization problems with the following associated cost function: where x ∈ {0, 1} N , and Q is an N -by-N square matrix.The mapping of Eq. ( 28) to the Ising Hamiltonian in Eq. ( 27), x = 1 2 (1 + s), leads to where We consider Q matrices where the non-zero coefficients are randomly drawn from a standard normal distribution N (µ = 0, σ 2 = 1).Note that the optimization of the function in Eq. ( 29) is then equivalent to that of Eq. ( 28).
The cost function in Eq. ( 29) exhibits no global symmetries.In such non-degenerate situations, the eigenstates appear to be ordered from low to excited states, developing a unique probability distribution p(H, E; x) that is centered at the center of the spectrum and shows a remarkable correlation between H xx ′ and E x ′ .Such behavior can be well observed in Fig. 3.Here we plot the distribution in the Hamming distance-energy plane of 25000 samples drawn from a continuous approximation to p(H, E; x) of a single instance of QUBO.In specific, we first estimate the continuous probability density function by a Kernel Density Estimation (KDE) on the actual discrete instance data (all H xx ′ and E x ′ pairs).Then we sample from such density function to obtain the plotted distributions.We also perform a fitting of the final distribution to a Gaussian mixture of one or two Gaussians using Variational Inference and plot the corresponding confidence ellipsoid.This fit was obtained with the BayesianGaussianMixture method of the Python package scikit-learn (Pedregosa et al., 2011).We show the result for x being the ground state (a) and the highest energy state (c).For the sake of clarity, let us consider the upper-left plot, which illustrates the probability distribution p(H, E; x) for the reference state x as the ground state.All alternative states x ′ possess higher energy values compared to E x .The difference E x ′ − E x shows a potential correlation with the number of spin flips required to transition from x to x ′ , i.e the Hamming distance H xx ′ .This positive correlation is evident through the presence of a Gaussian function oriented along a diagonal.On the contrary, opting for x as the highest energy state (bottom-left plot) reveals an opposite trend: the greater the number of spin flips, the lower the energy level.The Gaussian orientation depends on the covariance between H and E, which is highly correlated with the energy of the reference state E x , as shown in Fig. 4 of Díez-Valle et al. (2023).In Eq. ( 30) E[•] denotes the expected value of the variable.
Our simulations confirm that, in the non-degenerate scenario, the probability distribution p(H, E; x) resembles a bivariate Gaussian with a correlation between variables defined by x and its rank in the energy spectrum: where The spin model defines all the parameters of the bivariate Gaussian except the correlation parameter ρ, which encapsulates the whole dependence of the distribution on x.
To test the hypothesis that p(H, E; x) in Eq. ( 25) is well approximated by the continuous bivariate Gaussian distribution (31) we perform a graphical multivariate normality test on the (H xx ′ ,E x ′ ) data of 1000 independent instances of QUBO.As in Fig. 3, we analyze both the low and the high energy regimes with x the ground state and the highest energy state respectively.We use a technique based on the squared Mahalanobis distances between the sample points s = (H xx ′ , E x ′ ) and their averages µ where Σ = σ E σ EH σ EH σ H is the covariance matrix, and ρ = σ EH σ E σ H .The Mahalanobis distance is a multivariate measure to quantify the distance between a point and a distribution (Mahalanobis, 1936).Moreover, it is a useful tool to check when multidimensional data was sampled from a normal distribution, since the probability density p of a set of normally distributed samples s in any dimension is entirely determined by the Mahalanobis distance, Therefore, showing that the Mahalanobis distance D M (s), with s = (H xx ′ , E x ′ ), follows the chi-squared distribution with 2 degrees of freedom, would demonstrate that p(H, E) is compatible with a bivariate Gaussian distribution sampling (see Eq. ( 31)).Indeed, Fig. 4a and Fig. 4c shows the perfect agreement for at least 99.8% of the spectrum.M and the same theoretical quantiles of the χ 2 2 distribution.Over the scatter plot, we plot the straight line that would follow the points if D 2 M and χ 2 2 were described by the same distribution.We can see how about 99.8% of the D 2 M distribution is well described by the bivariate Gaussians, with the exception of ∼ 0.2% of outliers in the tail.
Hence, returning to the probability amplitude of single-layer QAOA on non-degenerate Ising models, Eq. ( 26) together with the probability distribution in Eq. ( 31) leads to the following interference amplitude,

Degenerate Ising Models
The MaxCut problems are a family of combinatorial optimization problems consisting in minimizing the following objective function: with x ∈ {0, 1} N , and Q is an N -by-N square matrix.Finding the minimum or an approximate solution very close to the minimum is known to be NP-Hard (Håstad, 2001).Again, we can map the binary cost function optimization in Eq. ( 37) to an equivalent spin Hamiltonian optimization (Eq.( 27)), where J = Q, and the magnetic field h cancels out.As in the non-degenerate case, the non-zero coefficients of Q are drawn from N (µ = 0, σ 2 = 1).This class of problems includes the Sherrington-Kirkpatrick model (Sherrington and Kirkpatrick, 1975) when all vertices of the graph are connected, i.e. when the coupling matrix has no null coefficients.
In contrast to QUBO, the cost function in Eq. ( 38) exhibits a global Z 2 symmetry that keeps invariant the energy under a global spin flip ẼMaxCut (s) = ẼMaxCut (−s), or what it is the same ẼMaxCut (x) = ẼMaxCut (1 − x).In the non-degenerate scenario (Eq.( 29)) the presence of the magnetic field h breaks the symmetry.The existence of such global symmetries in the Hamiltonian results in the division of the Hilbert space into two or more distinct hierarchies of eigenstates.Note that if x a = (x a 1 , x a 2 , . ..) represents a ground state, then there is another ground state in the opposite sector x b = 1 − x a .Consequently, a single excited state x ′ can be seen now as arising from flipping spins in either of the ground states x a or x b .Accordingly, an identical spin configuration exhibits significantly different Hamming distances to both ground states, while its energy remains unchanged.In specific, the Hamming distance reaches its maximum when considering the separation between the two degenerate ground states, where H x a x b = N .Then, this symmetry implies that given any two states x and x ′ , there is a unique alternative state x ′′ such that This phenomenon results in the separation of the probability distribution into a sum of two distributions, each measured with respect to one of the eigenstate hierarchies, where with A, B two complementary subspaces in the bit configurations space that represent the two eigenstate hierarchies.Since each distribution p ± (H x,x ′ , E x ′ ; x) defines itself a non-degenerate Hilbert subspace, it is natural to expect that they individually behave similarly to the probability distribution of the non-degenerate scenario explained in the previous section.Indeed, as shown in Fig. 3, these distributions resemble two shifted bivariate Gaussian distributions, where µ E , µ H and σ H are the same as Eq. ( 32), h 0 > 0 is a constant shift, and we have two separate correlation factors ρ + (x) = −ρ − (x) ≡ ρ.As in QUBO, we prove the Gaussian distribution of the eigenstates by a multivariate normality test based on the Mahalanobis distance.In this case, we first need to group the states in two clusters that represent A and B hierarchies.As previously explained, we do this by the fit of the sample points s = (H xx ′ , E x ′ ) to a mixture of two Gaussians so that we identify which states should belong to p + or p − .Then we calculate the Mahalanobis distance of all samples using their corresponding mean and covariance estimated from the Gaussian mixture fit.Again, we find a good correspondence between the Mahalanobis distance and the χ 2 2 distribution for the vast majority of the energy spectrum, concluding that p ± (H, E; x) are well approximated by the Gaussian expresions in Eqs. ( 42) and (44).
As in the non-degenerate case, the correlation factors ρ ± encapsulate the entire dependency in x, but we now have two functions that influence one another.Eq. ( 22) along this p + , p − interference leads to Therefore, we observe that the interference translates into a mixture of two exponentials together with an oscillatory term, with β ′ ≡ −γπσ E σ H . Nevertheless, except in the regime when r = − log(tan θ) ≈ 0 (θ ≈ π/4), one exponential is clearly dominant over the other due to the presence of the h 0 -shift, translates into an evident σ EH − E x correlation that can be numerically observed in Fig. 4 of (Díez-Valle et al., 2023).This dependence can be expressed as the sum of a linear function and a stochastic value ω with zero mean where c ∈ R > 0. In spite of the presence of the random term ω, the trend in Eq. ( 55) is noticeable.Therefore, introducing Eq. ( 55) in Eq. ( 54) leads to a thermal-like probability amplitude distribution for the Ising models, where β = cπγ for non-degenerate problems, and β = sgn(r)cπγ for degenerate ones, plays the role of the inverse temperature.Then the probability distribution in energy space can be expressed as, where d(E) is the density of states that in Ising models resembles a normal distribution centered in intermediate energies.
The Boltzmann distribution in Eq. ( 56) together with the random fluctuations ω manifests in regimes of optimized parameters, but the distribution may be modified by manipulating the angles θ ∈ (0, π/2) and γ (see Figs. 5 and 6).By the analysis of Eq. ( 56) and previous expressions, additional conclusions can be drawn about the shape of the distribution.
• The angle γ controls the direction of the optimization.For specific γ regimes, the single-layer QAOA probability amplitude distribution on Ising models resembles a thermal distribution at temperature T = β −1 = 1 cπ|γ| such that the ground state of the system aligns with the peak amplitude.Switching the sign of this angle is the same as changing the sign of the energy, and therefore of the temperature T = −1 cπ|γ| .These negative temperature states increase the probability of the highest excited state.• The degenerate Ising models exhibit antisymmetric behavior at the angle θ since T = sgn(r) cπ|γ| , with r = − log(tan θ).Therefore, for θ < π/4 and θ > π/4 the pure QAOA state resembles a thermallike distribution with positive/negative temperature.Eq. ( 45) shows that when θ ≈ π/4 the thermal type distribution disappears and we find a mixture of two Boltzmann exponentials with opposite temperatures.This behavior can be well observed in Fig. 6b.For non-degenerate Hamiltonians, in the r ≪ 0 regime the r 2 σ 2 H and −2rµ H terms of Eq. ( 52) increase the noise and the Boltzmann distribution vanishes.
• The Boltzmann distribution is apparent for a finite interval of angles γ ∈ (0, γ c ].The lowest temperature T is reached near γ ≈ γ opt , where γ opt is the angle that minimizes the mean energy in the variational principle.When |γ| > γ c , the Boltzmann term γπρσ E σ H becomes marginal in comparison with −γ 2 σ 2 E in Eq. ( 52), so that |F (x)| 2 ≈ cte.• For γ < γ c , we observe that β grows linearly with the angle γ, consistent with the numerical results of (Díez-Valle et al., 2023).

OUTLOOK
Sampling from complex probability distributions is a valuable computational task, both for its difficulty and its broad applicability.Quantum states projected on the computational basis are in essence classical probability distributions and measuring them is the same as sampling from such distributions.Therefore, a quantum computer can be roughly seen as a machine capable of creating exotic or suitable probability distributions thanks to quantum phenomena.Indeed, the first claim of quantum supremacy (Preskill, 2018) was performed on a sampling problem (Arute et al., 2019), and this task has been also envisioned as one of the candidates to show a practical quantum advantage in the near term (Wu et al., 2021;Zhong et al., 2021;Layden et al., 2023).
In this context, the Boltzmann distributions of Ising models could be good candidates to show such quantum advantage for two reasons.First, these distribtuions are classically intractable at low-temperature regimes.The most popular strategy to simulate thermal distributions is the use of Markov Chain Monte Carlo algorithms (MCMC) which guarantee convergence to the Boltzmann distribution at a given temperature T .While these methods work very well for relatively high temperatures, the number of iterations needed to converge scales exponentially when the temperature is too low.Much effort has been done to determine the range of temperatures that ensure rapid mixing and therefore polynomial convergence of MCMC.To the best of our knowledge, the state-of-the-art theoretical bound indicates that MCMC always converges in polynomial time to the thermal distribution of an Ising model at a temperature higher than ∥J∥ (Eldan et al., 2021), although practical realizations and state-of-the-art methods (e.g.parallel tempering, population annealing ...) might overcome this threshold.(Díez-Valle et al., 2023) showed that the single-layer QAOA already approximates Boltzmann distributions at temperatures beyond this theoretical MCMC bound.Second, there is a wide range of fields that would be impacted by an improvement in Boltzmann distribution sampling.In statistical mechanics, sampling from this distribution is crucial for simulating physical systems at thermal equilibrium and for computing observables such as magnetization in Ising models.Furthermore, machine learning makes use of this distribution in unsupervised learning techniques known as Boltzmann Machines (Ackley et al., 1985).Combinatorial optimization is another field of interest since some algorithms employ Boltzmann distribution sampling at decreasing temperatures as a subroutine to find the minimum of a cost function (Kirkpatrick et al., 1983).
The thermal-like distributions of the single-layer Quantum Approximate Optimization algorithm reveal a nice connection between quantum algorithmics and statistical physics that can help to gain a better understanding of the behavior of these ansätze.But in a broader perspective, the QAOA sampling from approximate Boltzmann distributions, also known as Gibbs sampling, may be extended to the multi-layer scenario with an improvement in the achievable temperature (Lotshaw et al., 2022), or to QAOA mixedstate ansätze to train unsupervised learning models as implemented in (Verdon et al., 2019).Furthermore, since QAOA resembles a trotterized approximation to an adiabatic quantum evolution, this picture might be expanded to the infinite depth scenario p → ∞ and the engineering of general time-dependent adiabatic passages.
Additionally, the approximate Boltzmann distributions of single-layer QAOA states present collateral implications in quantum information theory.For example, this behavior makes them suitable as warm initial states for more complex ansätze as has been demonstrated in (Leontica and Amaro, 2023).Another recent work (Sud et al., 2022) shows how a tight dependence between the energy distribution of the spin model and the final probability amplitude of the QAOA states allows a more efficient classical heuristic optimization of the QAOA parameters.The Boltzmann distribution unambiguously connects the energy of the eigenstates with their amplitudes, thus providing further arguments and explanations on these heuristics.
We are confident that this study, in combination with earlier work by (Díez-Valle et al., 2023), are just two initial works in a new field in which variational and other types of circuits are analyzed from a physical perspective, understanding not only their computational power but potentially significant physical insight behind the quantum computer's dynamics.

Figure 1 .
Figure 1.Interferometric interpretation of the Quantum Approximate Optimization Algorithm circuit.

Figure 3 .
Figure 3. Continuous approximation to the probability distribution p(H, E; x) estimated by 25000 samples drawn from a kernel density estimation (KDE) on a single instance of a non-degenerate problem (a and c), and of a degenerate one (b and d).We plot the extreme cases when the reference state x) is the ground state (a and b), and the highest energy state (c and d).We also fit the obtained distribution to a Gaussian mixture of one (a and c) or two (b and d) bivariate Gaussians, showing the clustering of the samples and the confidence ellipsoid over the distributions.

Figure 4 .
Figure4.Multivariate normality test to show that the structure of the energy levels of the Ising models studied resembles probability distributions p(H, E; x) that can be defined by one or two continuous bivariate Gaussians.We plot the results for 1000 14-qubit random instances of QUBO (a,c) and MaxCut (b,d) when the reference state x is the ground state (a,b) and the highest energy state (c,d).For each instance of the problem, we calculated the Mahalanobis distances, D M , of all (H xx ′ ,E x ′ ) samples, and display the results of the obtained distributions.If such samples are compatible with a set of samples drawn from a bivariate Gaussian, the obtained Mahalanobis distances must follow the χ 2 2 distribution.In the MaxCut case (b,d), a prior fit to a Gaussian Mixture was performed, so that the x ′ states were clustered into two Gaussians with corresponding means and covariance matrices.The smaller plots show the probability density function of D 2 M as a histogram along with the theoretical χ 2 2 distribution.The larger plots show quantile-quantile plots displaying 500 quantiles of D 2 M and the same theoretical quantiles of the χ 2 2 distribution.Over the scatter plot, we plot the straight line that would follow the points if D 2 M and χ 2 2 were described by the same distribution.We can see how about 99.8% of the D 2 M distribution is well described by the bivariate Gaussians, with the exception of ∼ 0.2% of outliers in the tail.