Linear regression

Formulation

Let $\mathcal{X}\subseteq {\mathbb{R}}^d$ and $\mathcal{Y}\subseteq \mathbb{R}$ be the (input) feature and target space respectively. A linear regression model takes the form

$$\begin{array}{rl}y& ={\beta }_0+\sum _{i=1}^d{\beta }_i{x}_i+\epsilon \\ & ={\beta }^{\top }\mathbf{x}+\epsilon \end{array}$$

where

$$\begin{array}{rl}\beta & ={\left[\begin{array}{cccc}{\beta }_0& {\beta }_1& \cdots & {\beta }_d\end{array}\right]}^{\top }\in {\mathbb{R}}^{d+1},\\ \mathbf{x}& ={\left[\begin{array}{cccc}1& x_1& \cdots & x_d\end{array}\right]}^{\top }\in {\mathbb{R}}^{d+1}.\end{array}$$

Linear regression can also be formulated in a matrix form. Suppose there are $N$ total instances, i.e., $\mathcal{X}=\{{\mathbf{x}}^{(i)}{\}}_{i=1}^N$ and $\mathcal{Y}=\{y^{(i)}{\}}_{i=1}^N$. Then,

$$\mathbf{y}=\mathbf{X}\beta +\epsilon$$

where

$$\begin{array}{rl}\mathbf{y}& ={\left[\begin{array}{cccc}{y}^{(1)}& y^{(2)}& \cdots & y^{(N)}\end{array}\right]}^{\top }\in {\mathbb{R}}^N,\\ \mathbf{X}& =\left[\begin{array}{c}{{\mathbf{x}}^{(1)}}^{\top }\\ {{\mathbf{x}}^{(2)}}^{\top }\\ \cdots \\ {{\mathbf{x}}^{(N)}}^{\top }\end{array}\right]=\left[\begin{array}{ccccc}1& x_1^{(1)}& x_2^{(1)}& \cdots & x_d^{(1)}\\ 1& x_1^{(2)}& x_2^{(2)}& \cdots & x_d^{(2)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_1^{(N)}& x_2^{(N)}& \cdots & x_d^{(N)}\end{array}\right]\in {\mathbb{R}}^{N\times (d+1)},\\ \beta & ={\left[\begin{array}{cccc}{\beta }_0& {\beta }_1& \cdots & {\beta }_d\end{array}\right]}^{\top }\in {\mathbb{R}}^{d+1},\\ \epsilon & ={\left[\begin{array}{cccc}{\epsilon }^{(1)}& {\epsilon }^{(2)}& \cdots & {\epsilon }^{(N)}\end{array}\right]}^{\top }\in {\mathbb{R}}^N.\end{array}$$

The values $\epsilon$ are called error terms, or sometimes noise, and captures all other factors that influence $y$ other than $\mathbf{x}$. It is common to assume that it follows a Gaussian distribution:

$$\epsilon \sim \mathcal{N}(0,{\sigma }^2).$$

Then, the likelihood function can be written as

$$p(y|\mathbf{x},\beta )=\mathcal{N}(y|\beta \mathbf{x},{\sigma }^2)$$

or equivalently

$$p(\mathbf{y}|\mathbf{X},\beta )=\prod _{i=1}^N\mathcal{N}\left(y^{(i)}|{\beta }^{\top }{\mathbf{x}}^{(i)},{{\sigma }^{(i)}}^2\right)$$

in multidimensional cases, assuming the data points are drawn independently from the distribution.

Basis Functions

The linear regression model assumes a linear relationship between $\mathbf{x}$ and $y$. However, in some cases the relationship may be nonlinear. To accommodate such cases, it is possible to apply a nonlinear transformation (called a basis function) $\varphi$ to the input $\mathbf{x}$, transforming the input into some other form $\varphi (\mathbf{x})$. As long as the parameters of $\varphi$ are fixed, the model remains linear in the parameters, even if it is not linear in the inputs. Common choices of basis functions include:

• Polynomial: ${\varphi }_i(\mathbf{x})=x_i^i$
• Gaussian: ${\varphi }_i(\mathbf{x})=\exp \left(-\frac{(x-{\mu }_i{)}^2}{2s^2}\right)$
• Sigmoidal: ${\varphi }_i(\mathbf{x})=\sigma \left(\frac{x-{\mu }_i}{s}\right)$

Parameter Estimation

Analytic

Maximum Likelihood

To maximise the likelihood, we can alternatively minimise the negative log likelihood.

$$\begin{array}{rl}\text{NLL}(\beta ,{\sigma }^2)& =-\log p(\mathbf{y}|\mathbf{X},\beta )\\ & =-\log \left(\prod _{i=1}^N\mathcal{N}\left(y^{(i)}|{\beta }^{\top }{\mathbf{x}}^{(i)},{{\sigma }^{(i)}}^2\right)\right)\\ & =-\sum _{i=1}^N\log \mathcal{N}\left(y^{(i)}|{\beta }^{\top }{\mathbf{x}}^{(i)},{{\sigma }^{(i)}}^2\right)\\ & =-\sum _{i=1}^N\log \left(\sqrt{\frac{1}{2\pi {{\sigma }^{(i)}}^2}}\cdot \exp \left(-\frac{1}{2{{\sigma }^{(i)}}^2}(y^{(i)}-{\beta }^{\top }{\mathbf{x}}^{(i)}{)}^2\right)\right)\\ & =\frac{1}{2}\sum _{i=1}^N\left[(y^{(i)}-{\beta }^{\top }{\mathbf{x}}^{(i)}{)}^2\cdot \frac{1}{2{{\sigma }^{(i)}}^2}+\log \left(2\pi {{\sigma }^{(i)}}^2\right)\right]\end{array}$$

The MLE is the point where ${\nabla }_{\beta }\text{NLL}(\beta ,{\sigma }^2)=0$.

$$\begin{array}{rl}{\nabla }_{\beta }\text{NLL}(\beta ,{\sigma }^2)& ={\nabla }_{\beta }\frac{1}{2}\sum _{i=1}^N\left[(y^{(i)}-{\beta }^{\top }{\mathbf{x}}^{(i)}{)}^2\cdot \frac{1}{2{{\sigma }^{(i)}}^2}+\log \left(2\pi {{\sigma }^{(i)}}^2\right)\right]\\ & \propto \frac{1}{2}{\nabla }_{\beta }\sum _{i=1}^N{\left(y^{(i)}-{\beta }^{\top }{\mathbf{x}}^{(i)}\right)}^2\\ & =\frac{1}{2}{\nabla }_{\beta }{\Vert \mathbf{X}\beta -\mathbf{y}\Vert }^2\\ & =\frac{1}{2}{\nabla }_{\beta }(\mathbf{X}\beta -\mathbf{y}{)}^{\top }(\mathbf{X}\beta -\mathbf{y})\\ & =\frac{1}{2}{\nabla }_{\beta }\left({\beta }^{\top }{\mathbf{X}}^{\top }\mathbf{X}\beta -2{\beta }^{\top }{\mathbf{X}}^{\top }\mathbf{y}+{\mathbf{y}}^{\top }\mathbf{y}\right)\\ & =\frac{1}{2}\left(2{\beta }^{\top }({\mathbf{X}}^{\top }\mathbf{X})-2{\mathbf{y}}^{\top }\mathbf{X}\right)\\ & ={\beta }^{\top }({\mathbf{X}}^{\top }\mathbf{X})-{\mathbf{y}}^{\top }\mathbf{X}\end{array}$$

Solving the equation ${\nabla }_{\beta }\text{NLL}(\beta ,{\sigma }^2)=0$ yields

$$\begin{array}{r}({\mathbf{X}}^{\top }\mathbf{X}){\beta }^{\ast }={\mathbf{X}}^{\top }\mathbf{y}\end{array}$$

and therefore,

$${\beta }^{\ast }=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}\mathbf{Xy}.$$

