图嵌入（拉普拉斯特征映射Laplacian Eigenmaps）

图嵌入（Graph embedding）的意义

Graph广泛存在于真实世界的多种场景中，即节点和边的集合。比如社交网络中人与人之间的联系，生物中蛋白质相互作用以及通信网络中的IP地址之间的通信等等。除此之外，我们最常见的一张图片、一个句子也可以抽象地看做是一个图模型的结构，图结构可以说是无处不在。通过对它们的分析，我们可以深入了解社会结构、语言和不同的交流模式，因此图一直是学界研究的热点。
图分析任务可以大致抽象为以下四类:
( a )节点分类：节点分类旨在基于其他标记的节点和网络拓扑来确定节点的标签(也称为顶点)
( b )链接预测：预测缺失链路或未来可能出现的链路的任务
( c )聚类：用于发现相似节点的子集，并将它们分组在一起
( d )可视化：有助于深入了解网络结构

为什么要使用图形嵌入

需要使用图形嵌入有以下几个原因:

1.机器学习图形是有限的。图由边和节点组成。这些网络关系只能使用数学、统计和机器学习的特定子集，而向量空间有更丰富的方法工具集。

2.嵌入是压缩的表示。
邻接矩阵描述图中节点之间的连接。它是一个 $\times|V|$ 矩阵 ( $∣ V ∣$ 是图中节点的个数)。矩阵中的每一列和每一行表示一个节点。矩阵中的非零值表示两个节点相连。但是使用邻接矩阵作为大型图的特征空间几乎是不可能的。

3.假设一个图有 $M$ 个节点和一个 $\times M$ 的邻接矩阵。嵌入比邻接矩阵更实用，因为它们将节点属性打包到一个维度更小的向量中。向量运算比图形上的可比运算更简单、更快。

4.真实的图（网络）往往是高维、难以处理的。

20世纪初，研究人员发明了图形嵌入算法，作为降维技术的一部分。

他们首先根据实际问题构造一个 $D$ 维空间中的图。然后将图的节点嵌入到 $d （ d < < D ）$ 维向量空间中。嵌入的思想是在向量空间中保持连接的节点彼此靠近。拉普拉斯特征映射（Laplacian Eigenmaps）和局部线性嵌入（Locally Linear Embedding ，LLE）是基于这一原理的算法的例子。然而，可伸缩性是这种方法的一个主要问题，它的时间复杂度是 $O (| V| ^2)$ 。

自2010年以来，关于图嵌入的研究已经转移到解决网络稀疏性的可伸缩图嵌入技术上。
例如，图分解（Graph Factorization）使用邻接矩阵的近似分解作为嵌入。LINE扩展了这种方法，并试图保持一阶和二阶近似。HOPE通过使用广义奇异值分解( SVD )分解相似性矩阵而不是邻接矩阵来扩展LINE以试图保持高阶邻近性。SDNE 使用自动编码器嵌入图形节点并捕捉高度非线性的依赖关系。这些新的可扩展方法的时间复杂度为 $O (∣ E ∣)$ 。

什么是图嵌入

通俗点讲：图嵌入（Graph Embedding，也叫Network Embedding）是一种将图数据（通常为高维稠密的矩阵）映射为低微稠密向量的过程，能够很好地解决图数据难以高效输入机器学习算法的问题。

节点的分布式表示
节点之间的相似性表示链接强度
编码网络信息并生成节点表示

图嵌入是将属性图转换为向量或向量集。嵌入应该捕获图的拓扑结构、顶点到顶点的关系以及关于图、子图和顶点的其他相关信息。更多的属性嵌入编码可以在以后的任务中获得更好的结果。我们可以将嵌入式大致分为两类:

顶点嵌入: 每个顶点(节点)用其自身的向量表示进行编码。这种嵌入一般用于在顶点层次上执行可视化或预测。比如，在2D平面上显示顶点，或者基于顶点相似性预测新的连接。
图嵌入: 用单个向量表示整个图。这种嵌入用于在图的层次上做出预测，或者想要比较或可视化整个图。例如，比较化学结构。

嵌入的挑战

如前所述，图嵌入的目标是发现高维图的低维向量表示，而获取图中每个节点的向量表示是十分困难的，并且具有几个挑战，这些挑战一直在推动本领域的研究：

属性选择：节点的“良好”向量表示应保留图的结构和单个节点之间的连接。第一个挑战是选择嵌入应该保留的图形属性。考虑到图中所定义的距离度量和属性过多，这种选择可能很困难，性能可能取决于实际的应用场景。它们需要表示图拓扑、节点连接和节点邻居。预测或可视化的性能取决于嵌入的质量。
可扩展性：大多数真实网络都很大，包含大量节点和边。嵌入方法应具有可扩展性，能够处理大型图。定义一个可扩展的模型具有挑战性，尤其是当该模型旨在保持网络的全局属性时。好的嵌入方法需要在大型图上高效。
嵌入的维数：实际嵌入时很难找到表示的最佳维数。
较高的维数可能会提高重建精度，但具有较高的时间和空间复杂性。
较低的维度虽然时间、空间复杂度低，但无疑会损失很多图中原有的信息。用户需要根据需求做出权衡。在一些文章中，他们通常报告说嵌入大小在128到256之间对于大多数任务来说已经足够了。在Word2vec方法中，他们选择了嵌入长度300。

图嵌入的方法

在过去的十年里，在图形嵌入领域已经有了大量的研究，重点是设计新的嵌入算法。发展到现在，大体上可以将这些嵌入方法分为三大类：
( 1 )基于因子分解的方法
( 2 )基于随机游走的方法
( 3 )基于深度学习的方法。

这里主要介绍基于因子分解的方法的Laplacian Eigenmaps （拉普拉斯的特征映射）
在介绍之前先说明一下背景知识

谱聚类

谱聚类（spectral clustering）是广泛使用的聚类算法，比起传统的K-Means算法，谱聚类对数据分布的适应性更强，聚类效果也很优秀，同时聚类的计算量也小很多，更加难能可贵的是实现起来也不复杂。在处理实际的聚类问题时，个人认为谱聚类是应该首先考虑的几种算法之一。下面我们就对谱聚类的算法原理做一个总结。

1. 谱聚类概述

谱聚类是从图论中演化出来的算法，后来在聚类中得到了广泛的应用。它的主要思想是把所有的数据看做空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。

乍一看，这个算法原理的确简单，但是要完全理解这个算法的话，需要对图论中的无向图，线性代数和矩阵分析都有一定的了解。下面我们就从这些需要的基础知识开始，一步步学习谱聚类。

2. 谱聚类基础之一：无向权重图

由于谱聚类是基于图论的，因此我们首先温习下图的概念。对于一个图 $G$ ，我们一般用点的集合 $V$ 和边的集合 $E$ 来描述。即为 $G (V, E)$ 。
其中 $V$ 即为我们数据集里面所有的点 $v_1,v_2,...,v_n)$ 。对于 $V$ 中的任意两个点，可以有边连接，也可以没有边连接。定义权重 $w_{ij}$ 为点 $v_i$ 和点 $v_j$ 之间的权重。由于是无向图，所以 $w_{ij}=w_{ji}$ 。

对于有边连接的两个点 $v_i$ 和 $v_j$ ， $w_{ij}>0$ ,对于没有边连接的两个点 $v_i$ 和 $v_j$ ， $w_{ij}>0$ 。对于图中的任意一个点 $v_i$ ，它的度 $d_i$ 定义为和它相连的所有边的权重之和，即

$di=∑j=1nwijd_i = \sum_{j=1}^n w_{ij}$

