The main task in a Linear Regression problem is to find the coefficients which, in 1-dimension case, are able to create a line where the difference between the datapoints and the created line is as small as possible. These differences are called "residuals" and they need to be minimized. This problem can be solved by using a procedure called Ordinary Least Squares (OLS) where the sum of all the squared residuals mentioned before are minimized. In case of $n$ datapoints, the $i^{th}$ datapoint $x_i$ is associated to the output $y_i$ and the line formula is $y_i=m{x}_i+b+{ϵ}_i$ where ${ϵ}_i$ is an an error term which takes into account all the variations of the output which are not associated to the known input variables. Mathematically, the OLS in 1-dimension can be written as: $$\begin{array}{}\text{(1)}& \underset{\hat{m},\hat{b}}{min}\sum _{i=1}^n(y_i-\hat{m}{x}_i-\hat{b}{)}^2\end{array}$$ Note that the $m$ and $b$ coefficients found are estimates so we call them respectively $\hat{m}$ and $\hat{b}$ however for the rest of this explanation they simply be called $m$ and $b$ to simplify the notation. In order to find these coefficients we can do a derivative of this quadratic function and set it equal to 0. This is a common way of finding the minima and maxima of a function. For a quadratic function, like in this case, there is only a minimum or maximum point. However, in this case the function is also convex so the derivative leads to the minimum point.

Find $b$

Let's start with the intercept $b$. We can do the derivative with respect to $b$ remembering that the derivative of a quadratic function $f(x{)}^2$ is 2 multiplied by the function and then multiplied by the derivative of the function (i.e. $2f(x)f^{\prime }(x)$). The first step is to set the derivative to 0: $$\begin{array}{}\text{(2)}& \frac{\partial }{\partial b}[\sum _{i=1}^n(y_i-m{x}_i-b{)}^2]=0\end{array}$$ Remembering that, for a datapoint $i$, $\frac{\partial }{\partial b}(y_i-m{x}_i-b{)}^2=2(y_i-m{x}_i-b)(-1)$, you can see that the derivative of the terms within the brackets is simply -1 because $y$ and $mx$ do not have terms in $b$ and the derivative of $-b$ with respect to $b$ is simply -1.

Equation 3 shows that the derivative can be done for all the $n$ points and, given that they are not interacting with each other, all these derivatives can simply be summed together. Also, given that we set this function equal to 0 we can remove the moltiplicative factor -2 because it is not going to change the point where the function is 0. $$\begin{array}{}\text{(3)}& \frac{\partial }{\partial b}[\sum _{i=1}^n(y_i-m{x}_i-b{)}^2]=-2\sum _{i=1}^n(y_i-m{x}_i-b)=0\end{array}$$

Equation 4 shows that we can split the sum in the three different terms. In particular, in the second term we can put $m$ outside the sum given that does not depend on $i$ and the third term is simply the coefficient $-b$ summed $n$ times so we can simply write $-nb$. $$\begin{array}{}\text{(4)}& \sum _{i=1}^n(y_i-m{x}_i-b)=\sum _{i=1}^n{y}_i-m\sum _{i=1}^n{x}_i-nb=0\end{array}$$

Equation 5 shows that if we put the first 2 terms on the other side of the equation and divide by $-n$ then we obtain an equation in $b$. In particular $b$ is equal to the sum of all the $y$ values divided by $n$ which is the mean of the $y$ terms, minus $m$ multiplied by the mean of the $x$ terms. $$\begin{array}{}\text{(5)}& b=\frac{\sum _{i=1}^n{y}_i}{n}-m\frac{\sum _{i=1}^n{x}_i}{n}\end{array}$$ We can calculate the mean of the $x$ terms and $y$ terms so the only missing value is $m$.

Find $m$

