Linear Regression

Linear regression adopts a linear-in-parameter regression function model $\{f_w{\}}_{w\in {\mathbb{R}}^n}$ such that

$$\hat{y}=f_w(x)=w^{T}x.$$

In this note, we focus on real-valued output $y\in \mathbb{R}$. We can think of $w$ as a set of weights that determine how each feature affects the prediction.

Linear regression is the correct model for data generated by a linear model:

$$Y_i=w^T{X}_i+{ϵ}_i,$$

where ${ϵ}_i$ is a zero-mean random noise. If the noise follows a Gauss distribution, we get a Gaussian Linear Model, whose estimator of $w$ have nice properties.

Suppose we have collected $n$ sample points $(X_i,Y_i)$, the matrix form of the model is

$$Y=Xw+ϵ,$$

where $X=(X_1,\dots ,X_n{)}^T\in {\mathbb{R}}^{n\times d}$ is called the design matrix, $Y=(Y_1,\dots ,Y_n{)}^T\in {\mathbb{R}}^n$ is the output vector, and $ϵ=({ϵ}_1,\dots ,{ϵ}_n{)}^T\in {\mathbb{R}}^n$ is the noise vector. Regression focuses on fitting the collected (training) data, corresponding to a fixed design problem; Prediction focuses on predicting the new (test) data, corresponding to a random design problem.

Ordinary Least Squares

We use Mean Squared Error as the performance measure

$${\mathrm{MSE}}_{train/test}=\frac{1}{m}\sum _i({\hat{y}}_i^{(train/test)}-y_i^{(train/test)}{)}^2=\frac{1}{m}\Vert {\hat{\mathbf{y}}}^{(train/test)}-{\mathbf{y}}^{(train/test)}{\Vert }_2^2$$

Ordinary least squares linear regression with a fixed design minimizes the $\mathrm{MS}{\mathrm{E}}_{\mathrm{train}}$. From now on we drop the $(\mathrm{train})$ superscript. The solution is just the least squares estimator (LSE) of $w$ (M-Estimator with Risk being MSE). Since the MSE is convex, we have

$$\begin{array}{rl}& {\nabla }_w{\mathrm{MSE}}_{\mathrm{train}}=0\\ \Rightarrow & {\nabla }_w\Vert Xw-y{\Vert }_2^2=0\\ \Rightarrow & {\nabla }_w{\left(Xw-y\right)}^T\left(Xw-y\right)=0\\ \Rightarrow & 2X^{T}Xw-2X^{T}y=0\\ \Rightarrow & w={\left(X^{T}X\right)}^{-1}{X}^{T}y\end{array}$$

The last set of equations is known as normal equations. Clearly, the least squares estimator requires $X$ to be full column rank, meaning that there are more samples than features. When there are more features than samples, which is common in the high-dimensional setting, we have an Underdetermined Linear System and there are infinitely many solutions with the lest square error, and we usually seek the least norm solution.

• 💡 The geometric perception of the solution is $X{w}_{\mathrm{LS}}={\mathrm{proj}}_{\mathrm{span}X}y$. Therefore, $w_{\mathrm{LS}}$ achieves the optimality. ^69df70

Misspecified Model

Consider a general non-parametric statistical model, $(X_i,Y_i)\stackrel{\text{iid}}{\sim }P\in \mathcal{P}\{P:P\in \mathrm{dist}({\mathbb{R}}^{d+1})\text{ w/ finite fourth moment}\}$ . The least square solution is targeting $w^{\ast }={\mathrm{argmin}}_w\mathbb{E}[(Y_i-X_i^{\top }w{)}^2]$, i.e., best linear model of $Y_i$ given $Z_i$. By LLN and CMT, ${\hat{w}}^{\mathrm{LS}}\stackrel{P}{\to }{w}^{\ast }$. Further ${\hat{w}}^{\mathrm{LS}}$ is asymptotically normal with some asymptotic variance $\Sigma$, which consists of the sample noise and misspecification noise.

More specifically, let’s look at the expression of $w^{\ast }$. Similar to the derivation of the normal equation, we have

$$\begin{array}{rl}& {\nabla }_w\mathbb{E}[(Y-X^{T}w{)}^2]=0\\ ⟹& {\nabla }_w\left(\mathbb{E}{Y}^2-2w^T\mathbb{E}XY+w^T(\mathbb{E}X{X}^T)w\right)=0\\ ⟹& w^{\ast }={\left(\mathbb{E}X{X}^T\right)}^{-1}\mathbb{E}XY.\end{array}$$

To incorporate the intercept, we can add a constant $1$ to the first dimension of each observation $X^{\prime }=(1;X)$. Then $w^{\ast }=(a^{\ast },b^{\ast })$, where $a^{\ast }\in \mathbb{R}$ is the intercept and $b^{\ast }\in {\mathbb{R}}^d$ is the slope. And using the Schur Complement, we have

