PCA、SVD與譜聚類

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可以通過下式得到:

SVD

SVD也是對矩陣進行分解，但是和特徵分解不同，SVD並不要求要分解的矩陣為方陣。假設我們的矩陣A是一個m×n的矩陣，那麼我們定義矩陣A的SVD為： 其中U是一個m×m的矩陣，Σ是一個m×n的矩陣，除了主對角線上的元素以外全為0，主對角線上的每個元素都稱為奇異值，V是一個n×n的矩陣。U和V都是酉矩陣（實數域上也就是正交矩陣），即滿足U^T

U=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，實現了樣本壓縮。

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)。

譜聚類

主要思想：把所有的資料看做空間當中的點，這些點之間可以用邊連線起來。距離較遠的兩個點之間的邊權重值比較低，較近的點權重值比較高，目的是：進行切圖，讓切圖後不同的子圖之間邊權重和儘可能低，而子圖內部的邊權重和儘可能的高，從而達到聚類的目的。

無向權重圖

由點集合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內的點其所有邊的權重和

相似矩陣

基本思想：距離較遠的兩個點之間的權重值較低，距離較近的兩個點之間權重值較高。 構造相似矩陣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核函式。

拉普拉斯矩陣

L=D-W（D是對稱陣，W也是對稱陣）（L也是對稱陣），由於L是對稱陣，所以特徵值都是實數。對於任意向量f有： 拉普拉斯矩陣是半正定的，且對應的n個實數特徵值0=λ1<=λ2<=λn，對應的特徵項向量的模長是1。

無向圖切圖

對於無向圖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。 與我們前面定義的樣本協方差矩陣，只差了1/(m-1)這個係數。所以，我們在PCA、LDA、譜聚類裡面，稱X^TX為協方差矩陣 兩個變數之間的皮爾遜相關係數定義為兩個變數之間的協方差和標準差的商（-1 ~ 1） 上面為總體相關係數，當這兩個變數分別採用兩個樣本集來觀測時， 其中n為某個樣本集中的樣本數目

採用numpy計算時，預設的也是把一行作為一個隨機變數對應的樣本集，如果想把一列作為隨機變數，則設定rowvar=False(0), numpy.corrcoef(X,rowvar=0)與樣本協方差矩陣類似，返回一個n x n的矩陣，分別是各個特徵變數（共n個）之間的相關係數。