利用每个点度的定义，我们可以得到一个 $n×nn\times n$ 的度矩阵 $D$ ,它是一个对角矩阵，只有主对角线有值，对应第i行的第 $i$ 个点的度数，定义如下：
$\begin{bmatrix} d_i & \cdots & \cdots \\ \cdots & d_2 & \cdots \\ \vdots & \vdots & \ddots \\ \cdots & \cdots & d_n \end{bmatrix}$
利用所有点之间的权重值，我们可以得到图的邻接矩阵 $W$ ，它也是一个 $\times n$ 的矩阵。
第 $i$ 行的第 $j$ 个值对应我们的权重 $w_{ij}$ 。

除此之外，对于点集 $V$ 的的一个子集 $A \subset V$ ，我们定义：
$∣ A ∣ : = 子集 A 中点的个数$
$\sum_{i \in A} d_i$

3. 谱聚类基础之二：相似矩阵

在上一节我们讲到了邻接矩阵 $W$ ，它是由任意两点之间的权重值 $w_{ij}$ 组成的矩阵。通常我们可以自己输入权重，但是在谱聚类中，我们只有数据点的定义，并没有直接给出这个邻接矩阵，那么怎么得到这个邻接矩阵呢？

基本思想是，距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，不过这仅仅是定性，我们需要定量的权重值。一般来说，我们可以通过样本点距离度量的相似矩阵 $S$ 来获得邻接矩阵 $W$ 。

构建邻接矩阵W的方法有三类。ϵ-邻近法，K邻近法和全连接法。

对于ϵ-邻近法，它设置了一个距离阈值ϵ，然后用欧式距离sij度量任意两点xi和xj的距离。即相似矩阵的 $s_ij=||x_i−x_j||_{2}^2$ , 然后根据 $s_{ij}$ 和ϵ的大小关系，来定义邻接矩阵 $W$ 如下：

