Week3

Recap on Singular Value Decomposition

Orthogonal Matrix

Let $V\in {\mathbb{R}}^{p\times k}$, given by

$$V=\left[\begin{array}{ccc}\mid & & \mid \\ v_1& \cdots & v_k\\ \mid & & \mid \end{array}\right].$$

Then $V$ has orthogonal columns if

$$v_i^T{v}_j=\{\begin{array}{cc}& 1,\text{ if }i=j,\\ & 0,\text{ if }i\ne j.\end{array}$$

For such $V$, we have $V^{T}V=I_k$. Note that it is often NOT the case that $V{V}^T=I_p$.

A square matrix $V\in {\mathbb{R}}^{k\times k}$ with orthogonal columns is called an orthogonal matrix. Note that $V^{-1}=V^T$.

Definition Thin SVD

For a matrix $A\in {\mathbb{R}}^{p\times q}$ ($p\ge q$), the thin SVD of $A$ is a representation of the form

$$A=VD{U}^T,V\in {\mathbb{R}}^{p\times q},D\in {\mathbb{R}}^{q\times q},U\in {\mathbb{R}}^{q\times q},$$

whrere $V$ and $U$ have orthogonal columns and $D$ is diagonal.

The column vectors of $V$ and $U$ are respectively called the left and right singular vectors. The singular values of $A$ are the diagonal elements of $D$.

Remarkable results: Every matrix $A\in {\mathbb{R}}^{p\times q}$ has an SVD.

Spectral Decomposition

A spectral decomposition of a symmetric matrix $A\in {\mathbb{R}}^{p\times p}$ is a representation of $A$ as

$$A=VD{V}^T,V\in {\mathbb{R}}^{p\times p},D\in {\mathbb{R}}^{p\times p},$$

where $V=\left[\begin{array}{cccc}{v}_1& v_2& \dots & v_p\end{array}\right]$ is orthogonal and $D=diag\{d_1,d_2,\dots ,d_p\}$ is diagonal.

The columns of $V$ are the eigenvector of $A$ and $d_i$ is the associated eigenvalue:

$$A{v}_i=VD{V}^T{v}_i=d_i{v}_i.$$

Remarkable results (spectral theorem): Every symmetric matrix $A\in {\mathbb{R}}^{p\times p}$ has a spectral decomposition.

Recap on Multivariate Statistics

Notations

Let $\mathit{X}$ be a $k$-dimensional random vector (a special case of random matrices)

$$\mathit{X}=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_k\end{array}\right],$$

or a $p\times q$ random matrix

$$\mathit{X}=\left[\begin{array}{ccc}{X}_{11}& \cdots & X_{1q}\\ ⋮& \ddots & ⋮\\ X_{p1}& \cdots & X_{pq}\end{array}\right],$$

where each $X_i\in$ or $X_{ij}$ is a random variable. $E(\mathit{X})$ is defined componentwise, i.e.

$$E(\mathit{X})=\left[\begin{array}{ccc}E(X_{11})& \cdots & E(X_{1q})\\ ⋮& \ddots & ⋮\\ E(X_{p1})& \cdots & E(X_{pq})\end{array}\right],$$

and $E(\mathit{X})$ has linearity: for constant matrices $A\in {\mathbb{R}}^{d\times p}$, $B\in\mathbb{R}^{q\times q}, $C\in {\mathbb{R}}^{d\times q}$,

$$E(A\mathit{X}B+C)=AE(\mathit{X})B+C.$$

For random vectors $\mathit{X}\in {\mathbb{R}}^p$, its covariance matrix is defined as the $p\times p$ symmetric matrix $Cov(\mathit{X})$, given by

$$[Cov(\mathit{X}){]}_{ij}=Cov(X_i,X_j).$$

Similar to the covariance of random variable, we have

$$\begin{array}{rl}Cov(\mathit{X})& =E\left((\mathit{X}-E(\mathit{X})){(\mathit{X}-E(\mathit{X}))}^T\right)\\ & =E(\mathit{X}{\mathit{X}}^T)-E(\mathit{X})E(\mathit{X}{)}^T.\end{array}$$

From the above, we have

$$Cov(A\mathit{X})=A\cdot Cov(X)\cdot A^T.$$

Multivariate Normal Distribution

The multivariate normal distribution of a $k$-dimensional random vector $\mathit{X}={(X_1,\dots ,X_k)}^T$ can be written in the following notation:

$$\mathit{X}\sim \mathcal{N}(\mathit{\mu },\mathbf{\Sigma })\text{ or }\mathit{X}\sim {\mathcal{N}}_k(\mathit{\mu },\mathbf{\Sigma }),$$

with $k$-dimensional mean vector

$$\mathit{\mu }=E(\mathit{X})=\left[\begin{array}{c}E\left(X_1\right)\\ E\left(X_2\right)\\ ⋮\\ E\left(X_k\right)\end{array}\right]$$

and $k\times k$ covariance matrix

$${\mathrm{\Sigma }}_{i,j}=E\left((X_i-{\mu }_i)(X_j-{\mu }_j)\right)=Cov(X_i,X_j),\mathrm{\forall }i,j=1,\dots ,k.$$

$\mathbf{\Sigma }$ is assumed to be positive definite (i.e. non-degenerate) and therefore, ${\mathbf{\Sigma }}^{-1}$ is also positive definite. In this case, the density of $\mathit{X}$ is given by

$$p(\mathit{z})=\frac{1}{\sqrt{(2\pi {)}^{k}det(\mathbf{\Sigma })}}\exp (-\frac{1}{2}(\mathit{z}-\mathit{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathit{z}-\mathit{\mu })).$$

Fact: for a full-rank matrix $A\in {\mathbb{R}}^{p\times q}$ ($q\le p$) and $b\in {\mathbb{R}}^q$, we have

$$A\mathit{X}+b\sim {\mathcal{N}}_k(A\mathit{\mu }+b,A\mathbf{\Sigma }{A}^T).$$

Linear Regression

The basic idea of linear regression is to assume that

$$y\approx {\beta }_0+\sum _{i=1}^p{\beta }_i{x}_i,$$

where $\mathit{X}=[x_1,x_2,\dots ,x_p{]}^T$ (for now we assume $\mathit{x}\in {\mathbb{R}}^p$). More precisely, we make a assumption about $p(y\mid \mathit{x})$ as follows:

$$y={\beta }_0+\sum _{i=1}^p{\beta }_i{x}_i+\epsilon ,$$

where $\epsilon$ is noise with $\epsilon \sim \mathcal{N}(0,{\sigma }^2)$ and $\mathit{x}$ is independent of the noise. Thus, we have

$$y\mid \mathit{x}\sim \mathcal{N}({\beta }_0+\sum _{i=1}^p{\beta }_i{x}_i,{\sigma }^2).$$

We know that (under squared error loss) the oracle predictor is

$$E(y\mid \mathit{x})={\beta }_0+\sum _{i=1}^p{\beta }_i{x}_i=f^{\ast }(\mathit{x}).$$

The goal is to find ${\beta }_0,\dots ,{\beta }_p$ and thus $f^{\ast }$.

Rewrite Training Data in Matrix Form

Assume training data $\{({\mathit{x}}^{(i)},y^{(i)}){\}}_{i=1}^n$ were i.i.d. generated via the assumption about $p(y\mid \mathit{x})$. To simplify the notations, we define

$$\mathit{Y}=\left[\begin{array}{c}{y}^{(1)}\\ ⋮\\ y^{(n)}\end{array}\right],\mathit{\epsilon }=\left[\begin{array}{c}{\epsilon }^{(1)}\\ ⋮\\ {\epsilon }^{(n)}\end{array}\right],\mathit{\beta }=\left[\begin{array}{c}{\beta }^{(0)}\\ ⋮\\ {\beta }^{(n)}\end{array}\right],$$

and design matrix $\mathit{X}\in {\mathbb{R}}^{n\times (p+1)}$ given by

$$\mathit{X}=\left[\begin{array}{cc}1& {\left({\mathit{x}}^{(1)}\right)}^T\\ 1& {\left({\mathit{x}}^{(2)}\right)}^T\\ ⋮\\ 1& {\left({\mathit{x}}^{(n)}\right)}^T\end{array}\right].$$

Then training data can written as

$$\mathit{Y}=\mathit{X}\mathit{\beta }+\mathit{\epsilon }.$$

Now we want to use training data to estimate ${\beta }_0,\dots ,{\beta }_p$ and thus $f^{\ast }$.

Recap on Likelihood Function

The likelihood function $L(\theta \mid \mathit{x})$ illustrates the probability of $\theta$ given data set $\mathit{x}=(x_1,\dots ,x_n{)}^T$.

More precisely, assume we have $n$ samples $(x_1,\dots ,x_n)$ observed from a distribution $p_{\theta }(x)$ with an unknown parameter $\theta$, and our goal is to estimate $\theta$ using the observed data.

We view the observed samples are realizations of some random variables $X_1,X_2,\dots ,X_n$, which has a joint density function $p(X_1,\dots ,X_n\mid \theta )$. Given $X_1=x_1,\dots ,X_n=x_n$, we may consider the probability of this event being observed, which is the likelihood function (a function of $\theta$) defined by:

$$L(\theta )=L(\theta \mid x_1,\dots ,x_n)=p(X_1=x_1,\dots ,X_n=x_n\mid \theta ).$$

Note that the likelihood function is NOT a probability density function. It measures the support provided by the data for each possible value of the parameter. If we compare the likelihood function at two parameter points and find that

$$L({\theta }_1\mid \mathit{x})>L({\theta }_2\mid \mathit{x})$$

then the sample we actually observed is more likely to have occurred if $\theta ={\theta }_1$ than if $\theta ={\theta }_2$. This can be interpreted as ${\theta }_1$ is a more plausible value for $\theta$ than ${\theta }_2$. Therefore, one approach to estimate $\theta$ is to choose the value of $\theta$ which gives you the highest probability among all possible values.

With the assumption that samples are drawn i.i.d., the joint probability is given by multiplied probabilities, i.e.

$$p(X_1,\dots ,X_n\mid \theta )=\prod _{i=1}^n{p}_{\theta }\left(X_i\right),$$

hence taking the log transforms this into a summation, which is usually easier to maximize analytically. Thus, we often write the log-likelihood rather than the likelihood.

Estimate the Coefficients ${\beta }_i$

Here we find an estimate of $\mathit{\beta }$ by maximum likelihood estimation:

By $\mathit{X}$ is independent of the noise (and vice versa), we have

$$\mathit{Y}\mid \mathit{X}\sim {\mathcal{N}}_n(\mathit{X}\mathit{\beta },{\sigma }^2{I}_n),$$

and thus, the (log-)likelihood function is given by the density function of the above multivariate normal distribution.

By maximizing the likelihood function, we can estimate $\mathit{\beta }$:

$$\begin{array}{rl}\hat{\mathit{\beta }}& =\underset{\beta }{\arg \max }\ln (p(\mathit{Y}\mid \mathit{X}))\\ & =\underset{\beta }{\arg \max }(\frac{1}{2{\sigma }^2}(\mathit{Y}-\mathit{X}\mathit{\beta }{)}^T(\mathit{Y}-\mathit{X}\mathit{\beta }))\\ & =\underset{\beta }{\arg \max }\Vert \mathit{Y}-\mathit{X}\mathit{\beta }{\Vert }^2,\end{array}$$

which becomes a least squares problem. We take a gradient with respect to $\mathit{\beta }$ and set the gradient for 0, and after solving for $\mathit{\beta }$, we have

$$\hat{\mathit{\beta }}={\left({\mathit{X}}^T\mathit{X}\right)}^{-1}{\mathit{X}}^T\mathit{Y},$$

where we assume $\mathit{X}$ is a full-rank matrix and thus ${\mathit{X}}^T\mathit{X}$ has an inverse

.

Furthermore, we can find the distribution of $\hat{\mathit{\beta }}$ conditional on the training features $\mathit{X}$:

$$E(\hat{\mathit{\beta }}\mid \mathit{X})=E({\left({\mathit{X}}^T\mathit{X}\right)}^{-1}{\mathit{X}}^T\mathit{Y}\mid \mathit{X})=\mathit{\beta },$$

which indicates that $\hat{\mathit{\beta }}$ is unbiased;

$$\begin{array}{rl}Cov(\hat{\mathit{\beta }}\mid \mathit{X})& =Cov({\left({\mathit{X}}^T\mathit{X}\right)}^{-1}{\mathit{X}}^T\mathit{Y}\mid \mathit{X})\\ & ={\left({\mathit{X}}^T\mathit{X}\right)}^{-1}{\mathit{X}}^T\cdot Cov(\mathit{Y}\mid \mathit{X})\cdot \mathit{X}{\left({\left({\mathit{X}}^T\mathit{X}\right)}^{-1}\right)}^T\\ & ={\sigma }^2({\mathit{X}}^T\mathit{X}{)}^{-1};\end{array}$$

and therefore,

$$\hat{\mathit{\beta }}\mid \mathit{X}\sim {\mathcal{N}}_{p+1}(\mathit{\beta },{\sigma }^2({\mathit{X}}^T\mathit{X}{)}^{-1}).$$

Evaluate OLS Prediction Estimate

The above estimation is called ordinary least squares (OLS) in statistics. Now we evaluate our estimation when receiving a new data pair $({\mathit{x}}_{\mathbf{\ast }},y_{\ast })$. We make our prediction by

$$\hat{f}({\mathit{x}}_{\mathbf{\ast }})={\tilde{\mathit{x}}}^T\hat{\mathit{\beta }},\tilde{\mathit{x}}=\left[\begin{array}{c}1\\ {\mathit{x}}_{\mathbf{\ast }}\end{array}\right].$$

Consider the bias of $\hat{f}$ (conditional on $\mathit{X},{\mathit{x}}_{\mathbf{\ast }}$),

$$E(\hat{f}({\mathit{x}}_{\mathbf{\ast }})\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})={\tilde{\mathit{x}}}^{T}E(\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})={\tilde{\mathit{x}}}^T\mathit{\beta }.$$

