El Karoui and Purdom wrote a mathematically solid paper on how the bootstrap as a general resampling technique fails in high dimensions: https://arxiv.org/abs/1608.00696.

I think it is a very important piece of work and I wonder why it has so few citations. This led me to imagine that perhaps most practitioners are either unaware of this issue (honestly, the derivations use random matrix theory, which are already quite sophisticated) or maybe in real life datasets, the bootstrap can be used.

I would like to understand if this paper is the current state of the art in our understanding of the bootstrap. Do applied researchers still use the bootstrap when they analyze high dimensional data, and if so, do they trust it?

I personally am not aware of any general purpose method that solves this issue of bootstrapping in high dimensions. May be the real datasets really have a lot of signal in them so that the bootstrap can be applied without much trouble?

What is your take on this issue? Do you agree with the message of the paper, or do you think otherwise?

Acknowledgement: Thanks to @Silverfish for suggesting the question framing!

  • 10
    $\begingroup$ (+1) I like this question but think it is liable to get closed as too broad/unfocussed even if not opinion-based. You present several (related) questions here and narrowing its scope could help avoid that. $\endgroup$
    – mkt
    Commented May 12, 2023 at 6:11
  • 1
    $\begingroup$ Thanks. Maybe then I will just talk about the bootstrap. It’s very widely used. It’s night here and I am in bed, so I will edit the question in the morning. $\endgroup$ Commented May 12, 2023 at 6:15
  • 4
    $\begingroup$ I'd suggest one problem with the current state of this (interesting!) question is that it's very much posed as "here's my (evidence-based and not entirely subjective) opinion - what do you think/do, can you change my mind?" Almost as if it is a Q trying to be an A at the same time? I think this might be a case where "less is more" - people do usually lengthen Qs to try to get them reopened, but this one might best be cut down to 2-3 paragraphs. And by "best" I mean "most likely to elicit useful responses". Refocusing the Q would probably help crystallise what the central point of the Q is too. $\endgroup$
    – Silverfish
    Commented May 12, 2023 at 18:59
  • 2
    $\begingroup$ There are several different angles you could take here, but I think however you play it, you probably want to step back a bit from using the question body to posit your own answer. Eg if you wanted to ask about the literature more broadly, you might reframe along the lines of "What's the current state of the literature on high-dimensional bootstrap? I've read this paper, which strikes me as very solid, that casts a lot of doubt on the viability of the method. Is this regarded as the definitive take? Has the paper been challenged or suggestions for more trustworthy procedures been put forward?" $\endgroup$
    – Silverfish
    Commented May 12, 2023 at 20:21
  • 2
    $\begingroup$ Okay, I have edited the question. It is also open now, so thanks moderators! I wish to get insightful answers now! $\endgroup$ Commented May 13, 2023 at 2:48

1 Answer 1


TL;DR The article that you refer to makes things look more worse than they actually are. Their bootstrapping procedure is not a good way to apply bootstrapping. In the case of OLS there shouldn't be big problems with high dimensionality if the sample size is large. If you can not get correct results with OLS, where a correct confidence interval can be easily computed analytically, then something must be wrong with the implementation of the bootstrapping method.

It is good though to be reminded that the residuals are not the same as the errors and that we can use simulations with OLS to test (potentially wrong) implementations of bootstrapping.

Simple reproduction of the article results

The article that you refer to is performing simulations of errors by bootstrapping/resampling of the residuals. Below is a simple example that reproduces this.

The model is a linear regression with $n=500$ samples (or 250 pairs) and $p=125$ parameters. The distributions that are plotted here are just for the first parameter estimate $\hat{\beta}_1$.

Discrepancy in estimated sample variance

The third image, resampling the true errors, gives a correct indication of the sample distribution of the coefficient.

The first and second images, resampling all residuals, or resampling the pairs, have distributions with a different variance. They lead to errors in the estimates of standard errors and confidence intervals.

The reason for the discrepancy is that bootstrapping only works when the bootstrapped samples are a good representation of the true distribution. This is not the case when $p/n$ is large.

  • resampling residuals The bootstrap samples are created by simulating errors by sampling from the residuals, however the variance of the residuals is lower than the variance of the errors $$\text{Var}(r_i) \approx \left(1-\frac{p}{n}\right) \text{Var}(\epsilon)$$

  • pairwise resampling in the case of pairwise resampling the distribution is effectively a scaled binomial distribution. The variance will be

    $$\text{Var}(r_{i,paired}) \approx \frac{1}{2} \text{Var}(\epsilon)$$

    The factor $0.5$ stems from the Bernoulli variable having variance $0.25$ and the scaling is the difference of the two pairs which has variance $2\text{Var}(\epsilon)$


