Generalized Probabilistic Approximate Optimization Algorithm

Domain: physics
Published: 2025-12-08
Editorial outcome: minor_revision

Abstract
We introduce the generalized Probabilistic Approximate Optimization Algorithm (PAOA), a classical variational Monte Carlo framework that extends and formalizes the recently introduced PAOA, enabling parameterized and fast sampling on present-day Ising machines and probabilistic computers. PAOA operates by iteratively modifying the couplings of a network of binary stochastic units, guided by cost evaluations from independent samples. We establish a direct correspondence between derivative-free updates and the gradient of the full Markov flow over the exponentially large state space, showing that PAOA admits a principled variational formulation. Simulated annealing emerges as a limiting case under constrained parameterizations, and we implement this regime on an FPGA-based probabilistic computer with on-chip annealing to solve large 3D spin-glass problems. Benchmarking PAOA against QAOA on the canonical 26-spin Sherrington–Kirkpatrick model with matched parameters reveals superior performance for PAOA. We show that PAOA naturally extends simulated annealing by optimizing multiple temperature profiles, leading to improved performance over SA on heavy-tailed problems such as SK–Lévy.

Introduction
Monte Carlo algorithms remain a central tool for exploring complex energy landscapes, especially in the context of combinatorial optimization and statistical physics. Classical methods such as simulated annealing (SA) have been widely applied across these domains, but their reliance on slowly equilibrating processes limits their performance on rugged energy landscapes1–6. New approaches are needed to construct non-equilibrium strategies that retain algorithmic simplicity while improving solution quality.

Inspired by the quantum approximate optimization algorithm (QAOA)7–17, Weitz et al.18 proposed a classical variational protocol based on the direct parameterization of low-dimensional Markov transition matrices. Their work introduced the term Probabilistic Approximate Optimization Algorithm (PAOA), raising the possibility of classical, variational analogs to QAOA within probabilistic architectures.

Building on this foundational work, we formalize and generalize PAOA. We move beyond the original proposal’s edge-local matrices to derive a global 2N×2N Markov-flow formulation applicable to any k-local Ising Hamiltonian of size N. This framework unifies a wide spectrum of variational ansätze, from global and local schedules to fully-parameterized couplings, and allows them to be stacked to an arbitrary depth p. Crucially, we connect this variational theory to practice by building on the p-computing framework19, where networks of binary stochastic units (p-bits) sample from Boltzmann-like distributions through asynchronous dynamics. This model has a demonstrated record of success in optimization and inference tasks20–23 and provides a natural substrate for our work. By implementing these dynamics on an FPGA, we demonstrate for the first time a scalable, hardware-based path for executing the PAOA.

We show that PAOA admits a broad class of parameterizations, including global, local, and edge-specific annealing schedules. Within this framework, SA emerges as a limiting case under constrained schedules. We implement this regime on an FPGA-based p-computer to demonstrate large-scale, high-throughput sampling. In contrast to standard SA, PAOA’s flexible parameterization enables the discovery of non-equilibrium heuristics that exploit structural features of the problem. One such example is demonstrated in heavy-tailed Sherrington–Kirkpatrick (SK) models, where PAOA assigns higher effective temperatures to strongly coupled nodes and lower temperatures to weakly connected ones, enabling annealing profiles with better performance over vanilla SA. This adaptive scheduling capability provides a powerful framework for the automated discovery of novel annealing heuristics.

We also benchmark PAOA against QAOA on the SK model using matched parameter counts and observe superior performance. This may not be surprising, as the practical performance advantages of QAOA over classical and probabilistic alternatives remain uncertain. A clear quantum advantage is expected to emerge in problems where interference plays an explicit role, such as in the synthetic examples constructed by Montanaro et al.12. Our work builds on the promising results of Weitz et al.18 on the max-cut problem by expanding PAOA into a scalable, hardware-compatible, and general-purpose optimization algorithm.

The remainder of this paper is organized as follows: we first present the theoretical formulation of PAOA, including its parameterization, Markov structure, and sampling-based approximation with p-bits. In order to illustrate the correspondence between analytical gradients and derivative-free optimization, we study a simple majority gate problem. We then demonstrate the recovery of SA as a special case and implement this regime on hardware, followed by a description of the FPGA architecture used to accelerate sampling. Using matched parameter counts, we benchmark PAOA against QAOA on the SK model. Finally, we study a heavy-tailed SK variant to show how PAOA can discover multiple annealing schedules that outperform standard single-schedule SA.

Results
Probabilistic Approximate Optimization Algorithm overview
The power of the PAOA comes from a conceptual departure from traditional optimization methods, such as SA. While SA explores a single, fixed energy landscape defined by a model Hamiltonian, PAOA treats the landscape itself as a variational object to be optimized. The goal is to dynamically reshape this landscape to make the ground state of the original problem easier to find.

This creates a two-level optimization loop, as illustrated in Fig. 1. In an outer loop, a classical optimizer adjusts a set of variational parameters, θ. In an inner loop, a probabilistic computer samples states from the variational landscape currently parameterized by θ. In the most general case, θ specifies the full set of couplings and fields {J(k), h(k)} at each layer k, so that the landscape itself is directly reconfigured by the optimizer. Simpler ansätze, such as global or node-wise schedules introduced in the Spectrum of variational ansätze subsection, can be viewed as constrained subsets of this general J, h parameterization (for example, a global inverse temperature amounts to scaling all entries of J by a single parameter). The sampler acts on M spins in the ansatz, while the original problem has N visible spins. In general M≥N because hidden spins may be introduced to enrich the representation, in direct analogy to neural quantum states24. Unless otherwise stated, we set M=N in the remainder of the paper, for simplicity, and leave the possibility of introducing hidden spins for better representation to future work.

During training, samples generated under θ are evaluated with a task-dependent cost function L(θ). In problems with a known target set of states (e.g., majority logic gate, discussed in the Representative example: majority gate subsection), L is typically chosen as a negative log-likelihood. In optimization problems such as 3D spin glasses, L is instead taken as the average energy over sampled configurations. In either case, feedback from the samples updates θ, decoupling the problem energy from the sampling dynamics and enabling PAOA to discover non-equilibrium pathways to solutions.

Spectrum of variational ansätze
The power and flexibility of PAOA lie in the choice of the variational parameters θ that define the search landscape. We refer to this parameterization as the variational ansatz. Let N denote the total number of spins in the problem. To explore the tradeoffs between expressiveness and generalization, in this paper, we study four representative ansätze, defined by the number of free parameters, Γ:

Γ = 1, 2, N, N(N-1)/2

When Γ=1, the landscape is scaled by a single, learned inverse temperature β(k) at each layer k. This is the global annealing schedule, which is the closest analog to SA, with the crucial difference that the schedule is learned rather than predefined.
When Γ=2, two independent schedules are assigned to distinct parts of the graph, typically divided based on node properties like weighted degree.
When Γ=N, each node is assigned a local schedule βi(k), allowing the algorithm to anneal different parts of the problem at different rates.
When Γ=N(N-1)/2, the entire symmetric coupling matrix J(k)ij is used as the variational ansatz, giving the algorithm maximum freedom to reshape the landscape at each layer.

