Bayesian Inference for Consistent Predictions
in Overparameterized Nonlinear Regression

Bayesian GLM has been studied intensively for high dimensions. The first paper to study high-dimensional Bayesian GLM [19] realized posterior concentration for the parameter without assuming any sparsity structure by imposing a restricted growth condition on the parameter dimension. This underparameterized setting differs from the present work. Subsequently, several studies explored the GLM setting where the number of parameters may exceed the sample size [25, 35, 24]. These works have demonstrated the asymptotic properties of the posterior distribution; however, they assume sparsity, implying that the parameters are inherently low-dimensional. The properties of Bayesian variable selection in GLMs have also been investigated in similar domains [9, 10, 43, 31, 23]. However, the assumption of parameter sparsity can be stringent, and the current study is conducted in a setting that does not rely on this assumption, thus expanding the scope of research in this area. The methods and theory of Bayesian inference in high-dimensional models are summarized in [2].

Numerous theoretical results suggest that overparameterization leads to good generalization performance. [6] demonstrated that the excess risk converges to zero in overparameterized linear models. This theory has been widely extended to other linear estimators, such as regularized regression [53, 30, 59, 37, 12], kernel regression [32, 33], linear neural networks [11, 51], and regression for dependent data [41, 54]. [14, 7, 27] ventured beyond these linear frameworks to investigate the benign nature of overparameterization, but their focus was limited to classification problems. [48] tackled nonlinear regression and demonstrated the absence of benign overfitting, but their assumptions differ from the present work mainly in terms of distribution of the covariate.

The remainder of this paper is structured as follows. In Section 2, the problem setting for one-parameter generalized linear models under the overparameterization assumption is introduced, providing a detailed description of the proposed methodology. Section 3 presents the theoretical results on posterior contraction for GLMs. Section 4 extends the results to the two-parameter case. Section 5 reports on the simulation studies conducted to validate the theoretical findings and assess the uncertainty estimation. Section 6 showcases a real data application, comparing the proposed approach with baseline methods in terms of classification performance. Section 7 concludes the paper, discussing the contribution and limitations of this work. The proof and discussion of the main theorems and details of the computational algorithms used in the analysis are provided in the Appendix.

For $q\in\mathbb{N}$ and any vector $\bm{v}$, $\|\bm{v}\|_{q}$ denotes the $\ell_{q}$-norm of $\bm{v}$. For any vector $\bm{v}$ and symmetric matrix of the same dimension $A$, we define $\|\bm{v}\|_{A}=\sqrt{\bm{v}^{\top}A\bm{v}}$. For sequences $\{a_{n}\}_{n\in\mathbb{N}}$ and $\{b_{n}\}_{n\in\mathbb{N}}$, $a_{n}\gtrsim b_{n}$ denotes the existence of a constant $c>0$ such that $a_{n}\geq cb_{n}$ for all $n\geq\overline{n}$ with some finite $\overline{n}\in\mathbb{N}$, and $a_{n}\lesssim b_{n}$ denotes the opposite; $a_{n}=o(b_{n})$ indicates that $a_{n}/b_{n}\to 0$ as $n\to\infty$, whereas $a_{n}=\omega(b_{n})$ indicates that $a_{n}/b_{n}\to\infty$ as $n\to\infty$; $a_{n}\gtrapprox b_{n}$ indicates that $a_{n}\gtrsim b_{n}(\log n)^{c}$ holds for every $c>0$, and $a_{n}\lessapprox b_{n}$ denotes the opposite. For natural numbers $a,b\in\mathbb{N}$ with $a<b$, we define $[a:b]:=\{a,a+1,a+2,...,b\}$. For any square matrix $A$, let $A^{-1}$ be the inverse matrix of $A$, $A^{\dagger}$ be the Moore–Penrose pseudoinverse matrix, and $\mathrm{tr}(A)$ be the sum of the diagonals of $A$. For a pseudo-metric space $(S,d)$ and a positive value $\delta$, $\mathcal{N}(\delta,S,d)$ denotes the covering number, that is, the minimum number of balls that cover $S$ with a radius $\delta$ in terms of $d$. For two probability measures $\Pi_{1}$ and $\Pi_{2}$ on the same measurable space $(\mathcal{X},\mathcal{B})$, the total variation distance is given by $\|\Pi_{1}-\Pi_{2}\|_{\mathrm{TV}}:=\sup_{B\in\mathcal{B}}|\Pi_{1}(B)-\Pi_{2}(B)|$, and the Kullback–Leibler (KL) divergence and the KL variation (also referred to as the centered second-order KL divergence) between them are respectively defined as $K(\Pi_{1},\Pi_{2})=\int_{\mathcal{X}}\log(\mathrm{d}\Pi_{1}/\mathrm{d}\Pi_{2})\mathrm{d}\Pi_{1}$ and $V(\Pi_{1},\Pi_{2})=\int_{\mathcal{X}}\{\log(\mathrm{d}\Pi_{1}/\mathrm{d}\Pi_{2})-K(\Pi_{1},\Pi_{2})\}^{2}\mathrm{d}\Pi_{1}$ if $\Pi_{1}$ is absolutely continuous with respect to $\Pi_{2}$; otherwise, both are $\infty.$ The Hellinger distance between two probability distributions $\Pi_{1}$ and $\Pi_{2}$ is defined as $H(\Pi_{1},\Pi_{2})=\sqrt{1-\int_{\mathcal{X}}\sqrt{\mathrm{d}\Pi_{1}\mathrm{d}\Pi_{2}}}$.

