Gaussian Optimization Protocols

• Gaussian Optimization Protocols are a family of techniques that use Gaussian representations to transform complex discrete or constrained spaces into continuous, tractable forms.
• They employ mechanisms like Gaussian parameterization, Gaussian-process uncertainty modeling, and Gaussian sampling to drive optimization in scheduling, black-box, quantum, and geometric domains.
• These protocols balance tractability with expressivity by smoothing non-differentiable states, preserving constraints via geometric updates, and enabling efficient numerical optimization.

“Gaussian Optimization Protocol” is not a single standardized algorithm in the literature surveyed here. Rather, the expression is used for a family of optimization procedures in which Gaussian structure is the central computational object: continuous Gaussian relaxations for discrete scheduling, isotropic Gaussian mutation in derivative-free search, Gaussian-process surrogates for Bayesian and policy optimization, Gaussian operations for refining non-Gaussian quantum resources, Gaussian boson sampling for proportional sampling, and symplectic or geometric optimization directly on Gaussian-state manifolds (Cai et al., 23 Feb 2026, Wong, 2023, Jing et al., 19 Sep 2025, Arrazola et al., 2018, Willby et al., 28 Jan 2026). Across these settings, the common pattern is to replace an intractable search space or constrained state space by Gaussian variables, Gaussian posteriors, Gaussian transformations, or Gaussian sampling rules that preserve enough structure to make optimization numerically feasible.

1. Terminological scope and recurring design principles

In the cited literature, the term denotes optimization protocols built around one of four recurring Gaussian mechanisms. The first is Gaussian parameterization, in which decision variables are represented by Gaussian means and variances rather than categorical or combinatorial assignments. The second is Gaussian-process uncertainty modeling, where optimization is performed through posterior means, variances, and acquisition functions. The third is Gaussian transformation of physical states, especially in continuous-variable quantum information and pure Gaussian-state geometry. The fourth is Gaussian sampling as a search primitive, either through isotropic mutations in derivative-free search or through hardware-native sampling distributions such as Gaussian boson sampling (Cai et al., 23 Feb 2026, Huang et al., 2019, Jing et al., 19 Sep 2025, Arrazola et al., 2018).

A second recurring principle is that Gaussian structure is used to balance tractability and expressivity. In scheduling, Gaussian start-time variables capture the ordinal nature of time while reducing parameter count from $O(|V|\cdot D)$ to $O(|V|)$. In Bayesian optimization, Gaussian processes provide closed-form posteriors and acquisition functions. In Gaussian-state optimization, symplectic or orthogonal group actions preserve physical constraints exactly, converting constrained variational problems into unconstrained optimization over structured generators (Cai et al., 23 Feb 2026, Candelieri et al., 18 May 2025, Windt et al., 2020, Willby et al., 28 Jan 2026).

A third common feature is the use of explicit post-processing or geometric preservation. Some protocols round and legalize continuous Gaussian optima into executable schedules; some preserve feasibility by conditioning a Gaussian process on boundary or state constraints; others remain inside the physical manifold by updating along symplectic or orthogonal transformations rather than in Euclidean parameter space (Cai et al., 23 Feb 2026, Khatab et al., 16 Mar 2026, Windt et al., 2020).

2. Gaussian relaxation and stochastic search in discrete optimization

The most explicit protocol-style formulation under this name appears in operator scheduling. GauS models start times on a DAG $G=(V,E)$ by independent Gaussian variables,

$s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$

with reparameterization

$$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$

This replaces categorical time assignments by continuous Gaussian variables, captures the ordinal nature of time, reduces the parameter space to exactly two parameters per operator, and enables low-variance, unbiased gradient estimation of expected objectives and penalties. The protocol smooths non-differentiable quantities such as makespan, peak resource usage, and peak memory with LogSumExp surrogates, uses Gaussian CDF bucket probabilities to align the continuous relaxation with discrete rounding, and combines objective terms with differentiable penalties and augmented Lagrangian updates. It also provides a differentiable formulation of modulo or pipelined scheduling through phase probabilities $P_i^{\mathrm{mod}}(t)$, modulo resource usage, modulo memory footprint, and recurrence penalties. After optimization, means are rounded by $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$ and then legalized by topological or fixed-point passes. On EPFL and large synthetic workloads, GauS is reported to achieve Pareto-optimal trade-offs, a 71.8% geometric mean improvement over GS-Schedule on shared feasible instances for a resource+communication objective, near-100% GPU utilization, and $10\text{–}100\times$ faster convergence when objectives coincide, while categorical baselines often OOM on large graphs (Cai et al., 23 Feb 2026).

