Book Notes on Pattern Recognition and Machine Learning (PRML)
在之前的模型中，我们利用训练集的数据来训练我们的模型 \(\beta\) ，训练结束后就可以丢弃掉训练数据，只保留模型用于预测。而还有一种方式是保留全部或部分训练集数据，并在预测阶段使用到这些数据。这样一种方式称为 memory-based method，而 kernel method 就是这样一种方式。
我们从线性模型中引入 kernel 的概念。和之前的假设一致，我们假设输入为 n 个数据，每个数据为 p 维；输出为 n 个一维数据。在之前介绍的线性模型中，对于训练集 \(X\) 中的第 i 个数据 \(X_i = [x_{i0},x_{i1},...,x_{ip}]^T\) ，我们直接将其输入到线性模型中训练模型 \(\beta\) ，此时该数据对应的 loss 可写为
\[L_{i}(\beta) = \left[ \sum_{j=0}^p (x_{ij}\beta_{j}) \right] - y_i^*\]但有时我们会对输入数据进行预处理，会引入非线性的因素。比如我们想对变换后的数据 \(X_i^{'} = [x_{i0},x_{i1}^2,...,\sqrt{x_{ip}}]^T\) 进行回归分类。此时我们仍然可以使用线性模型，只需要将变换后的数据 \(X^{'}\) 看成新的输入即可。我们用函数 \(\phi(x)\) 来表征这样一种变换，即 \(\phi(X_i) = [\phi_0(x_{i0}),\phi_1(x_{i1}),...,\phi_p(x_{ip})]^T\) 。我们写出ridge正则下基于所有训练数据的loss
\[L(\beta) = ||Y-\phi(X)\beta||^2 + \lambda||\beta||^2 =\sum_{i=1}^n(\phi(X_i)^T\beta-y_i^*)^2 + \lambda\beta^T\beta\]令损失函数导数为零可得
\[\beta = -\frac{1}{\lambda}\sum_{i=1}^n(\phi(X_i)^T\beta-y_i^*)\phi(X_i) = \sum_{i=1}^n c_i\phi(X_i) = \phi(X)^Tc\]和之前 X 的定义相同，这里的 \(\phi(X) =[\phi(X_1)^T,...,\phi(X_n)^T]^T\) 为 n*(p+1) 的矩阵。而这里 \(c = [c_1,...,c_n]^T\) ，我们用 \(c_i\) 来表示 \(\phi(X_i)\) 前面的系数，即
\[c_i = -\frac{1}{\lambda}(\phi(X_i)^T\beta-y_i^*)\]我们将 \(\beta = \phi(X)^Tc\) 代入到损失函数中，此时自变量将由 \(\beta\) 转换为 \(c\) ，我们有
\[L(c) = c^T\phi\phi^T\phi\phi^Tc-2c^T\phi\phi^TY+Y^TY+\lambda c^T\phi\phi^Tc\]我们用 \(K=\phi\phi^T\) 来代替，那么我们就引出了kernel的概念。显然，这里K是一个 n*n 的对称矩阵，其中的某个元素 \(K_{ij} = \phi(X_i)^T\phi(X_j)\) 。对于kernel，我们用kernel function来刻画，用来描述两个变量的关系，即 \(k(x_i,x_j) =\phi(x_i)^T\phi(x_j)\) 。所以此时的 K 矩阵也称为 kernel matrix。
此时的损失函数为
\[L(c) = ||Y-Kc|| + \lambda c^TKc\]此时预测值 \(y_{pred}(x) =\sum_{i=1}^nc_ik(x,X_i)\)
新的模型 c 的解析解为 \(c = (K+\lambda I_n)^{-1}Y\)
我们可以看到，此时我们需要求一个 N*N 矩阵的逆，这在大数据集中几乎是不可能的。但我们可以发现，在解析式和预测方程中，可以直接引入kernel function，而不需要明确 \(\phi(X)\) 。这个性质可以让我们隐式地使用高维甚至无穷维的特征空间。
为了利用 kernel method 的好处，我们需要构造出合法的kernel function。最简单的方式是我们先构造映射函数 \(\phi(x)\) ，然后利用 \(k(x_1,x_2) = \phi(x_1)^T\phi(x_2)\) 。而另外一种方法是我们直接构造 kernel function，但在此时我们怎样确定这是合法的kernel function 呢？
如果我们可以把 kernel function 拆分成两个映射函数的乘积，那么毫无疑问这种 kernel function 是合法的。但这种拆分并不容易，而且没有普适的拆分方法。不过还是有判断有效kernel的充要条件，那就是看kernel matrix是否是半正定的。这里的kernel matrix和上一节一致，即 \(K_{ij} =k(x_i,x_j)\) 。如果 K 是半正定的，那么对应的kernel是合法的。
接下来介绍几种常用的kernel
首先是最基础的 polynomial kernel，\(k(x,x^{'}) = (x^Tx^{'}+c)^M\) ，这里 c 是常数而 M 为阶数。如果该kernel计算式中不包含 c 这一项，则对应的映射函数 \(\phi(x)\) 只包含所有 M 阶的项。在加上 c 这一项后，包含所有阶数小于等于 M 的项。
其次是最常用的 Gaussian kernel，\(k(x,x^{'}) = exp(\frac{-\vert\vert x-x^{'}\vert \vert ^2}{2\sigma^2})\) 。广泛使用的一个原因是这个kernel对应的特征向量有无穷的维度。我们可以用高斯展开来解释这个原因。
\[k(x_1,x_2) = exp(-x_1^Tx_1/2\sigma^2)\{1+\frac{x_1^Tx_2}{\sigma^2}+\frac{(\frac{x_1^Tx_2}{\sigma^2})^2}{2!} +...\}exp(-x_2^Tx_2/2\sigma^2)\]而且，高斯核并不仅仅局限于欧氏距离，我们甚至可以用另一个非线性核代替高斯核内的 \(x_1^Tx_2\)，即
\[k(x_1,x_2) = exp\{-\frac{1}{2\sigma^2}(k'(x_1,x_1)+k'(x_2,x_2)-2k'(x_1,x_2)) \}\]同时，kernel method 也可以作用与离散变量。我们假设有两个集合 \(A_1\) 和 \(A_2\) ，此时我们可以定义kernel
\[k(A_1,A_2)=2^{|A_1\cap A_2|}\]最后，我们可以用生成模型来生成kernel，比如 Fisher kernel，或者是
\[k(x_1,x_2)=\sum_i p(x_1\vert i)p(x_2\vert i)p(i)\]radial basis function (rbf) 是指只依赖于径向距离(通常是欧氏距离)的basis function，例如 \(\phi_i(x)=h(||x-\mu_i||)\)
在历史上，rbf的引入是为了完全拟合训练集所有数据。我们通过 rbf 的线性组合来得到对应的函数 \(f(x)\) 即
\[f(x) = \sum_{i=1}^nw_ih(||x-x_i||)\]这里的系数 \(w_i\) 通过最小二乘确定。但显然这种方式在模式识别中是不合适的，一方面是计算量要求大，每个训练集数据都有一个待优化系数对应；另一方面是过拟合现象严重。
为此我们的解决方法是只选取有限个rbf进行线性组合，这就涉及到 rbf 中心点的选取（训练数据子集的选取）。我们可以用 k-means 聚类来选数据点，或者在每次迭代中选择使loss下降最大的那个数据
接下来我们从 Nadaraya-Watson model 的角度分析kernel回归模型。
对于训练集 \(\{X,Y\}\) ，我们用密度函数来建模联合分布 \(P(X,Y)\)，即
基于这个联合分布，我们可以表示出回归函数 \(y(x)\) ，即
\[y(x) = E[y\vert x] = \int_{-\infty}^{\infty}yP(y\vert x)dy = \frac{\int y P(x,y)dy}{\int P(x,y)dy} = \frac{\sum_n g(x-x_n)y_n}{\sum_m g(x-x_m)}\]这里我们用到了一个假设，假设密度函数的均值为零，即 \(\int_{-\infty}^{\infty}f(x,y)dy=0\) 基于这个假设，我们定义 \(g(x) = \int_{-\infty}^{\infty}f(x,y)dy\) 得到上式。此时我们定义 kernel function 为
\[k(x,x_i) = \frac{g(x-x_i)}{\sum_{j=1}^ng(x-x_j)}\]就可以得到 Nadaraya-Watson model（kernel regression）的形式
\[y(x) = \sum_{i=1}^n k(x,x_i)y_i^*\]从高斯过程的角度来看，我们不需要求解具体的模型，而是直接定义函数的先验分布。（ we can obtain the predictive distribution either by taking a parameter space viewpoint and using the linear regression result or by taking a function space viewpoint and using the Gaussian process result.）我们以回归任务为例，\(y_i^* = y_{pred}(x_i) + \epsilon_i\) ，我们假设噪声 \(\epsilon\) 服从高斯分布，则有
\[p(y^*_i\vert y_i) = N(y^*_i\vert y_i,\alpha^{-1})\]这里 \(\alpha\) 是控制噪声的超参，我们写成矩阵的形式为 \(P(Y^\star\vert Y) = N(Y^\star\vert Y,\alpha^{-1}I_n)\) 。而根据高斯过程的定义，我们同样可以写出边缘分布 \(P(Y) = N(Y\vert 0,K)\) （这源于我们对模型的先验分布为高斯分布的假设）。这里的 K 就是kernel matrix，我们需要选择合适的 kernel function 使得相似的两个数据 \(X_i\) 和 \(X_j\) 对应的预测值接近。在高斯回归中最常使用的 kernel function 是
\[k(x_i,x_j) = \theta_0exp\{-\frac{\theta_1}{2}||x_i-x_j||^2\} + \theta_2 + \theta_3 x_i^Tx_j\]接下来我们就可以求 ground-truth 的边缘分布 \(P(Y^\star)\) 。这用到了一条结论
\[If \qquad p(x) = N(x\vert \mu, M^{-1}),\quad p(y\vert x) = N(y\vert Ax+b,N^{-1})\] \[Then \qquad p(y) = N(y\vert A\mu +b,N^{-1}+AM^{-1}A^T)\]基于这个结论，我们也可以写出 ground-truth 的边缘分布，\(P(Y^{\star})=N(Y^{\star}\vert 0,C)\) ，这里的 C 代表协方差矩阵，\(C(x_i,x_j) = k(x_i,x_j)+\alpha^{-1}I_{ij}\)
以上都是针对训练集中的数据进行的变换。我们更关心的是在未知新数据上的表现，也就是说，对于一个新数据 \(x_{new}\) ，我们想知道它对应的 ground-truth 的分布情况。我们利用训练集中的数据进行预测，即确定 \(p(y^{\star}_{new}\vert Y^\star)\)
我们在这里还需要引入一个结论
\[If \qquad p(x_1) \sim N(\mu_1, \Sigma_{11}),\quad p(x_2) \sim N(\mu_2,\Sigma_{22})\] \[Then \qquad p(x_2\vert x_1) \sim N(\mu_2+\Sigma_{21}\Sigma_{11}^{-1}(x_1-\mu_1), \Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12})\]结合上面的变换，我们已经知道了
\[P(Y^*) \sim N(0,C) \qquad P(y^*_{new})\sim N(0,k(x_{new},x_{new})+\alpha^{-1})\]那么代入结论中的通式就可以得到
\[p(y^*_{new}\vert Y^*) \sim N(\mathbf{k}^TC^{-1}Y^*,c-\mathbf{k}^TC^{-1}\mathbf{k})\]这里 \(\mathbf{k} = (k(x_1,x_{new}),...,k(x_n,x_{new}))^T\) ，而 \(c = k(x_{new},x_{new})+\alpha^{-1}\)
自此就得到了针对未知新数据的预测分布，下一步就是超参的学习。我们用 \(\theta\) 表示高斯过程模型中的超参，我们可以写出极大似然对应的 log-likelihood 形式
\[\ln P(Y\vert \theta) = -\frac{1}{2}\ln |C_N| - \frac{1}{2} Y^TC_N^{-1}Y -\frac{N}{2} \ln (2\pi)\]这里同样考虑利用梯度来优化我们的目标，我们可以写出上式对应的梯度为
\[\frac{\partial}{\partial \theta_i}\ln P(Y\vert \theta) = -\frac{1}{2}Tr(C_N^{-1}\frac{\partial C_N}{\partial \theta_i}) + \frac{1}{2}Y^TC_N^{-1}\frac{\partial C_N}{\partial \theta_i}C_N^{-1}Y\]以上都基于我们对误差的方差为常数的假设，如果误差也取决于输入变量，我们需要引入二次高斯过程来建模方差对输入的依赖。