This is also called the normal solution.

Mean Square Error

Alternatively, we can define a loss function and find the optimal point that minimises the loss. A common choice for the loss function is the mean square loss (MSE), which is given as follows:

$$\begin{array}{rl}\text{MSE}(\beta )& =\frac{1}{N}\sum _{i=1}^N\frac{1}{2}{\left({\beta }^{\top }{\mathbf{x}}^{(i)}-y^{(i)}\right)}^2\\ & \propto \frac{1}{2}{\Vert \mathbf{X}\beta -\mathbf{y}\Vert }_2^2\end{array}$$

Now it suffices to find the point such that ${\nabla }_{\beta }\text{MSE}(\beta )=0$.

$$\begin{array}{rl}{\nabla }_{\beta }\text{MSE}(\beta )& \propto {\nabla }_{\beta }{\Vert \mathbf{X}\beta -\mathbf{y}\Vert }_2^2\\ & =\frac{1}{2}{\nabla }_{\beta }(\mathbf{X}\beta -\mathbf{y}{)}^{\top }(\mathbf{X}\beta -\mathbf{y})\\ & =\frac{1}{2}{\nabla }_{\beta }\left({\beta }^{\top }{\mathbf{X}}^{\top }\mathbf{X}\beta -2{\beta }^{\top }{\mathbf{X}}^{\top }\mathbf{y}+{\mathbf{y}}^{\top }\mathbf{y}\right)\\ & =\frac{1}{2}\left(2{\beta }^{\top }({\mathbf{X}}^{\top }\mathbf{X})-2{\mathbf{y}}^{\top }\mathbf{X}\right)\\ & ={\beta }^{\top }({\mathbf{X}}^{\top }\mathbf{X})-{\mathbf{y}}^{\top }\mathbf{X}\end{array}$$

Notice that the equation is the same as in the MLE scenario. The optimal parameters are given as:

$${\beta }^{\ast }=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}\mathbf{Xy}.$$

Numerical

Solving the normal solution includes computing the inverse of a large matrix ${\mathbf{X}}^{\top }\mathbf{X}$ of size $(d+1)\times (d+1)$. This can be a very expensive computation if $d$, the input dimension, becomes large. In such cases, finding the optimal solution numerically can be an alternative.

• Gradient descent
  • Stochastic gradient descent

References

• Wikipedia contributors. (2024b, January 15). Linear regression. Wikipedia. https://en.wikipedia.org/wiki/Linear_regression
• Murphy, K. P. (2022). Probabilistic Machine Learning: An Introduction. MIT Press.
• Bishop, C. M. (2016). Pattern recognition and machine learning. Springer.
• Deisenroth, M. P., Faisal, A. A., & Ong, C. S. (2020). Mathematics for machine learning. Cambridge University Press.
• James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013). An introduction to statistical learning: with Applications in R. Springer Science & Business Media.