Recent work in QAOA has shown that increasing the number of variational parameters per layer can improve convergence and reduce the required circuit depth8,9. A similar principle applies to PAOA: more expressive schedules allow the system to reach useful non-equilibrium distributions in fewer layers (e.g., a single layer for a fully parameterized AND gate, as we show in Supplementary Figs. S2, 3), but may also increase the risk of overfitting. The balance between expressiveness and generalization remains a key consideration when selecting an ansatz.

Implementation via probabilistic computers
The practical implementation of the sampling inner loop is achieved using a probabilistic computing (p-computing) framework19. This framework uses networks of binary stochastic units, p-bits, that sample from Boltzmann-like distributions through asynchronous dynamics. A single p-bit updates its state mi according to

mi = sgn[tanh(βIi) - randu(-1, 1)] (1)

where randu(-1, 1) is a random variable uniformly distributed in the range [-1, 1], β is the inverse temperature from the variational ansatz, and the local input Ii is

Ii = Σj Jij mj + hi (2)

As long as connected p-bits are updated sequentially, in any random order, this update rule generates samples from a Boltzmann distribution P({m}) ∝ exp[-βE({m})] over time, associated with the energy

E({m}) = - Σi<j Jij mi mj - Σi hi mi (3)

Formalism and connection to Markov chains
While in practice, PAOA relies on direct MCMC sampling, the process has a rigorous foundation in the theory of Markov chains. The evolution of the probability distribution ρ over the state space can be described by a sequence of transition matrices. To formalize this, we consider a p-process where the distribution after k layers, ρk, is given by

ρk = W(k)(θ(k)) ... W(1)(θ(1)) ρ0 (4)

where ρ0 ∈ R2N is the initial distribution and each W(k) is a 2N × 2N transition matrix parameterized by the variational parameters θ(k) for that layer. For a sequential update scheme, each W(k) can be factorized as a product of single-site update matrices, W(k) = w(k)1 ⊗ ... ⊗ w(k)N. The exact construction of these matrices from the p-bit update rule is detailed in Supplementary Section 3.

This layered application of stochastic matrices to a probability vector ρ is the direct classical analog of applying unitary operators U(k) to a wavefunction ψ in QAOA. The crucial distinction, however, is that each W(k) is a norm-one-preserving stochastic matrix that mixes non-negative probabilities, whereas each U(k) is a norm-two-preserving unitary matrix that rotates complex vectors.

The absence of complex phases in the classical evolution means that PAOA proceeds without the possibility of quantum interference. Consequently, for optimization problems where a constructive interference mechanism for QAOA is unclear, its advantage over classical alternatives is not guaranteed. Indeed, as we will demonstrate in the PAOA versus QAOA: Sherrington-Kirkpatrick model subsection, our benchmarking shows that PAOA consistently reaches higher approximation ratios on the Sherrington-Kirkpatrick model when compared to QAOA using an identical number of variational parameters.

The ultimate objective of PAOA is to find the states {m} that minimize the problem’s energy function (Eq. (3)). This is achieved through the two-level optimization loop described earlier. The outer loop minimizes a cost function (different from the problem’s energy function) over the variational parameters θ, thereby reshaping the sampling landscape to make the optimal states {m} more probable and easier to find through MCMC sampling in the inner loop.

For a target set of states X, this cost function can be the negative log-likelihood:

L(θ) = - Σ{m}∈X ln(ρp({m}; θ)) (5)

In practice, as the exact computation of ρp is intractable for large N, it is replaced by an empirical distribution ρb p estimated from NE independent MCMC samples. Because the samples are generated from a shallow Markov process (small p), the system typically remains out of equilibrium. This non-equilibrium character distinguishes PAOA from classical annealing methods and may unlock new features such as initial condition dependence and faster convergence to the solution. However, this may also become a liability for shallow circuits, as strong dependence on the initial distribution ρ0 or overfitting may prevent the learned heuristics from generalizing. Therefore, the same ρ0 used during training should also be used during inference. Deep PAOA circuits do not exhibit any initial condition dependence, especially with a small initial β that randomizes spins (see the Discovering SA with PAOA subsection).

Finally, the PAOA framework is general in two key respects. First, the cost function being minimized is not restricted to the two-local form, such as the Ising energy of Eq. (3); it can be any computable function over the states {m}, including higher-order k-local Hamiltonians or likelihood functions for machine learning tasks. Second, the variational landscape used for sampling, while implemented here with an Ising-like p-computer, is also not fundamentally limited. The PAOA approach is compatible with any model from which one can efficiently draw samples via MCMC, such as Potts models.

Representative example: majority gate
To illustrate the behavior of PAOA under both exact and approximate dynamics, we solve a small optimization problem involving a four-node majority gate. This problem serves as a tractable testbed for understanding the role of the classical optimizer.

The majority gate is defined over four binary variables [A, B, C, Y], where the output Y is given by

Y = MAJ(A, B, C) = A ∨ B, if C=1; A ∧ B, if C=0 (6)

The eight correct input-output combinations, shown in the truth table in Fig. 2a, define the target set of states X. The goal of the optimization is to find variational parameters that cause the final distribution, ρp, to be concentrated on this set.

As a starting graph, we consider a fully connected Ising graph with Jij=+1 and no biases. In this example, the variational parameters are node-specific schedules βi(p), corresponding to Γ=N. To validate the role of the numerical optimizer, we compare two optimization strategies: a full Markov chain with gradient-based optimization and an MCMC-based approximation using the COBYLA algorithm25.

For the gradient-based approach, we use the formulation laid out in the Formalism and connection to Markov chains subsection and optimize over two layers (p=2) with a uniform initial distribution ρ0 = 1/16. The cost is the negative log-likelihood over the eight correct states, and the variational parameters βi(p) are updated using gradient descent with a fixed learning rate η=0.004. Optimization stops when either the gradient norm falls below 10-7 or a fixed iteration budget is reached.

In the MCMC-based approach, the same initial condition is used, but sampling is performed using 107 independent runs. COBYLA is used to update the parameters based on the empirical distribution ρb. The optimization stops based on parameter convergence or a maximum number of function evaluations. Results are shown in Fig. 2b, where the two curves converge to nearly identical minima.

Figure 2c and d show the final distributions under both methods. The agreement between the exact Markov flow and the sampled histogram confirms that derivative-free optimization closely tracks the true gradient.

We also test the same procedure on a 5-node full-adder circuit using the fully-parameterized Jij(p) ansatz. The resulting distributions (see Supplementary Fig. S1) once again demonstrate tight correspondence between exact and approximate optimization at each layer. Although analytical optimization is feasible for small N, the 2N×2N transition matrix becomes intractable for larger systems. The success of COBYLA in this setting confirms that PAOA can be implemented efficiently on general p-computers using standard MCMC.

Interestingly, PAOA can explore parameter regimes not accessible to traditional Ising formulations. For example, optimal β values may be negative, corresponding to negative temperatures without a clear interpretation in the statistical physics-based context. Since equilibrium Boltzmann-sampling is not necessarily the starting point, alternative and more hardware-friendly activation functions could also be explored with PAOA, such as replacing the hyperbolic tangent in Eq. (1) with the error function erf, or with saturating linear functions.