We know that $E(y_{\ast }\mid {\mathit{x}}_{\mathbf{\ast }})={\tilde{\mathit{x}}}^T\mathit{\beta }$, and thus the OLS prediction is unbiased

.

Now, we consider the variance:

$$Var(\hat{f}({\mathit{x}}_{\mathbf{\ast }})\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})={\tilde{\mathit{x}}}^T\cdot Var(\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})\tilde{\mathit{x}}={\sigma }^2{\tilde{\mathit{x}}}^T({\mathit{X}}^T\mathit{X}{)}^{-1}\tilde{\mathit{x}}.$$

Roughly speaking, if ${\mathit{x}}_{\mathbf{\ast }}$ is similar to the training data features, the variance will be small; otherwise, the variance will be large.

Interval Estimate

We want to make a prediction interval for $y_{\ast }$ in the form of

$$[C_{low}(\mathit{Y},\mathit{X},{\mathit{x}}_{\mathbf{\ast }}),C_{high}(\mathit{Y},\mathit{X},{\mathit{x}}_{\mathbf{\ast }})].$$

We hope $y_{\ast }$ to be contained in this with probability of at least $1-\alpha$ (usually $\alpha =0.05$ or $\alpha =0.01$), i.e.

$$Pr(y_{\ast }\in [C_{low},C_{high}]\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})\ge 1-\alpha .$$