$$\left(\begin{array}{c}{a}^{\ast }\\ b^{\ast }\end{array}\right)=\left(\begin{array}{cc}1+\mathbb{E}{X}^T\mathrm{Cov}(X{)}^{-1}\mathbb{E}X& -\mathbb{E}{X}^T\mathrm{Cov}(X{)}^{-1}\\ -\mathrm{Cov}(X{)}^{-1}\mathbb{E}X& \mathrm{Cov}(X{)}^{-1}\end{array}\right)\left(\begin{array}{c}\mathbb{E}Y\\ \mathbb{E}XY\end{array}\right)=\left(\begin{array}{c}\mathbb{E}Y+\mathbb{E}{X}^T\mathrm{Cov}(X{)}^{-1}\mathrm{Cov}(X,Y)\\ \mathrm{Cov}(X{)}^{-1}\mathrm{Cov}(X,Y)\end{array}\right).$$

We can see that the slope measures the correlation between $Z$ and $Y$.

Let’s define the residual ${ϵ}_i=Y_i-X_i^T{w}^{\ast }$. Then, with the additional constant component, we have $\mathbb{E}{ϵ}_i=0$ and

$$\mathrm{Cov}(X_i,{ϵ}_i)=\mathrm{Cov}(X_i,Y_i)-\mathrm{Cov}(X_i)w^{\ast }=0.$$

Therefore, we can see that although the noise $ϵ$ now consists of both the sample noise and the misspecification noise, it is still uncorrelated with the input $X$ and has zero mean.

• ❗️ Importantly, zero correlation does not imply independence. Clearly, for a misspecified model, $Y-X^T{w}^{\ast }|X=x$ can have a larger variance and bias for some $x$ than others, and the noise $ϵ$ is not independent of $X$. The intercept balances out the overall bias and the slope turns perpendicular to the residual, which is consistent with the geometric perception of LSE.

Bias

We can see that it is more natural and general to use affine functions instead of linear functions as regression functions:

$$\hat{y}=w^{T}x+b$$

where $b\in \mathbb{R}$, is called the bias parameter or the intercept term, for the prediction is biased toward being $b$ in the absence of any input. This term is important to achieve an zero-mean residual (Misspecified Model).

Attaching a constant 1 to the first dimension of each observation $x_i$ is a simple way to adopt this bias term: $\hat{y}=(b;w{)}^T(1;x)$. Thus, linear regression encompasses affine functions.

A Probabilistic View: Maximum Likelihood Estimation

Suppose our residual noise follows a Gaussian distribution $ϵ\sim \mathcal{N}(0,{\sigma }^{2}I)$, then our statistical model is $y\mid x\sim \mathcal{N}(w^{T}x,{\sigma }^2)$ parameterized by $w$. Or, given a fixed design, $Y\sim \mathcal{N}(Xw,{\sigma }^{2}I)$.

Then, the Maximum Likelihood Estimation of $w$ is

$$\begin{array}{rl}w& =\underset{w}{\mathrm{argmax}}\sum _i\log P(Y_i|X_i;w)\\ & =\underset{w}{\mathrm{argmax}}\left(-m\log \sigma -\frac{m}{2}\log (2\pi )-\sum _i\frac{\Vert {\hat{Y}}_i-Y_i{\Vert }^2}{2{\sigma }^2}\right)\\ & =\underset{w}{\mathrm{argmin}}\sum _i\Vert {\hat{Y}}_i-Y_i{\Vert }^2\\ & =\underset{w}{\mathrm{argmin}}\Vert Y-Xw{\Vert }^2\\ & =\underset{w}{\mathrm{argmin}}\mathrm{MSE}\end{array}$$

Therefore, if the data-generating model is a Gaussian Linear Model, i.e., the noise is an independent/multivariate Gaussian noise, LSE is equivalent to MLE. Further, under this assumption, $w_{\mathrm{MLE}}$ is an unbiased estimate of $w$:

$$\mathbb{E}[w_{\mathrm{MLE}}]=\mathbb{E}[(X^{T}X{)}^{-1}{X}^{T}y]=(X^{T}X{)}^{-1}X\mathbb{E}[y]=(X^{T}X{)}^{-1}{X}^{T}Xw=w.$$

The variance of MLE is

$$\begin{array}{rl}\mathrm{Var}[w_{\mathrm{MLE}}]=& \mathbb{E}[w_{\mathrm{MLE}}{w}_{\mathrm{MLE}}^T]-\mathbb{E}[w_{\mathrm{MLE}}]\mathbb{E}[w_{\mathrm{MLE}}{]}^T\\ =& \mathbb{E}\left[{\left(X^{T}X\right)}^{-1}{X}^{T}Y{Y}^{T}X{\left(X^{T}X\right)}^{-1}\right]-w{w}^T\\ =& {\left(X^{T}X\right)}^{-1}{X}^T\mathbb{E}\left[Y{Y}^T\right]X{\left(X^{T}X\right)}^{-1}-w{w}^T\\ =& {\left(X^{T}X\right)}^{-1}{X}^T(\mathrm{Var}(Y)+\mathbb{E}[Y]\mathbb{E}[Y{]}^T)X{\left(X^{T}X\right)}^{-1}-w{w}^T\\ =& {\left(X^{T}X\right)}^{-1}{X}^T\left({\sigma }^{2}I+Xw{w}^T{X}^T\right)X{\left(X^{T}X\right)}^{-1}-w{w}^T\\ =& {\left(X^{T}X\right)}^{-1}{X}^T{\sigma }^{2}IX{\left(X^{T}X\right)}^{-1}+\\ & {\left(X^{T}X\right)}^{-1}{X}^{T}Xw{w}^T{X}^{T}X{\left(X^{T}X\right)}^{-1}-w{w}^T\\ =& {\sigma }^2{\left(X^{T}X\right)}^{-1}\\ =& {\sigma }^{2}V{\Sigma }^{-2}{V}^T,\end{array}$$