$wij={0sij>ϵϵsij≤ϵw_{ij} = \begin{cases} 0 & s_{ij} \gt ϵ \\ ϵ & s_{ij} \leq ϵ \end{cases}$
从上式可见，两点间的权重要不就是ϵ,要不就是0，没有其他的信息了。距离远近度量很不精确，因此在实际应用中，我们很少使用ϵ-邻近法。

第二种定义邻接矩阵 $W$ 的方法是 $K$ 邻近法，利用 $K N N$ 算法遍历所有的样本点，取每个样本最近的k个点作为近邻，只有和样本距离最近的k个点之间的 $wij>0w_{ij} \gt 0$ 。但是这种方法会造成重构之后的邻接矩阵 $W$ 非对称，我们后面的算法需要对称邻接矩阵。为了解决这种问题，一般采取下面两种方法之一：
第一种K邻近法是只要一个点在另一个点的K近邻中，则保留 $s_{ij}$

$wij=wji={0xi∉KNN(xj)andxj∉KNN(xi)exp(−∣∣xi−xj∣∣222σ2xi∈KNN(xj)orxj∈KNN(xi)w_{ij}= w_{ji} = \begin{cases} 0 & x_i \not\in KNN(x_j) & and & x_j \not\in KNN(x_i) \\ exp(-\frac{||x_i - x_ j ||_{2}^2} {2 \sigma^2} & x_i \in KNN(x_j) & or & x_j \in KNN(x_i) \end{cases}$
第二种K邻近法是必须两个点互为K近邻中，才能保留 $s_{ij}$
$wij=wji={0xi∉KNN(xj)orxj∉KNN(xi)exp(−∣∣xi−xj∣∣222σ2)xi∈KNN(xj)andxj∈KNN(xi)w_{ij}= w_{ji} = \begin{cases} 0 & x_i \not\in KNN(x_j) & or & x_j \not\in KNN(x_i) \\ exp(-\frac{||x_i - x_ j ||_{2}^2} {2 \sigma^2}) & x_i \in KNN(x_j) & and & x_j \in KNN(x_i) \end{cases}$
第三种定义邻接矩阵 $W$ 的方法是全连接法，相比前两种方法，第三种方法所有的点之间的权重值都大于0，因此称之为全连接法。可以选择不同的核函数来定义边权重，常用的有多项式核函数，高斯核函数和Sigmoid核函数。最常用的是高斯核函数RBF，此时相似矩阵和邻接矩阵相同：
$wij=sij=exp(−∣∣xi−xj∣∣222σ2)w_{ij}=s_{ij}= exp(-\frac{||x_i - x_ j ||_{2}^2} {2 \sigma^2})$
在实际的应用中，使用第三种全连接法来建立邻接矩阵是最普遍的，而在全连接法中使用高斯径向核RBF是最普遍的。

4. 谱聚类基础之三：拉普拉斯矩阵

单独把拉普拉斯矩阵(Graph Laplacians)拿出来介绍是因为后面的算法和这个矩阵的性质息息相关。它的定义很简单，拉普拉斯矩阵 $L = D - W$ 。 $D$ 即为我们第二节讲的度矩阵，它是一个对角矩阵。而 $W$ 即为我们第二节讲的邻接矩阵，它可以由我们第三节的方法构建出。

拉普拉斯矩阵有一些很好的性质如下：

1）拉普拉斯矩阵是对称矩阵，这可以由 $D$ 和 $W$ 都是对称矩阵而得。

2）由于拉普拉斯矩阵是对称矩阵，则它的所有的特征值都是实数。

3）对于任意的向量 $f$ ,我们有
$fTLf=12∑i,j=1nwij(fi−fj)2f^{T} L f = \frac{1}{2} \sum_{i,j=1}^{n} w_{ij}(f_i - f_j)^2$
这个利用拉普拉斯矩阵的定义很容易得到如下：

$fTLf=fTDf−fTWf=∑i=1ndifi2−∑i,j=1nwijfifj=12(∑i=1ndifi2−2∑i,j=1nwijfifj+∑j=1ndjfj2)由于：di=∑j=1nwij=12(∑i,j=1nwijfi2−2∑i,j=1nwijfifj+∑j,i=1nwijfj2)=12∑i,j=1nwij(fi−fj)2\begin{aligned} f^{T}Lf & = f^{T}Df − f^{T}Wf \\ &= \sum_{i=1}^n d_if_{i}^2 - \sum_{i,j=1}^n w_{ij}f_if_j \\ &= \frac{1}{2} ( \sum_{i=1}^n d_if_{i}^2 - 2\sum_{i,j=1}^n w_{ij}f_if_j + \sum_{j=1}^n d_jf_{j}^2) \\ 由于：d_i = \sum_{j = 1}^{n} w_{ij} \\ &= \frac{1}{2} ( \sum_{i,j=1}^n w_{ij}f_{i}^2 - 2\sum_{i,j=1}^n w_{ij}f_if_j + \sum_{j,i=1}^n w_{ij}f_{j}^2) \\ &=\frac{1}{2} \sum_{i,j=1}^{n} w_{ij}(f_i - f_j)^2 \end{aligned}$
4) 拉普拉斯矩阵是半正定的，且对应的n个实数特征值都大于等于0，即 $0=λ_1≤λ_2≤...≤λ_n$ ，且最小的特征值为0，这个由性质3很容易得出。

5. 谱聚类基础之四：无向图切图

对于无向图G的切图，我们的目标是将图 $G (V, E)$ 切成相互没有连接的k个子图，每个子图点的集合为： $A_1,A_2,..,A_k$ ，它们满足 $A_i∩A_j=∅$ ,且 $A_1∪A_2∪...∪A_k=V$ .

对于任意两个子图点的集合 $A, B \subset V$ ， $A \cap B = \emptyset$ ,

$我们定义A和B之间的切图权重为：W(A,B)=∑i∈A,j∈Bwij我们定义A和B之间的切图权重为：W(A,B)=\sum_{i \in A, j \in B} w_{ij}$
那么对于我们 $k$ 个子图点的集合： $A1,A2,.⋯,AkA_1,A_2,. \cdots,A_k$ ，我们定义切图cut为：

$cut(A1,A2,.⋯,Ak)=12∑i=1kW(Ai,Ai‾)cut(A_1,A_2,. \cdots,A_k)=\frac{1}{2} \sum_{i=1}^k W(A_i, \overline{A_i})$
其中 $Ai‾\overline{A_i}$ 为 $A_i$ 的补集，意为除 $A_i$ 子集外其他 $V$ 的子集的并集。

那么切图可以让子图内的点权重高，子图间的点权重低呢？一个自然的想法就是最小化 $cut(A1,A2,.⋯,Ak)cut(A_1,A_2,. \cdots,A_k)$ ，但是最小化这个只是达到了子图间的点权重低的效果。

如图：

我们选择一个权重最小的边缘的点，比如 $C$ 和 $H$ 之间进行 $c u t$ ，这样可以最小化 $cut(A1,A2,.⋯,Ak)cut(A_1,A_2,. \cdots,A_k)$ ，但是却不是最优的切图（明显在单独的左下角子图中，只有一个点，必然子图内的点权重就小），如何避免这种切图，并且找到类似图中"Best Cut"这样的最优切图呢？我们下一节就来看看谱聚类使用的切图方法。

6. 谱聚类之切图聚类

为了避免最小切图导致的切图效果不佳，我们需要对每个子图的规模做出限定，一般来说，有两种切图方式，第一种是RatioCut，第二种是Ncut。下面我们分别加以介绍。

6.1 RatioCut切图
在数学模型中能表示：两个之间相反变化（一个变大一个变小）同时还能达到同样的效果的，其中一个就是分数表示形式：
$ab\frac{a}{b}$
当a变小，b变大的时候，整个数学形式是在变小。

RatioCut切图为了避免第五节的最小切图，对每个切图，不光考虑最小化 $cut(A1,A2,.⋯,Ak)cut(A_1,A_2,. \cdots,A_k)$ ，它还同时考虑最大化每个子图点的个数。即：
$RatioCut(A1,A2,.⋯,Ak)=12∑i=1kW(Ai,Ai‾)∣Ai∣RatioCut(A_1,A_2,. \cdots,A_k) = \frac{1}{2} \sum_{i=1}^k \frac{W(A_i, \overline{A_i})}{|A_i|}$
那么怎么最小化这个RatioCut函数呢？牛人们发现，RatioCut函数可以通过如下方式表示。

我们引入指示向量 $hj∈{h1,h2,..hk},j=1,2,...k,h_j \in \{h_1,h_2,..h_k\}, j=1,2,...k,$ 对于任意一个向量 $h_j$ , 它是一个n维向量（n为样本数），我们定义 $h_{ij}$ 为：
$hij={0xi∉Aj1∣Aj∣xi∈Ajh_{ij}= \begin{cases} 0 & x_i \not\in A_j \\ \frac{1} {\sqrt{|A_j|}} & x_i \in A_j \end{cases}$

1xi∈Ajxi∈Aj

那么对于 $h_{i}^{T} Lh_i$ 就有：
$hiTLhi=12∑m=1∑n=1wmn(fim−fjn)2=12∑m∈Ai,n∉Aiwmn(1∣Ai∣−0)2+∑m∉Ai,n∈Aiwmn(0−1∣Ai∣)2=12∑m∈Ai,n∉Aiwmn1∣Ai∣+∑m∉Ai,n∈Aiwmn1∣Ai∣=12(cut(Ai,Ai‾)1∣Ai∣+cut(Ai‾,Ai)1∣Ai∣)=cut(Ai,Ai‾)∣Ai∣\begin{aligned} h_{i}^{T} Lh_i & = \frac{1}{2} \sum_{m=1}\sum_{n=1} w_{mn}(f_{im}- f_{jn})^2 \\ &= \frac{1}{2} \sum_{m\in A_i,n \not\in A_i} w_{mn}(\frac{1} {\sqrt{|A_i|}} - 0)^2 + \sum_{m \not\in A_i,n \in A_i} w_{mn}( 0 - \frac{1} {\sqrt{|A_i|}} )^2\\ & = \frac{1}{2} \sum_{m\in A_i,n \not\in A_i} w_{mn}\frac{1} {|A_i|} + \sum_{m \not\in A_i,n \in A_i} w_{mn}\frac{1} {\sqrt{|A_i|}}\\ & = \frac{1}{2} (cut(A_i, \overline{A_i})\frac{1} {|A_i|} +cut(\overline{A_i},A_i)\frac{1} {|A_i|}) \\ & = \frac{cut(A_i, \overline{A_i})}{|A_i|} \end{aligned}$

1−0)2+m∈Ai,n∈Ai∑wmn(0−∣Ai∣

1)2=21m∈Ai,n∈Ai∑wmn∣Ai∣1+m∈Ai,n∈Ai∑wmn∣Ai∣

1=21(cut(Ai,Ai)∣Ai∣1+cut(Ai,Ai)∣Ai∣1)=∣Ai∣cut(Ai,Ai)
上述第（1）式用了上面第四节的拉普拉斯矩阵的性质3. 第二式用到了指示向量的定义。可以看出，对于某一个子图i，它的

R a t i o C u t

对应于

h_{i}^{T}Lh_i

,那么我们的k个子图呢？对应的

R a t i o C u t

函数表达式为：

RatioCut(A1,A2,⋯,AK)=∑i=1kW(Ai,Ai‾)∣Ai∣=∑i=1khiTLhi=∑i=1k(HTLH)ii=tr(HTLH)RatioCut(A_1,A_2, \cdots, A_K) = \sum_{i=1}^k \frac{W(A_i, \overline{A_i})}{|A_i|} \\ = \sum_{i=1}^{k}h_{i}^{T} Lh_i = \sum_{i=1}^{k}(H^{T}LH)_{ii} = tr(H^{T}LH)

其中

tr(H^{T}LH)

为矩阵的迹。也就是说，我们的 $R a t i o C u t$ 切图，实际上就是最小化我们的 $tr(H^{T}LH)$ 。注意到

H^TH=I

,则我们的切图优化目标为:

argmintr(HTLH)s.t.HTH=I⏟H\underbrace{arg \quad min \quad tr(H^{T}LH) \quad \quad s.t.H^TH = I }\\ H

注意到我们H矩阵里面的每一个指示向量都是n维的，向量中每个变量的取值为0或者

1∣Aj∣\frac{1} {\sqrt{|A_j|}}

，就有

2^n

种取值，有k个子图的话就有k个指示向量，共有

k2^n

种

H

，因此找到满足上面优化目标的H是一个NP难的问题。那么是不是就没有办法了呢？

注意观察 $tr(H^TLH)$ 中每一个优化子目标 $hiTLhih^T_iLh_i$ ,其中 $h$ 是单位正交基， $L$ 为对称矩阵，此时 $hiTLhih^T_iLh_i$ 的最大值为L的最大特征值，最小值是 $L$ 的最小特征值。如果你对主成分分析PCA很熟悉的话，这里很好理解。在 $P C A$ 中，我们的目标是找到协方差矩阵（对应此处的拉普拉斯矩阵L）的最大的特征值，而在我们的谱聚类中，我们的目标是找到目标的最小的特征值，得到对应的特征向量，此时对应二分切图效果最佳。也就是说，我们这里要用到维度规约的思想来近似去解决这个NP难的问题。

对于 $hiTLhih^T_iLh_i$ ，我们的目标是找到最小的L的特征值，而对于 $tr(HTLH)=∑i=1khiTLhitr(H^TLH)=\sum_{i=1}^k h^T_iLh_i$ ，则我们的目标就是找到k个最小的特征值，一般来说，k远远小于n，也就是说，此时我们进行了维度规约，将维度从n降到了k，从而近似可以解决这个NP难的问题。

通过找到L的最小的k个特征值，可以得到对应的k个特征向量，这k个特征向量组成一个 $\times k$ 维度的矩阵，即为我们的 $H$ 。一般需要对 $H$ 矩阵按行做标准化，即:
　　　　 $Hij∗=hij(∑t=1khit2)12H_{ij}^{*} = \frac{h_{ij}}{ (\sum_{t=1}^ {k} h_{it} ^2)^{ \frac{1}{2}} }$
　　　　
由于我们在使用维度规约的时候损失了少量信息，导致得到的优化后的指示向量h对应的H现在不能完全指示各样本的归属，因此一般在得到 $\times k$ 维度的矩阵 $H$ 后还需要对每一行进行一次传统的聚类，比如使用 $K - M e a n s$ 聚类.

6.2 PCA（主成分析）以及KPCA（核主成分分析）数学原理
PCA（Principal Component Analysis）是一种常用的数据分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据的降维。网上关于PCA的文章有很多，但是大多数只描述了PCA的分析过程，而没有讲述其中的原理。这篇文章的目的是介绍PCA的基本数学原理，帮助读者了解PCA的工作机制是什么。

当然我并不打算把文章写成纯数学文章，而是希望用直观和易懂的方式叙述PCA的数学原理，所以整个文章不会引入严格的数学推导。希望读者在看完这篇文章后能更好的明白PCA的工作原理。

数据的向量表示及降维问题
一般情况下，在数据挖掘和机器学习中，数据被表示为向量。例如某个淘宝店2012年全年的流量及交易情况可以看成一组记录的集合，其中每一天的数据是一条记录，格式如下：

(日期, 浏览量, 访客数, 下单数, 成交数, 成交金额)

其中“日期”是一个记录标志而非度量值，而数据挖掘关心的大多是度量值，因此如果我们忽略日期这个字段后，我们得到一组记录，每条记录可以被表示为一个五维向量，其中一条看起来大约是这个样子：
$500,240,25,13,2312.15)^T$
注意这里我用了转置，因为习惯上使用列向量表示一条记录（后面会看到原因），本文后面也会遵循这个准则。不过为了方便有时我会省略转置符号，但我们说到向量默认都是指列向量。

我们当然可以对这一组五维向量进行分析和挖掘，不过我们知道，很多机器学习算法的复杂度和数据的维数有着密切关系，甚至与维数呈指数级关联。当然，这里区区五维的数据，也许还无所谓，但是实际机器学习中处理成千上万甚至几十万维的情况也并不罕见，在这种情况下，机器学习的资源消耗是不可接受的，因此我们必须对数据进行降维。

降维当然意味着信息的丢失，不过鉴于实际数据本身常常存在的相关性，我们可以想办法在降维的同时将信息的损失尽量降低。

举个例子，假如某学籍数据有两列M和F，其中M列的取值是如何此学生为男性取值1，为女性取值0；而F列是学生为女性取值1，男性取值0。此时如果我们统计全部学籍数据，会发现对于任何一条记录来说，当M为1时F必定为0，反之当M为0时F必定为1。在这种情况下，我们将M或F去掉实际上没有任何信息的损失，因为只要保留一列就可以完全还原另一列。

当然上面是一个极端的情况，在现实中也许不会出现，不过类似的情况还是很常见的。例如上面淘宝店铺的数据，从经验我们可以知道，“浏览量”和“访客数”往往具有较强的相关关系，而“下单数”和“成交数”也具有较强的相关关系。这里我们非正式的使用“相关关系”这个词，可以直观理解为“当某一天这个店铺的浏览量较高（或较低）时，我们应该很大程度上认为这天的访客数也较高（或较低）”。后面的章节中我们会给出相关性的严格数学定义。

这种情况表明，如果我们删除浏览量或访客数其中一个指标，我们应该期待并不会丢失太多信息。因此我们可以删除一个，以降低机器学习算法的复杂度。

上面给出的是降维的朴素思想描述，可以有助于直观理解降维的动机和可行性，但并不具有操作指导意义。例如，我们到底删除哪一列损失的信息才最小？亦或根本不是单纯删除几列，而是通过某些变换将原始数据变为更少的列但又使得丢失的信息最小？到底如何度量丢失信息的多少？如何根据原始数据决定具体的降维操作步骤？

要回答上面的问题，就要对降维问题进行数学化和形式化的讨论。而PCA是一种具有严格数学基础并且已被广泛采用的降维方法。下面我不会直接描述PCA，而是通过逐步分析问题，让我们一起重新“发明”一遍PCA。

向量的表示及基变换

既然我们面对的数据被抽象为一组向量，那么下面有必要研究一些向量的数学性质。而这些数学性质将成为后续导出PCA的理论基础。

内积与投影
下面先来看一个高中就学过的向量运算：内积。两个维数相同的向量的内积被定义为：
$a1,a2,⋯,an)^T⋅(b1,b2,⋯,bn)^T=a_1b_1+a_2b_2+⋯+a_nb_n$
内积运算将两个向量映射为一个实数。其计算方式非常容易理解，但是其意义并不明显。下面我们分析内积的几何意义。假设A和B是两个n维向量，我们知道n维向量可以等价表示为n维空间中的一条从原点发射的有向线段，为了简单起见我们假设A和B均为二维向量，则 $A=(x_1,y_1)$ ， $B=(x_2,y_2)$ 。则在二维平面上A和B可以用两条发自原点的有向线段表示，见下图：

好，现在我们从A点向B所在直线引一条垂线。我们知道垂线与B的交点叫做A在B上的投影，再设A与B的夹角是a，则投影的矢量长度为 $∣ A ∣ c o s (a)$ ，其中 $\sqrt{x^{2}_{1}+y^2_1}$

是向量A的模，也就是A线段的标量长度。

注意这里我们专门区分了矢量长度和标量长度，标量长度总是大于等于0，值就是线段的长度；而矢量长度可能为负，其绝对值是线段长度，而符号取决于其方向与标准方向相同或相反。

到这里还是看不出内积和这东西有什么关系，不过如果我们将内积表示为另一种我们熟悉的形式:
$\cdot B = |A||B|cos(a)$

现在事情似乎是有点眉目了：A与B的内积等于A到B的投影长度乘以B的模。再进一步，如果我们假设B的模为1，即让 $∣ B ∣ = 1$ ，那么就变成了:
$\cdot B = |A|cos(a)$

基
下面我们继续在二维空间内讨论向量。上文说过，一个二维向量可以对应二维笛卡尔直角坐标系中从原点出发的一个有向线段。例如下面这个向量：

在代数表示方面，我们经常用线段终点的点坐标表示向量，例如上面的向量可以表示为(3,2)，这是我们再熟悉不过的向量表示。

不过我们常常忽略，只有一个(3,2)本身是不能够精确表示一个向量的。我们仔细看一下，这里的3实际表示的是向量在x轴上的投影值是3，在y轴上的投影值是2。也就是说我们其实隐式引入了一个定义：以x轴和y轴上正方向长度为1的向量为标准。那么一个向量(3,2)实际是说在x轴投影为3而y轴的投影为2。注意投影是一个矢量，所以可以为负。

更正式的说，向量(x,y)实际上表示线性组合：
$x(1,0)^T+y(0,1)^T$

不难证明所有二维向量都可以表示为这样的线性组合。此处 $(1, 0)$ 和$(0,1)4叫做二维空间中的一组基。

所以，要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，就可以了。只不过我们经常省略第一步，而默认以 $(1, 0)$ 和 $(0, 1)$ 为基。

我们之所以默认选择 $(1, 0)$ 和 $(0, 1)$ 为基，当然是比较方便，因为它们分别是x和y轴正方向上的单位向量，因此就使得二维平面上点坐标和向量一一对应，非常方便。但实际上任何两个线性无关的二维向量都可以成为一组基，所谓线性无关在二维平面内可以直观认为是两个不在一条直线上的向量。

例如， $(1, 1)$ 和 $(- 1, 1)$ 也可以成为一组基。一般来说，我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了！实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。例如，上面的基可以变为 $(12,−12)和(−12,12)(\frac{1}{ \sqrt{2}},-\frac{1}{ \sqrt{2}}) 和(-\frac{1}{ \sqrt{2}},\frac{1}{ \sqrt{2}})$

1,−2

1)和(−2

1,2

1)。

现在，我们想获得 $(3, 2)$ 在新基上的坐标，即在两个方向上的投影矢量值，那么根据内积的几何意义，我们只要分别计算 $(3, 2)$ 和两个基的内积，不难得到新的坐标为 $(52,−12)(\frac{5}{ \sqrt{2}},-\frac{1}{ \sqrt{2}})$

5,−2

1)。下图给出了新的基以及