Recall the assumption that $y_{\ast }\mid \stackrel{\mathbf{~}}{\mathit{x}}\sim \mathcal{N}({\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }},{\sigma }^2)$. By

$$E(y_{\ast }-{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})=0,$$

and $Cov(y_{\ast },{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})=0$

and thus

$$\begin{array}{rl}& Var(y_{\ast }-{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})\\ & =Var(y_{\ast }\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})-2Cov(y_{\ast },{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})+Var({\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})\\ & ={\sigma }^2+{\sigma }^2{\stackrel{\mathbf{~}}{\mathit{x}}}^T({\mathit{X}}^T\mathit{X}{)}^{-1}\tilde{\mathit{x}},\end{array}$$

we know that $y_{\ast }-{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}$ is still normally distributed, given by

$$y_{\ast }-{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }}\sim \mathcal{N}(0,{\sigma }^2(1+{\stackrel{\mathbf{~}}{\mathit{x}}}^T({\mathit{X}}^T\mathit{X}{)}^{-1}\tilde{\mathit{x}})).$$

To obtain the interval estimation, we also need to get rid of ${\sigma }^2$ as it is unknown. Recall the assumption that the noise is normally distributed. We can estimate the variance ${\sigma }^2$ by the scaled residuals, i.e.

$${\hat{\sigma }}^2=\frac{\Vert \mathit{Y}-\mathit{X}\widehat{\mathit{\beta }{\Vert }^2}}{n-p}.$$

Then we have

$$\frac{y_{\ast }-{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}}{\sqrt{{\hat{\sigma }}^2(1+{\stackrel{\mathbf{~}}{\mathit{x}}}^T({\mathit{X}}^T\mathit{X}{)}^{-1}\tilde{\mathit{x}})}}\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }}$$

