4.1 Regression

4.1.1 Fundamentals

In general, we can define the “loss” of guessing \(\hat {y}\) when the true value was is \(y\) as \[ \mathscr {L}(\hat {y}, y) = \mathscr {L}(f({\mathbf{x}} ), y) \] We want to find \(f(\cdot )\) so that \(\mathscr {L}(f({\mathbf{x}} ), y)\) is as small as possible. But the loss for a particular \(f\) may depend on \({\mathbf{x}} _{\mathrm {new}}\), since different \({\mathbf{x}} _{\mathrm {new}}\) are associated with different distributions of \(y_{\mathrm {new}}\). It follows that it doesn’t make sense to minimize the loss — instead we want to minimize the risk, \[ \mathscr {L}(f) := \mathbb {E}\left [\mathscr {L}(f({\mathbf{x}} ), y)\right ] \qquad \text {(Risk)}. \] So we would like to find \[ f^{\star }(\cdot ) := \operatorname *{argmin}_{f} \mathscr {L}(f) = \operatorname *{argmin}_{f} \mathbb {E}\left [\mathscr {L}(f({\mathbf{x}} ), y)\right ]. \]

We’ll consider square loss first: \[ l({\mathbf{x}} , y) = \frac 1 2(y - f({\mathbf{x}} ))^2 \quad \hat {\mathcal {L}}(f) = \frac {1}{N} \sum _{i=1}^N l({\mathbf{x}} _i,y_i). \] We cannot actually compute \(f^*\), a general question is how to calculate the gap between \(f^*\) and \(\hat f\). \[ 0 \leq \mathscr {L}(\hat {f}) - \mathscr {L}(f^*) = \underbrace {\mathscr {L}(\hat {f}) - \hat {\mathscr {L}}(\hat {f})}_{\text {Difficult term!}} + \underbrace {\hat {\mathscr {L}}(\hat {f}) - \hat {\mathscr {L}}(f^*)}_{\text {Negative term}} + \underbrace {\hat {\mathscr {L}}(f^*) - \mathscr {L}(f^*)}_{\text {LLN term}}. \] It is called generalization error and the term \(\mathcal {L}(\hat f) - \hat {\mathcal {L}}(\hat f)\) can be bounded \(\sup _{f} [\mathcal {L}(f) - \hat {\mathcal {L}}(f)]\), which will be discussed in the later sections.

We begin with linear regression model: \[ y_n = \beta _1 x_{n1} + \beta _2 x_{n2} + \ldots + \beta _P x_{nP} + \varepsilon _n = \underbrace {{\mathbf{x}} _n^\top {\mathbf{\beta}} }_{f({\mathbf{x}} _n)} + \epsilon _n, \quad \text {For } n = 1, \ldots , N. \] where we use \(N\) to denote the number of data and \(P\) to denote the number of features. To get the optimal estimator \(\hat {{\mathbf{\beta}} }\), we let \[ 0 = \frac {\partial }{\partial {\mathbf{\beta}} } \mathcal {L}({\mathbf{\beta}} ) \implies {\mathbf{X}} ^\top ({\mathbf{Y}} -{\mathbf{X}} {\mathbf{\beta}} )=0 \implies \hat {{\mathbf{\beta}} } = ({\mathbf{X}} ^\top {\mathbf{X}} )^{-1} {\mathbf{X}} ^\top {\mathbf{Y}} . \]

In the general case, the optimal function is \(f^*({\mathbf{x}} ) = {\mathbb{E}} [y \mid {\mathbf{x}} ]\), and the linear model approximates this conditional expectation as \({\mathbb{E}} [y \mid {\mathbf{x}} ] = {\mathbf{\beta}} ^\top {\mathbf{x}} \). To expand the function space and improve expressiveness, we sometimes apply a feature mapping to \({\mathbf{x}} \). A simple example of a feature mapping is the polynomial mapping.

Splines For a partition \(\zeta _0, \ldots , \zeta _K\) of the \(x\)-axis, define the indicator regressors \[ \mathbf {z}_n^0 = \begin {pmatrix} \mathrm {I}(\zeta _0 \leq x < \zeta _1) \\ \vdots \\ \mathrm {I}(\zeta _{K-1} \leq x < \zeta _K) \end {pmatrix}. \] Regressing \(y_n \sim \mathbf {z}_n^0\) produces a piecewise constant fit, returning the average of training responses within each interval. Indicator regressors give discontinuous fits with zero derivative. To obtain a piecewise linear fit, augment with \[ \mathbf {z}_n^1 = \begin {pmatrix} \mathrm {I}(\zeta _0 \leq x < \zeta _1)(x - \zeta _0) \\ \vdots \\ \mathrm {I}(\zeta _{K-1} \leq x < \zeta _K)(x - \zeta _{K-1}) \end {pmatrix}. \] Regressing \(y_n \sim \mathbf {z}_n^0 + \mathbf {z}_n^1\) yields a separate first-order Taylor approximation within each segment. Generalizing to degree \(p\), \[ \mathbf {z}_n^p = \begin {pmatrix} \mathrm {I}(\zeta _0 \leq x < \zeta _1)(x - \zeta _0)^p \\ \vdots \\ \mathrm {I}(\zeta _{K-1} \leq x < \zeta _K)(x - \zeta _{K-1})^p \end {pmatrix}, \] gives a \(p\)-th order polynomial fit within each bucket.