Consider the pair $(\bm{X},Y)$, where $\bm{X}$ is an $\mathbb{R}^{p}$-valued centered random covariate, and $Y$ is an $\mathbb{R}$-valued univariate response. The GLM is defined as follows:

where $\bm{\beta}^{*}\in\mathbb{R}^{p}$ is the true unknown $p$-dimensional regression coefficient, and $g:\mathbb{R}\to\mathbb{R}$ is a known Lipschitz continuous function. In this setting, the response variable $Y$ follows a general distribution belonging to the one-parameter exponential family with a density function given by

where $\mu_{x,\bm{\beta}}\in\mathbb{R}$ is the parameter; $a:\mathbb{R}\to\mathbb{R}$ and $b:\mathbb{R}\to\mathbb{R}$ are known continuously differentiable functions, with $a$ having a nonzero derivative. The function $c:\mathbb{R}\to\mathbb{R}$ is chosen such that $L(y;\mu_{x,\bm{\beta}})$ is a proper probability density.

This formulation encompasses several well-known models: Normal linear models: $Y\mid\bm{X}$ follows a $g(\bm{x}^{\top}\bm{\beta}^{*})$-mean Gaussian distribution with fixed variance, and $g$ is the identity function; Logistic regression model: $Y\mid\bm{X}$ follows a Bernoulli distribution with parameter $g(\bm{x}^{\top}\bm{\beta}^{*})$, and $g$ is the logistic function; Poisson regression model: $Y\mid\bm{X}$ follows a Poisson distribution with parameter $g(\bm{x}^{\top}\bm{\beta}^{*})$, and $g$ is the exponential function. Note that the link function $g$ is not restricted to being canonical, thereby providing flexibility in modeling the relationship between the linear predictor and response variable.

This study focuses on a non-sparse overparameterization setting, where $n$ independent and identically distributed copies of $(\bm{X},Y)$ are generated, and their realizations $\mathcal{D}:=\{(\bm{x}_{1},y_{1}),\ldots,(\bm{x}_{n},y_{n})\}$ are observed. From these observations, a $p$-dimensional estimate $\widehat{\bm{\beta}}$ is constructed for prediction, with $p$ being larger than the sample size $n$. In the asymptotic setting, where $n$ diverges, $p$ possesses a larger order of divergence, i.e., $p=\omega(n)$.

The primary objective of this study is to examine the generalization performance of a predictor trained on the observations $\mathcal{D}$ in the overparameterized domain. We adopt as the predictor the mean $g(\bm{x}^{\top}\bm{\beta})$, with the estimated parameter $\widehat{\bm{\beta}}$ plugged in. Then, the notion of excess risk is introduced to quantify the generalization performance, measuring the difference between the oracle predictor $g(\bm{X}^{\top}\bm{\beta}^{*})$ and the predictor $g(\bm{X}^{\top}\widehat{\bm{\beta}})$:

where $\|\cdot\|_{P_{\bm{X}},2}$ denotes the $L_{2}(P_{\bm{X}})$ norm with respect to the distribution $P_{\bm{X}}$ of the covariates. As the later analysis shows, the convergence of the excess risk is related to the distributional convergence of the estimated parameter $\widehat{\bm{\beta}}$.

Predictor construction follows a systematic approach: first, the data are divided into two sets; then, a prior distribution of $\bm{\beta}$ is developed based on the empirical spectral information of $\Sigma$ from one set; next, the likelihood is computed in the other set; finally, the posterior distribution of $\bm{\beta}$ is computed based on the prior and likelihood.

Initially, we divide the dataset $\mathcal{D}$ into two equal parts, $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, each containing $n/2$ observations. Although the split ratio is arbitrary, it should not depend on the sample size $n$. The data splitting allows the use of $\mathcal{D}_{1}$ for constructing the prior distribution and $\mathcal{D}_{2}$ for evaluating the likelihood function, ensuring that the data used in these two steps are independent.

Next, we consider the empirical covariance matrix of the covariates $\bm{X}$, defined as:

where $\bar{\bm{x}}\in\mathbb{R}^{p}$ is the sample mean of the covariates from $\mathcal{D}_{1}$. The factor $2/n$ ensures that $\widehat{\Sigma}$ is an unbiased estimator of the population covariance matrix, whereby a spectral decomposition of $\widehat{\Sigma}$ is performed:

where $\{\widehat{\lambda}_{j}\}_{j=1}^{p}$ and $\{\widehat{\bm{v}}_{j}\}_{j=1}^{p}$ are the eigenvalues and eigenvectors of $\widehat{\Sigma}$, respectively, with $\widehat{\lambda}_{1}\geq\widehat{\lambda}_{2}\geq\cdots\widehat{\lambda}_{p}\geq 0$. Using the empirical spectra, we define a low-rank approximation of $\widehat{\Sigma}$ as:

This approximation captures the most significant $k$ eigenvalues and eigenvectors of $\widehat{\Sigma}$, effectively reducing the dimensionality of the covariate space while retaining the most informative directions of variation. In choosing $k$ as a tuning parameter to manage the computational cost, it has to be fixed to an interval of $[L_{\kappa}:U_{\kappa}]$ so that it satisfies Assumption 2.

We then employ the following prior distribution for $\bm{\beta}$.

Although this prior distribution resembles a $p$-variate Gaussian distribution with mean $\bm{0}$ and covariance $\widehat{\Sigma}_{1:k}$, there are some important distinctions. First, the distribution is truncated, with the indicator function restricting the domain to a centered sphere of radius $R$. This truncation is necessary because the low-rank approximation $\widehat{\Sigma}_{1:k}$ is singular, and defining the density function without support restrictions would cause the integrals over $\mathbb{R}^{p}$ to diverge. Second, unlike [43, 24, 55], we impose a non-sparse prior distribution, even though some parameters may be redundant. This choice is motivated by the desire to retain all potentially relevant covariates in the data and avoid the challenges associated with sparse priors, such as the need to specify sparsity-inducing hyperparameters.

Next, we consider the posterior distribution, which combines the prior distribution and the likelihood function using Bayes’ Theorem.

The Bayesian framework allows the quantification of events of interest, such as the probability of predictions being close to the true values through probability measures defined in this manner. Moreover, as demonstrated later, the distributional information provided by the posterior is valuable in practice.

Note that the information about the response variable $y_{i}$ in $\mathcal{D}_{1}$ is discarded since this subset is only used to compute the sample covariance matrix of the covariates. If excluding this information seems inefficient, the roles of $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ can be switched to construct a new estimator (posterior distribution), and the two estimators can be averaged in a manner similar to the sample splitting scheme [8]. The theoretical results presented next preserve this averaging.