A different use of the same phrase appears in Gaussian Crunching Search, a derivative-free optimizer for unconstrained black-box minimization. GCS samples a single candidate from an isotropic Gaussian centered at the current iterate,

$$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$

accepts it if $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$, resets the variance on success, and expands it by $O(|V|)$0 on failure. The “crunching” step is the hard reset $O(|V|)$1 after improvement; otherwise $O(|V|)$2. The paper recommends $O(|V|)$3 and $O(|V|)$4, and positions the method for plateau-like or noisy landscapes where objective differences are extremely small relative to numerical precision. On the benchmark $O(|V|)$5, all tested SciPy methods failed from $O(|V|)$6, whereas GCS achieved a passing objective $O(|V|)$7; the reported 100-iteration failure probability remained $O(|V|)$8 from $O(|V|)$9 to $G=(V,E)$0, rising to $G=(V,E)$1 at $G=(V,E)$2 (Wong, 2023).

These two protocols illustrate a sharp split within the term. In GauS, Gaussianity is a differentiable relaxation of a discrete structured problem. In GCS, Gaussianity is the mutation law of a stochastic hill-climber. The shared feature is not objective class or solver architecture, but the decision to make Gaussian sampling or Gaussian parameterization the primary search mechanism.

3. Gaussian-process surrogate optimization

A large fraction of the literature uses “Gaussian optimization” to mean optimization mediated by Gaussian-process posteriors. In these protocols, the search is driven by posterior mean, posterior variance, and acquisition rules rather than by direct optimization of the original expensive objective.

In multi-objective design optimization, the multi-objective multi-generation Gaussian process optimizer fits an independent GP for each objective, uses a squared exponential kernel with ARD lengthscales, and evaluates candidate offspring by lower confidence bounds $G=(V,E)$3. Candidates are then Pareto-ranked in LCB space, after which only the best are truly evaluated. The reported default schedule decays $G=(V,E)$4 from 2 by $G=(V,E)$5, uses $G=(V,E)$6, and screens $G=(V,E)$7 candidates per generation. On ZDT benchmarks with $G=(V,E)$8 and $G=(V,E)$9, the method is reported to converge faster and more stably than NSGA-II, MOPSO, MMOPSO, and WOF-SMPSO, with substantially better HV and IGD across most checkpoints (Huang et al., 2019).

In reinforcement learning, Gaussian Process Policy Optimization augments PPO with a GP defined over policy parameters $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$0, using the GP posterior to compute an expected-improvement bonus added to the clipped PPO surrogate. The combined loss is

$s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$1

and the paper reports training with no entropy bonus, a two-hidden-layer $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$2-unit tanh MLP, GP memory size $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$3, and fixed kernel settings of length-scale $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$4 and noise $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$5. On MuJoCo tasks, GPPO is described as comparable to PPO and TRPO and superior on Ant-v2 and HalfCheetah-v2 (Rao et al., 2020).

Several works modify the GP model itself rather than only the acquisition. WBGP-BO replaces single-hyperparameter GP fitting by a fixed ensemble of GPs and combines their univariate predictive Gaussians pointwise via the Wasserstein barycenter,

$s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$6

before computing LCB or other acquisitions. The motivation is that MLE instability in standard GP fitting can mislead Bayesian optimization; the reported empirical result is that the barycenter model converges on “tricky” one-dimensional problems where vanilla GP-BO does not (Candelieri et al., 18 May 2025). CGLO and AGLGP address a different weakness by combining a sparse inducing-point global GP with region-specific local GPs, alternating between a global density-penalized expected-improvement step and a local modified expected-improvement step, then allocating replications by OCBA in the current region (Meng et al., 2021).

Other variants adapt GP-based optimization to nonstandard computational settings. Hyper-optimization with Gaussian Process and Differential Evolution uses kernel coercion inside GP covariance evaluations to handle integer, categorical, and boolean variables and then uses differential evolution to optimize EI in the resulting non-smooth acquisition landscape; in the BlackBox 2020 challenge, its best reported public-test-set score under LH priming and quasi-random meta-initialization was 97.871 (Klus et al., 2021). Deterministic global optimization of trained Gaussian processes instead treats the GP posterior and acquisition functions as explicit factorable programs, propagates McCormick relaxations through reduced-space formulations, and globally optimizes EI, PI, or LCB via branch-and-bound; the reported runtime scaling is near-linear in the number of training points in reduced space, in contrast to much heavier full-space formulations (Schweidtmann et al., 2020).