follows a t-distribution with $n-p$ degree of freedom. Let $t_{n-p,\alpha /2}$ and $t_{n-p,1-\alpha /2}$ be the $\alpha /2$ and $1-\alpha /2$ quantile of the this distribution. Denote $\sqrt{{\hat{\sigma }}^2(1+{\stackrel{\mathbf{~}}{\mathit{x}}}^T({\mathit{X}}^T\mathit{X}{)}^{-1}\tilde{\mathit{x}})}$ as $C$. Then, we have

$$Pr(\frac{y_{\ast }-{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}}{C}\in [t_{n-p,\alpha /2},t_{n-p,1-\alpha /2}]\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})\ge 1-\alpha ,$$

that is

$$Pr(y_{\ast }\in [C{t}_{n-p,\alpha /2}+{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }},C{t}_{n-p,1-\alpha /2}+{\stackrel{\mathbf{~}}{\mathit{x}}}^T\hat{\mathit{\beta }}]\mid \mathit{X},{\mathit{x}}_{\mathbf{\ast }})\ge 1-\alpha .$$

Examples and Discrete Features

So far we assume $\mathit{x}$ lies in ${\mathbb{R}}^p$. Here we will introduce discrete features and mixed features by giving some examples.

Example Polynomial Regression ($\mathcal{X}=\mathbb{R}$):

The model takes $x$ and transforms it to a new feature

$$x\mapsto (1,x,x^2,\dots ,x^d),$$

and the model becomes

