Leverage as a Free Lunch: Four Properties You Already Paid For

What Makes It Free?

When you fit a linear regression, you compute a matrix decomposition of your design matrix (the table of features). This decomposition — the SVD — is the engine that produces the model's coefficients. But hidden inside that same decomposition, at no extra cost, is a number for each data point called its leverage score. This number tells you how unusual that point's features are relative to the rest of the data.

The leverage score is like a receipt stapled to the back of your model. You already paid for it when you fit the model. You just need to flip the receipt over and read it.

In concrete terms, this means the following. Suppose you have already fit a linear regression on a dataset with 1,000 rows and 30 features. To get adaptive prediction intervals using leverage, you need to store one $30 \times 30$ matrix (from the SVD you already computed) and perform one small matrix-vector product per test point. That is it — a few microseconds of additional work per prediction. No retraining, no new datasets, no validation loops.

Every other method for making prediction intervals adaptive requires additional work beyond fitting the base model: training extra models, selecting hyperparameters, or running expensive per-point computations. Leverage asks for nothing extra. This distinction matters in practice: in production systems where prediction latency is constrained, or in scientific settings where computational budgets are tight, the difference between "zero additional cost" and "train two more models" is often the difference between adopting adaptive intervals or not.

The Four Properties

Leverage scores have four properties that, taken together, distinguish them from every other adaptation mechanism:

1. Closed-form. There is a fixed formula for the weight function. You plug in the leverage score, and you get the answer. Nothing to tune, nothing to search over, nothing to cross-validate. The formula is: divide by the square root of one plus the leverage.
2. Model-free. Leverage is computed from the data (the feature matrix), not from any particular model. It measures how unusual a point's features are, regardless of whether you are using linear regression, a random forest, or a neural network. The geometric content — "this point is in a sparse region of feature space" — does not depend on your prediction method.
3. Computationally negligible. One matrix decomposition, which you probably already computed when fitting the model. After that, each new point's leverage costs a handful of multiplications — far less than making a prediction with a random forest, let alone training a new model.
4. Immune to the sign flip. In Part 6, we saw that training residuals are compressed at high-leverage points while prediction errors are amplified. Any method that estimates uncertainty from training residuals gets the direction backwards at high-leverage points. Leverage does not look at residuals at all. It is computed from the feature matrix. There is no sign to flip.

Among these methods, only leverage-based weighting has all four properties. Each alternative lacks at least two.

How It Works, Without Any Math

Here is the complete pipeline, from raw data to adaptive prediction intervals, using leverage:

The leverage pipeline: the SVD is already computed, so the extra cost is nearly zero.

Compare this to what the alternatives require:

• CQR needs two extra quantile regression models. Each requires choosing a model class, setting hyperparameters, and spending significant compute on training. In practice, this at least doubles the engineering effort and computational cost of the uncertainty estimation step.
• Studentized CP needs one extra scale-estimation model, and it inherits the sign flip problem: it learns the wrong variance at high-leverage points. The method is most inaccurate precisely where accurate uncertainty estimates matter most — at unusual feature combinations far from the training distribution.
• Localized CP is expensive at prediction time: it must compare each test point to every calibration point, with a kernel function and bandwidth that need to be chosen. For large test sets, this cost grows rapidly and can become a bottleneck.

When Does Leverage Actually Help?

Leverage-based adaptation helps when different points have meaningfully different leverage scores. If every point in the data has roughly the same leverage, the weights are roughly constant, and there is nothing to adapt. The method still works — it just reduces to the vanilla constant-width intervals. This is a graceful degradation, not a failure: you lose nothing by trying.

There is a simple diagnostic you can compute before committing to any method. Compute the leverage dispersion ratio: the standard deviation of the leverage scores divided by their mean. If this ratio is greater than 1, there is substantial variation, and leverage-based adaptation will noticeably improve interval quality. If the ratio is between 0.5 and 1, there is moderate variation, and the improvement will be real but not dramatic. If the ratio is below 0.5, leverage variation is mild, and the intervals will be similar to constant-width ones. Importantly, this diagnostic itself is free — it requires only the leverage scores you have already computed.

