Mathematics Behind Linear Regression

These are my notes that cover the mathematical foundations of Linear Regression, taken while attending CSCI-UA 9473 - Foundations of Machine Learning at NYU Paris. They make use of probability theory, statistics, calculus, optimization, and linear algebra to formalize the concept of linear regression.

Introduction

Linear regression is one of the most widely used statistical techniques in data analysis and machine learning. It provides a simple and intuitive linear model for modeling the relationship between a response variable and a set of explanatory variables.

Simple Linear Regression

Least Squares Criterion

where $\mathbf{w} = \begin{pmatrix} w_0 \\ w_1 \end{pmatrix}$. We want to find $\mathbf{w}$ such that $J_{\tiny N}(\mathbf{w})$ is minimized.

Therefore, lets compute $\nabla J_{\tiny N}(\mathbf{w})$:

where:

Setting them equal to zero, we get the critical values $(w_0, w_1)$ that minimize $J_{\tiny N}(\mathbf{w})$:

Similarly:

Replacing the value of $w_0$, we get:

However, it is commonly written in terms of Pearson's sample correlation coefficient $r_{xy}$.

For that, a slight modification is required in equation $(6)$:

Even though we have added the part in green, we can show that the second term is zero as:

Therefore, from equation $(9)$, we get:

Here, $\sigma_y$ and $\sigma_x$ are standard deviations of $x$ and $y$.

$r_{xy}$ is the Pearson's sample correlation coefficient between $x$ and $y$ defined as:

$$\begin{equation} r_{xy} = \frac{\displaystyle\sum_{i=1}^N \left( y_i - \overline{y} \right) \left(x_i - \overline{x}\right)}{\displaystyle\sqrt{\sum_{i=1}^N \left( y_i - \overline{y} \right)^2} \displaystyle\sqrt{\sum_{i=1}^N \left( x_i - \overline{x} \right)^2}} \end{equation}$$

Multiple Linear Regression

In multiple linear regression, there are multiple independent variables in the equation rather than just one. The mathematics behind multiple linear regression involves using matrix operations to solve for the coefficients of the regression equation.

The regression equation is now represented as:

where $x_1, x_2, \cdots, x_d$ are the independent variables and $w_0, w_1, \cdots, w_d$ are the coefficients of the equation.

We define $\mathbf{w} = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_d \end{pmatrix}$ and augment each data point as $\mathbf{x} = \begin{pmatrix} 1 \\ x_1 \\ \vdots \\ x_d \end{pmatrix}$, so that the equation can be written as:

Least Squares Criterion

where:

$\mathbf{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_{\tiny N} \end{pmatrix}$ is the vector of target values.

$\Phi \mathbf{w} = \begin{pmatrix} \langle \mathbf{x_1}, \mathbf{w} \rangle \\ \langle \mathbf{x_2}, \mathbf{w} \rangle \\ \vdots \\ \langle \mathbf{x_{\tiny N}}, \mathbf{w} \rangle \end{pmatrix} = \begin{pmatrix} \mathbf{x_1}^T \mathbf{w} \\ \mathbf{x_2}^T \mathbf{w} \\ \vdots \\ \mathbf{x_{\tiny N}}^T \mathbf{w} \end{pmatrix} = \begin{pmatrix} \mathbf{x_1}^T \\ \mathbf{x_2}^T \\ \vdots \\ \mathbf{x_{\tiny N}}^T \end{pmatrix}\mathbf{w}$ is the vector of predictions.

$\Phi = \begin{pmatrix} \mathbf{x_1}^T \\ \mathbf{x_2}^T \\ \vdots \\ \mathbf{x_{\tiny N}}^T \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1d} \\ 1 & x_{21} & x_{22} & \cdots & x_{2d} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{\tiny{N}1} & x_{\tiny{N}2} & \cdots & x_{\tiny{N}d} \end{pmatrix}$ is called the design matrix. It is of size $N \times (d+1)$.

$\mathbf{w} = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_d \end{pmatrix}$ is the vector of weights.

We want to find the minimizer for $J_{\tiny N}(\mathbf{w})$, so we solve:

Note that we are implicitly assuming that $\Phi^T \Phi$ is an invertible matrix. If either the number of linearly independent examples is less than the number of features, or if the features are not linearly independent, then $\Phi^T \Phi$ is not invertible.

This least squared solution is an estimator for $\mathbf{w}$. It represents the vector normal to the hyperplane that minimizes the sum of squared errors between the target values and the predictions.

However, this is not enough, we need to be sure of how confident we are in our predictions.

---

Homoscedastic Model

Assume that target values $y_i$ are represented by random variables $Y_1, Y_2, \cdots, Y_{\tiny N}$ that are independent and identically distributed (i.i.d.). Then,

where $\mathbf{w^{*}}$ is the vector of true weights that generates the data and $\mathbf{\epsilon}_i \sim \mathcal{N}(0, \sigma^2)$ is an independent gaussian noise term.

We can write $\mathbf{Y}$ as:

Here, $\Large \mathbf{\epsilon}$ $\sim \mathcal{N}(\mathbf{0}, \sigma^2 I_{\tiny N})$ where $I_{\tiny N}$ is the $N \times N$ identity matrix.

We can show this as:

$\mathbb{E} \left[\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_{\tiny N} \end{pmatrix} \right] = \left[\begin{array}{c} \mathbb{E} \left( \epsilon_1 \right) \\ \mathbb{E} \left( \epsilon_2 \right) \\ \vdots \\ \mathbb{E} \left(\epsilon_{\tiny N}\right) \end{array}\right] = \mathbf{0}$

