Paper deep dive

From Shallow Bayesian Neural Networks to Gaussian Processes: General Convergence, Identifiability and Scalable Inference

Gracielle Antunes de Araújo, Flávio B. Gonçalves

Year: 2026Venue: arXiv preprintArea: stat.MLType: PreprintEmbeddings: 100

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%

Last extracted: 7/20/2026, 10:35:25 AM

Summary

This paper establishes a general convergence result from shallow Bayesian Neural Networks (BNNs) to Gaussian Processes (GPs) under relaxed assumptions. It introduces a novel covariance function constructed as a convex mixture of components induced by four activation functions (tanh, sigmoid, ReLU, Leaky-ReLU). The study characterizes the positive definiteness and identifiability of this mixed kernel. Furthermore, it proposes a scalable Maximum A Posteriori (MAP) inference procedure using Nyström approximation to manage the computational cost-accuracy trade-off, demonstrating effectiveness on tabular datasets.

Entities (10)

Gaussian Processes · model · 98%Shallow Bayesian Neural Networks · model · 98%Maximum A Posteriori · method · 95%Mixed Covariance Function · method · 95%Nyström approximation · algorithm · 95%Leaky-ReLU · activationfunction · 92%ReLU · activationfunction · 92%Sigmoid · activationfunction · 92%Tanh · activationfunction · 92%Identifiability · property · 90%

Relation Signals (7)

Shallow Bayesian Neural Networks → convergesto → Gaussian Processes

confidence 95% · We first establish a general convergence result from BNNs to GPs by relaxing assumptions used in prior formulations

Nyström approximation → enables → Scalable MAP Training

confidence 92% · develop a scalable maximum a posterior (MAP) training and prediction procedure using a Nyström approximation

Mixed Covariance Function → composedof → Sigmoid

confidence 90% · proposed a new covariance function defined as a convex mixture of components induced by four widely used activation functions... sigmoid

Mixed Covariance Function → composedof → Leaky-ReLU

confidence 90% · proposed a new covariance function defined as a convex mixture of components induced by four widely used activation functions... Leaky-ReLU

Mixed Covariance Function → composedof → ReLU

confidence 90% · proposed a new covariance function defined as a convex mixture of components induced by four widely used activation functions... ReLU

Mixed Covariance Function → composedof → Tanh

confidence 90% · proposed a new covariance function defined as a convex mixture of components induced by four widely used activation functions... tanh

Mixed Covariance Function → hasproperty → Identifiability

confidence 88% · characterize key properties including positive definiteness and both strict and practical identifiability under different input designs

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:In this work, we study scaling limits of shallow Bayesian neural networks (BNNs) via their connection to Gaussian processes (GPs), with an emphasis on statistical modeling, identifiability, and scalable inference. We first establish a general convergence result from BNNs to GPs by relaxing assumptions used in prior formulations, and we compare alternative parameterizations of the limiting GP model. Building on this theory, we propose a new covariance function defined as a convex mixture of components induced by four widely used activation functions, and we characterize key properties including positive definiteness and both strict and practical identifiability under different input designs. For computation, we develop a scalable maximum a posterior (MAP) training and prediction procedure using a Nyström approximation, and we show how the Nyström rank and anchor selection control the cost-accuracy trade-off. Experiments on controlled simulations and real-world tabular datasets demonstrate stable hyperparameter estimates and competitive predictive performance at realistic computational cost.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

98,981 characters extracted from source content.

Expand or collapse full text

