本文讲解SVD和PCA
Singular Value Decomposition
abbr. SVD
The Singular Value Decomposition is a highlight of linear algebra.
A
is any
” A is diagonalized”
The singular vectors
v1…vr
are in the row space of A. The singular values
σ1…σr
are all positives number. When the
v′s
and
u′s
go into the columns of
V
and
proof:
v
是
那么
v
在
ATAvi=σ2ivi
ATA
一定是半正定的(
XTATAX≥0
)
所以存在
r
个eigenvalue 大于0 (r为rank)
row space
(proof: row space是
ATy
的集合与
nullspace
中的
x
的dot product
V
在
目标是
AV=UΣ
VT=V−1
A=UΣVT
ATAvi=σ2ivi
vTiATAvi=σ2ivTivi
∥Avi∥2=σ2i
so that
∥Avi∥=σi
AATAvi=σ2iAvi
gives
ui=Avi/σi
as a unit eigenvector of
AAT
ui
在A的column space上,所以存在
r
个正交基,与
将SVD用于降维:
Am×nVn×r=Um×rΣr×r
就是将
A
投影到
Principal Component Analysis
abbr. PCA
SVD是PCA的一种方法,还有一种常用的方法就是将SVD降维中所用到的
V
换成了
PCA的目的是使投影后的矩阵各向量之间的协方差为0,既去除了各个向量之间的线性关系,并且使得投影后的矩阵的在各个方向上的方差达到最大化。可以参考主成分分分析这篇文章所讲的最大方差理论。简单的对于分类来说,方差越大那么区分度越高。
算法步骤:
- 求矩阵 Am×n 中的每一列的均值,可以得到一个 n×1 的vector m 。其中A的每一行代表一个数据点,而每一列就是一个feature;
- 对
A 的每一行都减去 m ; - 求
A 的转置的协方差矩阵 C ; - 求
C 的eigenvalue和eigenvector; - 根据eigenvalue对eigenvector进行排序,可根据需求舍弃较小的eigenvalue所对应的eigenvector;
- 将eigenvector所组成的矩阵进行转置(
r×n
)与转置的
A
(
n×m )进行dot product,将结果转置既为所求结果。
因为用到了协方差矩阵,协方差矩阵中每一项度量的是两个随机变量的相关性,
Σij=cov(Xi,Xj)=E[(Xi−ui)(Xj−uj)T]
,协方差定义中的
Xi
是行向量,对应
Am×n
中的一列,既对一个random variable的一组观测值。所以需要对
A
进行转置,使得A的行向量成为feature。
假设投影后的矩阵为
Y
,那么
那么
cov(YT)=E((YT)(YT)T)
=E((PTAT)(PTAT)T)
=E(PT(AT(AT)T)P)
=PTcov(AT)P
由于
cov(AT)
是一个对称矩阵,所以可以将其对角化,而对角化所用的矩阵就是
cov(AT)的
eigenvetor所组成的矩阵(spectral theorem)。所以
pi
就是
cov(AT)
的eigenvector。在推导的过程中,要注意A是已经减去了均值的,要不不能这么推,这个过程反推也是有意义的,也能说明第一步为什么不需要减去均值。
接着就可以按照eigenvalue对eigenvector进行排序,根据需要cover多少variance来选几个eigenvector了。为什么这么选,可以想一下coviance matrix对角线上的值的意义。
import numpy as np
def pca(data,normalise=1):
# centre data
m = np.mean(data, axis=0)
data -= m
# Covariance matrix
C = np.cov(np.transpose(data))
# Compute eigenvalues and sort into descending order
evals, evecs = np.linalg.eig(C)
indices = np.argsort(evals)
indices = indices[::-1]
evecs = evecs[:, indices]
evals = evals[indices]
evalsNorm = evals / np.sum(evals)
cumEvals = np.cumsum(evalsNorm)
evecs = evecs[:,:(np.where(cumEvals > 0.9)[0][0] + 1)]
if normalise:
for i in range(np.shape(evecs)[1]):
evecs[:, i] = evecs[:, i] / np.linalg.norm(evecs[:, i])
# Produce the new data matrix
# x = np.transpose(np.dot(np.transpose(evecs), np.transpose(data)))
x = np.dot(data, evecs)
# Compute the original data again
y = np.transpose(np.dot(evecs, np.transpose(x))) + m
return x,y,evals,evecs
iris = np.genfromtxt(r"/home/real/pca/iris.data", delimiter=",", dtype="str")
iris = iris[:, range(np.shape(iris[0])[0] - 1)]
iris = iris.astype(np.float)
x,y,evals,evecs = pca(iris)
np.savetxt("result.csv", x, delimiter=",")
np.savetxt("resulty.csv", y, delimiter=",")