Neural Kernel Gaussian Process Regression

• NKGPR is a family of techniques that embeds neural networks into Gaussian process priors, enabling adaptive nonstationary kernel parameterization.
• It employs neural parameterization to modulate kernel hyperparameters, achieving flexible modeling of local variations and heteroscedastic noise.
• The method spans multiple formulations—from neural-GSM to NNGP/NTK—and has demonstrated competitive empirical performance on large, complex datasets.

Neural Kernel Gaussian Process Regression (NKGPR) denotes a family of Gaussian process regression methods in which neural networks determine some aspect of the GP prior, most often the covariance structure and, in some formulations, the mean as well. In the recent literature, the label encompasses at least four closely related constructions: neural parameterization of input-dependent kernel hyperparameters in nonstationary GPs, neural parameterization of generalized spectral mixture kernels, GP regression with kernels induced by infinite-width neural networks such as NNGP or NTK, and a matrix-parameterized cosmological variant that learns both mean and covariance by separate networks (James et al., 16 Jul 2025, Remes et al., 2018, Muyskens et al., 2024, Luo et al., 9 Aug 2025). This suggests that NKGPR is best understood as an umbrella term for GP models in which neural computation is embedded into kernel design or prior specification, rather than as a single standardized architecture.

1. Conceptual scope and taxonomy

The common structural idea is to preserve Gaussian-process inference while replacing some manually specified component of the prior by a learned neural map. In one line of work, the kernel family is fixed and the neural network outputs input-dependent hyperparameters such as local signal standard deviation, local noise standard deviation, or local lengthscale. In another, the neural network outputs the functions that define a nonstationary spectral kernel. In a third, the kernel itself is the GP limit of an infinitely wide network. A fourth variant, developed for cosmological reconstruction, uses one network for the mean function and another network for the covariance matrix directly (James et al., 16 Jul 2025, Remes et al., 2018, Muyskens et al., 2024, Luo et al., 9 Aug 2025).