(3, 2)

在新基上坐标值的示意图：

另外这里要注意的是，我们列举的例子中基是正交的（即内积为0，或直观说相互垂直），但可以成为一组基的唯一要求就是线性无关，非正交的基也是可以的。不过因为正交基有较好的性质，所以一般使用的基都是正交的。

基变换的矩阵表示

下面我们找一种简便的方式来表示基变换。还是拿上面的例子，想一下，将 $(3, 2)$ 变换为新基上的坐标，就是用 $(3, 2)$ 与第一个基做内积运算，作为第一个新的坐标分量，然后用 $(3, 2)$ 与第二个基做内积运算，作为第二个新坐标的分量。实际上，我们可以用矩阵相乘的形式简洁的表示这个变换：
$(1212−1212)(32)=(52−12)\begin{pmatrix} \frac{1}{ \sqrt{2}} & \frac{1}{ \sqrt{2}} \\ -\frac{1}{ \sqrt{2}} & \frac{1}{ \sqrt{2}} \end{pmatrix} \begin{pmatrix} 3\\ 2 \end{pmatrix} = \begin{pmatrix} \frac{5}{ \sqrt{2}} \\ -\frac{1}{ \sqrt{2}} \end{pmatrix}$

1−2

1)(32)=(2

5−2

太漂亮了！其中矩阵的两行分别为两个基，乘以原向量，其结果刚好为新基的坐标。可以稍微推广一下，如果我们有m个二维向量，只要将二维向量按列排成一个两行m列矩阵，然后用“基矩阵”乘以这个矩阵，就得到了所有这些向量在新基下的值。例如 $(1, 1)$ ， $(2, 2)$ ， $(3, 3)$ ，想变换到刚才那组基上，则可以这样表示：
$(1212−1212)(123123)=(224262000)\begin{pmatrix} \frac{1}{ \sqrt{2}} & \frac{1}{ \sqrt{2}} \\ -\frac{1}{ \sqrt{2}} & \frac{1}{ \sqrt{2}} \end{pmatrix} \begin{pmatrix} 1 & 2 & 3\\ 1 & 2 & 3 \end{pmatrix} = \begin{pmatrix} \frac{2}{ \sqrt{2}} & \frac{4}{ \sqrt{2}} & \frac{6}{ \sqrt{2}}\\ 0 & 0 & 0 \end{pmatrix}$

1−2

1)(112233)=(2

202

402

60)