where $X=U\Sigma V^T$ is the Singular Value Decomposition of $X$. And thus the final MSE is

$$\mathrm{MSE}(w_{\mathrm{MLE}})=\mathbb{E}\Vert w_{\mathrm{MLE}}-w{\Vert }^2=\mathbb{E}\mathrm{tr}((w_{\mathrm{MLE}}-w)(w_{\mathrm{MLE}}-w{)}^T))=\mathrm{tr}(\mathrm{Var}(w_{\mathrm{MLE}}))={\sigma }^2\mathrm{tr}({\left(X^{T}X\right)}^{-1}).$$

• ❗️ When the data are highly correlated, there will be small singular values, the values of $w_{\mathrm{MLE}}$ are very sensitive to the measured data $y$ and give unstable predictions for new data. This is bad if we want to analyze and predict using $w_{\mathrm{MLE}}$.

Ridge Regression can help with this problem.

Why "squares"?

The nice properties of least squares estimation and its connection to Gaussian distribution are not coincidences. The hidden player is the Euclidean geometry. Recall that an alternative definition of the Gaussian distribution is that it’s some affine transformation of some rotation-invariant distribution; and the geometric perception of LSE is that it is doing orthogonal projection. Here, both rotation and projection are with respect to the Euclidean geometry. More specifically, both Gaussian distributions and the ${\ell }_2$ norm in MSE are defined using the natural Inner Product of a self-dual vector space, which is the Euclidean space, giving natural quadratic (bilinear) forms.

Polynomial Regression

We can also use a polynomial function of degree n to fit the 2-D data. Then our estimation becomes

$$\hat{y}=w_n{x}^n+\cdots +w_{1}x+b$$

• ❗️ In polynomial regression, $n$ is often a Hyperparameter

The polynomial regression is the generalization of linear regression with bias parameter; however, it can be turned into a linear regression:

$$\hat{y}={\tilde{w}}^T\tilde{x}\triangleq (b,w_1,\dots ,w_n)(x^0,x^1,\dots ,x^n{)}^T$$

• 💡 Thus we can solve the polynomial regression in the same way as solving the linear regression. Actually, by the definition of linear regression, polynomial regression is still linear.

High-Dimensional Data

Polynomial regression can also apply to high-dimensional data. Similarly, we can extend an observation $x=(x_1,\dots ,x_d)$ to

$$\tilde{x}=(x_1,\dots ,x_d,x_1^2,\dots ,x_d^2,\dots ).$$

We can also add features like $x_1{x}_2$.

Generalized Linear Regression

By the definition of linear regression, the general form of linear regression is

$$y_i\approx f\left(x_i,w\right)=\sum _{s=1}^S{g}_s\left(x_i\right)w_s=w^{T}G(x_i).$$

For example,

$$\begin{array}{rl}& g_s\left(x_i\right)=x_{ij}^2\\ & g_s\left(x_i\right)=\log x_{ij}\\ & g_s\left(x_i\right)=\mathbb{I}\left(x_{ij}<a\right)\\ & g_s\left(x_i\right)=\mathbb{I}\left(x_{ij}<x_{i{j}^{\prime }}\right)\end{array}$$

As long as the function is linear in $w$, we can construct the matrix $X$ by putting the transformed $x_i$ on row $i$, and solve $w_{\mathrm{LS}}=(X^{T}X{)}^{-1}{X}^{T}y$.

• ❗️ As the number of functions increases, we need more data to avoid overfitting.

The indicator functions are called dummy variables, and the functions that depend on more than one factor (like the last example) capture the interaction effects. The whole model is called a generalized linear regression model, and $\{g_s{\}}_{s=1}^S$ are called basis functions.

Choosing indicator functions as basis functions gives step function regression. Another example of a generalized linear regression model is Splines.

Dummy Variables for Qualitative Predictors

We can use dummy variables (indicator functions) to represent qualitative predictors. One problem is, if we introduce a dummy variable for each level of a qualitative predictor, the matrix $X$ will be singular. To avoid this, we can introduce $k-1$ dummy variables for a qualitative predictor with $k$ levels. The left level can be regarded as the control/baseline/reference category. Then, the slope of each dummy variable represents the difference between the corresponding level and the baseline category. And the difference between different slopes represents the difference between different levels.