The shallow nature of the Markov chain used in PAOA introduces a dependency on the initial distribution ρ0. In contrast to standard Boltzmann training methods such as Contrastive Divergence26, which seek to approximate the equilibrium distribution through extensive sampling, PAOA optimizes the evolution of the distribution under a fixed initialization. This shortcut avoids the computational burden of long mixing times but restricts the learned dynamics to a specific initialization at shallow depth. While this dependence can hinder generalization at shallow depth, as we demonstrate in the Discovering SA with PAOA subsection, for deep PAOA circuits, initialization dependence is not an important concern.

Discovering simulated annealing with PAOA
An intriguing question is whether PAOA, when restricted to a global annealing schedule with one parameter per timestep, can recover the well-known structure of SA. In this section, we show that the answer is yes. Using a minimal parameterization with a single global inverse temperature, β(k), for each one of the p total layers, PAOA is able to discover SA-like profiles that optimize the average energy on a large three-dimensional (3D) spin-glass instance.

We apply PAOA to a 3D cubic lattice of size L3=63, which yields a sparse bipartite interaction graph. To leverage this structure for massive parallelism on our FPGA-based p-computer, we use chromatic updates21,22. We employ the simplest Γ = 1 global annealing ansatz, where all nodes share the same scalar inverse temperature β(k) at each layer k, forming a schedule of total length p ∈ {5, 10, 15}.

Sampling is performed on hardware using 720 Monte Carlo sweeps (MCS) per layer. We choose a uniform initial distribution, achieved by appending β=0 to the beginning of the schedule. After sampling, configurations are transferred to the CPU, where the CPU evaluates the average energy 〈E〉 using Eq. (3), and COBYLA updates the schedule β(k). In this training loop, the cost function is the average sampled energy, 〈E〉, estimated from NE independent MCMC runs. After training, however, each candidate's schedule is evaluated on a fresh batch of 105 runs and assigned a success probability, defined as the fraction of runs that reach the putative ground-state energy at the final layer. When referring to the best schedule or best β in Fig. 3, we mean the schedule that achieves the highest success probability out of this independent batch. Insets in Fig. 3c show these success probabilities for all 100 optimizer runs, illustrating both the variability across runs and the improvement with circuit depth.

The loop continues until either convergence or a maximum number of energy evaluations is reached. The algorithm is outlined in Algorithm 1, and parameters are listed in Table 1.

The results of this optimization are summarized in Fig. 3. The optimized schedules shift the energy distribution toward lower values (Fig. 3b), increasing the probability of reaching the putative ground state at E=−360, which was independently found using a linear SA (from β=0.1–5) schedule with 106 MCS. The schedules responsible for this improvement are shown in Fig. 3c. Notably, without any constraints on monotonicity, PAOA consistently learns schedules that begin at high temperature (low β) and gradually cool (increase β), resembling classical SA despite starting from a flat β=2 schedule.

As depth p increases, our schedules both cool further and spend more total sweeps at lower temperatures, so the energy histograms naturally put more weight near the ground state. In majorization language, the deeper runs appear to majorize the shallower ones, but we only show energy marginals; therefore, this is suggestive rather than proof. Similar cooling-with-depth trends were observed for QAOA27.

The elevated β at the first step (p=1) for the 15-layer case may be a result of the substantial degeneracy among near-optimal schedules. Many distinct β profiles achieve nearly identical success probabilities (shaded blue lines in Fig. 3c). The inset histograms show that while individual schedules vary across 100 optimizer runs, performance remains stable (even for non-monotonic schedules) with the smallest variance in success probabilities for deeper schedules (p=15).

Constraints could be added to restrict the class of allowed schedules, e.g., enforcing linear, exponential, or monotonic cooling profiles using inequality constraints in COBYLA. We did not use such constraints to avoid biasing the solution toward SA. These results demonstrate that SA emerges as a special case of variational Monte Carlo, recoverable from the global schedule ansatz. Moreover, PAOA proves to be more general, as it points to the possibility of alternative non-monotonic schedules with comparable performance, indicating the flexibility of the algorithm.

PAOA on FPGA with on-chip annealing
To support long annealing schedules in 3D spin-glass problems, we designed a p-computer architecture that performs annealing entirely on-chip, unlike earlier implementations that require off-chip resources for annealing (e.g., ref. 21). The annealing schedule β(k) is preloaded to the FPGA and indexed via a layer counter that advances after a fixed number of MCS. At each step, the current β value is used to scale the synaptic input Ii via a digital signal processing multiplier (see Supplementary Fig. S4).

The full p-layer annealing process runs uninterrupted on the FPGA, and the final spin configurations are read out only after the last layer is complete.

The FPGA operates in fixed-point arithmetic with s{4}{5} precision (where s denotes the sign bit, followed by 4 integer and 5 fractional bits) for both β(k) and Jij. Due to the size of the FPGA we use, we instantiate 10 independent replicas of the L3=63 spin-glass graph, enabling parallel updates of 2160 p-bits per cycle. Our FPGA architecture with on-chip annealing capability achieves about approximately an 800-fold reduction in wall-clock time in the FPGA compared to an optimized CPU implementation with graph-colored Gibbs sampling (Table 2).

PAOA is parallelizable at multiple levels. First, the cost estimator 〈E〉 is computed from NE independent MCMC runs, which can be executed independently on separate hardware resources. Second, the sampling process itself is parallelizable. Consequently, high-throughput implementations are possible using CPUs, GPUs, FPGAs, custom accelerators or nano devices. Our FPGA design is one point in this space and was chosen to accelerate the experiments in this work. The wall-clock numbers in Table 2 are an implementation case study rather than a universal hardware comparison.

All simulations use a uniform initialization, implemented by prepending a β=0 layer to the schedule. This yields random initial conditions and ensures randomized spin states for each replica. Further hardware design details are provided in Supplementary Section 4.

PAOA versus QAOA: Sherrington-Kirkpatrick model
The Sherrington-Kirkpatrick (SK) model defines a mean-field spin glass with all-to-all random couplings, and serves as a canonical benchmark for optimization algorithms. The classical energy function is defined as

E({m}) = - 1/sqrt(N) Σi<j Jij mi mj (7)

where mi ∈ {-1, 1} and Jij ∼ N(0, 1). The exact ground-state energy per spin in the large-N limit is known28.

QAOA has been studied extensively on the SK model7,15,29,30, with recent numerical work suggesting it can approximate the ground state in the infinite-size limit31. To enable a direct comparison, we evaluate PAOA using a two-schedule ansatz with exactly 2p parameters, matching the parameter count of depth-p QAOA. Although PAOA scales readily to much larger N, we stick to a fixed size of N=26, a typically studied SK-model, to enable a direct comparison with QAOA.

Each schedule is assigned to one half of the graph, split uniformly at random. For N=26, we randomly generate 30 training instances and optimize PAOA separately for each instance and layer depth p ∈ {1, ..., 17}. The resulting schedules are then averaged and applied without further adjustment to a disjoint set of 30 test instances. QAOA results plotted in Fig. 4 use optimal parameters from prior work7,29.