Variant	Neuralized object	Representative formulation
Neural hyperparameter GP	Input-dependent kernel hyperparameters	$\tilde\theta(x)=g_\phi(x)$
Neural-GSM	$w_i(x), \mu_i(x), \ell_i(x)$ in GSM kernels	$k_{\text{GSM}}(x,x')$ with NN outputs
NNGP/NTK-based GPR	Kernel induced by infinite-width NN architecture	$k_{(M,\sigma_a,\sigma_b)}(x,x')$
Mean–kernel matrix NKGPR	Mean vector and covariance matrix	Mean Net + Kernel Net with $K=LL^\top$

A recurring distinction in this literature is between neural parameterization of kernel hyperparameters and deep kernel learning. In the former, the input space is not necessarily warped; instead, a known positive-semidefinite kernel construction is retained and its local parameters become functions of $x$. In deep kernel learning, by contrast, a neural feature map transforms inputs and a usually stationary kernel is applied in latent space. The 2025 nonstationary GP framework explicitly classifies its construction as a “neural-parameterized nonstationary kernel GP” rather than a deep kernel learning model (James et al., 16 Jul 2025).

2. Neural parameterization of nonstationary kernels

A canonical NKGPR formulation starts from GP regression with independent Gaussian noise,

$$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$

and a GP prior

$$f\sim\mathcal{GP}\big(\mu(\cdot),k(\cdot,\cdot\mid \theta,\tilde\theta(\cdot))\big).$$

Here $\theta$ denotes stationary kernel parameters and $\tilde\theta(x)$ denotes nonstationary parameters that vary over feature space. The neural parameterization sets

$w_i(x), \mu_i(x), \ell_i(x)$0

so that the network outputs quantities such as $w_i(x), \mu_i(x), \ell_i(x)$1, $w_i(x), \mu_i(x), \ell_i(x)$2, or $w_i(x), \mu_i(x), \ell_i(x)$3 (James et al., 16 Jul 2025).

The simplest example is nonstationary variance modulation:

$w_i(x), \mu_i(x), \ell_i(x)$4

or, in neural form,

$w_i(x), \mu_i(x), \ell_i(x)$5

with $w_i(x), \mu_i(x), \ell_i(x)$6 implemented as a feed-forward network. This preserves positive-definiteness because it is a multiplicative modulation of a positive-definite base kernel. The same framework extends to multiple nonstationary parameters; the paper gives a one-dimensional example based on a kernel from Heinonen et al. (2016) in which both $w_i(x), \mu_i(x), \ell_i(x)$7 and $w_i(x), \mu_i(x), \ell_i(x)$8 are neural outputs (James et al., 16 Jul 2025).

The Neural-GSM construction generalizes the same principle to nonstationary spectral kernels. Building on the generalized spectral mixture kernel of Remes et al. (2017), it replaces GP-parameterized functions by neural networks:

$w_i(x), \mu_i(x), \ell_i(x)$9

yielding

$k_{\text{GSM}}(x,x')$0

In this form, amplitude nonstationarity is carried by $k_{\text{GSM}}(x,x')$1, local smoothness by $k_{\text{GSM}}(x,x')$2 through the Gibbs kernel, and local frequency variation by $k_{\text{GSM}}(x,x')$3. Because the kernel formula remains that of a valid GSM construction, replacing latent GP parameter functions by deterministic neural functions does not destroy positive semidefiniteness, provided the outputs satisfy the required domain constraints such as $k_{\text{GSM}}(x,x')$4 (Remes et al., 2018).

Across both formulations, positivity constraints are enforced at the network output layer. The nonstationary GP framework explicitly suggests softplus or exponential activations, for example

$k_{\text{GSM}}(x,x')$5

while Neural-GSM uses softplus outputs for $k_{\text{GSM}}(x,x')$6, $k_{\text{GSM}}(x,x')$7, and $k_{\text{GSM}}(x,x')$8 (James et al., 16 Jul 2025, Remes et al., 2018).

3. Training, inference, and implementation

The dominant training pattern is joint optimization of GP parameters and neural parameters by differentiating through the kernel. For neural hyperparameter GPs, the prior can be written as

$k_{\text{GSM}}(x,x')$9

and the approximate log marginal likelihood is

$k_{(M,\sigma_a,\sigma_b)}(x,x')$0

with gradients propagated to both $k_{(M,\sigma_a,\sigma_b)}(x,x')$1 and $k_{(M,\sigma_a,\sigma_b)}(x,x')$2 through automatic differentiation. The paper describes this explicitly as chain-rule or backpropagation through the GP kernel and trains with Adam (James et al., 16 Jul 2025).

Scalability is handled by standard sparse or inducing-point approximations. The 2025 nonstationary GP framework uses inducing point / low-rank approximation with 100 inducing points on large UCI datasets, implemented in GPyTorch through custom kernels, InducingPointKernel, VariationalStrategy, PyTorch autograd, and Krylov solvers. Neural-GSM instead adopts sparse stochastic variational GPs as in Hensman et al. (2015), with inducing points $k_{(M,\sigma_a,\sigma_b)}(x,x')$3, a Gaussian variational posterior $k_{(M,\sigma_a,\sigma_b)}(x,x')$4, and joint ELBO optimization over $k_{(M,\sigma_a,\sigma_b)}(x,x')$5, $k_{(M,\sigma_a,\sigma_b)}(x,x')$6, $k_{(M,\sigma_a,\sigma_b)}(x,x')$7, neural-network parameters, and noise variance. In that implementation, GPflow/TensorFlow provides automatic differentiation, mini-batches of size 64 or 128 are used, and L2 weight regularization is applied to neural parameters (James et al., 16 Jul 2025, Remes et al., 2018).

The cosmological mean–kernel matrix formulation departs from analytic kernel families but retains closed-form GP conditioning. A Mean Net learns $k_{(M,\sigma_a,\sigma_b)}(x,x')$8, while a Kernel Net maps a feature vector of length 200 to the $k_{(M,\sigma_a,\sigma_b)}(x,x')$9 entries of a lower-triangular matrix $K=LL^\top$0, from which the covariance is constructed as

$K=LL^\top$1

To mitigate scale imbalance across indices, the model applies a symmetric PSD scaling matrix $K=LL^\top$2 with entries

$K=LL^\top$3

and uses the Hadamard-scaled covariance $K=LL^\top$4. Kernel Net is trained against an observed covariance matrix with the LogDet divergence,

$K=LL^\top$5

while Mean Net uses MSE. Training employs Optuna-based hyperparameter tuning, a dynamic learning-rate schedule, and early stopping (Luo et al., 9 Aug 2025).

Despite these differences, posterior prediction remains the standard GP conditional Gaussian calculation. In the cosmological formulation, for example,

$K=LL^\top$6

$K=LL^\top$7

with $K=LL^\top$8, $K=LL^\top$9, and $x$0 extracted from the network-generated covariance matrix (Luo et al., 9 Aug 2025).

4. Infinite-width neural kernels and spectral theory

A distinct branch of NKGPR identifies the kernel with the infinite-width limit of a neural architecture. For fully connected ReLU networks with Gaussian weights and biases, the NNGP kernel is generated by a layerwise covariance recursion. If the input-layer covariance is

$x$1

then higher layers satisfy

$x$2

and the kernel for depth $x$3 is

$x$4

This kernel enters GP regression exactly like any other covariance function (Muyskens et al., 2024).

The 2024 comparison with Matérn kernels emphasizes two practical constraints. First, normalization is necessary: the authors embed each input coordinate on a unit circle via trigonometric features and then $x$5-normalize, which places the data on a unit hypersphere and stabilizes the correlation recursion. Second, not all $x$6 choices yield numerically valid or well-conditioned covariance matrices. The valid region shrinks with depth, and in practice NNGP-based GPR behaves very similarly to a Matérn-$x$7 kernel under the studied conditions. The paper therefore recommends Matérn kernels over this particular neural kernel in practical applications, despite the theoretical interest of NNGP (Muyskens et al., 2024).

A complementary theoretical perspective comes from learning-curve analysis for GP regression and kernel ridge regression with power-law eigenspectra. If the kernel eigenvalues satisfy $x$8 and the target coefficients satisfy $x$9, then the generalization error of GPR with a neural kernel, including NNGP or NTK cases, decays as a power law in sample size. In the natural $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$0 regime summarized in the paper, the decay exponent can be written as

$$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$1

so that the error behaves like $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$2. This provides a spectral explanation for power-law learning curves in infinite-width neural models and frames NKGPR as ordinary GP or KRR with a neural kernel whose sample efficiency is controlled by the interaction of kernel spectral decay and target smoothness (Jin et al., 2021).

5. Empirical behavior and applications

The most extensive nonstationary-kernel benchmark in the supplied literature evaluates 12 UCI regression datasets, with up to 53,437 observations and 385 features, plus a spatial climate dataset with $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$3 locations. The compared models are a stationary GP, the proposed variance-only neural-parameterized model, the proposed variance-plus-noise model, and hierarchical baselines in which the nonstationary parameters themselves have GP priors. Using a Matérn $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$4 base kernel and 100 inducing points, the variance-plus-noise model achieves the best MSE on 9/12 datasets and the best log-score on 8/12 datasets; on the spatial dataset, the variance-only model attains the best log-score and recovers regions of high local variance near the poles and in some mid-latitude regions. The same experiments report that the neural models are competitive with or faster than hierarchical models, especially on larger datasets (James et al., 16 Jul 2025).

The Neural-GSM experiments span solar irradiance, sunspots, motion capture, power, and protein data. Neural-GSM achieves the best log likelihood and MAE/MSE on Power and Protein, the best log likelihood and MSE on Sunspots, and competitive performance on Irradiance; GP-GSM is best on Motion Capture, but Neural-GSM is reported as less variable across runs. Runtime comparisons show similar iteration counts to convergence for Neural-GSM and GP-GSM, while Neural-GSM is approximately twice as fast per iteration in wall-time because it avoids GP-based parameter interpolation (Remes et al., 2018).

The cosmological mean–kernel formulation applies NKGPR to reconstruct the apparent magnitude–redshift relation of Pantheon+ Type Ia supernovae and combine it with DESI DR2 baryon acoustic oscillation measurements in a test of the cosmic distance duality relation. Relative to a GaPP baseline with zero mean and RBF kernel, the NKGPR reconstruction is described as smoother, less oscillatory at high redshift, and equipped with uncertainty bands that widen more gradually. In the resulting marginalized likelihood analysis, no significant departure from the expected CDDR relation is found (Luo et al., 9 Aug 2025).

By contrast, the NNGP-versus-Matérn study argues for a negative practical result. On Friedman, MODIS, and Borehole benchmarks, a flexible Matérn family with tuned $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$5 and $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$6 attains substantially lower best-case RMSE than NNGP, while NNGP exhibits severe worst-case numerical instability for some hyperparameter settings. The central empirical conclusion is not that neural kernels are ineffective in general, but that this particular infinite-width fully connected ReLU kernel is numerically delicate and practically close to a Matérn-$$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$7 model (Muyskens et al., 2024).

6. Limitations, misconceptions, and research directions

A common misconception is that NKGPR necessarily means arbitrary “direct kernel learning.” The surveyed constructions do not support that interpretation. The successful models maintain positive semidefiniteness either by inserting neural outputs into an already valid kernel formula, as in multiplicatively modulated Matérn or GSM kernels, or by parameterizing the covariance as $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$8 and optionally applying PSD-preserving transformations such as a Hadamard product with a PSD scaling matrix (James et al., 16 Jul 2025, Remes et al., 2018, Luo et al., 9 Aug 2025).

A second source of confusion is the relation to adjacent GP–NN hybrids. Neural hyperparameter GPs are not the same as deep kernel learning, because they do not merely warp inputs and retain a stationary kernel in latent space. They are also not the same as deep GPs, which place GP transformations in multiple layers and typically require harder variational inference. NNGP- and NTK-based methods are yet another category: there the kernel is induced by the architecture itself, rather than by a neural network that outputs local kernel parameters (James et al., 16 Jul 2025, Remes et al., 2018, Muyskens et al., 2024).

The limitations reported across the literature are substantial. Neural parameterizations may require enough data, can be less interpretable than hand-designed basis functions or explicit covariates, and do not automatically provide uncertainty quantification over the nonstationary parameter functions themselves. Extrapolation far from training points can be unreliable because both the mean and the kernel parameter maps depend on neural-network behavior outside the observed region. Neural-GSM highlights non-convex optimization, overfitting risk, and sensitivity to the number of mixture components, architecture, inducing points, learning rate, and batch size. The cosmological formulation adds a fixed $$y_i=f(x_i)+\epsilon_i,\qquad \epsilon_i\sim\mathcal N(0,\tau_i^2),$$9 covariance dimension, $$f\sim\mathcal{GP}\big(\mu(\cdot),k(\cdot,\cdot\mid \theta,\tilde\theta(\cdot))\big).$$0 conditioning cost, possible systematics from high-redshift data augmentation, and an explicitly acknowledged ad hoc scaling matrix $$f\sim\mathcal{GP}\big(\mu(\cdot),k(\cdot,\cdot\mid \theta,\tilde\theta(\cdot))\big).$$1. The NNGP formulation requires hypersphere normalization and careful restriction of hyperparameters to avoid nearly constant or near-singular covariance matrices (James et al., 16 Jul 2025, Remes et al., 2018, Luo et al., 9 Aug 2025, Muyskens et al., 2024).

The proposed research directions are correspondingly diverse. For neural hyperparameter GPs, the literature points to CNNs for images, graph neural networks for network data, deeper feature extractors for structured nonstationarity, and neural parameterization of richer objects such as local lengthscales, full anisotropy matrices, and multi-output cross-covariances. Bayesian uncertainty over neural parameter functions, for example through priors on network weights, is identified as a natural extension. Neural-GSM suggests higher-dimensional and more structured kernels, while the cosmological model highlights low-rank or sparse approximations, more principled alternatives to the scaling matrix, and broader use in model-independent reconstruction problems. The learning-curve analysis adds a theoretical agenda: estimating the spectral exponent $$f\sim\mathcal{GP}\big(\mu(\cdot),k(\cdot,\cdot\mid \theta,\tilde\theta(\cdot))\big).$$2 of realistic neural kernels and the effective target exponent $$f\sim\mathcal{GP}\big(\mu(\cdot),k(\cdot,\cdot\mid \theta,\tilde\theta(\cdot))\big).$$3 in order to predict sample-efficiency regimes for NKGPR beyond toy settings (James et al., 16 Jul 2025, Remes et al., 2018, Luo et al., 9 Aug 2025, Jin et al., 2021).

Taken together, these works portray NKGPR not as a single model class but as a design principle: retain GP posterior calculus, and use neural networks to supply either local kernel parameters, spectral kernel functions, an induced infinite-width kernel, or even the mean and covariance of a discretized GP prior. The central technical challenge is to enlarge expressivity without abandoning PSD structure, numerical stability, or scalable inference; the central empirical lesson is that neuralization is most compelling when it targets specific forms of nonstationarity, heteroscedasticity, or structured covariance that are difficult to encode by stationary kernels alone.