$$y={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+\cdots +{\beta }_d{x}^d+\epsilon ,$$

use this new feature in a regression. For instance, assume we have 4 data pairs for training and use a 2-degree polynomial in regression: training data is given by

$$\left[\begin{array}{c}{x}^{(1)}\\ x^{(2)}\\ x^{(3)}\\ x^{(4)}\end{array}\right]=\left[\begin{array}{c}2\\ -1\\ 3\\ 5\end{array}\right].$$

We first map the following training data to the design matrix $\mathit{X}$, given by

$$\left[\begin{array}{ccc}1& 2& 4\\ 1& -1& 1\\ 1& 3& 9\\ 1& 5& 25\end{array}\right].$$

The OLS estimate is still $\hat{\mathit{\beta }}={\left({\mathit{X}}^T\mathit{X}\right)}^{-1}{\mathit{X}}^T\mathit{Y}$. This allows us to fit $y$ wit ha non-linear function of $x$.

Example Categorical Feature:

A categorical feature could be, for instance 3 categories as

$$x=\{\begin{array}{cc}1,& \text{ if Age}<20,\\ 2,& \text{ if }20\le \text{Age}<40,\\ 3,& \text{ if Age}\ge 40.\end{array}$$

We transform categorical $x$'s via indicator functions.

First way is the dummy variable method: By choosing a baseline class ( say class '3'), we transform

$$x\mapsto (1,I(x=1),I(x=2)),$$

and the regression model becomes

$$y={\beta }_0+{\beta }_{1}I(x=1)+{\beta }_{2}I(x=2)+\epsilon ,$$

where we have $E(y\mid x=3)={\beta }_0$, $E(y\mid x=1)={\beta }_0+{\beta }_1$, and $E(y\mid x=2)={\beta }_0+{\beta }_2$.

Another way is more symetric -- we have the model

$$y={\beta }_0+{\beta }_{1}I(x=1)+{\beta }_{2}I(x=2)+{\beta }_{3}I(x=3)+\epsilon ,$$

where we have $E(y\mid x=1)={\beta }_0+{\beta }_1$, $E(y\mid x=2)={\beta }_0+{\beta }_2$, and $E(y\mid x=3)={\beta }_0+{\beta }_3$. Since ${\beta }_0$ is redundant, it can be dropped.

Assume we have training data:

$$\left[\begin{array}{c}{x}^{(1)}\\ x^{(2)}\\ x^{(3)}\\ x^{(4)}\end{array}\right]=\left[\begin{array}{c}1\\ 2\\ 3\\ 1\end{array}\right],$$

and the corresponding design matrix is

$$\mathit{X}=\left[\begin{array}{ccc}1& 1& 0\\ 1& 0& 1\\ 1& 0& 0\\ 1& 1& 0\end{array}\right].$$

Example Mixed Features

We can combine continous and categorical features. For instance, if $x\in \{1,2,3\}\times \mathbb{R}$, we can do similar mapping:

$$x\mapsto (1,I(x_1=1),I(x_1=2),x_2,x_2^2).$$

Assume the training data is

$$\left[\begin{array}{c}{x}^{(1)}\\ x^{(2)}\\ x^{(3)}\end{array}\right]=\left[\begin{array}{cc}1& 3.5\\ 2& 7\\ 3& -2.1\end{array}\right],$$

then the design matrix is given by

$$\mathit{X}=\left[\begin{array}{ccccc}1& 1& 0& 3.5& 12.25\\ 1& 0& 1& 7& 49\\ 1& 1& 0& -2.1& 4.41\end{array}\right].$$

Example Intersactions between Features

We can also have intersactions between features. For instance, if $\mathcal{X}\in {\mathbb{R}}^2$, then we can consider thier multiplication in the mapping:

$$x\mapsto (1,x_1,x_2,x_1{x}_2,x_1^2,x_2^2).$$

Remark Why don't we just use lots of features? I.e. in polynomial regression, why don't we take degree $d$ to be large? That is because problems owith overfitting! This is mentioned in the variance-bias tradeoff.

Variable Selection

Assume we have many features, i.e. $\mathit{X}$ has many columns. Given an index set $I=(i_1,i_2,\dots ,i_k)$, where $i_1\le i_2\le \cdots \le i_k$, define

$${\mathit{X}}_I=\left[\begin{array}{ccc}\mid & & \mid \\ x_{i_1}& \cdots & x_{i_k}\\ \mid & & \mid \end{array}\right],{\mathit{\beta }}_I=\left[\begin{array}{c}{\beta }_{i_1}\\ ⋮\\ {\beta }_{i_k}\end{array}\right].$$

Then we can obtain a submodel with features in $I$ give nby

$$\mathit{Y}={\mathit{X}}_I{\mathit{\beta }}_I+\epsilon .$$

Because of the problem of overfitting, we may not want to use empirical risks to evaluate the submodel (the empirical risk will always decrease when adding new features). A better idea is to add a penalty for including new features:

$$\begin{array}{rl}& \mathrm{AIC}(I)=\frac{1}{n}(\Vert \mathit{Y}-{\mathit{X}}_I{\hat{\mathit{\beta }}}_I{\Vert }^2+2k{\hat{\sigma }}^2),\\ & \mathrm{BIC}(I)=\frac{1}{n}(\Vert \mathit{Y}-{\mathit{X}}_I{\hat{\mathit{\beta }}}_I{\Vert }^2+\ln (n)k{\hat{\sigma }}^2),\end{array}$$

where if $\hat{\mathit{\beta }}$ is estimated by using the full model, $\hat{\sigma }$ is given by

$${\hat{\sigma }}^2=\frac{\Vert \mathit{Y}-{\mathit{X}}_I\hat{\mathit{\beta }}{\Vert }^2}{n-p}.$$

We want to choose proper index set $I$ to maximize AIC or BIC. For simplicity, we may view $I$ as a string in length $p$ consists of 0 and 1, where 0 indicates the feature is not included and 1 indicates the feature is included: for instance,

$$(0,1,0,0,1)$$

is an index set of 5 features where the second and fifth feature are included. In this way, $I$ is a point in $\{0,1{\}}^p$ and we may visualize the choice of $I$ as picking a vertex of a hypercube. Minimizing AIC is the same as searching for the optimal vertex on the cube.

For $p>40$, search becomes infeasible. If we enumerate the AIC for every $I$, this called best subset selection. Instead, we can do greedy search over the cube.

Forward selection:

1. Start with trivial null model $y={\beta }_0+\epsilon$.

2. Add a feature by searching for the model $y={\beta }_0+{\beta }_i{x}_i+\epsilon$ over $i$ that has the smallest AIC/BIC. Call this best $i$ to be $i_1$.

3. Continue by find an $j$ where the model $y={\beta }_0+{\beta }_{i_1}{x}_{i_1}+{\beta }_j{x}_j+\epsilon$ maximizes AIC/BIC. Call this $j$ to be $i_2$.

4. Keep on doing this until we reach the full model.

5. Choose the model in our sequence with the smallest AIC/BIC.

Backward selection: start at full model and remove variables until reaching the null model.

Best subset selection needs to fit $2^p$ linear models. On the other hand, for forward (backward) selection needs to fit

$$p+(p-1)+\cdots +1=\frac{p(p+1)}{2}$$

models, which is an order of $O(p^2)$ computations.