The initial focus is on the boundedness of the operator norm of $\Sigma$ and concentration of the empirical covariance matrix $\widehat{\Sigma}$ around $\Sigma$.

This assumption ensures that the empirical covariance matrix concentrates around the population covariance matrix. The sequence $\rho_{n}$ quantifies the rate of this concentration, whereas $\ell_{n}$ controls the probability of the concentration event. Notably, this assumption is mild and includes not only the sub-Gaussian distributions commonly imposed in high-dimensional statistics but also log-concave distributions, which encompass most unimodal parametric distributions [57, 58] that allow for heavy-tailed distributions such as the Laplace distribution. In contrast, other related literature impose more restrictive conditions such as sub-Gaussianity [6] or strong log-concavity [14]. In a high-dimensional setting, where many covariates are incorporated, the mildness of the distributional assumption is particularly important, as it allows for a wider range of applicable scenarios.

The next assumption concerns the behavior of the eigenvalues of $\Sigma$.

This assumption imposes certain regularity conditions on the eigenvalues of $\Sigma$. Intuitively, conditions (i) and (ii) control the rate at which the bulk of the eigenvalues (those indexed from $L_{\kappa}$ to $U_{\kappa}$) are allowed to decay moderately. Condition (iii) requires that the eigenvalues beyond the bulk (those indexed larger than $U_{\kappa}$) are sufficiently small, whereas condition (iv) constrains the size of $U_{\kappa}$ relative to sample size $n$ and rate $\varepsilon_{n}$. These conditions are crucial for deriving the posterior contraction rates in Theorem 1.

Finally, we assume a certain local Lipschitz condition of the link function $g$ in terms of the Kullback-Leibler (KL) divergence and variational form between the corresponding probability distributions.

This assumption limits the amount of variation in a distribution when different parameters are used, to control the complexity of the model and is satisfied by several common GLMs, as illustrated in the following examples.

These examples demonstrate that Assumption 3 is a reasonable condition that holds for several widely used GLMs.

We now present the main result, which establishes the posterior contraction for the GLM under the assumptions introduced in the previous subsection.

This theorem shows that under conditions suitable for the covariance matrix $\Sigma$ and prior distribution, the posterior distribution of the predictor concentrates around the oracle predictor at a rate $\varepsilon_{n}$ with respect to the $L_{2}(P_{\bm{X}})$ distance. The rate $\varepsilon_{n}$ depends on the interplay between sample size $n$, concentration rate $\rho_{n}$ of the empirical covariance matrix, and eigenvalues of $\Sigma$ in the bulk region.

Here, we consider a binary classification scenario where the response variable $Y_{i}$ takes values in $\{0,1\}$. The data are generated according to the following model:

where $g_{ls}(r):=1/(1+e^{-r})$ is the logistic function. Specifically, for $n=50,100,\ldots,500$, to ensure an overparameterized setting, $\bm{\beta}\sim N(\bm{0},I_{p})$ are generated with $p=n^{4/3}$. Two settings are considered for the covariate: (A) Gaussian covariate, $\bm{X}_{i}\sim N(\bm{0},\Sigma)$; (B) Laplacian covariate, $\bm{X}_{i}\sim L(\bm{0},\Sigma)$, where the eigenvalues of $\Sigma$ follow $\lambda_{j}(\Sigma)=10\exp(-j/8)+n\exp\{-\sqrt{n}\}/p$. Then, each $Y_{i}$ corresponding to $X_{i}$ and $\beta$ is sampled from (20).

For every $n$, $20$ different datasets were repeatedly generated through the above procedure and each was analyzed to determine the behavior of the proposed method in response to an increase in sample size. In each case, a test dataset of $1000$ instances was generated, to evaluate the performance of the trained classifier. The performance of the Bayesian predictions was evaluated, employing four widely used metrics: 0-1 loss, which quantifies the misclassification rate; area under the ROC curve (AUC), which assesses the discriminative ability of the classifier; rate of unconfident misclassifications (UM), which measures the proportion of unconfident predictions among misclassified instances; and confident correct (CC) predictions rate, which calculates the proportion of accurate predictions among confident ones. We refer to a confident prediction as one whose $95\%$ predictive interval does not cross the threshold $(0.5)$, whereas an unconfident one does. The latter two metrics are helpful when building systems that defer action when there is no confidence in the prediction. Additionally, principal component regression (PCR) [5] was performed as a baseline method, with the number of principal components selected by 5-fold cross-validation from $\{1,2,\ldots,30\}$.