To find $m$, now we do the derivative with respect to $m$ and set it to 0: $$\begin{array}{}\text{(6)}& \frac{\partial }{\partial m}[\sum _{i=1}^n(y_i-m{x}_i-b{)}^2]=0\end{array}$$ In this case, the derivative of the terms within the brackets for a datapoint $i$ is $-x_i$ because $y$ and $b$ do not have terms in $m$ but the derivative of $-m{x}_i$ with respect to $m$ is $-x_i$. So, $\frac{\partial }{\partial m}(y_i-m{x}_i-b{)}^2=2(y_i-m{x}_i-b)(-x_i)$. Equation 7 is similar to the previous case, we can again do the sum of all the derivatives and remove -2. $$\begin{array}{}\text{(7)}& \frac{\partial }{\partial m}[\sum _{i=1}^n(y_i-m{x}_i-b{)}^2]=-2\sum _{i=1}^n(y_i-m{x}_i-b)x_i=0\end{array}$$ Equation 8 shows that we can multiply $x_i$ for all the terms within the brackets and then split the sums similarly to what we did before. In this case note that the second term is a squared term while the third term is not anymore $-nb$ given that it is multiplied by $x_i$ which depends on $i$. $$\begin{array}{}\text{(8)}& \sum _{i=1}^n(y_i-m{x}_i-b)x_i=\sum _{i=1}^n{y}_i{x}_i-m\sum _{i=1}^n{x}_i^2-b\sum _{i=1}^n{x}_i=0\end{array}$$ Equation 9 shows that we can replace the $b$ value using the equation found before and now this equation contains only terms in $m$, $x$ and $y$. $$\begin{array}{}\text{(9)}& \sum _{i=1}^n{y}_i{x}_i-m\sum _{i=1}^n{x}_i^2-(\frac{\sum _{i=1}^n{y}_i}{n}-m\frac{\sum _{i=1}^n{x}_i}{n})\sum _{i=1}^n{x}_i=0\end{array}$$ In equation 10 we move all the terms in $m$ on the other side of the equation. $$\begin{array}{}\text{(10)}& \sum _{i=1}^n{y}_i{x}_i-(\frac{\sum _{i=1}^n{y}_i\sum _{i=1}^n{x}_i}{n})=m\sum _{i=1}^n{x}_i^2-m\frac{(\sum _{i=1}^n{x}_i{)}^2}{n}\end{array}$$ Then, in equation 11 we divide for all the terms that are not $m$ and we can obtain a formula to find $m$ which only depends on the $x$ and $y$ data. This is what we are looking for given that only by knowing the datapoints we can find the optimal $m$ value. $$\begin{array}{}\text{(11)}& m=\frac{\sum _{i=1}^n{y}_i{x}_i-(\frac{\sum _{i=1}^n{y}_i\sum _{i=1}^n{x}_i}{n})}{\sum _{i=1}^n{x}_i^2-\frac{(\sum _{i=1}^n{x}_i{)}^2}{n}}\end{array}$$ The formulas are a bit long but we can simplify them. In particular we can substitute the mean values of $x$ and $y$ respectively with $\overline{x}$ and $\overline{y}$ (i.e. $\overline{x}=\frac{\sum _{i=1}^n{x}_i}{n}$ and $\overline{y}=\frac{\sum _{i=1}^n{y}_i}{n}$) and this simplifies the two equations for $b$ and $m$. $$\begin{array}{}\text{(12)}& m=\frac{\sum _{i=1}^n{y}_i{x}_i-\overline{x}\sum _{i=1}^n{y}_i}{\sum _{i=1}^n{x}_i^2-\overline{x}\sum _{i=1}^n{x}_i}=\frac{\sum _{i=1}^n{y}_i(x_i-\overline{x})}{\sum _{i=1}^n{x}_i(x_i-\overline{x})}\end{array}$$ $$\begin{array}{}\text{(13)}& b=\overline{y}-m\overline{x}\end{array}$$ Equations 12 and 13 show the calculations that need to be performed, using $x$ and $y$ datapoints, to obtain the optimal coefficient values for $m$ and $b$. Click the button below to see the implementation of these formulas in Python.