The leverage dispersion diagnostic: a quick check for whether leverage-based adaptation will help.

When does leverage variation tend to be large? Two situations stand out. First, when the number of features $p$ is a non-negligible fraction of the number of training points (say, $p/n > 0.05$). In this regime, individual data points have enough influence on the fit that their geometric position matters. Second, when the features have a "spiked" structure — a few dominant directions with most features being near-redundant. This is common in datasets with correlated features, where a handful of principal components explain most of the variance and a point that is unusual along those components will have high leverage. In modern datasets with many features and moderate sample sizes, one or both conditions are often present.

What the SVD Already Gives Us

After fitting OLS on training data $\mathcal{D}_1$, we compute the thin SVD: $\mathbf{X} = \mathbf{U\Sigma V}^\top$. From this single decomposition, we have:

1. The model: $\hat{f}(x) = x^\top \hat{\beta}$, where $\hat{\beta} = \mathbf{V\Sigma}^{-1}\mathbf{U}^\top \mathbf{Y}$.
2. Leverage scores: $h(x) = \|\mathbf{\Sigma}^{-1}\mathbf{V}^\top x\|^2$ for any point $x$, computable in $O(p)$ time.
3. The prediction variance formula: $\text{Var}(Y_{\text{new}} - \hat{Y}(x)) = \sigma^2(1 + h(x))$.

Everything that follows is a consequence of these three facts. Nothing additional needs to be estimated or learned.

Property 1: Closed-Form Weight

The variance formula $\sigma^2(1+h(x))$ means that prediction errors are not identically distributed. Larger leverage means larger variance. To make them comparable, we divide by the scale factor:

This is the leverage-based weight function. It is a known formula of a known quantity. There is no estimation step, no tuning parameter, no model selection. The formula is derived from the exact variance decomposition (Part 6), not from an approximation.

Compare: CQR requires choosing two quantile regression models, each with its own architecture and hyperparameters. Studentized CP requires choosing a scale estimation model. Both involve cross-validation or other tuning procedures. The leverage weight involves none of this.

Property 2: Depends Only on the Design Matrix

The hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$ is a property of $\mathbf{X}$ alone. It does not depend on the response $\mathbf{Y}$, the model $\hat{f}$, or the noise distribution. The leverage score $h(x) = x^\top(\mathbf{X}^\top\mathbf{X})^{-1}x$ is similarly model-free.

This has a powerful consequence: you can use leverage-based weights even when the prediction model is not OLS. If you fit a random forest, a neural network, or any other model, you can still compute leverage scores from the raw feature matrix $\mathbf{X}$ and use them to weight the conformal scores. The coverage guarantee of weighted conformal prediction holds for any weight function and any model (Part 2). The leverage weights will correctly identify where the feature-space geometry makes prediction harder, regardless of the prediction method.

Property 3: Computational Cost

The total cost of leverage-based adaptation breaks down as follows:

• SVD of $\mathbf{X}$: $O(n_1 p^2)$. This is typically already computed as part of fitting OLS.
• Leverage for $n_2$ calibration points: $O(n_2 \cdot p)$. One matrix-vector product per point.
• Leverage for $n_{\text{test}}$ test points: $O(n_{\text{test}} \cdot p)$. Same operation.

Total additional cost beyond model fitting: $O((n_2 + n_{\text{test}}) \cdot p)$.

Computational cost comparison for $n_1 = 1000$, $n_2 = 500$, $p = 30$, $n_{\text{test}} = 1000$. Leverage is one to two orders of magnitude cheaper.

Property 4: No Residuals, No Sign Flip

The sign flip (Part 6) is the observation that training residuals have variance $\sigma^2(1 - h_i)$ while prediction errors have variance $\sigma^2(1 + h(x))$. Any method that estimates the scale function from training residuals — including Studentized CP — learns a function that decreases at high leverage, when the truth increases.

Leverage scores bypass this entirely. They are computed from $\mathbf{X}$, not from residuals. The formula $\sigma^2(1+h)$ describes prediction-time variance directly. At a high-leverage point with $h = 0.5$:

• Leverage-based method: width $\propto \sqrt{1 + 0.5} = \sqrt{1.5} \approx 1.22$.
• Residual-based method: learns width $\propto \sqrt{1 - 0.5} = \sqrt{0.5} \approx 0.71$.

The residual-based method underestimates the uncertainty by a factor of $1.22 / 0.71 \approx 1.73$. The leverage-based method gets it exactly right.

The Leverage Dispersion Diagnostic

Define the leverage dispersion ratio:

computed over the calibration set. Since $\text{mean}(h_i) = p/n_1$ for training leverages (and approximately $p/n_1$ for new points), this ratio measures how spread out the leverages are relative to their average.

• $\hat{\eta} > 1$: substantial leverage variation. Some points have leverage several times the average. Leverage weighting will produce noticeably better conditional coverage than vanilla CP.
• $0.5 < \hat{\eta} < 1$: moderate variation. Leverage weighting improves over vanilla, but the gap may be modest.
• $\hat{\eta} < 0.5$: mild variation. Leverages are concentrated near their mean, and leverage-weighted intervals will be similar to constant-width intervals.

Beyond Homoscedastic Noise

The weight $w(h) = (1+h)^{-1/2}$ is derived under homoscedastic noise: $\text{Var}(\varepsilon \mid X) = \sigma^2$. What if the noise is itself heteroscedastic, with variance depending on leverage?

If $\text{Var}(\varepsilon \mid X) = \sigma^2 g(h(X))$ for some function $g$, then the total prediction error variance becomes:

The first term is the noise variance; the second is the estimation variance amplified by leverage. The optimal weight is:

Special cases:

• Homoscedastic ($g = 1$): $w^*(h) = (1 + h)^{-1/2}$. This is the standard leverage-based weight.
• Leverage-proportional noise ($g(h) = 1 + h$): $w^*(h) = (1 + 2h)^{-1/2}$. Both noise and estimation variance grow with leverage; the weight compensates for both.

The key observation is that $w^*$ is still a function of leverage alone. No auxiliary model is needed — just the right functional form applied to a quantity we already have.

Full Comparison

Formal Statement of the Four Properties