Figure 1 illustrates the results, where the reported values are averaged over the $20$ datasets. As sample size increases, the 0-1 loss exhibits a consistent decrease in both settings, indicating an improvement in classification accuracy. This trend persists even in the overparameterized setting, aligning with the theoretical findings. The AUC metric also shows a steady increase with sample size, suggesting an enhanced ability to discriminate between the two classes. Our method outperformed PCR, and these results underscore the importance of overparameterization for prediction. The UM and CC metrics provide insights into the uncertainty estimation of the proposed method. UM consistently presents a high value. This indicates that the method can warn us of a lack of confidence when misclassifying. The CC also presents a similarly high value, suggesting that confident predictions are trustworthy. These observations suggest that the proposed method effectively captures the uncertainty associated with predictions, thereby leveraging it to make more informed decisions.

In this subsection, the scenario of true responses taking values along the positive real line $\mathbb{R}_{+}$ is considered. Specifically, the data are generated as follows:

where $g_{sp}(r):=\log(1+e^{r})$ is the softplus function. For $n=50,100,\ldots,500$, we set $p=n^{4/3}$ and sample the coefficient $\bm{\beta}\in\mathbb{R}^{p}$ and covariate vectors $X_{i}\in\mathbb{R}^{p}$ in the same manner as in the previous experiment. Then, the scalar response is generated by (21).

For every $n$, $20$ distinct datasets were generated using the above procedure and, on each one, Bayesian prediction performance changes were examined with an increase in training data size. For a newly generated test dataset $\{(\bm{x}_{i},g_{sp}(\bm{x}_{i}^{\top}\bm{\beta}))\}_{i=1}^{1000}$, the accuracy of Bayesian predictions was assessed using three widely adopted indicators: root-mean-squared error (RMSE), measuring the average deviation between predicted and true values; coverage probability (CP), evaluating the reliability of the $95\%$ prediction intervals by calculating the proportion of true values falling within the confidence intervals; and average length (AL), quantifying the width of the $95\%$ prediction intervals, providing insight into the precision of the predictions. Moreover, PCR was implemented, selecting the number of principal components using 5-fold cross-validation from $\{1,2,\ldots,30\}$.

The results, averaged over $20$ datasets, are shown in Figure 2. First, the RMSE decreases monotonically as sample size increases. Even when $p$ is large and the setting is clearly overparameterized, the fact that accuracy continues to improve supports the theoretical results experimentally. Also, our method excels PCR, providing advantages of our method for prediction problems. Next, we focus on CP and AL, which are indicators related to prediction intervals, with CP representing the reliability of the prediction interval and AL representing its width. In the figure, the coverage is not sufficient when $n=100$; however, after that, CP stabilizes at a high value. AL also slightly increases but not by much, implying that the interval is not trivial. (Regardless of $n$, the mean of $\{g_{sp}(\bm{x}_{i}^{\top}\bm{\beta})\}_{i=1}^{1000}$ in the test data is around $3.5$, and the variance is around $5$.) These results suggest that the proposed method properly quantifies the uncertainty in constructing prediction intervals.

We apply the proposed method to the ARCENE dataset, which was originally curated for the NIPS 2003 feature selection challenge for the binary classification problem. The dataset consists of mass spectrometric measurements from cancer patients and healthy individuals, with the goal of distinguishing between the two groups based on these measurements. The ARCENE dataset, accessible through the UC Irvine Machine Learning Repository [17], was constructed by merging three separate datasets: National Cancer Institute (NCI) ovarian cancer data, NCI prostate cancer data, and Eastern Virginia Medical School (EVMS) prostate cancer data. In total, the public dataset contains $200$ samples, comprising $88$ cancer samples and $112$ normal samples. These data are categorized as cancer and control samples, with $10000$ continuous covariates, with 7000 dimensions corresponding to real mass-spectrometry peaks and $3000$ dimensions representing random noise variables (probes). The primary challenges posed by the dataset lie in its high dimensionality and the presence of overlapping classes. The objective here is to discriminate between cancerous and normal samples.

