1 Introduction↩︎

We consider predictive regression models in settings where reliable uncertainty estimates are required alongside (i) computational efficiency, (ii) strong predictive performance, and (iii) interpretability. This requirement is particularly important in safety-critical and risk-sensitive applications [1], where deep neural networks (NNs) are known to produce poorly calibrated predictions that can be overconfident [2]–[4], especially under distribution shift. Moreover, regression models with quantified uncertainty serve as core components in broader learning systems, including reinforcement learning [5]–[8] and physics-informed neural networks [9]. Goals (i)–(iii) are inherently in tension. In what follows, we examine several widely used paradigms that address these requirements to varying degrees while exposing important trade-offs.

In terms of predictive performance and principled uncertainty estimation, Gaussian Processes (GPs) provide a canonical Bayesian framework for nonlinear regression with exact posterior inference [10]. However, their cubic computational complexity and sensitivity to kernel choice fundamentally limit scalability, particularly in high-dimensional settings where distance-based kernels suffer from the curse of dimensionality. Deep Kernel Learning (DKL) partially mitigates this limitation by learning data-dependent covariance functions using neural networks, but exact inference remains computationally intensive, and scalable implementations typically rely on approximations [11].

Fully Bayesian neural networks (BNNs) represent a different point in the performance–computation trade-off by placing probability distributions over network parameters [12], [13]. While they provide a flexible probabilistic framework for modeling uncertainty, exact inference is generally intractable, and practical implementations rely on approximate inference methods such as mean-field variational inference [14]. These approximations can impose restrictive assumptions on the posterior distribution, which may lead to miscalibrated uncertainty estimates [15], [16].

Bayesian last-layer (BLL) methods address the computational limitations of fully Bayesian neural networks by restricting Bayesian inference to the final linear layer of a neural network. This design enables efficient inference while retaining a probabilistic interpretation of predictions [17]. However, because uncertainty is placed only on a single learned representation, the representation itself is treated as fixed. As a result, uncertainty estimates can become overconfident when the learned representation fails to capture variability in the data, particularly outside the training distribution [4]. Extensions such as functional regularization have been proposed to mitigate this limitation, but they introduce additional computational overhead without fundamentally addressing the reliance on a single representation [4].

Deep ensembles (DE) provide a simple and scalable approach to uncertainty quantification by training multiple neural networks independently and uniformly averaging their predictive distributions [3]. These methods typically yield strong empirical performance while keeping computational costs moderate. However, they lack an explicit probabilistic interpretation, as uncertainty is quantified implicitly through model diversity rather than through a well-defined probability measure over functions. Moreover, uniform averaging is non-diagnostic: it does not provide probabilistic evidence for which ensemble members are supported by the data or how individual members contribute to predictions. As a result, performance can degrade when member quality varies, increasing the risk of overconfident predictions.

Taken together, these approaches illustrate a recurring pattern: methods that provide principled uncertainty estimation often incur substantial computational cost or restrict model flexibility, while scalable approaches typically rely on heuristic aggregation rules that offer limited interpretability. Designing models that simultaneously achieve computational efficiency, strong predictive performance, and interpretable uncertainty estimation, therefore, remains an open challenge.

1.1 Main contributions↩︎

In this work, we propose Bayesian Deep Kernel Networks, a regression framework that combines independently trained neural networks with closed-form Bayesian aggregation. The method trains \(H\) neural networks on a dataset \(\mathcal{D} = \{(\mathbf{x}_i,y_i)\}_{i=1}^N\), producing multiple learned feature representations that serve as basis functions. A Bayesian linear regression layer is then applied to these ensemble outputs, yielding an analytic posterior over aggregation weights and a closed-form posterior predictive distribution (Figure 1). This design enables principled uncertainty estimation while retaining the computational efficiency and predictive performance of standard deep ensembles.

The framework leverages the scalability of deep ensembles while replacing uniform averaging with Bayesian aggregation. By performing exact inference at the aggregation stage, the method produces interpretable posterior weights and calibrated uncertainty estimates without relying on approximate inference. Taken together, this design addresses the competing requirements (i)–(iii). We summarize the resulting contributions below.

Scalable Bayesian aggregation with learned neural predictors and finite-rank GP interpretation. We aggregate independently trained neural networks using closed-form Bayesian linear regression, yielding an analytic posterior predictive distribution with computational cost comparable to standard deep ensembles. The resulting model can be interpreted as a finite-rank Gaussian Process whose basis functions are learned directly through ensemble training, enabling principled Bayesian aggregation without kernel design or cubic-cost GP inference.

Separation of representation learning and probabilistic inference. Our framework maintains the standard deep-ensemble training pipeline and performs Bayesian inference only at the aggregation stage. Unlike Bayesian last-layer (BLL) methods that place uncertainty on a single latent representation, our method defines uncertainty directly over ensemble-member contributions. This separation maintains implementation simplicity while providing interpretable probabilistic weights for individual predictors.

Posterior-driven aggregation robust to heterogeneous ensemble member quality. Instead of uniform averaging, the aggregation weights are inferred from the posterior distribution. Ensemble members weakly supported by the data are automatically downweighted, improving robustness when model quality varies due to optimization variability or training conditions. This robustness is achieved without manual model selection or additional tuning heuristics.

To evaluate the practical implications of these design choices and their impact on computational efficiency, predictive performance, and interpretability, we conduct experiments on standard UCI regression benchmarks across a range of architectural and optimization settings.

2 Bayesian deep kernel networks↩︎

In this section, we introduce our framework while highlighting each improvement in comparison to existing methods. We start by introducing BLL and proceed to build upon it.

2.1 Learning the feature space representation by ensemble members↩︎

