Quantum Gaussian Processes Overview
- Quantum Gaussian Processes are frameworks that unify quantum state dynamics and machine learning by utilizing Gaussian channels and quantum kernels.
- They span continuous-variable Gaussian channels, quantum-assisted GP regression, and kernel methods, enabling efficient experimental protocols and robust inference.
- These models drive advancements in quantum state characterization, variational quantum algorithms, and scalable quantum inference for complex systems.
Quantum Gaussian Processes is a polysemous term spanning several research programs. In continuous-variable quantum information, it denotes Gaussian quantum processes or Gaussian channels: evolutions that map Gaussian states to Gaussian states and are naturally generated by Hamiltonians that are at most quadratic in the mode operators (Grove et al., 17 Mar 2025). In quantum machine learning, the same phrase is used for Gaussian process models whose covariance structure is induced by quantum circuits, for Gaussian process inference or training accelerated by quantum linear-algebra subroutines, and, more recently, for Bayesian learning frameworks defined directly on quantum states and quantum transformations (Rapp et al., 2023, Jäger et al., 30 Apr 2026). A further, related usage appears in many-body variational physics, where the “quantum Gaussian process state” is a kernel-inspired wavefunction ansatz built from quantum support states rather than classical support configurations (Rath et al., 2021).
1. Scope and nomenclature
The term covers distinct objects that share Gaussian-process language but differ in ontology, data model, and mathematics.
| Usage | Core object | Representative papers |
|---|---|---|
| Continuous-variable Gaussian process | Gaussian channel on bosonic modes | (Grove et al., 17 Mar 2025, Mehboudi et al., 2019, Braccini et al., 1 Oct 2025) |
| Quantum-assisted GP regression | Quantum algorithm for GP linear algebra or training | (Zhao et al., 2015, Zhao et al., 2018, Hu et al., 22 Mar 2025) |
| Quantum-kernel GP | Classical GP with covariance from a quantum circuit | (Rapp et al., 2023, Dai et al., 2022, Gandhi et al., 16 Feb 2026) |
| QGP for quantum data | GP prior over unknown quantum transformations | (Jäger et al., 30 Apr 2026) |
| GP limit of QNNs | Random or trained QNN outputs converge to a GP | (García-Martín et al., 2023, Girardi et al., 2024) |
| Quantum Gaussian process state | Kernel-inspired variational wavefunction ansatz | (Rath et al., 2021) |
A recurrent source of confusion is that only the first usage refers directly to physical Gaussian quantum dynamics in the standard continuous-variable sense. In the kernel and linear-solver literature, the GP itself usually remains the standard classical Gaussian process, while the quantum component enters through the kernel, the matrix inversion, or the hyperparameter optimization pipeline. This suggests that “quantum Gaussian process” is best treated as a family resemblance term rather than a single formalism.
2. Continuous-variable Gaussian processes as physical quantum dynamics
In an -mode bosonic system with annihilation and creation operators , the quadratures are defined by
and collected into
A Gaussian state is completely specified by its first and second moments,
A Gaussian quantum process maps Gaussian states to Gaussian states, and at the level of mean values is represented by
where is a real symplectic matrix satisfying and (Grove et al., 17 Mar 2025).
This class includes passive transformations such as linear optics and active transformations such as squeezing, parametric amplification, and other quadratic nonlinear optics. For passive linear optics, the same physical transformation can be written in terms of an 0 unitary 1, related to the symplectic matrix by
2
The symplectic description is more general, because it covers arbitrary Gaussian transformations rather than only number-preserving ones (Grove et al., 17 Mar 2025).
A 2025 characterization protocol reconstructs this symplectic matrix directly using only Gaussian resources: coherent-state probes and quadrature measurements. The basic procedure is to inject a coherent state into one input mode at a time, keep the remaining inputs in vacuum, and infer columns of 3 from output first moments. By choosing the probe phase to make either 4 or 5 vanish, one isolates one column of 6. Repeating this over all input-output mode pairs reconstructs all 7 real entries. The method requires 8 settings, is resilient to uniform loss, and in simulation heterodyne measurements outperform homodyne measurements for general Gaussian process reconstruction, while both are comparable for purely passive linear-optical transformations. The simulations were performed in Strawberry Fields for random symplectic matrices with 9 to 0, with either 1 or 2 uniform loss, repeated 500 times for each 3; the reconstruction error was reported as essentially independent of 4, and the method remained effective even with 5 uniform loss (Grove et al., 17 Mar 2025).
The same work models uniform loss by fictitious beam splitters of transmissivity 6, giving an experimentally reconstructed matrix 7. Since 8, the transmissivity is inferred directly from 9, so no extra calibration step is needed. This feature is significant because it ties experimental accessibility to the intrinsic symplectic structure, rather than to auxiliary non-Gaussian probes or separate detector characterization (Grove et al., 17 Mar 2025).
3. Response theory and hybrid superpositions of Gaussian dynamics
Beyond static characterization, Gaussian quantum processes also support a dynamical response theory formulated directly at the level of covariance matrices. For an 0-mode bosonic system with quadratures 1, a general Gaussian channel acts as
2
In the zero-mean setting 3, the channel is fully characterized by 4, with complete positivity condition
5
For a one-parameter family of Gaussian channels 6, the stationary covariance matrix satisfies
7
A linear-response theory for such systems gives, in continuous time,
8
where 9 is the static response (Mehboudi et al., 2019).
The same framework yields a fluctuation-dissipation theorem for Gaussian covariance matrices,
0
with 1 the symmetric logarithmic derivative. For thermal equilibrium under a quadratic Hamiltonian, the response reduces to a Gaussian Kubo-type formula,
2
and for quadratic Hamiltonians with linear Lindblad operators the covariance obeys
3
with steady state determined by the Lyapunov equation
4
This places non-equilibrium steady states, time-dependent Gaussian channels, and time-dependent Lindbladian master equations in a single Gaussian response formalism (Mehboudi et al., 2019).
A separate generalization extends Gaussian dynamics to hybrid continuous-variable–qubit or qudit systems, producing what are called Gaussian-Branched Cat States. For a single qubit,
5
and each branch is assigned a branched characteristic function
6
The state is then fully characterized by generalized covariance matrices, generalized first moments, and the qubit reduced density matrix. The formalism treats unitary and open dynamics exactly, without truncation, and includes ideal and noisy measurements. In this setting the system no longer follows a single Gaussian branch; instead, discrete degrees of freedom control a superposition of Gaussian processes, with interference encoded in generally complex off-diagonal branch data (Braccini et al., 1 Oct 2025).
This hybrid extension is structurally important because it preserves the calculational advantages of Gaussian phase-space methods while moving beyond globally Gaussian states. The reported examples—a squeezed, leaking, and measured resonator entangling two qubits, and a levitated nanoparticle undergoing Stern-Gerlach interferometry in a diffusive environment—show that Gaussian-process language in quantum physics can refer not only to Gaussian channels but also to controlled superpositions of Gaussian evolutions (Braccini et al., 1 Oct 2025).
4. Quantum algorithms for Gaussian process inference, training, and quadrature
In supervised learning, a standard GP regression model with covariance matrix 7 and noise variance 8 uses
9
Classically, exact inference or training is dominated by matrix factorization or inversion with 0 cost. The earliest quantum-assisted formulation applies the quantum linear systems algorithm to these two scalar quantities rather than to recovery of the full inverse, so the output model is already tailored to GP prediction (Zhao et al., 2015).
That line of work makes the usual conditional claims. Exponential improvement requires the regularized kernel matrix to be sparse and well conditioned, together with efficient QRAM-style state preparation of the relevant vectors; otherwise the advantage degrades to polynomial. The same access assumptions recur across later proposals. This is a central caveat rather than a secondary detail, because the complexity gain depends more on matrix structure and data access than on GP regression alone (Zhao et al., 2015).
Training by log marginal likelihood (LML) was developed further by decomposing
1
A 2018 algorithm estimates 2 by eigenvalue sampling via phase estimation and computes 3 using a modified quantum linear system routine based on 4. Under sparse access, the single-run complexity is
5
and the paper states an exponential speedup over classical 6 GP training in favorable sparse settings, with polynomial improvement in denser regimes (Zhao et al., 2018).
A later proposal moves from derivative-free LML evaluation to direct quantum optimization of GP hyperparameters. The LML gradient
7
is rewritten as the real part of an inner product, then estimated by block-encoding, quantum matrix inversion, and amplitude-estimation-style subroutines. The paper presents both a hybrid quantum-classical gradient descent and a fully quantum version in which the parameter register is updated in circuit. Its sparse-access runtime claim is
8
together with an exponential speedup claim for sparse and well-conditioned kernels and a cubic speedup claim in the dense worst case 9 (Hu et al., 22 Mar 2025).
Other algorithmic programs replace direct inversion by low-rank approximations. Quantum-Assisted Hilbert-Space Gaussian Process Regression rewrites the posterior mean and variance in terms of a basis-expansion matrix 0, then combines qPCA, conditional rotations, Hadamard tests, and Swap tests. Its reported total complexity is
1
with a stated polynomial computational complexity reduction over the classical Hilbert-space approximation route (Farooq et al., 2024). A closely related Gaussian process quadrature construction uses a Hilbert-space kernel approximation, qPCA, QPE, Hadamard tests, and SWAP tests to estimate the Bayesian quadrature mean and variance, with claimed complexity
2
and a stated polynomial advantage over classical Gaussian process quadrature methods (Galvis-Florez et al., 20 Feb 2025).
Near-term variants trade asymptotic claims for NISQ compatibility. A 2025 hybrid solver reformulates GP posterior computation as a set of linear systems solved by the Variational Quantum Linear Solver. For a one-dimensional Snelson test function with 3, the combination of a tapered Matérn 4 kernel and the MUHEA ansatz yielded classical GP MT MSE 5 and VQLS-GP MT MSE 6; the same paper reports a reduction from 41 Pauli strings for the RBF kernel to 23 for the MT kernel, about a 44% reduction (Bükrü et al., 17 Oct 2025). A separate multi-output GP for electrical-grid line-parameter estimation used HHL together with Approximate Quantum Compiling to implement up to a 13-qubit HHL circuit for a 7 kernel matrix on IBM Auckland; the result was presented as a hardware proof of concept rather than a competitive classical replacement (Ganeshamurthy et al., 2024).
5. Quantum kernels, Bayesian optimization, and hybrid Gaussian process models
In quantum-kernel GP models, the GP posterior equations remain classical, but the covariance function is computed from quantum states. A standard construction defines a feature map
8
with kernel
9
or, for pure states,
0
Once the Gram matrix is estimated, the predictive mean and covariance are the usual GP quantities (Rapp et al., 2023).
A central issue in this literature is preservation of posterior variance. In one quantum Gaussian process regression framework for Bayesian optimization, the Gram matrix estimated from finite shots and hardware noise may fail to be positive definite, so the paper uses eigenvalue-cutoff regularization: negative eigenvalues are truncated to zero and the matrix is reconstructed from the truncated spectrum and original eigenvectors. The stated reason is to preserve the variance structure as much as possible, since Bayesian optimization depends directly on the surrogate variance through expected improvement (Rapp et al., 2023).
The same study benchmarked one-dimensional regression on 1 using a hardware-efficient feature map with 2 qubits and 3 layers. It reported, for ideal statevector simulation, 4 and MSE 5; for sample-based simulation with 10,000 measurements, 6 and MSE 7; and on ibmq_montreal, 8 and MSE 9. The same model was used as a surrogate for Bayesian optimization on the Branin-Hoo function and for hyperparameter optimization of a gradient boosting model, with the stated result that the quantum Bayesian optimization version can match the performance of its classical counterpart (Rapp et al., 2023).
Quantum kernels have also been used for molecular potential energy surfaces. For H0O1, a six-qubit circuit assigns one qubit to one input dimension and defines the kernel by state overlap. The paper compares entangled and unentangled quantum kernels, learns the circuit parameters by maximizing a stabilized objective
2
and uses Bayesian optimization initialized with 20 random points, 3, with convergence typically in roughly 30–100 iterations. On a 6D global PES with 31,124 ab initio points, a representative model using 1000 training points, an entangled quantum kernel, and 72 BO iterations achieved an RMSE of about 4; with larger training sizes, the entangled quantum GP and the optimized classical RBF GP both approached RMSE values around 5. The paper further reports that entanglement improves extrapolation and that, for low training-energy thresholds, the entangled quantum GP can outperform the classical RBF GP (Dai et al., 2022).
A distributed variant extends quantum-kernel GPs to multi-agent settings. In that framework, local agents train quantum-kernel models and aggregate parameters via a Distributed consensus Riemannian Alternating Direction Method of Multipliers on the torus manifold. The quantum kernel may be a fidelity kernel or, primarily, a projected quantum kernel built from quantum expectation values followed by a classical outer kernel. On NASA Shuttle Radar Topography Mission elevation data and synthetic data generated by Quantum Gaussian Processes, the method reports, on average over SRTM datasets, 6 lower NRMSE than FACT-GP, 7 lower NRMSE than apxGP, and about 8 improvement over Full-GP in NRMSE; the paper also notes that DQGP can become somewhat overconfident or less calibrated in some larger-network settings (Gandhi et al., 16 Feb 2026).
Across these studies, quantum-kernel GPs are generally presented as hybrid models rather than fully quantum predictors. A common misconception is that the “quantum” label implies a uniformly superior surrogate. The published results are more restrained: some benchmarks show parity with strong classical baselines, some show improvements in specific extrapolation or non-stationary settings, and uncertainty calibration remains an explicit design constraint rather than an automatic consequence of quantum feature spaces (Rapp et al., 2023, Dai et al., 2022).
6. Quantum data, Gaussian-process limits of quantum models, and related state ansätze
A more recent line of work defines quantum Gaussian processes directly on quantum systems. In this framework, the unknown object is a quantum transformation 9, the data are quantum states 0, and the stochastic process is
1
A quantum Gaussian process arises when the induced random variables are jointly Gaussian under a suitable prior on 2. For matchgate, or free-fermionic, evolutions, the covariance kernel is
3
and the paper identifies this as the first scalable family in which the unknown unitary acts non-trivially on all qubits. It further reports a direct Bell-basis estimation protocol with sample complexity 4, compared with 5 for indirect coefficient-wise estimation, and gives demonstrations of accurate long-range extrapolation, phase-diagram learning, about 6 classification accuracy in a bond-alternating XXX-model phase task, and Bayesian optimization in a quantum sensing problem that found the maximally sensitive probe state in only 7 evaluations (Jäger et al., 30 Apr 2026).
Another notion of QGP emerges from random and trained quantum neural networks. For deep Haar-random QNNs, the outputs
7
converge in distribution to Gaussian processes in the large Hilbert-space limit 8, with covariance asymptotically proportional to the overlap kernel,
9
The same paper shows that, because these covariances are only 00, GP prediction becomes trivial under only polylogarithmically many shots: the predictive distribution reduces to the prior 01. It also proves stronger concentration for expectation values and gradients than earlier second-moment arguments, with tails of order 02 (García-Martín et al., 2023).
A complementary infinite-width result concerns parameterized quantum circuits with local readout. There the untrained network output converges to a GP whenever each measured qubit is correlated only with few others and
03
Under gradient descent with square loss, and as long as the model is not affected by barren plateaus, the trained network fits the training set, remains in the lazy regime, and still converges in distribution to a GP. The paper further proves that a polynomial number of measurements is sufficient for these statements to survive sampling noise (Girardi et al., 2024).
The “quantum Gaussian process state” is related in name but distinct in purpose. It is a variational many-body wavefunction ansatz defined by
04
where the support points are unentangled product states rather than discrete classical configurations. The ansatz is presented as systemically improvable to exactness as 05 is increased, can be viewed as an exponentiated matrix product state with diagonal, commuting site matrices, and admits a deterministic Bayesian sweep algorithm based on closed-form Gaussian posteriors over site weights. On frustrated 06-07 systems, the paper reports improved generalization from sparse training data, including learning from only 08 or 09 of the full Hilbert space of amplitudes for a 10 ground state (Rath et al., 2021).
Taken together, these results make clear that the phrase “quantum Gaussian process” now spans three technical levels: Gaussian quantum dynamics in continuous-variable physics, Gaussian-process machinery augmented by quantum kernels or quantum linear algebra, and genuinely quantum Bayesian models in which the prior is over unknown quantum transformations or quantum network functions. A plausible implication is that future work will continue to separate these meanings more sharply, because each has different scalability claims, different access assumptions, and different notions of what counts as “quantum” in the model.