Basis Expansion¶
Table of Contents
Warning
We can use the established framework for linear models to move beyond linearity just by using non-linear transforms on the features.
We define the transforms as functions \(h:\mathbb{R}^d\mapsto\mathcal{\Phi}\) where \(\mathcal{\Phi}\) can be finite or infinite dimensional.
Since these methods expand the basis from \(d\)-dimensions to \(\dim\mathcal{\Phi}\), these methods are called basis expansion methods.
Tip
For most of the methods, we consider the case where \(d=1\).
Finite Dimensional Expansion¶
Note
We define a finite set of transforms \(x\mapsto h_m(x)\) where \(m=1,\cdots,M\).
Each observation is mapped to the transformed vector \(\mathbf{h}=(h_1(x),\cdots,h_M(x))^T\)
We use the linear regression or logistic regression frameworks using these transforms
Regression: \(\hat{y}=f(x)=\sum_{i=1}^M \beta_mh_m(x)=\boldsymbol{\beta}^T\mathbf{h}\)
Classification: \(\log\frac{\mathbb{P}(G=k|X=x)}{\mathbb{P}(G=K|X=x)}=\boldsymbol{\beta}^T\mathbf{h}\)
Polynomial¶
Tip
[Taylor’s series] An infinite power series can arbitrarily approximate any real-valued infinitely differentiable function around a point \(a\).
\[f(x)=f(a)+\frac{f'(a)}{2!}(x-a)+\frac{f''(a)}{3!}(x-a)^2+\cdots\][Lagrange Polynomial] A polynomial of degree \(n\) can be fit to pass through exactly \(n\) points.
Say, we have a dataset \(\mathbf{X}=[x_1,\cdots,x_N]^T\) and target \(\mathbf{y}=[y_1,\cdots,y_N]^T\).
For each point \((x_k,y_k)\), we create the Lagrange coefficients
\[l_k(x)=\frac{(x-x_1)(x-x_2)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_{N-1})(x-x_N)}{(x_k-x_1)(x_k-x_2)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_{N-1})(x_k-x_N)}\]The polynomial passing through all these points (interpolation) is given by
\[f(x)=\sum_{i=1}^N y_1 l_i(x)\]
Warning
A polynomial of degree \(n\) has \(n-1\) turns.
With real data, often not all the turns are utilised for the points in the interior regions.
With no other constraints applied, polynomials of such higher degrees oscilate crazily towards the extremes where the data density is usually low.
Therefore, in those regions, it provides very poor generalisation of the model.
The situation is often worse if we move to higher dimensions, as the proportion of points around the exterior increases in higher dimension.
Note
With these in mind, a polynomial of smaller degree often works best, such as cubic polynomials.
\[x\overset{h}\mapsto[1,x,x^2,x^3]\]However, they often show higher error rate due to increased bias.
Piece-wise Polynomials¶
Note
We can split the input domain into \(M\) regions using knots, \(\xi_1,\cdots,\xi_{M-1}\).
We choose a function family \(\mathcal{H}\) and fit it separately as \(h_i\) in each region \((\xi_{i-1},\xi_i)\).
Warning
The knot-selection problem has to be taken care of separately.
In all such cases, the functions often becomes discontinuous at the knot points which is undesirable for generalisation.
One-hot encoding for regions¶
Note
In this case, the basis expansion works with indicator functions
\(h_1=I(x < \xi_1)\)
For \(i=2,\cdots,M-1, h_i(x)=I(\xi_{i-1} < x < \xi_i)\)
\(h_M(x)=I(\xi_M < x)\)
For regression problems, the fitted co-efficients are the average of the observation in the current piece.
\[\hat{\beta_i}=\frac{\sum_{x_k\in(\xi_{i-1},\xi_i)}y_k}{|x_k\in(\xi_{i-1},\xi_i)|}\]
Polynomial function for regions¶
Note
We can also design polynomials of any degree \(\mathcal{P}(n)\) to each of these regions.
\(h_1=I(x < \xi_1)\mathcal{P}_1(n)\)
For \(i=2,\cdots,M-1, h_i(x)=I(\xi_{i-1} < x < \xi_i)\mathcal{P}_i(n)\)
\(h_M(x)=I(\xi_M < x)\mathcal{P}_M(n)\)
Fitting this model with MSE for regression fits a separate polynomial in each region.
The fitted function is often discontinuous at the knot points.
Polynomial Spline¶
Note
Here we design the functions for each region in such a way so that the function becomes continuous at each of the knot points.
The key idea is to define additional polynomials of the same target degree in such a way that, for any given knot-point, the function stays 0 to the left of it but continuously becomes non-zero on the right of it.
For cubic splines, the functions are defined as:
\(h_1(x)=1\)
\(h_2(x)=x\)
\(h_3(x)=x^2\)
\(h_4(x)=x^3\)
\(h_{k+4}(x)=(x-\xi_k)^3_+\), for \(k=1,\cdots,M\) where
\[\begin{split}(x-\xi_k)^3_+=\begin{cases}0 & \text{if } x < \xi_k\\ (x-\xi_k)^3 & \text{if } x \ge \xi_k\end{cases}\end{split}\]
Natural Spline¶
Note
Since each region of a polynomial spline is fit with less data, often they show crazier behaviour near the boundaries than global polynomials.
To alleviate these problems, natural splines model the function as a linear function for the left of the leftmost and the right of the rightmost knot points.
We use the notation \(N_i(x)\) instead of \(h_i(x)\) to emphasize that we’re working with natural splines.
Tip
[TODO] Note on the number of parameters and degrees of freedom.
Smoothing Spline¶
Tip
For each of the piece-wise fitting approaches, knot selection remains a key-issue.
Smoothing splines address this by implicitly allowing a knot at every single data-point.
Since this approach can potentially create much higher degree polynomials, the complexity of the model is controlled via regularisation.
Note
The functions are restricted to be twice-differentiable equipped with norm (e.g. Sobolev space \(W^{k=2,p=2}\)).
Any \(f\in W^{k,p}\) and its (weak) derivatives up to order \(k\) have a finite \(L_p\) norm.
The objective function is defined as
\[\hat{f}=\min_{f\in W^{k=2,p=2}}\left[\sum_{i=1}^N(y-f(x))^2+\lambda\int\left(f''(z)\right)^2\mathop{dz}\right]\]The smoothness is captured in the double-derivative since it represents curvature.
\(\lambda\in[0,\infty]\) is a smoothing parameter which controls the model complexity
\(\lambda=0\): Rough fit, equivalent to interpolation using a Lagrange polynomial.
\(\lambda=\infty\): Linear fit since it reduces to an OLS problem with MSE loss.
Note
[TODO: Proof?] The solution for this is Natural splines
\[f(x)=\sum_{i=1}^N \beta_iN_i(x)\]We have \(f''(z)=\sum_{i=1}^N \beta_iN''_i(z)\) and
\[\int\left(f''(z)\right)^2\mathop{dz}=\sum_{i=1}^N \sum_{j=1}^N\beta_i\beta_j\int N''_i(z)N''_j(z)\mathop{dz}=\boldsymbol{\beta}^T\boldsymbol{\Omega}_N\boldsymbol{\beta}\]Here \(\{\boldsymbol{\Omega}_N\}_{i,j}=\int N''_i(z)N''_j(z)\mathop{dz}\)
The objective function can therefore be written as a generalised ridge regression
\[\min_\beta\left[(\mathbf{y}-\mathbf{N}\boldsymbol{\beta})^T(\mathbf{y}-\mathbf{N}\boldsymbol{\beta})+\lambda\boldsymbol{\beta}^T\boldsymbol{\Omega}_N\boldsymbol{\beta}\right]\][TODO: Write the final solution]
Non-linear Classification¶
Note
TODO
Moving Beyond 1 Dimension¶
Note
TODO Tensor product
Infinite Dimensional Expansion¶
Note
More general regression problems can be formulated using a similar framework for the smoothing splines
\[\hat{f}=\min_{f\in\mathcal{H}}\left[\sum_{i=1}^NL(y,f(x))+\lambda J(f)\right]\]Here \(\mathcal{H}\) is a function class (hypothesis space) (e.g. Sobolev space for smoothing splines).
\(J(f)\) is a regulariser which penalises functions for being too complex (to avoid overfitting).
\(\lambda\) is the regulariser parameter which controls the trade-off between the bias and the variance.
A Point Mapping to a Function¶
Tip
We usually think of points on the number line or in a higher dimensional space (\(x\in\mathbb{R}^d\)) as a discrete entity.
Instead, a point can be thought of as an impulse like Dirac-delta function, which has infinitely high probability of being found at it’s location.
\[\begin{split}\delta_x(t)=\begin{cases}+\infty & t=x \\ 0 & t\neq x\end{cases}\end{split}\]It has an infinitely sharp peak at \(x\) but that dies off immediately to \(0\) everywhere else.
[Side-node] This is a special case of generalised functions or distributions.
We can think of this impulse as a limit of a sequence of functions, \(\lim_\limits{\gamma\downarrow 0} f_\gamma\), such as Gaussian bumps around the point \(x\).
\[f_\gamma=\exp\left(-\gamma||x-x'||^2\right)\]While thinking of a map from a point to a function, we’re essentially going backwards from the impulse limit to a bump with a non-zero width.
Intuitively, this allows for some uncertainty about the exact location of the point.
Reproducing Kernel Hilbert Space¶
Note
We choose RKHS (\(\mathcal{H}_K\)) to be our function class.
We define the basis expansion as a mapping from the original domain to this function space
\[x\in\mathbb{R}^d\overset{h}\longmapsto (f:\mathbb{R}^d\mapsto\mathbb{R})\in\mathcal{H}_K\]For each point \(x_i\in\mathbf{X}\), we map it to a function \(f_i(\cdot)\in\mathcal{H}_K\) with a free parameter.
Warning
Let \(K\) be the reproducing kernel of \(\mathcal{H}_K\) with Dirac evaluation functional \(\delta_x:\mathcal{H}_K\mapsto\mathbb{R}\).
This means, for any \(f\in \mathcal{H}_K\),
\[\delta_x(f)=\langle K(\cdot,x), f\rangle_{\mathcal{H}_K}=f(x)\]We choose a specific mapping \(\varphi:\mathbb{R}^d\mapsto\mathcal{H}_K\) such that it maps each point \(x\) to this specific type of functions (canonical feature maps)
\[\varphi(x)=f_{\delta_x}=K(\cdot,x)\]We note that the kernel also provides the tool to calculate the inner products (and, hence, similarity measure via a metric) between the points.
\[\langle \varphi_i(\cdot), \varphi_j(\cdot)\rangle_{{\mathcal{H}}_K}=\langle K(\cdot,x_i), K(\cdot,x_j)\rangle_{{\mathcal{H}}_K}=K(x_i,x_j)\]
Note
\(\mathcal{H}_K\) contains (uncountably) infinite number of basis functions \(K(\cdot,x_i)\) for \(i\in\mathcal{I}\).
Every function \(f\in\mathcal{H}_K\) can be expressed using those basis
\[f(\cdot)=\sum_{i\in\mathcal{I}}\alpha_i K(\cdot,x_i)\]
Mercer Kernels: Eigen-decomposition of Linear Space Spanned by Kernel Functions¶
Note
From a potentially uncountable number of basis functions, eigen-decomposition reduces the basis functions to a countable case.
Requires added structural assumption on the domain of \(x\) and on the continuity of \(K\).
Assuming that the kernel has an eigen-decomposition with eigenfunctions \((\phi_i)_{i=1}^\infty\in\mathcal{H}_K\), the kernel also be written in term of these
\[K(x,y)=\sum_{i=1}^\infty \gamma_i\phi_i(x)\phi_i(y)\]We note that \(\varphi(x)\neq \phi(x)\), and therefore, this produces a different feature map.
Since kernels are symmetric and positive definite, the eigenvalues
are positive, i.e. \(\gamma_i\ge 0\), and
have bounded sum, i.e. \(\sum_{i=1}^\infty \gamma_i < \infty\)
Any function \(f\in\mathcal{H}_K\) can be expressed as a linear combination of the countable eigenfunctions
\[f(x)=\sum_{i=1}^\infty c_i\phi_i(x)\]
Kernel Ridge Regression¶
Note
We use the function norm as the regulariser as this captures how vigorously the function oscilate along the direction of each eigenfunctions.
\[||f||^2_{\mathcal{H}_K}\overset{\text{def}}=\sum_{i=1}^\infty c_i^2/\gamma_i < \infty\]
Note
The ridge regression problem using functions from kernel family can be expressed as
\[\hat{f}=\min_{f\in\mathcal{H}}\left[\sum_{i=1}^N L(y,f(x_i))+\lambda ||f||^2_{\mathcal{H}_K}\right]\]This reduces to
\[\hat{f}=\min_{(c_k)_{k=1}^\infty}\left[\sum_{i=1}^N L(y,\left(\sum_{k=1}^\infty c_k\phi_k(x_i)\right))+\lambda \sum_{k=1}^\infty c_i^2/\gamma_i\right]\]On a face-value, it seems that we’d need to estimate an infinite number of parameters.
[TODO: Proof?] However, the solution reduces to a finite dimensional one from a potentially uncountable one
\[f(x)=\sum_{i=1}^N\alpha_iK(x,x_i)=\mathbf{K}\boldsymbol{\alpha}\]The regulariser term is defined as
\[\lambda||f||_{\mathcal{H}_K}^2=\lambda\langle f,f\rangle_{\mathcal{H}_K}=\sum_{i=1}^N\sum_{j=1}^N K(x_i,x_j)\alpha_i\alpha_j=\lambda\boldsymbol{\alpha}^T\mathbf{K}\boldsymbol{\alpha}\]The minimization problem therefore becomes
\[\hat{f}=\min_{\boldsymbol{\alpha}}L(\mathbf{y}, \mathbf{K}\boldsymbol{\alpha})+\lambda\boldsymbol{\alpha}^T\mathbf{K}\boldsymbol{\alpha}\]We note that for MSE loss, this reduces to a generalised ridge regression problem.
Connection to Smoothness¶
Warning
[TODO] Intuition: Lesser weights to more wiggly component functions in the solution.
See also
For intuitive understanding, here is an example:
We have 1-dimensioanl data matrix \(\mathbf{X}=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}1.1\\1.5\\2.1\end{bmatrix}\).
We define the map \(x_i\overset{h}\mapsto f_i(\cdot)\) using Gaussian RBF kernel.
We define our 3 basis functions as
\[\mathbf{h}=[K(\cdot,x_1),K(\cdot,x_2),K(\cdot,x_3)]=\begin{bmatrix}\exp(-\frac{||x-1.1||^2}{0.05}),\exp(-\frac{||x-1.5||^2}{0.05}),\exp(-\frac{||x-2.1||^2}{0.05})\end{bmatrix}\]One linear combination of these functions is shown here in yellow.
\[f(x)=3h_1(x)+2h_2(x)-h_3(x)\]
Each datapoint \(x_i\) is the mapped to \(\mathbf{h_i}=[K(x_i,x_1),K(x_i,x_2),K(x_i,x_3)]\).