Let \(\mathcal{D} = \{(\boldsymbol{x}_i, y_i)\}_{i=1}^N\) denote the dataset, where \(\boldsymbol{x}_i \in \mathbb{R}^P\) are input vectors and \(y_i \in \mathbb{R}\) are scalar responses. We consider scalar responses throughout this work for simplicity and without loss of generality. Denote by \(\phi(\cdot;\boldsymbol{\theta}): \mathbb{R}^{P} \to \mathbb{R}^H\) a feature map with parameter \(\boldsymbol{\theta}.\) In BLLs the latent function \(f\) is modeled as \[\begin{align} y_i = f(\boldsymbol{x}_i;\boldsymbol{\theta})=\boldsymbol{\phi}(\boldsymbol{x}_i;\boldsymbol{\theta})^\top \boldsymbol{\beta}+ \varepsilon_i ,\quad \varepsilon_i \sim \mathcal{N}(0, \sigma^2_{\varepsilon}) \end{align}\] where \(\boldsymbol{\beta}\sim\mathcal{N}(\mu_{\beta},\Sigma_{\beta})\) denotes the feature weights.

Stacking the responses in the vector \(\mathbb{R}^N\) yields the compact representation \[\begin{align} \label{eq:lin95reg} \mathbf{y} = \boldsymbol{\Phi}(\mathbf{x};\boldsymbol{\theta})\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \end{align}\tag{1}\] where \(\boldsymbol{\Phi}(\mathbf{x};\boldsymbol{\theta}) = [\boldsymbol{\phi}(\mathbf{x}_1;\boldsymbol{\theta}),\ldots,\boldsymbol{\phi}(\mathbf{x}_N;\boldsymbol{\theta})\boldsymbol{]}^{\top} \in \mathbb{R}^{N\times H}.\) Then, training is performed by finding \(\boldsymbol{\theta},\,\sigma_{\varepsilon}^2,\) and \(\boldsymbol{\Sigma}_{\beta}\) which requires maximizing an exact or approximate marginal likelihood and \(\boldsymbol{\beta}\) is obtained through its posterior \(p(\boldsymbol{\beta}\mid \mathcal{D}, \theta)\) in closed form [4]. Observe that the feature map \(\phi(\mathbf{x};\boldsymbol{\theta})\) is learned once, and then treated as fixed during Bayesian inference. This implies that the variability employed to capture uncertainty stems only from the posterior over the weights \(\boldsymbol{\beta}\) with the representation itself remaining deterministic. Therefore, uncertainty is restricted only to linear combinations within a single learned feature space.

In order to improve uncertainty behavior in out-of-distribution regions, we combine ideas from DE-style representation with a BLL to obtain better calibrated uncertainty. Unlike standard BLL, which trains a single neural network to construct \(\boldsymbol{\Phi}.\) Thus, we define a collection of feature maps \(\{\phi_h(\cdot;\boldsymbol{\theta}_h)\}_{h=1}^H,\) where each element is an individual map parameterized by its own \(\boldsymbol{\theta}_h,\) yielding \[\begin{align} \boldsymbol{\Phi}(\mathbf{x};\boldsymbol{\theta}) = [\phi_1(\mathbf{x}_1;\boldsymbol{\theta}_1),\ldots, \phi_h(\mathbf{x}_h;\boldsymbol{\theta}_h)]^{\top} \end{align}\] in 1 . Each feature map of the \(H\) feature maps is implemented by an independently trained neural network. Training is performed by maximizing the Gaussian likelihood of \(y_i \sim \mathcal{N}(\phi_h(\mathbf{x}_i;\theta_h),\sigma^2_h(\mathbf{x}_i;\theta_h))\)

where both the mean function \(\phi_h(\cdot)\) and the variance function \(\sigma_h^2(\cdot)\) are outputs of the network. Note that \(\sigma^2_h(\cdot)\) is used only during training to induce a heteroscedastic objective with the goal of preventing all networks from fitting the data with the same smoothness profile. It reweights the contribution of each training point to the loss: points with larger predicted variance incur a smaller penalty, while those with smaller variance are fitted more tightly. This sample-dependent reweighting leads to heterogeneous learned feature maps across independently trained ensemble members.

The independent training yields the following benefits for our model. First, feature maps can be trained fully in parallel, allowing the ensemble size to scale without increasing training time beyond that of a single neural network. Second, each feature map represents a complete predictive model of the data, rather than a latent basis. In this way, the collection of independently learned feature maps captures multiple distinct functional representations of the data, which is particularly relevant when data are scarce compared to model size.

2.2 Bayesian inference in the induced function space↩︎

Once \(\boldsymbol{\Phi}\) has been obtained, the posterior distribution over \(\boldsymbol{\beta} \in \mathbb{R}^H\) is inferred via Bayesian linear regression in the function space 1 .

We place a conjugate prior over the linear coefficients \(\boldsymbol{\beta}\) and the observation-noise variance \(\sigma^2_{\varepsilon}\). We adopt Jeffreys’ prior for the noise variance to improve generalizability and uncertainty quantification, \[p(\sigma^2_{\varepsilon}) \propto \frac{1}{\sigma^2_{\varepsilon}},\] which is invariant under reparameterization and introduces no additional scale hyperparameters (See Appendix 6.1 for derivation). Conditioned on \(\sigma^2_{\varepsilon}\), the regression coefficients are assigned a Gaussian prior, \[\boldsymbol{\beta} \mid \sigma^2_{\varepsilon} \sim \mathcal{N}\!\big(\boldsymbol{\mu}_0,\; \sigma^2_{\varepsilon} \boldsymbol{\Lambda}_0^{-1}\big),\] where \(\boldsymbol{\Lambda}_0 \in \mathbb{R}^{H \times H}\) is positive definite.

This prior specification can be viewed as a special case of the Normal-Inverse-Gamma prior, where \(\sigma^2_{\varepsilon} \sim \mathrm{Inv\text{-}Gamma}(a_0,b_0)\), obtained in the limit \(a_0 \to 0\) and \(b_0 \to 0\). Jeffreys’ prior is a common choice when no prior information on variance is provided [18], [19], while maintaining analytical tractability and eliminating the need to estimate additional hyperparameters for the variance.