For each instance, exact ground states are computed by exhaustive search.

Our results show that PAOA has highly competitive performance under these iso-parametric conditions. As shown in Fig. 4a (left, PAOA), the average energy per spin consistently decreases with layer depth, and schedules trained on N=26 generalize well to much larger instances up to N=500 (see Supplementary Fig. S5). To formalize this comparison, we report the approximation ratio in Fig. 4b, defined as

Approx. Ratio = E({m}) / E({m})sol (8)

where E({m})sol is the cost of the exact ground state. PAOA achieves a consistently higher approximation ratio than standard QAOA across all depths. Notably, the performance of vanilla PAOA is comparable to that of recently proposed hybrid quantum-classical methods that use classical techniques to improve QAOA’s performance32, suggesting that a definitive quantum advantage on this problem has yet to be firmly established.

We emphasize that a clear quantum advantage over classical algorithms is expected to emerge in problems where interference plays an explicit role, such as in the synthetic examples constructed by Montanaro et al.12. In these cases, the superposition of many computational paths, enabled by either the problem structure or the quantum algorithm, creates interference patterns that amplify the correct solution while suppressing others, an effect that classical probabilistic samplers struggle to reproduce due to the well-known “sign problem33”. As shown in ref. 34, a probabilistic formulation of generic unitary evolutions yields complex (often purely imaginary) weights, so sampling effectively becomes randomized. While post-phase corrections recover unbiased estimates, the number of samples grows exponentially.

These interference-dominated regimes are precisely where we do not expect PAOA, or any classical sampler, to compete. On the other hand, unless such interference-driven advantages are clearly identified and demonstrated more broadly across practical problem classes, we do not expect any significant quantum advantage from QAOA. Meanwhile, PAOA offers a rigorous, scalable and high-performing benchmark for variational sampling at problem sizes beyond the reach of current quantum devices.

Multi-Schedule annealing in the Lévy SK model
To a heuristic solver, the SK model presents no intrinsic structure beyond a dense random coupling graph. To explore whether PAOA can adaptively exploit structure when it might exist, we consider a variant where the coupling weights Jij follow a heavy-tailed Lévy distribution35:

P(J) = |J|-(1+α)/2, |J|>1 (9)

We use α=0.9, which produces a distribution with diverging variance, dominated by rare, large couplings.

In this regime, a natural question arises: can PAOA discover principles that treat strongly and weakly coupled nodes differently? Prior work on Boltzmann machines and annealing schedules36,37 suggests that high-degree or strongly coupled nodes benefit from slower annealing (more time for high temperature) to avoid freezing too early. In these prior works, the heuristics are hand-designed by the ingenuity of the human algorithm designers. Here, we test whether the principle of annealing subgraphs with different temperatures can be discovered by PAOA through black-box optimization.

We begin by introducing a coupling strength per node:

λi = Σj |Jij| (10)

We then sort nodes by λi and split into a heavy group (top 20%) and a light group (bottom 80%), as shown in Fig. 5a. The split here is arbitrary and we confirmed that a 50–50% partition following the same procedure yields similar results (see Supplementary Fig. S6).

To test if PAOA can leverage this structure, we first identify a near-optimal single-schedule using SA at shallow depth (p=17). We then let COBYLA optimize two-schedules, corresponding to heavy and light nodes, freely. After convergence, we average the schedules across all instances. As shown in Fig. 5b, the optimized schedules exhibit a qualitative separation: heavy nodes are assigned a higher-temperature (lower β) profile, and light nodes anneal faster. We verified that the different schedules for heavy and light nodes lead to higher success probability for the double-schedule over the single-schedule. Due to the shallow depth (p=17) however the differences are not appreciably large.

However, when extrapolating the double-schedules discovered by PAOA to much deeper layers (Fig. 5c), we observe a clear difference against single-schedule SA. Our extrapolation is relatively straightforward but we report the details in the Methods Section. The results, shown in Fig. 5d for N=50 and 500 unseen instances, show that the two-schedule ansatz outperforms the single schedule, with a clear gap in success probabilities between the two approaches, confirming a robust benefit from structured, heterogeneous annealing. Details on determining the success probability for single and double schedules are in the Methods Section.

This behavior suggests that PAOA can discover problem-aware heuristics computationally. While such adaptive schedules have long been used manually in structured settings38, we are not aware of prior variational algorithms for non-equilibrium sampling that systematically learn them from data. In this sense, PAOA serves as a tool for discovering new algorithmic strategies in disordered systems.

Outlook
We have introduced the PAOA, a variational Monte Carlo framework that generalizes SA and supports a wider range of parameterizations compatible with existing Ising hardware. By treating the energy landscape itself as a variational object and updating it through classical feedback from Monte Carlo samples, PAOA enables efficient sampling from non-equilibrium distributions using shallow Markov chains.

We demonstrated that PAOA captures SA as a limiting case and, without constraints, can rediscover SA-like behavior when optimized for energy minimization. Using a minimal global schedule ansatz, PAOA efficiently learned monotonic cooling profiles that converge to low-energy states on 3D spin-glass problems. The approach was implemented on an FPGA-based p-computer, achieving hardware-accelerated annealing at rates exceeding CPU-based implementations by several orders of magnitude.

Beyond global schedules, we explored richer ansätze with multiple temperature profiles. On the Sherrington-Kirkpatrick (SK) model, PAOA exceeded the performance of QAOA under iso-parametric conditions, while scaling to system sizes beyond those accessible to quantum devices. In a heavy-tailed variant of the SK model, PAOA automatically learned to assign slower annealing schedules to strongly coupled nodes, reproducing hand-crafted heuristics from prior work. As such, PAOA can be a tool for discovering new algorithmic strategies for sampling in disordered systems.

Methods
Numerical simulation uses double precision (64-bit) in C++ to generate results for PAOA in Figs. 2d and 4a, and Supplementary Fig. S1d. QAOA results were reproduced from ref. 39.

Data transfer between FPGA and CPU
A PCIe interface was used to communicate between FPGA and CPU through a MATLAB interface for the “read/write” operations. A global “disable/enable” signal broadcast from MATLAB to the FPGA was used to freeze and resume all p-bits. Before a “read” instruction, the p-bit states were saved to the local block memory (BRAM) with a snapshot signal. Then the data were read once from the BRAM using the PCIe interface and sent to MATLAB for post-processing, that is, updating the annealing schedule (β(k)). For the “write” instruction, the “disable” signal was sent from MATLAB to freeze the p-bits before sending the updated schedule. After the “write” instruction was given, p-bits were enabled again with the “enable” signal sent from MATLAB.

Measurement of wall-clock time and flips per nanosecond
To measure the annealing wall-clock time per replica in FPGA along with the reading and writing overhead, we use MATLAB’s built-in tic and toc functions and average the total annealing time over 100 independent experiments. Adding all flip attempts for the entire system with all replicas, 2160 p-bits, the total flips per nanosecond (flips/ns) are computed. For CPU measurements, MATLAB’s built-in tic and toc functions were used to measure the total annealing time taken to perform the MCS in Table 2 with 1000 runs each. The average time per MCS is reported for a single replica for a network size of L3=63.