You should only use bootstrapping when the sampling is done from a distribution that represents the population distribution. This is not the case with paired resampling, which is sampling a distribution with half the variance, and with the residual resampling, which is sampling a distribution with smaller variance if $p/n$ is large.

The bootstrapping is often performed when a distribution of is difficult to compute. This is either the case when 1) the assumptions about the error distribution are false, or the error distribution is unknown 2) the propagation of errors is difficult to compute.

For ordinary linear regression, the second case is not an issue. The statistic is a linear sum of the data and it's sampling distribution will often approximate a normal distribution. With different cost functions the behaviour might not be too far. The problem is just to estimate the variance, and the residuals are often a good indication for this. But, one has to apply the right corrections.

The problem is more difficult in the situation where the error distribution has large tails and the variance is not easily estimated with a small sample. In this case the typical remedy is to simply gather more data. Potentially one could do an advanced semi-parametric bootstrapping by combining the residuals with a normal distribution that relates to the residuals being the errors with the estimate subtracted (the estimate being some correlated normal distribution).

Plot of reproduction

example of different variance


n = 500 # data samples
p = 125 # parameters 
m = 1000 # times resampling

### create paired data
X = matrix(rnorm(n*p/2),ncol =p)
X = rbind(X,X)
Y = rnorm(n)

solve(t(X)%*%X)[1,1] ### this is the theoretic variance

### compute main model
mod = lm(Y~X+0)

### variables used for resampling
Y_m = predict(mod)
res = mod$residuals
err = Y

### perform resampling of residuals 
b_residuals = sapply(1:m, FUN = function(i) {
   Y_s = Y_m + sample(res,n)

### perform resampling of errors 
b_errors = sapply(1:m, FUN = function(i) {
   Y_s = Y_m + sample(err,n)

### perform paired resampling 
b_paired = sapply(1:m, FUN = function(i) {
   selection = rep(1:(n/2),2)+rbinom(n,1,0.5)*n/2
   Y_s = Y_m + res[selection]

### plot histograms 
hist(b_residuals, breaks = seq(-0.5,0.5,0.02), 
     freq = 0, ylim = c(0,10), main = "resampling of residuals", xlab = expression(beta[1]))
lines(mod$coefficients[1]*c(1,1),c(0,10), lty = 2, col = 2)
hist(b_paired, breaks = seq(-0.5,0.5,0.02),
     freq = 0, ylim = c(0,10), main = "resampling of residual pairs", xlab = expression(beta[1]))
lines(mod$coefficients[1]*c(1,1),c(0,10), lty = 2, col = 2)
hist(b_errors, breaks = seq(-0.5,0.5,0.02),
     freq = 0, ylim = c(0,10), main = "resampling of true errors", xlab = expression(beta[1]))

lines(mod$coefficients[1]*c(1,1),c(0,10), lty = 2, col = 2)

  • $\begingroup$ In the article, the Figure 1 also shows that sampling corrected residuals already gives the correct coverage rates. I can not reproduce their results with the paired resampling. Resampling the pairs should underestimate the variance. Whereas they get "the average variance of $\hat\beta^\star_1$ roughly overestimates the true variance of $\hat\beta_1$". It is unclear to me what they did in order to get this overestimated variance. In the supplementary material they speak about using the boot package in R, but not how they did it exactly. $\endgroup$ Commented May 14, 2023 at 8:50
  • $\begingroup$ It is also surprising that their jackknife method did not work. If you can not get correct results with OLS, where a correct confidence interval can be easily computed analytically, then something must be wrong with the implementation of the method. $\endgroup$ Commented May 14, 2023 at 9:09
  • $\begingroup$ Might it be possible that you are looking at just a single coordinate while they are looking at the full coefficient beta? $\endgroup$ Commented May 14, 2023 at 9:13
  • $\begingroup$ @LandonCarter The sample distributions and confidence intervals are computed for each single $\beta_i$ seperately. For presentation of the example here, I picked out a single one, $\beta_1$, because the sample distributions of the $\beta_i$ have different mean and variance. Without some rescaling and shifting you can't treat them together. In another not posted code that I wrote, computing the CI error rates, I looked at all beta components. In the end it doesn't matter because the effect for $\beta_1$ is the same as for the other $\beta_i$. I believe that the article also just uses $\beta_1$. $\endgroup$ Commented May 14, 2023 at 9:39
  • 1
    $\begingroup$ I see now why their re-sampling of pairs is causing an overestimation of the variance. The reason is because resampling the pairs will result in several pairs being excluded (and in the place of it some pairs have their weight being increased). Effectively this is reducing the sample size and with a reduced sample size you get a larger variance. $\endgroup$ Commented May 15, 2023 at 7:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.