Bias and Variance Trade-off Let’s decompose this target \(\mathbb {E}\left [(\hat {f}(\mathbf {x}_{\mathrm {new}}) - y_{\mathrm {new}})^2\right ]\). \begin{align*} \mathbb {E}\left [(\hat {f}(\mathbf {x}_{\mathrm {new}}) - y_{\mathrm {new}})^2\right ] &= \mathbb {E}\left [(\hat {f}(\mathbf {x}_{\mathrm {new}}) - \bar {f}(\mathbf {x}_{\mathrm {new}}) + \bar {f}(\mathbf {x}_{\mathrm {new}}) - f^{\star }(\mathbf {x}_{\mathrm {new}}) + f^{\star }(\mathbf {x}_{\mathrm {new}}) - y_{\mathrm {new}})^2\right ] \\ &= \underbrace {\mathbb {E}\left [(\hat {f}(\mathbf {x}_{\mathrm {new}}) - \bar {f}(\mathbf {x}_{\mathrm {new}}))^2\right ]}_{\text {``variance''}} + \underbrace {\mathbb {E}\left [(\bar {f}(\mathbf {x}_{\mathrm {new}}) - f^{\star }(\mathbf {x}_{\mathrm {new}}))^2\right ]}_{\text {``bias''}} + \underbrace {\mathbb {E}\left [(f^{\star }(\mathbf {x}_{\mathrm {new}}) - y_{\mathrm {new}})^2\right ]}_{\text {``irreducible error''}}. \end{align*}

Define \(M_{xx} = {\mathbb{E}} [{\mathbf{x}} ^\top {\mathbf{x}} ]\) and \(M_{xy} = {\mathbb{E}} [{\mathbf{x}} ^\top y]\). Then \({\mathbf{\beta}} = M_{xx}^{-1}M_{xy}\). By CLT, \[ \hat {{\mathbf{\beta}} } - {\mathbf{\beta}} ^* = \left (\frac {1}{N} \sum _{n=1}^{N} \mathbf {x}_n \mathbf {x}_n^\top \right )^{-1} \frac {1}{N} \sum _{n=1}^{N} \mathbf {x}_n (y_n - \mathbf {x}_n^\top \mathbf {\beta }^*), \quad \sqrt {N}(\hat {{\mathbf{\beta}} } - {\mathbf{\beta}} ^*) \to \mathcal {N}(0, M_{xx}^{-1} U M_{xx}^{-1}), \] where \(U := \operatorname {Cov}({\mathbf{x}} _n \epsilon _n)\). We have \[ N \cdot \text {Variance} \approx \operatorname {tr}(M_{xx}^{-1/2} U M_{xx}^{-/2}) = \sigma ^2\operatorname {tr}(I_P) = P\sigma ^2 \implies \text {Variance} \approx \frac {P\sigma ^2}{N}. \]

4.1.2 Ridge (\(L_2\)) and Lasso (\(L_1\)) Regression

When \(P \gg N\), \({\mathbf{X}} ^\top {\mathbf{X}} \) is non-invertible, and the model is prone to overfitting—a common phenomenon in this regime. A simple remedy is to introduce a penalty term. We begin with the \(L_2\) penalty. \[ \mathcal {L}_{ridge}({\mathbf{\beta}} , \lambda ) = \|{\mathbf{Y}} -{\mathbf{X}} ^\top {\mathbf{\beta}} \|_2^2 + \lambda \|{\mathbf{\beta}} \|_2^2 \implies \hat {\mathbf{\beta}} _{ridge} = ({\mathbf{X}} ^\top {\mathbf{X}} + \lambda {\mathbf{I}} )^{-1} {\mathbf{X}} ^\top {\mathbf{Y}} . \] Assume \({\mathbf{Y}} = {\mathbf{X}} ^\top {\mathbf{\beta}} ^* + {\mathbf{\varepsilon}} \) with \(\mathbb {E}[{\mathbf{\varepsilon}} ] = {\mathbf{0}} \) and \(\mathrm {Cov}({\mathbf{\varepsilon}} ) = \sigma ^2 {\mathbf{I}} \). Let the singular value decomposition give \({\mathbf{X}} ^\top {\mathbf{X}} = {\mathbf{V}} {\mathbf{\Sigma}} ^2 {\mathbf{V}} ^\top \) with \({\mathbf{\Sigma}} ^2 = \mathrm {diag}(d_1, \ldots , d_P)\), so that \(({\mathbf{X}} ^\top {\mathbf{X}} + \lambda {\mathbf{I}} )^{-1} = {\mathbf{V}} ({\mathbf{\Sigma}} ^2 + \lambda {\mathbf{I}} )^{-1} {\mathbf{V}} ^\top \).

