Machine Learning第八周笔记：K means和降维
刚刚完成了Machine Learning第八周的课程，这一周主要介绍了K-means和降维，现将笔记整理在下面。
Unsupervised Learning
Clustering
Unsupervised Learning: Introduction
今天我们开始介绍无监督式学习（unsupervised learning）。下面两张图分别代表典型的有监督式学习和无监督式学习。一个典型的有监督式学习是从一个有标记的训练数据集出发（图中的两类数据点分别用圈圈和叉叉表示，圈圈一类，叉叉一类），目标是找到一个decision boundary能成功的将两类分开，并且得来的模型能够很好的对test set中的数据进行分类。而在无监督式学习中，训练数据集是没有标记的（图中的数据点全都用黑点表示），目标是找到一个算法能够发现数据集中的一些特点（图中使用聚类算法将数据点分成了两类，当然还有其他的算法）。聚类算法在科研和实际生产中有很多应用，如市场细分（market segmentation）、社交网络分析（social network analysis）、组织运算集群（organize computing clusters）、天文数据分析（Astronomical data analysis）等。
K-Means Algorithm
在聚类问题中，我们希望能有一个算法能够将给定数据集自动的分成不同的类或簇（同一个类或簇中的数据点存在某种相关）。K-means算法是目前使用最广泛的聚类算法，现在我们就来介绍什么是聚类算法以及它是怎样工作的。
我们希望将上图中的数据分成两类，下面使用K-means来完成这一工作。首先随机的选出两个点作为质心（cluster centroid），如图中两个叉号所示。K-means是一个迭代算法，在每一次迭代过程中主要完成两个工作：第一步是类分配（cluster assignment），第二步是移动质心（move centroid）。类分配是说遍历数据集中的每一个数据点，判断数据点到质心的距离，距离哪一个质心近就将数据点分配到哪一类，类分配结果见下图。移动质心是说根据前面刚刚完成的类分配，我们求出每一个类中数据点的中心点（这个点距离类中每一个点的平均距离或距离之和最小），将该点更新为该类新的质心，移动结果对比下面两张图。重复这样的迭代过程，直至数据集中每一个点的类别不再发生变化为止，此时我们能够得到两个类。
下面给出K-means的伪算法：
需要补充的是如果在一次类分配之后，某一个质心$\mu_i$对应的类别中没有一个数据点，我们就要将质心$\mu_i$删掉。
前面我们使用的都是很明显就可以分得开的数据（如下图左侧坐标系中的数据很明显可以分成3类），但是实际中处理的数据大多是直观上看起来不易分成多个类的数据（这种情况叫做non-separated clusters），如下图右侧中的数据。
Optimization Objective
我们用$\mu_k$来表示质心，共有K个质心。给定数据点$x^{(i)}$，我们用$c^{(i)}$表示该数据点所属的类别，用$\mu_{c^{(i)}}$表示该数据点所属类别对应的质心。下面给出了K-means的优化目标：
现在我们来分析一下迭代过程中的类分配和移动质心。类分配过程实际上就是一个在保持质心$\mu_1$到$\mu_K$不变的前提下最小化损失函数，而移动质心过程则是在保持$c^{(1)}$到$c^{(m)}$不变的前提下去最小化损失函数。
Random Initialization
如何初始化K个质心呢？首先，K应该小于数据集的大小m，否则不符合常理。在数据集中随机的选出K个数据点作为初始的质心。基于不同的初始质心，K-means会有不同的收敛结果。下图给出了同一组数据的3个不同结果，很容易看出只有右上是全局最优（global optimal），而另外两个都是局部最优。为了得到好的分类结果，我们就需要设置多组初始质心，然后在多组结果中选择最优的一组。
当K较小（2~10）时，由不同的初始质心得到的结果可能相差很大；当K较大时，由不同质心得到的结果差别可能就不是特别明显了。
Choosing the Number of Clusters
无监督式学习中数据没有标记，所以也就没有一个明确的答案（a clear cut answer）。一个数据集分成几类比较合适是比较模糊的，下图中的数据集可以分成2、3或者4类都是可以的。
定量的分析损失函数随K值得变化情况，可以得到下面的图像，可以看到图像好像人的胳膊，开始一段迅速下降，之后变化缓慢，这被称作elbow method。至于到底怎样选择K，还要结合实际的应用。
Dimensionality Reduction
这一节介绍数据的降维。
Motivation
Motivation Ⅰ: Data Compression
数据压缩（data compression）不仅能减少数据的存储量，节省了内存和磁盘，还会大大加快学习算法的运行速度。我们首先来介绍降维（dimensionality reduction）。假设我们已经收集了大量数据，其中含有特别多的特征，这里我们只选择其中的两个特征$x_1$和$x_2$，二者都用来表示某一事物的长度，$x_1$的单位是厘米，$x_2$的单位是英寸，如下图所示。很明显数据存在冗余，我们可以使用一个特征来表示事物的长度，数据从二维降到一维。当然，实际生产中我们可能不会出现这么低级的问题，这里仅以此来做说明。工业上的一个问题中可能会有上千个特征，很容易出现冗余。联系到直升机飞行员的例子，这里用$x_1$表示飞行员的驾驶技巧，用$x_2$表示飞行员对工作的喜爱程度，这两个特征时高度相关的，我们要从这两个特征中提取出一个新的特征$z_1$（如pilot attitude）来代替它们。如下图所示，我们找到这样一条直线，每一个数据点在这条直线上的投影作为新的数据点，我们就可以实现2维（原来的数据点是二维的）到一维的降维（原来的数据点与原点的连线在这条直线上的投影的长度是一个数字，是一维的）。通过类似这样的方法，我们就实现了降维。
再给出一个三维降到二维的例子。如下图所示，我们可以将三维坐标系中的数据点投影到一个二维平面上，实现三维到二维的降维。
Motivation Ⅱ: Visualization
在很多学习算法的应用中，对数据的理解可以帮助我们更好的改善模型，而降维恰好可以帮助我们实现数据的可视化，增进我们队数据的理解。当然，降维还有其他的许多作用。下图中列出了多个国家的多个指标，我们可以将这多个指标叫做衡量一个国家的个特征。衡量一个国家的特征可能会有几十个之多，我们如何去直观的理解这些数字呢，这里可能就需要将数据可视化。我们设法将下面表格中的多个特征降维到2个特征（如下面第二图所示），然后在直角坐标系中表示出它们。
上图中坐标系的横轴可以表示一个国家的经济总量，纵轴代表人均经济量，图中的每一个点代表一个国家。如右上侧的那个绿点可能是美国——经济总量世界第一，人均经济量也达到很高的一个数值；最上方的那个绿点可能代表新加坡——经济总量不大（国家太小），但人均量很高。了解了一个国家在这张图上的位置，我们也就知道了这个国家的经济状况。当然我们也可以将横纵表座变成其他的指标，实现对其他特征的可视化。
Principal Component Analysis
Principal Component Analysis Problem Formulation
这一节介绍主成分分析（Principal Component Analysis，PCA）。下图中叉号代表数据点，我们先给出了一条很好的直线（图中红线），通过将数据点投影到直线上来实现降维。为什么说这样的一条直线是一个很好的选择呢，下面给出答案。图中用蓝色的线段画出了数据点到投影的距离，称为投影误差（projection error），我们可以发现这些蓝色线段都是很短的。PCA的目标就是找到一个lower dimensional surface(对于二维的数据点来说，lower dimensional surface就是一条直线）并将数据点投影到上面，使得这些蓝色线段的长度之和尽可能小。作为比较，图中给出了另外一个lower dimensional surface（品红色线），我们可以看到数据点到其投影的距离都很大，所以这不会是一个好的选择。
对于二维降到一维，PCA的目标是找到一个二维向量（有了这个二维向量，也就有了上图中的那条红线，等价的），并且最小化projection error。对于从n维降到k维，我们则需要找到k个n维向量来实现。下一小节会介绍PCA的工作流程。
从前面二维降到一维的例子中，我们发现PCA有些像线性回归，下面我们就指出二者的不同。下图中的左右两个坐标系分别表示线性回归和PCA。两个算法的目的是找到一条直线，使得图中蓝色线段之和最小。但要注意的是，线性回归中的蓝色线段是数据点在y轴方向上到直线的距离，表示hypothesis的预测值与实际值的误差；而PCA中的蓝色线段是数据点到直线的垂直距离，可以理解为新的低维特征与原来的高维特征之间的信息差。（我是这样理解的，这个距离越小，说明在降维的过程中保留下来的信息越多，扭曲的信息越少，最后得到的hypothesis越可靠）这两个距离的含义截然不同，所以PCA和线性回归没有什么关系。
Principal Component Analysis Algorithm
这一节给出PCA的工作流程。
- 求解每个特征的均值 $$\mu_j=\frac{1}{m}\sum_{i=1}^{m}x_j^{(i)}$$
- 对数据进行标准化，即对每个$x_j^{(i)}$有（注意，图中还有一个$s_j$，没听明白怎么来的）:
$$x_j^{(i)}=x_j^{(i)}-\mu_j$$ - 求解协方差矩阵：
$$\sum=\frac{1}{m}\sum_{i=1}^{n}(x^{(i)})(x^{(i)})^T$$ - 根据svd求解特征值和特征向量：
$$[U, S, V] = svd(\sum)$$ - 将特征向量按对应特征值大小从左到右按列排列成矩阵U，取前k列组成矩阵$U_{reduce}$
- $Z=U_{reduce}^TX$即为降维到k维后的数据
以上是PCA的工作流程，我会在最后的补充中给出一些数学上的推导。
Applying PCA
Reconstruction from Compressed Representation
前面介绍了使用PCA对数据进行降维，已达到压缩数据提高计算速率的目标，现在我们尝试着从压缩后的数据出发还原数据，注意这里得到的数据是原始数据的近似。我们沿用前面二维降到一维的例子，现在我们尝试着从一维升维到二维。具体的公式是：
$$X_{approx}=U_{reduce}Z$$
如果project error能够取到很小很小的值，则会有：
$$X\approx X_{approx}$$
Choosing the Number of Principal Components
这一小节给出选择K值的一个指标。首先，我们给出一个计算公式：
$$err=\frac{\frac{1}{m}\sum_{i=1}^{m}\left\|x^{(i)}-x_{approx}^{(i)}\right\|^2}{\frac{1}{m}\sum_{i=1}^{m}\left\|x^{(i)}\right\|^2}$$
我们用这个公式来量化原始数据与降维后数据的差距。若该值不超过0.01，我们就可以说“99% of variance is retained”。我们可以把err作为选择k值得指标：取k=1，完成一系列降维计算，判断err是否满足要求，若满足则确定该k值为最终的k；若不满足，则k加1，继续前面的计算直至满足要求。
这里给出了计算err的一种简便方法。PCV求解过程中有这样一个公式
$$[U, S, V] = svd(\sum)$$
其中的S是一个对角矩阵，对角线上式要求特征值，从上到下依次递减。err可由下面的公式求得：
$$err=\frac{\sum_{i=1}^{k}S_{ii}}{\sum_{i=1}^{n}S_{ii}}$$
Advice for Applying PCA
假设一条数据中包含10000个特征，通过PCA降低到1000个特征，无论是在有监督式学习还是在无监督式学习中，这种降维对学习速度的提升是显著的。需要注意的一点是PCA只能在training set上进行，得到的$U_{reduce}$同样能应用在cross validation set和test set上。
在前面的学习中，我们知道过多的特征提取会导致overfitting，此时采取的措施是适当的去掉一些特征。而现在的降维恰好减小了特征的数量，那么我们是否可以通过降维来改善overfitting呢？答案是否定的。（个人理解，这种降维仍然是在尽量保持原有数据信息的前提下进行的，并不能达到前面直接去掉某些特征的效果。）
需要注意的一点是，在进行机器学习时，我们首先使用的是原始数据。只有当算法的速度过慢、内存或磁盘不够，或者你有足够的把握确定问题出在原始数据上（此处的问题当然不是指数据的真假问题），这时你你才可以考虑降维。
Supplement
查了一下Pattern Recognition And Machine Learning，先将有关PCA的内容摘抄在下面：
There are two commonly used definitions of PCA that give rise to the same algorithm. PCA can be defined as the orthogonal projection of the data onto a lower dimensional linear space, known as the principal subspace, such that the variance of the projected data is maximized(Hotelling, 1933). Equivalently, it can be defined as the linear projection that minimizes the average projection cost, defined as the mean squared distance between the data points and their projections(Pearson, 1901).
PCA有两个不同的定义，两个定义最终的算法是同一个。PCA可以定义成数据在一个低维线性空间上的垂直投影，这个空间称作主元子空间（principal subspace），这样得到的投影数据的方差达到最大值（Hotelling, 1993）。等价地，PCA可以定义成最小化投影成本的线性投影，还可以定义成数据点与其投影的平方距离的平均。
这里给出书中的第一个定义最大方差公式（maximum variance formulation）的解释，以后有时间再给出第二个。
先是摘抄，后是翻译。
Consider a data set of observation {$x_n$} where n = 1,…, N, and $x_n$ is a Euclidean variable with dimensionality D. Our goal is to project the data onto a space having dimensionality M < D while maximizing the variance of projected data.
To begin with, consider the projection on a one-dimensional space (M=1). We can define the direction of this space using a D-dimensional vector $u_1$, which for convenience (and without loss generality) we shall choose to be a unit vector so that $u_1^Tu_1=1$ (note that we only interested in the direction defined by $u_1$, not in the magnitude of $u_1$ itself). Each data point $x_n$ is then projected onto a scalar value $u_1Tx_n$. The mean of the projected data is $u_1^T\bar{x}$ where $\bar{x}$ is the sample set mean given by$$\bar{x}=\frac{1}{N}\sum_{n=1}^{N}x_n$$
and the variance of the projected data is given by
$$\frac{1}{N}\sum_{n=1}^{N}(u_1^Tx_n - u_1^T\bar{x})^2 = u_1^TSu_1$$
where S is the data covariance matrix defined by
$$S=\frac{1}{N}\sum_{n=1}^{N}(x_n - \bar{x})(x_n - \bar{x})^T$$.
We now maximize the projected variance $u_1^TSu_1$ with respect to $u_1$. Clearly, this has to be a constrained maximization to prevent $\left|u_1\right| \xrightarrow{} \infty$. The appropriate constraint comes from the normoalization condition $u_1^Tu_1 = 1$. To enforce this constraint, we introduce a Lagrange multiplier that we shall denote by $\lambda_1$, and then make an unconstrained maximization of
$$u_1^TSu_1 + \lambda_1(1 - u_1^Tu_1)$$
By setting the derivative with respect to u_1 equal to zero, we see that this quantity will have a stationary point when
$$Su_1 = \lambda_1u_1$$
which says that $u_1$ must be an eigenvector of S. If we left-multiply by $u_1^T$ and make use of $u_1^Tu_1 = 1$, we see that variance is given by
$$u_1^Tu_1 = \lambda_1$$
and so the variance will be a maximum when we set $u_1$ equal to the eigenevector having the largest eigenvalue $\lambda_1$. The eigenvector is known as the first principal component.
We can define additional principal components in an incremental fashion by choosing each new direction to be that which maximizes the projected variance amongst all possible directions orthogonal to those already considered. If we consider the general case of an M-dimensional projection space, the optimal linear projection for which the variance of the projected data is maximized is now defined by the M eigenvectors $u_1$,…, $u_M$ of the data covariance matrix S corresponding to the M largest eigenvalues $\lambda_1$,…, $\lambda_M$. This is easily shown using proof by induction.
To summarize, principal component analysis involves evaluating the mean $\bar{x}$ and the covariance matrix S of the data set and then finding the M eigenvectors of S corresponding to the M largest eigenvalues. Algorithms for finding the eigenvectors and eigenvalues, as well as additional theorems related to eigenvector decomposition, can be found in Golub and Van Loan (1996). Note that the computational cost of computing the full eigenvector decomposition for a matrix of size $ D \times D$ is $O(D^3)$.
If we plan to project our data onto the first M principal components, then we only need to find the first M eigenvalues and eigenvectors. This can be done with more efficient techniques, such as the power method (Golub and Van Loan, 1996), that scale like $O(MD^3)$, or alternatively we can make sure of EM algorithm.
考虑数据集{$x_n$}，其中n = 1,…,N，且$x_n$是D维欧几里得向量。我们的目标是将数据投影到一个M（< D）维空间上使得投影数据的方差最大。
首先，我们来考虑在一维空间上的投影。我们可以用一个D维向量$u_1$来定义这个空间的方向向量，为求方便且不失一般性，我们可以选择一个单位向量，故有$u_1^Tu_1=1$。（注意，我们的重点在由$u_1$定义的这个方向向量，而不是$u_1$本身）每个数据点$x_n$投影到上面得到一个标量$u_1^Tu_1=1$。投影数据的均值为$u_1^T\bar{x}$，其中$\bar{x}$由样本均值给出： $$\bar{x}=\frac{1}{N}\sum_{n=1}^{N}x_n$$而投影数据的方差则由下式给出：
$$\frac{1}{N}\sum_{n=1}^{N}(u_1^Tx_n - u_1^T\bar{x})^2 = u_1^TSu_1$$其中S是由下式定义的协方差矩阵：
$$S=\frac{1}{N}\sum_{n=1}^{N}(x_n - \bar{x})(x_n - \bar{x})^T$$我们现在来最大化投影方差$u_1^TSu_1$，其与$u_1$有关。很显然，这是约束最优问题。引入Lagrange因子$\lambda_1$，求解下式的最大值：
$$u_1^TSu_1 + \lambda_1(1 - u_1^Tu_1)$$
令上式对$u_1$的求导等于零，我们会得到
$$Su_1 = \lambda_1u_1$$
这说明$u_1$必须是S的一个特征向量。等式两边同时左乘$u_1^T$，又因为$u_1^Tu_1 = 1$，此时我们发现投影数据的方差变成了下面这样：
$$u_1^Tu_1 = \lambda_1$$
所以当我们将$u_1$设置成最大特征值对应的特征向量时，我们就可以得到投影数据方差的最大值。（好神奇，是不是！！！）这个特征向量被称作第一主成分（the first principal component）。
我们可以以一种递增的顺序选择每一个新的方向向量，使得投影方差最大化，且每一个新的方向向量都与已选择了的正交。考虑M维投影空间的一般情况，最佳的投影线性空间则由协方差矩阵S的前M大的特征值$\lambda_1$,…, $\lambda_M$对应的特征向量$u_1$,…,$u_M$构成，以使得投影数据的方差最大。（这通过推导很容易证明，还不懂的可以去查高等代数，虽然我也快忘光了。）
如果我们计划将数据投影到前M个主成分上，那我们只需要找到前M大的特征值及其特征向量。注意对一个$ D \times D$的矩阵进行特征分解的算法，其时间复杂度为$O(D^3)$。
PCA的算法流程也包括在上面的解释当中了，一举两得。
另外，在网上查阅资料时，找到一篇相关博文，感觉不错，见PCA数学原理。
想要了解更多，请查阅Pattern Recognition And Machine Learning等其他资料。