于是一组向量的基变换被干净的表示为矩阵的相乘。

一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果

数学表示为：

$(p1p2⋮pR)(a1a2⋯aM)=(p1a1⋯p1aM⋮⋱⋮pRa1⋯pRaM)\begin{pmatrix} p_1 \\ p_2\\ \vdots \\ p_R \end{pmatrix} \begin{pmatrix} a_1 & a_2 & \cdots &a_M\\ \end{pmatrix} = \begin{pmatrix} p_1a_1 & \cdots & p_1a_M\\ \vdots & \ddots & \vdots \\ p_Ra_1 & \cdots & p_Ra_M \end{pmatrix}$
其中 $p_i$ 是一个行向量，表示第 $i$ 个基， $a_j$ 是一个列向量，表示第 $j$ 个原始数据记录。

特别要注意的是，这里R可以小于N，而R决定了变换后数据的维数。也就是说，我们可以将一N维数据变换到更低维度的空间中去，变换后的维度取决于基的数量。因此这种矩阵相乘的表示也可以表示降维变换。

最后，上述分析同时给矩阵相乘找到了一种物理解释：两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。更抽象的说，一个矩阵可以表示一种线性变换。很多同学在学线性代数时对矩阵相乘的方法感到奇怪，但是如果明白了矩阵相乘的物理意义，其合理性就一目了然了。

协方差矩阵及优化目标

上面我们讨论了选择不同的基可以对同样一组数据给出不同的表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。但是我们还没有回答一个最最关键的问题：如何选择基才是最优的。或者说，如果我们有一组N维向量，现在要将其降到K维（K小于N），那么我们应该如何选择K个基才能最大程度保留原有的信息？

要完全数学化这个问题非常繁杂，这里我们用一种非形式化的直观方法来看这个问题。

为了避免过于抽象的讨论，我们仍以一个具体的例子展开。假设我们的数据由五条记录组成，将它们表示成矩阵形式：
$(1124213344)\begin{pmatrix} 1 & 1 & 2 & 4 & 2 \\ 1 & 3 & 3 & 4 & 4 \end{pmatrix}$
其中每一列为一条数据记录，而一行为一个字段。为了后续处理方便，我们首先将每个字段内所有值都减去字段均值，其结果是将每个字段都变为均值为0（这样做的道理和好处后面会看到）。

我们看上面的数据，第一个字段均值为2，第二个字段均值为3，所以变换后：

$(−1−1020−20011)\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}$

我们可以看下五条数据在平面直角坐标系内的样子：

现在问题来了：如果我们必须使用一维来表示这些数据，又希望尽量保留原始的信息，你要如何选择？

通过上一节对基变换的讨论我们知道，这个问题实际上是要在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。这是一个实际的二维降到一维的问题。

那么如何选择这个方向（或者说基）才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。

以上图为例，可以看出如果向x轴投影，那么最左边的两个点会重叠在一起，中间的两个点也会重叠在一起，于是本身四个各不相同的二维点投影后只剩下两个不同的值了，这是一种严重的信息丢失，同理，如果向y轴投影最上面的两个点和分布在x轴上的两个点也会重叠。所以看来x和y轴都不是最好的投影选择。我们直观目测，如果向通过第一象限和第三象限的斜线投影，则五个点在投影后还是可以区分的。

下面，我们用数学方法表述这个问题。

