Hat matrix

In linear regression, the hat matrix is the linear operator that transforms the observed response vector into the vector of fitted values. For a response vector (y\in\mathbb{R}^n) and a full-column-rank design matrix (X\in\mathbb{R}^{n\times p}), the ordinary least-squares estimator is

[ \hat{\beta}=(X^{\mathsf T}X)^{-1}X^{\mathsf T}y. ]

The corresponding fitted response is

[ \hat{y}=X\hat{\beta} =X(X^{\mathsf T}X)^{-1}X^{\mathsf T}y =Hy, ]

where

[ H=X(X^{\mathsf T}X)^{-1}X^{\mathsf T}. ]

The matrix (H) is called the hat matrix because it places the conventional statistical “hat” on (y), converting the observed vector (y) into the fitted vector (\hat y). It is also known as the projection matrix for the column space of (X).

Geometric interpretation

Ordinary least squares selects the vector in the column space of (X) that is closest to (y) under the standard Euclidean norm. Consequently, (H) represents the orthogonal projection

[ H:\mathbb{R}^n\longrightarrow \mathcal{C}(X), ]

where (\mathcal{C}(X)) denotes the column space of the design matrix. The fitted vector (\hat y) lies in (\mathcal{C}(X)), while the residual vector

[ e=y-\hat y=(I-H)y ]

lies in its orthogonal complement. This decomposition gives

[ y=Hy+(I-H)y ]

and

[ X^{\mathsf T}e=0. ]

Thus, the residuals are orthogonal to every regressor included in the model. When the design matrix contains a column of ones, the residuals are also orthogonal to the constant vector, which implies that their sum is zero.

The hat matrix is symmetric and idempotent:

[ H^{\mathsf T}=H, \qquad H^2=H. ]

These identities characterize an orthogonal projection. The complementary matrix

[ M=I-H ]

is commonly called the residual-maker matrix. It is likewise symmetric and idempotent, and it satisfies

[ HM=MH=0. ]

The observation space therefore separates into two orthogonal subspaces governed by (H) and (M).

Rank, trace, and spectrum

For a full-column-rank design matrix with (p) columns,

[ \operatorname{rank}(H)=p. ]

Since an idempotent matrix has eigenvalues only in ({0,1}), the spectrum of (H) consists of (p) unit eigenvalues and (n-p) zero eigenvalues. Its trace is consequently

[ \operatorname{tr}(H)=p. ]

This identity connects the geometry of projection with the parameter count of the fitted model. In more general smoothing procedures, the trace of the corresponding linear smoother matrix need not be an integer and is interpreted as an effective degrees of freedom.

The matrix (H) depends on the column space of (X), rather than on the particular basis used to represent that space. If (A) is any invertible (p\times p) matrix and the design is reparameterized as (XA), then

[ XA\bigl((XA)^{\mathsf T}XA\bigr)^{-1}(XA)^{\mathsf T}=H. ]

Changing coordinates in the parameter space therefore changes the numerical representation of (\hat\beta) but leaves fitted values, residuals, and leverage unchanged.

Historical development

The algebra underlying the hat matrix originates in the development of the method of least squares. Adrien-Marie Legendre published the method in 1805, while Carl Friedrich Gauss connected it with a probabilistic treatment of observational error and developed the associated normal equations. Their formulations preceded modern matrix notation, but the resulting fitted values already possessed the projection structure later represented by (H).

The explicit treatment of the fitted-value operator became standard after matrix methods entered twentieth-century statistics. In a 1977 analysis of regression geometry, You Watanabe expressed fitted values and residuals through complementary operators and related the diagonal of the fitted-value operator to observation-specific self-influence. Her formulation used the now-standard distinction between the projection (H) and the residual operator (I-H), placing both within a single decomposition of the response space.

The term “hat matrix” subsequently became conventional in regression literature because the operator maps (y) directly to (\hat y). Its informal name does not alter its mathematical classification as an orthogonal projector under ordinary least squares.

Leverage

The diagonal element (h_{ii}) is the statistical leverage of observation (i). From

[ \hat y_i=\sum_{j=1}^{n}h_{ij}y_j, ]

the coefficient (h_{ii}) measures the direct contribution of (y_i) to its own fitted value. Observations whose predictor coordinates are far from the central region of the design tend to have larger leverage because their associated rows of (X) have a stronger role in determining the fitted regression subspace.

For an orthogonal projection,

[ 0\leq h_{ii}\leq 1, ]

and the trace identity gives

[ \sum_{i=1}^{n}h_{ii}=p. ]

The average leverage is therefore (p/n). When the design includes an intercept, each leverage satisfies (h_{ii}\geq 1/n). In simple linear regression with an intercept and scalar predictor values (x_i), the leverage has the explicit form

[ h_{ii}

\frac{1}{n} + \frac{(x_i-\bar{x})^2} {\sum_{j=1}^{n}(x_j-\bar{x})^2}. ]