Posterior updates. Given the design matrix \(\boldsymbol{\Phi}\) and targets \(\mathbf{y}\), conjugacy is maintained and the posterior factorizes as \(p(\boldsymbol{\beta},\sigma^2_{\varepsilon}\mid \mathbf{y},\boldsymbol{\Phi}) = p(\boldsymbol{\beta}\mid\sigma^2_{\varepsilon},\mathbf{y},\boldsymbol{\Phi})\; p(\sigma^2_{\varepsilon}\mid \mathbf{y},\boldsymbol{\Phi}).\) The conditional posterior over coefficients is Gaussian, \(\boldsymbol{\beta}\mid\sigma^2_{\varepsilon},\mathbf{y},\boldsymbol{\Phi} \sim \mathcal{N}\!\big(\boldsymbol{\mu}_n,\; \sigma^2_{\varepsilon}\boldsymbol{\Lambda}_n^{-1}\big),\) with \(\boldsymbol{\mu}_n = \boldsymbol{\Lambda}_n^{-1} \big(\boldsymbol{\Lambda}_0\boldsymbol{\mu}_0 + \boldsymbol{\Phi}^\top \mathbf{y}\big)\) and \(\boldsymbol{\Lambda}_n = \boldsymbol{\Lambda}_0 + \boldsymbol{\Phi}^\top\boldsymbol{\Phi}.\) Under the Jeffreys prior on \(\sigma^2_{\varepsilon}\), the marginal posterior over \(\sigma^2_{\varepsilon}\) is \(\sigma^2_{\varepsilon}\mid \mathbf{y},\boldsymbol{\Phi}\sim \mathrm{Inv\text{-}Gamma}(a_n,b_n),\) where \(a_n = \frac{N}{2}\) and \(b_n = \frac{1}{2}\!\left( \mathbf{y}^\top\mathbf{y} + \boldsymbol{\mu}_0^\top\boldsymbol{\Lambda}_0\boldsymbol{\mu}_0 - \boldsymbol{\mu}_n^\top\boldsymbol{\Lambda}_n\boldsymbol{\mu}_n \right).\)

Posterior predictive. For a new input with feature vector \(\boldsymbol{\phi}_* = \boldsymbol{\phi}(\mathbf{x}_\star)\), the predictive distribution integrates out both \(\boldsymbol{\beta}\) and \(\sigma^2_{\varepsilon}\), yielding a Student-\(t\) distribution, \[\begin{align} p(y_*\mid \boldsymbol{\phi}_*,\mathbf{y},\boldsymbol{\Phi}) = \mathrm{St}\!\Big( y_*;\; \nu = 2a_n,\; m = \boldsymbol{\phi}_*^\top\boldsymbol{\mu}_n, s^2 = \frac{b_n}{a_n} \big(1+\boldsymbol{\phi}_*^\top\boldsymbol{\Lambda}_n^{-1}\boldsymbol{\phi}_*\big) \Big). \label{eq:posterior32prediction} \end{align}\tag{2}\]

In 2 , \(\boldsymbol{\mu}_0\) and \(\boldsymbol{\Lambda}_0\) denote hyperparameters that must be specified. In the Bayesian linear regression and BLL literature, these are commonly set to \(\boldsymbol{\mu}_0 = \mathbf{0}\) and \(\boldsymbol{\Lambda}_0 = \sigma^2_{\beta}\mathbf{I}\), where \(\sigma^2_{\beta}\) is estimated via empirical Bayes, which we also adopt [4]. Specifically, \(\sigma^2_{\beta}\) is obtained by maximizing the closed-form marginal likelihood. See Appendix 6.3 and Appendix 6.4 for more details. Algorithm 4 in Appendix 6.4 outlines the main steps of BDKN.

2.3 Connection to Bayesian model↩︎

Our method has a direct interpretation in terms of Gaussian Processes. Since \(f(\boldsymbol{x})\) is an affine transformation of a Gaussian random vector, it follows that, conditional on \(\sigma^2_{\varepsilon}\), \(f(\boldsymbol{x})\) is a Gaussian Process, \[f(\cdot)\mid\sigma^2_{\varepsilon} \sim \mathcal{GP}\!\big(m(\cdot),k(\cdot,\cdot\mid\sigma^2_{\varepsilon})\big).\] where the mean function is \[m(\boldsymbol{x})=\mathbb{E}[f(\boldsymbol{x})]=\boldsymbol{\phi}(\mathbf{\boldsymbol{x}})^\top\boldsymbol{\mu}_0,\] and the covariance function follows as \[\begin{align} k(\boldsymbol{x},\boldsymbol{x}'\mid\sigma^2_{\varepsilon}) =\mathrm{Cov}\!\big(f(\boldsymbol{x}),f(\boldsymbol{x}')\mid\sigma^2_{\varepsilon}\big) =\sigma^2_{\varepsilon}\,\boldsymbol{\phi}(\mathbf{x})^\top\Lambda_0^{-1}\boldsymbol{\phi}(\mathbf{x'}). \end{align}\]

With \(\Lambda_0 = \sigma_{\beta}^{-2} I\), the kernel reduces to \(k(\boldsymbol{x},\boldsymbol{x}'\mid\sigma^2_{\varepsilon}) = \sigma^2_{\varepsilon}\sigma_{\beta}^2 \sum_{h=1}^H \phi_h(\boldsymbol{x}) \phi_h(\boldsymbol{x}'),\) which corresponds to an isotropic Gaussian prior over the coefficients. See Appendix 6.2 and Appendix 6.3 for more details and derivations. Conditioned on the learned feature maps \(\{\phi_h\}_{h=1}^H\), the hypothesis space of the induced Gaussian process is the finite-dimensional linear span \(\mathcal{H}_H = \mathrm{span}\{\phi_1,\dots,\phi_H\}.\) Marginalizing the noise variance \(\sigma^2_{\varepsilon}\) under the Jeffreys prior yields a Student-\(t\) process [20] with the same finite-dimensional support \(\mathcal{H}_H\), which is consistent with the closed-form Student-\(t\) posterior predictive distribution presented in 2 .

2.4 Computational complexity↩︎

Computationally, BDKN has nearly the same practical cost as a standard deep ensemble. A DE trains \(H\) independent neural networks and aggregates their predictions by uniform averaging, so its dominant cost is the cost of training the \(H\) ensemble members. BDKN uses the same ensemble training stage, but replaces uniform averaging with a Bayesian aggregation layer over the ensemble-member outputs. Given the ensemble-induced design matrix \(\boldsymbol{\Phi}\in\mathbb{R}^{N\times H}\), this aggregation requires forming an \(H\times H\) matrix and solving a small \(H\)-dimensional Bayesian linear regression problem, giving an additional cost of \(O(NH^2 + H^3)\) time and \(O(H^2)\) memory. The \(O(NH^2)\) term comes from forming \(\boldsymbol{\Phi}^\top\boldsymbol{\Phi}\), which sums pairwise interactions between ensemble members over all \(N\) training points, while the \(O(H^3)\) term comes from solving the resulting \(H\times H\) linear system. Since H is small and fixed in our proposed method, this additional cost is negligible relative to neural network training and remains linear in \(N\).

3 Experiments↩︎

We evaluate the performance of the proposed model, namely BDKN, on a set of standard UCI regression benchmarks [21] and compare it against a diverse set of baseline methods, including Gaussian Processes (GP) [10], Deep Ensembles (DE) [3], Deep Kernel Learning (DKL) [11], Bayesian Last-Layer (BLL), Latent Derivative Bayesian Last-Layer (LD-BLL), and Mean-Field Variational Inference Bayesian Neural Networks (MFVI-BNN) [4], and Variational Bayesian Last-Layer (VBLL) [22]. The UCI datasets are widely used benchmarks [3], [4], [22] for probabilistic nonlinear regression, covering a range of regression problems with varying input dimensionalities and sample sizes.

We report the root mean squared error (RMSE) to evaluate predictive accuracy and the negative log-likelihood (NLL) to assess uncertainty calibration. RMSE is defined as \(\mathrm{RMSE} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_i - \hat{y}_i)^2}\), where \(\hat{y}_i = \mathbb{E}[y_i \mid \mathbf{x}_i]\) denotes the predictive mean. NLL is computed as \(\mathrm{NLL} = -\frac{1}{N}\sum_{i=1}^{N}\log p(y_i \mid \mathbf{x}_i)\), where \(p(y_i \mid \mathbf{x}_i)\) denotes the predictive distribution of the model.

