Reproducing kernel Hilbert space

A reproducing kernel Hilbert space (RKHS) ${\displaystyle {\mathcal {H}}_{k}}$ is a Hilbert space of functions defined by a symmetric, positive-definite function ${\displaystyle k:{\mathcal {X}}\times {\mathcal {X}}\rightarrow \mathbb {R} }$ called the reproducing kernel such that the function ${\displaystyle k(\mathbf {x} ,\cdot )}$ belongs to ${\displaystyle {\mathcal {H}}_{k}}$ for all ${\displaystyle \mathbf {x} \in {\mathcal {X}}}$.^[6][7][8] There are three main properties make an RKHS appealing:

1. The reproducing property, which gives name to the space,

${\displaystyle f(\mathbf {x} )=\langle f,k(\mathbf {x} ,\cdot )\rangle _{k},\quad \forall \ f\in {\mathcal {H}}_{k},}$

where ${\displaystyle \langle \cdot ,\cdot \rangle _{k}}$ is the inner product in ${\displaystyle {\mathcal {H}}_{k}}$.

2. Functions in an RKHS are in the closure of the linear combination of the kernel at given points,

${\displaystyle f(\mathbf {x} )=\sum _{i}k(\mathbf {x} _{i},\mathbf {x} )c_{i}}$.

This allows the construction in a unified framework of both linear and generalized linear models.

3. The squared norm in an RKHS can be written as

${\displaystyle \|f\|_{k}^{2}=\sum _{i,j}k(\mathbf {x} _{i},\mathbf {x} _{j})c_{i}c_{j}}$

and could be viewed as measuring the complexity of the function.

The regularized functional

The estimator is derived as the minimizer of the regularized functional

${\displaystyle {\frac {1}{n}}\sum _{i=1}^{n}(f(\mathbf {x} _{i})-y_{i})^{2}+\lambda \|f\|_{k}^{2},}$

(2)

where ${\displaystyle f\in {\mathcal {H}}_{k}}$ and ${\displaystyle \|\cdot \|_{k}}$ is the norm in ${\displaystyle {\mathcal {H}}_{k}}$. The first term in this functional, which measures the average of the squares of the errors between the ${\displaystyle f(\mathbf {x} _{i})}$ and the ${\displaystyle y_{i}}$, is called the empirical risk and represents the cost we pay by predicting ${\displaystyle f(\mathbf {x} _{i})}$ for the true value ${\displaystyle y_{i}}$. The second term in the functional is the squared norm in a RKHS multiplied by a weight ${\displaystyle \lambda }$ and serves the purpose of stabilizing the problem^[3][5] as well as of adding a trade-off between fitting and complexity of the estimator.^[2] The weight ${\displaystyle \lambda }$, called the regularizer, determines the degree to which instability and complexity of the estimator should be penalized (higher penalty for increasing value of ${\displaystyle \lambda }$).

Derivation of the estimator

The explicit form of the estimator in equation (1) is derived in two steps. First, the representer theorem^[9][10][11] states that the minimizer of the functional (2) can always be written as a linear combination of the kernels centered at the training-set points,

${\displaystyle {\hat {f}}(\mathbf {x} ')=\sum _{i=1}^{n}c_{i}k(\mathbf {x} _{i},\mathbf {x} ')=\mathbf {k} ^{\top }\mathbf {c} ,}$

(3)

for some ${\displaystyle \mathbf {c} \in \mathbb {R} ^{n}}$. The explicit form of the coefficients ${\displaystyle \mathbf {c} =[c_{1},\ldots ,c_{n}]^{\top }}$ can be found by substituting for ${\displaystyle f(\cdot )}$ in the functional (2). For a function of the form in equation (3), we have that

${\displaystyle {\begin{aligned}\|f\|_{k}^{2}&=\langle f,f\rangle _{k},\\&=\left\langle \sum _{i=1}^{N}c_{i}k(\mathbf {x} _{i},\cdot ),\sum _{j=1}^{N}c_{j}k(\mathbf {x} _{j},\cdot )\right\rangle _{k},\\&=\sum _{i=1}^{N}\sum _{j=1}^{N}c_{i}c_{j}\langle k(\mathbf {x} _{i},\cdot ),k(\mathbf {x} _{j},\cdot )\rangle _{k},\\&=\sum _{i=1}^{N}\sum _{j=1}^{N}c_{i}c_{j}k(\mathbf {x} _{i},\mathbf {x} _{j}),\\&=\mathbf {c} ^{\top }\mathbf {K} \mathbf {c} .\end{aligned}}}$

We can rewrite the functional (2) as

${\displaystyle {\frac {1}{n}}\|\mathbf {y} -\mathbf {K} \mathbf {c} \|^{2}+\lambda \mathbf {c} ^{\top }\mathbf {K} \mathbf {c} .}$

