繼續教自己 Linear Algebra 。最近在 Stanford 的 Machine Learning 課,學了如何用 Linear Algebra 的方法進行 Principal Component Analysis (PCA) 。以前寫過一篇文,是講如何用 R 進行 PCA ,今次想只用 R 的 Linear Algebra 運算來複製上次的計算結果。的確,這是很無聊的行為,但卻有教育意義。

[code]
#< -
morpho <- read.csv("http://www.chainsawriot.com/data_epb/pca_example.csv", header = TRUE, sep = ",", quote="\"", dec=".",fill = TRUE, comment.char="")[,2:7]
morpho <- as.matrix(morpho)
# dimension of morpho, 28 x 6
m <- dim(morpho)[1]
n <- dim(morpho)[2]
# mean normalization
column.mean <- apply(morpho, 2, mean)
# tedious, can use sweep function
mat.mean <- t(t(matrix(1, m, n)) * column.mean)
# the normalization is only x - mean
adj.morpho <- morpho - mat.mean
# covariance matrix, remember the averaging factor is 1/(m-1), instead of 1/m
# it is the same as cov(morpho)
mat.cov <- (1/(m-1)) *(t(adj.morpho) %*% adj.morpho)
# mat.cov can be covert to correlation matrix using cov2cor function.
# But i want to do it using linear algebra
# conversion of covariance to correlation = cov(x,y) / (sd(x)*sd(y))
column.sd <- apply(morpho, 2, sd)
# colculation the outer product of column.sd. Can be done using
# column.sd %o% column.sd
outer.sd <- column.sd %*% t(column.sd)
mat.cor <- mat.cov / outer.sd
# svd the mat.cor
svd.cor <- svd(mat.cor)
eigen.value <- svd.cor$d
#screeplot
plot(eigen.value, type="b", xlab="Component", ylab="Variances")
abline(h=1, lty=3)

# we are going to take the first two components and compute the two component scores
u.reduce <- svd.cor$u[,1:2]

# For compatibility with score from Princomp()
# sds is the idiosyncractic scaling factor used by Princomp
# not the same as column.sd

sds <- sqrt(diag(mat.cov*(1 - 1/m)))
mat.sds <- t(t(matrix(1, m, n)) * sds)
mat.zscore <- adj.morpho / mat.sds
comp.score <- mat.zscore %*% u.reduce
[/code]

講講每一行在做甚麼。

L2-3:

讀入上次的數據。由於 read.csv 只會回報 data frame ,要用 as.matrix 的方法轉成 matrix 。此 matrix 的 dimension 應是 28 x 6 。

L5-6:

儲存起 morpho 此一 matrix 的 row 數(m,即有幾多 set data) 及 column 數。 (n, 即有幾多個 variable )

L7:

計算 column mean (1x6 vector),以作 normalization / scaling 之用。

L 10-12:

其實只想將 morpho 每一 row 減去 column mean 。我只想出用此方法。需注意的是及後計算的 covariance matrix ,是無需要除以 sd 。

L 15:

計算出 adj.morpho 的 covariance matrix 。
t(adj.morpho) 是 6 x 28 的 matrix ,乘 adj.morpho 是 28 x 6 Matrix ,結果是 6 x 6 matrix 。在 R ,如要計算與 cov() 算出來的 covariance matrix 一樣的結果,是要除 m-1 。我在 Stanford 的 Machine Learning 課是教除以 m 的。 ((為何是除 m-1 ,是因為除 m 的是 biased estimator 。如果 m 是超大的話,減不減那 1 分別不大。))

L 19-23:

其實是可以用 covariance matrix 進行 PCA 的,但是學術界又流行用 correlation matrix ,故此要將 covariance matrix 轉成 correlation matrix 。由於 correlation(x,y) = covariance(x,y) / (sd(x) * sd*(y)) ,可用此一方法轉換。注意的是 22 行,計算了 sd vector 的 outer product 。

L 25-32:

Singular Value Decomposition (svd) 方法計算 Eigen values 和 Loadings (Eigen Vectors) 。再用 Eigen Value 繪製 Scree Plot 。由 Scree plot 可見只有 Component 1 和 Component 2 的 Eigen Values 高於 1 。故此,只取頭兩個 component 。將頭兩個 component 的 loadings (u, Eigen Vector) 存為 u.reduce 。

L38-41:

最後是要計出 Component Score 。在 Stanford 的課,教做的方法是 t(u.reduce) %*% morpho ,即是 loading * x 的 linear combination 。經過不同的方法 scaling ,可得到類似的結果,但都不能 reproduce princomp 的結果。最後的方法是看 princomp 的 source code , (stats:::princomp.default) 才知道 princomp 所用的 scaling 方法,是用一個叫 sds 的物體。但此物又不是 column.sd