Property 1 (Closed-Form). Let $Y = X^\top\beta + \varepsilon$ with $\varepsilon \sim (0, \sigma^2)$ and $\hat{\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{Y}$ the OLS estimator. For a new point $x$ independent of the training data:

This is exact (not asymptotic) under the stated assumptions. The variance-stabilizing weight $w(h) = (1+h)^{-1/2}$ follows immediately: after dividing by $\sqrt{1+h(x)}$, the normalized residuals have constant variance $\sigma^2$.

Proof sketch. Write $Y_{\text{new}} - x^\top\hat{\beta} = \varepsilon_{\text{new}} - x^\top(\hat{\beta} - \beta)$. The two terms are independent (new noise vs. estimation error). $\text{Var}(\varepsilon_{\text{new}}) = \sigma^2$. $\text{Var}(x^\top(\hat{\beta} - \beta)) = \sigma^2 x^\top(\mathbf{X}^\top\mathbf{X})^{-1}x = \sigma^2 h(x)$. Summing gives $\sigma^2(1 + h(x))$.

Property 2 (Model-Free). The leverage $h(x) = x^\top(\mathbf{X}^\top\mathbf{X})^{-1}x$ depends only on the design matrix $\mathbf{X}$ and the test point $x$. It is invariant to:

• The response vector $\mathbf{Y}$.
• The fitted model $\hat{f}$ (OLS, ridge, random forest, etc.).
• The noise distribution beyond its variance.

Property 3 (Computational Cost). The cost accounting is as follows:

Total additional cost: $O(n_1 p^2 + n_2 p + n_{\text{test}} \cdot p)$, dominated by the SVD if not already computed, or $O((n_2 + n_{\text{test}}) \cdot p)$ if it is.

For comparison, CQR with gradient-boosted trees (100 trees, depth 6) costs approximately $O(n_1 \cdot 100 \cdot p \cdot 6) = O(600 \cdot n_1 \cdot p)$ per quantile model, for a total of $O(1200 \cdot n_1 \cdot p)$. With $n_1 = 1000$ and $p = 30$, this is $\sim 3.6 \times 10^7$ operations, versus $\sim 4.5 \times 10^4$ for leverage computation of 1500 points. The ratio is approximately $800 \times$.

Property 4 (Sign-Flip Immunity). For a training point $i$, the residual is $e_i = Y_i - \hat{Y}_i = (I - H)_i \cdot \varepsilon$, so $\text{Var}(e_i) = \sigma^2(1 - h_i)$. For a new test point $x$, $\text{Var}(Y_{\text{new}} - \hat{Y}(x)) = \sigma^2(1 + h(x))$. The sign of the leverage contribution flips between training and prediction.

Any method that regresses $|e_i|$ on $h_i$ to estimate $\sigma(x)$ learns the function $\sigma\sqrt{1-h}$ instead of $\sigma\sqrt{1+h}$. At the point $h = p/n_1$ (the average leverage), the relative error is:

For $p/n_1 = 0.1$, this is $\approx 10\%$. For $p/n_1 = 0.3$, this is $\approx 35\%$. The error grows with $p/n_1$ and is worst at the highest-leverage points, where accurate uncertainty estimation matters most.

Leverage-based weights use $h(x)$ directly from the design matrix, never passing through the training residuals. The formula $\sigma^2(1+h)$ is the prediction-time variance, not an estimate of it via training residuals.

Randomized SVD

For very large datasets where even the standard SVD is expensive, randomized SVD (Halko, Martinsson, and Tropp, 2011) computes an approximate thin SVD in $O(n_1 p \log p)$ time — nearly linear in the data size. The approximation error is controlled and can be made arbitrarily small by oversampling a few extra dimensions. This means that leverage computation scales to datasets with millions of rows and thousands of features.

The randomized SVD produces $\tilde{\mathbf{U}}, \tilde{\mathbf{\Sigma}}, \tilde{\mathbf{V}}$ such that $\|\mathbf{X} - \tilde{\mathbf{U}}\tilde{\mathbf{\Sigma}}\tilde{\mathbf{V}}^\top\|$ is close to optimal. The approximate leverage scores $\tilde{h}(x) = \|\tilde{\mathbf{\Sigma}}^{-1}\tilde{\mathbf{V}}^\top x\|^2$ satisfy $|\tilde{h}(x) - h(x)| \leq \epsilon$ for a user-chosen tolerance $\epsilon$, provided a modest oversampling parameter (typically 5-10 extra dimensions).

General Heteroscedastic Case

Suppose $\text{Var}(\varepsilon \mid X) = \sigma^2 g(h(X))$, where $g: [0,1] \to \mathbb{R}_+$ is a known or estimable function. Then:

The optimal variance-stabilizing weight is:

Important special cases and their optimal weights:

In every case, the optimal weight depends only on $h$. The entire adaptation mechanism is a univariate function of leverage.

When Heteroscedasticity Is Not Captured by Leverage

If $\text{Var}(\varepsilon \mid X) = \sigma^2(X) $ depends on features in ways that are not a function of $h(X)$, leverage-based adaptation handles only the geometric (estimation-variance) component. The total prediction error variance is then:

where $\sigma_0^2$ is the average noise level and $\sigma^2(x)$ is the feature-dependent noise. The leverage component $\sigma_0^2 \cdot h(x)$ is still handled exactly by leverage weighting. For the non-leverage component $\sigma^2(x)$, one can combine leverage weighting with a lightweight scale estimator. A natural hybrid approach is to use a small random forest (say, 10 trees) to estimate $\hat{\sigma}(x)$ from training residuals, then define the weight as $w(x) = (\hat{\sigma}^2(x) \cdot (1 + h(x)))^{-1/2}$. The leverage component handles the geometry; the scale estimator handles the feature-dependent noise.

Even in this hybrid approach, the leverage component provides a substantial portion of the adaptation for free. The scale estimator needs to capture only the residual heteroscedasticity not explained by leverage, which is a simpler estimation problem than capturing all heteroscedasticity from scratch.

Full Mathematical Pipeline