如果你的職業定位是數據分析師/計算生物學家，那么不懂PCA、t-SNE的原理就說不過去了吧。跑通軟件沒什么了不起的，網上那么多教程，copy一下就會。關鍵是要懂其數學原理，理解算法的假設，適合解決什么樣的問題。
學習可以高效，但卻沒有捷徑，你終將為自己的思維懶惰和行為懶惰買單。
PCA的原理和普通實現
PCA原理
2019年05月16日
用了這么久的PCA，看了很多人的講解，基本上都是一上來就講協方差矩陣、特征值、特征向量和奇異值分解，其實這對新手是非常不友好的。
PCA只是一種思想，核心就是線性變化，線性代數里的工具只是一種高效的實現PCA的手段，但並不是唯一的工具。
PCA的核心思想：假設我們的數據有D1, D2,···,Dn個維度，PCA就是要構造一個線性變換，PCi = W1D1 + W2D2 +···+ WnDn，Wj就是第j維度在第i個PC中的權重。找PC有先后順序，我們總是先找總方差最大的PC，方差解釋度是統計里最重要的一個概念。然后我們再找與前一個PC線性無關的能解釋最多方差的下一個PC，以此類推，知道得到所有的n個PC。最終的結果就是原先的n個維度通過線性變換，變成了新的n個線性無關的按方差解釋度排序的PC。主成分中的“主principal”針對的就是方差解釋度。
先從一個玩具的案例開始，熟悉的R數據mtcars，取其中1、3列，二維可視化方便些：
mpg disp Mazda RX4 21.0 160 Mazda RX4 Wag 21.0 160 Datsun 710 22.8 108 Hornet 4 Drive 21.4 258 Hornet Sportabout 18.7 360 Valiant 18.1 225
我們先用R的prcomp包跑一下，熟悉一下PCA的輸入輸出：
my_data <- mtcars[,c(1,3)] head(my_data) plot(my_data, main="raw data") my_data <- scale(my_data, center = T, scale = T) plot(my_data, main = "scaled data") # mtcars.pca <- prcomp(mtcars[,c(1:7,10,11)], center = TRUE,scale. = TRUE) mtcars.pca <- prcomp(my_data, center = F,scale. = F) summary(mtcars.pca) plot(mtcars.pca$x, main = "after PCA")
Importance of components:
PC1 PC2
Standard deviation 1.3592 0.39045
Proportion of Variance 0.9238 0.07622
Cumulative Proportion 0.9238 1.00000
PCA的輸入很簡單，就是一個矩陣，行為樣本，列為維度/特征；
這里我把scale這一步單獨列出來了，大家可以看下scale的效果是什么，PCA的原始數據如果沒有scale，PCA的結果就是錯的，因為PCA非常看重方差占比，我們的數據的單位往往千差萬別，必須scale，使得不同維度的方差可比。
PCA的輸出，sdev、rotation、x；sdev標准差是每個PCA的方差解釋度的平方根，肯定是遞減的；rotation旋轉，很形象，PCA只是把原維度做了旋轉而已，rotation是個矩陣里面就是原維度如何通過線性變換，得到PC的；x就是PCA之后原來數據在新的PC下的坐標。
sdev the standard deviations of the principal components (i.e., the square roots of the eigenvalues of the covariance/correlation matrix, though the calculation is actually done with the singular values of the data matrix). rotation the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings. x if retx is true the value of the rotated data (the centred (and scaled if requested) data multiplied by the rotation matrix) is returned. Hence, cov(x) is the diagonal matrix diag(sdev^2).
For the formula method, napredict() is applied to handle the treatment of values omitted by the na.action.
輸入輸出明確了，接下來我們開始做一些個性化的事情：
1. 如果我有新的數據，如何看它在原PCA里的位置？
看下PC是如何得到的，從rotation里我們可以得到線性變換的方程。
PC1 PC2 mpg -0.7071068 0.7071068 disp 0.7071068 0.7071068
PC1 = -0.7071068 * mpg + 0.7071068 * disp
PC2 = 0.7071068 * mpg + 0.7071068 * disp
把新數據代入到這個方程里就可以求新數據的PC score了，在R里，一般用矩陣乘法實現：Newdata %*% rotation
2. 如何找對我們某個PC貢獻最大或top10的維度？取絕對權重最大的top維度集合
在對RNA-seq數據進行PCA分析時非常實用，我們發現了感興趣的PC，往往要找哪些基因貢獻最大，或者直接可視化我們某個基因在前兩個PC中的位置。
3. 如何在原數據里可視化我們的投射/旋轉方法？
藍線就代表了我們原坐標下的PC1方向，如何反向找出藍線的斜率呢？
假設藍線為y=kx，那么向量[1, k]在rotation的作用下必然變為[c, 0]，c的大小和符號是不確定的，但我們有一個確定的0，有點線性代數基礎的就知道如何解k了。
plot(my_data, main = "scaled data") abline(v=0) abline(h=0) abline(coef = c(0, -0.7071068/0.7071068), col="blue") abline(coef = c(0, 0.7071068/0.7071068), col="red")
那么我們就知道了第一個方差最大的方向是y = -0.7071068/0.7071068 * x，第二個方差最大的就是y = 0.7071068/0.7071068.
普通實現
如何用普通方法，不用線性代數和矩陣來求呢？
先明確總方差這個概念：
單一維度的方差好理解，直接求var就行，比如我們scale后，mpg和disp的方差就都為1了。
但是我們的數據PCA之后，PC1的方差就為1.847425，PC2的方差為0.1524486。
Importance of components:
PC1 PC2
Standard deviation 1.3592 0.39045
Proportion of Variance 0.9238 0.07622
Cumulative Proportion 0.9238 1.00000
如何理解總方差呢？其實很簡單，就是所有PC的方差之和。PC1的Proportion of Variance是這么算的：1.847425/(1.847425 + 0.1524486)
現在我們手動來做PCA，先求第一個PC，就是要找藍線y = kx，如何求k？
我們要把所有點投射到藍線上，然后要讓所有的投影點的在藍線這個維度方差最大化。
這里沒辦法了，點太多不可能手算，所以只能用最大似然估計了。
最大似然估計用的mle函數，此函數需要一個似然函數，變得就是A，就是Ax+By+C=0，我們就是不停地試A，然后看什么情況下我們所有的點在此線上投影的方差最大。
最后的結果是：1.009254*x + y = 0，可以看到與原來的線幾乎重合，我們手動算出來的PC和prcomp算出來的一模一樣。
這里的這個似然函數的構造有一點技巧：需要知道怎么算點到線的距離；根據勾股定理算投影長度，需要判斷投影的正負。
PC1定了，其實PC2也就定了，因為只有一條線與PC1正交，就是垂直的線。所以也就不重復算了。
似然函數和最大似然估計不清楚的可以看我的這個文章：似然函數 | 最大似然估計 | likelihood | maximum likelihood estimation | R代碼
##################
nll <- function(A) {
# A=1
B=1; C=0
m <- my_data[,1]
n <- my_data[,2]
if (A!=0) {
raw_sign <- as.integer(n > -m/A)
raw_sign[raw_sign==0] <- -1
} else {
raw_sign <- 1
}
p <- sqrt(((A*m+B*n+C)^2 / (A^2 + B^2)) + m^2 + n^2) * raw_sign
var(p)
}
est <- stats4::mle(minuslog=nll, start=list(A=1), )
summary(est)
plot(my_data, main = "scaled data")
abline(v=0)
abline(h=0)
abline(coef = c(0, -0.7071068/0.7071068), col="blue")
abline(coef = c(0, -1.009254), col="red")
參考：
Principal component analysis - wiki
Principal Component Analysis in R
主成分分析（Principal components analysis）-最大方差解釋
PCA的線性代數實現
自此，我們已經了解了PCA的思想內涵，普通實現。
這個時候再來講PCA的線性代數實現才是合適的，為什么線性代數里的工具能高效的實現PCA？
R的PCA包，prcomp和princomp，The function princomp() uses the spectral decomposition approach. The functions prcomp() and PCA()[FactoMineR] use the singular value decomposition (SVD).
待續~
2019年04月25日
不該先說covariacne matrix協方差矩陣的，此乃后話，先從直覺理解PCA。先看一個數據實例，明顯的兩個維度之間有一個相關性，大部分的方差可以被斜對角的維度解釋，少數的noise則被虛線解釋。
我們的任務就是通過線性變化，找出一個維度能最大的解釋數據的方差；然后繼續找正交的第二個維度，直至最后一個正交維度。
不用矩陣思維我們都可以做，構造一個線性組合，優化方案就是方差最大，然后尋找另一個正交維度。
熟悉了普通解法，我們再來矩陣的解法。
PCA里最重要的就是covariacne matrix協方差矩陣，它完美的包含了我們需要的信息，變量內的方差，變量間的相關程度協方差；
我們想要干什么呢？我們的目的是將covariacne matrix對角化，也就是變量內的方差最大，變量間的相關程度協方差為0；
PCA的靈魂是什么？特征值和特征向量。Ax=λx. A是covariacne matrix，λ是特征值，x是特征向量。這個線性變化意味着什么呢？意味着在經過x處理（線性變換/移動）通過正交矩陣正交化后，A變成了特征值對角矩陣，也就是變量內的方差最大，變量間的相關程度協方差為0。這也解釋了為什么covariacne matrix的特征值其實就是每個PC的variance。按variance排列也就自然而然了。同時也解釋了為什么轉化后的PC彼此線性無關，因為它是正交矩陣特征值構成的是對角矩陣。
特征向量的矩陣表征的就是線性變化，所以原先的數據矩陣在施加x運動后，會投射到新的PC里，這也從直覺上理解了PCA的精髓，線性投射。
PCA本身沒有降維的功能，只是PCA告訴我們了每個PC的variance解釋度，我們選取top 99%的PCs就能解釋絕數據里的大部分信息，所以PCA才有了降維的功能！！！
大誤區：沒有明白特征向量的意義，搞混了特征向量、原數據矩陣和協方差矩陣之間的關系。沒有完全理解對角化和行列式求解之間的聯系。如何直觀理解特征向量矩陣/投影矩陣的由來還是有點難度的！！！
矩陣特征值是對特征向量進行伸縮和旋轉程度的度量，實數是只進行伸縮，虛數是只進行旋轉，復數就是有伸縮有旋轉。我們一般求得的是實數，所以特征值最大表示在這個方向的運動最大。
如何求實對稱矩陣的特征值和特征向量？
看來不把線性代數學通，想徹底搞懂PCA是不可能的了。
講得人很多，但沒有真正講清楚的。從特征分解到協方差矩陣：詳細剖析和實現PCA算法
求協方差矩陣的特征值和特征向量已經很明確了，但對角化的目的，以及為什么協方差矩陣的特征向量具有投影的功能沒講清楚。
協方差矩陣是實對稱矩陣，不具有運動特性，所以不該從運動的角度去解釋。
從定義出發，Ax=cx：A為矩陣，c為特征值，x為特征向量。矩陣A乘以x表示，對向量x進行一次轉換（旋轉或拉伸）（是一種線性轉換），而該轉換的效果為常數c乘以向量x（即只進行拉伸）。我們通常求特征值和特征向量即為求出該矩陣能使哪些向量（當然是特征向量）只發生拉伸，使其發生拉伸的程度如何（特征值大小）。這樣做的意義在於，看清一個矩陣在那些方面能產生最大的效果（power），並根據所產生的每個特征向量（一般研究特征值最大的那幾個）進行分類討論與研究。知乎Rex
talia有個PPT講得非常嚴謹。
不得不驚嘆數學的神奇，一個復雜問題用這么精簡巧妙的方式解決了。
對於基因表達（>30000 gene），表觀（十萬以上enhancer）等超多特征的數據，降維是繞不開的，尤其是在計算樣本或細胞間的關系時，降維能有效的減少計算量。
降維的另一個作用就是可視化，所以分析單細胞數據總能見到PCA、t-SNE和diffusionMap。人最多可以看三維的數據，還得移動看，學術上嚴謹些，只用二維數據。
PCA，以基因表達數據為例，每個基因在樣本間都是在變得，而且有很多是相互相關的，舉個例子，如果兩個/一群基因是完全相關的，那我們完全只需要一個變量/維度就可以解釋樣本里的所有變異。PCA就是通過線性變化，把原來的所有變量（n個）轉化為新的n個PC，這些PC彼此間線性無關，而且可以根據變異的解釋量來排序PC，每個變量對新的PC都有一個簡單的線性轉換公式。
PCA的數學基礎：要有線性代數基礎。PCA和SVD的聯系。
問題：
1.為什么說PCA有降維的功能？明明線性變化后有與變量數相等的PC。
2.PCA適合給什么樣的數據降維？什么情況下用其他降維方法會更好？
3.能不能把PCA的大致計算過程寫下來？
4.PCA的應用，我有一個很感興趣的維度，請把對該維度貢獻最大的變量找出來；我有一個新的樣品，能否給我一個新的PC打分
5.variance，如何理解特征向量表征了變異比例？
6.如何直觀理解協方差，covariance？
7.PCA里是如何利用特征值和特征向量的？
8.如何理解“主成分分析就是把協方差矩陣做一個奇異值分解，求出最大的奇異值的特征方向”？
9.PCA和ICA有什么關系？
10.特征值和特征向量單獨表示很好理解，如何寫成對角矩陣形式？矩陣和向量基本運算要熟練。這里糾結了好久
11.PCA在圖像壓縮和圖像識別方面也有應用。
12. 為什么特征向量的矩陣可以用於空間轉換？
學習是一個思維發散和收斂的過程，不斷深入，不斷回歸。
為了防止我又忘了我的目的，以下梳理一下我的任務：
1. 如何對RNA-seq的數據做PCA分析？
2. PCA原理；如何做到重新分配variance的？
3. 如何理解得到的PC；
4. 有個感興趣的PC，如何提取出關鍵基因；
5. 如何把一個新的樣本投射到已經建立的PCA坐標里，打分預測？
好了，開始發散，先到線性代數里走一圈。
PCA計算的起點是協方差矩陣，方差描述的是一群點偏離均值的程度（偏離度（平方）的均值/期望），標准差是方差的標准化，具有與原數據一樣的單位/量綱。協方差描述的是兩個變量之間的協同關系（兩個變量與其各自均值的協變度），是相關性系數的前體（除以兩個變量的標准差后就是相關性系數）。
協方差公式簡單，但卻有許多驚奇的特性，它居然可以描述兩個變量的協同變化程度。協方差矩陣更是計算生物學的利器（協方差矩陣有什么意義？），很多工具分析的起點。
問題：為什么有了標准差還要方差，兩者不是重復嗎？為什么方差的度量選擇了平方？為什么描述性統計的兩個核心指標是均值和方差，而不是均值和標准差？參考鏈接， 可導。
題外話：為什么樣本的方差要除以n-1而不是n？參考鏈接：無偏估計。
接下來就是求協方差矩陣covariance matrix的特征值和特征向量了eigenvalue and eigenvector。
直觀理解，A就是我們的矩陣，它可以看作是運動，v是特征向量，特征值和特征向量表征的是一個矩陣的特征（如何理解矩陣特征值？），在可視化工具里，給定一個矩陣，我們可以找到它的一些特征值和特征向量，意思就是特征向量在運動后（乘以A），方向不變，長度變了，變化的長度就是λ。
線性代數就是解方程，前面解的方程是Ax=b，現在解的方程是Ax=λx.
說多了感覺脫離實際了，show me your code
需要注意的細節：scale會影響結果
head(mtcars) res <- prcomp(mtcars, scale. = F, center = F, retx = T) # res2 <- princomp(mtcars) plot(res$x[,1:2]) res$x[1:5,1:5] res$sdev
res$rotation[1:5,1:5] recover <- as.matrix(mtcars) %*% as.matrix(res$rotation) recover[1:5,1:5]
> head(mtcars)
mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
一般的主成分分析返回的信息有： sdev標准差，平方后就是每一個PC所能解釋的方差，也就是特征向量；rotation對應了特征向量的矩陣；x就是最終的投射后的矩陣，就是原矩陣乘以特征向量rotation。
有個無法還原的就是Ax=λx，協方差矩陣乘以特征向量矩陣不等於右邊，不知道哪里出錯了。需要看下源碼
The Mathematics Behind Principal Component Analysis 數學基礎，步驟很清楚，怎么做講得很清楚，但沒講出緣由，為什么可以這么做
Principal Component Analysis in R: prcomp vs princomp R代碼
Principal Component Analysis (PCA) - THE MATH YOU SHOULD KNOW! 數學部分完全不懂
Principal Components Analysis: Theory and Application to Gene Expression Data Analysis 正統paper
再談協方差矩陣之主成分分析 - 傳神的講了協方差矩陣的神奇作用
PCA的數學原理 非常不錯了
總結一下：做PCA，首先我們會構建特征/變量的協方差矩陣，然后求其特征值和特征向量，top特征向量構建的矩陣就能將原始數據投射到新的相互獨立的且排好序的PC空間里。
新的問題：
1. 為什么這樣就可以投射到新的PC了？特征向量可以理解為對原來的所有維度進行重構了，特征向量相當於是對應的線性權重，這就是PCA核心的線性變換過程。為什么特征向量就是權重，如何理解？
2. 為什么求解的特征值就是對應特征向量的方差，這些是怎么聯系起來的？必須理解協方差矩陣的特征值和特征向量意味着什么！
深刻理解到，線性代數和概率論不學好，搞統計研究就會舉步維艱。
Dimensionality reduction for visualizing single-cell data 需要總結一下
番外篇：
1. PCA與聚類的關系，PCA不是聚類，它只是降維，只是在RNA-seq當中，好的replicates往往會聚在一起，才會誤以為PCA可以做聚類。PCA是降維，是聚類的准備工作，最常見的聚類是k-means聚類，為了降低計算復雜度，我們可以在PCA的結果里做聚類。
2. 協方差矩陣為什么比相關性矩陣用途更廣？因為協方差里包含了兩種信息，變量間的相關性以及變量內部的方差，相關性標准化后丟失了太多的信息，所以用得少了。
例子：以小學生基本生理屬性為案例，分享下R語言的具體實現，分別選取身高（x1）、體重(x2)、胸圍(x3)和坐高(x4)。具體如下：
> student<- data.frame( + x1=c(148,139,160,149,159,142,153,150,151), + x2=c(41 ,34 , 49 ,36 ,45 ,31 ,43 ,43, 42), + x3=c(72 ,71 , 77 ,67 ,80 ,66 ,76 ,77,77), + x4=c(78 ,76 , 86 ,79 ,86 ,76 ,83 ,79 ,80) + ) > student.pr <- princomp(student,cor=TRUE) > summary(student.pr,loadings=TRUE) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Standard deviation 1.884603 0.57380073 0.30944099 0.152548760 Proportion of Variance 0.887932 0.08231182 0.02393843 0.005817781 Cumulative Proportion 0.887932 0.97024379 0.99418222 1.000000000 Loadings: Comp.1 Comp.2 Comp.3 Comp.4 x1 -0.510 0.436 -0.139 0.728 x2 -0.513 -0.172 -0.741 -0.398 x3 -0.473 -0.761 0.396 0.201 x4 -0.504 0.448 0.524 -0.520 > screeplot(student.pr,type="lines")
standard deviation：標准差
Proportion of Variance：方差的占比
Cumulative Proportion：累計貢獻率
由上圖可見四項指標做分析后，給出了4個成分，他們的重要性分別為：0.887932、0.08231182、0.02393843、0.005817781，累積貢獻為：0.887932、0.97024379、0.99418222 1.000000000各個成分的碎石圖也如上，可見成份1和成份2的累積貢獻已經達到95%，因此采用這兩個成份便可充分解釋學生的基本信息。
可以算出 Z1 和 Z2的公式
> temp<-predict(student.pr) > plot(temp[,1:2])
參考鏈接：R語言與數據分析之五：主成分分析（還有很多系列，慢慢看）
主成分分析（Principal Component Analysis，PCA）的目標是用一組較少的不相關的變量代替大量相關變量，同時盡可能保留原始變量的信息，推導所得的變量就成為主成分，是原始變量的線性組合。也就是將N個變量（N維），通過線性組合，降維到K個綜合變量（K維，K <N）來歸納性解釋某一個現象。
先舉個簡單例子幫助理解吧：
某籃球俱樂部有40名男同學，同學之間各個指標存在或大或小的差異，包括身高、體重、視力、百米速度、肺活量、每天練球時長、睡眠時間等指標。在季度選拔中，40名同學進球得分數（成績）存在差異，那么這些指標是否均與成績相關？或者相關性有多大？
此時可能就要用到PCA，分析方法簡述如下：
1選擇初始變量
比如以上7個指標作為變量（a1 - a7），40名同學作為樣本。
2對原始數據矩陣進行標准化，做相關系數矩陣
（1）原始數據矩陣：每行為40名男同學的各項指標值，每列為各項指標在40名同學中的體現；
（2）因為各個指標度量單位不同，取值范圍不同，不宜直接由協方差矩陣出發，因此選擇相關系數矩陣。
3計算特征值及相應特征向量
4判斷主成分的個數
最常見的方式是根據特征值判斷，一般選擇特征值大於1的變量數作為PC個數，假如，此項分析中特征值大於1的有兩個，則最終可以有2個主成分（具體主成分個數可以根據實際研究調整）。
5得到主成分表達式
通過分析得到第一主成分（PC1）和第二主成分（PC2），假設表達式如：（a1*表示a1標准化后的數值）
6結合數據的實際意義開展分析
PC1比PC2更能解釋樣本間差異的原因（如下圖中橫縱軸的百分比）。PC1的線性組合中a1、a2、a3、a6貢獻度較大（前面的系數較大），PC2的線性組合中a5貢獻度較大。以PC1為橫軸，各個樣本根據成績大小有明顯區分，說明以上7個指標中，身高、體重、視力、每天練球時長這4個指標與同學的成績相關性更強。為了迎合線性組合的概念，應該找個更合適的詞語來綜合描述身高、體重、視力、每天練球時間這4個指標，以涵蓋這4個指標的意義。理解了以上的概念，再將主成分分析用於RNA-seq也很容易理解。
參考鏈接：
StatQuest: Principal Component Analysis (PCA) clearly explained (2015) 必須一看
How to perform dimensionality reduction with PCA in R 具體實現
主成分分析（Principal components analysis）-最大方差解釋
Principal component analysis for clustering gene expression data
主成分分析 - stanford
主成分分析法 - 智庫
主成分分析（Principal Component Analysis）原理
主成分分析及R語言案例 - 文庫
主成分分析法的原理應用及計算步驟 - 文庫
【機器學習算法實現】主成分分析(PCA)——基於python+numpy
主成分分析（PCA）原理詳解（推薦）
主成分分析-生物信息
主成份分析（PCA）在生物芯片樣本篩選中的應用及在R語言中的實現