Last reviewed
May 26, 2026
Sources
No citations yet
Review status
Needs citations
Revision
v3 · 4,370 words
Add missing citations, update stale details, or suggest a clearer explanation.
Regression is a family of statistical and machine learning methods for modelling the relationship between a numeric outcome variable and one or more explanatory variables (often called predictors, features, or covariates). The goal is to estimate a function that maps inputs to a continuous output and to quantify how the inputs influence that output, both for prediction on unseen data and for inference about the underlying relationship. The term traces to Francis Galton, who in 1886 described how the heights of adult children tend to "regress" toward the population mean of their parents' generation, a phenomenon now known as regression toward the mean.[1] The mathematical machinery, however, predates Galton: Adrien-Marie Legendre published the method of least squares in 1805 and Carl Friedrich Gauss gave a probabilistic derivation in 1809, both motivated by the calculation of planetary and cometary orbits.[2][3] Regression today spans simple straight-line fits, generalized linear models, penalized estimators such as the lasso and ridge, and large nonparametric and tree-based predictors used throughout the empirical sciences and machine learning.
The earliest published method for fitting a line through noisy observations is Legendre's "Sur la méthode des moindres carrés", an appendix to Nouvelles méthodes pour la détermination des orbites des comètes (Paris, 1805). Legendre proposed choosing parameter values that minimize the sum of squared residuals, calling the technique a "very simple" device that "establishes a kind of equilibrium among the errors".[2] Gauss, in Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium (1809), supplied a probabilistic justification by showing that the least squares estimate is the maximum likelihood solution when errors are normally distributed; he also claimed prior use of the method since 1795, triggering a long priority dispute with Legendre.[3] A later result, now called the Gauss-Markov theorem after Gauss's 1821 Theoria Combinationis Observationum and Andrei Markov's 1900 textbook reformulation, states that among linear unbiased estimators of the regression coefficients, ordinary least squares has the smallest variance whenever the errors have zero mean, constant variance, and no correlation.[4]
Galton supplied the name and the conceptual frame. In "Regression towards Mediocrity in Hereditary Stature", published in the Journal of the Anthropological Institute of Great Britain and Ireland in 1886, he reported measurements of 928 adult children and their parents and observed that tall parents produced children who, on average, were tall but closer to the mean stature than their parents; short parents produced offspring closer to the mean in the opposite direction.[1] Galton called this drift "regression towards mediocrity" and quantified it through what he termed the coefficient of regression, the slope relating mid-parent stature to filial stature. His student Karl Pearson, in collaboration with George Udny Yule, generalized the procedure into the multiple linear regression framework familiar today; Sewall Wright introduced the coefficient of determination R-squared in 1921 as part of his work on path analysis in genetics.[5]
Twentieth-century developments extended the framework along several axes. David R. Cox's 1958 paper "The Regression Analysis of Binary Sequences" introduced the modern formulation of logistic regression as a regression model for binary outcomes.[6] John Nelder and Robert Wedderburn unified normal, binomial, Poisson, and gamma regression as instances of a single family in their 1972 paper "Generalized Linear Models", along with a common fitting algorithm (iteratively reweighted least squares) and a single deviance-based goodness of fit measure.[7] Andrey Tikhonov proposed L2 regularization for ill-posed integral equations in 1943; Arthur Hoerl and Robert Kennard reintroduced the same penalty into multiple regression as ridge regression in their 1970 Technometrics paper "Ridge Regression: Biased Estimation for Nonorthogonal Problems".[8][9] Robert Tibshirani's 1996 paper introduced the lasso (least absolute shrinkage and selection operator), which uses an L1 penalty to perform simultaneous shrinkage and variable selection;[10] Hui Zou and Trevor Hastie's 2005 elastic net combined L1 and L2 penalties to handle correlated predictors.[11] Nonparametric estimators including the kernel regression of Èlizbar Nadaraya and Geoffrey Watson (1964) and the regression trees of Leo Breiman, Jerome Friedman, Richard Olshen, and Charles Stone (1984) opened the way to flexible function fitting without a fixed parametric form.[12][13]
In supervised learning the predictive task is conventionally split into regression, in which the response is a numeric quantity that can in principle take any value on a continuum, and classification, in which the response is a categorical label drawn from a finite set.[14] Regressing housing price on square footage, estimating blood pressure from age and weight, or forecasting next quarter's electricity demand are regression problems; predicting whether a borrower will default, whether an email is spam, or which of ten digits is shown in an image are classification problems. The two tasks share most of their algorithmic machinery (decision trees, neural networks, kernel methods, ensemble learners) but differ in the loss function (squared or absolute error versus cross-entropy or zero-one loss), in evaluation metrics (RMSE versus accuracy or area under the ROC curve), and in the form of the output (a real number versus a probability or label).[14]
The line between the two is porous. logistic regression is technically a regression model in that it estimates a continuous quantity (the log-odds of the positive class) and reports a probability, but it is almost always used as a classifier by thresholding that probability. Conversely, ordinal and count outcomes, where the response is discrete but the categories carry numerical meaning, are often handled by regression models such as ordered probit or Poisson regression. Ranking, survival, and structured prediction problems borrow tools from both camps. The label "regression" in modern machine learning usually means simply that the response is continuous and that performance is measured by an error magnitude rather than a classification rate.
The simple linear regression model assumes that a scalar response y depends on a scalar predictor x through y = beta_0 + beta_1 x + epsilon, where epsilon is a random error term with mean zero. Multiple linear regression generalizes this to y = beta_0 + beta_1 x_1 + ... + beta_p x_p + epsilon and admits a compact matrix form y = X beta + epsilon, with X the n-by-(p+1) design matrix that includes a column of ones for the intercept.[15] The ordinary least squares estimator chooses beta to minimize the residual sum of squares (y - X beta)'(y - X beta); when X has full column rank, the closed-form solution is beta hat = (X'X)^(-1) X' y. Equivalently, beta hat is the orthogonal projection of y onto the column space of X.[15]
OLS rests on a standard set of assumptions. Linearity says that the conditional expectation of the response is a linear function of the parameters (not necessarily of the original variables, since polynomial and interaction terms are still linear in their coefficients). Strict exogeneity, or zero conditional mean of the errors, requires E[epsilon | X] = 0; this rules out simultaneity and omitted-variable confounding. Independence rules out autocorrelation across observations; homoscedasticity requires constant error variance; and no perfect multicollinearity requires that X have full column rank.[4] Under these assumptions the Gauss-Markov theorem guarantees that OLS is the best linear unbiased estimator. Adding the further assumption that the errors are jointly normal makes OLS coincide with the maximum likelihood estimator and supplies exact t and F distributions for hypothesis tests on the coefficients.[3][15]
Polynomial regression, in which polynomial powers and cross terms of the original predictors are added as columns of X, remains a special case of multiple linear regression: it is linear in the unknown coefficients even though it can describe curved relationships between y and the raw inputs. Stepwise procedures, best-subset selection, and information criteria such as Akaike's AIC (Akaike, 1974) and Schwarz's BIC (1978) supply tools for choosing which terms to include in the design matrix.[16][17] Weighted least squares and generalized least squares relax the homoscedasticity and uncorrelated-errors assumptions by replacing OLS's identity weighting with an explicit covariance matrix of the errors.[15]
Real outcomes are often not continuous and Gaussian. Counts of accidents, binary disease indicators, durations, and proportions all violate the normal-error assumption. Nelder and Wedderburn's 1972 framework, the generalized linear model or generalized linear model, reorganizes a large class of regression problems around three components: a response distribution from the exponential family, a linear predictor eta = X beta, and a link function g that connects the mean of the response to the linear predictor through g(mu) = eta.[7] Ordinary linear regression is the case in which the response is Gaussian and the link is the identity; other choices yield familiar named models.
logistic regression, used when the response is binary, takes the response to be Bernoulli and the link to be the logit function logit(p) = log(p / (1 - p)). The model fits the log-odds as a linear function of the predictors and the resulting coefficients can be interpreted as log odds ratios. Cox's 1958 paper introduced the procedure in essentially this form, and it remains the workhorse of medical statistics, credit scoring, and many other applications.[6] Poisson regression handles counts with a log link, log(mu) = X beta, modelling expected event counts as a multiplicative function of the predictors; it underlies much epidemiology and insurance pricing.[7] Probit regression replaces the logit link with the inverse standard normal cumulative distribution function and was historically dominant in bioassay before logistic models overtook it. Gamma and inverse Gaussian GLMs handle positive continuous outcomes such as costs and durations.
All GLMs are typically fit by maximum likelihood (mle), and Nelder and Wedderburn showed that the score equations can be solved by an iteratively reweighted least squares routine that reuses the same code path as OLS at each iteration.[7] Goodness of fit is measured by the residual deviance, a likelihood-ratio statistic against the saturated model, while individual coefficients are tested with Wald or likelihood-ratio statistics.
When the number of predictors p is large relative to the number of observations n, or when predictors are correlated, OLS estimates become unstable: small changes in the data produce large swings in the fitted coefficients, the matrix X'X becomes nearly singular, and out-of-sample prediction degrades. Regularized regression introduces a penalty on the coefficient vector that biases estimates toward zero in exchange for reduced variance, sitting deliberately on the favorable side of the bias variance tradeoff.
ridge regression adds a quadratic penalty lambda * sum(beta_j^2) to the residual sum of squares, giving the closed-form solution beta hat = (X'X + lambda I)^(-1) X' y. The estimator was introduced as a numerical device by Tikhonov for inverse problems in 1943 and reintroduced into statistics by Hoerl and Kennard in 1970 to stabilize multiple regression in the presence of multicollinearity.[8][9] Ridge regression shrinks all coefficients toward zero but rarely sets them exactly to zero, so it does not perform automatic variable selection. It is equivalent to maximum a posteriori estimation under a Gaussian prior on the coefficients and corresponds to weight decay in neural network training.[14]
The lasso, proposed by Tibshirani in 1996, replaces the quadratic penalty with an absolute-value penalty lambda * sum(|beta_j|).[10] Because the L1 ball has corners on the coordinate axes, the solution path of the lasso often hits these corners exactly, producing coefficient vectors with many entries equal to zero. The lasso therefore performs continuous variable selection along with shrinkage. It has no closed-form solution but can be computed efficiently with coordinate descent or the least angle regression algorithm of Bradley Efron, Trevor Hastie, Iain Johnstone, and Tibshirani (2004).[14] In the orthonormal design case the lasso reduces to soft-thresholding the OLS estimates.
The elastic net of Zou and Hastie (2005) combines the two penalties as alpha * lambda * sum(|beta_j|) + (1 - alpha) * lambda * sum(beta_j^2).[11] The L2 component preserves the grouping behavior of ridge regression (correlated predictors are kept together with similar coefficients) while the L1 component preserves the variable selection of the lasso. Elastic net is widely used in genomics, text classification, and other settings where predictors come in correlated blocks. Other penalty families in active use include the adaptive lasso, smoothly clipped absolute deviation (SCAD) of Jianqing Fan and Runze Li, the minimax concave penalty of Cun-Hui Zhang, the group lasso of Ming Yuan and Yi Lin, and the fused lasso of Tibshirani and collaborators for ordered features.[14]
A different penalty form, motivated by extreme distributions of the residuals rather than by collinearity, is the L1 regression or least absolute deviation estimator, which minimizes the sum of absolute residuals instead of the sum of squares. Quantile regression, developed by Roger Koenker and Gilbert Bassett Jr. in their 1978 Econometrica paper "Regression Quantiles", generalizes this idea by introducing a tilted absolute-value loss whose minimizer estimates a chosen conditional quantile (median, tenth percentile, ninetieth percentile) of the response.[18] Quantile regression makes no assumption that the conditional mean is the object of interest and is robust to heavy-tailed errors.
The methods so far assume a parametric form for the regression function. Nonparametric and machine learning approaches let the data determine the shape of the function.
Kernel regression, proposed independently by Nadaraya and Watson in 1964, estimates the conditional mean E[y | x] as a weighted average of nearby observed y values, with weights given by a kernel function centered at x and a bandwidth that controls the width of the smoothing window.[12] The Nadaraya-Watson estimator is the simplest member of a family that includes local linear and local polynomial regression. These methods are flexible but suffer from the curse of dimensionality: smoothing windows shrink rapidly as the dimension of x grows.
Spline regression fits piecewise polynomials joined at knot points with continuity constraints. Smoothing splines and penalized regression splines penalize integrated squared curvature, producing the function that minimizes residual sum of squares plus a wiggliness penalty. Generalized additive models (GAMs), introduced by Trevor Hastie and Tibshirani in 1986, extend this to multiple predictors by writing the regression function as a sum of smooth univariate functions of each predictor.[14]
Regression trees, formalized in the CART (Classification and Regression Trees) framework of Breiman, Friedman, Olshen, and Stone in their 1984 monograph, recursively partition the predictor space into rectangular regions and fit a constant in each region.[13] Single trees are easy to interpret but high variance; ensembles dramatically improve accuracy. random forest, introduced by Breiman in his 2001 paper of the same name, averages many decorrelated trees grown on bootstrap samples with random feature subsets at each split, and is now a default workhorse for tabular data.[19] gradient boosting, introduced in Friedman's 2001 paper "Greedy Function Approximation: A Gradient Boosting Machine", grows trees sequentially, with each tree fit to the negative gradient of a chosen loss function evaluated at the current predictions.[20] xgboost (Tianqi Chen and Carlos Guestrin, 2016), LightGBM (Microsoft, 2017), and CatBoost (Yandex, 2018) are widely used scalable implementations of gradient boosting that have repeatedly topped Kaggle leaderboards on tabular regression tasks.[21]
Bayesian linear regression places a prior distribution on the coefficients (commonly Gaussian) and derives the posterior given the data; under conjugate priors the posterior is again Gaussian and gives credible intervals for both coefficients and predictions. Gaussian process regression, treated systematically by Carl Edward Rasmussen and Christopher Williams in their 2006 book Gaussian Processes for Machine Learning, places a prior directly on the regression function: any finite collection of function values is jointly Gaussian with a covariance specified by a kernel function.[22] gaussian process regression produces a full posterior over functions and is widely used in geostatistics (where it is called kriging), Bayesian optimization, and emulation of expensive simulators.
neural network regression replaces the linear function X beta with a parameterized nonlinear function f(x; theta), typically a multilayer perceptron trained by stochastic gradient descent on squared-error loss. Deep networks dominate regression on unstructured data such as images, audio, and text, where they can exploit features that hand-engineered designs cannot easily encode. On tabular data the picture is mixed: gradient boosted decision trees often match or beat deep networks of comparable size, and the comparative advantage of deep learning over gradient boosting depends heavily on dataset size, feature type, and effort spent on tuning.[14] Hybrid models such as TabNet, FT-Transformer, and neural additive models attempt to combine the strengths of both.
Several optimization principles produce regression estimates from data. Ordinary least squares minimizes squared residuals and is the canonical choice for linear regression with Gaussian-like errors; under the Gauss-Markov assumptions it is the best linear unbiased estimator.[4] Maximum likelihood estimation (mle) maximizes the joint probability of the observed responses under a specified distributional model and is the standard fitting principle for GLMs and most parametric regression models; under the normal-error linear model it coincides with OLS.[3][7] Maximum a posteriori (MAP) estimation combines the likelihood with a prior on the parameters and is the Bayesian counterpart of regularized estimation: a Gaussian prior gives ridge, a Laplace prior gives the lasso, and a horseshoe prior gives an aggressive shrinkage estimator. Full Bayesian regression integrates over the posterior rather than maximizing it, producing distributions over coefficients and predictions.
For large datasets and nonlinear models, closed-form solutions are unavailable and the parameters are found by numerical optimization. gradient descent, stochastic gradient descent, and modern variants (Adam, AdamW, RMSprop) are universal default choices, especially for neural network regression and for fitting regularized linear models on big data. Coordinate descent is the dominant algorithm for fitting lasso and elastic net at scale; iteratively reweighted least squares is standard for GLMs. The least angle regression algorithm produces the full lasso solution path in a number of operations comparable to one OLS fit.[14]
Predictive performance of a regression model on continuous outcomes is most often summarized by mean squared error (MSE), root mean squared error (RMSE), or mean absolute error (MAE). MSE = (1/n) sum((y_i - y_i hat)^2) is the natural counterpart of the squared-error loss minimized during fitting and is sensitive to large outliers because of the square. RMSE = sqrt(MSE) is reported in the same units as the response and is intuitive to compare across studies. MAE = (1/n) sum(|y_i - y_i hat|) is less sensitive to outliers and corresponds to median rather than mean regression.[14] Mean absolute percentage error and symmetric MAPE are common in business forecasting; the log-cosh and Huber losses interpolate between squared and absolute error to combine smooth gradients with robustness.
The coefficient of determination R-squared, introduced by Sewall Wright in 1921, reports the fraction of total variance in y explained by the model: R-squared = 1 - SS_residual / SS_total.[5] For OLS on the same training data it lies between zero and one, but for arbitrary nonlinear models or out-of-sample evaluation it can be negative. Adjusted R-squared penalizes R-squared for the number of predictors, attempting to compensate for the mechanical increase that follows from adding any feature. Information criteria provide alternative measures of fit that penalize complexity directly: Akaike's AIC (1974) is -2 log L + 2k for log-likelihood L and parameter count k; Schwarz's BIC (1978) is -2 log L + k log n.[16][17] AIC is justified as an estimator of expected out-of-sample Kullback-Leibler divergence; BIC is justified as a large-sample approximation to a Bayes factor.
cross-validation, systematized by Mervyn Stone in his 1974 paper "Cross-Validatory Choice and Assessment of Statistical Predictions", estimates prediction error by repeatedly fitting the model on subsets of the data and evaluating on the held-out portion.[23] K-fold cross-validation partitions the data into k chunks and rotates which chunk is held out; leave-one-out cross-validation is the extreme case k = n; nested cross-validation separates the selection of hyperparameters from the estimation of final performance. Cross-validation is the standard tool for choosing regularization strengths in ridge, lasso, and elastic net and for tuning hyperparameters in tree ensembles and neural networks.[14]
Diagnostic plots are used to check that the assumptions underlying linear regression are met. Plots of residuals against fitted values should show no systematic pattern; curvature suggests missing nonlinear terms, fan-shaped scatter suggests heteroscedasticity. Quantile-quantile plots of residuals against the normal distribution check the normality assumption. Leverage and Cook's distance identify observations with disproportionate influence on the fit. Variance inflation factors quantify multicollinearity among predictors. Durbin-Watson and Ljung-Box statistics detect autocorrelation in time-ordered residuals.
Beyond prediction, regression supports inference: hypothesis tests and confidence intervals for the coefficients, predictions intervals for new observations, and decompositions of variance attributable to different predictors. Under the classical normal linear model the coefficient vector beta hat has a multivariate normal distribution with mean beta and covariance sigma squared (X'X)^(-1), and individual t statistics and overall F statistics have exact null distributions.[15] For non-Gaussian GLMs the same machinery applies asymptotically. Bootstrap resampling provides nonparametric confidence intervals when distributional assumptions are doubtful.
Selecting predictors based on the data (stepwise procedures, lasso, screening on p-values) and then performing standard inference on the selected variables produces invalid p-values and confidence intervals, an issue known as post-selection inference. Methods including the selective inference framework of Tibshirani, Jonathan Taylor, and collaborators, the knockoff filter of Rina Foygel Barber and Emmanuel Candes, and various sample-splitting schemes provide valid inference after model selection but are still less widely deployed than the older, technically invalid procedures.
Regression methods are pervasive across the empirical sciences and engineering. In economics and finance, the linear regression model underpins econometrics: estimation of demand curves and Phillips curves, returns of stocks on market factors (CAPM and Fama-French models), and treatment effects from observational data with instrumental variables and difference-in-differences identification strategies. Logistic regression is the workhorse of credit scoring, marketing response models, and customer churn prediction.
In epidemiology and medicine, logistic regression models case-control studies and binary outcomes; Cox proportional hazards regression (a separate method due to David R. Cox in 1972 building on his 1958 framework) models survival times; Poisson regression handles incidence counts; mixed-effects models handle clustered and longitudinal data. Regression underlies risk scoring systems such as the Framingham heart score, QRISK, and ASCVD calculators.
In machine learning, regression appears throughout supervised learning pipelines. Continuous targets such as image attributes, pose estimates, audio features, depth maps, regression heads on object detectors, value functions in reinforcement learning, and reward models in alignment work all use regression objectives. Many language model post-training procedures use regression-style fits of reward or value functions over candidate completions.
In time series analysis and forecasting, autoregressive models (AR, ARMA, ARIMA, VAR) treat the past values of a series as predictors for its future, augmented in modern forecasting systems with exogenous variables, splines, neural networks, and gradient boosting. State-space models combine regression-style measurement equations with dynamic latent processes.
In physical sciences and engineering, regression fits calibration curves, response surfaces in design of experiments, and emulators for expensive simulations. Gaussian process regression has become standard in geostatistics for spatial interpolation under the name kriging and in computer experiments for surrogate modelling.