Experiments are conducted following [3], using 20 repeated random training/test splits. For each split, 90% of the data are used for training, and the remaining 10% for testing. We report the mean and standard deviation of the evaluation metrics over the 20 splits; for the Protein dataset, which is substantially larger, 5 repeated splits are used, and for the Song dataset, a single split is used. More details about the datasets are provided in Appendix 6.5.1. All methods are evaluated on the benchmark under two experimental settings: an aggressive optimization configuration and a more conservative configuration. In both settings, all neural-network-based approaches use two fully connected hidden layers with 50 neurons each and are trained with a batch size of 32. The only exception is the low-learning-rate setting on Sarcos, where BLL, LDBLL, and MFVI-BNN use three hidden layers with widths \(50\), \(200\), and \(200\), respectively. Our proposed model uses \(H=5\) feature maps (ensemble members). Additional details for each experimental setting and method, including activation functions, are provided in Appendix 6.5.2.

In the first setting, we use a learning rate of \(0.1\) and train all models for 40 epochs, following [3], corresponding to a commonly used lightweight configuration in deep-ensemble studies and is widely adopted in practice. Results for this setting are reported in Table ¿tbl:tab:setting195rmse95nll?. BDKN achieves competitive performance across datasets in both RMSE and NLL, and consistently outperforms DE. While DE exhibits large variance and unstable behavior on several datasets (e.g., Power, Sarcos), BDKN remains stable, demonstrating robustness to variability in optimization quality. GP remains a strong baseline on several datasets. However, due to its cubic computational scaling and prohibitive memory requirements, we could not evaluate GP on large datasets such as Protein, Sarcos, and Song. Compared to other BLL-based methods (BLL, LDBLL) and variational baselines (MFVI, VBLL), BDKN generally provides a better balance between predictive accuracy and uncertainty calibration.

4pt

To also consider scenarios in which more conservative optimization is used, we introduce a second experimental setting with lower learning rates. In this setting, we follow the benchmark configurations proposed by [4], which use dataset-specific learning rates and training schedules. This setting allows us to demonstrate the benefits of our method even when optimization is less aggressive, and to assess robustness to changes in optimization hyperparameters. The exact configurations, including learning rates, numbers of epochs, and training procedures for each dataset, are provided in Appendix 6.5.2.

The empirical performance of all methods under this second experimental setting is summarized in Table ¿tbl:tab:setting295rmse95nll?, where methods exhibit increased stability under this more conservative optimization regime. BDKN remains competitive in both RMSE and NLL and continues to provide well-calibrated uncertainty across datasets. DE performs strongly in this setting, with reduced variations in optimization quality, while GP again achieves strong performance on several benchmarks. In contrast, BLL-based and variational methods exhibit dataset-dependent behavior. The results obtained for LDBLL and VBLL are consistent with those reported in [22]. MFVI also shows strong and stable performance as well, but at a substantially higher computational cost. Runtime comparisons for both experimental settings are reported in Appendix 6.5.3.

4pt

Comparing the two experimental settings, BDKN exhibits robust behavior across both high and low learning-rate regimes. Figure 2 compares test RMSE and NLL over training epochs on the Boston dataset under the two settings. In the high learning-rate setting, all models exhibit stable training dynamics and rapid convergence, except for LDBLL, which degrades after a few epochs (Figure 2 (a)). Moreover, MFVI continues to converge slowly despite the aggressive learning rate. Under the low learning-rate setting, optimization is more stable across methods, with smoother convergence and reduced variance. Figure 2 (b) shows the test-data evaluation over the first 500 training epochs for all methods. GP demonstrates good performance, although with slower convergence. BLL exhibits limited improvement, converging quickly but to higher RMSE and NLL values. LDBLL initially degrades after the first few epochs and then begins to converge very slowly. DKL also shows degradation after a few epochs. In contrast, BDKN and DE maintain fast and stable convergence throughout training. Across both settings, BDKN maintains stable training dynamics and consistent performance in RMSE and NLL, further supporting its robustness.

