Linear regression is used to find linear relationships in data. In a set of data, it is the line (or plane) of best fit. You have probably seen linear regression in a scientific or financial graph where there is a set of plotted points and a line through the points.

To use linear regression, you need data with one or more features that are being used to predict a \(1\)-dimensional outcome. For example, you could use location (2 coordinates), square feet and number of bathrooms to predict the price of a house. You will need at least one data point for each feature.

Math behind Linear Regression

Suppose you have data with \(n\) features that you use to predict a real valued outcome. Suppose that you also have \(m \geq n\) data points. Data point \(d_i\) will have feature values \((a_{1i}, a_{2i}, \dots, a_{ni})\) and outcome \(b_i.\) To approximate the data points with a linear function, you need values \(x_1, x_2, \dots, x_n\) such that \begin{align} & a_{11}x_1 + a_{12}x_2 + \dots + a_{1n}x_n = b_1 \\ & a_{21}x_1 + a_{22}x_2 + \dots + a_{2n}x_n = b_2 \\ & \vdots \\ & a_{m1}x_1 + a_{m2}x_2 + \dots + a_{mn}x_n = b_m \\ \end{align} Once you have found \(x_1, x_2, \dots, x_n,\) given a data point \(d = (c_1, c_2, \dots, c_n),\) you can predict the outcome is \[c_1x_1 + c_2x_2 + \dots + c_nx_n\] Rewriting the formula for the known data in matrix notation, we want to solve for \(x_1, x_2, \dots, x_n\) in the formula \[ \begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & & & \\ a_{m1} & a_{m2} & \dots & a_{mn} \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{bmatrix} \] Writing the \((a_{ij})\) matrix as \(A,\) the \((x_i)\) vector as \(\overrightarrow{x},\) and the \((b_i)\) vector as \(\overrightarrow{b},\) we can write this as \[A\overrightarrow{x}=\overrightarrow{b}\] Since \(m \geq n\) there will be \(0\) or \(1\) solutions. Most likely there will be \(0\) solutions if \(m > n.\)

In the case that there is \(1\) solution \(\overrightarrow{x},\) that is the solution you should use to estimate outcomes.

In the case that there are \(0\) solutions, the linear regression algorithm says you should use the vector \(\hat{x}\) which minimizes the square error. In this case, given a solution \(\overrightarrow{x},\) the error term for the \(i\)th data point is \(e_i(\overrightarrow{x})\) defined by \[e_i(\overrightarrow{x}) = (a_{i1}x_1 + a_{i2}x_2 + \dots + a_{in}x_n) - b_i\] and the squared error of the vector \(\overline{x} = (x_1, x_2, \dots, x_n)\) is the sum of the squared errors for each data point: \[e_1(\overrightarrow{x})^2 + e_2(\overrightarrow{x})^2 + \dots + e_m(\overrightarrow{x})^2\] Since the squared error is positive quadratic there is exactly one minimum, \(\hat{x}.\) It is known from linear algebra that this minimum is \[\hat{x} = (A^\intercal A)^{-1} A^\intercal \overrightarrow{b} \]

2D Case