$\mathrm{Cov}(\mathbf{\Large{\epsilon}}) = \begin{pmatrix} \mathrm{Cov}(\epsilon_i, \epsilon_j) \end{pmatrix}_{1 \leq~i,~j~\leq N} = \begin{pmatrix} \sigma^2 & 0 & 0 & \cdots & 0 \\ 0 & \sigma^2 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \sigma^2 \end{pmatrix} = \sigma^2 I_{\tiny N}$

$\because \forall ~i \neq j \mathrm{: } \mathrm{Cov}(\epsilon_i, \epsilon_j) = 0$ and $\forall ~i = j \mathrm{~: Cov}(\epsilon_i, \epsilon_j) = \mathrm{Var}(\epsilon_i) = \sigma^2$

$\mathbf{\hat{w}}$ is our estimator for $\mathbf{w^{*}}$. We can write:

Substituting value of $\mathbf{Y}$ from equation $(17)$, we can write:

Then taking expectation of both sides of the equation $(22)$, we get:

We can see that our least-squared estimator $\mathbf{\hat{w}}$ is unbiased in probability with the true $\mathbf{w^{*}}$.

---

Lemma

Any affine transformation of a gaussian vector $\mathbf{z} \sim \mathcal{N} (\mathbf{0}, \Sigma)$ is also a gaussian vector:

Then, if $A$ is a matrix (that represents a linear transformation):

Then, if $A = \left( \Phi^T \Phi \right)^{-1} \Phi^T$:

---

Confidence Interval

To compute confidence interval on the weights i.e. on the $i^{th}$ component of $\mathbf{\hat{w}} = \begin{pmatrix} \hat{w_0} \\ \hat{w_1} \\ \vdots \\ \hat{w_d} \end{pmatrix} \sim \mathcal{N}( \mathbf{w^{*}}, \sigma^2 \left( \Phi^T \Phi \right)^{-1} )$ , we define:

$S = \left( \Phi^T \Phi \right)^{-1} = \begin{pmatrix} s_{00} & s_{11} & \cdots & s_{1d} \\ s_{10} & s_{11} & \cdots & s_{1d} \\ \vdots & \vdots & \ddots & \vdots \\ s_{d0} & s_{d1} & \cdots & s_{dd} \end{pmatrix}$

Then, we can write, $\text{Var}(\hat{w}) = \sigma^2 S$, which implies:

Now,

Then, we can write:

For example, say that we want to compute the $95\%$ confidence interval for the $i^{th}$ component of $\mathbf{\hat{w}}$, then $\alpha = 0.05$ and $z_{\alpha/2} = 1.96$:

Here, $\left[ \hat{w_i} \pm 1.96 \sigma \sqrt{s_{ii}} ~\right]$ is the $95\%$ confidence interval for the $i^{th}$ component of $\mathbf{\hat{w}}$.

---

Practical Problem

In practice, we don't know the true $\sigma^2$ and we have to estimate it from the data. We can estimate it using the following formula:

Here, $d$ is the number of features and $N$ is the number of samples.

The term $N - d - 1$ rather than $N$ in the denominator is used to make $\hat{\sigma}^2$ an unbiased estimate of $\sigma^2$.

Now, the empirical variance of the residual is nearly equal to $\sigma^2$. The term $\left( \frac{\hat{w_i} - w^{*}_i}{\hat{\sigma} \sqrt{s_{ii}}} \right)$ is now asymptotically gaussian with mean $0$ and variance $1$.

Notes

• If the value of $N$ is small, it follows a t-distribution (Student distribution). The exact Confidence Intervals can even be derived using Cochran's theorem.
• Our math for confidence intervals relies on the assumption that our data is homoscedastic. In the case that it is not, we can compute the confidence intervals empirically using the bootstrap method, which is based on resampling and the law of large numbers.

---

Bias-Variance Tradeoff

Regularization

Currently, we are using the following loss function:

We can also add a regularization term in the loss function, to prevent overfitting:

This is called Ridge or L2 Regularization. Here, $\lambda$ is the regularization parameter. It is a hyperparameter and can be tuned using cross-validation.

Note: We can also use other regularization methods like lasso (L1) regularization, square-root lasso, etc.

Now, to find a value of $\mathbf{w}$ that minimizes the loss function, we will solve:

Note: Now, even if $d > N$, $\left( \Phi^T \Phi + N \lambda \mathbb{I_d} \right)^{-1}$ is well defined because of the shift.

---

Prediction Error

Let's now look at the impact $\lambda$ on the generalization error.

Then,

Substituting this value in the above equation, we get:

Here:

• $\textcolor{#C060A1}{\langle N \lambda \left( \Phi^T \Phi + N \lambda \mathbb{I_d} \right)^{-1} \mathbf{w}^*, \mathbf{x_i} \rangle}$ is the non-zero bias term, which is deterministic. (It is $0$ if $\lambda = 0$.)
• $\textcolor{#03C988}{\langle \left( \Phi^T \Phi + N \lambda \mathbb{I_d} \right)^{-1} \Phi^T \epsilon, \mathbf{x_i} \rangle}$ is the variance term.
• $\textcolor{8B9A46}{\epsilon_i}$ is the noise term which is uncontrollable.

Here, if we take the expectation of $y_i - \hat{y_i}$, we get a non-zero term, which is the bias term. It's a pity because our non-regularized model had an unbiased least squared estimator (expectation was the ground truth). Therefore, by introducing a bias in our model, we are either underestimating or overestimating the ground truth.

The interest of introducing a bias in our model is that it enables us to control the variance term. Therefore, with a small $\lambda$, we have a high variance and a low bias. On the other hand, with a large $\lambda$, we have a low variance and a high bias. We call this tradeoff the bias-variance tradeoff.