A key advantage of BDKN is that its feature maps can be trained fully in parallel as ensemble members. Although we use five ensemble members following [3] for simplicity and scalability, the ensemble size can be increased with parallel computing without substantially increasing training time. Figure 3 illustrates the effect of ensemble size on performance for the Yacht dataset under Setting 1 with 50 training epochs. As the number of ensemble members increases, both RMSE and NLL improve consistently, and larger ensembles exhibit faster convergence across splits. This behavior is evident throughout training (Figure 3 (a)) and at the final epoch, where the median RMSE and NLL decrease monotonically with ensemble size (Figure 3 (b)). At the final epoch, Figure 3 (b) shows a clear decrease in the median RMSE and NLL as the ensemble size increases. These results indicate that increasing the ensemble size yields a richer and more expressive feature space, improving predictive accuracy and uncertainty calibration without a significant increase in training time when sufficient computational resources are available.

4 Related work↩︎

Gaussian process. Gaussian Processes (GPs) provide a nonparametric Bayesian framework for regression by placing a prior directly over functions [10]. However, their practical use is often limited by computational complexity and dependence on kernel choice [23], [24], motivating a large body of approximation methods such as sparse and inducing-point approaches [25], [26]. Beyond scalability, numerous extensions have been proposed to improve expressiveness, including neural-network-based kernels and deep GP variants [11], [27]–[29]. Related theoretical work has also established connections between infinitely wide neural networks and GPs, characterizing the function-space priors induced by neural architectures [30]–[32].

Bayesian neural networks. BNNs place prior distributions over network parameters and perform posterior inference to capture predictive uncertainty [12]. While fully Bayesian inference via Monte Carlo sampling methods [33]–[35] provides a principled treatment of uncertainty, such approaches are typically computationally expensive and difficult to scale to modern architectures [36]. As a result, much of the recent work on BNNs has focused on scalable approximate inference techniques, most notably variational inference methods, including mean-field approximations, which trade posterior fidelity for computational tractability [37]–[40].

Bayesian last layer networks. BLLs restrict Bayesian inference to a linear head on top of learned deterministic features, offering a practical trade-off between scalability and uncertainty estimation. Building on this formulation, several extensions refine training and uncertainty behavior. In particular, an efficient reformulation of BLL enables backpropagation-based marginal likelihood training and improves extrapolative uncertainty quantification [41]. To address overconfident predictions and improve uncertainty estimation outside the training distribution, LD-BLL networks introduce functional priors over network derivatives with respect to the inputs to encourage more expressive and calibrated uncertainty [4]. More recently, VBLL has been proposed to further improve scalability and uncertainty estimation through a sampling-free, single-pass training objective [22].

Deep ensembles. DE combines independently trained neural networks, typically via uniform averaging, and is widely used in practice, including reinforcement learning [42]–[47]. Several studies have provided theoretical and analytical perspectives on standard deep ensembles and their uncertainty behavior [48]–[50]. Prior work has shown that promoting ensemble diversity can improve uncertainty estimation in DE [51]–[53]. For instance, ensemble diversity can be encouraged by augmenting the loss function with an auxiliary regularization term that maximizes ensemble predictive variance on inputs sampled uniformly from the input space, improving uncertainty behavior away from the training data [54]. Diversity can also be induced structurally through architectural heterogeneity across ensemble members, such as by randomizing activation functions across networks [55], [56].

5 Discussion↩︎

We introduced BDKN, a scalable framework that bridges the gap between the empirical success of DE and the statistical rigor of GPs. By interpreting independently trained neural networks as induced basis functions, BDKN constructs a finite-rank deep kernel model that enables exact Bayesian inference without the cubic computational costs associated with traditional GPs. Since the ensemble members in BDKN are trained independently, the feature maps can be learned in parallel. This allows the number of feature maps to be scaled to increase model expressivity without a significant increase in training time, which we empirically show leads to improved performance. Our approach offers a clear advantage over standard weight-space BNNs and variational approximations: it yields closed-form posterior predictive distributions together with an interpretable posterior over ensemble members. BDKN also addresses a key limitation of DE, namely their non-diagnostic uniform averaging, by replacing it with a Bayesian aggregation layer that automatically downweights ensemble members weakly supported by the data. This aggregation mechanism improves robustness by reducing the impact of poorly optimized or low-quality ensemble members, which can otherwise degrade performance under aggressive training regimes. While baselines such as DE and LD-BLL show sensitivity to learning-rate selection and can become unstable in high learning-rate regimes, BDKN consistently maintains competitive predictive accuracy and well-calibrated uncertainty across both conservative and aggressive training configurations.

This study focuses on regression; extending the BDKN framework to classification tasks is a natural next step. A key open challenge is maintaining analytic tractability of the posterior predictive distribution under non-Gaussian likelihoods while retaining the robustness and interpretability observed in the regression setting. Further improvements in uncertainty calibration may be achieved by enhancing variability in the learned feature maps or by incorporating principles from latent-derivative modeling. Finally, extending BDKN to other domains, such as deep reinforcement learning and physics-informed machine learning, represents a promising direction for future work.

6 Appendix↩︎

6.1 Jeffreys prior for the noise variance↩︎