In the case that you have one input and one output, you are trying to find the line of best fit. You need \(m > 1,\) or \(m \geq 2\) data points. Data point \(d_i\) will have feature value \((a_i\) and output \(b_i.\)

The slope-intercept form of the equation of a line is \[y = mx + b\] To match the notation used above, we will replace \(m\) by \(x_1\) and \(b\) by \(x_2.\) So, to find an outcome that best fits the data, we need to find a slope \(x_1\) and a \(y\)-intercept \(x_2.\) To incorporate the constant into the equation, we can add a constant feature that has no effect on the outcome. This feature will always have input \(1.\) So, data point \(d_i\) will have features \((a_i, 1)\) and outputs \(b_i.\) (You can use this trick to incorporate a constant in any number of dimensions.\)

We want to find values \(x_1\) and \(x_2\) such that \begin{align} & x_1a_1 + x_2 = b_1 \\ & x_1a_2 + x_2 = b_2 \\ & \vdots \\ & x_1a_m + x_2 = b_m \\ \end{align} Written in matrix form, this is \[ \begin{bmatrix} a_1 & 1 \\ a_2 & 1 \\ \vdots & \\ a_m & 1 \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{bmatrix} \] The least squares solution is \(\hat{x} = (A^\intercal A)^{-1} A^\intercal \overrightarrow{b}.\) In this case, \[A^\intercal A = \begin{bmatrix} a_1 & a_2 & \dots & a_m \\ 1 & 1 & \dots & 1\\ \end{bmatrix} \begin{bmatrix} a_1 & 1 \\ a_2 & 1 \\ \vdots & \\ a_m & 1 \\ \end{bmatrix} = \begin{bmatrix} a_1^2+a_2^2+\dots+a_m^2 & a_1+a_2+\dots+a_m \\ a_1+a_2+\dots+a_m & m \\ \end{bmatrix} \] We shall multiply and divide by \(m\) so that we can write the notation more concisely. Namely, the mean of the squares of the data is \[\overline{a^2} = \frac{a_1^2+a_2^2+\dots+a_m^2}{m}\] and the mean of the data is \[\overline{a} = \frac{a_1+a_2+\dots+a_m}{m}\] Using this notation, we have \[A^\intercal A = m\begin{bmatrix} \overline{a^2} & \overline{a} \\ \overline{a} & 1 \\ \end{bmatrix}\] Therefore the inverse is \[(A^\intercal A)^{-1} = \frac{1}{m(\overline{a^2}-\overline{a}^2)} \begin{bmatrix} 1 & -\overline{a} \\ -\overline{a} & \overline{a^2} \\ \end{bmatrix} \] Let \(\sigma_a^2\) be the variance in the data. So, \[\sigma_a^2 = \overline{a^2}-\overline{a}^2\] Then \[(A^\intercal A)^{-1} = \frac{1}{m\sigma_a^2} \begin{bmatrix} 1 & -\overline{a} \\ -\overline{a} & \overline{a^2} \\ \end{bmatrix} \] The rest of the expression is \[A^\intercal \overrightarrow{b} = \begin{bmatrix} a_1 & a_2 & \dots & a_m \\ 1 & 1 & \dots & 1\\ \end{bmatrix} \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{bmatrix} = \begin{bmatrix} a_1 b_1 + a_2 b_2 + \dots a_m b_m \\ b_1 + b_2 + \dots b_m \\ \end{bmatrix} = m \begin{bmatrix} \overline{ab} \\ \overline{b} \\ \end{bmatrix} \] Therefore, \begin{align} \hat{x} & = (A^\intercal A)^{-1} A^\intercal \overrightarrow{b} \\ & = \frac{1}{m\sigma_a^2} \begin{bmatrix} 1 & -\overline{a} \\ -\overline{a} & \overline{a^2} \\ \end{bmatrix} \cdot m \begin{bmatrix} \overline{ab} \\ \overline{b} \\ \end{bmatrix} \\ & = \frac{1}{\sigma_a^2} \begin{bmatrix} \overline{ab}-\overline{a}\cdot\overline{b} \\ -\overline{a}\cdot\overline{ab} + \overline{a^2}\cdot\overline{b} \\ \end{bmatrix} \\ & = \frac{1}{\sigma_a^2} \begin{bmatrix} \text{Cov}(a,b) \\ -\overline{a}(\text{Cov}(a,b)+\overline{a}\cdot\overline{b}) + (\sigma_a^2 + \overline{a}^2)\cdot\overline{b} \\ \end{bmatrix} \\ & = \begin{bmatrix} \text{Cov}(a,b)/\sigma_a^2 \\ \overline{b} - \overline{a}\text{Cov}(a,b)/\sigma_a^2 \\ \end{bmatrix} \\ \end{align} This solves for the slope and \(y\)-intercept. \begin{align} & \text{Slope: } x_1 = \text{Cov}(a,b)/\sigma_a^2 & \text{y-intercept: } x_2 = \overline{b} - \overline{a}\text{Cov}(a,b)/\sigma_a^2 \end{align}

Pearson's Correlation Coefficient - \(r\) Score

One way to measure how well the line fits the data is with Pearson's Correlation Coefficient, which is represented with an \(r.\) The value \(r\) measures the linear correlation between the \(x\) and \(y\) values, and can be computed as \[r = \frac{n\sum x_i y_i - \overline{x}\overline{y}}{\sqrt{\left(n\sum x_i^2 - \overline{x}^2\right)\left(n\sum y_i^2 - \overline{y}^2\right)}}\]