方差
上文说到，我们希望投影后投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。此处，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，即：
$\frac{1}{m} \sum_{i=1}^{m} (a_i - \mu)^2$

由于上面我们已经将每个字段的均值都化为0了，因此方差可以直接用每个元素的平方和除以元素个数表示：

$\frac{1}{m} \sum_{i=1}^{m} a_i ^2$

于是上面的问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

协方差

对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。

如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。

数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：
$\frac{1}{m} \sum_{i=1}^{m} a_i b_i$

可以看到，在字段均值为0的情况下，两个字段的协方差简洁的表示为其内积除以元素数m。

当协方差为0时，表示两个字段完全独立。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。

至此，我们得到了降维问题的优化目标：将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差）。

协方差矩阵

上面我们导出了优化目标，但是这个目标似乎不能直接作为操作指南（或者说算法），因为它只说要什么，但根本没有说怎么做。所以我们要继续在数学上研究计算方案。

我们看到，最终要达到的目的与字段内方差及字段间协方差有密切关系。因此我们希望能将两者统一表示，仔细观察发现，两者均可以表示为内积的形式，而内积又与矩阵相乘密切相关。于是我们来了灵感：

假设我们只有a和b两个字段，那么我们将它们按行组成矩阵X：
$\begin{pmatrix} a_1 & a_2 & \cdots & a_m \\ b_1 & b_2 & \cdots & b_m \end{pmatrix}$

然后我们用X乘以X的转置，并乘上系数1/m：

$1mXXT=(1m∑i=1mai21m∑i=1maibi1m∑i=1maibi1m∑i=1mbi2)\frac{1}{m}XX^T= \begin{pmatrix} \frac{1}{m} \sum_{i=1}^m a_i^2 & \frac{1}{m} \sum_{i=1}^m a_ib_i\\ \frac{1}{m} \sum_{i=1}^m a_ib_i & \frac{1}{m} \sum_{i=1}^m b_i^2 \end{pmatrix}$

奇迹出现了！这个矩阵对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。两者被统一到了一个矩阵的。

根据矩阵相乘的运算法则，这个结论很容易被推广到一般情况：
设我们有 $m$ 个 $n$ 维数据记录，将其按列排成n乘m的矩阵X，设 $C=1mXXTC=\frac{1}{m}XX^T$ ，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第 $i$ 行 $j$ 列和 $j$ 行i列元素相同，表示 $i$ 和 $j$ 两个字段的协方差。

方法一
设原始数据矩阵 $X$ 对应的协方差矩阵为 $C$ ，而 $P$ 是一组基按行组成的矩阵，设 $Y = P X$ ，则 $Y$ 为 $X$ 对 $P$ 做基变换后的数据。设 $Y$ 的协方差矩阵为 $D$ ，我们推导一下 $D$ 与 $C$ 的关系：

$D=1mYYT=1m(PX)(PX)T=1mPXXTPT=P(1mXXT)PT=PCPT\begin{aligned} D &= \frac{1}{m} YY^T \\ &= \frac{1}{m}(PX)(PX)^T \\ &= \frac{1}{m}PXX^TP^T \\ &= P( \frac{1}{m} XX^T)P^T \\ &= PCP^T \end{aligned}$

现在事情很明白了！我们要找的 $P$ 不是别的，而是能让原始协方差矩阵对角化的 $P$ 。换句话说，优化目标变成了寻找一个矩阵 $P$ ，满足 $PCP^T$ 是一个对角矩阵，并且对角元素按从大到小依次排列，那么 $P$ 的前 $K$ 行就是要寻找的基，用 $P$ 的前 $K$ 行组成的矩阵乘以 $X$ 就使得 $X$ 从 $N$ 维降到了 $K$ 维并满足上述优化条件。

至此，我们离“发明”PCA还有仅一步之遥！

现在所有焦点都聚焦在了协方差矩阵对角化问题上，有时，我们真应该感谢数学家的先行，因为矩阵对角化在线性代数领域已经属于被玩烂了的东西，所以这在数学上根本不是问题。

由上文知道，协方差矩阵C是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：

1）实对称矩阵不同特征值对应的特征向量必然正交。

2）设特征向量λ重数为r，则必然存在r个线性无关的特征向量对应于λ，因此可以将这r个特征向量单位正交化。

由上面两条可知，一个 $n$ 行 $n$ 列的实对称矩阵一定可以找到n个单位正交特征向量，设这 $n$ 个特征向量为 $e_1,e_2,⋯,e_n$ ，我们将其按列组成矩阵：
$E=(e_1,e_2,⋯,e_n)$
则对协方差矩阵C有如下结论：
$ETCE=Λ=(λ1λ2⋱λn)E^TCE= \Lambda = \begin{pmatrix} \lambda_1& & &\\ & \lambda_2 & &\\ & &\ddots &\\ & & & \lambda_n \end{pmatrix}$
其中 $Λ$ 为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。

以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。

到这里，我们发现我们已经找到了需要的矩阵 $P$ ：
$P=E^T$

$P$ 是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是 $C$ 的一个特征向量。如果设 $P$ 按照 $Λ$ 中特征值的从大到小，将特征向量从上到下排列，则用P的前K行组成的矩阵乘以原始数据矩阵 $X$ ，就得到了我们需要的降维后的数据矩阵 $Y$ 。

除了矩阵推导之外还有利用拉格朗日乘数法——方法二

至此我们完成了整个 $P C A$ 的数学原理讨论。在下面的一节，我们将给出 $P C A$ 的一个实例。

想法分析：一个数学模型，求取带有限条件函数（矩阵）的极值。自然想到是拉格朗日乘数法。故我们把它转换成为最优化问题利用拉格朗日乘数法来推导：

样本点 $x_i$ 在基 $W$ 下的坐标为： $x_i, w) = x_i^T w$ , 于是我们有方差：
$D(x)=1m∑i=1m(xiTw)2=1m∑i=1m(xiTw)T(xiTw)=1m∑i=1mwxixiTw=wT(1m∑i=1mxixiT)w\begin{aligned} D(x) &= \frac{1}{m} \sum_{i=1}^{m}(x_i^Tw)^2 \\ &=\frac{1}{m} \sum_{i=1}^{m}(x_i^Tw)^T(x_iTw) \\ &=\frac{1}{m} \sum_{i=1}^{m} wx_ix_i^Tw \\ &=w^T(\frac{1}{m}\sum_{i=1}^{m} x_ix_i^T)w \end{aligned}$

