因子分析、主成分分析（PCA）、獨立成分分析（ICA）——斯坦福CS229機器學習個人總結（六）

因子分析是一種數據簡化技術，是一種數據的降維方法。
因子分子可以從原始高維數據中，挖掘出仍然能表現衆多原始變量主要信息的低維數據。此低維數據可以通過高斯分佈、線性變換、誤差擾動生成原始數據。
因子分析基於一種概率模型，使用EM算法來估計參數。

主成分分析（PCA）也是一種特徵降維的方法。
學習理論中，特徵選擇是要剔除與標籤無關的特徵，比如「汽車的顏色」與「汽車的速度」無關；
PCA中要處理與標籤有關、但是存在噪聲或者冗餘的特徵，比如在一個汽車樣本中，「千米/小時」與「英里/小時」中有一個冗餘了。
PCA的方法比較直接，只要計算特徵向量就可以降維了。

獨立成分分析（ICA）是一種主元分解的方法。
其基本思想是從一組混合的觀測信號中分離出獨立信號。比如在一個大房間裏，很多人同時在說話，樣本是這個房間裏各個位置的一段錄音，ICA可以從這些混合的錄音中分離出每個人獨立的說話的聲音。
ICA認爲觀測信號是若干個統計獨立的分量的線性組合，ICA要做的是一個解混過程。

因爲因子分析、PCA、ICA都是對數據的處理方法，就放在這同一份總結裏了。

1、因子分析（Factor analysis）

1.1、因子分析的直觀理解

因子分析認爲高維樣本點實際上是由低維樣本點經過高斯分佈、線性變換、誤差擾動生成的。讓我們來看一個簡單例子，對低維數據如何生成高維數據有一個直觀理解。

假設我們有m=5個2維原始樣本點如下：

圖一

那麼按照因子分析的做法，原始數據可以由以下過程生成：
①在一個低維空間（此處是1維）中，存在着由高斯分佈生成的 m 個點 z(i) ， z(i) ~ N(0,I) ：

圖二
②使用某個 Λ=(a,bT) 將1維的 z 映射到2維的空間中：

圖三
③加上 μ(μ1,μ2)T ，讓直線過 μ ——實際上是將樣本點橫座標加 μ1 ，縱座標加 μ2 ：

圖四
④對直線上的點做一定的擾動，其擾動爲 ε ~ N(0,ψ) ：

圖五
黑點就是圖一中的原始數據。

1.2、因子分析的一般過程

因子分析認爲m個n維特徵的訓練樣例 (x(1),x(2),⋯,x(m)) 的產生過程如下：
①在一個 k 維空間中，按照多元高斯分佈生成m個 z(i) （ k 維向量， k<n ），即
z(i) ~ N(0,I)
②存在一個變換矩陣 Λ∈Rn∗k ，將 z(i) 映射到 n 維空間中，即
Λz(i)
③將 Λz(i) （ n 維）加上一個均值 μ （ n 維），即
μ+Λz(i)
④對每個點加上符合多元高斯分佈的擾動 ε ~ N(0,ψ) （ n 維向量），即
x(i)=μ+Λz(i)+ε

1.3、因子分析模型

模型與參數概述

由上面的分析，我們定義因子分析的模型爲：

z ~ N(0,I)
ε ~ N(0,ψ)
x=μ+Λz+ε(1)
其中 z 和 ε 是相互獨立的。並且由上面的分析過程，我們可以直觀地感受到我們的 參數是 μ∈Rn 、 Λ∈Rn∗k 、 ψ∈Rn∗n 。

另一個等價的假設是， (x,z) 聯合分佈如下，其中 z∈Rk 是一個隱藏隨機變量：

x∣z ~ N(μ+Λz,ψ) (2)
這個假設會在使用EM算法求解因子分析參數，E步中迭代 Q 分佈的時候用到。

接下來的課程，是使用高斯模型的矩陣表示法來對模型進行分析。矩陣表示法認爲 z 與 x 聯合符合多元高斯分佈，即：

[zx] ~ N(μzx,Σ)
多元高斯分佈的原始模型是：

f(x)=12πk|Σ|−−−−−−√exp(−12(x−μ)TΣ−1(x−μ))(3)