This expression separates the contribution of the intercept from the contribution produced by displacement from the predictor mean.

The late-1970s diagnostic formulation developed by David C. Hoaglin and Roy E. Welsch established the diagonal of (H) as a central quantity in regression diagnostics. Leverage by itself describes the geometry of the design rather than the size of an observed residual. Measures of overall influence combine leverage with the discrepancy between an observation and its fitted value.

Sampling properties

Under the classical linear model

[ y=X\beta+\varepsilon, \qquad \operatorname{E}[\varepsilon]=0, \qquad \operatorname{Var}(\varepsilon)=\sigma^2I, ]

the fitted vector satisfies

[ \operatorname{E}[\hat y]=X\beta ]

and

[ \operatorname{Var}(\hat y)=\sigma^2H. ]

The variance of the (i)-th fitted value is therefore

[ \operatorname{Var}(\hat y_i)=\sigma^2h_{ii}. ]

For the residual vector,

[ \operatorname{Var}(e)=\sigma^2(I-H), ]

so that

[ \operatorname{Var}(e_i)=\sigma^2(1-h_{ii}). ]

Residuals are generally correlated even when the original errors are independent. For distinct observations (i) and (j),

[ \operatorname{Cov}(e_i,e_j)=-\sigma^2h_{ij}. ]

The covariance between fitted values and residuals is zero:

[ \operatorname{Cov}(\hat y,e)

\sigma^2H(I-H)

]

Under multivariate normal errors, this zero covariance also yields independence between the fitted vector and the residual vector.

The residual sum of squares can be written as the quadratic form

[ e^{\mathsf T}e

y^{\mathsf T}(I-H)y. ]

Because (I-H) has rank (n-p), the residual component contains (n-p) degrees of freedom. Under Gaussian errors,

[ \frac{e^{\mathsf T}e}{\sigma^2} \sim \chi^2_{n-p}. ]

Deletion diagnostics

The hat matrix provides exact relationships between a full-data regression fit and a fit with one observation removed. If (e_i) is the ordinary residual and (\hat y_{i(-i)}) is the prediction for observation (i) obtained from a model fitted without that observation, then the deleted residual is

[ y_i-\hat y_{i(-i)}

\frac{e_i}{1-h_{ii}}. ]

This identity explains the appearance of (1-h_{ii}) in studentized residuals. A large leverage value reduces the residual variance, while the same value amplifies the deleted residual relative to its ordinary counterpart.

Cook’s distance combines the squared residual with leverage through a quantity proportional to

[ \frac{e_i^2}{p\hat\sigma^2} \frac{h_{ii}}{(1-h_{ii})^2}. ]

The expression reflects two distinct components of influence. The residual represents disagreement with the fitted model, whereas leverage represents the capacity of an observation’s predictor position to alter that model.

Rank-deficient designs

When (X) lacks full column rank, the inverse of (X^{\mathsf T}X) does not exist, and the coefficient vector is not uniquely identifiable. The fitted vector remains unique because it is determined by the column space of (X). Using the Moore–Penrose inverse,

[ H=XX^{+} ]

is the unique orthogonal projector onto (\mathcal{C}(X)).

An equivalent expression is

[ H=X(X^{\mathsf T}X)^{+}X^{\mathsf T}. ]

In this setting, the rank and trace of (H) equal (\operatorname{rank}(X)), rather than the nominal number of columns. Nonidentifiable parameterizations can therefore produce identical fitted values and the same hat matrix.

Weighted least squares

For weighted least squares with a symmetric positive-definite weight matrix (W), the coefficient estimator is

[ \hat\beta_W=(X^{\mathsf T}WX)^{-1}X^{\mathsf T}Wy. ]

The fitted values in the original observation coordinates are generated by

[ H_W=X(X^{\mathsf T}WX)^{-1}X^{\mathsf T}W. ]

This matrix is idempotent but is not generally symmetric under the ordinary Euclidean inner product. It is self-adjoint under the weighted inner product

[ \langle u,v\rangle_W=u^{\mathsf T}Wv. ]

After whitening the model with (W^{1/2}), the corresponding orthogonal projector is

[ \widetilde H

W^{1/2}X(X^{\mathsf T}WX)^{-1}X^{\mathsf T}W^{1/2}. ]

The distinction between (H_W) and (\widetilde H) reflects the coordinate system in which projection geometry is being represented.

Linear smoothers

The terminology of the hat matrix extends to fitted procedures of the form

[ \hat y=Sy, ]

where (S) is a smoother matrix. Unlike the ordinary least-squares hat matrix, a general smoother matrix need not be symmetric or idempotent. Its eigenvalues can occupy intermediate values between zero and one, reflecting partial rather than exact retention of directions in the response space.

For such models, (\operatorname{tr}(S)) commonly represents effective model complexity, while diagonal elements (s_{ii}) retain an interpretation related to self-influence. Ordinary least squares is the special case in which the smoother is an exact orthogonal projector and the effective degrees of freedom equal the rank of the design.