看到 $1m∑i=1mxixiT\frac{1}{m}\sum_{i=1}^{m} x_ix_i^T$ 就是原样本的协方差，我们另一个矩阵为 $Λ\Lambda$ , 于是我们就有了：
${maxwTΛws.t.wTw=1\begin{cases} max{w^T \Lambda w} \\ s.t.w^Tw = 1 \end{cases}$
此时就可以用拉格朗日乘数法——构造拉格朗日函数：
$w^T \Lambda w= \lambda(1-w^Tw)$

对 $w$ 求导得：
$\Lambda w - \lambda w=0 \\ \Lambda w= \lambda w$
此时方差为：
$w^T \Lambda w = \lambda w^Tw =\lambda$
于是我们发现， $x$ 投影后的方差就是协方差矩阵的特征值。我们要找到最大方差也就是协方差矩阵最大的特征值，最佳投影方向就是最大特征值所对应的特征向量，次佳就是第二大特征值对应的特征向量，以此类推。

至此我们完成了基于最大可分性的 $P C A$ 数学证明。

PCA算法
总结一下PCA的算法步骤：

设有 $m$ 条 $n$ 维数据。

1）将原始数据按列组成 $n$ 行 $m$ 列矩阵 $X$

2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值

3）求出协方差矩阵 $C=1mXXTC=\frac{1}{m} XX^T$
4）求出协方差矩阵的特征值及对应的特征向量

5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P

6） $Y = P X$ 即为降维到 $k$ 维后的数据

KPCA（核主成分分析）

KPCA，中文名称”核主成分分析“，是对PCA算法的非线性扩展，言外之意，PCA是线性的，其对于非线性数据往往显得无能为力，例如，不同人之间的人脸图像，肯定存在非线性关系，自己做的基于ORL数据集的实验，PCA能够达到的识别率只有88%，而同样是无监督学习的KPCA算法，能够轻松的达到93%左右的识别率（虽然这二者的主要目的是降维，而不是分类，但也可以用于分类），这其中很大一部分原因是，KPCA能够挖掘到数据集中蕴含的非线性信息。

PCA使用的方式：线性可分得数据可以直接使用PCA，通过降维来进行数学向量表示。
KPCA使用的方式：线性不可分的数据再使用PCA就得不到很好的效果，甚至会产生负面效果。故把此类数据先放到高维度下，使其变得线性可分，再使用PCA。

例如：

如左图，红色部分的点为一类数据，黑色部分的点为另一类，在一维空间中，你不可能通过一刀切将两类数据分开，至少需要两刀。OK，这就说明数据分布是非线性的，我们采用高维映射，当然了，例子中只是映射到了二维空间，但已经足够说明问题了，在右图中，完全可以通过沿着X轴方向的一刀切将两类数据分开，说明在二维空间中，数据已经变成线性可分的了。

例如：

在原始数据中，明显没法使用线性的切割方式来把2类点分开，强行使用PCA会适得其反，而采用KPCA会起到非常好的效果。

代码如下：

from sklearn.datasets import make_moons
from sklearn.decomposition import PCA
from sklearn.decomposition import KernelPCA
import matplotlib.pyplot as plt
import numpy as npx2, y2 = make_moons(n_samples=100, random_state=123)if __name__ == "__main__":pca = PCA(n_components=2)x_spca = pca.fit_transform(x2)fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(14, 6))ax[0].scatter(x_spca[y2 == 0, 0], x_spca[y2 == 0, 1], color='red', marker='^', alpha=0.5)ax[0].scatter(x_spca[y2 == 1, 0], x_spca[y2 == 1, 1], color='blue', marker='o', alpha=0.5)ax[1].scatter(x_spca[y2 == 0, 0], np.zeros((50, 1)) + 0.02, color='red', marker='^', alpha=0.5)ax[1].scatter(x_spca[y2 == 1, 0], np.zeros((50, 1)) + 0.02, color='blue', marker='o', alpha=0.5)ax[0].set_xlabel('PC1')ax[0].set_ylabel('PC2')ax[1].set_ylim([-1, 1])ax[1].set_yticks([])ax[1].set_xlabel('PCA')kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)x_kpca = kpca.fit_transform(x2)ax[2].scatter(x_kpca[y2 == 0, 0], x_kpca[y2 == 0, 1], color='red', marker='^', alpha=0.5)ax[2].scatter(x_kpca[y2 == 1, 0], x_kpca[y2 == 1, 1], color='blue', marker='o', alpha=0.5)ax[2].set_xlabel('kPCA')plt.show()

由以上例子看出KPCA的优势。接下来正式介绍KPCA。

理论

KPCA的公式推导和PCA十分相似，只是存在两点创新：

为了更好地处理非线性数据，引入非线性映射函数 $Φ\Phi$ ，将原空间中的数据映射到高维空间，注意，这个 $Φ\Phi$ 是隐性的，我们不知道，也不需要知道它的具体形式是啥。
引入了一个定理：空间中的任一向量（哪怕是基向量），都可以由该空间中的所有样本线性表示，这点对KPCA很重要，我想大概当时那个大牛想出KPCA的时候，这点就是它最大的灵感吧。话说这和”稀疏“的思想比较像。

假设中心化后的样本集合 $X$ （ $\times N$ ，N个样本，维数d维，样本”按列排列“），现将 $X$ 映射到高维空间，得到，假设在这个高维空间中，本来在原空间中线性不可分的样本现在线性可分了，然后呢？想啥呢！果断上PCA啊！

于是乎！假设 $D （ D > > d ）$ 维向量 $wi(i=1,⋯,d)w_i(i=1, \cdots, d)$ 为高维空间中的特征向量， $λi(i=1,⋯,d)\lambda _{i}(i=1, \cdots, d)$ 为对应的特征值，根据PCA的数学模型 $XXTwi=λiwiXX^Tw_i = \lambda _i w_i$ ，高维空间中的PCA如下：

$Φ(X)Φ(X)Twi=λiwi\Phi (X) \Phi (X)^T w_i = \lambda_iw_i$
我们根据基的概念得到：一个空间中的基向量完全可以用所有的训练样本线性表示出来，即：
$wiΦ=∑k=1Nαkϕ(xk)=Φ(X)αiw_i^{\Phi} = \sum_{k=1}^{N} \alpha_k \phi(x_k)=\Phi(X) \alpha_i$

带入到高维PCA的表示式中得：
$Φ(X)Φ(X)TΦ(X)αi=λiΦ(X)αi\Phi (X) \Phi (X)^T \Phi(X) \alpha_i = \lambda_i\Phi(X) \alpha_i \\$
变换,等号两侧同时左乘一个东西
$Φ(X)TΦ(X)Φ(X)TΦ(X)αi=λiΦ(X)TΦ(X)αi\Phi(X)^T\Phi (X) \Phi (X)^T \Phi(X) \alpha_i = \lambda_i\Phi(X)^T\Phi(X) \alpha_i$

此时定义核函数使得下式成立;
$k(xi，xj)=Φ(xi)Φ(xj)k(x_i，x_j)= \Phi(x_i) \Phi(x_j)$

请注意，这个公式表达的意思只是高维空间中向量的内积等于原空间中向量之间的核函数值，根核函数的具体形式没有一毛钱的关系。我们继续，一个样本是这样的，那么矩阵化的形式如下：
$\Phi(X)^T\Phi(X)$

那么高维PCA就变成了：
$K2αi=λiKαiKαi=λiαiK^2 \alpha_i = \lambda_i K \alpha_i \\ K \alpha_i = \lambda_i \alpha_i$

求解公式的含义就是求K最大的几个特征值所对应的特征向量，由于K为对称矩阵，所得的解向量彼此之间肯定是正交的。

但是，请注意，这里的 $α\alpha$ 只是 $K$ 的特征向量，但是其不是高维空间中的特征向量，高维空间中的特征向量 $w$ 应该是由 $α\alpha$ 进一步求出。

其中它借用核函数的性质，将高维映射 $w$ 成功的抹掉了。所以我们根本就不需要关心高维映射的基向量。

解决另一个问题又带来另一个问题，那么什么样的函数可以被称作和函数呢？这里大牛们给出了一个著名的定理，称作mercer定理。

Mercer 定理：任何半正定的函数都可以作为核函数。所谓半正定的函数 $f(x_i,x_j)$ ，是指拥有训练数据集合 $x_1,x_2,...x_n)$ ，我们定义一个矩阵的元素 $a_{ij} = f(x_i,x_j)$ ，这个矩阵式 $n×nn\times n$ 的，如果这个矩阵是半正定的，那么 $f(x_i,x_j)$ 就称为半正定的函数。