其中 x 是 k 維向量， μ 是 k 維向量， Σ 是 k∗k 協方差矩陣。
很明顯在多元高斯分佈模型下，參數是 μzx,Σ ——它們是由 x,z 的聯合分佈生成的，所以我們可以用我們的原始參數 μ,Λ,ψ 來表示 μzx,Σ ，求得 x 的邊緣分佈，再把相關參數帶入式（3），這就得到了關於我們參數的概率分佈，然後就可以通過最大似然估計來求取我們的參數。

求取 μzx

μzx 是 x,z 聯合分佈的期望值（期望的定義：所有結果*相應概率的總和）：

μzx=E[zx]=[E(z)E(x)](4)

由 z ~ N(0,I) 我們可以簡單獲得 E(z)=0 。
類似地由 ε ~ N(0,ψ) ， x=μ+Λz+ε ， μ 是一個常數，我們有：

E[x]=E[μ+Λz+ε]=E[μ]+ΛE[z]+E[ε]=μ+0+0=μ(5)

所以：

μzx=[0⃗ μ](6)

求取 Σ

Σ 是 x,z 聯合分佈的協方差矩陣。
方差，度量隨機變量與期望之間的偏離程度，定義如下：

Var(X)=E((X−E(X))2)=E(X2)−(E(X)2)(7)

協方差，兩個變量總體誤差的期望，定義如下：

Cov(X,Y)=E((X−E(X))(Y−E(Y)))(8)

協方差、方差、期望之間的一些相互關係如下：

Cov(X,X)=Cov(X)=Var(X)=E(XXT)=σ2(9)

下面開始求取 Σ 。

Σ=Cov[zx]=[ΣzzΣxzΣzxΣxx]=E[(z−E(z))(z−E(z))T(x−E(x))(z−E(z))T(z−E(z))(x−E(x))T(x−E(x))(x−E(x))T](10)

由 z ~ N(0,I) ，可以簡單得到：

Σzz=Cov(z)=σ2=I(11)

由 ε ~ N(0,ψ) ， x=μ+Λz+ε ， E(x)=μ ，並且 z 和 ε 是相互獨立，有：

Σzx=E[(z−E(z))(x−E(x))T]=E[(z−0)(μ+Λz+ε−μ)T]=E[zzT]ΛT+E[zεT]=IΛT+0=ΛT(12)

類似地，我們可以得到：

Σxx=E[(x−E(x))(x−E(x))T]=E[(μ+Λz+ε−μ)(μ+Λz+ε−μ)T]=ΛE[zzT]ΛT+E[εεT]=ΛIΛT+ψ=ΛΛT+ψ(13)

用最大似然估計法求解參數

經過上面的步驟，我們就把 μzx,Σ 用我們的參數 μ,Λ,ψ 表示出來了：

[zx] ~ N(μzx,Σ) ~ N([0⃗ μ],[IΛΛTΛΛT+ψ])
然後我們可以求得 x 的邊緣分佈：
x ~ N(μ,ΛΛT+ψ)
因此，給定一個訓練集 {x(i);i=1,2,⋯,m} ，把參數帶入式（3），我們可以寫出下面的似然函數：

l(μ,Λ,ψ)=log∏i=1m12πn∣∣ΛΛT+ψ∣∣−−−−−−−−−−−√exp(−12(x(i)−μ)T(ΛΛT+ψ)−1(x(i)−μ))(14)

對此似然函數做最大似然估計，就能求得我們的參數。

1.4、因子分析的EM求解

可以感受到，直接對這個似然函數求解是很困難的，在這個情況下，用EM算法就登場了——當一個似然函數難以直接求解其最大值的時候，可以通過EM算法不斷建立下界、最大化下界的方式不斷逼近該似然函數真實的最大值，當EM算法收斂，我們就認爲已經求得了此最大值。

E-step

對於EM算法的E-step，我們有：

Qi(z(i)):=p(z(i)∣x(i);μ,Λ,ψ)(15)

進一步地：

Qi(z(i))=12πk∣∣Σz(i)∣x(i)∣∣−−−−−−−−−−√exp(−12(z(i)−μz(i)∣x(i))TΣ−1z(i)∣x(i)(z(i)−μz(i)∣x(i)))(16)

