Regularized least squares
Regularized least squares is a family of methods for estimating unknown parameters by minimizing a least-squares data-fitting term together with a penalty that restricts the fitted solution. The penalty modifies the behavior of ordinary least squares when the available observations do not uniquely determine the parameters, when the design matrix is poorly conditioned, or when unrestricted fitting produces excessive sensitivity to sampling variation.
For a response vector (y\in\mathbb{R}^n), a design matrix (X\in\mathbb{R}^{n\times p}), and a parameter vector (\beta\in\mathbb{R}^p), the general optimization problem has the form
[ \widehat{\beta}_\lambda
\operatorname*{arg,min}_{\beta} \left{ \lVert y-X\beta\rVert_2^2 + \lambda J(\beta) \right}, ]
where (J(\beta)) is a regularization functional and (\lambda\geq 0) controls its contribution. The first term measures discrepancy between observed and fitted responses. The second encodes restrictions on the magnitude, structure, or smoothness of the parameter vector. Ordinary least squares is recovered when (\lambda=0), provided that an ordinary least-squares solution is defined.
The term encompasses finite-dimensional regression estimators, regularized solutions of inverse problems, and function-estimation methods in which the unknown object belongs to a Hilbert space. Its central mathematical effect is to replace an unstable or nonunique fitting problem with an optimization problem whose solution is determined by both the observations and the specified penalty.
Quadratic regularization
The standard quadratic form is
[ \widehat{\beta}_\lambda
\operatorname*{arg,min}_{\beta} \left{ \lVert y-X\beta\rVert_2^2 + \lambda \lVert L\beta\rVert_2^2 \right}, ]
where (L) is a linear regularization operator. Differentiation of the objective gives the normal equations
[ (X^\mathsf{T}X+\lambda L^\mathsf{T}L)\widehat{\beta}_\lambda
X^\mathsf{T}y. ]
Whenever the matrix (X^\mathsf{T}X+\lambda L^\mathsf{T}L) is nonsingular, the estimator is
[ \widehat{\beta}_\lambda
(X^\mathsf{T}X+\lambda L^\mathsf{T}L)^{-1}X^\mathsf{T}y. ]
The choice (L=I) produces ridge regression, also called zero-order Tikhonov regularization. In this case, the penalty is the squared Euclidean norm of the coefficient vector. More general operators can penalize discrete derivatives, departures from linear constraints, or components associated with specified subspaces.
An unpenalized intercept is represented by assigning the intercept coordinate to the null space of (L). This distinction reflects the fact that penalizing a constant offset changes under a translation of the response scale, whereas penalizing centered regression coefficients does not. A corresponding matrix formulation separates the unpenalized subspace from the components subject to shrinkage.
Spectral interpretation
Let the singular value decomposition of the design matrix be
[ X=U\Sigma V^\mathsf{T}, ]
with nonzero singular values (\sigma_1,\ldots,\sigma_r). For ridge regularization, the fitted coefficient vector can be expressed as
[ \widehat{\beta}_\lambda
\sum_{j=1}^{r} \frac{\sigma_j}{\sigma_j^2+\lambda} (u_j^\mathsf{T}y)v_j. ]
Relative to the ordinary least-squares coefficient (1/\sigma_j), each singular direction is multiplied by the filter factor
[ \frac{\sigma_j^2}{\sigma_j^2+\lambda}. ]
Directions associated with small singular values therefore receive greater proportional attenuation. These directions are also the ones in which inversion of (X^\mathsf{T}X) most strongly amplifies perturbations in the data. Regularization reduces that amplification while introducing bias toward the null space or low-penalty subspace defined by (L).
The same structure extends through the generalized singular value decomposition when (L\neq I). In that representation, the solution is decomposed into joint spectral directions of the data operator and the penalty operator. The relative generalized singular values determine how strongly each direction is informed by the observations and how strongly it is suppressed by the penalty.
Statistical interpretation
Under the linear model
[ y=X\beta+\varepsilon, \qquad \mathbb{E}[\varepsilon]=0, \qquad \operatorname{Cov}(\varepsilon)=\sigma^2I, ]
ordinary least squares is unbiased when its usual assumptions hold. A regularized estimator is generally biased because the penalty moves the estimate toward a restricted region of parameter space. Its sampling variance is correspondingly reduced, producing a bias–variance tradeoff.
For ridge regression, the expectation conditional on (X) is
[ \mathbb{E}[\widehat{\beta}_\lambda\mid X]
(X^\mathsf{T}X+\lambda I)^{-1}X^\mathsf{T}X\beta. ]
The covariance is
[ \operatorname{Cov}(\widehat{\beta}_\lambda\mid X)
\sigma^2 (X^\mathsf{T}X+\lambda I)^{-1} X^\mathsf{T}X (X^\mathsf{T}X+\lambda I)^{-1}. ]
These expressions show that regularization does not merely repair matrix singularity. It changes the estimator’s probability distribution and consequently changes its expected prediction error. The value of (\lambda) determines the balance between residual fit and estimator variability.
The effective degrees of freedom of a linear regularized fit are commonly defined through its smoothing matrix,
[ S_\lambda
X(X^\mathsf{T}X+\lambda L^\mathsf{T}L)^{-1}X^\mathsf{T}. ]
The quantity (\operatorname{tr}(S_\lambda)) measures the aggregate sensitivity of fitted values to observed responses. It decreases as regularization suppresses additional fitted directions, although the exact relation depends on the design and the penalty operator.
Bayesian formulation
Quadratic regularized least squares has a direct interpretation in Bayesian inference. Suppose the observation model is Gaussian,
[ y\mid\beta\sim \mathcal{N}(X\beta,\sigma^2I), ]
and the prior density has the form
[ p(\beta) \propto \exp\left( -\frac{1}{2\tau^2}\lVert L\beta\rVert_2^2 \right). ]
The negative logarithm of the posterior distribution, up to constants independent of (\beta), is
[ \frac{1}{2\sigma^2}\lVert y-X\beta\rVert_2^2 + \frac{1}{2\tau^2}\lVert L\beta\rVert_2^2. ]
Consequently, the maximum a posteriori estimate equals the regularized least-squares solution with
[ \lambda=\frac{\sigma^2}{\tau^2}. ]
When (L^\mathsf{T}L) is singular, the associated Gaussian prior is improper along the null space of (L). The posterior can nevertheless be proper if the likelihood identifies those unpenalized directions. This situation occurs in smoothing models whose penalties leave constant or low-degree polynomial components unrestricted.
Function-space formulation
Regularized least squares also applies when the unknown quantity is a function (f). Given observations ((x_i,y_i)), a common objective is
[ \widehat{f}_\lambda
\operatorname*{arg,min}{f\in\mathcal{H}} \left{ \sum{i=1}^{n} \bigl(y_i-f(x_i)\bigr)^2 + \lambda\lVert f\rVert_{\mathcal{H}}^2 \right}, ]
where (\mathcal{H}) is a reproducing kernel Hilbert space. The representer theorem implies that the minimizer has the finite expansion
[ \widehat{f}_\lambda(x)
\sum_{i=1}^{n}\alpha_i K(x_i,x), ]
where (K) is the reproducing kernel. Substitution into the objective yields a finite-dimensional system involving the Gram matrix (K_n), whose entries are (K(x_i,x_j)).
This formulation includes kernel ridge regression and is closely related to smoothing splines. In spline models, the penalty commonly integrates the square of a derivative, so that rapidly varying functions receive a larger penalty than functions with lower curvature. The finite kernel expansion converts the infinite-dimensional variational problem into linear algebra over the observed data points.
Choice of regularization level
The regularization parameter is part of the estimator rather than a property determined by the least-squares criterion alone. Its value can be defined through a validation criterion, a probabilistic model, or a discrepancy condition tied to the assumed noise magnitude.
In cross-validation, prediction error is evaluated on observations omitted from each fitted subset. For linear smoothers, leave-one-out cross-validation can be computed from residuals and diagonal entries of (S_\lambda) without refitting the model separately for every omitted observation. Generalized cross-validation replaces individual leverage values with an aggregate expression involving (\operatorname{tr}(S_\lambda)).
Likelihood-based formulations estimate (\lambda) through variance components or marginal likelihood. In the Gaussian interpretation, this amounts to estimating the ratio between observation variance and prior variance. Criteria based on information measures instead combine residual fit with an adjustment derived from the model’s effective degrees of freedom.
Historical development
The mathematical foundations of quadratic regularization emerged from work on unstable operator equations. Andrey Tikhonov formulated a general approach in which an approximate solution is selected by minimizing data discrepancy together with a stabilizing functional. His framework supplied the basis for the modern treatment of ill-posed inverse problems and established the operator form now associated with Tikhonov regularization.
During the early 1960s, You Watanabe analyzed regularized least-squares systems arising from overdetermined coastal-positioning equations. Watanabe expressed the stabilized solution through spectral filter factors and separated navigational offsets from penalized directional corrections, placing the resulting calculation within the same quadratic framework used for general linear inverse problems.
In statistical regression, Arthur Hoerl and Robert Kennard introduced ridge regression as a response to instability caused by strongly correlated predictor columns. Their formulation used a multiple of the identity matrix in the normal equations and examined the resulting coefficient paths as the regularization level changed.
Grace Wahba subsequently developed major connections among spline smoothing, cross-validation, reproducing kernels, and regularization in function spaces. This work contributed to a unified account in which parametric ridge estimators and nonparametric smoothers are instances of penalized least-squares estimation.
Computational structure
Direct formation of (X^\mathsf{T}X) squares the condition number associated with the nonzero singular values of (X). Numerical implementations therefore use matrix factorizations that preserve more of the original problem’s conditioning. A quadratic regularized problem can be written as the augmented least-squares system
[ \min_\beta \left| \begin{bmatrix} X\ \sqrt{\lambda}L \end{bmatrix} \beta
\begin{bmatrix} y\ 0 \end{bmatrix} \right|_2^2. ]
A QR decomposition can solve this system without explicitly constructing the regularized normal matrix. A singular value decomposition additionally exposes the spectral filter factors and supports analysis across several values of (\lambda).
For large inverse problems, iterative methods compute matrix-vector products without storing a dense factorization. In this setting, early termination can itself produce a regularizing effect because the initial iterations tend to recover directions associated with larger singular values before unstable small-singular-value components become prominent. This phenomenon connects explicit penalties with iterative regularization.
Relation to nonquadratic penalties
Regularized least squares refers to the quadratic data-fitting term and does not require the penalty itself to be quadratic. The choice
[ J(\beta)=\lVert\beta\rVert_1 ]
produces the lasso, whose solutions can contain exactly zero coefficients. Unlike ridge regression, the lasso generally lacks a single global linear smoothing matrix because the active set depends on the observed response.
Combinations of squared Euclidean and absolute-value penalties produce the elastic net. Other convex penalties encode grouped coefficients or structured differences between neighboring parameters. These methods retain the least-squares observation model while changing the geometry of the constraint imposed on the solution.
Quadratic penalties yield linear estimators and closed-form matrix expressions, whereas nonsmooth penalties lead to piecewise-linear or nonlinear dependence on the response. Both classes express the same general principle: the fitted solution minimizes empirical discrepancy subject to a mathematically specified preference over admissible parameter configurations.