We derive the Jeffreys prior for the Gaussian noise variance \(\sigma^2_{\varepsilon}\) in the Bayesian linear model used in Step 2. Consider the observation model \[\mathbf{y}\mid \boldsymbol{\beta},\sigma^2_{\varepsilon} \sim \mathcal{N}\!\big(\boldsymbol{\Phi}\boldsymbol{\beta},\;\sigma^2_{\varepsilon}\mathbf{I}_N\big),\] where \(\boldsymbol{\Phi}\in\mathbb{R}^{N\times H}\) is the ensemble-induced design matrix. The log-likelihood of \(\sigma^2_{\varepsilon}\) is \[\ell(\sigma^2_{\varepsilon}) = -\frac{N}{2}\log\sigma^2_{\varepsilon} -\frac{1}{2\sigma^2_{\varepsilon}} \|\mathbf{y}-\boldsymbol{\Phi}\boldsymbol{\beta}\|^2 + C,\] where \(C\) denotes a constant independent of \(\sigma^2_{\varepsilon}\). Differentiating twice with respect to \(\sigma^2_{\varepsilon}\) yields \[\frac{\partial^2 \ell}{\partial (\sigma^2_{\varepsilon})^2} = \frac{N}{2}\,\frac{1}{(\sigma^2_{\varepsilon})^2} - \frac{\|\mathbf{y}-\boldsymbol{\Phi}\boldsymbol{\beta}\|^2}{(\sigma^2_{\varepsilon})^3}.\] Taking the negative expectation under \(\mathbf{y}\sim\mathcal{N}(\boldsymbol{\Phi}\boldsymbol{\beta},\sigma^2_{\varepsilon}\mathbf{I}_N)\) gives the Fisher information \[\mathcal{I}(\sigma^2_{\varepsilon}) = \frac{N}{2}\,\frac{1}{(\sigma^2_{\varepsilon})^2}.\] By definition, the Jeffreys prior is proportional to the square root of the Fisher information, yielding \[p(\sigma^2_{\varepsilon})\;\propto\;\sqrt{\mathcal{I}(\sigma^2_{\varepsilon})}\;=\;\frac{1}{\sigma^2_{\varepsilon}}.\]

6.2 GP connection↩︎

Recall \[f(x) = \boldsymbol{\phi}(x)^\top \boldsymbol{\beta}, \qquad \boldsymbol{\beta}\mid\sigma^2_{\varepsilon} \sim \mathcal{N}\!\big(\boldsymbol{\mu}_0,\;\sigma^2_{\varepsilon}\boldsymbol{\Lambda}_0^{-1}\big).\]

The conditional mean function is \[\begin{align} m(x) &= \mathbb{E}\!\big[f(x)\mid\sigma^2_{\varepsilon}\big] \\ &= \mathbb{E}\!\big[\boldsymbol{\phi}(x)^\top\boldsymbol{\beta}\mid\sigma^2_{\varepsilon}\big] \\ &= \boldsymbol{\phi}(x)^\top \boldsymbol{\mu}_0. \end{align}\]