其中：

μz(i)∣x(i)Σz(i)∣x(i)=ΛT(ΛΛT+ψ)−1(x(i)−μ)=I−ΛT(ΛΛT+ψ)−1Λ(17)

μz(i)∣x(i),Σz(i)∣x(i) 是講義與課上直接給出的，這裏也不進行推導。

M-step

在M-step中，我們需要最大化如下公式來求取參數 μ,Λ,ψ ：

∑i=1m∫z(i)Qi(z(i))logp(x(i),z(i);μ,Λ,ψ)Qi(z(i))dz(i)(18)

視爲期望，打開log

在這裏，因爲 z 是連續的，所以使用積分；如果是離散的，則使用累加。
並且，積分部分可以當成 z 服從 Q 分佈時，函數 logp(x(i),z(i);μ,Λ,ψ)Qi(z(i)) 的期望，這裏將會用E表示，省略 z(i) ~ Qi 的下標；對於函數中 x,z 的聯合分佈，我們可以用貝葉斯公式把它打開 p(x,z)=p(x∣x)p(z) ；爲了方便計算我們還要把log函數打開——經過這些分析，我們有如下推導：

∑i=1m∫z(i)Qi(z(i))logp(x(i),z(i);μ,Λ,ψ)Qi(z(i))dz(i)=∑i=1mE[logp(x(i),z(i);μ,Λ,ψ)Qi(z(i))]=∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)p(z(i))Qi(z(i))]=∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)+logp(z(i))−logQi(z(i))](19)

去掉無關項後帶入具體分佈

這就比較清爽了，然後，記住我們的目標是求得參數 μ,Λ,ψ ，但是它們不能一起求解，所以下面以參數 Λ 爲例，對公式進行求解——在式（19）中，對參數 Λ 求偏導。另外式（19）中的 p(z(i) 與 Qi(z(i)) 與 Λ 無關，可以忽略掉，所以實際上就是對下式求偏導：

∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)](20)

在對式（20）求偏導之前，還可以對其進行一些處理——由式（2），並且 x∣z 服從多元高斯分佈，所以有：

∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)]=∑i=1mE[log12πn|ψ|−−−−−−√exp(−12(x(i)−(μ+Λz(i)))Tψ−1(x(i)−(μ+Λz(i))))]=∑i=1mE[−12log|ψ|−n2log(2π)−12(x(i)−μ−Λz(i))Tψ−1(x(i)−μ−Λz(i)))](21)

去掉無關項後求偏導

同樣地，我們的目標是與 Λ 有關的項，所以忽略掉前面的無關項之後，我們實際上是對下式求偏導並求解：

∇Λ∑i=1mE[−12(x(i)−μ−Λz(i))Tψ−1(x(i)−μ−Λz(i)))]=∑i=1m∇Λ−E[12((x(i)T−μT)ψ−1A−z(i)TΛTψ−1B)((x(i)−μ)C−Λz(i)D)]=∑i=1m∇Λ−E[12(AC−AD−BC+BD)](22)

打開後我們可以發現， AC 這一項是與 Λ 無關的，把這一項忽略掉，所以由式（22）繼續推導有：

∑i=1m∇Λ−E[12(−AD−BC+BD)]=∑i=1m∇Λ−E[12(−(x(i)T−μT)ψ−1Λz(i)E−z(i)TΛTψ−1(x(i)−μ)F+z(i)TΛTψ−1Λz(i))](23)

因爲期望是一個常數，又因爲 a=tr(a) ，所以可以直接對上式求跡；
因爲 trA=trAT ，可以對E求轉置，又因爲對角矩陣的轉置是它本身—— (ψ−1)T=ψ−1 ，所以有 trE=trET=trF ，對式（23）繼續推導有：

∑i=1m∇Λ−E[12(−(x(i)T−μT)ψ−1Λz(i)E−z(i)TΛTψ−1(x(i)−μ)F+z(i)TΛTψ−1Λz(i))]=∑i=1m∇Λ−E[tr12(−z(i)TΛTψ−1(x(i)−μ)ET−z(i)TΛTψ−1(x(i)−μ)F+z(i)TΛTψ−1Λz(i))]=∑i=1m∇ΛE[−tr12z(i)TΛTψ−1Λz(i)+trz(i)TΛTψ−1(x(i)−μ)](24)