请注意，这个mercer定理不是核函数必要条件，只是一个充分条件，即还有不满足mercer定理的函数也可以是核函数。所谓半正定指的就是核矩阵K的特征值均为非负。

证明过程：https://www.cnblogs.com/jerrylead/archive/2011/03/18/1988406.html

其中一些常见的核函数在：https://blog.csdn.net/wsj998689aa/article/details/47027365

至此PCA以及KPCA介绍完成。

6.3 Ncut切图

$N c u t$ 切图和 $R a t i o C u t$ 切图很类似，但是把 $R a t i o c u$ t的分母 $A_i|$ 换成 $vol(A_i)$ . 由于子图样本的个数多并不一定权重就大，我们切图时基于权重也更合我们的目标，因此一般来说 $N c u t$ 切图优于 $R a t i o C u t$ 切图。
$NCut(A1,A2,⋯,Ak)=12∑i=1kW(Ai,Ai‾)vol(Ai)NCut(A_1, A_2, \cdots, A_k) = \frac{1}{2} \sum_{i=1}^{k} \frac{W(A_i, \overline{A_i})}{vol(A_i)}$

对应的， $N c u t$ 切图对指示向量h做了改进。注意到 $R a t i o C u t$ 切图的指示向量使用的是 $1vol(Ai)\frac{1}{\sqrt{vol(A_i)}}$

1标示样本归属，而Ncut切图使用了子图权重

1vol(Ai)\frac{1}{\sqrt{vol(A_i)}}

来标示指示向量

h

，定义如下:

hij={0vi∉Ai1vol(Ai)vi∈Aih_{ij} = \begin{cases} 0 & v_i \not\in A_i \\ \frac{1}{\sqrt{vol(A_i)}} & v_i \in A_i \end{cases}

那么我们对于

h_i^T L h_i

就有 (根据前面的推导）:

hiTLhi=cut(Ai,Ai‾)vol(Ai)h_i^T L h_i = \frac{cut(A_i, \overline{A_i})}{vol(A_i)}

同样的我们的优化目标推到方式和 $R a t i o C u t$ 切图一样:
$NCut(A1,A2,⋯,Ak)=∑i=1kcut(Ai,Ai‾)vol(Ai)=tr(HTLH)NCut(A_1, A_2, \cdots, A_k) = \sum_{i=1}^{k} \frac{cut(A_i, \overline{A_i})}{vol(A_i)} = tr(H^TLH)$

有一点不同的是此时， $HTH≠IH^TH \ne I$ 而是 $H^TDH = I$ .，推导如下：
$hiTDhi=∑j=1khij2dj=1vol(Ai)∑j=1kdjh_i^T D h_i = \sum_{j=1}^{k} h_{ij}^2 d_j = \frac{1}{vol(A_i)} \sum_{j=1}^{k} d_j$

根据 $vol(Ai)=∑j=1kdjvol(A_i) = \sum_{j=1}^k d_j$ :

则：
$\frac{1}{vol(A_i)} vol(A_i)=1$

也就是说，此时我们的优化目标最终为：

$argmintr(HTLH)s.t.HTDH=I⏟H\underbrace{arg \quad min \quad tr(H^{T}LH) \quad \quad s.t.H^TDH = I }\\ H$

argmintr(HTLH)s.t.HTDH=IH

此时我们的 $H$ 中的指示向量 $h$ 并不是标准正交基，所以在 $R a t i o C u t$ 里面的降维思想不能直接用。怎么办呢？其实只需要将指示向量矩阵 $H$ 做一个小小的转化即可。

令 $D^{-\frac{1}{2}} F$ , 则： $HTLH=FTD−12LD−12F,HTDH=FTF=IH^TLH = F^TD^{- \frac{1}{2}} L D^{-\frac{1}{2}}F, H^TDH = F^TF=I$ 也就是说我们的优化目标：
$argmintr(HTLH)s.tFTF=I⏟H\underbrace{arg \quad min \quad tr(H^{T}LH) \quad \quad s.t F^TF = I }\\ H$

argmintr(HTLH)s.tFTF=IH
可以发现这个式子和

R a t i o C u t

基本一致，只是中间的L变成了

D−12LD−12D^{- \frac{1}{2}} L D^{-\frac{1}{2}}

。这样我们就可以继续按照

R a t i o C u t

的思想，求出

D−12LD−12D^{- \frac{1}{2}} L D^{-\frac{1}{2}}

的最小的前k个特征值，然后求出对应的特征向量，并标准化，得到最后的特征矩阵

F

,最后对F进行一次传统的聚类（比如K-Means）即可。

一般来说， $D−12LD−12D^{- \frac{1}{2}} L D^{-\frac{1}{2}}$ 相当于对拉普拉斯矩阵L做了一次标准化即: $Lijdi∗dj\frac{L_{ij}}{\sqrt{d_i * d_j}}$

Lij

7. 谱聚类算法流程

铺垫了这么久，终于可以总结下谱聚类的基本流程了。一般来说，谱聚类主要的注意点为相似矩阵的生成方式（参见第二节），切图的方式（参见第六节）以及最后的聚类方法（参见第六节）。

最常用的相似矩阵的生成方式是基于高斯核距离的全连接方式，最常用的切图方式是Ncut。而到最后常用的聚类方法为K-Means。下面以Ncut总结谱聚类算法流程。

输入：样本集 $D=(x_1,x_2,...,x_n)$ ，相似矩阵的生成方式, 降维后的维度 $k_1$ , 聚类方法，聚类后的维度 $k_2$
输出：簇划分 $C(c_1,c_2,...c_{k_2})$ .

1）根据输入的相似矩阵的生成方式构建样本的相似矩阵 $S$

2）根据相似矩阵 $S$ 构建邻接矩阵 $W$ ，构建度矩阵 $D$

3）计算出拉普拉斯矩阵 $L$

4）构建标准化后的拉普拉斯矩阵 $D−12LD−12D^{- \frac{1}{2}} L D^{-\frac{1}{2}}$

5）计算 $D−12LD−12D^{- \frac{1}{2}} L D^{-\frac{1}{2}}$ 最小的 $k_1$ 个特征值所各自对应的特征向量 $f$

6）将各自对应的特征向量f组成的矩阵按行标准化，最终组成 $n×k1n\times k_1$ 维的特征矩阵 $F$

7）对 $F$ 中的每一行作为一个 $k_1$ 维的样本，共 $n$ 个样本，用输入的聚类方法进行聚类，聚类维数为 $k_2$ 。

8）得到簇划分 $C(c_1,c_2,...c_{k_2})$ .

8. 谱聚类算法总结
　　　　谱聚类算法是一个使用起来简单，但是讲清楚却不是那么容易的算法，它需要你有一定的数学基础。如果你掌握了谱聚类，相信你会对矩阵分析，图论有更深入的理解。同时对降维里的主成分分析也会加深理解。

下面总结下谱聚类算法的优缺点。

谱聚类算法的主要优点有：

1）谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效。这点传统聚类算法比如K-Means很难做到

2）由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好。

谱聚类算法的主要缺点有：

1）如果最终聚类的维度非常高，则由于降维的幅度不够，谱聚类的运行速度和最后的聚类效果均不好。

2) 聚类效果依赖于相似矩阵，不同的相似矩阵得到的最终聚类效果可能很不同。

参考：

https://www.cnblogs.com/pinard/p/6221564.html
https://blog.csdn.net/wsj998689aa/article/details/40398777
http://blog.codinglabs.org/articles/pca-tutorial.html
https://zhuanlan.zhihu.com/p/85977023
https://zhuanlan.zhihu.com/p/100586855