Chapter 4 Linear Methods for Classification
前一章节介绍了针对 quantitative 变量的线性回归 Linear Regression,本章节讲适用于 qualitative 变量的线性分类 Linear Classification。这两个章节是本书所有内容的基础。根据我的理解,regression 相当于是优化问题中的 linear programming,而 classification则对应 integer programming。根据优化理论,两者的理论基础是共通的,前者如果是 convex problem的话就可以直接求得唯一最优解,后者没有统一的最优解求解方法,并且最优解不唯一。优化问题的很多难题都是在 integer programming 领域,那么对应的,我推测 learning接下来的很多困难都来自classification。
假设有编号为 $1,2,\cdots,K$ 的 $K$ 个分类,线性模型匹配后得到 $$ \begin{equation} \hat{f}_k(x) = \hat{\beta}_{k0}+\hat{\beta}_k^Tx. \end{equation} $$ 那么决定输入 $x$ 属于 $k$ 还是 $l$ 的分类边界条件为 $\hat{f}_k(x)=\hat{f}_l(x)$。该条件是一个 affine 集合或者 hyperplane,即为 $$ \begin{align} \left\{x:\left(\hat{\beta}_{k0}-\hat{\beta}_{l0}\right)+\left(\hat{\beta}_{k}-\hat{\beta}_{l}\right)^Tx=0\right\}. \end{align} $$
该分类方法其实是一种回归方法,其为每一类定义判别函数 discriminant functions $\delta_k(x)$,或者后验概率 posterior probabilities $Pr\left(G=k|X=x\right)$(这里分类值用G,以区分回归值Y),然后选择函数值或者概率最大的分类。当判别函数或者后验概率在$x$上是线性的时候,分类的决策边界就是线性的。
边界条件 $\hat{f}_k(x)=\hat{f}_l(x)$ 的判别通常有两种方法:减法 $\hat{f}_k(x)-\hat{f}_l(x)=0$ 和除法 $\frac{\hat{f}_k(x)}{\hat{f}_l(x)}=1$。在计算上,我们通常偏向于除法(理由以后有兴趣再去考究,感觉是除法可以消除好多的公因子,简化计算)。除法运算后得到的值为 1,通常我们希望和 0 进行比较,于是就引入了对数变换 logit transformation, $\log{\frac{\hat{f}_k(x)}{\hat{f}_l(x)}}=0$。
然后,最大的匹配值即为对应的分类 $$ \begin{align} \hat{G}(x)=\arg \max_{k\in G}\hat{f}_k(x). \end{align} $$ 这个方法显然是行不通的,不然就不用研究 classification了。其本质类似于用 linear programming 的方法来求解 integer programming 的问题,很多情况下得不到最优解。只有小部分的 integer programming 问题可以直接等价与 linear programming 问题,其它的都不可行,甚至连次优解都得不到。
假设$f_k(x)$满足独立标准分布(高斯分布)$(\mu_k,\mathbf{\Sigma})$, $$ \begin{align} f_k(x) = \frac{1}{(2\pi)^{p/2}\left|\mathbf{\Sigma}_k\right|^{1/2}}e^{-\tfrac{1}{2}(x-\mu_k)^T\mathbf{\Sigma}_k^{-1}(x-\mu_k)}. \end{align} $$ 可以得到线性判别函数 $$ \begin{align} \delta_k(x)=x^T\mathbf{\Sigma}^{-1}\mu_k - \tfrac{1}{2}\mu_k^T\mathbf{\Sigma}^{-1}\mu_k+\log \pi_k. \end{align} $$ 然后,最大的匹配值即为对应的分类 $$ \begin{align} \hat{G}(x)=\arg \max_{k}\delta_k(x). \end{align} $$
where $\alpha \in [0,1]$。另外,上式中的 $\hat{\mathbf{\Sigma}}$ 可以进一步 regularized 趋近于单位矩阵, $$ \begin{align} \hat{\mathbf{\Sigma}}(\gamma) = \gamma \hat{\mathbf{\Sigma}} + (1-\gamma)\hat{\sigma}^2\hat{\mathbf{I}}, \end{align} $$ for $\gamma \in [0,1]$.
目标是最大化 log-likelihood 函数 $$ \begin{align} \ell (\beta) = \sum_{i=1}^{N}\log \Pr\left(G=k|X=x;\beta\right) \end{align} $$ 最优化的 $\beta$ 可以通过 Newton-Raphson 算法求得 $$ \begin{align} \beta^{new} = \beta^{old} - \left(\frac{\partial^2 \ell (\beta)}{\partial \beta \partial \beta^T} \right)^{-1}\frac{\partial \ell (\beta)}{\partial \beta} \end{align} $$
其中 $M$ 是错误归类的点集合。算法更新如下 $$ \begin{align} \begin{pmatrix}\beta \\ \beta_0\end{pmatrix} \leftarrow \begin{pmatrix}\beta \\ \beta_0\end{pmatrix} + \rho \begin{pmatrix}y_ix_i \\ y_i\end{pmatrix}. \end{align} $$ where $\rho$ is the learning rate, equaling to 1 by default.
这里应该介绍下该算法什么时候可行,以及为什么失败。但是我暂时还没弄明白,需要参考其它资料,以后再补充吧。
学过统计的同学会说,LDA就是Fisher算法的特例,增加了正态分布,covariance一致$\mathbf{\Sigma}_1=\mathbf{\Sigma}_2=\cdots=\mathbf{\Sigma}$且满秩等条件。Fisher Linear Discriminant定义问题为
Find the linear combination $Z=a^T X$ such that the between-class variance is maximized relative to the within-class variance.
其数学表达形式为 $$ \begin{align} \max_{a}S=\frac{\delta^2_{between}}{\delta^2_{within}}=\frac{a^T\mathbf{B}a}{a^T\mathbf{W}a}, \end{align} $$ 其中,$\mathbf{B}$和 $\mathbf{W}$分别是 between-class variance 和 within-class variance。另外,$\mathbf{B}+\mathbf{W}=\mathbf{T}$ 就是整个输入 $X$ 的 total covariance matrix。Fisher 的最优化问题等价于 $$ \begin{align} \max_a a^T\mathbf{B}a \quad \textrm{ subject to } \quad a^T\mathbf{W}a=1. \end{align} $$ 最优解a为$\mathbf{W}^{-1}\mathbf{B}$的最大特征值.
当$x_k$的 mean 和 variance 分别是 $\mu_k$ 和 $\mathbf{\Sigma}_k$ 时,$Z=a^T X$ 的 mean 和 variance 分别是 $a^T\mu_k$ 和 $a^T\mathbf{\Sigma}_ka$。如果仅有两个变量的话,我们有(参考的wiki) $$ \begin{align} S=\frac{\delta^2_{between}}{\delta^2_{within}}=\frac{\left(a^T\mu_1 - a^T\mu_2\right)^2}{a^T\mathbf{\Sigma}_1a + a^T\mathbf{\Sigma}_2a}=\frac{\left(a^T(\mu_1 -\mu_2)\right)^2}{a^T(\mathbf{\Sigma}_1+\mathbf{\Sigma}_2)a}. \end{align} $$ 从而可以得到最优解 $$ \begin{align} a \propto (\mathbf{\Sigma}_1+\mathbf{\Sigma}_2)^{-1}(\mu_1 -\mu_2). \end{align} $$
对于通信背景的我来说,Fisher的思路其实是最大化 SNR。在通信系统中,输入变量 $x$ 和输出变量 $y$ 的对应关系为 $y=h x$,其中 $h$ 为信道,对应此处的 $a$ 映射。为了成功解调出 $y$ ,其实就是信号分类,对于输入变量$x_1$,$x_2$为干扰噪声,所以我们要使变量自身的 variace (即 $\mathbf{\Sigma}_k$ 或者 within-class variance $\mathbf{W}$) 尽量小,集中在均值点周围,这样$x_2$对应的$y_2$才不会跑到$y_1$的地盘;同时我们尽可能使$y_1$ 和 $y_2$ 相隔足够远,即 between-class variance $\mathbf{B}$ 最大化,这样在 $y_1$ 和 $y_2$ 之间划一条分割线可以使两个集合的交集最小,误差最小。因此,我们在设计通信系统 $h$ 时,尽力满足正交性,即与 $x_1-x_2$ 正交。
在 Fisher 的算法中,分子 $a^T\mathbf{B}a$ 比较好理解,为了更好的区分两个类,当然需要让两个类的中心相隔最远,最大化 between-class variance。而分母 $a^T\mathbf{W}a$ 的解释,见下图4.9。由于每个变量的 within-class variance 不容,变量分布在各个维度上会有倾向性(这个词是我的感官,可能不专业)。如果只考虑分子不除以分母,会得到图4.9左边的图,效果不好;除以分母后,得到右边的图。
我的理解,在除以 within-class variance 后,变量在各个维度上的倾向性被消除了,变成了 variance 为 $\mathbf{I}$ 的变量。即标准化了
LDA 包括 QDA 的计算先对角化其 covariance matrix,$\hat{\mathbf{\Sigma}}$ or $\hat{\mathbf{\Sigma}}_k$。
在实际计算过程中,先计算 $\mathbf{\Sigma}_k^T(x-\hat{\mu}_k)$ 部分。
由于 Linear Regression 简称为 LR,那么这个就叫 logit 吧,参考的其他博客。
以后验概率为例,为了保证 logit 之后的边界条件满足线性要求, $$ \begin{align} \log{\frac{\Pr\left(G=k|X=x\right)}{\Pr\left(G=K|X=x\right)}}=\beta_{k0}+\beta_k^Tx, \end{align} $$ 通常会对概率进行指数化处理,$Pr\left(G=k|X=x\right)=\exp(\beta_{k0}+\beta_k^Tx)$,然后归一化得到 $$ \begin{align} \Pr\left(G=k|X=x\right)=&\frac{\exp(\beta_{k0}+\beta_k^Tx)}{\sum_{l=1}^{K}\exp(\beta_{l0}+\beta_l^Tx)}. \end{align} $$ 如果任意选取一组为参照组,默认为第 $K$ 组,考虑归一化后得到 $$ \begin{align} \Pr\left(G=k|X=x\right)=&\frac{\exp(\beta_{k0}+\beta_k^Tx)}{1+\sum_{l=1}^{K-1}\exp(\beta_{l0}+\beta_l^Tx)}, \\ \Pr\left(G=K|X=x\right)=&\frac{1}{1+\sum_{l=1}^{K-1}\exp(\beta_{l0}+\beta_l^Tx)}. \end{align} $$ 目标是最大化 log-likelihood 函数 $$ \begin{align} \ell (\theta) = \sum_{i=1}^{N}\log p_{g_i}(x_i;\theta) \end{align} $$ where $p_k(x_i;\theta)=\Pr\left(G=k|X=x;\theta\right)$。
以 two-class 为例。为什么是 two?虽然思路是一样,算法还是可以实现,但是多了就太复杂了,没办法写出 expression。令当 $G_i = 1$ ($G_i = 2$) 时 $y_i = 1$ ($y_i = 0$)。同时,令 $p_1(x;\beta)= \frac{\exp{\beta^Tx}}{1+\exp{\beta^Tx}} = p(x;\beta)$,那么 $p_2(x;\beta)= 1- p(x;\beta)$。我们可以得到 log-likelihood 函数 $$ \begin{align} \ell{\beta} = \sum_{i=1}^{N}\left\{y_i\log p(x_i;\beta) + (1-y_i)\log (1-p(x_i;\beta)) \right\} = \sum_{i=1}^{N}\left\{y_i \beta^T x_i - \log (1+ e^{\beta^Tx_i}) \right\} \end{align} $$ 一阶导数,或者 score equation 为 $$ \begin{align} \frac{\partial \ell (\beta)}{\partial \beta} = \sum_{i=1}^{N}x_i(y_i-p(x_i;\beta)) =0 \end{align} $$ 二阶导数,或者 Hessian matrix 为 $$ \begin{align} \frac{\partial^2 \ell (\beta)}{\partial \beta \partial \beta^T} = -\sum_{i=1}^{N}x_i x_i^Tp(x_i;\beta)(1-p(x_i;\beta)) =0. \end{align} $$ 通过 Newton-Raphson 算法求得 $$ \begin{align} \beta^{new} = \beta^{old} - \left(\frac{\partial^2 \ell (\beta)}{\partial \beta \partial \beta^T} \right)^{-1}\frac{\partial \ell (\beta)}{\partial \beta} \end{align} $$
将以上推到写成矩阵形式, $$ \begin{align} \frac{\partial \ell (\beta)}{\partial \beta} = & \mathbf{X}^T(\mathbf{y} - \mathbf{p}) \\ \frac{\partial^2 \ell (\beta)}{\partial \beta \partial \beta^T} =& -\mathbf{X}^T\mathbf{W}\mathbf{X} \end{align} $$ 其中 $\mathbf{X}$ 是 $N\times(p+1)$ 维矩阵, $$ \begin{align} \mathbf{P} =& \left[p(x_0;\beta),\ p(x_1;\beta),\ p(x_2;\beta),\ \cdots,\ p(x_N;\beta) \right]^T \\ \mathbf{W} =& diag \left[p(x_0;\beta)(1-p(x_0;\beta)) \quad p(x_1;\beta)(1-p(x_1;\beta)) \quad p(x_2;\beta)(1-p(x_2;\beta))\quad \cdots \quad p(x_N;\beta)(1-p(x_N;\beta)) \right] \end{align} $$ 最后牛顿迭代可以写成 $$ \begin{align} \beta^{new} = \beta^{old} + (\mathbf{X}^T\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^T(\mathbf{y}-\mathbf{p}) = (\mathbf{X}^T\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}\mathbf{z} \end{align} $$ where $$ \begin{align} \mathbf{z} = \mathbf{X} \beta^{old} + \mathbf{W}^{-1}(\mathbf{y}-\mathbf{p}). \end{align} $$ 上面的迭代算法被称为 iteratively reweighted least squares (IRLS),因为每一步相当于求解一个 weighted least squares 问题: $$ \begin{align} \beta^{new}\leftarrow \arg \min_{\beta} (\mathbf{z} - \mathbf{X}\beta )^T\mathbf{W} (\mathbf{z} - \mathbf{X}\beta ). \end{align} $$
该方法依旧是基于线性条件,定义 Hyperplane or affine set $L$ 为 $$ \begin{align} f(x) = \beta_0 + \beta^T x = 0. \end{align} $$ 基于向量原理,任意一点 $x$ 到 $L$ 的有向距离为 $$ \begin{align} \frac{\beta^T}{\|\beta\|} (x-x_0) = \frac{1}{\|\beta\|} (\beta^T x+\beta_0) = \frac{1}{\|f'(x)\|}f(x) \end{align} $$ where $x_0$ is any point in $L$ and it satisfies $$ \begin{align} \beta^T x_0 = -\beta_0. \end{align} $$
Rosenblatt 在 1958 年提出的算法,让错误归类的点到 $L$ 的距离最小。该方法常用于神经网络分析,对滴,就是信号传输走最短距离。 $$ \begin{align} D(\beta, \beta_0) = - \sum_{i\in M}y_i(x_i^T\beta + \beta_0). \end{align} $$ 其中 $M$ 是错误归类的点集合。
分别对 $\beta$ 和 $\beta_0$ 求导,然后采用随机梯度下降法 stochastic gradient descent 求解。 $$ \begin{align} \partial \frac{D(\beta, \beta_0)}{\partial \beta} =& - \sum_{i\in M}y_ix_i\\ \partial \frac{D(\beta, \beta_0)}{\partial \beta_0} =& - \sum_{i\in M}y_i \end{align} $$ 算法更新如下 $$ \begin{align} \begin{pmatrix}\beta \\ \beta_0\end{pmatrix} \leftarrow \begin{pmatrix}\beta \\ \beta_0\end{pmatrix} + \rho \begin{pmatrix}y_ix_i \\ y_i\end{pmatrix}. \end{align} $$
该方法存在3个问题(Ripley, 1996):
解决了 PLA 算法答案不唯一的问题。最大化任意一点到 Hyperplane 的最小距离 $$ \begin{align} &\max_{\beta,\beta_0, \|\beta\|=1} M \\ \textrm{subject to } &y_i(x_i^T\beta+\beta_0) \geq M, \ \ i=1,\cdots,N. \end{align} $$ 由于 $\beta$ 的 scale 可以调整,问题等价于 $$ \begin{align} &\min_{\beta,\beta_0} \frac{1}{2}\|\beta\|^2 \\ \textrm{subject to } &y_i(x_i^T\beta+\beta_0) \geq 1, \ \ i=1,\cdots,N. \end{align} $$ 该 convex 优化问题对应的拉格朗日函数为, $$ \begin{align} L_p = \frac{1}{2}\|\beta\|^2 - \sum_{i=1}^{N} \alpha_i[y_i(x_i^T\beta+\beta_0) -1]. \end{align} $$ 对 $\beta$ 和 $\beta_0$求导可得 $$ \begin{align} \beta =& \sum_{i=1}^{N} \alpha_i y_ix_i \\ 0 =& \sum_{i=1}^{N} \alpha_i y_i \end{align} $$ 同时,由 KKT 条件得最优解满足 $$ \begin{align} \alpha_i[y_i(x_i^T\beta+\beta_0) -1] = 0. \end{align} $$ 要么 $\alpha_i = 0$;或者 $y_i(x_i^T\beta+\beta_0) =1$,即 $x_i$ 在 $L$ 的平行线上。我们在比较下 $\beta$ 的公式,有意思了,出去了所有 $\alpha_i = 0$ 的相式, $\beta$ 是由所有这些在 $L$ 的平行线上的 $x_i$ 决定的。这些 $x_i$ 被称为 support points。可以相信,这个概念和以后的 support vector machine 有关。
LDA,logit,和 PLA 都是线性分类算法,大家强调不同输入变量。
根据不同的分布,可以选择不同的分类算法。
如果预先不知道输入分布?我觉得 logit 更适合。因为 logit 和贝叶斯概率函数有关。即使我们不知道输入变量的分布,但是只要数据量足够大,一定符合大数定理,概率就准了。