A direct computation gives \[ \mathbb {E}[{\hat{\mathbf{\beta}}} _{\text {ridge}}] - {\mathbf{\beta}} ^* = -\lambda {\mathbf{V}} ({\mathbf{\Sigma}} ^2 + \lambda {\mathbf{I}} )^{-1} {\mathbf{V}} ^\top {\mathbf{\beta}} ^*, \qquad \mathrm {Cov}({\hat{\mathbf{\beta}}} _{\text {ridge}}) = \sigma ^2 {\mathbf{V}} \, \mathrm {diag}\!\left (\frac {d_j}{(d_j+\lambda )^2}\right ) {\mathbf{V}} ^\top , \] which yields \[ \text {Bias: }\| \mathbb {E}[{\hat{\mathbf{\beta}}} _{\text {ridge}}] - {\mathbf{\beta}} ^* \|_2^2 = \sum _{j=1}^{P} \frac {\lambda ^2}{(d_j + \lambda )^2} ({\mathbf{v}} _j^\top {\mathbf{\beta}} ^*)^2, \qquad \text {Variance: }\mathrm {tr}\,\mathrm {Cov}({\hat{\mathbf{\beta}}} _{\text {ridge}}) = \sigma ^2 \sum _{j=1}^{P} \frac {d_j}{(d_j + \lambda )^2}. \] As \(\lambda \) increases, the variance decreases while the bias increases.

Another choice is \(L_1\) penalty, which we called LASSO, \[ \mathcal {L}_{lasso}({\mathbf{\beta}} , \lambda ) = \|{\mathbf{Y}} -{\mathbf{X}} ^\top {\mathbf{\beta}} \|_2^2 + \lambda \|{\mathbf{\beta}} \|_1. \]

4.1.3 Kernel Methods

Consider the Ridge Regression, by some Mathematics tricks \[ \hat {\mathbf{\beta}} _{ridge} = ({\mathbf{X}} ^\top {\mathbf{X}} + \lambda {\mathbf{I}} _P)^{-1} {\mathbf{X}} ^\top {\mathbf{Y}} = {\mathbf{X}} ({\mathbf{X}} {\mathbf{X}} ^\top + \lambda {\mathbf{I}} _N)^{-1} {\mathbf{Y}} . \] Note that we have replaced an \(P\times P\) inverse with an \(N\times N\) inverse, and may reduce computation (in some case). \({\mathbf{X}} {\mathbf{X}} ^\top = (x_i^\top x_j)_{i,j}\), replace \(x_i^\top x_j\) into \(k(x_i,x_j)\), and we have kernel function.

For a kernel \(k : {\mathcal{X}} \times {\mathcal{X}} \to \mathbb {R}\), we define an integral operator \({\mathcal{T}} _k : L^2_\mu ({\mathcal{X}} ) \to L^2_\mu ({\mathcal{X}} )\) as follows \[ {\mathcal{T}} _k f(x) = \mathbb {E}_{x' \sim \mu }[k(x, x') f(x')]. \]

This theorem ensures the existence of an eigendecomposition of a kernel \(k\), i.e., the corresponding integral operator \({\mathcal{T}} _k\). Note that \((\lambda _j)_{j \geq 1}\) and \((e_j)_{j \geq 1}\) are the eigenvalues and eigenfunctions of the integral operator \({\mathcal{T}} _k\) in the sense that

• \({\mathcal{T}} _k e_j = \lambda _j e_j\), i.e., \(\mathbb {E}_{x' \sim \mu }[k(x, x') e_j(x')] = \lambda _j e_j(x)\).
• \(\langle e_i, e_j \rangle _{L^2_\mu ({\mathcal{X}} )} = \mathbb {E}_{x \sim \mu }[e_j(x) e_i(x)] = \delta _{i,j}\).

Feature map. Mercer’s theorem gives a feature map for the kernel \(k\). Let \[ \Phi : {\mathcal{X}} \to \ell ^2, \qquad \Phi (x) = \left (\sqrt {\lambda _1} e_1(x), \sqrt {\lambda _2} e_2(x), \ldots , \sqrt {\lambda _j} e_j(x), \ldots \right )^\top . \] Then, \[ k(x, x') = \sum _{j=1}^{\infty } \sqrt {\lambda _j} e_j(x) \sqrt {\lambda _j} e_j(x') = \langle \Phi (x), \Phi (x') \rangle _{\ell ^2}. \]