Linear Regression Code

The linear regression problem can also be defined for more than one indepedent variable and the equations are typically written in matrix form. For a generic point with output $y_i$ then there are $p$ different input variables and each of them is associated to a specific $\beta$ coefficient and then there is also the coefficient ${\beta }_0$ not associated to any variable because this is the intercept. The equation for a generic $i$ data point is: $$\begin{array}{}\text{(14)}& y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_p{x}_{ip}+{ϵ}_i\end{array}$$

This can be written in matrix form, taking into account all $n$ data points, in this way: $$\begin{array}{}\text{(15)}& \mathit{y}=\mathit{X}\mathit{\beta }+\mathit{ϵ}\end{array}$$ Where: $$\begin{array}{}\text{(16)}& \mathit{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right]\mathit{\beta }=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right]\mathit{ϵ}=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_n\end{array}\right]\end{array}$$ $$\begin{array}{}\text{(17)}& \mathit{X}=\left[\begin{array}{cccc}1& x_{11}& \cdots & x_{1p}\\ 1& x_{21}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{n1}& \cdots & x_{np}\end{array}\right]\end{array}$$

Given that there are $n$ data points, $\mathit{y}$ is a vector of $n$ by 1 dimensions and also the error term $\mathit{ϵ}$ is a vector of $n$ by 1 dimensions. On the other hand, $\mathit{\beta }$ as previously mentioned is a vector of $p+1$ by 1 dimensions where $p$ is the number of input variables and then we add the intercept. For the data, this is a matrix of $n$ by $p+1$ dimensions given that there are $n$ data points but for each data point there is a multiplication for the $\mathit{\beta }$ vector of coefficients which is $p+1$ dimensions. In the matrix one column is all ones, that is the intercept column which multiplies the ${\beta }_0$ coefficient and this is a compact way to include the intercept in the input matrix.

The Ordinary Least Squares (OLS) defined in matrix form is similar to the 1-dimension case. The equation below shows how the problem is formally posed and that is essentially the same equation used previously where the sum of all the squared residuals is calculated and then minimized but in this case the aim is to find a vector $\mathit{\beta }$ with all its $p+1$ coefficients. $$\begin{array}{}\text{(18)}& \underset{\mathit{\beta }}{min}\Vert \mathit{y}\mathbf{-}\mathit{X}\mathit{\beta }{\Vert }^2\end{array}$$ The two vertical bars at the left and right of the equation, mean that the euclidean norm of the vector is calculated. Given a vector $\mathit{a}$, the euclidean norm of $\mathit{a}$ is defined as $\Vert \mathit{a}\Vert =\sqrt{a_1^2+a_2^2+\cdots +a_p^2}=\sqrt{{\mathit{a}}^{\mathit{T}}\mathit{a}}$. So, the euclidean norm can be written as the multiplication between the vector $\mathit{a}$ transpose and the original vector $\mathit{a}$. The euclidean norm has the square root but the function that we are trying to minimize is squared so that we are simply considering the sum of the squared terms without the need of calculating the square root.

