Linear regression models a continuous target as a weighted sum of features plus noise: y = β0 + β1x1 + … + βpxp + ε. Fit by ordinary least squares (OLS), it has a closed-form solution, a clean geometric interpretation, and a mature statistical theory for confidence intervals and hypothesis tests. It is the starting point for every regression problem, the foundation of polynomial and regularized regression, and — as a single linear layer — a building block of every neural network.
Part of Regression within Supervised Learning.
In matrix form, with a design matrix X (n × (p+1), including a column of ones for the intercept), the model is y = Xβ + ε. OLS picks the coefficients that minimize the residual sum of squares:
RSS(β) = Σ (yi − xi·β)2 = ‖y − Xβ‖2
Minimizing squared error is equivalent to maximum likelihood under Gaussian noise, and it estimates the conditional mean E[y | x]. If you care about the median or are worried about outliers, absolute-error (QuantileRegressor) or Huber loss (HuberRegressor) are the robust alternatives. Standard evaluation metrics are RMSE, MAE and R2 = 1 − RSS/TSS, the fraction of variance explained relative to predicting the mean.
Setting the gradient −2XT(y − Xβ) to zero gives the normal equations XTXβ = XTy. Geometrically, the fitted vector ŷ = Xβ̂ is the orthogonal projection of y onto the space spanned by the feature columns, and the residual is perpendicular to every column. When XTX is invertible the solution is unique.
Numerically, libraries solve with a QR decomposition or SVD (numpy.linalg.lstsq, scikit-learn's LinearRegression), which is stable even when features are nearly collinear; the cost is about O(n p2). For very large or streaming datasets, stochastic gradient descent (SGDRegressor) updates β ← β + η(yi − xi·β)xi one sample or mini-batch at a time; it requires scaled features and a learning-rate schedule.
Predictions only require that the relationship is approximately linear in the chosen features. Valid standard errors, p-values and intervals require more (the Gauss–Markov conditions plus normality for small samples):
y, weighted least squares, or robust (HC) standard errors.βj is the expected change in y for a one-unit increase in xj holding all other features fixed. That "holding fixed" clause is where interpretation breaks: with correlated features, a coefficient can be small or even flip sign while the feature is strongly related to y on its own. Coefficients describe association, not causation, unless the design supports a causal reading. Standardizing features makes magnitudes comparable; a log-transformed target turns coefficients into approximate percentage effects (100·βj% per unit for small β).
import numpy as np
import statsmodels.api as sm
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
X, y = fetch_california_housing(return_X_y=True, as_frame=True)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=0)
# Prediction-focused: scikit-learn
lr = LinearRegression().fit(X_tr, y_tr)
pred = lr.predict(X_te)
print("RMSE:", round(mean_squared_error(y_te, pred) ** 0.5, 3), " R^2:", round(r2_score(y_te, pred), 3))
# Inference-focused: statsmodels gives standard errors, t-stats, CIs
ols = sm.OLS(y_tr, sm.add_constant(X_tr)).fit(cov_type="HC3") # heteroscedasticity-robust SEs
print(ols.summary().tables[1])
# Closed form via least squares (what LinearRegression does internally)
A = np.c_[np.ones(len(X_tr)), X_tr.to_numpy()]
beta, *_ = np.linalg.lstsq(A, y_tr.to_numpy(), rcond=None)
print(np.allclose(beta[1:], lr.coef_))
Use scikit-learn when the goal is prediction inside a pipeline, and statsmodels when you need the inferential table. HC3 standard errors remain valid under heteroscedasticity, which this housing dataset clearly has.
| Strength | Limitation |
|---|---|
| Closed-form, fast, deterministic. | Only captures relationships that are linear in the features. |
| Coefficients and confidence intervals are interpretable. | Sensitive to outliers (squared loss). |
| Extrapolates in a predictable (linear) way. | Unstable with collinear or many features — regularize. |
| Strong baseline; well understood statistical theory. | Usually less accurate than tree ensembles on complex tabular data. |
X singular; drop one level.p > n) — OLS has infinitely many solutions; use ridge or lasso.With one feature the normal equations reduce to two familiar formulas: slope β1 = Sxy / Sxx and intercept β0 = ȳ − β1x̄, where Sxy = Σ(xi − x̄)(yi − ȳ) and Sxx = Σ(xi − x̄)2.
For the four points above, x̄ = 2.5 and ȳ = 4. The centered products are 3, 0.5, 0.5 and 3, so Sxy = 7; the centered squares sum to Sxx = 5. The slope is 1.4 and the intercept 0.5. Fitted values are 1.9, 3.3, 4.7 and 6.1, residuals 0.1, −0.3, 0.3, −0.1. Two properties of OLS are visible: the residuals sum to zero (because there is an intercept) and they are uncorrelated with x.
The residual sum of squares is 0.2 against a total sum of squares of 10, so R2 = 0.98. The noise variance estimate is RSS/(n − 2) = 0.1, and the slope's standard error is √(0.1 / Sxx) ≈ 0.14, giving a t-statistic near 10. With only two residual degrees of freedom, though, that inference rests heavily on the normality assumption; small samples are exactly where the assumptions in Section 3 stop being a formality.
"Linear" means linear in the coefficients, not in the raw inputs, so much of the craft is building the right columns of X:
x1x2 lets the slope of one feature depend on another; always keep the main effects alongside it.SplineTransformer) capture curvature with less wild extrapolation than high-degree polynomials.Each added column costs variance. Once the column count grows, switch to ridge or lasso and choose the penalty by cross-validation.
Two intervals are easy to confuse. A confidence interval for the mean response at x0 reflects uncertainty in the estimated line only and shrinks toward zero as n grows. A prediction interval for a new observation adds the irreducible noise σ2, so it never shrinks below roughly ±2σ. Forecasting a single house price needs the prediction interval; estimating the average price of houses with given features needs the confidence interval. In statsmodels, results.get_prediction(X_new).summary_frame(alpha=0.05) returns both (mean_ci_* and obs_ci_*). Both widen as x0 moves away from the bulk of the training data, which is the formal version of "do not extrapolate".
In production, monitor the residual distribution over time, not just RMSE: a drifting mean residual signals that the relationship has shifted, and the model should be refit on recent data.