The conditional covariance (kernel) is given by \[\begin{align} k(x,x'\mid\sigma^2_{\varepsilon}) &= \mathrm{Cov}\!\big(f(x), f(x') \mid \sigma^2_{\varepsilon}\big) \\ &= \mathbb{E}\!\Big[ \big(f(x)-m(x)\big)\big(f(x')-m(x')\big) \,\Big|\, \sigma^2_{\varepsilon} \Big] \\ &= \mathbb{E}\!\Big[ \big(\boldsymbol{\phi}(x)^\top\boldsymbol{\beta} - \boldsymbol{\phi}(x)^\top\boldsymbol{\mu}_0\big) \big(\boldsymbol{\phi}(x')^\top\boldsymbol{\beta} - \boldsymbol{\phi}(x')^\top\boldsymbol{\mu}_0\big) \,\Big|\, \sigma^2_{\varepsilon} \Big] \\ &= \mathbb{E}\!\Big[ \boldsymbol{\phi}(x)^\top(\boldsymbol{\beta}-\boldsymbol{\mu}_0) (\boldsymbol{\beta}-\boldsymbol{\mu}_0)^\top \boldsymbol{\phi}(x') \,\Big|\, \sigma^2_{\varepsilon} \Big] \\ &= \boldsymbol{\phi}(x)^\top \mathbb{E}\!\Big[ (\boldsymbol{\beta}-\boldsymbol{\mu}_0) (\boldsymbol{\beta}-\boldsymbol{\mu}_0)^\top \,\Big|\, \sigma^2_{\varepsilon} \Big] \boldsymbol{\phi}(x') \\ &= \boldsymbol{\phi}(x)^\top \mathrm{Cov}(\boldsymbol{\beta}\mid\sigma^2_{\varepsilon}) \boldsymbol{\phi}(x') \\ &= \boldsymbol{\phi}(x)^\top \big(\sigma^2_{\varepsilon} \boldsymbol{\Lambda}_0^{-1}\big) \boldsymbol{\phi}(x') \\ &= \sigma^2_{\varepsilon}\,\boldsymbol{\phi}(x)^\top \boldsymbol{\Lambda}_0^{-1} \boldsymbol{\phi}(x'). \end{align}\]

Writing \(\boldsymbol{\phi}(x) = [f_1(x),\dots,f_H(x)]^\top\) and \(\boldsymbol{\Lambda}_0^{-1} = \big[(\boldsymbol{\Lambda}_0^{-1})_{hh'}\big]_{h,h'=1}^H,\) the kernel can be expressed as \[k(x,x'\mid\sigma^2_{\varepsilon}) = \sigma^2_{\varepsilon} \sum_{h=1}^H \sum_{h'=1}^H (\boldsymbol{\Lambda}_0^{-1})_{hh'} \, f_h(x)\, f_{h'}(x').\]

6.2.0.1 Inner-product (finite-rank) form.

Let \(\boldsymbol{\Lambda}_0^{-1/2}\) denote any symmetric square root of \(\boldsymbol{\Lambda}_0^{-1}\) and define the transformed feature map \(\tilde{\boldsymbol{\phi}}(x)\triangleq \boldsymbol{\Lambda}_0^{-1/2}\boldsymbol{\phi}(x)\). With this notation, the kernel has the inner-product representation \[k(x,x'\mid\sigma^2_{\varepsilon}) = \sigma^2_{\varepsilon}\,\tilde{\boldsymbol{\phi}}(x)^\top \tilde{\boldsymbol{\phi}}(x').\] Consequently, for any collection of inputs \(\{x_i\}_{i=1}^N\), the associated kernel matrix can be written as \(K=\sigma^2_{\varepsilon}\tilde{\Phi}\tilde{\Phi}^\top\) with \(\tilde{\Phi}\in\mathbb{R}^{N\times H}\), implying that \(\mathrm{rank}(K)\le H\) independently of \(N\). The induced Gaussian Process is therefore finite-rank, with sample paths that are almost surely confined to the \(H\)-dimensional linear subspace \(\mathrm{span}\{f_1,\dots,f_H\}\). In contrast to nondegenerate (full-rank) kernels such as the squared-exponential or Matern kernels [10], the function-space complexity of the proposed model is explicitly controlled by the ensemble size \(H\).

6.3 Isotropic prior over ensemble coefficients↩︎

In our experiments, we adopt an isotropic Gaussian prior over the ensemble coefficients by setting \[\boldsymbol{\Lambda}_0 = \sigma_{\beta}^{-2}\mathbf{I}_H,\] which corresponds to \(\boldsymbol{\beta}\mid\sigma^2_{\varepsilon} \sim \mathcal{N}(\boldsymbol{\mu}_0,\;\sigma^2_{\varepsilon}\sigma_{\beta}^2 \mathbf{I}_H).\) Under this choice, all ensemble-induced basis functions are a priori treated symmetrically and assigned the same variance, reflecting the absence of additional structural assumptions about the relative importance of individual ensemble members.

Substituting \(\boldsymbol{\Lambda}_0^{-1} = \sigma_{\beta}^2 \mathbf{I}_H\) into the kernel expression derived in the previous section yields \[\begin{align} k(x,x'\mid\sigma^2_{\varepsilon}) &= \sigma^2_{\varepsilon}\,\boldsymbol{\phi}(x)^\top \boldsymbol{\Lambda}_0^{-1} \boldsymbol{\phi}(x') \\ &= \sigma^2_{\varepsilon} \sigma_{\beta}^2 \boldsymbol{\phi}(x)^\top \boldsymbol{\phi}(x') \\ &= \sigma^2_{\varepsilon} \sigma_{\beta}^2 \sum_{h=1}^H f_h(x)\,f_h(x'). \end{align}\] This isotropic prior offers several practical and conceptual advantages. First, it avoids introducing additional biases that would favor specific ensemble members, allowing uncertainty to be driven primarily by the diversity of the learned functions themselves. Second, the single scale parameter \(\sigma_{\beta}^2\) controls the overall strength of the prior in the ensemble-induced function space, enabling efficient empirical Bayes tuning via the marginal likelihood without increasing model complexity. Finally, because the dimensionality of the induced feature space is small (equal to the ensemble size \(H\)), inference remains computationally inexpensive while retaining a principled Bayesian interpretation equivalent to finite-rank Gaussian Process regression.

6.4 BDKN inference↩︎

This section discuss the posterior updates, marginal likelihood optimization, and inference algorithm used in the Bayesian inference procedure of BDKN.

Posterior updates↩︎

Given the design matrix \(\boldsymbol{\Phi}\) and targets \(\mathbf{y}\), conjugacy is maintained and the posterior factorizes as \(p(\boldsymbol{\beta},\sigma^2_{\varepsilon}\mid \mathbf{y},\boldsymbol{\Phi}) = p(\boldsymbol{\beta}\mid\sigma^2_{\varepsilon},\mathbf{y},\boldsymbol{\Phi})\; p(\sigma^2_{\varepsilon}\mid \mathbf{y},\boldsymbol{\Phi}).\) The conditional posterior over coefficients is Gaussian, \(\boldsymbol{\beta}\mid\sigma^2_{\varepsilon},\mathbf{y},\boldsymbol{\Phi} \sim \mathcal{N}\!\big(\boldsymbol{\mu}_n,\; \sigma^2_{\varepsilon}\boldsymbol{\Lambda}_n^{-1}\big),\) with \[\begin{align} \label{eq:mu95n} \boldsymbol{\mu}_n &= \boldsymbol{\Lambda}_n^{-1} \big(\boldsymbol{\Lambda}_0\boldsymbol{\mu}_0 + \boldsymbol{\Phi}^\top \mathbf{y}\big),\\ \boldsymbol{\Lambda}_n &= \boldsymbol{\Lambda}_0 + \boldsymbol{\Phi}^\top\boldsymbol{\Phi}. \end{align}\tag{3}\] Under the Jeffreys prior on \(\sigma^2_{\varepsilon}\), the marginal posterior over \(\sigma^2_{\varepsilon}\) is \(\sigma^2_{\varepsilon}\mid \mathbf{y},\boldsymbol{\Phi}\sim \mathrm{Inv\text{-}Gamma}(a_n,b_n),\) where \[\begin{align} \label{eq:b95n} a_n &= \frac{N}{2},\\ b_n &= \frac{1}{2}\!\left( \mathbf{y}^\top\mathbf{y} + \boldsymbol{\mu}_0^\top\boldsymbol{\Lambda}_0\boldsymbol{\mu}_0 - \boldsymbol{\mu}_n^\top\boldsymbol{\Lambda}_n\boldsymbol{\mu}_n \right). \end{align}\tag{4}\]

Marginal likelihood↩︎

The hyperparameter \(\boldsymbol{\Lambda}_0\), specifically \(\sigma^2_{\beta}\), is optimized by an empirical Bayes estimate, maximizing the marginal likelihood in a closed form. Specifically, one maximizes the log marginal likelihood, written as \[\log p(\mathbf{y}\mid \boldsymbol{\Phi}) = -a_n\log b_n +\frac{1}{2}\log\lvert\boldsymbol{\Lambda}_0\rvert -\frac{1}{2}\log\lvert\boldsymbol{\Lambda}_n\rvert + C, \label{eq:log95marginal95likelihood}\tag{5}\] where \(C\) denotes a constant independent of the hyperparameters. Algorithm 4 outlines the main steps of BDKN.

6.5 Experiments details↩︎

6.5.1 Data↩︎

We evaluated all methods on standard UCI regression benchmarks. Table 1 summarizes the key characteristics of each dataset, including the total number of samples, input dimensionality, number of random splits, and the corresponding train-test sizes. Except for Song, experiments use a 90%/10% train-test split repeated across multiple random splits as reported in the table; Song uses the fixed train-test split shown in the table. In all experiments, input features are standardized using statistics computed on the training data, resulting in zero-mean and unit-variance inputs.

6.5.2 Settings↩︎

Table 2 summarizes the optimization hyperparameters used in the first experimental setting, which corresponds to a high learning-rate configuration. Following the lightweight training regime commonly adopted in DE benchmarks, all models are trained with a learning rate of 0.1 for the full budget of 40 epochs on all datasets. Activation functions follow standard choices for each model class (ReLU for BDKN, DE, and VBLL; Leaky ReLU for BLL, LD-BLL, and MFVI-BNN), GP uses an RBF kernel, and DKL uses a ReLU feature extractor with an RBF GP layer. Results corresponding to this configuration are referred to as Setting 1 throughout the paper.

3pt

Table 3 summarizes the optimization hyperparameters used in the second experimental setting, which corresponds to a lower learning-rate configuration. This setting follows the benchmark protocols proposed by Watson et al. [4], using dataset- and method-specific learning rates and training budgets as reported in the table. For BDKN, DE, DKL, and VBLL, models are trained for the full budget of 500 epochs on all datasets. For BLL, LD-BLL, and MFVI-BNN, the effective number of training epochs is selected via validation: for each split, 80% of the training data are used for optimization and the remaining 20% are used as a validation set to determine the epoch with the best validation performance. The model is then retrained on the full training set using this selected number of epochs. For GP regression, validation-based early stopping is computationally expensive; therefore, we employ a marginal-likelihood–based stopping criterion and terminate optimization when the relative decrease in the average marginal likelihood over the most recent \(\kappa\) iterations, compared to the average over the preceding \(\kappa\) iterations, falls below a threshold \(\rho\), following a procedure similar to that used in prior work [4]. For all datasets, we set \(\kappa = 11\) and \(\rho = 10^{-4}\). Apart from the learning rates and maximum epochs specified in Table 3, all other experimental details are unchanged. Results obtained under this configuration are referred to as Setting 2 throughout the paper.

3pt

For both experimental settings, all neural network-based methods use architectures with two fully connected hidden layers of width 50 and are trained with a batch size of 32 using the Adam optimizer, except for BLL, LD-BLL, and MFVI-BNN on Sarcos in the low-learning-rate setting, which use hidden dimensions \(50,200,200\). DKL uses the same ReLU feature extractor architecture as the other neural-network baselines, with a two-dimensional output representation following Wilson et al. [11], and places an RBF GP layer on top of the learned representation. For MFVI-BNN, we follow the priors and inference procedure described in [4], using independent \(\mathcal{N}(0, \omega/\sqrt{n_{\mathrm{in}}})\) priors for all network weights, where \(n_{\mathrm{in}}\) denotes the number of input features to the corresponding layer and \(\omega = 4\), and independent \(\mathcal{N}(0,1)\) priors for bias terms. Variational inference is implemented in Pyro using diagonal Gaussian variational distributions, optimized via stochastic variational inference (SVI) with the Trace_ELBO objective. For validation and evaluation, predictive distributions are estimated using 100 Monte Carlo samples from the variational posterior. VBLL is implemented using the official Python library (VBLL) provided in [22] with the default settings reported in the paper, and GP regression is implemented using GPyTorch.

For both experimental settings, we initialize the prior variance to \(\sigma_{\beta}^2 = 1\) and perform a small number of marginal-likelihood optimization steps (5 iterations) using the Adam optimizer with a learning rate of \(0.1\), following Algorithm 4. We adopt a zero-mean isotropic Gaussian prior over the linear coefficients, i.e., \(\boldsymbol{\beta} \sim \mathcal{N}(0, \sigma_{\beta}^2 I)\), which is a standard and commonly used choice in Bayesian last-layer models [4].

6.5.3 Runtime↩︎

All experiments were conducted on a workstation running Windows 10 Enterprise (64-bit), equipped with a 12th-generation Intel Core i7-12700 CPU (12 cores, 20 logical processors), 32 GB of RAM, and an NVIDIA RTX A2000 GPU with 12 GB of dedicated memory. The system uses NVMe SSD storage. GPU computations were performed using CUDA 12.8 with NVIDIA driver version 573.44. Tables 4 and 5 report the training time (in minutes) for all methods under the two experimental settings. The runtimes correspond to end-to-end training on standard regression benchmarks using the same hardware and implementation framework. Under Setting 1 (high learning rate), BDKN exhibits runtime comparable to standard DE, reflecting that its additional Bayesian aggregation step incurs negligible overhead relative to training the ensemble members. In contrast, fully Bayesian or variational approaches, such as MFVI-BNN and VBLL, incur substantially higher computational costs on larger datasets due to sampling-based inference or variational optimization. GP is efficient on smaller datasets but becomes expensive or infeasible on larger datasets, such as Protein, Sarcos, and Song. Under Setting 2 (low learning rate), overall runtimes increase across all neural network-based methods due to longer training schedules, with MFVI-BNN showing particularly large increases on high-dimensional or large-scale datasets. Despite this, BDKN continues to scale comparably to DE and remains significantly more efficient than fully Bayesian baselines. Entries marked as N/A indicate cases where a method could not be evaluated due to GPU memory constraints.

6.6 Impact Statement↩︎

This paper presents work whose goal is to advance the field of Machine Learning by developing a scalable and analytically tractable Bayesian framework for uncertainty quantification in regression. The proposed Bayesian Deep Kernel Networks (BDKN) are intended to improve the reliability and interpretability of predictive uncertainty relative to commonly used baselines (e.g., Gaussian Processes, Bayesian last-layer models, and deep ensembles), which can be beneficial for downstream decision-making in risk-sensitive settings.

Potential positive societal impacts include enabling more trustworthy model deployment by reducing overconfident predictions, supporting better uncertainty-aware decisions in applications such as scientific modeling, engineering design, and other domains where regression is a core component. By providing a transparent posterior over ensemble-member contributions and closed-form predictive uncertainty, this work may also facilitate diagnostic analyses that help practitioners detect brittle predictors and improve model accountability.

To the best of our knowledge, this work does not introduce direct negative societal impacts or high-risk misuse concerns.

References↩︎

---

1. Corresponding author: j.lee@tamu.edu↩︎