Quantum kernel surrogates extend the same template. Quantum Gaussian Process Regression for Bayesian Optimization replaces the classical kernel by a quantum overlap kernel $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$7, emphasizes PSD-preserving regularization of noisy Gram matrices by eigenvalue cutoff, and reports that statevector and sample-based quantum Bayesian optimization match classical GP-BO on Branin and on hyperparameter optimization of a gradient boosting regressor using EI with $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$8 (Rapp et al., 2023).

4. Quantum and photonic protocols

In continuous-variable quantum computation, a Gaussian Optimization Protocol is used not to optimize Gaussian states themselves, but to refine approximate non-Gaussian resources by task-specific Gaussian operations. The protocol applies only momentum displacement $s_i \sim \mathcal{N}(\mu_i,\sigma_i^2),$9 and squeezing $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$0 to approximate cubic phase resources, then optimizes $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$1 for either magic-state-based gate fidelity or nonlinear-quadrature variance. For MSBQC, the optimized cost is the gate fidelity

$$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$2

and for MBQC it is $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$3 computed from position-space expectation integrals. At $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$4, Gaussian preprocessing raises the gate fidelities of $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$5 and $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$6 to approximately $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$7, and reduces nonlinear-quadrature variance to $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$8 for $$s=\mu+\sigma\odot\epsilon,\qquad \epsilon\sim\mathcal{N}(0,I).$$9 and $P_i^{\mathrm{mod}}(t)$0 for $P_i^{\mathrm{mod}}(t)$1. A related genetic-algorithm state-preparation scheme based on finite Fock superpositions followed by Gaussian operations reaches gate fidelity $P_i^{\mathrm{mod}}(t)$2 and variance $P_i^{\mathrm{mod}}(t)$3 at $P_i^{\mathrm{mod}}(t)$4 (Jing et al., 19 Sep 2025).

A more hardware-centric variant appears in Gaussian boson sampling. There, the optimization problem is Max-Haf, and the Gaussian resource is not a surrogate model but the physical sampling distribution itself. For collision-free $P_i^{\mathrm{mod}}(t)$5-photon outputs $P_i^{\mathrm{mod}}(t)$6, a GBS device samples subsets with probability proportional to $P_i^{\mathrm{mod}}(t)$7, which realizes proportional sampling over candidate submatrices. The paper introduces GBS-Explore for seed generation and GBS-Tweak for local proposals, then plugs these into random search, simulated annealing, and greedy search. On a $P_i^{\mathrm{mod}}(t)$8 Haar-random instance with $P_i^{\mathrm{mod}}(t)$9, all three algorithms improved under GBS seeding, and GBS-enhanced random search performed best overall (Arrazola et al., 2018).

In broadband $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$0-type quantum memory, the Gaussian object is the temporal envelope. Signal and control are both restricted to Gaussian pulses,

$s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$1

and optimization is carried out over control pulse area $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$2, delay $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$3, and control duration $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$4. The protocol identifies resonant regimes by $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$5: an ATT region with $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$6, an ATS band $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$7 with $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$8, $s_i^{\mathrm{disc}}=\lfloor \mu_i+0.5\rfloor$9, $10\text{–}100\times$0, and an EIT regime with $10\text{–}100\times$1 and $10\text{–}100\times$2. For $10\text{–}100\times$3, the protocol-independent storage bound is $10\text{–}100\times$4, giving total efficiency $10\text{–}100\times$5, and the paper reports that Gaussian controls saturate this bound over a broad resonant region when $10\text{–}100\times$6 (Shinbrough et al., 2020).

These examples show that “Gaussian optimization” in physics often means optimization through Gaussian operations, Gaussian pulses, or Gaussian sampling, even when the target resource or objective is non-Gaussian.

5. Optimization on Gaussian-state manifolds and function spaces

A distinct line of work treats Gaussianity as the state manifold itself. In bosonic ground-state optimization, covariance matrices must satisfy the uncertainty principle $10\text{–}100\times$7. Symplectic Optimization on Gaussian States removes this constraint by parameterizing pure-state covariance matrices as positive-definite symplectic matrices through a unit-triangular factorization,

$10\text{–}100\times$8

with $10\text{–}100\times$9 and each $$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$0 symmetric. The variational problem becomes

$$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$1

which is globally unconstrained in the $$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$2. The method preserves physicality exactly, supports warm-start reuse across nearby Hamiltonians, and on weakly dipole-coupled lattices recovers ground-state energies, covariance matrices, and spectral gaps accurately. In the reported QDO lattice experiments, warm-starting from a nearby spacing reduced the number of optimization steps by about a factor of two (Willby et al., 28 Jan 2026).

A related but more geometric formulation works directly on the manifolds of pure bosonic and fermionic Gaussian states, $$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$3 and $$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$4. The optimization variables are group elements acting on a reference complex structure $$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$5, and the gradient is computed in a basis of generators orthogonal to the stabilizer. Because the metric components