Schedule generation for the Lévy SK model benchmark
The performance benchmark in the Multi-Schedule annealing in the Lévy SK model subsection required comparing two annealing strategies. The schedules for these strategies were generated via a grid search based on a geometric functional form:

β(k) = βinitial * (βfinal/βinitial)^(k/(p-1)), k=0, ..., p-1 (11)

The parameter ranges searched were β ∈ [0.005, 0.12] for the N=50 instances.

The two strategies were constructed as follows:

Single-Schedule (SA): We performed a grid search over βfinal to find the single geometric schedule that produced the lowest average energy across all instances. This served as the baseline.
Two-Schedule PAOA: The optimized baseline schedule was assigned to the heavy nodes, as suggested by PAOA. For light nodes, we scale up the baseline schedule by a scaling factor that reflects the separation proposed by PAOA. Then, both geometric schedules are extrapolated to large number of layers to study the effect of deep PAOA.

In order to compare the two approaches, we define success probability as the ratio of the number of runs that hit the putative ground state to the total number of runs (i.e., 100 independent runs). The putative ground state was found by running a long SA with one million sweeps and ten independent runs, then choosing the lowest energy among all ten million states.

Optimizer choice
For the outer optimization loop, we used the COBYLA algorithm as it is a robust and widely used method for derivative-free optimization problems with a moderate number of variables40–42. Our COBYLA implementation is based on the non-linear optimization package43. While the performance of PAOA may depend on the choice of classical optimizer, a detailed comparison of different optimization techniques is beyond the scope of this work.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors discuss the Probabilistic Approximate Optimization Algorithm (PAOA), which is a classical probabilistic\nvariational circuit inspirited by Quantum Approximate Optimization Algorithm (QAOA). PAOA was introduced by Weitz et al\nin a paper published in PRA earlier this year. The authors implement the algorithm on an FPGA-based probabilistic\ncomputer and demonstrated significant speedups over other approaches including QAOA for a small instance of\nSherrington–Kirkpatrick model. The authors note that their algorithm can discover annealing schedules that are better than\nsingle-schedule simulated annealing.\nThe authors generalize the original PAOA framework and realize it experimentally with FPGAs. They demonstrate the\nalgorithm on graphs with up to 500 vertices. While this looks impressive, I see several problems with the paper that would\nneed to be sorted out before it is accepted in Nature Communications.\nIn my opinion, the advantages that the authors claim are not yet fully validated in terms of scalability, generalization and\ncomparison fairness. The fundamental weakness of the approach as compared to simulated annealing is strong\ndependence on the initial distribution. This can limit the generalization, as the authors admit. Heuristics learned for one\ninitialization may not transfer to others. In the present manuscript, it is not clear to what extend the advantages of PAOA are\nreal or are they the artifacts of carefully chosen initial conditions.\nI find the QAOA comparison misleading. PAOA has already been compared to QAOA by Weitz et al, so why the authors\nrepeat similar discussion here? Did the authors find an example that shows QAOA superiority? It would be good to see\nexamples that point towards QAOA advantage as well. The message coming from the PAOA-vs-QAOA section seems to be\n“PAOA is a practical baseline” – this has been established before.\nI find parts of the presentation ambiguous. The authors repeatedly claim that simulated annealing emerges as a limiting case\nof their approach. Is PAOA then just a generalization of simulated annealing or is there a new class of problems that PAOA\ncan solve better? I am trying to understand whether this approach is fundamentally different or an incremental improvement\nover old techniques. The paper fails to convincingly establish the approach as a genuinely new paradigm. So far, the\npresentation and results suggest that the method is an optimized heuristic that may work on carefully chosen benchmarks\nand only when properly initialized.", "reviewer_2": "Referee Report on Abdelrahman et al. \"Probabilistic Approximate Optimization: A New Variational Monte Carlo Algorithm\"\nThe paper introduces the Probabilistic Approximate Optimization Algorithm (PAOA), a new variational Monte Carlo\nframework. PAOA was inspired by the Quantum Approximate Optimization Algorithm (QAOA) and thus is a \"dequantized\"\nversion of QAOA. The authors show that PAOA generalizes simulated annealing and can rediscover annealing schedules.\nThey demonstrate hardware acceleration on an FPGA p bit computer with large performance gains over CPUs. PAOA is\nbenchmarked against QAOA on the Sherrington–Kirkpatrick model and achieves higher approximation ratios under matched\nparameters. The work is significant because it provides a scalable classical baseline and a tool for discovering new non-\nequilibrium heuristics, with impact on optimization and probabilistic computing.\nI also wanted to point out, in the spirit of scientific reproducibility, the authors have open-sourced their GitHub repository,\nwhich is commendable. I took a cursory look at the code.\nI think the work is clear, careful, and highly novel. I strongly support publication of this work. However I would like the authors\nto respond to the comments/clarifications/questions below.\nComments, Clarifications, Questions\n------------------------------------\nMost of my comments are aimed at improving clarity and presentation. The authors should feel free to implement only those\nchanges they believe clarify or strengthen the paper. If a suggested change is marginal, it need not be implemented.\n1. Definition of Symbols and Notation\n- Define N in the abstract and again in Section II. From my understanding\n- N = number of spins/variables in the optimization problem.\n- 2^N = size of the full state space.\n- In Fig. 1, explicitly define symbols: ρ, hat ρ, p, M, m.\n- Clarify upfront why M (ansatz size, possibly with hidden variables) differs from N (problem size). This was explained later in\nSection II.C, but until I read that part I was confused.\n2. Section IV and Fig 3\n- In Section IV, is \"betst beta\" the optimization criterion that minimization of expected energy, i.e. best ≡ min ⟨E⟩.\n- In the inset figures in Fig 3 (c). State / clarify that ⟨p⟩ in plots refers to an expectation over independent runs. That if i was to\ncompute the average of the numbers shown that i'd get the average number.\n- Numerical Anomalies? In Fig. 3, the higher energy observed at p = 1 for 15 layers is unexplained. Is this just a numerical\nartifact? It may or may not Some comment would help the reader.\n- I noticed in Fig. 3(b) an apparent majorization ordering across circuit depths. Ordering probabilities from largest to smallest\ngives $P_{15} \\succ P_{10} \\succ P_{5}$, suggesting structured improvement beyond average energy. Since entropy and\npurity are Schur-convex/concave, one expects $H(P_{15}) < H(P_{10}) < H(P_{5})$. Entropy and purity from truncated data\nmay not show this, but the full distributions likely would. This indicates that adding layers reduces entropy and concentrates\nprobability on the correct answer. This may or may not be an interesting additional insight for the manuscript.\n3. The y-axis in Fig. 4.\n-The scale seems a little off. Is it logarithmic, or is there a label typo from illustrator editing? The spacing between |-0.8 -\n(-0.6)| \\neq |-0.6 - (-0.4)|. Regardless I wonder if cliping the y axis to -0.3 might allow for the data to look less compressed. Or\nrescaling (possibly log scale) to improve readability. [This also applies to Figure S5]\n4. Fig. 5 e.\nI'm trying to get a feel for the improvement of PAOA over SA. Are the differences shown meaningful? Or am I missing the\npoint of what is illustrated by the plots here?\n5. Reference [31] should be corrected/clarified: Make it clear that it is Chapter 10 in Simulated Annealing, Emile Aarts, Jan\nKorst, and Wil Michiels."}, "author_response": "We sincerely thank the reviewer for their thoughtful review and the open-minded assessment. We believe we have\naddressed each of the points in the revised manuscript, and a color-tracked diff file highlighting every change is\nappended for the reviewer’s convenience.\n(1) The fundamental weakness of the approach as compared to simulated annealing is strong dependence on the initial\ndistribution. This can limit the generalization, as the authors admit. Heuristics learned for one initialization may not\ntransfer to others.\nAUTHOR RESPONSE\nWethankthereviewerforhighlightingtheroleofinitialization. Ouroriginaltextdidnotclearlyseparatethebehavior\nof shallow versus deep PAOA circuits. We would like to make a few clarifying points:\n(i) PAOA strictly contains simulated annealing (SA) as a special case (global schedule with one parameter per layer),\nwhile also supporting richer, structured parameterizations.\n(ii)Initializationdependenceisnot aninherentlimitationofPAOA.Itonlyarisesinshallow,non-equilibriumsettings.\n• Shallow PAOA. For small depth p, PAOA implements a short Markov flow that can depend on the chosen\ninitialization ⇢ . This can be advantageous for targeted sampling, e.g., in the majority-gate example (Fig. 2) a\n0\nspecific initial state is steered into a distribution peaked on the correct truth-table entries.\n• Deep PAOA. With su cient layers and sweeps, the learned schedules start hot (small  ), which rapidly\nmixes the chain and erases memory of ⇢ before cooling. In our 3D spin-glass experiment (Fig. 3), PAOA\n0\nrecoversSA-likebehaviorwithout imposingmonotonicityconstraintson ,whichisknowntobeinitialcondition\nindependent.\n2\n(iii)PAOA’smainadvantageisexpressivity. Asanexample,intheSK–L´evymodel,PAOAdiscoversadouble-schedule\nansatz that assigns a slower (higher-temperature) annealing schedule to “heavy” nodes and a faster one to “light”\nnodes. Thislearned heterogeneityoutperformsahighlyoptimizedsingle-scheduleSA(revisedFig.5andSectionVII).\n(iv) The learned schedules already generalize across instances and sizes without retraining. In Fig. 4, schedules\ntrained on one set of n=26 SK instances are applied to a disjoint test set; the same schedules also generalize to\nn 400,500 (Fig. S5). Likewise, in Fig. 5, the SK–L´evy heuristics learned at small n are used at larger n without\n2{ }\nretraining.\nAUTHOR ACTION\nWe revised the second to last paragraph of Section IIB, where we discuss initialization, as follows:\nBecause the samples are generated from a shallow Markov process (small p), the system typically remains\nout of equilibrium. This non-equilibrium character distinguishes PAOA from classical annealing methods\nand may unlock new features such as initial condition dependence and faster convergence to the solution.\nHowever, this may also become a liability, for shallow circuits, as strong dependence on the initial distri-\nbution ⇢ or overfitting may prevent the learned heuristics from generalizing. Therefore, the same ⇢ used\n0 0\nduring training should also be used during inference. As we show in later Section IV, deep PAOA circuits\ndo not exhibit any initial condition dependence especially with small initial   that randomizes spins.\nPlease also see the last paragraph of Section III where we clarified the initial condition dependence between shallow\nand deep PAOA circuits.\n(2) In the present manuscript, it is not clear to what extend the advantages of PAOA are real or are they the artifacts\nof carefully chosen initial conditions.\nAUTHOR RESPONSE\nInourwork,wepresentedtworesultswherePAOAdemonstratesadvantagesoverexistingmethodsandneithershows\nany dependence on initial conditions. They are robust and real.\nThe first result is our comparison against QAOA in Fig. 4, which was explicitly designed to test for generalization.\nThe protocol involved optimizing the PAOA schedules on a dedicated training set of instances. We then applied the\nlearned average schedule to a fully disjoint set of 30 random test instances without any retraining or fine-\ntuning. The superior approximation ratio achieved by PAOA across all depths (Fig. 4b) is therefore a demonstrated\nadvantage on unseen problem instances, not an artifact of carefully chosen conditions.\nThe second result is our analysis of the SK-L´evy model in Fig. 5, now completely revised, which shows that these\nadvantages hold in deep circuits where initial condition dependence is not a factor:\n1. Irrelevance of initial conditions: For the deep circuits used in this test, both the SA and PAOA schedules\nbegin at very low inverse temperatures (  0), corresponding to an initial high-temperature phase. This\n⇡\nstandard annealing practice ensures the system is thoroughly randomized, erasing any memory of the specific\nstarting spin configuration and making the subsequent performance independent of the initial state.\n2. Generalization across scale: The advantage of the multi-schedule heuristic discovered by PAOA is not tied\nto a specific problem instance or size. The strategy of using distinct schedules for “heavy” and “light” nodes\nconsistently outperforms the best single-schedule SA. This demonstrates that the advantage is a robust feature\nof the learned strategy itself.\nInbothcases,theadvantagesofPAOAareshowntobegenuineandgeneralizable,notartifactsoftheinitialconditions.\nPlease see Author Action below to see the concrete changes we made to the manuscript to better clarify these points.\nAUTHOR ACTION\n• Please see the second to last paragraph of Section IIC\n3\n• Please see our completely revised Fig. 5, Section VII, as well as Fig. 4.\n(3) I find the QAOA comparison misleading. PAOA has already been compared to QAOA by Weitz et al, so why the\nauthors repeat similar discussion here?\nAUTHOR RESPONSE\nOur comparison of PAOA to QAOA is unique and is not a repeat of prior work. The comparison by Weitz et al.\nfocused on the MaxCut problem for specific graph types. To the best of our knowledge, our work represents the\nfirst direct comparison of PAOA and QAOA on the Sherrington-Kirkpatrick (SK) model, which is a canonical,\nmean-field spin glass problem and a celebrated benchmark in numerous foundational QAOA studies. A rigorous\nperformance comparison on this specific problem is therefore a critical and previously unaddressed question.\nFurthermore, our work is not a repetition but a substantial generalization that validates and extends the PAOA\nframework to new problem domains and into a scalable, hardware-realized algorithm, as the original authors of\nRef.[1], Weitz et al., explicitly acknowledge in their own publication when citing our work (boldface ours):\n“Second, Abdelrahman et al. [44] extend our theoretical framework and realize the PAOA exper-\nimentally. They developed hardware capable of executing the PAOA, including the classical optimizer,\non graphs with up to 500 vertices. This provides independent experimental validation of the PAOA and a\nsubstantial generalization of our theoretical results.”\n(4) Did the authors find an example that shows QAOA superiority? It would be good to see examples that point\ntowards QAOA advantage as well.\nAUTHOR RESPONSE/ACTION\nWhile our objective is to conduct a rigorous, iso-parametric benchmark of PAOA on a canonical (and not previously\nstudied) QAOA benchmark, we agree that a fair assessment with respect to QAOA is absolutely essential.\nWe now expand the last paragraph of Section VI in the revised manuscript to better illustrate the fundamental\npossible advantages of QAOA where we state that\n“...thataclearquantumadvantageoverclassicalalgorithmsisexpectedtoemergeinproblemswhereinterferenceplays\nanexplicitrole,suchasinthesyntheticexamplesconstructedbyMontanaroetal[12]. Inthesecases,thesuperposition\nof many computational paths, enabled by either the problem structure or the quantum algorithm, creates interference\npatternsthatamplifythecorrectsolutionwhilesuppressingothers,ane↵ectthatclassicalprobabilisticsamplersstruggle\nto reproduce due to the well-known “sign problem” [32]. As shown in [33], a probabilistic formulation of generic\nunitaryevolutionsyieldscomplex(oftenpurelyimaginary)weights, sosamplinge↵ectivelybecomesrandomized. While\npost-phase corrections recover unbiased estimates, the number of samples grows exponentially.\nTheseinterference-dominatedregimesarepreciselywherewedonotexpectPAOA,oranyclassicalsampler,tocompete.\nOn the other hand, unless such interference-driven advantages are clearly identified and demonstrated more broadly\nacross practical problem classes, we do not expect any significant quantum advantage from QAOA. Meanwhile, PAOA\no↵ers a rigorous, scalable and high-performing benchmark for variational sampling at problem sizes beyond the reach\nof current quantum devices.\n(5) The message coming from the PAOA-vs-QAOA section seems to be “PAOA is a practical baseline” – this has\nbeen established before.\nAUTHOR RESPONSE\nAs we discussed in our response to (3), PAOA vs QAOA on the canonical SK problem has not been established\nbefore. Moreover, the conclusion the reviewer is drawing, namely, PAOA is a practical baseline is incorrect, as we\ndiscuss in our response to (4): PAOA is a superior baseline (not just a practical one) for the SK problem.\n(6) I find parts of the presentation ambiguous. The authors repeatedly claim that simulated annealing emerges as a\nlimiting case of their approach. Is PAOA then just a generalization of simulated annealing or is there a new class of\nproblems that PAOA can solve better?\nSimulatedAnnealingisindeedasub-caseofthegeneralPAOAframework. However,PAOAisnotjust ageneralization\nof SA; it is a more powerful framework that can solve certain classes of problems better by learning heuristics that\nare inaccessible to standard SA. Our new results in FIG. 5 provide a clear demonstration of this advantage. To\nsummarize:\n• SAisaspecific,constrainedsub-casewithinthebroaderPAOAframework. WhenPAOAisrestrictedtoasingle\nglobal schedule and run for a su cient number of layers, it learns to recover the behavior of SA.\n• ThetrueadvantageofPAOAemergesinproblemswithunderlyingstructure,suchastheSherrington-Kirkpatrick\n(SK)modelwithheavy-tailedL´evybonds. StandardSAislimitedtoapplyingasingle,globalannealingschedule,\nwhich is a suboptimal strategy for such a heterogeneous problem. In contrast, PAOA can utilize a double-\nschedule ansatz. As shown in our new FIG. 5(b), PAOA automatically discovers a problem-aware heuristic\nby learning to assign a slower, higher-temperature schedule to the “heavy” nodes and a faster-cooling schedule\ntothe“light”nodes. Thislearnedstrategyresultsinaclearandstatisticallysignificantperformanceadvantage.\nAs shown in FIG. 5(d), the double-schedule PAOA achieves a consistently higher success probability than the\nbest-performing single-schedule SA, outperforming it on a clear majority of the 500 tested instances.\nWe believe the double-schedule example we show is the tip of the iceberg: one can imagine many other parameteri-\nzations of PAOA, such as private temperatures to each node, improving replica-based algorithms such as simulated\nquantum annealing and parallel tempering, or entirely new heuristics. The main contribution of our work is to\nshow PAOA can provide a systematic method of exploring these methods by approximating derivatives of very high\ndimensional Markov matrices. We hope that our foundational work will motivate such future studies.\nAUTHOR ACTION\nPlease see our revised Fig. 5 and Section VII.\n(7) I am trying to understand whether this approach is fundamentally di↵erent or an incremental improvement over\nold techniques. The paper fails to convincingly establish the approach as a genuinely new paradigm. The presentation\nand results suggest that the method is an optimized heuristic that may work on carefully chosen benchmarks and only\nwhen properly initialized.\nAUTHOR RESPONSE/ACTION\nAs discussed earlier in several of our responses (e.g., (2)-(3)), PAOA does not exhibit initial condition dependence in\ngeneral.\nMoreover, PAOA is an entirely new framework that generalizes Simulated Annealing enabling countless other param-\neterizations in a principled optimization framework.\nOur core contributions over Ref. 1 are:\n• AglobalMarkov-flowformulationandvariationalprinciplelinkingderivative-freeupdatestothegradientofthe\nfull (exponentially large) Markov process.\n• A scalable hardware realization (FPGA) with on-chip annealing and massively parallel sampling.\n• New iso-parametric SK benchmarks versus QAOA with generalization across instances and sizes.\n• Automaticdiscoveryofstructured,heterogeneousschedulesinheavy-tailedSK–L´evythatoutperformsingle-schedule\nSA.\nWe genuinely believe PAOA is a fertile ground. QAOA has a vast literature with thousands of papers; by contrast,\nPAOA is discussed in just a single one so far (excluding ours). We believe PAOA can be a serious contender and\nmerits broader exploration.\nOne final note: We were alerted by the Editor that our current title needs to be modified, as such we are proposing\n“Generalized Probabilistic Approximate Optimization Algorithm” as a new title of the paper.\n5\n==================================================================\nII. REVIEWER #2\n==================================================================\nThe paper introduces the Probabilistic Approximate Optimization Algorithm (PAOA), a new variational Monte Carlo\nframework. PAOAwasinspiredbytheQuantumApproximateOptimizationAlgorithm(QAOA)andthusisa“dequan-\ntized” version of QAOA. The authors show that PAOA generalizes simulated annealing and can rediscover annealing\nschedules. They demonstrate hardware acceleration on an FPGA p-bit computer with large performance gains over\nCPUs. PAOA is benchmarked against QAOA on the Sherrington–Kirkpatrick model and achieves higher approxima-\ntion ratios under matched parameters. The work is significant because it provides a scalable classical baseline and a\ntool for discovering new non-equilibrium heuristics, with impact on optimization and probabilistic computing.\nAUTHOR RESPONSE\nWe sincerely thank the reviewer for their thoughtful review and positive assessment. We believe we have addressed\neach of the points in the revised manuscript, and a color-tracked diff file highlighting every change is appended for\nthe reviewer’s convenience.\n(1) Definition of Symbols and Notation. Define N in the abstract and again in Section II. From my understanding\nN = number of spins/variables in the optimization problem, and 2N = size of the full state space.\nAUTHOR RESPONSE/ACTION\nThe reviewer is correct about the definition. We removed the explicit mention of N from the abstract for clarity and\nadded descriptions in the caption of Fig. 1 in Section II.\n(2) In Fig. 1, explicitly define symbols: ⇢, ⇢ˆ, p, M, m. Clarify upfront why M (ansatz size, possibly with hidden variables) di↵ers from N (problem size). This was explained later in Section II.C, but until I read that part I was confused.\nAUTHOR RESPONSE/ACTION\nWe agree. We have revised the Fig. 1 caption to explicitly define the symbols ⇢, ⇢ˆ, p, M, and m, and clarified why\nthe ansatz size M can di↵er from the problem size N. We have also added an explicit discussion at the beginning of\nSection II so that this distinction is clear upfront rather than deferred to Section II.C.\n(3) In Section IV, is “best beta” the optimization criterion defined as minimization of expected energy, i.e., best\nmin E ?\n⌘ h i\nAUTHOR RESPONSE\nWe agree. In Sec. IV the optimizer updates the schedules by minimizing the average sampled energy E over N\nE\nh i\nindependent runs. After training, each candidate schedule is evaluated on a separate batch of 105 inference runs.\n“Best  ” denotes the schedule whose inference runs yield the highest success probability, where success probability is\ndefinedasthefractionofrunsthatreachtheputativeground-stateenergyatthefinallayer. Thesuccessprobabilities\nshown in the insets of Fig. 3(c) are computed independently for each optimizer run under its corresponding schedule.\nAUTHOR ACTION\nWe have written the third paragraph of Section IV to better clarify these points.\n6\n(4) Numerical anomalies? In Fig. 3, the higher energy observed at p=1 for 15 layers is unexplained. Is this just a\nnumerical artifact? Some comment would help the reader.\nAUTHOR RESPONSE\nThe reviewer is likely referring to the inverse temperature schedule,  (p), rather than to the energy itself. One might\nreasonably expect an optimal   schedule must be monotonically increasing over the layer depth, p. In practice,\nhowever, we observe that there is a large degeneracy of   schedules. The shaded blue lines in Fig. 3 show all the\nschedules found by the optimizer and the insets report the corresponding success probabilities. The near-uniform\nsuccess probabilities indicate that there are many degenerate and near-equivalent schedules. Despite the apparent\noutlieratitssecondpoint,theredscheduleatp=15providesthehighestsuccessprobabilityamong100independent\noptimizer runs.\nAUTHOR ACTION\nWe have revised the discussion around Fig. 3 to include this observation.\n(5) In Fig. 3(b), I noticed an apparent majorization ordering across circuit depths. Ordering probabilities from largest\nto smallest gives P P P , suggesting structured improvement beyond average energy. Since entropy and purity\n15 10 5\n   \nare Schur-convex/concave, one expects H(P )<H(P )<H(P ). Entropy and purity from truncated data may not\n15 10 5\nshow this, but the full distributions likely would. This indicates that adding layers reduces entropy and concentrates\nprobability on the correct answer. This may or may not be an interesting additional insight for the manuscript.\nAUTHOR RESPONSE\nThereviewermakesadeepandinsightfulobservation. Weagreethatthefulldistributions(thatareoutofreach)likely\nwould show explicit majorization. It may also be possible to make mathematical statements about this majorizing\nbehavior in the appropriate limit. We hope that the small number of PAOA papers (starting with the one our work\nbuilt on Ref. [1]) will motivate such studies further.\nAUTHOR ACTION\nWe added the following paragraph to Section IV:\nAs depth p increases, our schedules both cool further and spend more total sweeps at lower temperatures,\nso the energy histograms naturally put more weight near the ground state. In majorization language, the\ndeeper runs appear to majorize the shallower ones, but we only show energy marginals, therefore, this is\nsuggestive rather than proof. Similar cooling-with-depth trends were observed for QAOA Ref[27].\n(6) The y-axis in Fig. 4. The scale seems a little o↵. Is it logarithmic, or is there a label typo from Illustrator editing?\nThe spacing between 0.8 ( 0.6) = 0.6 ( 0.4). Regardless, I wonder if clipping the y-axis to 0.3 might\n|      | 6 |      |  \nallow for the data to look less compressed. Or rescaling (possibly log scale) to improve readability. [This also applies\nto Figure S5]\nAUTHOR RESPONSE/ACTION\nWe agree. We revised the main and supplementary figures and used linear scaling for clarity. Following the excellent\nsuggestion by the reviewer, we also clipped y-axis appropriately for the main and supplementary figures.\n(7) Fig. 5 e. I’m trying to get a feel for the improvement of PAOA over SA. Are the di↵erences shown meaningful?\nOr am I missing the point of what is illustrated by the plots here?\nAUTHOR RESPONSE\nIn light of this comment, we have completely revised Fig. 5 and the associated Section VII. Even though the\nimprovement of double-schedule PAOA over standard (single schedule) SA is statistically significant in our original\nshallow depth (p) we used: the di↵erences between the anneals did not have enough “time” to grow.\nTo better show whether the double-schedule provided by PAOA improves SA at deeper layers, we extrapolate the\nschedules provided by the optimizer, using a straightforward procedure. With more depth, we observe a clearer\nincrease in success probability (from about 16% to 20%) for the double-schedule PAOA over the single SA, for this\nclass of heavy-tailed SK problem.\nAUTHOR ACTION\nPlease observe the new Fig. 5 in the marked up manuscript along with the discussion.\n(8) Reference [31] should be corrected/clarified: Make it clear that it is Chapter 10 in Simulated Annealing, Emile\nAarts, Jan Korst, and Wil Michiels.\nAUTHOR RESPONSE/ACTION\nWe thank the reviewer for catching a serious citational error: our original reference was to another book by Emile\nAarts and Jan Korst. We now corrected the citation in the revised manuscript.\nOne final note: We were alerted by the Editor that our current title needs to be modified, as such we are proposing\n“Generalized Probabilistic Approximate Optimization Algorithm” as a new title of the paper."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors answered my questions in a satisfactory way. They clarified their claims and improved the presentation of their\nresults. I now support the publication of the manuscript in its present form.", "reviewer_2": "The authors have substantially revised the manuscript and addressed my comments.\nI reread the paper and noticed a few minor typos and inconsistencies, that will likely be caught at the proof stage. I did catch\none issue. You compute an estimated distribution (cid:1)_f but never define f. Define it on first use (e.g.,f = “final layer”?).\nAdmittedly I got side tracked by an interesting observation about figure 5 (d). I spent 3 hours exploring that data. Clearly that\nwas research and not related to the report. I hope many other researchers get side tracked by the wonderful results in this\nmanuscript.\nI recommend publication without further delay. I am looking forward to reading the published version of the article."}, "author_response": "We thank the reviewer for the comprehensive review that has markedly improved the manuscript.\nI did catch one issue. You compute an estimated distribution ρ but never define f. Define it on first use (e.g.,f =\nf\n“final layer”?).\nAUTHOR RESPONSE/ACTION\nWe corrected this notational inconsistency by using ρ all throughout, where p represents the number of layers and ρ\np\nrepresents the estimated distribution."}]

Decision letter