然後利用 trAB=trBA ，把式（25）中的 z(i)T 放到它們自己的後面，再把求導切換到期望裏面——求導是針對 Λ ，期望是針對 z(i) ，所以是可以切換的：

∑i=1m∇ΛE[−tr12z(i)TΛTψ−1Λz(i)+trz(i)TΛTψ−1(x(i)−μ)]=∑i=1m∇ΛE[−tr12ΛTψ−1Λz(i)z(i)T+trΛTψ−1(x(i)−μ)z(i)T]=∑i=1m(E[−∇Λtr12ΛTψ−1Λz(i)z(i)T]+E[∇ΛtrΛTψ−1(x(i)−μ)z(i)T])(25)

對於式（25），先用矩陣的跡的性質 ∇ATf(A)=(∇Af(A))T 處理一下：

∑i=1m(E[−(∇ΛTtr12ΛTψ−1Λz(i)z(i)T)T]+E[(∇ΛTtrΛTψ−1(x(i)−μ)z(i)T)T])(26)

對式（26）的第一項使用 ∇AtrABATC=CAB+CTABT 的性質，對第二項使用 ∇AtrAB=BT 的性質：

∑i=1m⎛⎝⎜⎜E⎡⎣⎢⎢−⎛⎝⎜∇ΛTAtr12ΛTAψ−1BΛATz(i)z(i)TC⎞⎠⎟T⎤⎦⎥⎥+E⎡⎣⎢⎢⎛⎝⎜∇ΛTAtrΛTAψ−1(x(i)−μ)z(i)TB⎞⎠⎟T⎤⎦⎥⎥⎞⎠⎟⎟=∑i=1mE((−12z(i)z(i)TΛTψ−1−12(z(i)z(i)T)TΛT(ψ−1)T)T+((ψ−1(x(i)−μ)z(i)T)T)T)=∑i=1mE((−z(i)z(i)TΛTψ−1)T+ψ−1(x(i)−μ)z(i)T)=∑i=1mE(−ψ−1Λz(i)z(i)T+ψ−1(x(i)−μ)z(i)T)(27)

令式（27）=0，並化簡，就可以求得參數 Λ ：

⟹⟹⟹∑i=1mE(−ψ−1Λz(i)z(i)T+ψ−1(x(i)−μ)z(i)T)=0∑i=1m−ψ−1ΛE[z(i)z(i)T]+∑i=1m−ψ−1(x(i)−μ)E[z(i)T]=0∑i=1mΛE[z(i)z(i)T]=∑i=1m(x(i)−μ)E[z(i)T]Λ=(∑i=1m(x(i)−μ)E[z(i)T])(∑i=1mE[z(i)z(i)T])−1(28)

我們發現，這裏的公式與線性迴歸中最小二乘法的矩陣形式相似。
相似原因：在因子分析中， x 是 z 的線性函數，在E-step中給出 z 的 Q 分佈之後，在M-tep中尋找 x 與 z 的映射關係 Λ ；在線性迴歸的最小二乘中，也是尋找 x 與 y 的線性關係。
不同之處：最小二乘只用到了 z 的最優估計，因子分析還用到了 z(i)z(i)T 的估計。

對於參數 Λ ，這裏還有未知數 E[z(i)T] 與 E[z(i)z(i)T] ，並且此處的期望是在 z(i) 服從 Qi 前提下計算的，所以對於前者，通過式（17）我們有：

E[z(i)T]=μTz(i)∣x(i)(29)

對於後者，由式（7）~式（9）方差與協方差的性質，我們有：

Cov(z)⟹E[z(i)z(i)T]=E(zzT)−E(z)E(zT)=E(z(i))E(z(i)T)+Cov(z)=μz(i)∣x(i)μTz(i)∣x(i)+Σz(i)∣x(i)(30)

注意這裏的 E[z(i)z(i)T] 不僅僅等於 E(z(i))E(z(i)T) ，後面還有加上一個後驗概率 p(z∣x) 協方差，要特別注意。

到這裏，我們就可以把參數 Λ 的最終形式給出來了：

