Three Faces of Leverage Scores
11 min read
Statisticians, numerical linear algebraists, and spectral graph theorists each attach a number to every row of a matrix. It took me embarrassingly long to notice it is the same number.
A leverage score measures how much an observation influences its own fitted value in linear regression [HW, TAS ’78]. The beautiful connection is that the same number guides both matrix approximation and graph simplification.
Self-Influence in Regression
We can measure an observation’s influence by varying its response and tracking the fitted value. Fix a design matrix \(\mathbf{A}\in\mathbb{R}^{n\times d}\) with full column rank and rows \(\mathbf{a}_i^\top\), where \(\mathbf{a}_i\in\mathbb{R}^d\). The responses are \(\mathbf{y}\in\mathbb{R}^n\). Least squares chooses coefficients \(\boldsymbol{\beta}\in\mathbb{R}^d\) to minimize the squared error \(\|\mathbf{y}-\mathbf{A}\boldsymbol{\beta}\|_2^2\). Setting the gradient to zero gives the normal equations:
\[ \mathbf{A}^\top\mathbf{A}\,\boldsymbol{\beta} = \mathbf{A}^\top\mathbf{y}. \]
Substitute the solution back into the fitted values:
\[\begin{align} \hat{\mathbf{y}} &= \mathbf{A}\boldsymbol{\beta} \\ &= \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{y}. \end{align}\]
The matrix multiplying \(\mathbf{y}\) is the hat matrix \(\mathbf{H}=\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\in\mathbb{R}^{n\times n}\), named for the operation that takes \(\mathbf{y}\) to \(\hat{\mathbf{y}}\). The \(i\)-th fitted value is a linear combination of the responses, and its coefficient on \(y_i\) measures self-influence. This coefficient is the row’s leverage score:
\[\begin{align} \tau_i &=\frac{\partial\hat y_i}{\partial y_i} \\ &=[\mathbf{H}]_{i,i} \\ &=\mathbf{a}_i^\top(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{a}_i. \end{align}\]
The second line differentiates the linear fit; the third reads the diagonal entry of the hat matrix. The score depends only on \(\mathbf{A}\). Moving \(y_i\) changes the fit, but it does not change how strongly that observation influences its own fitted value.
With \(d = 1\) the matrix is a single column of numbers \(a_1, \dots, a_n\), and the formula reduces to row \(i\)’s share of the total squared mass:
\[ \tau_i = \frac{a_i^2}{\sum_{j=1}^{n} a_j^2} . \]
Ten identical rows split the mass ten ways, a share of \(1/10\) each. A row holding all the mass takes share \(1\). Every share sits in \([0, 1]\), and the shares sum to \(1\) no matter how many rows there are.
In the plot, the teal point steers the line from the far right: delete it and the fit swings from the solid line to the dashed one, though its residual is small. The amber point does the opposite, carrying the largest residual on the plot while barely moving the fit. Leverage is about being alone in a direction of \(\mathbf{A}\), and only the residual looks at the responses at all.
The Range of Leverage
The range of possible influence follows from the geometry of the fit: the hat matrix \(\mathbf{H}\) is a projection. Multiply it by itself and the middle factors \((\mathbf{A}^\top\mathbf{A})^{-1}(\mathbf{A}^\top\mathbf{A})\) collapse to the identity:
\[\begin{align} \mathbf{H}^2 &= \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top \\ &= \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top \\ &= \mathbf{H}. \end{align}\]
It is also symmetric, \(\mathbf{H}^\top = \mathbf{H}\), so it is the orthogonal projection onto the column space of \(\mathbf{A}\). Idempotence and symmetry together give \(\mathbf{H} = \mathbf{H}^\top\mathbf{H}\), which lets us read a diagonal entry as a squared row norm and bound it below by its own square:
\[\begin{align} \tau_i = [\mathbf{H}]_{i,i} &= [\mathbf{H}^\top\mathbf{H}]_{i,i} \\ &= \sum_{j} [\mathbf{H}]_{i,j}^2 \\ &\ge [\mathbf{H}]_{i,i}^2 = \tau_i^2. \end{align}\]
A number no smaller than its own square lies in \([0, 1]\), so \(0 \le \tau_i \le 1\).
Take the extreme \(\tau_i = 1\). The identity \(\mathbf{H} = \mathbf{H}^\top\mathbf{H}\) turns the diagonal entry into a squared length, with \(\mathbf{e}_i \in \mathbb{R}^n\) the \(i\)-th standard basis vector:
\[ \tau_i = \mathbf{e}_i^\top \mathbf{H}\mathbf{e}_i = \|\mathbf{H}\mathbf{e}_i\|_2^2 . \]
A projection preserves length only by preserving the vector, so \(\tau_i = 1\) makes \(\mathbf{e}_i\) a fixed point:
\[ \mathbf{H}\mathbf{e}_i = \mathbf{e}_i . \]
Push that through the fitted value, moving \(\mathbf{H}\) onto \(\mathbf{e}_i\) by symmetry:
\[\begin{align} \hat{y}_i &= \mathbf{e}_i^\top \mathbf{H}\mathbf{y} \\ &= (\mathbf{H}\mathbf{e}_i)^\top \mathbf{y} \\ &= \mathbf{e}_i^\top \mathbf{y} \\ &= y_i . \end{align}\]
The fit reproduces observation \(i\) exactly, for every \(\mathbf{y}\), however absurd it is.
Rebuilding a Row
The regression calculation assigned a score using only \(\mathbf{A}\). We can get the same score by trying to reconstruct the row, without choosing any responses.
Claim: the leverage score is the least energy needed to build row \(i\) out of all the rows:
\[ \tau_i = \min\big\{ \|\mathbf{w}\|_2^2 \;:\; \mathbf{w} \in \mathbb{R}^n,\ \mathbf{A}^\top\mathbf{w} = \mathbf{a}_i \big\}. \]
Write out the constraint one row at a time:
\[ \mathbf{A}^\top\mathbf{w} = \sum_{j=1}^{n} w_j\, \mathbf{a}_j = \mathbf{a}_i . \]
A feasible \(\mathbf{w}\) writes row \(i\) as a linear combination of all the rows, and \(\|\mathbf{w}\|_2^2 = \sum_j w_j^2\) is the energy of that combination. One feasible choice is \(\mathbf{w} = \mathbf{e}_i\): row \(i\) builds itself at energy \(1\), so \(\tau_i \le 1\) falls out again. The minimum is small when the other rows can reconstruct \(\mathbf{a}_i\) with coefficients of low squared norm. Ten parallel copies of one row split the coefficient ten ways for total energy \(10 \cdot (1/10)^2 = 1/10\), the picture of a redundant row. A row holding the only mass in some direction of \(\mathbb{R}^d\) cannot be reconstructed by the others and has \(\tau_i = 1\).
Proof
Consider the candidate \(\mathbf{w}^* = \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{a}_i\), the \(i\)-th column of the hat matrix \(\mathbf{H}\mathbf{e}_i\). It is feasible, since the leading \(\mathbf{A}^\top\mathbf{A}\) cancels the inverse:
\[ \mathbf{A}^\top\mathbf{w}^* = (\mathbf{A}^\top\mathbf{A})(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{a}_i = \mathbf{a}_i . \]
Cancelling the same middle factors makes its energy equal the leverage score:
\[\begin{align} \|\mathbf{w}^*\|_2^2 &= \mathbf{a}_i^\top(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{a}_i \\ &= \mathbf{a}_i^\top(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{a}_i \\ &= \tau_i . \end{align}\]
It is also optimal. For any feasible \(\mathbf{w}\), the difference satisfies \(\mathbf{A}^\top(\mathbf{w} - \mathbf{w}^*) = \mathbf{0}\) and is therefore orthogonal to the column space of \(\mathbf{A}\), while \(\mathbf{w}^*\) lies inside it. Pythagoras then splits the energy:
\[ \|\mathbf{w}\|_2^2 = \|\mathbf{w}^*\|_2^2 + \|\mathbf{w} - \mathbf{w}^*\|_2^2 \ge \|\mathbf{w}^*\|_2^2 = \tau_i . \]
The reconstruction constraint also has an electrical interpretation: for a graph incidence matrix, the coefficients describe a flow and their squared norm is its dissipated power.
Effective Resistance
Take a connected graph with \(d\) vertices and edge set \(E\) of size \(n\), orient each edge arbitrarily, and let the row of \(\mathbf{B} \in \mathbb{R}^{n \times d}\) for edge \((u, v)\) be \(\mathbf{e}_u - \mathbf{e}_v\), the difference of the endpoints’ indicator vectors in \(\mathbb{R}^d\). Multiplying out, the Gram matrix is the graph Laplacian \(\mathbf{L} \in \mathbb{R}^{d \times d}\):
\[ \mathbf{B}^\top\mathbf{B} = \mathbf{L}, \]
with degrees on the diagonal and \(-1\) in entry \((u, v)\) for each edge. Every row of \(\mathbf{B}\) sums to zero. Since the graph is connected, the all-ones vector spans the kernel and the rank is \(d-1\). The leverage formula for a rank-deficient matrix uses a pseudoinverse, so \(\mathbf{L}^+\) takes the place of the inverse Gram matrix.
Make every edge a one-ohm resistor, push one ampere in at \(u\), and pull it out at \(v\). Collect the vertex potentials in \(\mathbf{p} \in \mathbb{R}^d\). Ohm’s law makes an edge’s current equal the potential drop across it, so the net current leaving vertex \(z\) sums those drops over its neighbors \(z'\):
\[ \sum_{z'} (p_z - p_{z'}) = [\mathbf{L}\mathbf{p}]_z . \]
Conservation makes that net current \(+1\) at \(u\), \(-1\) at \(v\), and zero elsewhere, a linear system in \(\mathbf{p}\):
\[ \mathbf{L}\mathbf{p} = \mathbf{e}_u - \mathbf{e}_v, \qquad \text{solved by} \quad \mathbf{p} = \mathbf{L}^+(\mathbf{e}_u - \mathbf{e}_v) . \]
The effective resistance between \(u\) and \(v\) is the voltage this experiment requires, the potential gap \(p_u - p_v\):
\[ p_u - p_v = (\mathbf{e}_u - \mathbf{e}_v)^\top \mathbf{L}^+ (\mathbf{e}_u - \mathbf{e}_v) . \]
That is the definition of \(\tau_i\) with \(\mathbf{e}_u - \mathbf{e}_v\) as the row: the effective resistance of an edge is its leverage score in \(\mathbf{B}\).
A coefficient vector \(\mathbf{w}\) with one entry per edge is a flow, and the constraint \(\mathbf{B}^\top\mathbf{w} = \mathbf{e}_u - \mathbf{e}_v\) says one unit of flow enters at \(u\), exits at \(v\), and is conserved at every other vertex. The energy \(\sum_e w_e^2\) is the power the flow dissipates in the resistors, so the minimum-energy representation becomes Thomson’s principle: the electrical flow is the cheapest unit flow, and its energy is the effective resistance.
An edge has leverage one exactly when it is a bridge. If deleting edge \((u,v)\) disconnects the graph, every unit flow pushes the whole unit across that edge, so its energy is at least \(1\); an edge on a cycle sends a fraction around the detour and comes in under \(1\).
In the diagram, thickness tracks resistance: the clique’s edges have many alternate routes, while the bridge has none and carries the whole unit.
Sampling the Rows
The minimum-energy interpretation suggests sampling rows that are hard to replace more often [DMM, SODA ’06]. To check that choice, we need to see what a sampled row contributes to the approximation. Write the Gram matrix in the outer-product view, one rank-one term per row:
\[ \mathbf{A}^\top\mathbf{A} = \sum_{i=1}^{n} \mathbf{a}_i \mathbf{a}_i^\top . \]
Sampling rows is Monte Carlo on this sum. Discard any zero rows, which contribute nothing, so every remaining row has positive leverage. Draw \(m\) indices independently, taking index \(i\) with probability \(p_i\), and stack the sampled rows scaled by \(1/\sqrt{m\,p_i}\) into \(\tilde{\mathbf{A}} \in \mathbb{R}^{m \times d}\). With every \(p_i>0\), inverse-probability scaling gives an unbiased Gram matrix. Summing the expected contributions of the \(m\) draws gives:
\[ \mathbb{E}[\tilde{\mathbf{A}}^\top\tilde{\mathbf{A}}] =m\sum_{i=1}^n p_i\frac{\mathbf{a}_i\mathbf{a}_i^\top}{m p_i} =\mathbf{A}^\top\mathbf{A}. \]
Unbiasedness does not ensure concentration, and uniform probabilities \(p_i = 1/n\) can produce a singular sampled Gram matrix. If one direction of the column space lives in a single row (picture \(n - 1\) generic rows supported on the first \(d - 1\) coordinates, plus one needle \(\mathbf{e}_d^\top\)), uniform sampling needs \(m = \Omega(n)\) draws to see the needle once. Until it does, the sampled Gram matrix is singular in that direction: a relative error of one, when the target is a small relative error \(\varepsilon\). The isolated row has leverage \(1\): no combination of the other rows can reconstruct it.
Whitening shows why leverage is the right sampling probability for every matrix. Full column rank makes \(\mathbf{A}^\top\mathbf{A}\) positive definite, with an orthonormal eigenbasis \(\mathbf{v}_1, \dots, \mathbf{v}_d \in \mathbb{R}^d\) and positive eigenvalues \(\lambda_1, \dots, \lambda_d\), so halving the exponents defines its inverse square root:
\[ \mathbf{A}^\top\mathbf{A} = \sum_{k=1}^{d} \lambda_k\, \mathbf{v}_k \mathbf{v}_k^\top, \qquad (\mathbf{A}^\top\mathbf{A})^{-1/2} = \sum_{k=1}^{d} \lambda_k^{-1/2}\, \mathbf{v}_k \mathbf{v}_k^\top . \]
Conjugating by it replaces each row by its whitened version:
\[ \tilde{\mathbf{a}}_i = (\mathbf{A}^\top\mathbf{A})^{-1/2}\,\mathbf{a}_i . \]
Whitening turns the target into \(\mathbf{I}_d\), so “preserve the spectrum” becomes “keep every eigenvalue of the sampled sum near \(1\)”. The whitened term \(\tilde{\mathbf{a}}_i\tilde{\mathbf{a}}_i^\top\) has spectral norm equal to the squared whitened length, exactly the leverage score:
\[ \|\tilde{\mathbf{a}}_i\|_2^2 = \mathbf{a}_i^\top(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{a}_i = \tau_i . \]
To turn leverage scores into probabilities, we need their sum. For our full-column-rank matrix \(\mathbf{A}\), that sum is \(d\). Cycling the factors inside the trace cancels the inverse Gram matrix:
\[ \sum_{i=1}^n\tau_i =\mathrm{tr}\big((\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{A}\big) =d. \]
We can therefore sample a row in proportion to its leverage:
\[ p_i = \frac{\tau_i}{d} . \]
Whichever row is drawn, the scaled term has spectral norm exactly \(d\):
\[ \Big\|\,\frac{\tilde{\mathbf{a}}_i\tilde{\mathbf{a}}_i^\top}{p_i}\,\Big\| = \frac{\tau_i}{p_i} = d , \]
so no sampled term can dominate the sum. A matrix Chernoff bound (Tropp [Tropp, FoCM ’12]) turns that flatness into a two-sided guarantee from \(m = O(d \log d / \varepsilon^2)\) draws, with high probability:
\[ (1 - \varepsilon)\, \mathbf{A}^\top\mathbf{A} \;\preceq\; \tilde{\mathbf{A}}^\top\tilde{\mathbf{A}} \;\preceq\; (1 + \varepsilon)\, \mathbf{A}^\top\mathbf{A}, \]
where \(\preceq\) compares quadratic forms direction by direction. The sketch preserves the squared length \(\|\mathbf{A}\boldsymbol{\beta}\|_2^2\) for every coefficient vector \(\boldsymbol{\beta}\), with a sample size that does not depend on \(n\). To preserve all least-squares residuals \(\|\mathbf{A}\boldsymbol{\beta}-\mathbf{y}\|_2^2\), apply the same argument to the augmented matrix \([\mathbf{A}\;\mathbf{y}]\), using its column space if it is rank-deficient. (The \(\log d\) is not slack: on the identity matrix every row has leverage \(1\), so leverage sampling is uniform and collecting all \(d\) rows is coupon collecting.)
Computing \(\tau_i\) exactly needs \((\mathbf{A}^\top\mathbf{A})^{-1}\), about as expensive as the problem the sketch was meant to shrink. Constant-factor approximations are enough, at a constant factor in \(m\), and they need no inverse.
Applications
On the column space of a graph incidence matrix, the same argument applies with rank \(d-1\). The trace identity gives the sum of the effective resistances of the edges of any connected graph:
\[ \sum_{(u,v) \in E} \tau_{(u,v)} = d - 1 . \]
This is Foster’s theorem (1949), and it gives the normalizing constant for sampling edges in proportion to their resistance. Spielman & Srivastava [SS, STOC ’08] use it to sparsify graphs: sample \(O(d \log d / \varepsilon^2)\) edges of a \(d\)-vertex graph with probability proportional to effective resistance, and the Laplacian, hence every cut and eigenvalue, is preserved to \(1 \pm \varepsilon\). On the barbell graph, uniform sampling can miss the high-resistance bridge and disconnect the sparsifier, while resistance sampling keeps it with high probability.
Regression packages print leverage beside the residual, while randomized least-squares solvers use it to sample a smaller problem. Leverage-score sampling of coalitions is also the engine inside LeverageSHAP, which turns the KernelSHAP regression into a Shapley value estimator with provable guarantees (Musco & Witter [MW, ICLR ’25]).