This functional is convex in ${\displaystyle \mathbf {c} }$ and therefore we can find its minimum by setting the gradient with respect to ${\displaystyle \mathbf {c} }$ to zero,

${\displaystyle {\begin{aligned}-{\frac {1}{n}}\mathbf {K} (\mathbf {Y} -\mathbf {K} \mathbf {c} )+\lambda \mathbf {K} \mathbf {c} &=0,\\(\mathbf {K} +\lambda n\mathbf {I} )\mathbf {c} &=\mathbf {Y} ,\\\mathbf {c} &=(\mathbf {K} +\lambda n\mathbf {I} )^{-1}\mathbf {Y} .\end{aligned}}}$

Substituting this expression for the coefficients in equation (3), we obtain the estimator stated previously in equation (1),

${\displaystyle {\hat {f}}(\mathbf {x} ')=\mathbf {k} ^{\top }(\mathbf {K} +\lambda n\mathbf {I} )^{-1}\mathbf {Y} .}$

A review of Bayesian probability

As part of the Bayesian framework, the Gaussian process specifies the prior distribution that describes the prior beliefs about the properties of the function being modeled. These beliefs are updated after taking into account observational data by means of a likelihood function that relates the prior beliefs to the observations. Taken together, the prior and likelihood lead to an updated distribution called the posterior distribution that is customarily used for predicting test cases.

The Gaussian process

A Gaussian process (GP) is a stochastic process in which any finite number of random variables that are sampled follow a joint Normal distribution.^[12] The mean vector and covariance matrix of the Gaussian distribution completely specify the GP. GPs are usually used as a priori distribution for functions, and as such the mean vector and covariance matrix can be viewed as functions, where the covariance function is also called the kernel of the GP. Let a function ${\displaystyle f}$ follow a Gaussian process with mean function ${\displaystyle m}$ and kernel function ${\displaystyle k}$,

${\displaystyle f\sim {\mathcal {GP}}(m,k).}$

In terms of the underlying Gaussian distribution, we have that for any finite set ${\displaystyle \mathbf {X} =\{\mathbf {x} _{i}\}_{i=1}^{n}}$ if we let ${\displaystyle f(\mathbf {X} )=[f(\mathbf {x} _{1}),\ldots ,f(\mathbf {x} _{n})]^{\top }}$ then

${\displaystyle f(\mathbf {X} )\sim {\mathcal {N}}(\mathbf {m} ,\mathbf {K} ),}$

where ${\displaystyle \mathbf {m} =m(\mathbf {X} )=[m(\mathbf {x} _{1}),\ldots ,m(\mathbf {x} _{N})]^{\top }}$ is the mean vector and ${\displaystyle \mathbf {K} =k(\mathbf {X} ,\mathbf {X} )}$ is the covariance matrix of the multivariate Gaussian distribution.

Derivation of the estimator

In a regression context, the likelihood function is usually assumed to be a Gaussian distribution and the observations to be independent and identically distributed (iid),

${\displaystyle p(y|f,\mathbf {x} ,\sigma ^{2})={\mathcal {N}}(f(\mathbf {x} ),\sigma ^{2}).}$

This assumption corresponds to the observations being corrupted with zero-mean Gaussian noise with variance ${\displaystyle \sigma ^{2}}$. The iid assumption makes it possible to factorize the likelihood function over the data points given the set of inputs ${\displaystyle \mathbf {X} }$ and the variance of the noise ${\displaystyle \sigma ^{2}}$, and thus the posterior distribution can be computed analytically. For a test input vector ${\displaystyle \mathbf {x} '}$, given the training data ${\displaystyle S=\{\mathbf {X} ,\mathbf {Y} \}}$, the posterior distribution is given by

${\displaystyle p(f(\mathbf {x} ')|S,\mathbf {x} ',{\boldsymbol {\phi }})={\mathcal {N}}(m(\mathbf {x} '),\sigma ^{2}(\mathbf {x} ')),}$

where ${\displaystyle {\boldsymbol {\phi }}}$ denotes the set of parameters which include the variance of the noise ${\displaystyle \sigma ^{2}}$ and any parameters from the covariance function ${\displaystyle k}$ and where

${\displaystyle {\begin{aligned}m(\mathbf {x} ')&=\mathbf {k} ^{\top }(\mathbf {K} +\sigma ^{2}\mathbf {I} )^{-1}\mathbf {Y} ,\\\sigma ^{2}(\mathbf {x} ')&=k(\mathbf {x} ',\mathbf {x} ')-\mathbf {k} ^{\top }(\mathbf {K} +\sigma ^{2}\mathbf {I} )^{-1}\mathbf {k} .\end{aligned}}}$