Λ=(∑i=1m(x(i)−μ)μTz(i)∣x(i))(∑i=1mμz(i)∣x(i)μTz(i)∣x(i)+Σz(i)∣x(i))−1(31)

另外對於其他兩個參數 Λ,ψ ，使用相同的方法可以求得，這裏直接給出結果
：

μ=1m∑i=1mx(i)(32)

ϕ=1m∑i=1mx(i)x(i)T−x(i)μTz(i)∣x(i)ΛT−Λμz(i)∣x(i)x(i)T+Λμz(i)∣x(i)μTz(i)∣x(i)+Σz(i)∣x(i)ΛT(33)

注意這裏的參數是 ϕ 不是 ψ ，得到 ϕ 之後還需要將 ψii=ϕii ，因爲 ϕ 不是對角矩陣，所以只需要取 ϕ 對角線上的值即可。

2、主成分分析（Principal component analysis，簡稱PCA

PCA的意義

PCA技術的一大好處是對數據進行降維的處理。我們可以對新求出的「主元」向量的重要性進行排序，根據需要取前面最重要的部分，將後面的維數省去，可以達到降維從而簡化模型或是對數據進行壓縮的效果。同時最大程度的保持了原有數據的信息。

PCA將n個特徵降維到k個，可以用來進行數據壓縮，如果100維的向量最後可以用10維來表示，那麼壓縮率爲90%。同樣圖像處理領域的KL變換使用PCA做圖像壓縮。但PCA要保證降維後，還要保證數據的特性損失最小。

預處理

運行PCA算法之前，數據一般要進行預處理。預處理步驟如下：
①令 μ=1m∑mi=1x(i)
②用 (x(i)−μ) 代替 x(i)
③令 σ2=1m∑mi=1(x(i)j)2
④用 x(i)j/σ2 代替 x(i)j
步驟①與②將數據的均值變爲0；
步驟③與④將數據每個維度的方差變爲1，使每個維度都在同一個維度下被度量。
如果已知數據均值爲0，或者數據已在同樣一個維度下，就無需進行上面的步驟。

最大方差理論解釋PCA

Ng在課上說，PCA有9到10種解釋方法，這裏僅用其中的最大方差理論來解釋PCA
在信號處理中認爲信號具有較大的方差，噪聲有較小的方差，信噪比就是信號與噪聲的方差比，越大越好。如前面的圖，樣本在橫軸上的投影方差較大，在縱軸上的投影方差較小，那麼認爲縱軸上的投影是由噪聲引起的。
因此我們認爲，最好的k維特徵是將n維樣本點轉換爲k維後，每一維上的樣本方差都很大。

下面舉例來說明如何尋找主方向。

圖六
圖六左圖是經過預處理的樣本點，均值爲0，特徵方差歸一；
圖六的中圖和右圖是將樣本投影到某個維度上的表示， 該維度用一條過原點的直線表示，前處理的過程實質是將原點移到樣本點的中心點。

圖七
圖七中，原點到樣本在直線上的投影的距離 x(i)Tu 既是方差；樣本點到直線上的距離是平方誤差。
我們可以直觀感受到，圖六中間的圖方差和是最大的，平方誤差是最小的。下面給出用最大方差理論尋找該最佳方向的公式定義。

形式化

設 x(i) 是數據集中的點， u 是要求解的單位向量，那麼方差最大化可以形式化爲最大化下式：

1m∑i=1m(x(i)Tu)2=1m∑i=1muTx(i)x(i)Tu=uT(1m∑i=1mx(i)x(i)T)u(34)

因爲歸一化後的數據，投影值的均值也爲0，所以在方差計算中直接平方。
同時這個式子還有一個約束條件，即 ∥u∥2=1 。熟悉的最大化某個帶約束的函數，引入拉格朗日乘子來求解：

l=uT(1m∑i=1mx(i)x(i)T)u−λ(∥u∥2−1)=uTΣu−λ(uTu−1)(35)

對 u 求導：

∇ul=∇u(uTΣu−λ(uTu−1))=∇utr(uTΣu)−λ∇utr(uTu)=(∇uTtr(uTΣu))T−(λ∇uTtr(uTu))T=((Σu)T)T−(λuT)T=Σu−λu(36)

這裏主要運用了 ∇ATf(A)=(∇A