$$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$6

are constant in this basis, the geometric gradient can be evaluated efficiently and updates can be taken by exponentials or Cayley transforms while preserving purity and any subgroup constraints. The paper develops this framework for variational ground states, Gaussian circuit complexity, and entanglement of purification, and reports numerical and analytical evidence that Gaussian purifications suffice for computing the entanglement of purification of arbitrary mixed Gaussian states (Windt et al., 2020).

A further generalization moves from finite-dimensional Gaussian states to Gaussian-process sample functions in Sobolev spaces. The consensus-based optimization algorithm in Sobolev spaces uses GP priors whose sample paths lie in $$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$7, conditions them on boundary, initial, or state constraints, and evolves an ensemble of feasible functions by a consensus-based update,

$$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$8

where the noise samples $$x_{\mathrm{cand}}^{(t)}=x^{(t)}+\sigma_t z_t,\qquad z_t\sim\mathcal{N}(0,I_d),$$9 come from a homogeneous posterior GP so that constraints remain preserved. The protocol is demonstrated on linear and nonlinear PDE boundary-value problems and on a nonlinear ODE-constrained shepherd control problem; in the latter, the best $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$0 run achieved nearly a 30% reduction relative to gradient descent, with mean reductions of about 27% (Khatab et al., 16 Mar 2026).

6. Limitations, misconceptions, and open directions

A common misconception is that “Gaussian Optimization Protocol” names a single method. The literature surveyed here suggests the opposite: the same label is used for distinct procedures in scheduling, black-box search, Bayesian optimization, photonic sampling, continuous-variable quantum computation, Gaussian-state geometry, and Sobolev-space control. What unifies them is methodological rather than algorithmic: Gaussian structure is used as the primary representation of uncertainty, search, or state space.

A second misconception is that Gaussian alignment with a nominal target automatically predicts downstream performance. In continuous-variable quantum computation, the task-oriented resource paper explicitly shows that state fidelity to the ideal cubic phase state is not predictive of task performance: $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$1 has the highest fidelity to $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$2 yet the worst MSBQC gate fidelity, while $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$3 can retain gate fidelity near $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$4 even when its state fidelity is nearly zero (Jing et al., 19 Sep 2025). This suggests that in Gaussian optimization protocols, the operative criterion is often task-conditioned utility rather than geometric closeness to a reference object.

The limitations are equally heterogeneous. GauS uses independent $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$5 variables and therefore ignores correlations among operators, which can slow convergence or miss global improvements; full Gaussian covariances would cost $f(x_{\mathrm{cand}}^{(t)})<f(x^{(t)})$6 (Cai et al., 23 Feb 2026). GCS has no covariance adaptation and may struggle in narrow curved valleys or ill-conditioned basins (Wong, 2023). GP-based methods inherit cubic GP training costs, kernel-hyperparameter instability, and acquisition sensitivity; WBGP-BO is explicitly motivated by criticalities of MLE fitting in Bayesian optimization (Candelieri et al., 18 May 2025). Quantum-kernel surrogates require Gram-matrix PSD correction and are exposed to shot noise and hardware noise (Rapp et al., 2023). Symplectic Gaussian-state optimization, as formulated in the cited work, targets pure states; mixed-state extensions lie outside that parameterization (Willby et al., 28 Jan 2026).

The strongest controversy concerns regret theory. “Gaussian Process Optimization with Mutual Information” claimed cumulative-regret bounds improving by an exponential factor over GP-UCB and introduced GP-MI with even stronger bounds, but the later erratum states that the main theorem does not hold under the standard Bayesian-optimization observation model because a martingale measurability argument fails. The authors further note that GP-MI can become overconfident and incur linear regret in some scenarios (Contal et al., 2013). This episode is significant because it shows that Gaussian uncertainty quantification does not by itself guarantee valid exploration theory.

Open directions in the cited literature remain substantial. GauS proposes mixtures of Gaussians, learned constraint surrogates, joint optimization with reinforcement learning or Bayesian optimization for tuning weights and initiation interval, and hybrid handoff to CP-SAT or ILP polishers (Cai et al., 23 Feb 2026). Task-oriented CV protocols point toward multimode Gaussian layers and broader Gaussian operation sets (Jing et al., 19 Sep 2025). Symplectic optimization identifies tensor-network enhancements and large-scale approximate treatments of weakly non-quadratic interactions as natural extensions (Willby et al., 28 Jan 2026). Taken together, these directions suggest that the future of Gaussian optimization is less about a single canonical protocol than about increasingly specialized Gaussian representations matched to the geometry and constraints of particular problem classes.