# linear regression

*Entry of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org) (course edition), by Alexander Jung, Aalto University. CC BY 4.0. Updated 2026-11-08. [HTML](https://dictionaryofml.org/terms/linreg.html) · [typeset PDF](https://dictionaryofml.org/terms/linreg.pdf) · [LaTeX](https://dictionaryofml.org/terms/linreg.tex) · [Python demo](https://dictionaryofml.org/terms/linreg.py) · [notebook](https://dictionaryofml.org/terms/linreg.ipynb)*

Math is LaTeX with the book's macros expanded, so no notation package is needed to read it.

## Abstract

Linear regression is a regression method that learns a linear hypothesis map for predicting the numeric label of a data point from its features. Its least squares variant chooses a map that minimizes the average squared error loss on a training set. Linear regression is an instance of empirical risk minimization that is obtained by using the linear model and the squared error loss. The optimal model parameters are determined by the normal equations. This system can be solved in closed form or iteratively, for example, by gradient descent. Adding a penalty term to the average squared error loss yields regularized variants such as ridge regression and least absolute shrinkage and selection operator.

## Definition

Linear regression can be used to predict tomorrow's temperature from weather measurements recorded over the last five days, using them as features and tomorrow's temperature as the label. Linear regression methods use a linear model: they learn a linear hypothesis map that delivers a prediction of the numeric label of a data point. The prediction is based solely on the features of the data point, which are fed as input to the learned hypothesis map.

Formally, linear regression learns a linear hypothesis map $h^{({\mathbf w})}({\mathbf x}) = {\mathbf w}^{\top}{\mathbf x}$ to predict the numeric label $y \in \mathbb{R}$ of a data point from its feature vector ${\mathbf x} = \big(x_{1},\,\ldots,\,x_{d}\big)^{\top} \in \mathbb{R}^{d}$. Here, $d$ denotes the number of features of a data point. The model parameters ${\mathbf w}$ are learned from a training set $\mathcal{D}^{(\rm train)} = \big\{ \left( {\mathbf x}^{(r)},y^{(r)} \right) \big\}_{r=1}^{m}$ of data points. Appending a constant feature to ${\mathbf x}$ lets one entry of ${\mathbf w}$ act as an intercept, so ${\mathbf w}^{\top}{\mathbf x}$ represents an affine function of the original measurements. Below, ${\mathbf x}$ denotes whichever feature vector is fed to the linear hypothesis map, original or augmented, unless stated otherwise.

The least squares variant of linear regression measures the quality of a linear hypothesis map by its average loss, the squared error loss, on the training set. As an instance of empirical risk minimization, it learns the model parameters ${\mathbf w}$ by solving the optimization problem

$$
\min_{{\mathbf w} \in \mathbb{R}^{d}} \frac{1}{m}
\sum_{r=1}^{m}
\big( y^{(r)} - {\mathbf w}^{\top} {\mathbf x}^{(r)} \big)^{2} \text{.}
$$

Fig. 1 illustrates this linear regression problem for a training set with a constant scalar feature $x^{(r)}=1$ for $r=1,\ldots,m$.

![Figure 1 of the entry linreg](https://dictionaryofml.org/terms/fig/linreg-fig1.svg)

*Figure 1: For a linear model with $d=1$ and the constant feature $x=1$ for every data point (blue filled circles), linear regression reduces to computing the average $\widehat{w} = (1/m) \sum_{r=1}^{m} y^{(r)}$ of the labels: this average minimizes the average squared error loss. The horizontal dashed line sits at the height $\widehat{w}$, and the vertical gray segments are the errors $y^{(r)} - \widehat{w}$ entering that loss*

The optimization problem can be written more compactly using the feature matrix ${\mathbf X}$ and the label vector ${\mathbf y}$,

$$
{\mathbf X} = \big({\mathbf x}^{(1)},\,\ldots,\,{\mathbf x}^{(m)}\big)^{\top}
\in \mathbb{R}^{m \times d} \text{,} \qquad
{\mathbf y} = \big( y^{(1)},\,\ldots,\,y^{(m)} \big)^{\top}
\in \mathbb{R}^{m} \text{.}
$$

In terms of ${\mathbf X}$ and ${\mathbf y}$, the optimization problem reads

$$
\min_{{\mathbf w} \in \mathbb{R}^{d}} \underbrace{\frac{1}{m}
\left\lVert {\mathbf y} - {\mathbf X} {\mathbf w} \right\rVert_{2}^{2}}_{f({\mathbf w})} \text{.}
$$

The average squared error loss $f({\mathbf w})$ is a convex and differentiable function of ${\mathbf w}$. By the zero-gradient condition (Boyd and Vandenberghe, 2004, Sect. 4.2.3), a vector $\widehat{{\mathbf w}}$ solves this optimization problem if and only if it satisfies the system of linear equations (Strang, 2016)

\begin{equation}
\tag{1}
{\mathbf X}^{\top}{\mathbf X} \widehat{{\mathbf w}} = {\mathbf X}^{\top} {\mathbf y} \text{.}
\end{equation}

The normal equations Eq. (1) always have a solution because ${\mathbf X}^{\top}{\mathbf y}$ lies in the column space of ${\mathbf X}^{\top}$, which coincides with the column space of ${\mathbf X}^{\top}{\mathbf X}$ (Meyer, 2000, Sect. 4.5).

However, the solution of the normal equations Eq. (1) is unique only if the feature matrix ${\mathbf X}$ has full column rank. This requires at least as many data points as features, $m \geq d$. When this holds, and the columns of ${\mathbf X}$ are linearly independent, the matrix ${\mathbf X}^{\top}{\mathbf X}$ is invertible and the normal equations Eq. (1) have the unique closed-form solution $\widehat{{\mathbf w}} = \big({\mathbf X}^{\top}{\mathbf X}\big)^{-1} {\mathbf X}^{\top} {\mathbf y}$ (Golub and Loan, 2013, Sect. 5.3.1).

If instead $m < d$ (more features than data points), ${\mathbf X}$ cannot have full column rank, ${\mathbf X}^{\top}{\mathbf X}$ is singular, and the optimization problem has infinitely many solutions (Golub and Loan, 2013, Sect. 5.5). All of them incur the same minimal average squared error loss on the training set, but their predictions for data points outside the training set can differ arbitrarily. Choosing among them is therefore a matter of generalization (see overfitting). A unique solution can be selected, for example, by picking the minimum-norm solution ${{\mathbf X}}^{+} {\mathbf y}$ given by the pseudoinverse ${{\mathbf X}}^{+}$ (Golub and Loan, 2013, Sect. 5.5.2).

A common remedy is regularization: adding a penalty term to the average squared error loss. The penalty term can be interpreted as an estimate of how much higher the loss is on data points outside the training set than on it (Hastie et al., 2009, Sect. 7.4). This yields regularized variants of linear regression. Ridge regression uses the penalty term $\alpha \left\lVert {\mathbf w} \right\rVert_{2}^{2}$ and the least absolute shrinkage and selection operator uses $\alpha \left\lVert {\mathbf w} \right\rVert_{1}$. The hyperparameter $\alpha > 0$ sets the weight of the penalty term.

Linear regression also has a statistical interpretation. Consider a probabilistic model with a joint probability distribution $\mathbb{P}^{({\mathbf x}, y)}$ over the features and the label. Under the squared error loss, the Bayes estimator (i.e., the hypothesis with minimum risk) is the conditional expectation $\mathbb{E}\{y \mid {\mathbf x}\}$ (Lehmann and Casella, 1998, Ch. 4; Papoulis and Pillai, 2002). In this statistical interpretation, ${\mathbf x}$ denotes the original feature vector. Let the features and the label be jointly Gaussian RVs with zero mean and invertible covariance matrix $\mathbf{C}^{({\mathbf x})} = \mathbb{E}\{{\mathbf x} {\mathbf x}^{\top}\}$. This Bayes estimator is then (Bishop, 2006, Sect. 2.3.1, Eq. (2.81))

$$
h^{*}({\mathbf x}) = \big({\mathbf w}^{\star}\big)^{\top}{\mathbf x}
\text{, with } {\mathbf w}^{\star} = \big(\mathbf{C}^{({\mathbf x})}\big)^{-1}\, \mathbf{c}^{({\mathbf x},y)} \text{.}
$$

Here, $\mathbf{c}^{({\mathbf x},y)} = \mathbb{E}\{{\mathbf x} y\}$ is the covariance between the features ${\mathbf x}$ and the label $y$. This Bayes estimator is linear in ${\mathbf x}$. Without the zero-mean assumption, it is an affine function of ${\mathbf x}$, which a constant feature absorbs as before.

Least-squares linear regression is the empirical risk minimization counterpart of this predictor. It replaces $\mathbf{C}^{({\mathbf x})}$ and $\mathbf{c}^{({\mathbf x},y)}$ with the sample-based estimates $(1/m){\mathbf X}^{\top}{\mathbf X}$ and $(1/m){\mathbf X}^{\top}{\mathbf y}$. For invertible ${\mathbf X}^{\top}{\mathbf X}$, this recovers the unique solution $\widehat{{\mathbf w}} = \big({\mathbf X}^{\top}{\mathbf X}\big)^{-1}{\mathbf X}^{\top}{\mathbf y}$ of the normal equations Eq. (1).

Instead of solving the normal equations Eq. (1) directly (via the inverse or pseudoinverse), machine learning methods often solve it by an iterative optimization method, as implemented in widely used machine learning software libraries such as scikit-learn (Pedregosa et al., 2011) and PyTorch (Paszke et al., 2019). Starting from initial model parameters ${\mathbf w}^{(0)}$, such a method repeatedly applies an update operator $\mathcal{F}$,

$$
{\mathbf w}^{(t+1)} = \mathcal{F}\big({\mathbf w}^{(t)}\big)
\text{, for } t = 0,1,\ldots \text{.}
$$

The operator $\mathcal{F}$ is designed such that the sequence of model parameters ${\mathbf w}^{(t)}$ converges to a solution of the normal equations Eq. (1).

One such iterative optimization method is gradient descent, which uses the gradient descent step operator

$$
\mathcal{F}^{(\eta)}({\mathbf w}) = {\mathbf w} - \eta \nabla f({\mathbf w}) \text{.}
$$

Here, $\eta > 0$ is a step size that must be chosen small enough to ensure convergence (Bertsekas, 2016, Prop. 1.2.2). The superscript in $\mathcal{F}^{(\eta)}$ makes the dependence of the operator on $\eta$ explicit. Inserting the average training error into the general form of the gradient descent step yields the explicit update

$$
\begin{aligned}
\mathcal{F}^{(\eta)}({\mathbf w})
&= {\mathbf w} + \frac{2\eta}{m} \sum_{r=1}^{m}
\big( y^{(r)} - {\mathbf w}^{\top}{\mathbf x}^{(r)} \big)
{\mathbf x}^{(r)} \\
&= \Big(\mathbf{I} - \frac{2\eta}{m}{\mathbf X}^{\top}{\mathbf X}\Big){\mathbf w}
+ \frac{2\eta}{m}{\mathbf X}^{\top}{\mathbf y} \text{.}
\end{aligned}
$$

The fixed points of $\mathcal{F}^{(\eta)}$ are exactly the solutions of the normal equations Eq. (1), since $\mathcal{F}^{(\eta)}(\widehat{{\mathbf w}}) = \widehat{{\mathbf w}}$ holds if and only if ${\mathbf X}^{\top}{\mathbf X}\widehat{{\mathbf w}} = {\mathbf X}^{\top}{\mathbf y}$. The update operator is affine,

$$
\mathcal{F}^{(\eta)}({\mathbf w}) = {\mathbf M} {\mathbf w} + {\mathbf b} \text{,}
$$

with linear part ${\mathbf M} = \mathbf{I} - (2\eta/m){\mathbf X}^{\top}{\mathbf X}$ and offset ${\mathbf b} = (2\eta/m)\, {\mathbf X}^{\top}{\mathbf y}$.

When ${\mathbf X}$ has full column rank, so that ${\mathbf X}^{\top}{\mathbf X}$ is positive definite, $\mathcal{F}^{(\eta)}$ is a contractive operator with respect to the Euclidean norm for every step size $0 < \eta < m/\lambda_{\max}$. Here, $\lambda_{\max}$ is the largest eigenvalue of ${\mathbf X}^{\top}{\mathbf X}$. The eigenvalues of the linear part ${\mathbf M}$ are $1 - (2\eta/m)\,\lambda_{j}$, with the positive eigenvalues $\lambda_{j}$ of ${\mathbf X}^{\top}{\mathbf X}$. They lie strictly between $-1$ and $1$ exactly for such step sizes. The more general case of any smooth, strongly convex objective function and step sizes $\eta \leq m/(\lambda_{\min} + \lambda_{\max})$ is covered by Nesterov (2004, Thm. 2.1.15), with the smallest eigenvalue $\lambda_{\min}$ of ${\mathbf X}^{\top}{\mathbf X}$. The fixed-point iteration converges to the unique solution $\widehat{{\mathbf w}}$. Moreover, the convergence speed of gradient descent is governed by the ratio of the largest to the smallest eigenvalue of ${\mathbf X}^{\top}{\mathbf X}$, that is, by its condition number $\kappa\!\left({\mathbf X}^{\top}{\mathbf X}\right) = \lambda_{\max}/\lambda_{\min}$ (Boyd and Vandenberghe, 2004, Sect. 9.3).

The gradient descent step above uses the entire training set to compute the gradient. When the data points instead arrive sequentially at time instants $t=1,2,\ldots$, or the training set is too large to fit in memory, the model parameters can be learned by online gradient descent. At each time step $t$, it applies a single gradient descent step. That step uses only the arriving data point $\left( {\mathbf x}^{(t)},y^{(t)} \right)$,

$$
{\mathbf w}^{(t+1)} = {\mathbf w}^{(t)} + 2\eta \big(
y^{(t)} - \big({\mathbf w}^{(t)}\big)^{\top}
{\mathbf x}^{(t)} \big) {\mathbf x}^{(t)} \text{.}
$$

This is the gradient descent step on the loss function $\big(y^{(t)} - {\mathbf w}^{\top}{\mathbf x}^{(t)}\big)^{2}$ of the arriving data point alone, rather than on the average squared error loss $f({\mathbf w})$. It is the least mean squares update (Bishop, 2006, Sect. 3.1.3). It avoids storing the entire training set, and its per-iteration complexity is independent of the training set size $m$.

The optimality condition Eq. (1) allows one to study the stability of linear regression. Ideally, the learned model parameters are insensitive to limited perturbations of the training set. A perturbed training set yields a feature matrix $\widetilde{{\mathbf X}} = {\mathbf X} + \Delta {\mathbf X}$ and a label vector $\widetilde{{\mathbf y}} = {\mathbf y} + \Delta {\mathbf y}$, with perturbation matrix $\Delta {\mathbf X}$ and vector $\Delta {\mathbf y}$. This gives the perturbed normal equations

\begin{equation}
\tag{2}
\widetilde{{\mathbf X}}^{\top} \widetilde{{\mathbf X}} \widetilde{{\mathbf w}}
= \widetilde{{\mathbf X}}^{\top} \widetilde{{\mathbf y}} \text{.}
\end{equation}

Matrix perturbation theory quantifies how much $\widetilde{{\mathbf w}}$ deviates from a solution $\widehat{{\mathbf w}}$ of the clean normal equations Eq. (1) (Golub and Loan, 2013, Sect. 5.3.6).

For a label-only perturbation, which replaces a single label of the training set with an outlier or another corrupted value, $\Delta {\mathbf X} = \mathbf{0}$, and the perturbed and clean normal equations share the coefficient matrix ${\mathbf X}^{\top}{\mathbf X}$. Assume ${\mathbf X}$ has full column rank, so ${\mathbf X}^{\top}{\mathbf X}$ is invertible and the solution is unique. Subtracting the clean normal equations Eq. (1) from the perturbed normal equations Eq. (2) yields

$$
\widetilde{{\mathbf w}} - \widehat{{\mathbf w}} = \big({\mathbf X}^{\top}{\mathbf X}\big)^{-1} {\mathbf X}^{\top} \Delta {\mathbf y} = {{\mathbf X}}^{+} \Delta {\mathbf y} \text{,}
$$

with the pseudoinverse ${{\mathbf X}}^{+}$. The effect of perturbing the training set is thereby quantified,

$$
\left\lVert \widetilde{{\mathbf w}} - \widehat{{\mathbf w}} \right\rVert_{2} \leq \left\lVert {{\mathbf X}}^{+} \right\rVert_{2}\, \left\lVert \Delta {\mathbf y} \right\rVert_{2} \text{.}
$$

When ${\mathbf X}$ has full column rank, the pseudoinverse has the norm $\left\lVert {{\mathbf X}}^{+} \right\rVert_{2} = 1/\sqrt{\lambda_{\min}({\mathbf X}^{\top}{\mathbf X})}
= \sqrt{\kappa\!\left({\mathbf X}^{\top}{\mathbf X}\right)}/\left\lVert {\mathbf X} \right\rVert_{2}$, where $\kappa\!\left({\mathbf X}^{\top}{\mathbf X}\right) = \lambda_{\max}/\lambda_{\min}$ is the condition number of ${\mathbf X}^{\top}{\mathbf X}$ (see pseudoinverse). In the setting of Fig. 1, ${\mathbf X} = \mathbf{1}$ (the all-ones vector), so ${{\mathbf X}}^{+} \Delta {\mathbf y} =
(1/m) \sum_{r=1}^{m} \Delta y^{(r)}$. Fig. 2 illustrates such a label-only perturbation of the training set from Fig. 1: a single label perturbation shifts $\widehat{w}$ by one third of it, as Fig. 2 shows.

![Figure 2 of the entry linreg](https://dictionaryofml.org/terms/fig/linreg-fig2.svg)

*Figure 2: The perturbation $\Delta y^{(3)} = 6$ (red arrow) replaces the third label of the training set from Fig. 1 with an outlier: filled circles denote the original data points, as in Fig. 1, and the open circle the perturbed value. The outlier raises the average of the labels: the learned model parameter shifts by $\widetilde{w} - \widehat{w} = 6/3 = 2$, from the clean solution $\widehat{w} = 3$ (dashed line) to the perturbed solution $\widetilde{w} = 5$ (dotted line)*

The effect of perturbations can also be studied on the level of a specific optimization method. For example, the update of online gradient descent becomes

$$
\begin{aligned}
{\mathbf w}^{(t+1)}
&= {\mathbf w}^{(t)} + 2\eta \big(
\widetilde{y}^{(t)} - \big({\mathbf w}^{(t)}\big)^{\top}
\widetilde{{\mathbf x}}^{(t)} \big) \widetilde{{\mathbf x}}^{(t)} \\
&= {\mathbf w}^{(t)} + 2\eta \big(
y^{(t)} - \big({\mathbf w}^{(t)}\big)^{\top}
{\mathbf x}^{(t)} \big) {\mathbf x}^{(t)} + {\boldsymbol \varepsilon}^{(t)} \text{.}
\end{aligned}
$$

Here, ${\boldsymbol \varepsilon}^{(t)}$ is a perturbation term that depends on the data point perturbation $\left( \Delta {\mathbf x}^{(t)},\Delta y^{(t)} \right)$ and the current model parameters ${\mathbf w}^{(t)}$. A bounded perturbation term moves each iterate by a bounded amount. The perturbed online gradient descent update has no fixed point in general, since a single perturbed data point leaves $\widetilde{{\mathbf w}}$ unchanged only if $\widetilde{{\mathbf w}}$ interpolates it. The gradient descent step operator $\mathcal{F}^{(\eta)}$ computed on the perturbed training set instead has as its fixed points exactly the solutions of the perturbed normal equations Eq. (2), by the same computation as for the clean training set. Among them is the perturbed solution $\widetilde{{\mathbf w}}$, whose displacement from $\widehat{{\mathbf w}}$ was computed above.

**See also:** regression, linear model, least squares, squared error loss, [empirical risk minimization](https://dictionaryofml.org/terms/erm.html), [gradient descent](https://dictionaryofml.org/terms/gd.html), ridge regression, least absolute shrinkage and selection operator, pseudoinverse.

## References

- Boyd and Vandenberghe (2004). Convex Optimization. *Cambridge Univ. Press*. <https://doi.org/10.1017/CBO9780511804441>
- Strang (2016). Introduction to Linear Algebra. *Wellesley-Cambridge Press*.
- Meyer (2000). Matrix Analysis and Applied Linear Algebra. *SIAM*.
- Golub and Loan (2013). Matrix Computations. *The Johns Hopkins Univ. Press*.
- Hastie et al. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction. *Springer Science+Business Media*. <https://doi.org/10.1007/978-0-387-84858-7>
- Lehmann and Casella (1998). Theory of Point Estimation. *Springer-Verlag*. <https://doi.org/10.1007/b98854>
- Papoulis and Pillai (2002). Probability, Random Variables, and Stochastic Processes. *McGraw-Hill Higher Education*.
- Bishop (2006). Pattern Recognition and Machine Learning. *Springer Science+Business Media*. <https://doi.org/10.1007/978-0-387-45528-0>
- Pedregosa et al. (2011). Scikit-learn: Machine Learning in Python. *Journal of Machine Learning Research*. <https://www.jmlr.org/papers/v12/pedregosa11a.html>
- Paszke et al. (2019). PyTorch: An Imperative Style, High-Performance Deep Learning Library. *Adv. Neural Inf. Process. Syst.*. <https://proceedings.neurips.cc/paper_files/paper/2019/hash/bdbca288fee7f92f2bfa9f7012727740-Abstract.html>
- Bertsekas (2016). Nonlinear Programming. *Athena Scientific*.
- Nesterov (2004). Introductory Lectures on Convex Optimization: A Basic Course. *Kluwer Academic*. <https://doi.org/10.1007/978-1-4419-8853-9>

## Cite this entry

```bibtex
@incollection{dictml_linreg,
  author = {Jung, Alexander},
  editor = {Olioumtsevits, Konstantina and Schnoor, Ekkehard},
  title = {linear regression},
  booktitle = {Dictionary of Applied Machine Learning},
  publisher = {Springer Nature Singapore},
  year = {2026},
  doi = {10.1007/978-981-92-2012-0_13-1},
  isbn = {978-981-92-2012-0},
  url = {https://doi.org/10.1007/978-981-92-2012-0_13-1}
}
```