To find the optimal coefficients, Equation 19 shows that we can write the equation using the transpose and then we do the derivative and set it equal to 0. In this case the derivative is in respect of a vector $\mathit{\beta }$ and not single variables as previously shown. $$\begin{array}{}\text{(19)}& \frac{\partial }{\partial \mathit{\beta }}\Vert \mathit{y}\mathbf{-}\mathit{X}\mathit{\beta }{\Vert }^2=\frac{\partial }{\partial \mathit{\beta }}\mathbf{(}\mathit{y}\mathbf{-}\mathit{X}\mathit{\beta }{\mathbf{)}}^{\mathit{T}}\mathbf{(}\mathit{y}\mathbf{-}\mathit{X}\mathit{\beta }\mathbf{)}=0\end{array}$$ Equation 20 focuses on performing the multiplications of the terms within the brackets by taking into account that the transpose of two matrices (or a matrix and a vector) are equal to the multiplication of the matrices (or matrix and vector) both transposed but they are flipped (i.e. $\mathbf{(}\mathit{A}\mathit{B}{\mathbf{)}}^{\mathit{T}}\mathbf{=}{\mathit{B}}^{\mathit{T}}{\mathit{A}}^{\mathit{T}}$). $$\begin{array}{}\text{(20)}& \frac{\partial }{\partial \mathit{\beta }}\mathbf{(}\mathit{y}\mathbf{-}\mathit{X}\mathit{\beta }{\mathbf{)}}^{\mathit{T}}\mathbf{(}\mathit{y}\mathbf{-}\mathit{X}\mathit{\beta }\mathbf{)}=\frac{\partial }{\partial \mathit{\beta }}({\mathit{y}}^{\mathit{T}}\mathit{y}\mathbf{-}{\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{y}\mathbf{-}{\mathit{y}}^{\mathit{T}}\mathit{X}\mathit{\beta }\mathbf{+}{\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}\mathit{\beta })=0\end{array}$$ If you look at all these 4 terms and apply matrix multiplications you can see that dimensionally all of them are a single dimension so that we are essentially obtaining a single value out of this OLS equation and this value is the minimum error possible given the optimal coefficients that will be found. We can do the derivative of these 4 terms by remembering some of the derivative identities in respect of a vector.

• Given 2 vectors $\mathit{a}$ and $\mathit{x}$, $\frac{\partial {\mathit{x}}^{\mathit{T}}\mathit{a}}{\partial \mathit{x}}=\frac{\partial {\mathit{a}}^{\mathit{T}}\mathit{x}}{\partial \mathit{x}}={\mathit{a}}^{\mathit{T}}$
• Given a matrix $\mathit{A}$ and a vector $\mathit{x}$, $\frac{\partial {\mathit{x}}^{\mathit{T}}\mathit{A}\mathit{x}}{\partial \mathit{x}}=2{\mathit{x}}^{\mathit{T}}\mathit{A}$ ($\mathit{A}$ must be symmetric)

Now let's use these identities:

• The first term does not have $\mathit{\beta }$ terms, so the derivative is 0: $\frac{\partial {\mathit{y}}^{\mathit{T}}\mathit{y}}{\partial \mathit{\beta }}=0$
• The second term is ${\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{y}$. This can be seen as ${\mathit{\beta }}^{\mathit{T}}\mathbf{(}{\mathit{X}}^{\mathit{T}}\mathit{y}\mathbf{)}$, using the previous identity in this case the result is ${\mathit{y}}^{\mathit{T}}\mathit{X}$. : $\frac{\partial {\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{y}}{\partial \mathit{\beta }}={\mathit{y}}^{\mathit{T}}\mathit{X}$
• The third term is similar to the second but in this case is it can be seen as $\mathbf{(}{\mathit{y}}^{\mathit{T}}\mathit{X}\mathbf{)}\mathit{\beta }$ but the result is still the same: $\frac{\partial {\mathit{y}}^{\mathit{T}}\mathit{X}\mathit{\beta }}{\partial \mathit{\beta }}={\mathit{y}}^{\mathit{T}}\mathit{X}$
• The fourth term is ${\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}\mathit{\beta }$. As shown in the second identity, this derivative is equal to $2{\mathit{\beta }}^{\mathit{T}}$ and the symmetric matrix which in this case is ${\mathit{X}}^{\mathit{T}}\mathit{X}$: $\frac{\partial {\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}\mathit{\beta }}{\partial \mathit{\beta }}=2{\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}$

Equation 21 shows that the sum of all these terms should be equal to 0: $$\begin{array}{}\text{(21)}& 0\mathbf{-}{\mathit{y}}^{\mathit{T}}\mathit{X}\mathbf{-}{\mathit{y}}^{\mathit{T}}\mathit{X}\mathbf{+}\mathbf{2}{\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}\mathbf{=}0\end{array}$$ Equation 22 shows that we can remove the first term and sum together the second and third term which both are negative and then the equation has only 2 terms. $$\begin{array}{}\text{(22)}& \mathbf{-}\mathbf{2}{\mathit{y}}^{\mathit{T}}\mathit{X}\mathbf{+}\mathbf{2}{\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}\mathbf{=}0\end{array}$$ We can divide the both by 2 and then move the first term on the other side of the equation and this leads to Equation 23. $$\begin{array}{}\text{(23)}& {\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}\mathbf{=}{\mathit{y}}^{\mathit{T}}\mathit{X}\end{array}$$ We can transpose both sides as you can see from Equation 24. $$\begin{array}{}\text{(24)}& \mathbf{(}{\mathit{\beta }}^{\mathit{T}}{\mathit{X}}^{\mathit{T}}\mathit{X}{\mathbf{)}}^{\mathit{T}}\mathbf{=}\mathbf{(}{\mathit{y}}^{\mathit{T}}\mathit{X}{\mathbf{)}}^{\mathit{T}}\end{array}$$ Remember that ${\mathit{X}}^{\mathit{T}}\mathit{X}$ is symmetric so the transpose of this matrix has not effect and returns again ${\mathit{X}}^{\mathit{T}}\mathit{X}$. Also, as mentioned before, for two matrices $\mathit{A}$ and $\mathit{B}$, if they get transposed $\mathbf{(}\mathit{A}\mathit{B}{\mathbf{)}}^{\mathit{T}}\mathbf{=}{\mathit{B}}^{\mathit{T}}{\mathit{A}}^{\mathit{T}}$. This means that in Equation 24 we can swap the ${\mathit{\beta }}^{\mathit{T}}$ and ${\mathit{X}}^{\mathit{T}}\mathit{X}$ terms and $\mathit{\beta }$ is not anymore transposed while ${\mathit{X}}^{\mathit{T}}\mathit{X}$ remains the same. On the right hand of the equation we can swap ${\mathit{y}}^{\mathit{T}}$ and $\mathit{X}$ and in this case they become respectively $\mathit{y}$ and ${\mathit{X}}^{\mathit{T}}$. You can see this in Equation 25. $$\begin{array}{}\text{(25)}& {\mathit{X}}^{\mathit{T}}\mathit{X}\mathit{\beta }\mathbf{=}{\mathit{X}}^{\mathit{T}}\mathit{y}\end{array}$$ Now, remembering that the inverse of matrix $\mathit{A}$ (i.e ${\mathit{A}}^{\mathbf{-}\mathbf{1}}$) which multiplies $\mathit{A}$ is equal to the identity matrix $\mathit{I}$ which is composed only by ones on the diagonal. The identity matrix behaves like the number 1 in a single dimension, if you multiply a number by 1 you get the same number. In the same way, if you have an indentity matrix which multiplies a vector, you simply obtain the vector. So the last thing to do is to multiply both sides by $\mathbf{(}{\mathit{X}}^{\mathit{T}}\mathit{X}{\mathbf{)}}^{\mathbf{-}\mathbf{1}}$. In this way the left hand side is the vector $\mathit{\beta }$ while the right hand side is composed by only terms in $\mathit{X}$ and $\mathit{y}$ as you can see from Equation 26. $$\begin{array}{}\text{(26)}& \mathit{\beta }\mathbf{=}\mathbf{(}{\mathit{X}}^{\mathit{T}}\mathit{X}{\mathbf{)}}^{\mathbf{-}\mathbf{1}}{\mathit{X}}^{\mathit{T}}\mathit{y}\end{array}$$ This is a closed-form solution which you can apply to obtain an estimate of the $\mathit{\beta }$ coefficients. This method assumes that we can do the inverse of the matrix, this is not always the case but in this case we are assuming that we can compute it. Click the button below to see the implementation of these formulas in Python.

Linear Regression Code