From Shallow Bayesian Neural Networks to Gaussian Processes: General Convergence, Identifiability and Scalable Inference Gracielle Antunes de Araújo Universidade do Estado de Minas Gerais (UEMG), FaEnge – Faculdade de Engenharia. gracielle.araujo@uemg.br Flávio B. Gonçalves Universidade Federal de Minas Gerais (UFMG), ICEx – Instituto de Ciências Exatas. fbgoncalves@ufmg.br Abstract In this work, we study scaling limits of shallow Bayesian neural networks (BNNs) via their connection to Gaussian processes (GPs), with an emphasis on statistical modeling, identifiability, and scalable inference. We first establish a general convergence result from BNNs to GPs by relaxing assumptions used in prior formulations, and we compare alternative parameterizations of the limiting GP model. Building on this theory, we propose a new covariance function defined as a convex mixture of components induced by four widely used activation functions, and we characterize key properties including positive definiteness and both strict and practical identifiability under different input designs. For computation, we develop a scalable maximum a posterior (MAP) training and prediction procedure using a Nyström approximation, and we show how the Nyström rank and anchor selection control the cost-accuracy trade-off. Experiments on controlled simulations and real-world tabular datasets demonstrate stable hyperparameter estimates and competitive predictive performance at realistic computational cost. 1 Introduction Mobile devices and connected systems continuously generate large volumes of complex and weakly structured data, fostering the use of machine learning models in prediction and decision-making tasks. In this setting, neural networks (NNs) stand out for their strong empirical performance, yet their increasing complexity often reduces interpretability and hampers rigorous uncertainty quantification—an issue that becomes critical in applications requiring calibration, robustness, and reliability (????). A natural way to incorporate uncertainty is Bayesian modeling, which treats weights and biases as random variables and yields predictive distributions equipped with explicit measures of uncertainty. For shallow BNNs, classical results show that, under standard initialization and scaling assumptions, the infinite-width regime converges to a GP whose kernel is determined by the activation function and the prior distributions over the parameters (????). This connection provides a useful statistical bridge between N models and stochastic processes, but it leaves open practical choices in modeling and inference when one seeks both kernel flexibility and computational tractability at scale. From a modeling perspective, much of the literature still focuses on isolated activations or restricted families, which limits the diversity of limiting covariance structures and may constrain the ability to capture, in an interpretable way, combined effects of smoothness and angular behavior commonly observed in real data (????). From a computational perspective, even with an expressive kernel, exact GP inference scales cubically in time and quadratically in memory with the sample size, which makes direct use prohibitive for large datasets (??). Thus, exploiting the network–process equivalence brings a recurring tension between kernel expressiveness, interpretability/uncertainty quantification, and computational cost. This paper studies scaling limits of shallow BNNs via their connection to GPs, with emphasis on modeling, identifiability, and scalable inference. First, we establish a general convergence result from wide BNNs to GPs under relaxed assumptions. Second, we introduce a new covariance function based on a convex mixture of activation-induced components (tanh, sigmoid, ReLU, and Leaky-ReLU), yielding an interpretable kernel that combines smooth and angular behavior. Third, we study positive definiteness and both strict and practical identifiability as a function of the input design, highlighting regimes where different mixtures become difficult to distinguish from finite samples. Finally, we develop a scalable MAP training and prediction procedure via a Nyström approximation, making the cost–accuracy trade-off explicit through the rank and anchor-selection strategy. We evaluate the proposed methodology on simulations and real-world tabular datasets, emphasizing estimation stability, predictive performance, and scalability as a function of the Nyström rank and anchor selection. The remainder of the article is organized as follows. Section 1.1 reviews the connection between shallow BNNs and GPs and introduces notation. Section 2 establishes the infinite-width limit and the general convergence result. Section 3 derives the mixed kernel, discusses its properties, and presents identifiability results. Section 4 describes the Nyström-based MAP procedure for scalable training and prediction. Section 5 presents the experiments and discusses the empirical findings. Finally, Section 6 summarizes the conclusions and outlines directions for future work. 1.1 Background and Related Work Neural networks provide flexible function approximation and strong predictive performance, but increasing model complexity can hinder interpretability and rigorous uncertainty quantification. A Bayesian formulation addresses these issues by placing priors on weights and biases and producing predictive distributions with explicit measures of uncertainty (????). This perspective is particularly relevant in applied machine learning settings where decisions must be accompanied by calibrated notions of reliability. A central theoretical result is that shallow BNNs, under standard scaling and independence assumptions, converge in the infinite-width limit (in the sense of finite-dimensional distributions) to Gaussian processes (??). For specific architectures, closed-form expressions for the induced covariance function can be derived, completing the description of the limiting GP model (??). The broader theory and exact inference machinery for GPs are consolidated in (??), providing a unified statistical framework for probabilistic prediction. For deeper architectures, this correspondence extends through layer-wise recursions and is often described through the Neural Network Gaussian Process (NNGP∞NNGP_∞) and the Neural Tangent Kernel (NTK), which clarify how activation, scaling, and initialization shape inductive bias and regularity, and how training dynamics can be characterized in large-width regimes (??????????). From a modeling standpoint, however, much of the literature focuses on isolated activation choices or restricted kernel families, limiting the range of covariance structures available for capturing combined smooth and angular behaviors that can arise in practice (????). Moreover, practical deployment requires scalable inference: exact GP computations scale cubically in time and quadratically in memory with the sample size, motivating scalable approximations such as low-rank methods and sparse neighbor-based constructions (????????????). Two additional considerations become important when moving from limit theory to practice. First, finite-width networks exhibit residual dependencies that can affect effective kernels and uncertainty calibration (??????). Second, as one adopts richer kernel parameterizations, identifiability and numerical stability become central, since distinct components may be difficult to distinguish under realistic input designs and finite samples (????). These factors motivate methodological development that jointly addresses kernel expressiveness, identifiability, and scalable inference. In this work, we build on these foundations by establishing a general convergence result from shallow Bayesian neural networks to Gaussian processes under weaker assumptions (????), proposing an interpretable mixed covariance function as a convex combination of activation-induced components with theoretical guarantees, and developing a scalable MAP training and prediction pipeline based on a Nyström approximation, where rank and anchor selection transparently govern the cost–accuracy trade-off. 2 Infinite-width limit for shallow BNN The Bayesian approach to feedforward neural networks assigns prior distributions to model parameters, thereby inducing a prior distribution over functions. A classical result shows that, when the number of neurons in one (or more) hidden layers tends to infinity (H→∞H→∞), the induced prior over functions converges to a GP, whose covariance is jointly determined by the activation function and the parameter priors (??????). 2.1 Model: shallow BNN with one hidden layer Architecture 1 (Feedforward network with one hidden layer). Let ⊆ℝIX ^I and ⊆ℝKY ^K. Fix H∈ℕH (number of hidden units) and activation functions h:ℝ→ℝh:R (hidden layer) and g:ℝ→ℝg:R (output layer), both applied componentwise. We define the network :→f:X with parameters ∈ℝH×I,∈ℝH,∈ℝK×H,∈ℝK.U ^H× I,\ a ^H,\ V ^K× H,\ b ^K. Given an input ∈x , the forward pass is =+∈ℝH,=h()∈ℝH,()=g(+)∈ℝK,z=Ux+a ^H, =h(z) ^H, (x)=g\! (Vh+b ) ^K, where g(⋅)g(·) is applied elementwise in ℝKR^K. In indexed form, for k=1,…,Kk=1,…,K, fk()=g(bk+∑j=1Hvkjh(aj+∑i=1Iujixi)).f_k(x)=g\! (b_k+ _j=1^Hv_kj\,h\! (a_j+ _i=1^Iu_jix_i ) ). For regression, we take g(⋅)≡Id(⋅)g(·) (·). Hence, ()=+∈ℝK,fk()=bk+∑j=1Hvkjh(aj+∑i=1Iujixi).f(x)=Vh+b ^K, f_k(x)=b_k+ _j=1^Hv_kj\,h\! (a_j+ _i=1^Iu_jix_i ). (1) This architecture is a fully-connected shallow multilayer perceptron (MLP), i.e., a feedforward network with a single hidden layer. Figures 1 and 2 provide two complementary views of the same one-hidden-layer feedforward network, as defined in Architecture 1. In Figure 1, an input =(x1,…,xI)x=(x_1,…,x_I) feeds all H hidden units simultaneously: each neuron j computes the pre-activation zj=aj+∑i=1Iujixiz_j=a_j+ _i=1^Iu_jix_i and applies the nonlinearity h(⋅)h(·), yielding hj=h(zj)h_j=h(z_j). The hidden activations are then linearly combined at the output layer: for each output component k, the model sums the contributions vkjhjv_kjh_j and adds the bias bkb_k, resulting in fk=bk+∑j=1Hvkjhjf_k=b_k+ _j=1^Hv_kjh_j (in regression, with g(⋅)≡Idg(·) ). Figure 2 presents the same computation from a dataset perspective: one selects a column (p)x^(p) from the input matrix X, propagates this observation through the network to obtain ((p))∈ℝKf(x^(p)) ^K, and by repeating this for p=1,…,np=1,…,n, stacks the outputs to form =[((1))⋯((n))]∈ℝK×nF=[f(x^(1))\ ·s\ f(x^(n))] ^K× n. Figure 1: Shallow feedforward network (one-hidden-layer MLP): I inputs, H hidden units, and K outputs. The figure illustrates the regression case, where g(⋅)≡Id(⋅)g(·) (·).1x_1ix_iIx_I⋮ ⋮ 1h_1hjh_jhHh_H⋮ ⋮ z1=a1+∑i=1Iu1ixiz_1=a_1+ _i=1^Iu_1ix_i h1=h(z1)h_1=h(z_1) zj=aj+∑i=1Iujixiz_j=a_j+ _i=1^Iu_jix_i hj=h(zj)h_j=h(z_j) zH=aH+∑i=1IuHixiz_H=a_H+ _i=1^Iu_Hix_i hH=h(zH)h_H=h(z_H) f1f_1fkf_kfKf_K⋮ ⋮ f1=b1+∑j=1Hv1jhjf_1=b_1+ _j=1^Hv_1j\,h_j fk=bk+∑j=1Hvkjhjf_k=b_k+ _j=1^Hv_kj\,h_j fK=bK+∑j=1HvKjhjf_K=b_K+ _j=1^Hv_Kj\,h_j Figure 2: The column (p)∈ℝIx^(p) ^I is propagated through the network, producing the output vector ((p))∈ℝKf(x^(p)) ^K. Repeating this for p=1,…,np=1,…,n and stacking the outputs as columns yields =[((1))⋯((n))]∈ℝK×nF=[f(x^(1))\ ·s\ f(x^(n))] ^K× n. =[(1)⋯(p)⋯(n)]∈ℝI×nX= [\,x^(1)\ ·s\ x^(p)\ ·s\ x^(n)\, ] ^I× nChoose a column: (p)∈ℝIx^(p) ^I(p)x^(p)x1(p)x_1^(p)⋮ (p)x_i^(p)⋮ (p)x_I^(p) (p)x^(p)h1h_1⋮ _j⋮ _Hf1(p)f_1^(p)⋮ (p)f_k^(p)⋮ (p)f_K^(p)((p))=[f1(p)⋮fk(p)⋮fK(p)]∈ℝKf(x^(p))= bmatrixf_1^(p)\\ \\ f_k^(p)\\ \\ f_K^(p) bmatrix ^K=[((1))⋯((p))⋯((n))]∈ℝK×nF= [\,f(x^(1))\ ·s\ f(x^(p))\ ·s\ f(x^(n))\, ] ^K× n ((p))f(x^(p))Repeat for p=1,…,np=1,…,n and stack the resulting outputs.F results from stacking the output vectors column by column We assume that the observed responses (p)y^(p) are noisy measurements of the network outputs. In particular, we consider a nonlinear regression model with i.i.d. Gaussian observation noise: (p)=((p))+(p),(p)∼(,σϵ2K),p=1,…,n,y^(p)\;=\;f(x^(p))\;+\; ^(p), ^(p) (0,\, _ε^2I_K), p=1,…,n, (2) where, for k=1,…,Kk=1,…,K, fk()=bk+∑j=1Hvkjh(aj+∑i=1Iujixi).f_k(x)\;=\;b_k\;+\; _j=1^Hv_kj\,h\! (a_j\;+\; _i=1^Iu_ji\,x_i ). (3) Under the Bayesian approach, we place priors on the network parameters. Specifically, ujij,i∼i.i.d.(0,σu2),ajj=1H∼i.i.d.(0,σa2),vkjk,j∼i.i.d.(0,σv2H),bkk=1K∼i.i.d.(0,σb2).\u_ji\_j,i .i.d. N(0, _u^2),\;\;\;\a_j\_j=1^H .i.d. N(0, _a^2),\;\;\;\v_kj\_k,j .i.d. N\! (0, _v^2H ),\;\;\;\b_k\_k=1^K .i.d. N(0, _b^2). (4) These priors contribute to inference by regularizing the likelihood implied by the Gaussian nonlinear regression model. In particular, the choice vkj∼(0,σv2H)v_kj \! (0, _v^2H ) ensures that the output variance remains finite as the hidden-layer width H grows. Without the 1/H1/H scaling, each additional hidden unit would add variance to the output, and the variance of fk()f_k(x) would typically diverge as H→∞H→∞. Architecture 2 (Bayesian feedforward neural network (BFNN)). Let ⊆ℝIX ^I and ⊆ℝKY ^K. Fix H∈ℕH (number of hidden units) and a componentwise activation h:ℝ→ℝh:R . The observation model and the prior specification are given by (2), (3), and (4). Under mild moment conditions (see Section˜2.3), finite-dimensional distributions of the network outputs converge to Gaussian distributions, thus characterizing the limit as a GP. xh1h_1⋮ _j⋮ _HH→∞H→∞fkf_kOne-hidden-layer network with output k(fk(1),…,fk(n))⊤→H→∞(0,K(,)),fk(⋅)∼(0,K). aligned (f_k^(1),…,f_k^(n) ) & [d]H→∞N (0,\,K(X,X) ),\\[2.77501pt] f_k(·)& (0,\,K ). alignedGP associated with output k Figure 3: Limit of a wide one-hidden-layer network to a GP for output k. As H→∞H→∞, the vector (fk(1),…,fk(n))⊤ (f_k^(1),…,f_k^(n) ) converges in distribution to a multivariate normal with covariance K(,)K(X,X), and the collection fk()\f_k(x)\_x defines a GP (0,K)GP(0,K). 2.2 Limiting kernel and interpretation In the H→∞H→∞ regime, the model is no longer described inferentially by individual weights. Instead, it is fully characterized by a covariance function K(,′)K(x,x ), determined by h and the parameter priors. For (3), the limiting covariance takes the form K(,′)=σb2+σv2[h(Z)h(Z′)],K(x,x )\;=\; _b^2\;+\; _v^2\,E\! [h(Z)\,h(Z ) ], (5) where (Z,Z′)(Z,Z ) denotes the pair of pre-activations associated with (,′)(x,x ) induced by the priors on a and u (details in Appendix A.1). 2.3 Minimal regularity condition To ensure that (5) is well defined and that the output variance remains finite under (vkj)=σv2/HV(v_kj)= _v^2/H, it suffices to assume that [h(Z)2]<∞for Z distributed according to the implied prior on pre-activations.E[h(Z)^2]<∞ Z distributed according to the implied prior on pre-activations. (6) Common activations satisfy (6): tanh and sigmoid are bounded; ReLU and Leaky ReLU grow at most linearly and thus have a finite second moment under Gaussian Z. A more general sufficient condition (polynomial growth) and a proof are given in Appendix A.2. 2.4 General convergence theorem We state a general version of the convergence result (in the sense of finite-dimensional distributions) for shallow BNNs. Theorem 2.1 (Convergence of a shallow BNN to a GP). Consider the model (3) with the prior independent, centered parameters with finite variances and the scaling (vkj)=σv2/HV(v_kj)= _v^2/H. Assume h and the priors are such that Var(fk())<∞Var(f_k(x))<∞ for all x (e.g., [h(Z)2]<∞E[h(Z)^2]<∞). Then, for each k∈1,…,Kk∈\1,…,K\, as H→∞H→∞, fkf_k converges in law (in the sense of finite-dimensional distributions) to a Gaussian process fk(⋅)∼(0,K(⋅,⋅)),f_k(·) \! (0,\,K(·,·) ), with kernel given by (5), namely K(,′)=σb2+σv2[h(Z)h(Z′)],K(x,x )= _b^2+ _v^2\,E\! [h(Z)\,h(Z ) ], (7) where (Z,Z′)(Z,Z ) is the pre-activation pair associated with (,′)(x,x ) under the parameter priors. Proof sketch. For fixed inputs (1,…,n)(x_1,…,x_n), the vector (fk(1),…,fk(n)) (f_k(x_1),…,f_k(x_n) ) can be written as a sum of i.i.d. terms over hidden units j, with variance controlled by the 1/H1/H scaling. Under (6), the multivariate CLT (e.g. via Cramér–Wold) yields convergence to a multivariate normal with covariance K(,)K(X,X). Since this holds for any finite set of inputs, the limit defines a GP. The full proof is given in Appendix A.3. 3 Mixed kernel derived from a BNN In this section, we show how our mixed kernel arises directly from a BNN prior. The key point is that, under the standard scaling for wide hidden layers, the BNN prior induces a GP whose covariance function can be written explicitly as an expectation under a bivariate Gaussian distribution. Our mixed kernel is obtained by summing (and, optionally, adding) multiple kernels induced in this way, corresponding to different components of the BNN prior. 3.1 Kernels induced by nonlinear activations in BNNs By the infinite-width convergence result (Theorem 2.1) under the architecture 3, for fixed inputs ,′x,x the kernel induced by an activation h takes the form Kh(,′)=σb2+σv2[h(Z)h(Z′)],K_h(x,x )= _b^2+ _v^2\,E\! [h(Z)\,h(Z ) ], where (Z,Z′)(Z,Z ) is the pair of pre-activations induced by the prior and is a centered bivariate Gaussian. Writing σz2=Var(Z) _z^2=Var(Z), σz′2=Var(Z′) _z ^2=Var(Z ), and ρ=Corr(Z,Z′)ρ=Corr(Z,Z ), we obtain the integral representation [h(Z)h(Z′)]=∬ℝ2h(z)h(z′)ϕρ,σz,σz′(z,z′)zz′,E\! [h(Z)\,h(Z ) ]= _R^2h(z)\,h(z )\; _ρ, _z, _z (z,z )\,dz\,dz , where ϕρ,σz,σz′ _ρ, _z, _z denotes the bivariate Normal density with zero means, standard deviations σz,σz′ _z, _z , and correlation ρ. Next, we collect expressions available in the literature for the inner term [h(Z)h(Z′)]E[h(Z)h(Z )], in closed form when available or via standard approximations, for four reference activations: tanh. The integral does not admit a simple closed form; one uses the standard approximation tanh⁡(x)≈erf(π2x) (x) \! ( π2x ) (??), which yields the arcsine kernel: [tanh⁡(Z)tanh⁡(Z′)]≈2πarcsin⁡(π2ρσzσz′(1+π2σz2)(1+π2σz′2)).E[ (Z) (Z )]\;≈\; 2π\, \! ( π2\,ρ\, _z _z (1+ π2 _z^2 ) (1+ π2 _z ^2 ) ). (8) Sigmoid. Similarly (via an erferf approximation), one obtains: [sig(Z)sig(Z′)]≈14+12πarcsin⁡(π8ρσzσz′(1+π8σz2)(1+π8σz′2)).E[sig(Z)\,sig(Z )]\;≈\; 14+ 12π\, \! ( π8\,ρ\, _z _z (1+ π8 _z^2 ) (1+ π8 _z ^2 ) ). (9) ReLU. For ReLU(z)=max⁡(0,z)ReLU(z)= (0,z), a closed form is available: [ReLU(Z)ReLU(Z′)]=σzσz′2π[1−ρ2+ρ(π−arccos⁡(ρ))].E[ReLU(Z)\,ReLU(Z )]= _z\, _z 2π [ 1-ρ^2+ρ (π- (ρ) ) ]. (10) LeakyReLU. For LeakyReLUα(z)=max⁡(z,αz)LeakyReLU_α(z)= (z,α z), define S(ρ):=12π[1−ρ2+ρ(π−arccos⁡(ρ))],S(ρ):= 12π [ 1-ρ^2+ρ (π- (ρ) ) ], and obtain: [LeakyReLUα(Z)LeakyReLUα(Z′)]=σzσz′[αρ+(1−α)2S(ρ)].E[LeakyReLU_α(Z)\,LeakyReLU_α(Z )]= _z _z [α\,ρ+(1-α)^2\,S(ρ) ]. (11) Finally, substituting each of the expressions above into Kh(,′)=σb2+σv2[h(Z)h(Z′)]K_h(x,x )= _b^2+ _v^2\,E[h(Z)h(Z )], we obtain the corresponding induced kernels. The calculation steps and approximations are detailed Supplementary Material (Section S.2). 3.2 Mixed kernel from the GP limit of a BNN Each choice of activation function h determines, in the infinite-width limit, a specific covariance function through the term [h(Z)h(Z′)]E[h(Z)h(Z )]. In applications, it is often useful to combine features of different activations. In this work, we implement such a combination through an additive (block) construction, i.e., as a sum of M independent subnetworks, each using one activation hmh_m. This construction yields, in the GP limit, a kernel exactly in the form of an additive (and convex) combination of the component kernels, with no cross-terms. For simplicity, we first consider a single-output version (the extension to each output component k is immediate). Let zj()=aj+∑i=1Iujixi,j=1,…,H,z_j(x)=a_j+ _i=1^Iu_jix_i, j=1,…,H, and define the function prior f()=b+∑m=1Mwm∑j=1Hvj(m)hm(zj()),∑m=1Mwm=1,wm∈(0,1),f(x)=b+ _m=1^M w_m _j=1^Hv^(m)_j\,h_m\! (z_j(x) ), _m=1^Mw_m=1,\;\;w_m∈(0,1), (12) with vj(m)∼i.i.d.(0,σv2H),b∼(0,σb2),v^(m)_j i.i.d. N\! (0, _v^2H ), b (0, _b^2), and independence across the collections vj(m)j=1H\v^(m)_j\_j=1^H for different m. This independence (together with [vj(m)]=0E[v^(m)_j]=0) implies that, when expanding the covariance, all cross-terms with m≠m′m≠ m vanish, leading to an additive kernel in the GP limit. Moreover, to ensure that the CLT-based argument applies in the standard way, it suffices to assume that each component satisfies [hm(Z)2]<∞E[h_m(Z)^2]<∞ under the implied prior for the pre-activations; in that case, the output has finite variance and the GP limit is well defined. Proposition 3.1 (Convergence with mixed activations (additive blocks)). Consider the shallow network under the additive construction (12), with the same assumptions as in Theorem 2.1 for (aj,uji)(a_j,u_ji), and with output-weight blocks vj(m)\v^(m)_j\ independent across components m, each scaled as vj(m)∼(0,σv2/H)v^(m)_j (0, _v^2/H). Then, as H→∞H→∞, f converges in distribution (in the sense of finite-dimensional distributions) to a Gaussian process f(⋅)∼(0,K(⋅,⋅)),f(·) \! (0,K(·,·) ), with covariance function K(,′)=σb2+σv2∑m=1MwmKm(,′),Km(,′)=[hm(Z)hm(Z′)],K(x,x )= _b^2+ _v^2 _m=1^Mw_m\,K_m(x,x ), K_m(x,x )=E\! [h_m(Z)\,h_m(Z ) ], (13) where (Z,Z′)(Z,Z ) denotes the pre-activation pair associated with (,′)(x,x ) under the priors (e.g., a centered bivariate Normal under Gaussian (aj,uji)(a_j,u_ji)), as described in Appendix A.1. Proof idea. Write f()=b+∑m=1Mwmgm()f(x)=b+ _m=1^M w_m\,g_m(x), where gm()=∑j=1Hvj(m)hm(zj())g_m(x)= _j=1^Hv^(m)_j\,h_m(z_j(x)). Because the blocks vj(m)\v^(m)_j\ are independent across m and centered, cross-covariances vanish: Cov(gm(),gm′(′))=0Cov(g_m(x),g_m (x ))=0 for m≠m′m≠ m . For each fixed m, the standard wide-network argument (multivariate CLT with 1/H1/H scaling) yields Cov(gm(),gm(′))→σv2[hm(Z)hm(Z′)]Cov(g_m(x),g_m(x ))→ _v^2\,E[h_m(Z)h_m(Z )]. Summing over m and adding Var(b)=σb2Var(b)= _b^2 gives (13). A complete proof is given in Appendix A.5. 3.3 Identifiability of the mixed kernel The inner term [h(Z)h(Z′)]E[h(Z)\,h(Z )] fully determines the dependence structure of the process over the input space and, therefore, each choice of activation hmh_m induces a specific covariance term Km(,′)=[hm(Z)hm(Z′)].K_m(x,x )=E\! [h_m(Z)\,h_m(Z ) ]. In the infinite-width regime, the GP associated with the BNN has a mixture-form kernel K(,′)=σb2+σv2∑m=1MwmKm(,′),K(x,x )= _b^2+ _v^2 _m=1^Mw_m\,K_m(x,x ), (14) where we assume wm≥0w_m≥ 0 (and, if desired for a mixture interpretation, ∑mwm=1 _mw_m=1). Table 1: Pearson correlation between the shapes of the kernels (curves as functions of ρ with σz=σz′=1 _z= _z =1), after centering and normalization. Compared pairs Correlation tanh × sigmoid 0.9999 ReLU × LeakyReLU (α=0.1α=0.1) 0.9983 ReLU × LeakyReLU (α=0.3α=0.3) 0.9919 The expressions (8) to (11) report four inner terms KmK_m arising in the infinite-width limit. Although these kernels are distinct in principle, in practice they can be very close as functions of ρ(,′)ρ(x,x ) (for comparable scales of σz _z and σz′ _z ). As a consequence, different combinations of weights and hyperparameters may yield very similar covariance matrices over a finite set of inputs, leading to (near) non-identifiability and ill-conditioned estimation problems, a phenomenon that is typical in mixture models. In particular, monotone saturating activations such as tanh and sigmoid induce arc-sine family kernels (via the erf approximation) and produce smooth curves, monotone in ρ, with similar curvature around ρ=0ρ=0 and saturation as |ρ|→1|ρ|→ 1. Likewise, ReLU and LeakyReLU induce nearly colinear shapes: the LeakyReLU term can be viewed as a perturbation of ReLU by a linear term in ρ weighted by α, which for usual values of α leads to a numerically small difference. This proximity is reflected in Table 1, which reports very high Pearson correlations between the curves (after centering and normalization) as functions of ρ. In light of this, we adopt two representative components. For the smooth component, we take Ksmooth≡K1K_smooth≡ K_1, since K1≈K2K_1≈ K_2 under the conditions above. For the angular component, we choose Kangular≡K4K_angular≡ K_4 (LeakyReLU) as the representative of the pair, because K3≈K4K_3≈ K_4 and K4K_4 includes the parameter α that controls asymmetry: Ksmooth(,′) K_smooth(x,x ) ≡K1(,′)≈K2(,′), \;≡\;K_1(x,x )\;≈\;K_2(x,x ), Kangular(,′) K_angular(x,x ) ≡K4(,′)≈K3(,′). \;≡\;K_4(x,x )\;≈\;K_3(x,x ). Thus, we introduce a single mixing parameter w∈(0,1)w∈(0,1) and define the reduced mixed kernel Kmix(,′) K_mix(x,x ) =σb2+σv2[wKsmooth(,′)+(1−w)Kangular(,′)] = _b^2+ _v^2 [w\,K_smooth(x,x )+(1-w)\,K_angular(x,x ) ] (15) =σb2+σv2[w2πarcsin(π2ρσzσz′(1+π2σz2)(1+π2σz′2)) = _b^2+ _v^2 [w\, 2π\, \! ( π2\,ρ\, _z\, _z (1+ π2\, _z^2 )\, (1+ π2\, _z ^2 ) ) +(1−w)(αρσzσz′+(1−α)2σzσz′2π(1−ρ2+ρ(π−arccos(ρ))))]. +(1-w)\, (α\,ρ\, _z\, _z +(1-α)^2\, _z\, _z 2π\, ( 1-ρ^2+ρ (π- (ρ) ) ) ) ]. The kernel KmixK_mix is positive semidefinite because it is a convex combination (in w) of positive semidefinite kernels, plus the bias-noise term σb2 _b^2; in particular, linear combinations with nonnegative coefficients preserve positive semidefiniteness. Proposition 3.2 (Identifiability of the mixed kernel). Consider the zero-mean GP with covariance KmixK_mix in (15), with parameter space Θ=θ=(σa2,σu2,σb2,σv2,w,α):σa2,σu2,σb2,σv2>0,w∈(0,1),α∈(0,1). = \\,θ=( _a^2, _u^2, _b^2, _v^2,w,α)\;:\; _a^2, _u^2, _b^2, _v^2>0,\;w∈(0,1),\;α∈(0,1)\, \. Then the mapping Θ⟶,θ↦Kmix(⋅,⋅;θ), \; \;K, θ K_mix(·,·;θ), is injective in the interior of Θ . In particular, if for any θ1,θ2∈Θ _1, _2∈ we have Kmix(,′;θ1)=Kmix(,′;θ2),∀,′∈ℝd,K_mix(x,x ; _1)=K_mix(x,x ; _2), ∀\,x,x ^d, then necessarily θ1=θ2 _1= _2. Proof idea. The complete proof is given in Appendix A.5. The argument exploits families of pairs (,′)(x,x ) to separate parameters: (i) by varying the angle between vectors of fixed norm, the constant term forces σb,12=σb,22 _b,1^2= _b,2^2; (i) using orthogonal pairs with ‖=‖′‖=r\|x\|=\|x \|=r and varying r, one identifies (σa2,σu2)( _a^2, _u^2) via the induced expressions for ρ(r)ρ(r) and σz2(r) _z^2(r); (i) with these fixed, equality of the kernel implies equality of the remaining parameters (σv2,w,α)( _v^2,w,α), yielding θ1=θ2 _1= _2 in the interior of Θ . Boundary cases (w∈0,1w∈\0,1\ or α∈0,1α∈\0,1\) are excluded since they may remove parameters from the kernel. 3.4 Stratified input designs and practical identifiability Proposition 3.2 establishes theoretical identifiability of Kmix(⋅,⋅;θ)K_mix(·,·;θ) in the interior of Θ . In practice, however, identifiability can become numerically fragile in moderate-to-high dimension when inputs ii=1n\x_i\_i=1^n are sampled i.i.d.. The concentration of norms and angles around typical values, which is characteristic of high-dimensional spaces, implies that the covariance matrix (θ)K(θ) is evaluated only over a limited range of radial and angular variation (????). As a result, the effectively observable information for distinguishing kernel components is reduced, which can lead to numerically ill-conditioned estimation (??????). Let (θ)∈ℝn×nK(θ) ^n× n be the matrix with entries ij(θ)=Kmix(i,j;θ)K_ij(θ)=K_mix(x_i,x_j;θ). Information about θ is determined by two geometric aspects of the input set. The first is the dispersion of the norms ‖i‖\|x_i\|, which governs marginal variances, and the second is the dispersion of the angles between pairs (i,j)(x_i,x_j), which governs covariances. These two sources of variability feed, respectively, the diagonal and off-diagonal terms of (θ)K(θ), as follows. • Diagonal (radial information): For the models considered, the pre-activation variance satisfies σz2()=σu2‖2+σa2 _z^2(x)= _u^2\|x\|^2+ _a^2. Thus, separating (σu2,σa2)( _u^2, _a^2) depends on observing sufficient dispersion in ‖i‖2\|x_i\|^2 across inputs. • Off-diagonal (angular information): For i≠ji≠ j, the kernel depends on the pre-activation correlation ρ(i,j)ρ(x_i,x_j), which is a monotone function of the cosine similarity ρ~(i,j)=i⊤j/(‖i‖‖j‖) ρ(x_i,x_j)=x_i x_j/(\|x_i\|\|x_j\|). Therefore, distinguishing the smooth component A(ρ)A(ρ) from the angular component B(ρ;α)B(ρ;α) requires observing ρ over a sufficiently wide range. Under i.i.d. sampling in high dimension, both ‖2/I\|x\|^2/I and ρ~(,′) ρ(x,x ) concentrate around typical values, with fluctuations that decay as O(I−1/2)O(I^-1/2) (????). In this regime, different values of θ may induce very similar matrices (θ)K(θ) in finite samples, which hampers numerical separation between kernel components and can degrade estimation stability (??????). To make identifiability operational in practice, we compare two input designs, a uniform (i.i.d.) design and a stratified design. The stratified design preserves variation along the two relevant axes of the kernel, with radial stratification to increase dispersion in ‖\|x\| and angular stratification to increase the diversity of similarities between pairs of points. This choice allows us to study, in a controlled way, how the geometry of the input design affects numerical identifiability and estimation stability in the simulation scenarios, isolating this effect from other factors. This type of construction is consistent with recommendations for space-filling designs to avoid clustering, reduce voids, and increase coverage of the input domain (??). It is important to emphasize that this stratification is not proposed as a strategy for real datasets, where the inputs are already given and, in general, cannot be redesigned. Here, it is used only in the simulations as a tool to investigate, in high dimension and in a controlled manner, how the geometry of the input design affects numerical identifiability and estimation stability. Supporting calculations and implementation details are provided in Supplementary Section S2.2. 4 MAP estimation via Nyström and Bayesian prediction In this section, we describe MAP estimation of the hyperparameters of the GP induced by the mixed kernel, and we introduce the Nyström approximation to make the procedure scalable when n is large. 4.1 Optimization via MAP Let =(i,yi)i=1nD=\(x_i,y_i)\_i=1^n and let θ denote the set of model hyperparameters. MAP estimation is equivalent to minimizing the objective (log-posterior density) ℒMAP(θ)=−log⁡p(∣,θ)−log⁡p(θ),L_MAP(θ)=- p(y ,θ)\;-\; p(θ), (16) where the first term is the log-likelihood induced by the GP measure and the second term is the contribution of the log-prior density (acting as probabilistic regularization). Updates are performed by gradient descent: θ(t+1)=θ(t)−η∇θℒMAP(θ).θ^(t+1)=θ^(t)-η\, _θ\,L_MAP(θ). In practice, gradients are obtained by automatic differentiation applied to the numerical construction of ~θ K_θ (e.g., via a Cholesky factorization or low-rank approximations). 4.2 GP log-likelihood We assume a regression model with independent Gaussian noise: yi=f(i)+εi,εi∼(0,σϵ2),y_i=f(x_i)+ _i, _i (0, _ε^2), and, under a zero-mean GP with kernel Kmix(⋅,⋅;θ)K_mix(·,·;θ), we have ∣,θ∼(,θ),∣,θ∼(,~θ),f ,θ (0,K_θ), ,θ (0, K_θ), where ~θ=θ+σϵ2,[θ]ij=Kmix(i,j;θ). K_θ=K_θ+ _ε^2I, [K_θ]_ij=K_mix(x_i,x_j;θ). (17) Therefore, the log-likelihood takes the standard GP form: log⁡p(∣,θ)=−12⊤~θ−1−12logdet(~θ)−n2log⁡(2π). p(y ,θ)=- 12\,y K_θ^-1y- 12\, ( K_θ)- n2 (2π). (18) Mixed-kernel form (reduced model). In our case, θ=σb2 11⊤+σv2[wsmooth(σa2,σu2)+(1−w)angular(σa2,σu2,α)],K_θ= _b^2\,11 + _v^2 [w\,K_smooth( _a^2, _u^2)+(1-w)\,K_angular( _a^2, _u^2,α) ], (19) with θ=σϵ2,σa2,σu2,σb2,σv2,α,wθ=\ _ε^2, _a^2, _u^2, _b^2, _v^2,α,w\. 4.3 MAP loss and gradients (standard form) With the priors σq2∼Inv-Gamma(aq,bq),q∈ϵ,a,u,b,v,α∼Beta(aα,bα),w∼Beta(aw,bw), _q^2 -Gamma(a_q,b_q), q∈\ε,a,u,b,v\, α (a_α,b_α), w (a_w,b_w), the MAP loss (up to additive constants) is ℒMAP(θ) _MAP(θ) =12⊤~θ−1+12logdet(~θ)+∑q∈ϵ,a,u,b,v[(aq+1)log⁡σq2+bqσq2] = 12\,y K_θ^-1y+ 12\, ( K_θ)+ _q∈\ε,a,u,b,v\ [(a_q+1) _q^2+ b_q _q^2 ] −(aα−1)log⁡α−(bα−1)log⁡(1−α)−(aw−1)log⁡w−(bw−1)log⁡(1−w). -\,(a_α-1) α\;-\;(b_α-1) (1-α)-\,(a_w-1) w\;-\;(b_w-1) (1-w). (20) For the likelihood term in (18), the gradient with respect to a generic parameter θj _j can be written compactly as ∂θj(−log⁡p(∣,θ))=12tr[(~θ−1⊤~θ−1−~θ−1)∂~θ∂θj]. ∂ _j (- p(y ,θ) )= 12\,tr\! [ ( K_θ^-1yy K_θ^-1- K_θ^-1 )\; ∂ K_θ∂ _j ]. (21) The specific derivatives ∂~θ/∂θj∂ K_θ/∂ _j (e.g., for w, α, and the variance parameters) follow directly from (17)–(19) and are listed in Supplementary Material (Section S.3) for reference. From (19), ∂~θ∂w=σv2(smooth−angular), ∂ K_θ∂ w= _v^2 (K_smooth-K_angular ), and the corresponding gradient follows from (21). 4.4 Nyström approximation Let ∈ℝn×nK ^n× n denote the covariance matrix induced by a kernel Kθ(⋅,⋅)K_θ(·,·). In this work, we use the approach proposed in (??), which we refer to as the Nyström method, as it approximates K by a low-rank decomposition constructed from a subset S⊂1,…,nS⊂\1,…,n\ with |S|=r|S|=r (the “anchors”). As shown in Eq. 20, MAP optimization involves logdet(~θ) ( K_θ) and products with ~θ−1 K_θ^-1. Therefore, we apply the Nyström approximation directly to the noisy covariance matrix ~θ K_θ: ~θ(r)=~θ,:S~θ,SS−1~θ,S:,~θ=θ+σϵ2,Swith|S|=r. K_θ^(r)= K_θ,:S\, K_θ,S^-1\, K_θ,S:, K_θ=K_θ+ _ε^2I, S\;with\;|S|=r. (22) where S⊂1,…,nS⊂\1,…,n\ is a subset of size r. With this construction, the dominant computational cost (e.g., forming the required terms and solving the associated linear systems) typically reduces to (nr2)O(nr^2) instead of (n3)O(n^3). Anchors S can be selected by different strategies, such as random sampling, k-means++ initialization (or centroids obtained via k-means) (??), or randomized sampling/projection techniques that approximate the column space of K (????). In this work, we adopt two criteria depending on the experimental context: • first, where S=1,…,rS=\1,…,r\ (the first r points), used in simulations (via Vecchia approximation) due to its simplicity and reproducibility; or • k-means++, which selects the first centroid uniformly at random from the data and then selects each subsequent centroid with probability proportional to the squared distance to the nearest already-chosen centroid (??). In summary, the Nyström approximation avoids the exact factorization of an n×n× n matrix—the main source of the (n3)O(n^3) cost in GPs—by replacing it with low-rank operations. This makes it practical to compute logdet(~θ) ( K_θ) and solve ~θ= K_θu=y within MAP optimization on large-scale datasets. 4.5 Bayesian prediction In a fully Bayesian analysis, uncertainty about the model parameters θ is propagated to prediction through the posterior predictive distribution: p(y⋆∣,⋆)=∫p(y⋆∣⋆,θ)p(θ∣)θ.p(y_ ,x_ )= p(y_ _ ,θ)\,p(θ )\,dθ. Thus, prediction is not obtained by conditioning on a single point estimate of θ, but by averaging the conditional predictions p(y⋆∣⋆,θ)p(y_ _ ,θ) weighted by the posterior distribution p(θ∣)p(θ ). As a result, predictive intervals reflect both the intrinsic variability of the model and uncertainty about θ. Although conceptually appealing, this approach is typically infeasible in large-scale problems, since evaluating the posterior predictive distribution requires repeatedly recomputing quantities based on the covariance matrix ~θ K_θ, including costly terms involving ~θ−1 K_θ^-1. To balance uncertainty quantification and computational cost, we estimate the parameters via MAP, obtaining a point estimate θ θ, and perform prediction conditional on θ θ. Although this does not explicitly integrate over parameter uncertainty, predictions remain probabilistic because the GP provides a full predictive distribution for y⋆y_ . Accordingly, let ∈ℝn×IX ^n× I denote the training inputs and ∈ℝny ^n the observed responses. We define ~θ^=θ^+σ^ϵ2,[θ^]ij=Kmix(i,j;θ^). K_ θ=K_ θ+ σ_ε^2\,I, [K_ θ]_ij=K_mix(x_i,x_j; θ). For a new input ⋆x_ , we let ⋆=(Kmix(⋆,1;θ^),…,Kmix(⋆,n;θ^))⊤,k⋆=Kmix(⋆,⋆;θ^).k_ = (K_mix(x_ ,x_1; θ),…,K_mix(x_ ,x_n; θ) ) , k_ =K_mix(x_ ,x_ ; θ). Then, the posterior predictive distribution is Gaussian: y⋆∣,,⋆,θ^∼(μ⋆,σ⋆2),y_ ,X,x_ , θ \! ( _ , _ ^2 ), (23) with mean and variance given by μ⋆=⋆⊤~θ^−1,σ⋆2=k⋆−⋆⊤~θ^−1⋆+σ^ϵ2. _ =k_ K_ θ^-1y, _ ^2=k_ -k_ K_ θ^-1k_ + σ_ε^2. (24) Note that, given θ θ, GP prediction requires manipulating the covariance matrix ~θ K_ θ. Computing and manipulating this matrix exactly is impractical at large scale. Hence, although the predictive expressions are written in terms of ~θ^−1 K_ θ^-1, in practice we replace ~θ K_ θ by its Nyström approximation ~θ^(r) K_ θ^(r) (Eq. 22) and compute the products (~θ^(r))−1 ( K_ θ^(r) )^-1y and (~θ^(r))−1⋆ ( K_ θ^(r) )^-1k_ via low-rank algebra (equivalently, by solving linear systems in the approximated space), avoiding explicit matrix inversion and enabling scalable evaluation of predictive means and variances (??). The proposed GP model, conditioned on θ θ and using the Nyström approximation, offers important advantages over finite NNs. First, whereas in finite networks the number of parameters grows with the number of neurons H, the GP limit is governed by a small set of parameters, such as the variances (σa2,σu2,σb2,σv2,σϵ2)( _a^2, _u^2, _b^2, _v^2, _ε^2) and, when applicable, parameters indexing the activation function, which in turn determine the covariance function. Second, uncertainty quantification is more faithful and robust, especially outside the data domain. In deterministic NNs, prediction is purely pointwise and there is no predictive distribution. Moreover, in finite NNs with a nugget term, predictive variance tends to be approximately constant and dominated by the noise term, remaining small even when ⋆x_ is far from the training data. In contrast, in the proposed GP, predictive variance includes an epistemic component that naturally increases as ⋆x_ becomes less correlated with the training inputs, approaching the prior variance in regions poorly informed by data. Thus, even with parameters fixed via MAP, the GP consistently captures both aleatoric and epistemic uncertainty. Finally, the predictive mean μ⋆ _ and predictive variance σ⋆2 _ ^2 (Eqs. 23–24) provide the basis for computing the evaluation metrics: the Mean absolute error (MAE) and the Root mean squared error (RMSE) assess point-prediction accuracy using μ⋆ _ , whereas the Mean expected squared error (MESE) and the Standard deviation of the expected squared error (SDESE) jointly assess accuracy and the quality of uncertainty quantification by explicitly incorporating σ⋆2 _ ^2 (see Section 5.3 for a formal definition of those metrics). Under the Gaussian assumption, the expected squared error decomposes as [(Y⋆−y⋆)2∣,⋆,θ^]=(μ⋆−y⋆)2+σ⋆2,E\! [(Y_ -y_ )^2 ,x_ , θ ]=( _ -y_ )^2+ _ ^2, i.e., the sum of the predictive-mean error and the predictive uncertainty, reflecting both accuracy and the reliability of the reported uncertainty. 5 Experimental Setup We evaluate the proposed method on both simulated data and two public regression datasets. In the simulations, we vary the sample size n, the input dimension I, and the sampling design for X (uniform vs. stratified), and we compare anchor-selection strategies for the Nyström approximation. On real data, we report predictive performance and total runtime on Superconductivity and YearPredictionMSD, including a scalability study on the full YearPredictionMSD dataset by varying the Nyström rank r. 5.1 Simulated data: generation, scenarios, and replications We consider eight experimental scenarios. In each scenario, inputs i∈[0,1]Ix_i∈[0,1]^I are generated under two sampling designs: (i) uniform and (i) stratified over the unit hypercube. We then apply a centering transformation, ~i=i−0.5, x_i=x_i-0.5, ensuring ~i∈[−0.5,0.5]I x_i∈[-0.5,0.5]^I, consistent with the model parameterization. The latent function is defined as a Gaussian process f∼(0,Kmix)f (0,K_mix). Observations are obtained as y(~)=f(~)+ε(~),ε(~)∼iid(0,σε2),y( x)=f( x)+ ( x), ( x) iid N(0, _ ^2), with ε assumed independent of f. The nugget ε is calibrated per scenario as a fixed fraction η of the average marginal variance of the kernel evaluated at the input locations: σε2=ηK¯,K¯=1n∑i=1nKmix(~i,~i), _ ^2=η\, K, K= 1n _i=1^nK_mix( x_i, x_i), (25) with η=0.04η=0.04. In the implementation, K¯ K is computed in batches without forming the full n×n× n matrix. To enable scalable GP sampling in regimes with large n without incurring cubic cost, we generate f using a sequential Vecchia-type construction: we sample an initial exact block of size Ninit=500N_init=500 via Cholesky, and for each new point, we sample the univariate conditional distribution based on the Nviz=500N_viz=500 nearest neighbors among the previously simulated ones. This procedure is used only in the simulation stage (data generation). We vary n∈10.000,20.000,50.000n∈\10.000,20.000,50.000\, I∈20,80I∈\20,80\ and, for n=50.000n=50.000, the design of X (uniform vs. stratified), yielding eight scenarios (C1–C8) summarized in Table 2. To quantify variability induced by the randomness of the latent process and observation noise, we perform R=20R=20 replications for scenario C1, keeping the inputs fixed and varying only seed_y. A summary of the replications is reported in the Supplementary Material (Table S.1), and full details are provided in the Supplementary Material (Section S.5). Table 2: Simulation scenarios using a Vecchia-type construction with a nugget automatically calibrated from Kmix(x,x)¯ K_mix(x,x). Scenario Design of X I n (,)¯ K_mix(x,x) Nugget (σε)( _ ) _ ^2 Time (s) C1 Uniform 20 10000 2.132745 0.292078 0.085310 63.72 C2 Uniform 80 10000 3.769967 0.388328 0.150799 34.90 C3 Uniform 20 20000 2.131981 0.292026 0.085279 68.46 C4 Uniform 80 20000 3.769044 0.388281 0.150762 74.07 C5 Uniform 20 50000 2.131707 0.292007 0.085268 182.54 C6 Uniform 80 50000 3.768044 0.388229 0.150722 215.66 C7 Stratified 20 50000 2.132325 0.292050 0.085293 186.13 C8 Stratified 80 50000 3.772796 0.388474 0.150912 227.50 Table 2 defines the eight simulation regimes by varying n, the input dimension I, and the design of X. Automatically calibrating the nugget as a fixed fraction of the average marginal variance of Kmix(,)K_mix(x,x) standardizes the noise scale across scenarios, avoiding ad hoc choices when KmixK_mix changes with I or with the sampling design. 5.2 Real data and preprocessing We evaluate the method on the public regression benchmarks Superconductivity and YearPredictionMSD. For YearPredictionMSD, we analyze random subsamples with n∈10,000,20,000,50,000n∈\10,000,20,000,50,000\ and, additionally, the full dataset while varying the Nyström rank r, in order to characterize the cost–accuracy trade-off of the approximation. Let X∈ℝn×IX ^n× I denote the covariate matrix and y∈ℝny ^n the response vector. We rescale each column of X to [0,1]I[0,1]^I using min–max normalization (fit on the training set and applied to the test set), and we center the inputs via X~=X−0.5 X=X-0.5. We standardize the target on the training set via y~=(y−μy)/σy y=(y- _y)/ _y and de-standardize predictions to report all metrics on the original scale. All datasets have no missing values; hence no imputation was required. 5.3 Training and evaluation protocol To ensure direct comparability across scenarios, we adopt a fixed-budget training protocol: we train the model for 5050 epochs in all configurations, keeping all other optimization choices unchanged. Model fitting is performed via MAP estimation (Eq. 20), by minimizing ℒMAPGP(θ)=NLL(θ)−log⁡p(θ),NLL(θ)=−log⁡p(∣,θ),L_MAP^GP(θ)=NLL(θ)- p(θ), (θ)=- p(y ,θ), where p(θ)p(θ) denotes the prior. For scalability, we use the Nyström approximation applied to the mixed kernel KmixK_mix. In the simulation scenarios, we fix the rank at r=500r=500. For real data, we use r=500r=500 on the full Superconductivity dataset and on the YearPredictionMSD subsamples, so as to keep computational cost comparable across configurations. Only in the experiment on the full YearPredictionMSD dataset do we vary r, to characterize the cost–accuracy trade-off. Anchor selection is performed either by first or by k-means, depending on the experiment. In all runs, we optimize the parameters with Adam (learning rate 10−310^-3), computing objective terms with mini-batches. The prior term is always included in the objective. The nugget is estimated jointly with the remaining parameters, with a dedicated learning rate of 10−310^-3, and is automatically initialized using the rule in Eq. 25 (with η=0.04η=0.04). When available, kernel parameters are initialized with the values used in the simulation. Across all datasets, we use reproducible splits (shared random seeds) and reserve 10%10\% for testing. After training on trainD_train, we evaluate predictive performance on the held-out test set. For each test input ix_i, the fitted GP yields a Gaussian predictive distribution Yi∣train,i∼(y^i,σ^i2)Y_i _train,x_i ( y_i, σ_i^2). Let (y^i,σ^i2)( y_i, σ_i^2) denote the predictive mean and variance for each test point (i,yi)(x_i,y_i), i=1,…,ntesti=1,…,n_test. Then, we compute point-accuracy metrics, MAE=1ntest∑i=1ntest|y^i−yi|,MSE=1ntest∑i=1ntest(y^i−yi)2,RMSE=MSE.MAE= 1n_test _i=1^n_test y_i-y_i , = 1n_test _i=1^n_test( y_i-y_i)^2, = MSE. To jointly assess accuracy and uncertainty, we use the Expected squared error (ESE) under the predictive distribution, ESEi=[(Yi−yi)2∣train,i]=(y^i−yi)2+σ^i2,ESE_i=E\! [(Y_i-y_i)^2 _train,x_i ]=( y_i-y_i)^2+ σ_i^2, and report MESE=1ntest∑i=1ntestESEi,SDESE=1ntest−1∑i=1ntest(ESEi−MESE)2.MESE= 1n_test _i=1^n_testESE_i, = 1n_test-1 _i=1^n_test (ESE_i-MESE )^2. For real datasets, in addition to metrics on the original target scale, we also report normalized versions based on the training-target scale to enable comparisons across datasets with different units (Supplementary Material S5.4). Overall, the experimental design and protocol above define a reproducible procedure to compare, in a controlled manner, the effect of the simulation regime, the Nyström anchor-selection strategy, and, in the large-scale setting, the rank r. In the next section, we report MAP estimates, test-set predictive metrics, and total time, highlighting the cost–accuracy trade-off induced by the scalable approximation. 5.4 Results 5.4.1 Simulation results The Tables 3 and 4 disentangle two effects: (i) the stability of the MAP estimates of the kernel, and (i) the impact of the anchor-selection strategy on the quality of the Nyström approximation. Across all scenarios, the estimated structural parameters (σ^b2,σ^v2,σ^u2,σ^a2,α^,w^)( σ_b^2, σ_v^2, σ_u^2, σ_a^2, α, w) remain very similar under first and k-means, suggesting stable optimization under a fixed training budget, with no evidence of collapse to degenerate cases (e.g., w∈0,1w∈\0,1\). In contrast, the test-set metrics highlight a clear cost–accuracy trade-off. In more challenging settings (notably I=80I=80), k-means anchor selection tends to yield more noticeable reductions in MAE/RMSE, whereas in lower dimension the differences are small. This pattern is consistent with the role of anchors in Nyström, since better coverage of the input space improves the approximation of the correlated component of the process. As expected, these gains may come at increased computational cost due to clustering, which is reflected in the total runtime. Table 3: MAP estimates under two anchor-selection strategies in Nyström. Scenario Anchors σ^b2 σ_b^2 σ^v2 σ_v^2 σ^u2 σ_u^2 σ^a2 σ_a^2 α α w w σ^ϵ2 σ_ε^2 C1 first 1.025371 0.980080 0.981301 0.994972 0.508588 0.495693 8.56×10−28.56× 10^-2 C1 k-means 1.026216 0.987962 0.995776 0.987359 0.508575 0.496666 8.29×10−28.29× 10^-2 C2 first 1.015887 1.031211 1.030972 0.968986 0.487482 0.512434 1.64×10−11.64× 10^-1 C2 k-means 1.028323 1.032194 1.032033 0.968549 0.487544 0.511323 1.63×10−11.63× 10^-1 C3 first 1.028055 0.971542 0.972107 0.983965 0.511472 0.491441 8.64×10−28.64× 10^-2 C3 k-means 1.023113 0.974400 0.976894 0.996529 0.510373 0.490998 8.18×10−28.18× 10^-2 C4 first 0.999669 1.031572 1.031348 0.968746 0.487415 0.511870 1.64×10−11.64× 10^-1 C4 k-means 1.018052 1.031881 1.031831 0.968938 0.487784 0.509282 1.63×10−11.63× 10^-1 C5 first 1.028704 0.971421 0.971586 0.994188 0.511165 0.490527 8.34×10−28.34× 10^-2 C5 k-means 1.025100 0.969336 0.969129 1.001332 0.511568 0.489583 8.21×10−28.21× 10^-2 C6 first 1.001530 1.031871 1.031808 0.969091 0.487441 0.512053 1.64×10−11.64× 10^-1 C6 k-means 1.014398 1.031409 1.031123 0.969269 0.487602 0.511393 1.64×10−11.64× 10^-1 C7 first 1.025906 0.971412 0.971880 1.002389 0.511691 0.489485 8.32×10−28.32× 10^-2 C7 k-means 1.024638 0.971043 0.972473 1.013075 0.512253 0.489012 8.33×10−28.33× 10^-2 C8 first 0.983979 1.031795 1.031675 0.968993 0.487443 0.512017 1.64×10−11.64× 10^-1 C8 k-means 1.018084 1.031796 1.031743 0.969018 0.487622 0.511714 1.64×10−11.64× 10^-1 Table 4: Total time and test-set metrics under two Nyström anchor-selection strategies. Scenario Anchors Total time (s) MAE MSE RMSE MESE SDESE C1 first 1.334 0.265983 0.112261 0.335053 0.236112 0.165608 C1 k-means 2.155 0.263822 0.110102 0.331816 0.233065 0.161552 C2 first 1.284 0.573294 0.504148 0.710034 0.892681 0.665183 C2 k-means 2.514 0.534228 0.440660 0.663822 0.829185 0.585415 C3 first 1.881 0.265764 0.110175 0.331926 0.233302 0.152864 C3 k-means 2.137 0.264066 0.109563 0.331004 0.229287 0.153301 C4 first 2.101 0.571924 0.515084 0.717694 0.903406 0.718319 C4 k-means 2.334 0.562438 0.498106 0.705766 0.886464 0.707048 C5 first 1.623 0.269528 0.114101 0.337789 0.233814 0.160743 C5 k-means 2.408 0.267161 0.111680 0.334185 0.231012 0.159149 C6 first 1.550 0.593175 0.562680 0.750120 0.950702 0.807495 C6 k-means 3.511 0.566497 0.503485 0.709567 0.892490 0.718490 C7 first 1.389 0.270088 0.115020 0.339146 0.234563 0.165507 C7 k-means 2.592 0.270696 0.116113 0.340753 0.237124 0.168015 C8 first 1.335 0.583224 0.537355 0.733045 0.927608 0.777588 C8 k-means 1.395 0.581387 0.538477 0.733810 0.929836 0.774792 To isolate the effect of the input design, we compare pairs with fixed n and I that differ only in the sampling design of X. In low dimension (I=20I=20), the design effect is small, with C5 and C7 yielding very similar metrics. In high dimension (I=80I=80), the effect depends on the anchor-selection strategy. With first, the stratified design (C8) slightly improves over the uniform design (C6), whereas with k-means the opposite occurs, with C8 performing worse than C6. Overall, the stratified design does not yield systematic gains its effect is secondary at I=20I=20 and, at I=80I=80, it interacts with the choice of anchors. This pattern is consistent with the idea that, in more complex regimes, input-space coverage and the anchor-selection mechanism can interact to affect approximation quality. Table 5: Summary of the 20 replications of scenario C1 (r=500r=500, k-means anchors, ntrain=9000n_train=9000). Quantidade Média (dp) [min, max] MAE 0.273514 (0.004816) [0.263822, 0.282890] RMSE 0.342126 (0.005893) [0.331816, 0.353658] MESE 0.239397 (0.004587) [0.231336, 0.246951] SDESE 0.167666 (0.005495) [0.161552, 0.179720] total_s 2.486 (0.616) [1.829, 4.463] σε2 _ ^2 0.08341 (0.00142) [0.08033, 0.08700] w 0.499987 (0.002367) [0.494409, 0.504269] alpha 0.509978 (0.001317) [0.508133, 0.513169] σb2 _b^2 1.028096 (0.001465) [1.025079, 1.031158] σv2 _v^2 0.975803 (0.006260) [0.963308, 0.987962] σu2 _u^2 0.982022 (0.006519) [0.969905, 0.995776] σa2 _a^2 0.980596 (0.007583) [0.965003, 0.999919] The results (Table 5) show that test-set metrics exhibit non-negligible variability across replications (e.g., mean RMSE ≈0.342≈ 0.342 with sd ≈0.006≈ 0.006), while the MAP estimates of the parameters remain stable (e.g., w≈0.5w≈ 0.5 with low dispersion). This suggests that Nyström-based fitting is numerically stable in scenario C1 and reinforces the importance of reporting performance as mean± , avoiding conclusions drawn from a single realization. To contextualize the magnitude of this variability, note that the replications were conducted using k-means anchors. The comparison between anchor-selection strategies (first vs. k-means) is presented separately for scenario C1 (Tables 3–4), where k-means typically yields a slight error reduction at the cost of increased runtime due to the additional clustering step. Taken together, the replications indicate that the MAP estimates remain stable, while predictive performance varies in a manner consistent with the randomness of the latent process and observation noise. Next, we evaluate the method on real datasets and examine its behavior in the large-scale regime. 5.4.2 Real data To maintain comparability with the simulation scenarios and control computational cost, we use k-means anchor selection for both the Superconductivity dataset and the YearPredictionMSD subsamples. For the full YearPredictionMSD dataset, we instead adopt the first strategy, aiming for improved scalability and lower preprocessing overhead. The results in Tables 6 and 7 show that the MAP estimates are stable across the YearPredictionMSD settings, in particular with α^≈0.488 α≈ 0.488 and w^≈0.512 w≈ 0.512, indicating that the fit preserves a balanced mixture between the components of KmixK_mix and does not collapse to degenerate cases (w∈0,1w∈\0,1\). Since training is carried out on the standardized target scale, the nugget is estimated in that scale; after de-standardization, its magnitude reflects the problem scale in the original target units (Table 8). In particular, we obtain σ^ε2(z)≈4.4×10−4 σ_ ^2(z)≈ 4.4× 10^-4, corresponding to σ^ε(z)≈0.021 σ_ (z)≈ 0.021. Table 6: MAP estimates of the model parameters on the real-world datasets under k-means anchor selection. Case n I σ^b2 σ_b^2 σ^v2 σ_v^2 σ^u2 σ_u^2 σ^a2 σ_a^2 α α w w σ^ε2 σ_ ^2 S1_super_full_all 21263 79 0.983783 1.031285 1.031180 0.969412 0.487748 0.512268 5.1759e-01 R1_year_full_10k 10000 90 1.028279 1.032023 1.031995 0.968947 0.487489 0.512480 5.4418e-02 R2_year_full_20k 20000 90 1.021539 1.031795 1.031790 0.969213 0.487546 0.512372 5.3020e-02 R3_year_full_50k 50000 90 1.027704 1.032048 1.032012 0.968841 0.487509 0.512499 5.3630e-02 Table 7: Test-set performance metrics reported on the original target scale and total time on the real datasets. Case Dataset n I Total time (s) MAE RMSE MESE SDESE S1_super_full_all superconduct 21263 79 1.521 11.996 18.464 347.035 808.186 R1_year_full_10k yearpredictionmsd 10000 90 2.141 8.151 11.346 129.177 276.568 R2_year_full_20k yearpredictionmsd 20000 90 1.593 8.483 11.578 134.426 272.604 R3_year_full_50k yearpredictionmsd 50000 90 2.266 8.551 11.614 135.163 292.555 Table 8: Scaled metrics normalized by the training-target standardization, relative to Table 7. Case sd(ytrain)sd(y_train) MAEz RMSEz σ^ε2(z) σ_ ^2(z) σ^ε2(orig) σ_ ^2(orig) S1_super_full_all 34.231 0.350 0.539 4.417e-04 5.1759e-01 R1_year_full_10k 11.104 0.734 1.022 4.414e-04 5.4418e-02 R2_year_full_20k 10.957 0.774 1.057 4.416e-04 5.3020e-02 R3_year_full_50k 11.022 0.776 1.054 4.415e-04 5.3630e-02 Predictively, we observe only marginal changes as n increases in YearPredictionMSD while keeping r=500r=500: RMSE remains around 11.3511.35–11.6111.61 (original scale), suggesting that, under a fixed-rank approximation and a fixed training budget, part of the error becomes dominated by the low-rank Nyström approximation. Even so, computational cost remains low (total time ≈1.5≈ 1.5–2.32.3 s per setting), highlighting scalability and reproducibility on large tabular datasets. Finally, because the targets are measured in different natural units (years in YearPredictionMSD and Kelvin in Superconductivity), absolute MAE/RMSE values on the original scale are not directly comparable across datasets. Therefore, in addition to reporting metrics on the original scale, we also report versions normalized by the training-target scale (MAEz, RMSEz) in Table 8, which allow comparisons of relative error magnitudes across settings and datasets. Table 9: Full YearPredictionMSD: MAP estimates of the model parameters under Nyström, using first anchors, for different ranks r. r σ^b2 σ_b^2 σ^v2 σ_v^2 σ^u2 σ_u^2 σ^a2 σ_a^2 α α w w σ^ε2 σ^2_ 1000 0.999745 1.031698 1.031808 0.968802 0.487477 0.512471 5.280451×10−2\,5.280451× 10^-2\, 1500 1.001928 1.030765 1.031028 0.969650 0.487674 0.512086 5.281479×10−2\,5.281479× 10^-2\, 2000 1.000840 1.031140 1.031370 0.969235 0.487561 0.512327 5.281858×10−2\,5.281858× 10^-2\, 2500 1.000225 1.031369 1.031605 0.968976 0.487469 0.512463 5.281404×10−2\,5.281404× 10^-2\, 3000 0.999901 1.031510 1.031739 0.968839 0.487441 0.512496 5.282333×10−2\,5.282333× 10^-2\, 3500 0.999993 1.031638 1.031865 0.968724 0.487422 0.512520 5.282105×10−2\,5.282105× 10^-2\, 4000 1.000136 1.031792 1.032011 0.968610 0.487406 0.512538 5.282208×10−2\,5.282208× 10^-2\, 4500 0.999972 1.017239 1.018181 0.978938 0.492810 0.523093 5.280713×10−2\,5.280713× 10^-2\, 5000 1.000011 1.022378 1.022732 0.973577 0.490332 0.518135 5.280737×10−2\,5.280737× 10^-2\, At this stage, we evaluate the GP model with the mixed kernel KmixK_mix on the full YearPredictionMSD dataset (UCI), using a random split with 10%10\% held out for testing and the same preprocessing adopted in the other real-data experiments: min–max rescaling of X to [0,1]I[0,1]^I on the training set (and applied to the test set), centering X~=X−0.5 X=X-0.5, and standardizing the target on the training set, with predictions de-standardized to report metrics on the original scale. We fit the model for each value of the rank r, keeping the anchor selection fixed to the first strategy, in order to isolate the effect of r on the quality of the Nyström approximation. We choose first on the full dataset for scalability, since selecting anchors via k-means would introduce an additional preprocessing cost that becomes significant at this scale. Table 10: Full YearPredictionMSD: test-set performance on the original target scale and total time under Nyström with first anchors and different ranks r. r Total time (s) MAE RMSE MESE SDESE 1000 5.284 7.446 10.262 105.416 231.199 1500 10.218 7.179 10.112 102.348 247.560 2000 16.643 7.087 9.947 98.988 248.515 2500 28.760 7.032 10.057 101.158 276.705 3000 41.467 6.986 9.690 93.922 231.058 3500 61.748 6.843 9.648 93.231 236.175 4000 92.689 6.903 9.594 92.119 239.006 4500 139.815 7.778 10.057 101.152 283.829 5000 179.735 7.198 9.677 93.660 247.062 Table 11: Full YearPredictionMSD: metrics normalized by the training-target scale from Table 10. r sd(ytrain)sd(y_train) MAEz RMSEz σ^ε2(z) σ_ ^2(z) σ^ε2(orig) σ_ ^2(orig) 1000 10.936 0.681 0.938 4.414865×10−4\,4.414865× 10^-4\, 5.280451×10−2\,5.280451× 10^-2\, 1500 10.936 0.656 0.925 4.415726×10−4\,4.415726× 10^-4\, 5.281479×10−2\,5.281479× 10^-2\, 2000 10.936 0.648 0.910 4.416462×10−4\,4.416462× 10^-4\, 5.281858×10−2\,5.281858× 10^-2\, 2500 10.936 0.643 0.920 4.415373×10−4\,4.415373× 10^-4\, 5.281404×10−2\,5.281404× 10^-2\, 3000 10.936 0.639 0.886 4.416178×10−4\,4.416178× 10^-4\, 5.282333×10−2\,5.282333× 10^-2\, 3500 10.936 0.626 0.882 4.416175×10−4\,4.416175× 10^-4\, 5.282105×10−2\,5.282105× 10^-2\, 4000 10.936 0.631 0.877 4.416362×10−4\,4.416362× 10^-4\, 5.282208×10−2\,5.282208× 10^-2\, 4500 10.936 0.711 0.920 4.414865×10−4\,4.414865× 10^-4\, 5.280713×10−2\,5.280713× 10^-2\, 5000 10.936 0.658 0.885 4.415092×10−4\,4.415092× 10^-4\, 5.280737×10−2\,5.280737× 10^-2\, Tables 9–10 highlight two main findings. First, the MAP estimates of the KmixK_mix parameters are stable as r varies. In particular, α^≈0.487 α≈ 0.487 and w^≈0.512 w≈ 0.512 remain nearly constant, indicating that the model preserves a balanced mixture of the components and does not collapse to the degenerate cases (w∈0,1w∈\0,1\). Similarly, the estimated nugget remains essentially constant across the grid of r, showing variation only in the 5th–6th decimal place when reported with higher precision. On the full YearPredictionMSD dataset, σ^ε2(orig) σ_ ^2(orig) stays within [5.280451, 5.282333]×10−2[5.280451,\,5.282333]× 10^-2, which corresponds to a relative variation of approximately 0.036%0.036\% (Table 11). This stability is consistent with interpreting σε2 _ ^2 as i.i.d. noise: as r changes, it primarily affects the Nyström approximation’s ability to represent the correlated component of the latent process, while the estimated observational noise level remains essentially unchanged. Residual differences can be attributed to numerical/optimization effects and rounding. In addition to assessing the stability of the fit as r varies, Table 10 can be interpreted operationally as a direct guide for choosing the Nyström rank under a total time budget T (training + prediction). For a fixed budget, a simple rule is to select rmax(T)r_ (T), defined as the largest rank in the grid such that time(r)≤Ttime(r)≤ T, which makes explicit the cost of increasing the capacity of the approximator. For instance, under T=30T=30s, the largest feasible rank in the grid is r=2500r=2500 (28.760s), with MAE =7.032=7.032 and RMSE =10.057=10.057; under T=60T=60s, rmax(T)=3000r_ (T)=3000 (41.467s), with MAE =6.986=6.986 and RMSE =9.690=9.690; and under T=300T=300s, rmax(T)=5000r_ (T)=5000 (179.735s), with MAE =7.198=7.198 and RMSE =9.677=9.677. At the same time, the error is not strictly monotone in r: some intermediate values perform worse than smaller ranks (e.g., r=2500r=2500 compared with r=2000r=2000), and there is also a point with a marked degradation (r=4500r=4500). This behavior is consistent with numerical variability and sensitivity of training/optimization when combined with a low-rank approximation, and it reinforces the trade-off interpretation with diminishing returns: increasing r raises cost, but does not guarantee an immediate improvement in MAE/RMSE. Therefore, a second practical rule is to choose, within each budget T, the rank rbest(T)r_best(T) that minimizes the test RMSE. Table 12 summarizes both choices. In particular, for T=60T=60s, r=3000r=3000 emerges as a competitive point by combining low cost with reduced RMSE; for a more generous budget (T=300T=300s), the lowest RMSE in the grid occurs at r=4000r=4000 (92.689s), still well below the cost of r=5000r=5000, indicating that additional gains beyond r∈[3000,4000]r∈[3000,4000] tend to be marginal relative to the increase in runtime. Table 12: Choice of r on the full YearPredictionMSD dataset under total runtime budgets T. We report rmax(T)r_ (T) (the largest r such that time(r)≤Ttime(r)≤ T) and, additionally, the r that yields the lowest RMSE within the same budget. Roof T rmax(T)r_ (T) Total time (s) MAE RMSE rbest(T)r_best(T) RMSEbest_best 30s 2500 28.760 7.032 10.057 2000 9.947 60s 3000 41.467 6.986 9.690 3000 9.690 300s 5000 179.735 7.198 9.677 4000 9.594 In summary, the full-dataset study confirms the Nyström rank as the main driver of the cost–accuracy trade-off: larger r values tend to reduce error up to a saturation regime, while the kernel parameters and the nugget remain essentially stable. Considering accuracy and computational cost jointly, r=4000r=4000 yields the lowest RMSE among the evaluated values, although configurations such as r∈[3000,4000]r∈[3000,4000] are also competitive at a substantially lower cost than r=5000r=5000. Overall, on the real datasets, MAP estimates remain consistent across sample sizes, and performance confirms that the scalable approximation is the primary source of the time–error trade-off. In the large-scale study, the approximation complexity directly controls this trade-off, with gains that tend to saturate beyond an intermediate range and improvements concentrated in intermediate regimes. 6 Conclusion This work strengthens the connection between infinite-width limits of BNNs and GPs, and leverages this bridge to propose and evaluate a mixed kernel KmixK_mix, interpreted as a combination of two complementary process behaviors. From a theoretical standpoint, the formulation establishes conditions that guarantee fundamental properties of the proposed class and discusses regimes of non-identifiability in finite samples, motivating a reduced parameterization that is more stable and interpretable. From a computational standpoint, we develop a scalable procedure based on MAP estimation with a Nyström approximation, enabling efficient training and prediction with explicit uncertainty quantification. Empirically, simulation studies and experiments on real-world datasets support three main findings: (i) MAP fitting is stable across scenarios, with the mixing parameter w consistently estimated and showing no systematic collapse to degenerate cases; (i) anchor selection in the Nyström approximation affects the cost–accuracy trade-off, with more noticeable gains in more complex regimes; and (i) at scale, the approximation rank directly governs the time–performance trade-off, with diminishing returns beyond an intermediate range. These results indicate that the mixing in KmixK_mix provides practical robustness by avoiding losses associated with fixed kernel choices, while retaining a small number of parameters with clear probabilistic interpretation. Despite these advances, the study has limitations. First, the experimental evaluation was conducted over a finite set of scenarios and configuration choices, so other regimes may exhibit different behavior. Second, we focus on univariate regression with i.i.d. Gaussian noise, which limits immediate applicability to multivariate, heteroscedastic, or classification settings. Third, while Nyström provides scalability, performance depends both on anchor selection and on the complexity of the low-rank approximation, and a systematic sensitivity analysis of these factors was not explored exhaustively. As future work, we aim to extend the model to multivariate outputs and non-Gaussian likelihoods, including classification. We also plan to investigate hierarchical variants inspired by the BNN → GP bridge to capture local heterogeneities, broaden empirical validation across additional domains and scaling regimes with more systematic sensitivity analyses, and compare alternative scalable approximations and inference strategies beyond MAP, with the goal of improving calibration and robustness under controlled computational budgets. Acknowledgments We thank PPGEst–UFMG for a collaborative research environment and UEMG for institutional support. This work was supported by CNPq, FAPEMIG, and CAPES. Appendix A Technical details and proofs A.1 Pre-activation moments and dependence across inputs For a fixed hidden unit j, consider the pre-activations z()=a+⊤,z(′)=a+⊤′,z(x)=a+u x,z(x )=a+u x , where a and u are drawn from the parameter priors (with zero mean and finite variances). Under the assumptions in the main text, we have: Mean. [z()]=[a]+∑i=1Ixi[ui]=0,[z(′)]=0.E[z(x)]=E[a]+ _i=1^Ix_iE[u_i]=0, [z(x )]=0. Variance. Assuming independence between a and u, and finite variances, Var(z())=Var(a)+Var(⊤)=σa2+σu2⊤,Var (z(x) )=Var(a)+Var(u x)= _a^2+ _u^2\,x x, and similarly Var(z(′))=σa2+σu2′⁣⊤′Var (z(x ) )= _a^2+ _u^2\,x x . Covariance. Note that z()z(x) and z(′)z(x ) are not independent: they share the same a and the same u. Their covariance is Cov(z(),z(′))=Var(a)+Cov(⊤,⊤′)=σa2+σu2⊤′.Cov (z(x),z(x ) )=Var(a)+Cov(u x,u x )= _a^2+ _u^2\,x x . In particular, when a and u are Gaussian, (z(),z(′))(z(x),z(x )) is a centered bivariate Gaussian with covariance matrix determined by the quantities above. A.2 A general regularity condition We next state a simple sufficient condition ensuring [h(Z)2]<∞E[h(Z)^2]<∞ when Z is Gaussian. Proposition A.1 (Polynomial growth implies a finite second moment). If h:ℝ→ℝh:R satisfies |h(z)|≤C(1+|z|p)|h(z)|≤ C(1+|z|^p) for some C>0C>0 and p≥0p≥ 0, then for Z∼(0,σ2)Z (0,σ^2) we have [h(Z)2]<∞E[h(Z)^2]<∞. Proof. By assumption, h(Z)2≤2C2(1+|Z|2p)h(Z)^2≤ 2C^2(1+|Z|^2p). Since Gaussian random variables have moments of all orders, [|Z|2p]<∞E[|Z|^2p]<∞. Hence, [h(Z)2]≤2C2(1+[|Z|2p])<∞E[h(Z)^2]≤ 2C^2(1+E[|Z|^2p])<∞. ∎ A.3 Proof of Theorem 2.1 For a fixed k, the function fk()=bk+∑j=1Hvkjh(zj())f_k(x)=b_k+ _j=1^Hv_kj\,h\! (z_j(x) ) is the sum of the bias term bkb_k and H terms of the form vkjh(zj())v_kjh(z_j(x)). Under the the prior independence assumptions, the terms vkjh(zj())v_kjh(z_j(x)) are independent across j. Define vkj∗=Hvkjv_kj^*= H\,v_kj, so that vkj∗∼(0,σv2)v_kj^* (0, _v^2). Then, by the CLT, H∑j=1Hvkj∗h(zj())H⟹(0,σv2[h(zj())2]). H _j=1^H v_kj^*\,h(z_j(x))H\; \;N\! (0,\, _v^2\,E[h(z_j(x))^2] ). Independently of the specific distributions of aja_j, ujiu_ji, and vkjv_kj, as long as each term vkjh(zj())v_kjh(z_j(x)) has mean zero and finite variance, the sum ∑j=1Hvkjh(zj()) _j=1^Hv_kjh(z_j(x)) remains a sum of H independent terms, due to the independence of vkjv_kj and zjz_j. Thus, for each x, [fk()]=[bk]+∑j=1H[vkj][h(zj())] [f_k(x)]=E[b_k]+ _j=1^HE[v_kj]\,E[h(z_j(x))] =0+∑j=1H0⋅[h(zj())]=0, =0+ _j=1^H0·E[h(z_j(x))]=0, and Var[fk()]=Var[bk]+∑j=1HVar(vkjh(zj()))=σb2+∑j=1Hσv2H[h(zj())2]=σb2+σv2[h(zj())2]. [f_k(x)]=Var[b_k]+ _j=1^HVar\! (v_kjh(z_j(x)) )= _b^2+ _j=1^H _v^2H\,E[h(z_j(x))^2]= _b^2+ _v^2\,E[h(z_j(x))^2]. Therefore, fk()⇒(0,σb2+σv2[h(zj())2]).f_k(x)\; d\;N\! (0,\, _b^2+ _v^2\,E[h(z_j(x))^2] ). Note that one must verify that second-order moments are finite; if h satisfies Proposition A.1, then the theorem’s assumptions guarantee that Var[fk()]<∞Var[f_k(x)]<∞. Now consider n input points (1,2,…,n)(x_1,x_2,…,x_n). For each p∈1,…,np∈\1,…,n\, fk(p)f_k(x_p) is a sum of terms vkjh(zj(p))v_kjh(z_j(x_p)) plus a bias, so fk(p)∼(0,σb2+σv2[h(zj(p))2])f_k(x_p) \! (0, _b^2+ _v^2E[h(z_j(x_p))^2] ), i.e., ∼n(,)f _n(0,K). Note that for different px_p and p′x_p , the terms vkjh(zj(p))v_kjh(z_j(x_p)) and vkjh(zj(p′))v_kjh(z_j(x_p )) are correlated because they share the same vkjv_kj. Moreover, since vkjv_kj and h(zj(p′))h(z_j(x_p )) are independent, Cov[fk(p),fk(p′)] [f_k(x_p),f_k(x_p )] =[fk(p)fk(p′)]−[fk(p)][fk(p′)] =E[f_k(x_p)f_k(x_p )]-E[f_k(x_p)]E[f_k(x_p )] =[fk(p)fk(p′)]−0 =E[f_k(x_p)f_k(x_p )]-0 =Var[bk]+∑j=1HCov(vkjh(zj(p)),vkjh(zj(p′))). =Var[b_k]+ _j=1^HCov\! (v_kjh(z_j(x_p)),\,v_kjh(z_j(x_p )) ). Since vkjv_kj are independent across different j, ∑j=1HCov(vkjh(zj(p)),vkjh(zj(p′))) _j=1^HCov\! (v_kjh(z_j(x_p)),v_kjh(z_j(x_p )) ) =∑j=1H[vkj2][h(zj(p))h(zj(p′))] = _j=1^HE[v_kj^2]\,E[h(z_j(x_p))h(z_j(x_p ))] =∑j=1Hσv2H[h(zj(p))h(zj(p′))] = _j=1^H _v^2H\,E[h(z_j(x_p))h(z_j(x_p ))] =σv2[h(z(p))h(z(p′))]. = _v^2\,E[h(z(x_p))h(z(x_p ))]. (26) Therefore, Cov[fk(p),fk(p′)] [f_k(x_p),f_k(x_p )] =σb2+σv2[h(z(p))h(z(p′))] = _b^2+ _v^2\,E[h(z(x_p))h(z(x_p ))] =σb2+σv2∫h(z(p))h(z(p′))⋅p(zj(p),zj(p′))zz′ = _b^2+ _v^2 h(z(x_p))\,h(z(x_p ))· p\! (z_j(x_p),z_j(x_p ) )\,dz\,dz =σb2+σv2C(p,p′) = _b^2+ _v^2\,C(x_p,x_p ) =K(p,p′). =K(x_p,x_p ). The fact that, for any finite set of points 1,2,…,n\x_1,x_2,…,x_n\, the vector (fk(1),…,fk(n))(f_k(x_1),…,f_k(x_n)) converges to a multivariate Normal distribution implies that fkf_k converges to a GP with zero mean and the covariance defined above, namely =σb2 1nn⊤+σv2.K\;=\; _b^2\,1_n1_n \;+\; _v^2\,C. A.4 Proof of Proposition 3.1 Step 0 (setup and independence). We work under the additive (block) mixture construction described in the main text. For each component m=1,…,Mm=1,…,M, the network uses a separate collection of output weights (vj(m))j=1H(v^(m)_j)_j=1^H, and we assume: vj(m)are mutually independent in j,independent across m,vj(m)∼iid(0,σv2H).v^(m)_j\ are mutually independent in j,\ independent across m, v^(m)_j iid N\! (0, _v^2H ). The pre-activations are zj()=aj+j⊤,j=1,…,H,z_j(x)=a_j+u_j x, j=1,…,H, where (aj,j)(a_j,u_j) are i.i.d. in j and independent of all vj(m)v^(m)_j. We also assume b∼(0,σb2)b (0, _b^2) independent of all other parameters. Step 1 (network output). For a single-output network, the block-mixture prior is f()=b+∑m=1Mwm∑j=1Hvj(m)hm(zj()),wm∈(0,1),∑m=1Mwm=1.f(x)=b+ _m=1^M w_m _j=1^Hv^(m)_j\,h_m\! (z_j(x) ), w_m∈(0,1),\ \ _m=1^Mw_m=1. (For multiple outputs k, the argument applies componentwise to fkf_k.) Step 2 (expand the covariance). For two inputs i,i′x_i,x_i , by bilinearity of covariance, Cov[f(i),f(i′)] \! [f(x_i),f(x_i ) ] =Var(b) =Var(b) +Cov(∑m=1Mwm∑j=1Hvj(m)hm(zj(i)),∑m′=1Mwm′∑ℓ=1Hvℓ(m′)hm′(zℓ(i′))). +Cov\! ( _m=1^M w_m _j=1^Hv^(m)_jh_m(z_j(x_i)),\; _m =1^M w_m _ =1^Hv^(m )_ h_m (z_ (x_i )) ). Step 3 (cross terms vanish). All terms with m≠m′m≠ m vanish because the blocks vj(m)\v^(m)_j\ and vℓ(m′)\v^(m )_ \ are independent and centered. Likewise, for fixed m, all terms with j≠ℓj≠ vanish because vj(m)v^(m)_j are independent across j and centered. Therefore, Cov[f(i),f(i′)]=σb2+∑m=1Mwm∑j=1HCov(vj(m)hm(zj(i)),vj(m)hm(zj(i′))).Cov\! [f(x_i),f(x_i ) ]= _b^2+ _m=1^Mw_m _j=1^HCov\! (v^(m)_jh_m(z_j(x_i)),\,v^(m)_jh_m(z_j(x_i )) ). Step 4 (factor each remaining covariance). Since vj(m)v^(m)_j is independent of zj(⋅)z_j(·) and [vj(m)]=0E[v^(m)_j]=0, Cov(vj(m)hm(zj(i)),vj(m)hm(zj(i′)))=[(vj(m))2][hm(zj(i))hm(zj(i′))].Cov\! (v^(m)_jh_m(z_j(x_i)),\,v^(m)_jh_m(z_j(x_i )) )=E[(v^(m)_j)^2]\,E\! [h_m(z_j(x_i))\,h_m(z_j(x_i )) ]. With the scaling [(vj(m))2]=σv2/HE[(v^(m)_j)^2]= _v^2/H, we obtain Cov[f(i),f(i′)]=σb2+∑m=1Mwm∑j=1Hσv2H[hm(zj(i))hm(zj(i′))].Cov\! [f(x_i),f(x_i ) ]= _b^2+ _m=1^Mw_m _j=1^H _v^2H\,E\! [h_m(z_j(x_i))\,h_m(z_j(x_i )) ]. Step 5 (simplify the sum over hidden units). Because (aj,j)(a_j,u_j) are i.i.d. in j, the pair (zj(i),zj(i′)) (z_j(x_i),z_j(x_i ) ) has the same distribution for every j, so [hm(zj(i))hm(zj(i′))]=[hm(Z)hm(Z′)]=:Km(i,i′),E\! [h_m(z_j(x_i))\,h_m(z_j(x_i )) ]=E\! [h_m(Z)\,h_m(Z ) ]=:K_m(x_i,x_i ), where (Z,Z′)(Z,Z ) denotes the pre-activation pair induced by the priors at inputs (i,i′)(x_i,x_i ). Thus, ∑j=1Hσv2H[hm(zj(i))hm(zj(i′))]=H⋅σv2HKm(i,i′)=σv2Km(i,i′). _j=1^H _v^2H\,E\! [h_m(z_j(x_i))\,h_m(z_j(x_i )) ]=H· _v^2H\,K_m(x_i,x_i )= _v^2\,K_m(x_i,x_i ). Substituting back yields the kernel form Cov[f(i),f(i′)]=σb2+σv2∑m=1MwmKm(i,i′).Cov\! [f(x_i),f(x_i ) ]= _b^2+ _v^2 _m=1^Mw_m\,K_m(x_i,x_i ). Step 6 (GP limit and positive semidefiniteness). For any finite set of inputs, f()f(x) is a sum (over j) of independent terms with variance controlled by the 1/H1/H scaling. Under the same moment assumptions as in Theorem 2.1 (e.g., [hm(Z)2]<∞E[h_m(Z)^2]<∞), the multivariate CLT implies that the finite-dimensional distributions converge to a multivariate Gaussian with covariance matrix built from K(,′)=σb2+σv2∑m=1MwmKm(,′).K(x,x )= _b^2+ _v^2 _m=1^Mw_m\,K_m(x,x ). Hence, as H→∞H→∞, f(⋅)f(·) converges (in the sense of finite-dimensional distributions) to a GP with covariance K. Finally, since each KmK_m is symmetric positive semidefinite and wm≥0w_m≥ 0, the kernel K is also symmetric positive semidefinite. A.5 Proof of Proposition 3.2 Proof. Suppose there exist θ1≠θ2 _1≠ _2 in Θ such that Kmix(,′;θ1)≡Kmix(,′;θ2),∀(,′)∈ℝd×ℝd.K_mix(x,x ; _1)≡ K_mix(x,x ; _2), ∀\,(x,x ) ^d×R^d. Write Kmix(,′;θ)=σb2+σv2[wA(,′)+(1−w)B(,′;α)],K_mix(x,x ;θ)= _b^2+ _v^2 [w\,A(x,x )+(1-w)\,B(x,x ;α) ], where A denotes the smooth (arcsine) term and B(⋅;α)B(·;α) denotes the angular (ReLU/LeakyReLU) term, both depending on (,′)(x,x ) through ρ(,′)ρ(x,x ) and through σz(),σz′(′) _z(x), _z (x ), which are determined by (σa2,σu2)( _a^2, _u^2). Step 1: identification of σb2 _b^2. From the identity Kmix(⋅,⋅;θ1)≡Kmix(⋅,⋅;θ2)K_mix(·,·; _1)≡ K_mix(·,·; _2) it follows that σb,12−σb,22=Fθ2(,′)−Fθ1(,′),∀(,′), _b,1^2- _b,2^2=F_ _2(x,x )-F_ _1(x,x ), ∀(x,x ), where Fθ:=σv2[wA+(1−w)B(⋅;α)]F_θ:= _v^2 [wA+(1-w)B(·;α) ]. The left-hand side is constant. Fixing ‖=‖′‖=r\|x\|=\|x \|=r and choosing different angles between x and ′x (which yields different values of ρ), both A and B vary; hence the difference on the right-hand side cannot be constant over the whole domain unless Fθ1≡Fθ2F_ _1≡ F_ _2. Consequently, σb,12=σb,22 _b,1^2= _b,2^2 and Fθ1(,′)≡Fθ2(,′).F_ _1(x,x )≡ F_ _2(x,x ). Step 2: non-proportionality of A and B(⋅;α)B(·;α). Restrict attention to pairs with ‖=‖′‖=r\|x\|=\|x \|=r, so that σz=σz′ _z= _z is fixed and ρ varies in (−1,1)(-1,1). Within this family, A(ρ)A(ρ) is real-analytic on (−1,1)(-1,1) (since it involves arcsin ), whereas B(ρ;α)B(ρ;α) contains terms such as 1−ρ2 1-ρ^2 and arccos⁡(ρ) (ρ), with qualitatively different behavior as |ρ|→1|ρ|→ 1. Therefore, A and B(⋅;α)B(·;α) cannot be proportional as functions of ρ over the entire interval (−1,1)(-1,1). Step 3: identification of (σa2,σu2)( _a^2, _u^2). Now consider orthogonal pairs ⟂′x with ‖=‖′‖=r>0\|x\|=\|x \|=r>0. Then ρ(r)=σa2σa2+σu2r2,σz2(r)=σz′2(r)=σa2+σu2r2.ρ(r)= _a^2 _a^2+ _u^2r^2, _z^2(r)= _z ^2(r)= _a^2+ _u^2r^2. The identity Fθ1≡Fθ2F_ _1≡ F_ _2 for all r implies that the functions r↦ρ(r)r ρ(r) and r↦σz2(r)r _z^2(r) coincide under θ1 _1 and θ2 _2. Since σz2(r) _z^2(r) is affine in r2r^2 and ρ(r)ρ(r) depends only on the pair (σa2,σu2)( _a^2, _u^2) through a rational function of r2r^2, it follows that (σa,12,σu,12)=(σa,22,σu,22).( _a,1^2, _u,1^2)=( _a,2^2, _u,2^2). Step 4: identification of σv2,w,α _v^2,w,α. With (σa2,σu2)( _a^2, _u^2) identified, the quantities ρ(,′)ρ(x,x ), σz() _z(x), and σz′(′) _z (x ) match on both sides. Since A and B(⋅;α)B(·;α) are not proportional in ρ, the equality Fθ1≡Fθ2F_ _1≡ F_ _2 forces the coefficients to coincide: σv,12=σv,22,w1=w2,α1=α2. _v,1^2= _v,2^2, w_1=w_2, _1= _2. Together with σb,12=σb,22 _b,1^2= _b,2^2, we obtain θ1=θ2 _1= _2, contradicting θ1≠θ2 _1≠ _2. Hence, the mapping θ↦Kmix(⋅,⋅;θ)θ K_mix(·,·;θ) is injective in the interior of Θ . In boundary cases (w∈0,1w∈\0,1\ or α∈0,1α∈\0,1\), some parameters may drop out of the kernel and identifiability may fail; these cases are excluded by assumption. ∎

Instant research discovery