Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
128,835 characters · 11 sections · 75 citation commands
Debiased Bayesian Inference for High-dimensional Regression Models
\pdfbookmark[1]{Title}{title}
Applied researchers now routinely work with regression models that feature a large number of covariates. A primary inferential goal in econometrics is to estimate the ceteris paribus effect of a specific variable while controlling for other variables BelloniChernoHansen2013HD, BelloniChernozhukovCherverikovWei2018Many. The prevailing practice interprets the coefficient on a regressor as a causal effect, conditional on the included controls. As the plausibility of conditional unconfoundedness is often argued using a large set of covariates, practitioners have increasingly embraced high-dimensional regression models. This setting has been extensively studied, predominantly using frequentist methods.
Bayesian inference, on the other hand, has long been valued for its coherent framework for handling uncertainty in statistical analysis. As highlighted by Rubin1984Applied, Bayesian methods provide direct answers to many empirical questions by quantifying uncertainty about unknown parameters conditional on the observed data.\footnote{In discussing the commonly-used rule-of-thumb confidence interval, Rubin1984Applied states that “the interval is-at least in my experience-nearly always interpreted Bayesianly, that is, as providing a fixed observed interval in which the unknown (parameter of interest) $\mu$ lies with 95% probability.” As advocated by Imbens2021Pvalue, it may even be preferable to adopt Bayesian inference in cases where Bayesian and frequentist procedures lead to different conclusions.} This appeal has grown among applied researchers, who often seek probabilistic statements about particular parameters of interest given the specific dataset at hand. Bayesian posterior distributions conveniently encapsulate both sampling variation and parametric uncertainty, unifying estimation and inference in a way that aligns well with the needs of empirical work.
In addressing the challenges posed by high dimensionality, recent methodological advances have substantially expanded the Bayesian toolkit for regression models with many covariates. Notable progress includes the development of spike-and-slab priors MitchellBeauchamp_BayesianVariable_1988,GeorgeMcCulloch_VariableSelection_1993 and horseshoe priors CarvalhoPolsonScott2010Horseshoe, as well as scalable approximate Bayesian inference techniques RaySzabo2022Variational. These innovations have attracted a growing interest in economics, as illustrated by GiannoneLenzaPrimiceri2021Sparsity. However, while much of this literature focuses on model selection and estimation, relatively less attention has been devoted to delivering valid inference for a particular parameter of interest in the presence of high-dimensional controls---a central need in many empirical applications. This paper aims to fill this gap.
Specifically, we develop a novel inferential procedure for high-dimensional linear regression that introduces a debiasing step for the Bayesian posterior distribution of the parameter of interest. Building on the concept of debiased point estimators, which is well established in the frequentist literature, we instead debias the entire posterior distribution, yielding a new debiased Bayesian approach. Our method is tailored to high-dimensional settings and constructs credible sets from the debiased posterior, while remaining faithful to Bayesian principles by conditioning on the observed data.
The core theoretical contribution of our work is the establishment of new Bernstein-von Mises (BvM) results for the debiased Bayesian procedure. Our framework allows the number of covariates $p$ to exceed the sample size $n$, with $p$ growing exponentially in $n$. These results formally justify the asymptotic normality of the debiased posterior and ensure that credible sets achieve correct frequentist coverage under repeated sampling, thereby helping to bridge the gap between Bayesian and frequentist inference. Importantly, our results show that, after debiasing, Bayesian credible sets can match the performance of debiased frequentist methods---such as double machine learning ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018Double---in high-dimensional regimes. This effort also echoes an important point made by Rubin_BayesianBootstrap_1981 regarding the assessment of Bayesian procedures under repeated sampling. Our framework is notably general. It is agnostic to the choice of prior for the regression coefficients and permits flexible selection of pilot estimators for the precision matrix. The specific implementations we present in Section (ref) illustrate how the method can be practically adapted to high-dimensional settings. This flexibility allows the procedure to be tailored to a wide range of empirical applications without sacrificing computational scalability.
Another key advantage of our approach is its compatibility with computationally efficient approximate Bayesian inference, such as variational Bayes RaySzabo2022Variational. Traditional spike-and-slab priors, while theoretically appealing, often entail substantial computational cost due to the reliance on Markov chain Monte Carlo (MCMC) sampling. In contrast, variational Bayes delivers approximate posteriors within seconds, while still attaining optimal contraction rates in high-dimensional settings. However, without debiasing, credible sets constructed from such approximate posteriors are not guaranteed to achieve correct frequentist coverage. By incorporating our debiasing procedure, the variational Bayes posterior becomes a principled tool for valid uncertainty quantification.
Our approach parallels an important line of work in the frequentist literature, in particular the debiasing of LASSO-type estimators JavanmardMontanari2014CI,VandeGeer_OnAsymptotically_2014,ZhangZhang2014CI. The LASSO estimator itself can be viewed as the posterior mode under a Laplace prior, yet it typically exhibits bias and therefore requires post-estimation correction. Whereas debiased frequentist methods adjust a point estimator and then rely on its asymptotic properties for inference, our Bayesian debiasing method extends this idea to the entire posterior distribution. Because the debiasing step is applied at the posterior level, we develop new technical arguments to establish its large-sample properties. Monte Carlo simulations show that our procedure delivers substantial improvements over standard (uncorrected) Bayesian methods and remains competitive with debiased frequentist procedures in terms of coverage and estimation accuracy. In particular, when regression coefficients are large, our method offers pronounced advantages over debiased frequentist approaches, providing practitioners with a powerful tool for empirical analysis in high-dimensional settings.
Our work contributes to the active literature on the frequentist validity of Bayesian inference in high-dimensional regression. Early contributions by Ghosal1999HDL and Bontemps2011BvM established asymptotic normality of the posterior in Gaussian linear regression when the number of covariates grows slowly with the sample size. However, these results do not accommodate sparsity and require the ambient dimension to remain smaller than the sample size (see Bontemps2011BvM). More recently, Castilloetal_BayesianSparse_2015 analyzed Bayesian posteriors under spike-and-slab priors, establishing mixed normality under strong signal conditions; see also WuNarisettyYang2023ConditionalBayes for extensions. Yang2019HDL introduced a debiasing approach based on reparameterization and targeted prior modification, assigning a data-dependent Gaussian prior to the parameter of interest. Building on this idea, Castilloetal2024VB incorporated mean-field approximation strategies to distinguish between high-dimensional nuisance parameters and parameters of interest. In contrast, our approach attains debiasing via a post-processing step applied to the posterior, preserving standard Bayesian modeling and computation and allowing practitioners to directly leverage existing Bayesian methods for high-dimensional regression.
Our paper also connects with broader developments in semiparametric Bayesian inference and correction methods for posterior or prior distributions RayVaart2020Semi,BreunigLiuYu2022DR,BreunigLiuYu2024DiD,BreunigLiuYu2025Bart,YiuFongHolmesRousseau2023Semiparametric. These works depart from naive plug-in principles to achieve robust large-sample properties under weak conditions. We contribute to this literature by explicitly addressing high-dimensional regression with sparsity-inducing priors and by showing that debiasing can be performed on approximate posteriors. Theoretically, we also establish BvM results for a growing subset of coefficients, which is a novel contribution in this context. Recently, DiTragliaLiu2025BML proposed a Bayesian analog of double machine learning for linear regression when the number of covariates is relatively small, in the sense that $p = o(n)$. Their approach relies on priors for reduced-form covariance matrices to remain consistent with the likelihood principle Walker_Parametrization_2023. However, when $p$ is large relative to $n$, practical computation of posteriors for high-dimensional covariance matrices remains challenging.
The remainder of the paper is organized as follows. Section (ref) introduces our inferential procedure and discusses implementation details, while Section (ref) presents two empirical applications. In Section (ref), we develop our main theoretical results under high-level conditions, which are then verified under more primitive assumptions. Section (ref) examines the finite-sample performance of our procedure through Monte Carlo simulations. Section (ref) concludes. All proofs and supporting results are relegated to the appendices, which also contain additional simulation evidence.
Consider the prototypical linear regression model:
where $\{Y_i,X_i\}_{i=1}^n$ are independent and identically distributed (i.i.d.) and $\mathbb{E}_0[\cdot]$ is the population expectation. The covariate $X_i$ is a $p$-dimensional vector, with $p$ potentially exceeding the sample size $n$. The true regression coefficient vector $\beta_0$ is assumed to be sparse, in the sense that only a small subset of regressors have nonzero effects. We denote the population covariance matrix of $X_i$ by $\Omega_0:=\mathbb{E}_0[X_iX_i^\intercal]$, and its inverse, the precision matrix, by $\Theta_0:=\Omega_0^{-1}$.
Since the seminal work of Tibshirani_Lasso_1996, the LASSO has become the default method for estimating the high-dimensional linear regression model, owing to its ability to perform model selection and coefficient estimation simultaneously. It is well known that the LASSO estimator, as a penalized least-squares estimator, can be interpreted as a Bayesian point estimator---specifically, the posterior mode under independent Laplace priors on the regression coefficients (see Remark (ref)). However, when one considers the entire posterior distribution rather than only its mode, the Laplace prior alone fails to induce sparsity: the resulting posterior remains non-sparse and does not contract toward the true parameter at the optimal rate, unlike its mode. As shown in Theorem 7 of Castilloetal_BayesianSparse_2015, the full Bayesian posterior corresponding to the Laplace prior lacks desirable asymptotic properties.\footnote{Castilloetal_BayesianSparse_2015 remark that “the LASSO is essentially non-Bayesian, in the sense that the corresponding full posterior distribution is a useless object” (p.1988). In practice, sampled posterior distributions of regression coefficients under the Laplace prior are indeed non-sparse Baietal_ReviewSpikeSlab_2021.}
Our method is inspired by the frequentist literature on debiased inference ZhangZhang2014CI, which applies a debiasing step to the LASSO estimator. Let $\hat{\beta}^{\text{Pilot}}_n$ denote any pilot estimator of $\beta_0$, and let $\hat{\Theta}_n$ be a pilot estimator of the precision matrix $\Theta_0$. The celebrated debiased inference procedure is based on the following construction:
It is well established that the individual coordinates of the debiased estimator $\hat{\beta}_n^{\text{Debias}}$ are asymptotically normal, enabling valid statistical inference.
Our debiased Bayesian inference procedure departs from the frequentist approach in several key respects. The Bayesian modeling perspective treats the unknown regression coefficient as the random vector $\beta$. We begin by assigning a sparsity-inducing prior, such as a spike-and-slab prior or a horseshoe-type prior, and computing the (approximate) initial posterior distribution of $\beta$, denoted by $\Pi_{\beta}(\,\cdot\,|\{Y_i,X_i\}_{i=1}^n)$. Under sparsity, the resulting posterior achieves (near-)optimal contraction rates Castilloetal_BayesianSparse_2015,GaoVaartZhou2020General,SongLiang_NearlyOptimal_2023. We then debias the entire posterior distribution using a correction term analogous to the frequentist construction. In this step, however, we replace the uniform weights $1/n$ with Bayesian bootstrap weights. Specifically, the Bayesian bootstrap weights are defined by normalizing independent standard exponential random variables (independent of $\{Y_i,X_i\}_{i=1}^n$ and the prior of $\beta$):
In summary, we study the debiased posterior of the form:
where $\beta\sim \Pi_{\beta}(\,\cdot\,|\{Y_i,X_i\}_{i=1}^n)$ and $\hat{\Theta}_n$ is as introduced above. The pseudo-code for generating draws from the debiased posterior is provided in Algorithm (ref).
Several remarks are in order. First, the weight replacement is essential. To see this, consider the case $p<n$, where one may take $\hat{\Theta}_n =(\sum_{i=1}^{n} X_{i}X_{i}^\intercal/n)^{-1}$. Without the weight replacement, the expression simplifies to $\tilde{\beta} = (\sum_{i=1}^{n}X_{i}X_{i}^{\intercal})^{-1}\sum_{i=1}^{n}X_{i}Y_{i}$, which exhibits no posterior uncertainty. Second, the debiasing step introduces no additional computational burden. Since the Bayesian bootstrap weights are independent of the posterior draws of $\beta$, steps (a) and (b) of Algorithm (ref) can be parallelized efficiently to save time. Moreover, the one-time computation of $\hat{\Theta}_n$ can be carried out simultaneously, which remains computationally efficient as known in the frequentist literature (see the additional discussion in Section (ref)). Third, step (c) of Algorithm (ref) admits a natural interpretation as a Bayesian version of the residual bootstrap. Specifically, consider the case $p<n$ with $\hat{\Theta}_n =(\sum_{i=1}^{n} X_{i}X_{i}^\intercal/n)^{-1}$. For a given posterior draw $\beta^b$, define bootstrap residuals $\varepsilon_i^{b\ast}$ as $\varepsilon_i^{b}=Y_i-X_i^\intercal \beta^b$ weighted by the Bayesian bootstrap weights $W_{ni}^b$. Construct the pseudo responses $Y_{i}^{b\ast} = X_i^{\intercal}\beta^b + \varepsilon^{b\ast}$ and regress $Y_{i}^{b\ast}$ on $X_i$. The resulting estimator coincides with $\tilde{\beta}^{b}$. Thus, each debiased posterior draw $\tilde\beta^b$ can be viewed as the regression coefficient obtained from a “Bayesian residual bootstrap,” where residuals are reweighted via Bayesian bootstrap rather than resampled as in the classical residual bootstrap. While other bootstrap weights can be employed, the Bayesian bootstrap weights provide a natural Bayesian interpretation (see Section (ref)).
Our debiased Bayesian procedure provides simultaneous point estimation and uncertainty quantification. For $j=1,\ldots,p$, let $\beta_{0,j}$ and $\tilde{\beta}^{b}_j$ be the $j$th coordinate of $\beta_{0}$ and $\tilde{\beta}^{b}$ respectively. For $0<\alpha<1$, a $100(1-\alpha)\%$ credible set for the regression coefficient $\beta_{0,j}$ is defined as:
where $\hat{c}_{n,j}(\alpha,B)$ denotes the $\alpha$th quantile of the posterior draws $\{\tilde{\beta}_j^b : b=1,\ldots,B\}$. The Bayesian point estimator (posterior mean) is obtained by averaging the simulation draws: $\bar{\tilde \beta}_j = \sum_{b=1}^B \tilde{\beta}_j^b/B$ for each $j=1\ldots, p$.
Our paper aims not only to present a comprehensive theoretical development, but also to demonstrate the empirical benefits through concrete empirical applications. Before presenting the theoretical properties of our debiased Bayesian inference method, we first demonstrate its practical utility through two empirical exmaples adapted from GiannoneLenzaPrimiceri2021Sparsity. In both cases, our goal is not to determine the overall sparsity of the regressors, but rather to evaluate the significance of a particular regressor—an objective well suited to our inference procedure.
For each application, we begin by presenting standard Bayesian inference results for the coefficient of interest using both a spike-and-slab prior and a horseshoe-type prior. The prior specifications and hyperparameter settings follow those used in our Monte Carlo simulations. For the spike-and-slab prior, we obtain the posterior distribution via the variational Bayesian approximation described in RaySzabo2022Variational, drawing 8,000 samples from the resulting approximate posterior. The posterior under the horseshoe prior is generated using the MCMC algorithm of KimLeeGupta2020BayesianSC, with a total of 16,000 draws and the first 8,000 discarded as burn-in. In both applications, the spike-and-slab posterior collapses to a point mass at zero. Under the horseshoe prior, the posterior exhibits greater dispersion, though in the second application the density remains tightly concentrated near zero. In both cases, the 95% credible intervals include zero, implying that the effects are insignificant at the $5\%$ level.
To implement the debiasing adjustment, we estimate the precision matrix using the nodewise LASSO regression method of VandeGeer_OnAsymptotically_2014, adopting the default tuning parameters from the R package hdi DBMM2015hdi. For the Bayesian bootstrap, we generate 8,000 additional sets of bootstrap weights in parallel with the posterior draws of the regression coefficients. After applying the debiasing step, the posterior densities become approximately Gaussian, providing empirical support for our theoretical results. Overall, the debiased inference reveals a significant effect in the first application but no significant effect in the second.
The influential work of Barro1991Growth initiated a long-standing debate over what drives long-term economic growth across countries. Researchers have identified numerous potential predictors, many of which are included in the original data set compiled by Barro1991Growth. In line with BelloniChernozhukovHansen2013InferHigh, we use this data set to analyze average GDP growth between 1960 and 1985 for 90 countries. The data set contains 60 candidate predictors, covering pre-1960 measures of a wide range of socio-economic, institutional, and geographical factors. Using standard Bayesian inference methods, we find that all of the posterior results are statistically insignificant. However, applying our debiased Bayesian approach uncovers a significant negative association between the initial GDP level and subsequent growth. This suggests the presence of catch-up effect, everything else equal, which is in line with the neoclassical economic growth theory.
Using U.S. state-level data, DonohueLevitt2001Abortion identified a strong relationship between the legalization of abortion following the 1973 Roe v.\ Wade decision and the subsequent decline in crime rates. In their analysis, the dependent variable is the change in log per capita murder rates between 1986 and 1997 across states. This outcome is regressed on a measure of the effective abortion rate—which is always included as a predictor alongside twelve year dummy variables—and a collection of controls capturing alternative determinants of crime, such as the number of police officers and prisoners per 1,000 residents, among others. BelloniChernozhukovHansen2014HighATE extended the set of controls by incorporating these variables in various forms, including levels, differences, squared differences, cross-products, initial conditions, and interactions with linear and quadratic time trends, resulting in a database of 284 variables and 576 observations. When analyzing this expanded model, the coefficient on the abortion rate is consistently found to be insignificant across all methods. Note that the horseshoe induced posterior places almost all of its posterior mass at zero, which is close to the degenerate (at zero) posterior of the spike-and-slab. Therefore, these two debiased posteriors—whether based on the spike-and-slab or horseshoe priors—are largely indistinguishable in this application. In this example, standard Bayesian approaches place (nearly) all posterior mass at zero, while our debiased Bayesian procedure yields more nuanced posterior distributions that may better capture estimation uncertainty.
Various features of posterior distributions are used by empirical researchers for frequentist type inference. In particular, regions with high posterior probabilities, known as credible sets, are often interpreted as confidence sets. The frequentist validity of the debiased posterior distributions serves as a theoretical justification for the proposed Bayesian procedure. In this section, we present the main theoretical results, beginning with a generic BvM theorem formulated under a set of high-level conditions. The result applies broadly and is not restricted to specific prior choices.
For the reader’s convenience, we introduce notation that will be used throughout the paper. Let $I_p$ denote the $p\times p$ identity matrix, $e_j$ its $j$th column, and $E_J$ the matrix formed by any $J$ columns of $I_p$. For $1 \leq q < \infty$, we write $\|v\|_{q}$ for the standard $\ell_q$-norm of a vector $v$, that is, $\|v\|_{q}: = (\sum_j |v_j|^{q})^{1/q}$. We use $\|v\|_{0}$ to denote the number of nonzero entries of $v$, and $\|v\|_{\infty}$ for its sup-norm. For a matrix $A$, $\|A\|_{q}$ denotes the induced (operator) $\ell_q$-norm for $1\leq q\leq \infty$, and $\|A\|_{\max}$ the elementwise sup-norm. For two positive sequences $a_n$ and $b_n$, we write $a_n\lesssim b_n$ if $a_n\leq C b_n$ for some constant $C$ and all $n$, and $a_n\asymp b_n$ if $a_n\lesssim b_n$ and $b_n \lesssim a_n $. Let $Z^{(n)} := \{Y_i,X_i\}_{i=1}^n$ denote the data and $W^{(n)} := \{W_{ni}\}_{i=1}^n$ the Bayesian bootstrap weights. As previously, $\mathbb{E}_0[\cdot]$ indicates that the expectation is evalauted with respect to the distribution of $Z^{(n)}$ under $\beta=\beta_0$, while $\Pi_W(\,\cdot\,|Z^{(n)})$ denotes the conditional distribution of $W^{(n)}$ given $Z^{(n)}$. Since $W^{(n)}$ is independent of $Z^{(n)}$, this conditional distribution does not depend on $Z^{(n)}$. For a positive sequence $r_n$, the asymptotic symbols $o_{P_0}(r_n)$ and $O_{P_0}(r_n)$ are defined with respect to the same underlying probability measure $P_0$. The sub-Gaussian norm of a random variable $Z$, denoted by $\|Z\|_{\psi_2}$, is defined as $\|Z\|_{\psi_2}:=\inf\{t>0: \mathbb E_{0}[\exp(Z^2/t^2)]\leq 2\}$. For a random vector $Z$, its sub-Gaussian norm is defined as $\|Z\|_{\psi_2}:=\sup _{\|x\|_2=1}\| Z^{\intercal}x\|_{\psi_2}$.
We begin by establishing the result under high-level assumptions that specify the posterior contraction rate of $\Pi_{\beta}(\,\cdot\,| Z^{(n)})$, the convergence rate of the precision matrix estimator $\hat{\Theta}_n$, a frequentist point estimator that centers the debiased posterior distribution, and suitable moment and regularity conditions.
Assumption (ref) imposes a contraction rate on the initial posterior distribution of $\beta$. We later provide primitive conditions under which this assumption holds for both spike-and-slab priors and horseshoe-type priors Castilloetal_BayesianSparse_2015, SongLiang_NearlyOptimal_2023. Importantly, we allow for either the exact posterior or an approximate posterior for $\beta$. The exact posterior arises directly from the Bayes’ rule, given a prior and the likelihood function. In low-dimensional settings, standard MCMC algorithms can deliver fairly accurate approximations to this exact posterior. In contrast, in high-dimensional problems, it may no longer be feasible. For example, the point-mass spike-and-slab prior is often considered theoretically ideal for sparse Bayesian problems. However, exploring the full posterior over the entire model space using point-mass spike-and-slab priors can be computationally prohibitive, because of the combinatorial complexity of updating the discrete indicators whether to include each variable or not. Recently, the variational Bayesian approximation has become increasingly popular where one relies on an approximate posterior that minimizes the Kullback–Leibler divergence to the true posterior within a restricted family, such as the mean-field variational family. We provide sufficient conditions to ensure that Assumption (ref) holds for mean-field variational approximations under spike-and-slab priors RaySzabo2022Variational.
Assumption (ref) specifies the convergence rate of the precision matrix estimator. When the dimensionality $p$ is fixed, this condition is trivially satisfied by the plug-in estimator $\hat{\Theta}_n = \hat{\Omega}_n^{-1}$ with $\gamma_n = 0$ and $\delta_n = n^{-1/2}$ under mild moment assumptions. In high-dimensional settings ($p \to \infty$), we will present primitive conditions under which the assumption holds when $\hat{\Theta}_n$ is obtained from the nodewise LASSO regressions of VandeGeer_OnAsymptotically_2014 or the CLIME approach of Caietal_ConstrainedSparse_2011. Assumption (ref) requires the existence of a frequentist point estimator that serves as the center of the debiased posterior. It is satisfied by the OLS estimator when $p$ is fixed and by the debiased LASSO estimator in high-dimensional models VandeGeer_OnAsymptotically_2014 under the primitive conditions provided in Section (ref). Besides LASSO, other types of pilot estimators have also been studied in the literature; see Remark (ref). Finally, Assumption (ref) imposes sub-Gaussianity on both the regressors and the error terms. These regularity conditions are standard in the high-dimensional inference literature and ensure well-behaved concentration of sample quantities.
The debiased posterior distribution is jointly determined by the initial posterior distribution of $\beta$ and the distribution of the Bayesian bootstrap weights. Specifically, the debiased posterior distribution is determined by
where $\Pi_{\beta}(\,\cdot\,|Z^{(n)})$ and $\Pi_W(\,\cdot\,|Z^{(n)})$ are, respectively, the initial posterior distribution of $\beta$ and the distribution of the Bayesian bootstrap weights conditional on the data. For each $j=1,\ldots,p$, let $\mathcal{L}_{\Pi}\big(e_{j}^{\intercal}\sqrt{n}(\tilde\beta-\hat{\beta}_{n})\mid Z^{(n)}\big)$ denote the posterior law of $e_{j}^{\intercal}\sqrt{n}(\tilde\beta-\hat{\beta}_{n})$ given $Z^{(n)}$. We study the weak convergence of these posterior laws, which we measure using the bounded Lipschitz distance $d_{BL}$. For probability measures $P$ and $Q$ on $\mathbf{R}^k$, the bounded Lipschitz distance is defined as
where the bounded Lipschitz norm of a measurable function $f:\mathbf{R}^k\to\mathbf{R}$ is given by
Theorem (ref) establishes a coordinatewise BvM result for the debiased posterior. In particular, the asymptotic variance $\sigma^{2}_{0,j}$ for the $j$th coordinate of $\sqrt{n}(\tilde{\beta}-\hat{\beta}_{n})$ coincides with the asymptotic variance for the centering point estimator in Assumption (ref) under general, possibly heteroskedastic, errors. The rate conditions $\sqrt{n}\epsilon_n\gamma_{n}\to0$ and $\sqrt{n}(\epsilon_n\|\Theta_0\|_\infty+ \delta_n)(\sqrt{\log p/n} + \log^2 n\log^{2} p/n)\to 0$ ensure two key properties: (i) the bias stemming from the initial posterior and the error from estimating the precision matrix are asymptotically negligible; and (ii) the residual posterior uncertainty, after the debiasing correction, is asymptotically Gaussian. Importantly, the theorem accommodates a high-dimensional regime in which the dimensionality $p$ may grow exponentially with the sample size $n$. As will be shown later, under spike-and-slab or horseshoe priors for $\beta$ and nodewise LASSO or CLIME estimation for $\Theta_0$, the quantities $\epsilon_n$, $\gamma_n$, and $\delta_n$ can all be of order $\sqrt{\log p/n}$ if the sparsity levels of $\beta_0$ and $\Theta_0$ are bounded (implying $\|\Theta_0\|_\infty$ is bounded). Hence, the result remains valid as long as $\log^{5/2} p$ grows at most linearly with $n$, up to logarithmic factors.
The BvM theorem enables the construction of marginal credible intervals for any coordinate of the regression coefficient from the debiased posterior. For $0<\alpha<1$, let $q_{n,j}(\alpha)$ denote the $\alpha$th quantile of the debiased posterior distribution for the $j$th coefficient---that is, the $\alpha$th conditional quantile of $\tilde{\beta}_{j}$ (the $j$th coordinate of $\tilde{\beta}$) given $Z^{(n)}$. We study the asymptotic frequentist coverage of the $100(1-\alpha)\%$ credible intervals, defined for $j=1,\ldots,p$ as
In practice, $\mathcal{C}_{n,j}(\alpha)$ is infeasible but can be approximated by its simulated counterpart $\mathcal{C}_{n,j}(\alpha,B)$ defined in (ref) by taking a large $B$.
This corollary provides an important frequentist implication of Theorem (ref): the Bayesian credible intervals derived from the debiased posterior achieve asymptotically correct frequentist coverage for each regression coefficient. In other words, for large samples, the posterior-based uncertainty quantification coincides with the frequentist notion of confidence intervals, thereby unifying Bayesian and frequentist inference in high-dimensional settings.
We next consider simultaneous inference for a growing number of coefficients. Without loss of generality, we focus on the first $J$ coordinates. For this purpose, we strengthen Assumption (ref)(iv) to the following condition.
Assumption (ref)(i) is satisfied if, for example, $\mathbb{E}_0[\varepsilon_{i}^{2}|X_i]$ is constant and the eigenvalues of $\Omega_0$ are bounded away from zero. Assumption (ref) (ii) holds with $\nu_J = O(\log J)$ when the random vectors $\Theta_0X_{i}\varepsilon_i$ are sub-exponential.
Theorem (ref) extends the marginal BvM result in Theorem (ref) to simultaneous inference on a growing subset of regression coefficients. Specifically, the result establishes joint asymptotic normality of the debiased posterior for any subvector of $\tilde{\beta}$ of size $J$, provided that $J$ increases with the sample size $n$ at a suitable rate. This rate is derived from the strong approximation result for the exchangeable bootstrap in FangSantosShaikhTorgovitsky2023LP. The condition $\nu_J^{6}J^{3}\log^{9}(1+J)/n \to 0$ captures the complexity of the simultaneous inference problem, balancing the growth of the subset dimension against sample size and moment bounds. If we assume enough moment restrictions on $\|E_{J}^{\intercal}\Theta_0X_{i}\varepsilon_i\|_{\infty}$ (e.g., $\mathbb E_0[\|E_{J}^{\intercal}\Theta_0X_{i}\varepsilon_i\|^{4}_{\infty}]$ being bounded uniformly in $n$ and $p$), we can obtain Theorem (ref) under $J/n^2\to\infty$ (up to logarithmic factors) by Theorem (ref).
We now verify the high-level conditions in Assumptions (ref)-(ref) by introducing a set of primitive conditions. For concreteness, we focus on a benchmark specification in which the regression coefficients are assigned a spike-and-slab prior, the precision matrix is estimated via the nodewise LASSO procedure, and the centering point estimator is constructed by the deibased LASSO. The general theoretical framework, however, extends to global–local shrinkage priors, including the horseshoe family, and other estimation approaches for the precision matrix, including the CLIME approach. Detailed results for these additional cases are provided in the Appendix.
We follow the construction of Castilloetal_BayesianSparse_2015, who specify the spike-and-slab prior through the following hierarchical mixture for $\beta = (\beta_1,\ldots,\beta_p)^{\intercal}$,
where $\delta_0$ denotes a point mass at zero, and $\psi(\beta_j|\lambda) = \frac{\lambda}{2}\exp(-\lambda|\beta_j|)$ is the Laplace density with hyperparameter $\lambda>0$. The mixing weight $r$ controls the proportion of nonzero coefficients\footnote{If $r$ were taken to be deterministic, it would represent the expected number of nonzero coefficients a priori.} and thus governs the degree of model sparsity. Following Castilloetal_BayesianSparse_2015, we assign a Beta hyper-prior
with hyperparameter $u>1$, large values of which favor sparse models by placing most of the prior mass on small values of $r$.
The point mass in the spike-and-slab prior is designed for explicit variable selection. However, the posterior under this prior places mass over all $2^p$ candidate models, which is computationally prohibitive in high dimensions. We therefore adopt a variational Bayesian approximation scheme. The variational posterior is defined as the best approximation of the posterior within a given class. The crux is to choose this class both sufficiently rich to approximate the exact posterior well, while at the same time also simple enough so that the minimization problem can be efficiently solved. For this purpose, we consider the following mean-field family from RaySzabo2022Variational for $\mu = (\mu_1,\ldots,\mu_p)^{\intercal}$, $\sigma^2 = (\sigma^2_1,\ldots,\sigma^2_p)^{\intercal}$, and $\gamma = (\gamma_1,\ldots,\gamma_p)^{\intercal}$,
where $\gamma_j$ plays the role of a variational inclusion probability. The resulting variational Bayesian approximate posterior is defined as the minimizer of the Kullback-Leibler divergence with respect to the exact posterior:
where $\mathrm{KL}(\,P\,\|\,Q\,):=\int \log (dP/dQ) dP$ for two probability distributions $P$ and $Q$, and $\Pi_{\beta}(\,\cdot\,|Z^{(n)})$ is the exact posterior. Note that the above mean-field class enforces substantial independence in the approximated posterior. By doing so, it significantly reduces the model complexity, because there are only $p$ inclusion variables $\gamma_j$ to consider, rather than a total of $2^p$ models which the exact posterior puts mass on.
The mean-field class approximates the posterior distribution by a product-form distribution, thereby ignoring dependence among different coordinates. Nevertheless, it has been shown to achieve the desired rate of contraction to the true parameter RaySzabo2022Variational, which is essential for our high-level conditions to ensure the asymptotic normality of the debiased posterior. Without the debiasing step, however, the variational posterior generally fails to deliver asymptotically correct coverage, even in low-dimensional parametric settings WangBlei2019VB. This highlights the versatility of our debiasing proposal from a complementary perspective. The trade-off for relaxing the exactness of the posterior distribution lies in the substantial computational gains. Component-wise coordinate-ascent variational inference algorithms for $\tilde{\Pi}_{\beta}(\,\cdot\,|Z^{(n)})$ are given in RaySzabo2022Variational. Beyond the mean-field family, our arguments also extend to other closely related variational classes that allow dependence among the nonzero coordinates (see Equation (9) in RaySzabo2022Variational), yielding analogous theoretical guarantee.
Turning to Assumption (ref), we consider the nodewise LASSO regression following VandeGeer_OnAsymptotically_2014. Let $X_{j,i}$ denote the $j$th component of $X_{i}$ and let ${X}_{-j,i}\in\mathbf{R}^{p-1}$ be the subvector of $X_{i}$ excluding its $j$th entry. For each $j=1,\ldots, p$, we define the following LASSO regression:
with the corresponding residual variance given by
where $\lambda_{j}>0$ is a tuning parameter. Writing $\hat{\theta}_j=(\hat{\theta}_{j, 1},\ldots, \hat{\theta}_{j, j-1}, \hat{\theta}_{j, j+1}, \ldots, \hat{\theta}_{j, p})$, we construct the estimated precision matrix $\hat{\Theta}_n$ as
The estimator $\hat\Theta_n$ is motivated by noting that the restriction $\Omega_0^{-1}\Omega_0=I_p$ amounts to first order conditions of $p$ linear projections after parametrizing $\Omega_0^{-1}$ in the same structure as $\hat\Theta_n$.
We choose the centering point estimator as the debiased LASSO estimator. For $\hat{\Theta}_n$ given in (ref)-(ref), the debiased LASSO estimator is given by
where $\hat{\beta}_n^{\text{LASSO}}$ is a LASSO estimator given by
where $\rho>0$ is a tuning parameter.
We adopt the following standard notation in analyzing the high-dimensional regression. For a vector $\beta=(\beta_1,\cdots,\beta_p)^{\intercal}\in\mathbf{R}^p$ and a subset $S\subset\{1,\cdots,p \}$ of indices, let $\beta_S$ be the vector $(\beta_j)_{j\in S}\in \mathbf{R}^{|S|}$, where $|S|$ is the cardinality of $S$. Let $S_{\beta}:=\{ j\in \{1,\cdots,p \} :\beta_{j}\neq 0\}$ denote the index set of nonzero components of $\beta$, and $s_{\beta}:= |S_{\beta}|$ denote the sparsity level of $\beta$. Recall that $\hat{\Omega}_n= \sum_{i=1}^{n}X_iX_i^\intercal/n$. Define $\nu(\hat{\Omega}_n):=\max_{1\leq j\leq p}\sqrt{e_{j}^{\intercal}\hat{\Omega}_n e_{j}}$ as the square root of the largest diagonal element of $\hat{\Omega}_n$.
All the assumptions are standard in the literature. Following Castilloetal_BayesianSparse_2015, we assume $\varepsilon_i \sim N(0,1)$ for simplicity. In the case of an unknown variance, $\varepsilon_i \sim N(0,\sigma_0^2)$, one may rescale the data using an estimate of $\sigma_{0}^2$, thereby adopting an empirical Bayes approach. Alternatively, a fully Bayesian treatment can be employed by assigning a prior to $\sigma_{0}^2$, for instance, an inverse-Gamma prior. In contrast to Castilloetal_BayesianSparse_2015 and RaySzabo2022Variational, we provide primitive conditions for their high-level conditions on compatibility numbers. Since $\|\Theta_0\|_{\infty }=O(\max_{1\leq j\leq p}\sqrt{s_{\Theta_0,j}})$, Assumption (ref) involves the sparsity levels of $\beta_0$ and $\Theta_0$ and the dimension $p$. However, the requirements are slightly different from those in the frequentist inference literature. For example, VandeGeer_OnAsymptotically_2014 require $s_{\beta_0}\log p/\sqrt{n}\to 0$ and $\max_{1\leq j\leq p}s_{\Theta_0,j}\log p/\sqrt{n}\to 0$. When the sparsity levels of $\beta_0$ and $\Theta_0$ are bounded, Assumption (ref) requires $\log^{5/2}p = o(n/\log^2 n)$, which is slightly more retrictive than $\log p = o(\sqrt{n})$ required in VandeGeer_OnAsymptotically_2014. As in the literature on frequentist inference, $p$ is allowed to grow exponentially with $n$.
Our proposal can be seen as performing one-step update based on the identity $\beta_0 = \chi(\beta_0,\Theta_0,F_0)$, where $F_0$ denotes the distribution of $(Y_i,X_i)$, and $\chi$ is a mapping from the space of $(\beta_0,\Theta_0,F_0)$ to that of $\beta_0$ defined by
We consider the posterior law of $F$, denoted by $\Pi_{F}(\,\cdot\,|Z^{(n)})$, to be the Bayesian bootstrap law, which can be viewed as the posterior of $F$ in a nonparametric Bayesian model with a Dirichlet process prior having a degenerate (zero-mass) base measure Rubin_BayesianBootstrap_1981. Then, (ref) can be expressed as
where $(\beta,F)|Z^{(n)} \sim \Pi_{\beta}(\,\cdot\,|Z^{(n)}) \times \Pi_{F}(\,\cdot\,|Z^{(n)})$. Informallly speaking, the debiased posterior is obtained by plugging in the frequentist estimator $\hat{\Theta}_n$ for $\Theta_0$, the initial posterior for $\beta$, and the Bayesian bootstrap for $F$. Note that the evaluation of $\Pi_{F}(\,\cdot\,|Z^{(n)})$ boils down to $\Pi_{W}(\,\cdot\,|Z^{(n)})$, whose randomness involves the weighted exponential random variables. Since the two posteriors are conditionally independent, the debiased posterior distribution arises from the product of the initial posterior distribution of $\beta$ and the posterior distribution of $F$. Next, we discuss the comparison of our approach with a number of recent proposals in the literature.
Our debiased posterior aligns with the one-step posterior framework of YiuFongHolmesRousseau2023Semiparametric, which studies scalar parameters of interest in general semiparametric models. Nonetheless, our proposal is distinct in several important ways. First, we consider a high-dimensional regime where the number of covariates may exceed the sample size, and we allow for simultaneous inference for a growing number of coefficients. The one-step correction term in YiuFongHolmesRousseau2023Semiparametric involves the influence function perturbed by the Bayesian bootstrap weights. In a linear regression model, the $j$th regression coefficient $\beta_{0,j}$ has the influence function $(y,x^{\intercal})^{\intercal}\mapsto e_{j}^{\intercal}\Theta_{0} x(y-x^\intercal \beta_0)$ at the truth. This representation permits the application of their one-step posterior in fixed-dimensional settings. In contrast, our framework explicitly employs sparsity-inducing priors in high-dimensional regimes. Moreover, the debiasing construction applies to the entire vector of regression coefficients, allowing the number of covariates to exceed the sample size while both $n$ and $p$ diverge. The resulting general BvM theorems enable inference on subvectors of regression coefficients whose dimension can grow with $(n,p)$.
Second, we estimate $\Theta_0$ via a frequentist method rather than a Bayesian one, which substantially improves computational efficiency. This necessitates addressing non-trivial technical challenges and developing new theoretical results. Suppose the covariates $\{X_i\}_{i=1}^n$ are i.i.d. draws from ${N}(0,\Theta_0^{-1})$. The likelihood of $\{X_i\}_{i=1}^n$ as a function of the precision matrix $\Theta_0$ is
A natural Bayesian approach would place a sparsity-inducing prior, such as a spike-and-slab prior or a horseshoe-type prior, on the entries of $\Theta_0$ Atchade2019Graph. However, computing the posterior distribution of a high-dimensional precision matrix is prohibitive in practice. Our proposal instead adopts a more pragmatic strategy: we simply plug in a frequentist estimator. In particular, one may compute $\hat{\Theta}_n$ either by running $p$ nodewise LASSO regressions or by solving $p$ linear programming problems in CLIME, both of which are computationally efficient. This plug-in approach is feasible because our procedure only requires $\hat{\Theta}_n$ to satisfy certain convergence rate conditions. Since our primary inferential target is the regression coefficients $\beta_0$, and $\Theta_0$ only plays a supporting role, this frequentist plug-in approach is both theoretically valid and practically effective.
Third, we accommodate a broad class of bootstrap weights beyond the Bayesian bootstrap weights; while the Bayesian bootstrap weights provide a natural Bayesian interpretation, alternative weights can be incorporated. A closer look at proofs indicates that other common exchangeable weights may also be employed. The Bayesian bootstrap law corresponds to the posterior law when the empirical distribution of $\{Y_i,X_i\}_{i=1}^n$ is equipped with a Dirichlet process prior with degenerate (zero-mass) base measure. For more general base measures, however, the resulting posterior law becomes substantially more cumbersome both theoretically and computationally; see Equation (2.1) in RayVaart2021BvM for a representation of the posterior in this case.
In this section, we evaluate the finite-sample performance of the proposed debiased Bayesian inference method through a series of Monte Carlo experiments. Specifically, we examine the frequentist coverage of the Bayesian credible sets and assess estimation accuracy in terms of bias and root mean squared error (RMSE). Beyond the stylized linear model with standard Gaussian errors presented in Section (ref), we further demonstrate the robustness and effectiveness of our approach across a range of more complex data-generating processes. We highlight that current results are all based on the default choices of hyper-parameters and tuning parameters in the existing R packages, which do not require additional interventions from users.
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9956444 0.9379778 0.9522667 1 0.95 0.1120000 0.9060000 0.8870000 2 0.95 0.0750000 0.9050000 0.9100000 3 0.95 0.0500000 0.8940000 0.8820000 4 0.95 0.0670000 0.9070000 0.6290000 5 0.95 0.0820000 0.9080000 0.8400000 }\pfnoheter
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9947778 0.9392 0.9522222 1 0.95 0.0620000 0.9210 0.9480000 2 0.95 0.0530000 0.9090 0.9270000 3 0.95 0.0570000 0.8880 0.9110000 4 0.95 0.0600000 0.9020 0.6450000 5 0.95 0.0730000 0.9250 0.8710000 }\pfnohomochi
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9913333 0.9454444 0.9519111 1 0.95 0.3070000 0.8910000 0.9250000 2 0.95 0.4590000 0.8810000 0.8470000 3 0.95 0.4670000 0.8850000 0.7230000 4 0.95 0.5210000 0.8990000 0.6220000 5 0.95 0.5810000 0.9010000 0.7560000 }\pfnohomon
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9998632 0.9392316 0.9528421 1 0.95 0.0530000 0.9210000 0.9080000 2 0.95 0.0560000 0.9310000 0.9440000 3 0.95 0.0780000 0.8900000 0.8970000 4 0.95 0.0820000 0.9060000 0.6760000 5 0.95 0.0940000 0.9290000 0.8530000 }\ponoheter
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9999158 0.9371368 0.9528105 1 0.95 0.0390000 0.9400000 0.9590000 2 0.95 0.0460000 0.9130000 0.9290000 3 0.95 0.0470000 0.9020000 0.9310000 4 0.95 0.0400000 0.9020000 0.7100000 5 0.95 0.0530000 0.9250000 0.8700000 }\ponohomochi
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9995158 0.9480526 0.9601368 1 0.95 0.1520000 0.9210000 0.9330000 2 0.95 0.3130000 0.9000000 0.8700000 3 0.95 0.3330000 0.8840000 0.7330000 4 0.95 0.4080000 0.8970000 0.5900000 5 0.95 0.5460000 0.9140000 0.7280000 }\ponohomon
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9999333 0.9378872 0.953159 1 0.95 0.0390000 0.9270000 0.921000 2 0.95 0.0430000 0.9220000 0.944000 3 0.95 0.0570000 0.8740000 0.918000 4 0.95 0.0500000 0.9140000 0.728000 5 0.95 0.0530000 0.9290000 0.885000 }\ptonoheter
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.999959 0.9368205 0.9530615 1 0.95 0.018000 0.9210000 0.9540000 2 0.95 0.029000 0.9190000 0.9410000 3 0.95 0.046000 0.8840000 0.9340000 4 0.95 0.045000 0.9170000 0.7540000 5 0.95 0.042000 0.9280000 0.8780000 }\ptonohomochi
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9997385 0.9457128 0.9647282 1 0.95 0.1190000 0.9230000 0.9440000 2 0.95 0.2240000 0.8780000 0.8640000 3 0.95 0.2200000 0.8500000 0.7250000 4 0.95 0.3010000 0.8790000 0.5640000 5 0.95 0.4730000 0.9240000 0.7400000 }\ptonohomon
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9998444 0.9439556 0.9522 1 0.95 0.1900000 0.9210000 0.8590 2 0.95 0.2890000 0.9180000 0.9380 3 0.95 0.3330000 0.9330000 0.9500 4 0.95 0.3780000 0.9250000 0.9410 5 0.95 0.5790000 0.9450000 0.9680 }\CovpfnoheterSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9998889 0.9355111 0.9498444 1 0.95 0.0800000 0.9090000 0.9370000 2 0.95 0.1220000 0.9240000 0.9390000 3 0.95 0.1300000 0.9270000 0.9550000 4 0.95 0.1540000 0.9370000 0.9550000 5 0.95 0.2510000 0.9450000 0.9560000 }\CovpfnohomochiSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9991778 0.9604222 0.9544222 1 0.95 0.6210000 0.9200000 0.9470000 2 0.95 0.8080000 0.9120000 0.9240000 3 0.95 0.8710000 0.9230000 0.9320000 4 0.95 0.9050000 0.9210000 0.9360000 5 0.95 0.9280000 0.9450000 0.9510000 }\CovpfnohomonSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9998211 0.9430632 0.9519053 1 0.95 0.1470000 0.9070000 0.8520000 2 0.95 0.1970000 0.9230000 0.9500000 3 0.95 0.2510000 0.9220000 0.9500000 4 0.95 0.3170000 0.9270000 0.9470000 5 0.95 0.5290000 0.9490000 0.9590000 }\CovponoheterSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9999579 0.9342526 0.9499684 1 0.95 0.0460000 0.9280000 0.9580000 2 0.95 0.0700000 0.9310000 0.9560000 3 0.95 0.0790000 0.9140000 0.9370000 4 0.95 0.1110000 0.9340000 0.9530000 5 0.95 0.2010000 0.9550000 0.9610000 }\CovponohomochiSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9991053 0.9630316 0.9560211 1 0.95 0.4760000 0.8800000 0.9440000 2 0.95 0.6820000 0.9040000 0.9310000 3 0.95 0.8080000 0.9240000 0.9390000 4 0.95 0.8660000 0.9300000 0.9380000 5 0.95 0.9250000 0.9460000 0.9480000 }\CovponohomonSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9994667 0.9422359 0.9518769 1 0.95 0.1050000 0.8900000 0.8770000 2 0.95 0.1460000 0.9030000 0.9400000 3 0.95 0.1870000 0.9250000 0.9410000 4 0.95 0.2320000 0.9280000 0.9520000 5 0.95 0.4480000 0.9430000 0.9550000 }\CovptnoheterSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9998256 0.9352974 0.9508923 1 0.95 0.0280000 0.9030000 0.9450000 2 0.95 0.0430000 0.9160000 0.9470000 3 0.95 0.0590000 0.9200000 0.9530000 4 0.95 0.0680000 0.9360000 0.9580000 5 0.95 0.1570000 0.9540000 0.9650000 }\CovptnohomochiSigmaone
\pgfplotstableread{ beta alpha Bayes DebiasBayes DebiasLasso 0 0.95 0.9987333 0.9634718 0.9564872 1 0.95 0.3640000 0.8100000 0.9550000 2 0.95 0.6020000 0.8850000 0.9380000 3 0.95 0.7430000 0.9130000 0.9350000 4 0.95 0.8540000 0.9300000 0.9420000 5 0.95 0.9240000 0.9390000 0.9520000 }\CovptnohomonSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.001405958 0.0004518953 0.001585543 1 0.188266654 0.0198290763 0.006154626 2 0.379296830 0.0079773669 0.006147921 3 0.551779900 0.0223169900 0.005791037 4 0.687083060 0.0496947506 0.029428829 5 0.798797334 0.0089568229 0.008365539 }\BiaspfnoheterSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0001724629 0.021572592 0.022899255 1 0.2298244455 0.021302251 0.005174676 2 0.4558012007 0.018473234 0.002541299 3 0.6727092236 0.017276890 0.005215583 4 0.8796584906 0.036355861 0.024209296 5 1.4444430016 0.005467414 0.001221988 }\BiaspfnohomochiSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.001383397 0.005604515 0.003177875 1 0.123256509 0.019453508 0.011714049 2 0.084265179 0.003992994 0.008970469 3 0.064312110 0.008131654 0.019866299 4 0.066967318 0.014106107 0.030450883 5 0.055819965 0.006161213 0.022271943 }\BiaspfnohomonSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.002200622 0.0005885165 0.0004309773 1 0.209249474 0.0309213007 0.0058489507 2 0.433571940 0.0500124761 0.0268647271 3 0.606398851 0.0267364013 0.0079332456 4 0.727287633 0.0412250312 0.0239651555 5 0.883873170 0.0031275531 0.0058530058 }\BiasponoheterSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0004517858 0.017746747 0.014082731 1 0.2440684313 0.015509021 0.011653418 2 0.4768083808 0.025150244 0.004387045 3 0.6944520531 0.036664305 0.026167961 4 0.8928678258 0.009734311 0.024369349 5 1.5364042550 0.005458272 0.002653375 }\BiasponohomochiSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0003041703 0.003508030 0.00230409 1 0.1652696245 0.034634330 0.01369987 2 0.1777309113 0.028141826 0.02851929 3 0.1271038110 0.010364185 0.02419351 4 0.1082425540 0.013614436 0.03484569 5 0.0426113377 0.009365502 0.01871877 }\BiasponohomonSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0001799846 0.009249999 0.007015938 1 0.2165986597 0.041992712 0.011857832 2 0.4451585143 0.041864287 0.010924753 3 0.6271699454 0.043673796 0.018468610 4 0.7847656468 0.027481907 0.008362704 5 1.0399637874 0.027770277 0.028437195 }\BiasptnoheterSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.001088812 0.03254957 0.034389545 1 0.240153859 0.02675843 0.004916506 2 0.479685095 0.03712065 0.012825109 3 0.710580816 0.03511702 0.023680412 4 0.929041048 0.05068757 0.035217777 5 1.652809687 0.04848980 0.038397464 }\BiasptnohomochiSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.001279871 0.005058094 0.004698813 1 0.182678905 0.046673645 0.016240946 2 0.215813474 0.031674831 0.023118812 3 0.183858651 0.023055634 0.029313589 4 0.140857600 0.016156247 0.031925811 5 0.065958831 0.012642770 0.038299084 }\BiasptnohomonSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.01925207 0.01589084 0.09247683 1 0.22005008 0.05499430 0.04028210 2 0.41834090 0.08614162 0.12418467 3 0.64245754 0.14470681 0.16520591 4 0.79990089 0.12796052 0.32591465 5 0.40817257 0.08722685 0.18174314 }\BiaspfnoheterSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.01861526 0.02142069 0.10998463 1 0.23570318 0.05861051 0.03837186 2 0.44150229 0.07770014 0.11975520 3 0.65084910 0.13654482 0.14655078 4 0.83032609 0.12435710 0.33407180 5 0.43106131 0.09160899 0.18111704 }\BiaspfnohomochiSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.003968657 0.003084803 0.01835252 1 0.113172058 0.020778603 0.03020753 2 0.111124704 0.030065644 0.06907534 3 0.062742655 0.022823021 0.10995558 4 0.041286355 0.017870030 0.13628569 5 0.026227520 0.012266941 0.08989703 }\BiaspfnohomonSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0003854495 0.01096429 0.09335733 1 0.2331656821 0.06218585 0.03936348 2 0.4632173206 0.08336527 0.10991463 3 0.6692672165 0.15808925 0.15691825 4 0.8253575326 0.13633374 0.31972620 5 0.4103119204 0.09383610 0.18014970 }\BiasponoheterSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0000417 0.01837294 0.11257658 1 0.2421614 0.05504538 0.03360933 2 0.4672417 0.07278867 0.10560262 3 0.6924365 0.15432506 0.14150559 4 0.8859219 0.13005213 0.32324648 5 0.4421301 0.09837324 0.17942505 }\BiasponohomochiSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.003506978 0.000316633 0.01817540 1 0.184636804 0.037862667 0.03510348 2 0.248930128 0.055971252 0.07236356 3 0.212215880 0.047770395 0.12039333 4 0.123404768 0.030420800 0.15347303 5 0.041941850 0.016230821 0.10470971 }\BiasponohomonSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0007001274 0.01035257 0.09325036 1 0.2377599729 0.05398724 0.02732927 2 0.4640292922 0.07979414 0.09458665 3 0.6941966000 0.16645638 0.14223467 4 0.8925797685 0.14864868 0.30480988 5 0.4387484356 0.09187280 0.15887714 }\BiasptnoheterSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0010033 0.01765966 0.1096036 1 0.2443898 0.05233355 0.0271311 2 0.4785253 0.09218970 0.1101223 3 0.6974934 0.15740223 0.1260607 4 0.9078971 0.13117534 0.2983622 5 0.4587348 0.10312010 0.1683140 }\BiasptnohomochiSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.00540785 0.00046631 0.01902918 1 0.19907592 0.03883401 0.03311536 2 0.31594224 0.07267480 0.07249739 3 0.32860568 0.07351858 0.12647247 4 0.21528844 0.05249209 0.16446921 5 0.06759288 0.01927046 0.10592719 }\BiasptnohomonSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.08254837 0.5086305 0.5441416 1 0.25971430 0.2637163 0.2786553 2 0.46170014 0.3033568 0.3153846 3 0.65905194 0.3591929 0.3726307 4 0.84968007 0.4218454 0.4267199 5 1.29128758 0.4603910 0.4594298 }\RMSEpfnoheterSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.06209489 0.6611826 0.6778838 1 0.25160765 0.2675794 0.2839598 2 0.48733017 0.3771388 0.3906908 3 0.72631279 0.4635989 0.4665322 4 0.95027629 0.5102706 0.5193839 5 1.73324747 0.5847290 0.5860584 }\RMSEpfnohomochiSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.06815679 0.2451978 0.2804317 1 0.19176764 0.1098113 0.1133537 2 0.26271838 0.1660429 0.1713026 3 0.27672554 0.1957949 0.2049842 4 0.28410034 0.2180002 0.2322905 5 0.24572306 0.2394319 0.2590445 }\RMSEpfnohomonSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.04539301 0.5078626 0.5396947 1 0.25982750 0.2650866 0.2876093 2 0.47776810 0.2918635 0.3088945 3 0.68868893 0.3547176 0.3620187 4 0.88052546 0.4314152 0.4370849 5 1.36814845 0.4550914 0.4526235 }\RMSEponoheterSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.02811768 0.6416459 0.6581995 1 0.24828050 0.2525080 0.2779783 2 0.49103707 0.3650983 0.3790894 3 0.73751676 0.4876046 0.4917899 4 0.96280747 0.5483693 0.5549353 5 1.78887747 0.5720341 0.5695781 }\RMSEponohomochiSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0353110 0.2342876 0.2740181 1 0.2103085 0.1154194 0.1137970 2 0.3150441 0.1729335 0.1697419 3 0.3360473 0.2032149 0.2035512 4 0.3418442 0.2322440 0.2403927 5 0.2545867 0.2439282 0.2550367 }\RMSEponohomonSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.06266688 0.5109558 0.5367663 1 0.25404252 0.2516392 0.2783975 2 0.48008979 0.2982676 0.3109914 3 0.70624598 0.3689816 0.3798191 4 0.91121555 0.4225943 0.4320719 5 1.47325863 0.4628739 0.4641672 }\RMSEptnoheterSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.03904932 0.6827244 0.6942065 1 0.25422926 0.2530075 0.2775954 2 0.49637675 0.3846165 0.3996284 3 0.73710518 0.4641862 0.4659896 4 0.97600696 0.5270730 0.5407012 5 1.84099104 0.5693631 0.5718987 }\RMSEptnohomochiSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.04626209 0.2380254 0.2807367 1 0.22161738 0.1220298 0.1164194 2 0.34942491 0.1761957 0.1717340 3 0.39417550 0.2128253 0.2157517 4 0.38773882 0.2362682 0.2429289 5 0.27160848 0.2506018 0.2717188 }\RMSEptnohomonSigmaone
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.09168509 0.2163082 0.2046474 1 0.27463131 0.3054994 0.3104504 2 0.47637559 0.2479545 0.2600198 3 0.70733907 0.2712974 0.2602138 4 0.90924111 0.2656003 0.3858627 5 0.49621814 0.2501143 0.2675998 }\RMSEpfnoheterSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.08779023 0.2299180 0.2245792 1 0.26106128 0.2656319 0.2720877 2 0.48299388 0.2550444 0.2657801 3 0.71286718 0.2826652 0.2627743 4 0.92871305 0.2820502 0.4032482 5 0.50143339 0.2504466 0.2700973 }\RMSEpfnohomochiSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.02958598 0.09069449 0.08294056 1 0.20510245 0.11858631 0.11920002 2 0.23966992 0.12577262 0.13439824 3 0.21783138 0.12206385 0.15986872 4 0.19774494 0.11978397 0.17797528 5 0.11721226 0.10809705 0.14070302 }\RMSEpfnohomonSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.01113178 0.2225581 0.2142471 1 0.26707937 0.2972665 0.3013392 2 0.48454461 0.2314504 0.2482455 3 0.71034970 0.2735636 0.2609650 4 0.91087610 0.2646515 0.3805860 5 0.47610880 0.2380804 0.2650259 }\RMSEponoheterSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.00019168 0.2336082 0.2315057 1 0.25160764 0.2581977 0.2692475 2 0.48637188 0.2539778 0.2669972 3 0.72393887 0.2779722 0.2557865 4 0.94569961 0.2771307 0.3924041 5 0.48940920 0.2476144 0.2717588 }\RMSEponohomochiSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.02842996 0.09851303 0.08989901 1 0.22966045 0.11707309 0.11997834 2 0.33936494 0.12602572 0.13516457 3 0.38408620 0.13598380 0.16585410 4 0.30282318 0.12907162 0.19242741 5 0.15056226 0.11501320 0.15130191 }\RMSEponohomonSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.0154824 0.2169399 0.2146040 1 0.2563912 0.2867962 0.2976811 2 0.4870332 0.2355932 0.2488787 3 0.7215580 0.2788700 0.2581164 4 0.9449017 0.2617744 0.3629847 5 0.4860447 0.2361386 0.2545549 }\RMSEptnoheterSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.02227146 0.2315107 0.2322819 1 0.25263231 0.2660762 0.2803563 2 0.49126995 0.2557743 0.2762522 3 0.72526834 0.2825386 0.2575012 4 0.95687267 0.2664669 0.3690936 5 0.49710271 0.2577885 0.2759601 }\RMSEptnohomochiSigmatwo
\pgfplotstableread{ beta Bayes DebiasBayes DebiasLasso 0 0.03176429 0.09775338 0.09127577 1 0.23455020 0.11603385 0.11992921 2 0.38856441 0.13797691 0.13812807 3 0.48110861 0.15255231 0.17382285 4 0.40918779 0.14570092 0.20269544 5 0.19329238 0.12088637 0.15347705 }\RMSEptnohomonSigmatwo
We compare our debiased Bayes method with both standard Bayesian inference based on the spike-and-slab prior and a frequentist benchmark---specifically, the debiased inference procedure using a LASSO pilot estimator. The methods under comparison are summarized as follows.
\noindentSimulation design. We fix the sample size at $n = 100$ and vary the ambient dimensionality of the covariates by setting $p = 50, 100,$ and $200$, respectively. The true regression coefficient vector $\beta_0$ is sparse with $5$ nonzero entries. Specifically, the first five elements of $\beta_0$ are set to $(0.25, 0.5, 0.75, 1, 2)$, while the remaining coefficients are zero.
Data-generating processes. We consider six simulation scenarios designed to examine the robustness of each method:
The heteroskedastic specification in S3 is similar to that in Section 4.2 of HouMaWang2023CompositeQuantile. For S1–S3, the regressors are drawn independently as $X_i \sim N(0, \Theta_0^{-1})$, where $\Theta_0$ is a $p\times p$ diagonal matrix with entries $(1, 2, \ldots, p)$ on the diagonal.
Each simulation scenario is replicated 1,000 times.
Implementation details. For the standard Bayesian method, we employ the variational Bayes algorithm of RaySzabo2022Variational to approximate the posterior distribution. Each posterior sample consists of $B = 8{,}000$ draws from the variational Bayes approximated posterior. The spike-and-slab prior uses a Laplace slab with $\lambda = 1$ and a Beta$(1, p^u)$ hyperprior on the inclusion probability with $u = 1$. No additional tuning is performed. For the precision matrix estimator $\hat{\Theta}_n$, we adopt the nodewise LASSO regression of VandeGeer_OnAsymptotically_2014, as implemented in the R package hdi DBMM2015hdi. This estimator is used for both the Debiased Bayes and the Debiased LASSO methods to ensure comparability. Our numerical experiments suggest that the default tuning parameters in the package yield stable and robust performance. We anticipate that further optimization of hyper-parameters could potentially improve finite-sample performance, such exploration is beyond the scope of the present study.
Figures (ref) and (ref) illustrate the empirical coverage probabilities of the $95\%$ credible or confidence sets for the regression coefficients across the simulation scenarios. Specifically, Figure (ref) reports results under specifications S1 (first column), S2 (second column), and S3 (third column), with each row corresponding to a different dimensionality setting: $p = 50$, $p = 100$, and $p = 200$. On the horizontal axis, the label 0 represents the average coverage across all zero coefficients in $\beta_0$, whereas labels 1–5 correspond to the coverage probabilities for the nonzero coefficients $(0.25, 0.5, 0.75, 1, 2)$, respectively. Figure (ref) displays the analogous coverage results for specifications S4–S6, following the same layout and interpretation. We also compare the estimation accuracy of the posterior mean of the debiased posterior distribution against the standard Bayesian posterior mean and the frequentist debiased estimator. Figures (ref) and (ref) report the empirical bias of the three estimators, while Figures (ref) and (ref) present the corresponding RMSE.
We summarize our findings as follows. The standard Bayesian approach based on the spike-and-slab prior performs well in identifying the zero coefficients, as reflected by favorable average coverage, low bias, and small RMSE for these components. However, its performance deteriorates substantially for the nonzero coefficients,\footnote{By contrast, RaySzabo2022Variational considered much larger nonzero entries of $\beta_0$ (as large as 10), whereas ours are more moderate.} and becomes more fragile when the error distribution departs from the baseline homoscedastic Gaussian specification. In contrast, our proposed debiased Bayes approach consistently outperforms the standard method, achieving markedly better estimation accuracy and more reliable uncertainty quantification across all designs. Relative to the frequentist debiased LASSO benchmark, our method delivers comparable performance under the independent-covariate settings (S1–S3). Notably, in the more challenging correlated-covariate scenarios (S4–S6), the debiased Bayes procedure tends to attain substantially improved coverage, lower bias, and smaller RMSE for relatively larger coefficients. Overall, these results demonstrate the robustness and accuracy of the debaised Bayes approach across diverse conditions and the importance of debiasing standard Bayesian procedures in high-dimensional settings. It provides clear improvements over conventional Bayesian inference and exhibits complementary strengths relative to the leading frequentist benchmark.
In this paper, we develop a new debiased Bayesian inferential method for high-dimensional linear regression models. The construction resembles the frequentist debiasing step, whereas the key difference is that we correct the entire posterior distribution rather than the point estimator. Our approach is tailored to building confidence intervals based on sparsity-inducing priors such as spike-and-slab or horseshoe type priors of the regression coefficients. We establish the frequentist validity of our proposal in the general setup and also provide low-level conditions. It is straightforward to observe that our debiasing step easily extends to other parametric or semiparametric models. To demonstrate the versatility of our general methodology, we mention two possible extensions.
One may consider the generalized linear model for the conditional density function of the dependent variable given covariates as follows:
where $d$ is a convex link function. In this case, our debiased Bayesian method builds on
where $\hat{\Theta}_n$ stands for a regularized inverse of the Hessian matrix $n^{-1}\sum_{i=1}^nd^{\prime\prime} (X_{i}^{\intercal}\hat\beta_n)X_iX_i^\intercal$, given some pilot estimator $\hat{\beta}_n$.
One may also incorporate group structure which often occurs in additive models Baietal2022Group, multivariate outcome variables regression, and regressors collected over mixed frequencies MoglianiSimoni2021BMidas. For example, consider the linear regression specified by the following form:
where $\beta_{g,0}$ is an $m_g\times 1$ vector of coefficients, and $X_g$ is an $m_g\times 1$ vector of covariates corresponding to group $g = 1, \cdots,G$. The number of groups $G$ is potentially larger than the sample size $n$. In this case, our debiasing procedure is based on
where $\mathcal{E}_g$ here is the matrix formed by columns of $I_p$ such that $\mathcal{E}_g\beta=\beta_g$ with $\beta=(\beta_1^\intercal,\ldots,\beta_G^\intercal)^\intercal$ and $p=m_1+\cdots+m_G$, $X_i=(X_{1,i}^\intercal,\ldots,X_{G,i}^\intercal)^\intercal$, and the initial posterior for $\beta$ can be obtained by the group spike-and-slab prior from Baietal2022Group.
For both extensions, we believe our debiased Bayesian inference offers interesting new insights. We will defer the detailed theoretical development, as well as practical performance to some future work.