The performance of the six distinct methods were compared.

• •

  Logistic Regression with effective spectral (ES) prior distribution: The proposed Bayesian methodology, as described in Section 2.2, was applied to the logistic regression model to make predictions on unseen data.

• •

  Logistic Regression with the horseshoe (HS) prior: To handle the high-dimensionality of covariates, logistic regression was implemented using the Horseshoe prior [13]. The ”horseshoenlm” package in R [39] was used for the implementation.

• •

  Logistic Regression with the normal (NM) prior: We investigated the performance of logistic regression with a standard normal prior distribution, which is the most commonly used prior and serves as a ridge-like regularization [16]. This prior provides flexibility in modeling the coefficient distribution and can be easily implemented using the Pyro probabilistic programming language [3].

• •

  Logistic Regression with the Laplace (LP) prior: the performance of logistic regression was also explored using a standard Laplace prior distribution, which promotes sparsity in the coefficient estimates as a Bayesian Lasso [44]. This prior has proven effective in high-dimensional settings and can be efficiently implemented using Pyro, leveraging its support for non-Gaussian prior distributions [3].

• •

  LightGBM: As a baseline method for table data prediction tasks, LightGBM [28], which is an efficient and scalable gradient boosting framework, was implemented adopting the ”lightgbm” package in Python for the implementation and training of the model using Bayesian optimization [49, 15] via the ”Optuna” package [1].

• •

  Principal Component Regression: As a representative method for high-dimensional data, PCR was implemented with the number of bases selected from $1$ to $30$ by 5-fold cross-validation.

For all methods, performance evaluation was conducted using 5-fold cross-validation. In each fold, the model was trained on the training data, and predictions were made on the test data. As in the previous section, the evaluation metrics used to assess predictive performance were Accuracy, AUC, UM, and CC.

The results of the experiment are presented in Table 1. The entries represent the average values of the metrics over $5$-fold cross-validation. Overall, among all the methods, logistic regression with the ES prior distribution achieves the highest accuracy and AUC, demonstrating its effectiveness in capturing the underlying structure of the high-dimensional ARCENE dataset. The superior performance of the ES prior can be attributed to its ability to adapt to the intrinsic dimensionality of the data by incorporating the empirical spectral information of the covariance matrix [60].

Let us take a closer look at the comparison. The horseshoe prior, known for its flexibility and sparsity-inducing properties [13], performs well, with the second-highest accuracy and AUC. This suggests that the regression on ARCENE dataset may exhibit some degree of sparsity, which the HS prior is able to capture effectively. However, the ES prior outperforms the HS prior, indicating that the adaptive nature of the ES prior is more beneficial in this context. Also, the fact that the $3000$ covariates are essentially irrelevant to the regression and have zero coefficients but are not highly sparse, may be the reason why the proposed method prevails. The LP and NM priors, while still providing competitive results, do not perform as well as the ES and HS priors. Nevertheless, their performance is still superior to the LightGBM baseline, highlighting the advantages of the Bayesian approach in handling excess covariates. PCR effectively selects relevant components from a large number of covariates although the results are only moderate and do not match those of our method which retains all covariates without reduction.

A notable aspect of the results is the uncertainty estimation of the Bayesian methods. The ES prior achieves the highest UM and CC values, indicating that it is highly effective in identifying uncertainty and providing reliable confidence estimates. This is particularly valuable in practical applications, where the ability to defer decision-making when the model is uncertain can prevent the risk of misclassifications. The HS prior also demonstrates good uncertainty quantification, with the second-highest UM value. However, the LP and NM priors exhibit very low UM values, suggesting that they may be overconfident in their predictions. Note that the LightGBM method does not provide uncertainty estimation, as it is not a probabilistic model.