PCA、SVD、谱聚类

1. PCA

所谓降维，就是要把n维向量X(i)投射到k维的空间（k<n），在这个k维空间里面，样本点的投影长度尽可能大，这样就能保证这个新的空间保留了原来n维空间里面尽可能多的variance。下面是公式描述：
if x⁽ⁱ⁾ is a point(n x 1), then its projection onto a unite vector u (n x 1) is the distance：x^(i)Tu
Hence, to maximize the variance of the projections, we
would like to choose a unit-length u so as to maximize:


where X=(x⁽¹⁾ x⁽²⁾ ……x^(m)) is n x m matrix, and n is the number of features
如果用拉格朗日乘子法，对u求偏导并另其为0，我们可以发现：

最大值就是XX^T的最大特征值，该特征值对应的单位特征向就是我们要求的u。
上面的例子u是一维空间，下面扩展到k维空间U=(u₁ u₂ ………u_k)，这里每个u_j都是n x 1的单位向量，因此U为n x k矩阵。把x⁽ⁱ⁾投影到这k个向量张成的空间，得到k x 1维向量：

我们的最优化问题，就是最大化k个投影分量，也就是z⁽ⁱ⁾的模：

采用拉格朗日乘子法

根据偏导公式：

对U求偏导，得

化简后得：

也就是求XX^T的最大的k个特征根及对应的特征向量，需要补充的是，该矩阵是对称阵，因此一定有n个实数特征值。
最后如果要还原到n维向量空间

实际操作的时候，需要首先对所有的样本进行标准化。另外，我们不指定降维后的k的值，而是指定一个降维到的主成分比重阈值t。这个阈值t在（0,1]之间。假如我们的n个特征值为λ 1 ≥λ 2 ≥…≥λ n,则k可以通过下式得到:

2. SVD

SVD也是对矩阵进行分解，但是和特征分解不同，SVD并不要求要分解的矩阵为方阵。假设我们的矩阵A是一个m×n的矩阵，那么我们定义矩阵A的SVD为：

其中U是一个m×m的矩阵，Σ是一个m×n的矩阵，除了主对角线上的元素以外全为0，主对角线上的每个元素都称为奇异值，V是一个n×n的矩阵。U和V都是酉矩阵（实数域上也就是正交矩阵），即满足U^TU=I,V^TV=I。

A^T A是一个nxn的对称矩阵（事实上对于任何矩阵A,A^TA是半正定的，也就是说所有特征值大于等于0），所有特征值分别代入特征方程后，所有基础解系的总数为n，而且一定正交相似于对角矩阵，也就是说V矩阵是A^TA基础解系正交化后张成的。而对于nxn的矩阵Σ^T Σ=Σ^{每个对角线元素取平方}（如果m<n,则奇异值只有m个，维度不够的填零，即A^T A的0特征值），就是A^T A的特征值矩阵，也就是说，每个奇异值是A^TA的每个特征值开平方σ_i=sqrt(λ_i )

AA^T是一个mxm的对称矩阵（事实上对于任何矩阵A,AA^T是半正定的，也就是说所有特征值大于等于0），所有特征值分别代入特征方程后，所有基础解系的总数为m，而且一定正交相似于对角矩阵，也就是说U矩阵是AA^T基础解系正交化后张成的。而对于m x m的矩阵 ΣΣ^T =Σ^{每个对角线元素取平方}（如果m>n,则奇异值只有n个，维度不够的填0，即AA^T的0特征值），就是AA^T的特征值矩阵，也就是说，每个奇异值是AA^T的每个特征值开平方σ_i=sqrt(λ_i )。
AA^T和A^T A非零特征值相同，奇异值矩阵中，奇异值也是按照从大到小排列，而且奇异值的减少特别快，在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上的比例。也就是说，我们也可以用最大的k个的奇异值和对应的左右奇异向量来近似描述矩阵。也就是说：

k要比m和n小很多，也就是一个大的矩阵A可以用三个小的矩阵来表示。由于这个重要的性质，SVD可以用于PCA降维，来做数据压缩和去噪!
对于含有m个样本（每个样本有n个特征）的m x n样本集矩阵

可以通过SVD的右矩阵V得到X^T X (大小为n x n）最大的k个特征向量（列向量）张成的矩阵V^T_k×n（大小为n x k）用于特征降维，相当于是说我们把所有特征归纳成了k类；还可以用左矩阵U 得到XX^T（大小为m x m）的前d个特征值张成的矩阵U_m×d（大小为m x d）, 相当于把所有样本归纳成了d类

X’是一个d×n的矩阵,这个矩阵和我们原来的m×n的样本相比，行数从m减到了d，实现了样本压缩。

3. LDA

首先来看看瑞利商的定义。瑞利商是指这样的函数:

其中x为非零向量，而A为的Hermitan矩阵。所谓的Hermitan矩阵就是满足共轭转置矩阵和自己相等的矩阵， 如果我们的矩阵A是实矩阵，则对称矩阵即为Hermitan矩阵。
瑞利商有一个非常重要的性质，即它的最大值等于矩阵A最大的特征值，而最小值等于矩阵A的最小的特征值。
广义瑞利商是指这样的函数 :

其中x为非零向量，而A,B为n x n的Hermitan矩阵。B为正定矩阵。它的最大值和最小值是什么呢？其实我们只要通过将其通过标准化就可以转化为瑞利商的格式。B是对称矩阵，一定可以被一个 正交矩阵对角化，即存在正交矩阵C，使得B=CDC^T=CD^{对角线元素开根号}D^{对角线元素开根号}C^T=CD^{对角线元素开根号}C^TCD^{对角线元素开根号}C^T=(CD^{对角线元素开根号}C^T) (CD^{对角线元素开根号}C^T)。B^1/2=CD^{对角线元素开根号}C^T，对这个矩阵再求逆，B^-1/2= （CD^{对角线元素开根号}C^T）^-1=C^-TD^{对角线元素开根号再求倒数}C^-1。由于C是正交矩阵C^-1=C^T, C^-T=C，因此B^-1/2= CD^{对角线元素开根号再求倒数}C^T。B^-1/2也是对称矩阵。
B^-1/2 B B^-1/2= CD^{对角线元素开根号再求倒数}C^T CD C^T CD^{对角线元素开根号再求倒数}C^T= CD^{对角线元素开根号再求倒数}D D^{对角线元素开根号再求倒数}C^T=CC^T=E

B^-1/2 B^-1/2= CD^{对角线元素开根号再求倒数}C^T CD^{对角线元素开根号再求倒数}C^T= CD^{对角线元素开根号再求倒数}D^{对角线元素开根号再求倒数}C^T=CD^{对角线元素求倒数}C^T=B^-1
令x=B^(-1/2)x’
则分母转化为：

而分子转化为：

广义瑞利商转化成了

它的极值等于矩阵B^-1/2AB^-1/2的特征值。
对于特征值问题B^-1/2AB^-1/2 x^’=λx^’，左乘B^-1/2，得B^-1 AB^-1/2 x^’=λB^-1/2 x^’

运用x=B^(-1/2) x^’，上式得B^-1 Ax=λx
因此B^-1 A的特征值问题就是上述广义瑞利商的极值。
现在我们回到LDA的原理上，对于二类问题希望找到一个向量w,使得同一种类别数据的投影点尽可能的接近，而不同类别的数据的类别中心之间的距离尽可能的大。对任意一个样本本xⁱ,它在向量w上的投影(xⁱ)’=w^Txⁱ。它是实数。

假设我们的数据集D={(x¹,y¹),(x²,y²),…,((x^m,y^m))},其中任意样本xⁱ为n维向量，我们定义N_j(j=0,1)为第j类样本的个数，X_j(j=0,1)为第j类样本的集合。
μ_j(j=0,1)为第j类样本的均值向量(n x 1)

第j类样本的所有投影的均值（scaler）

定义Σ_j(j=0,1)为第j类样本的协方差矩阵（n x n）

第j类样本的所有投影的协方差为(scaler)

（如果A为n阶对角矩阵，由二次型知w^TAw=w₁²A₁₁+w₂ ²A₂₂+….w_n²A_nn为一个scaler）
因此我们 的最优化目标就是

我们定义类内散度S_w为(nx n)

我们再定义类间散度S_b为(n x n)

最优化目标(scaler)

最大值为矩阵S_w⁻¹S_b的最大特征值，要找的向量w为对应的特征向量。

而(μ₀-μ₁)^T w是一个标量，因此S_bw的方向恒为μ₀−μ₁。不妨令S_bw=λ(μ₀−μ₁)，将其带入：(S_w⁻¹S_b)w=λw，可以得到
w=S_w ^-1(μ₀−μ₁)
如果一共要分成N类。假设我们投影到的低维空间的维度为d，对应的基向量为(w₁,w₂,…w_d)，基向量组成的矩阵为W, 它是一个n×d的矩阵。
类内散度矩阵为(n x n)


类间散度矩阵为(n x n)

其中μ是所有示例的均值向量
下面我们来分析这个d x d的矩阵

(如果W为n阶对角矩阵，则W^TAW相当于对A的每一个元素作A_ij*W_ii*W_jj，i,j=1…n)
这个矩阵的每个对角元素，就相当于二类问题里面投影到一个向量w的情况。因此我们可以把最优化问题写为：

也就是说，转化为对d个(w_i^TS_bw_i)/(w_i^T S_w w_i )子问题最优化，每个子问题的最大值是矩阵S_w⁻¹S_b的最大特征值，因此要最大化J，就是求S_w⁻¹S_b的d个最大特征值的乘积, 此时要找的投影空间W为这d个特征值对应的特征向量张成的矩阵。
(μ_j-μ)为一个nx1的矩阵，For A ∈ R_m×n, rank(A) ≤ min(m, n)，因此对于一个非零矩阵，其秩为1。A ∈ R_m×n, B ∈ R_n×p, rank(AB) ≤ min(rank(A), rank(B))，因此(μ_j-μ)(μ_j-μ)^T的秩也为1，再由 rank(A + B) ≤ rank(A) + rank(B)得

为N项相加，因此秩最大为N,而μ_k可以由μ₁…… μ_k-1线性表出，也就是S_b可以化简为N-1项相加，因此S_b的秩最大为N-1。再由rank(AB) ≤ min(rank(A), rank(B))得S_w⁻¹S_b的秩最大为N-1。矩阵的秩就是非零特征值的数目，因此对于第i个最优子问题，如果i大于N-1，S_w⁻¹S_b的特征值为0，将其带入损失函数J，最大化目标就没有意义了。因此我们最多取N-1个特征值，也就是取N-1个特征向量，因此最多降维到N-1维空间W=(w₁,w₂,…w_N-1)。

4. 谱聚类

主要思想：把所有的数据看做空间当中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值比较低，较近的点权重值比较高，目的是：进行切图，让切图后不同的子图之间边权重和尽可能低，而子图内部的边权重和尽可能的高，从而达到聚类的目的。

4.1 无向权重图

由点集合V，边集合E来描述的。G（V，E）

定义权重wij为点vi，vj之间的权重，由于是无向图，所以wij=wji。如果是有向图，wij不等于wji（一般情况下）。对于有连接的点vi，vj，wij>0，如果没有连接，例如上图中的v2，v4，我们认为w24=w42=0。利用wij，我们可以建立图的邻接矩阵W（nn, n为样本数），第i行第j列对应wij。在W矩阵当中，我们定义对角线上的元素均为0.

定义：di为点vi和它相连的所有边的权重之和。d2=w21+w22+w23+w24+w25。这里w22=0。利用di，我们可以建立一个矩阵D（nn），一个对角阵，主对角线为dn，其他的位置均为0.

定义：对于点集V的一个子集 A属于V。|A|：是子集A当中点的个数。Vol（A）：定义了子集A内的点其所有边的权重和

4.2 相似矩阵

基本思想：距离较远的两个点之间的权重值较低，距离较近的两个点之间权重值较高。
构造相似矩阵S就是在对wij究竟如何取值，进行定义。

1. ε-邻近法：
   设置阈值ε，用欧氏距离sij来度量vi，vj之间的距离，然后sij=||vi-vj||^2,根据sij和ε的大小关系来定义W。
   如果wij={0 sij>ε
   ε，sij<=ε}
   这种定义不够精细，实际应用比较少；
2. K近邻法：利用KNN算法来遍历所有样本点，取每个样本最近的K个点来作为我们的近邻点。只有离i最近的几个点j才有wij，i的K近邻含j的时候，j的k近邻可能没有i，因此wij不等于wji。所以有三种权重定义方法：

• 只要一个点在另一个点的K近邻中，则得权重wij=wji=
  
• 两个点都必须在互为K近邻时才有权重wij=wji=
  
• 全连接法：保证了所有点的权重都大于0可以选择不同的核函数来定义边权重，常用的有多项式核函数、高斯核函数、sigmoid核函数。

4.3 拉普拉斯矩阵

L=D-W（D是对称阵，W也是对称阵）（L也是对称阵），由于L是对称阵，所以特征值都是实数。对于任意向量f有：




拉普拉斯矩阵是半正定的，且对应的n个实数特征值0=λ1<=λ2<=λn，对应的特征项向量的模长是1。

4.4 无向图切图

对于无向图G的切图，我们的目标是将图G(V,E)切成相互没有连接的k个子图，每个子图的点集为：A₁,A₂,…A_k，它们满足A_i∩A_j=∅,且A₁∪A₂∪…∪A_k=V。对于任意两个子图点的集合A,B⊂V, A∩B=∅, 我们定义A和B之间的切图权重为：

那么对于我们k个子图的点集：A₁,A₂,…A_k，我们定义切图cut为：

其中 A¯_i为Ai的补集，意为除A_i子集外其他V的子集的并集。

我们选择一个权重最小的边缘的点，比如C和H之间进行cut，这样得到的cut(A₁,A₂,…A_k)包含的w_ij项最少， 但是却不是最优的切图。我们需要对每个子图的规模做出限定，一般来说，有两种切图方式，第一种是RatioCut，第二种是Ncut。

1. ratiocut
   不光考虑最小化cut(A₁,A₂,…A_k)，它还同时考虑最大化每个子图点的个数，即：
   
   我们引入指示向量
   
   j=1,2,…k, n为样本数。对于任意一个向量h_j，它的各个元素：
   
   对于每个子图i：
   
   
   其中H=（h₁ h₂……h_k）为n x k的矩阵
   实际上就是最小化tr(H^TLH)，并且注意到H^TH=I,则我们的切图优化目标为:
   
   其中h是单位正交基， L为对称矩阵，目标tr(H^TLH)的每一个子目标h_i^TLh_i的最大值为L的最大特征值，最小值是L的最小特征值。我们的目标是找到k个最小的特征值，一般来说，k远远小于n，也就是说，此时我们进行了维度规约，将维度从n降到了k。另外，我们得到了对应的k个特征向量，这k个特征向量组成一个nxk维度的矩阵，即为我们要求的H。
   
   从上图看出，每一行对应一个样本，一个样本点只会属于一个子图，H就体现了各个样本的分类情况。如果某个样本在第j列为0，由指示函数hj的定义可知，该样本不属于第j个子图，也就是不属于第j类,j=1,2…k。由于我们只取另外一部分特征向量，导致有些样本无法确定归属（比如上图的第一行，即第一个样本）。因此一般在得到nxk维度的矩阵H后还需要对每一行进行一次传统的聚类，比如使用K-Means聚类。

2. Ncut
   
   子图样本的个数多并不一定权重就大，因此一般来说Ncut切图优于RatioCut切图。
   
   
   推导方式和RatioCut完全一致。也就是说，我们的优化目标仍然是
   
   但是此时我们的H^TH≠I,而是H^TDH=I。推导如下：
   
   此时我们的H中的各个指示向量h_j不是单位正交 的，所以在RatioCut里面的降维思想不能直接用。其实只需要将指示向量矩阵H做一个小小的转化即可。
   我们令H=D^−1/2F, 则：H^TLH=F^TD^−1/2LD^−1/2F，H^TDH=F^TF=I,也就是说优化目标变成了:
   
   这样我们就可以继续按照RatioCut的思想，求出D^1/2LD^−1/2的最小的前k个特征值，然后求出对应的特征向量，得到最后的特征矩阵F,最后对F进行一次传统的聚类（比如K-Means）即可，每一行作为一个k维的样本，共n个样本，用输入的聚类方法进行聚类，聚类维数为k2，得到簇划分C(c₁,c₂,…c_k2)。
   D是对角矩阵，因此D^-1/2就是对角元素取根号再求倒数，因此 D^-1/2 LD^-1/2就是对L的所有元素作如下操作：L_ij/sqrt(D_ii) /sqrt(D_jj), i, j=1…n

附录1 “秩”和“特征值”

n x n的方阵A, 它一定有n个特征值（算上重数）。矩阵的行列式为

如果特征值全部不相等（每个特征值的代数重数为1），那么一定有n个线性无关的特征向量（不相等的特征值对应的特征向量线性无关），该矩阵可以对角化。对于某个特征值λ，如果其有重数，带入特征方程（λI-A）x=0，其基础解系的个数为n-rank（λI-A），称为该特征值的几何重数（几何重数小于等于代数重数）。基础解系线性组合的任何一个结果都可以作为该特征值对应的特征向量，因此特征向量的方向不确定。如果n个特征值都不相同（代数重数全为1），也就是说每个特征值的几何重数都为1，（λI-A）x=0的基础解系只有1个，由它线性组合（拉伸）得到的特征向量的方向是确定的。对于任何一个特征值λ，如果其几何重数等于λ的重数（代数重数），几何重数之和为n(基础解系的总数为n), 则该方阵一定是可以对角化的，把所有n个基础解系排列成矩阵P, 并且P一定可逆(最大线性无关组个数为n，也就说满秩)，P-1AP=E。如果某些特征值的几何重数小于其代数重数，则所有基础解系的数目相加一定是小于n的，无法对角化。
对于实对称矩阵来说，不同特征值对应的特征向量一定也是正交的。而即便是特征值有重数，其几何重数也等于其代数重数，因此实对称矩阵一定可以对角化，要注意的是相同特征值对应的基础解系不一定是正交的。但是这所有n个基础解系一定是线性无关的，因此我们把所有n个基础解系排列成矩阵P, P一定可逆，P^-1AP=E。当然我们也可以把带重数的特征值的基础解系做正交化，这样重新得到n个基础解系，它们相互正交，把它们排列成P ’，P ’为正交矩阵。P’^TAP’=E，也就是说 实对称矩阵一定正交相似于对角矩阵。对称矩阵由于特征值可以为正数，0，负数，从相似对角矩阵就可以看出它不一定是满秩的。正定矩阵所有特征值一定都是大于0的，因此一定是满秩的。
强调文本 强调文本
rank(X_mxn)<= min(m,n)，如果rank(X_mxn)=m，称为行满秩，如果rank(X_mxn)=n，称为列满秩；既是行满秩又是列满秩就只能是方阵了。
rank(X)=rank(X^T)。如果 n>m则rank(X^TX _nxn) =rank(X)=m，因此方阵X^TX一定不满秩，行列式为0，存在0特征根，半正定；如果 n<m则rank(X^TX _nxn)=rank(X)<=n，如果X是列满秩的，则rank(X^TX _nxn)=rank(X) =n，XTX满秩，行列式不为0，无0特征根，正定。

附录2 协方差

对多维随机变量x=[x₁, x₂,…, x_n]^T，我们往往需要计算各维度之间的协方差，这样协方差就组成了一个n×n的矩阵，称为协方差矩阵。协方差矩阵是一个对角矩阵，对角线上的元素是各维度上随机变量的方差。我们定义矩阵内的元素Σ_ij (scaler)为

协方差矩阵（n x n）为

其中，E(X)=[E(x₁),E(x₂)…E(x_n)]^T
样本的协方差矩阵与上面的协方差矩阵相同，只是每个随机变量x_i以m个样本替换了。所有样本可以组成一个m×n的矩阵，这里每行代表一个样本。

c_i表示第i维的随机变量xi的样本集合，而由大数定理可知,随机变量的期望可由样本均值代替E(x_i)=c¯_i 。因此，样本的协方差矩阵(n x n)

numpy默认情况是将每一行作为一个随机变量，因此如果要把一列作为随机变量的话，要设置numpy.cov(X,rowvar=False)
如果每个维度上的随机变量已经进行了标准化，即 E(xi)= c¯_i =0。



与我们前面定义的样本协方差矩阵$\hat{\Sigma }$，只差了1/(m-1)这个系数。所以，我们在PCA、LDA、谱聚类里面，称X^TX为协方差矩阵
两个变量之间的皮尔逊相关系数定义为两个变量之间的协方差和标准差的商（-1 ~ 1）

上面为总体相关系数，当这两个变量分别采用两个样本集来观测时，

其中n为某个样本集中的样本数目

采用numpy计算时，默认的也是把一行作为一个随机变量对应的样本集，如果想把一列作为随机变量，则设置rowvar=False(0), numpy.corrcoef(X,rowvar=0)与样本协方差矩阵类似，返回一个n x n的矩阵，分别是各个特征变量（共n个）之间的相关系数。

附录3 卡方检验

一般用于离散型特征
以下为一个典型的四格卡方检验，我们想知道喝牛奶对感冒发病率有没有影响：

	感冒人数	未感冒人数	合计	感冒率
喝牛奶组	43	96	139	30.94%
不喝牛奶组	28	84	112	25.00%
合计	71	180	251	28.29%

通过简单的统计我们得出喝牛奶组和不喝牛奶组的感冒率为30.94%和25.00%，两者的差别可能是抽样误差导致，也有可能是牛奶对感冒率真的有影响。
为了确定真实原因，我们先假设喝牛奶对感冒发病率是没有影响的，即喝牛奶喝感冒时独立无关的，所以我们可以得出感冒的发病率实际是（43+28）/（43+28+96+84）= 28.29%
所以，理论的四格表应该如下表所示：

	感冒人数	未感冒人数	合计
喝牛奶组	=139*0.2829	=139*(1-0.2829)	139
不喝牛奶组	=112*0.2829	=112*(1-0.2829)	112

即下表（如果喝牛奶和感冒真的是独立无关的，那么四格表里的理论值和实际值差别应该会很小。）：

	感冒人数	未感冒人数	合计
喝牛奶组	39.3231	99.6769	139
不喝牛奶组	31.6848	80.3152	112
合计	71	180	251

卡方检验的计算公式为：

其中，A为实际值，T为理论值。
卡方 = (43 - 39.3231)平方 / 39.3231 + (28 - 31.6848)平方 / 31.6848 + (96 - 99.6769)平方 / 99.6769 + (84 - 80.3152)平方 / 80.3152 = 1.077
这里还需要用到一个自由度的概念，自由度等于V = (行数 - 1) * (列数 - 1)，对四格表，自由度V = 1。对V = 1，相关置信度为5%的临界卡方值是：3.84。即如果卡方值小于3.84，则认为相关置信度小于5%。
显然1.077<3.84，没有达到卡方分布的临界值，所以喝牛奶和感冒独立不相关的假设成立。