此文是自己的学习博客,主要参考Rachel-Zhang的Stanford机器学习—第十讲. 数据降维
博客中说到,感谢先行的数学家,让许多问题在数学上根本不是问题,也要感谢你们将自己的理解写成博客分享给后来学习者~thx
一、为什么要降维?
我们需要一组关于XXX的数据,定义就铺天盖地的来了,百万级个特征拿过来,我们怎么进行机器学习啊?!李航老师在他的博客《机器学习新动向:从人机交互中学习》中提到,学习精度越高,学习确信度越高,学习模型越复杂,所需要的样本也就越多。样本复杂度满足以下不等式
由此可见,feature太多会造成模型复杂,训练速度过慢,因此我们引入降维。
About Visualization:
多维数据很难进行可视化分析,因此我们需要降维分析。
二、Principal components analysis 主成分分析PCA
这里有一篇博客PCA的数学原理(非常值得阅读)!!!!讲解的非常非常清楚。
总结下就是:
1. 数据的特征向量表示
2. 基变换(我们降维目的就是是把原来的数据转换到基数比较少的基坐标上)
3. 基变换的矩阵表示(因为数据都是多维的,多个基构成矩阵A,多维向量构成矩阵B)
4. 如何选择基才是最优的,引入协方差矩阵
5. 很巧,协方差矩阵的对角线为数据投影到一个基上的方差,越大越好;其他元素表示不同基之间的协方差,相关性越小越好
6. 协方差矩阵
7. 寻找这样的基坐标构成的矩阵使得协方差矩阵对角化!(这个在线性代数里面考烂了,对称矩阵对角化)
有个疑问,怎么降到多少维呢??这个怎么确定,也就是矩阵A的维度,基坐标的个数怎么确定?
还有个疑问,基坐标构成的矩阵A应该是方阵吧,比如,那么矩阵A的列数等于矩阵B的个数,而矩阵B的行数就是数据的特征个数吧?这有点矛盾啊??
貌似理解了,之前基坐标个数确实与特征个数一致,然后选取其中的部分向量就满足降维的目的,但这个部分怎么取呢?博客中从二到一很简单,1000呢,降维到100还是200呢?
三、PCA算法流程
3.1 数据预处理
假设有m个样本,n维特征,即
第1、2步,去均值化
第3、4步,归一化:ensures that different attributes are all treated on the same“scale.”
3.2 PCA算法选取k个主分量
1.求协方差矩阵
2.奇异值SVD求取特征值和特征向量:
from numpy import *
A = array([[1,1],[0,0],[1,1]])
# print(Base)
U,S,V = linalg.svd(A)
print(U,S,V)
[[ -7.07106781e-01 -7.07106781e-01 -1.57009246e-16]
[ 0.00000000e+00 0.00000000e+00 1.00000000e+00]
[ -7.07106781e-01 7.07106781e-01 1.57009246e-16]]
[ 2.00000000e+00 3.35470445e-17]
[[-0.70710678 -0.70710678]
[ 0.70710678 -0.70710678]]
SVD为奇异值分解(singular value decomposition),在numpy库中有函数[U,S,V] = linalg.svd(A) 返回一个数组S(由Σ的奇异值组成),两个酉矩阵U和V,且满足A= U*S* V。若A为m×n阵,则U为m×m阵,V为n×n阵。奇异值在S的对角线上,非负且按降序排列。(奇异值等于sqrt(特征值))
协方差矩阵肯定是方阵,那么对于方阵Σn*n呢,就有
Σ = USV
ΣΣ’= USV*VS’U = U(ΣΣ’)U
ΣΣ = VS’U*USV = V(Σ’Σ)V
i.e. U是ΣΣ’的特征向量矩阵;V是Σ’Σ的特征向量矩阵,都是n*n的矩阵
由于方阵的SVD相当于特征值分解,所以事实上U = V, 即Σ = USU’, U是特征向量组成的正交矩阵
我们的目的是,从n维降维到k维,也就是选出这n个特征中最重要的k个,也就是选出特征值最大的k个~
3.选择k个分量
svd分解后,我们得到了一个n×n的协方差矩阵Σ和U,这时,我们就需要从U中选出k个最重要的分量;即选择前k个特征向量,即为Ureduce,矩阵大小为n*k
这样对于一个n维向量x,就可以降维到k维向量z了:
注意这里有个转置,是因为svd求出来的U矩阵,其每一列对对应特征值的特征向量,而上图中将U转置,行向量是对应的特征向量。
3.3 从压缩数据中恢复原数据
我们已经知道,可以根据z(i) = Ureduce’× x(i) 将n维向量x降维到k维向量z,那么有时我们需要恢复n维向量,怎么做呢?
由于Ureduce是正交矩阵(下面Ureduce简记为U),即U’ = U-1, 所以
xapprox = (U’)-1×z = (U-1)-1×z = Uz
(PS:这里的逆操作为伪逆操作)
注意:这里恢复出的xapprox并不是原先的x,而是向量x的近似值。
3.4 怎样决定降维个数/主成分个数
这个就是之前思考的问题呀~
首先从一个general一点的思路去想呢,我们是希望,选出主成分之后进行数据分析,不会造成大量特征的丢失,也就是说可以用下式的error ratio表示经过压缩后的性能如何~果然是需要一个度量的~~
从数学上可以证明,上面这个error ratio 可以表示为
所以,可以用下式进行k的合理选取:
3.5 应用PCA进行降维的建议
PCA可以降维,那让我们联想到了之前说过拟合问题是由维度过高或者参数过多造成的,那么可不可以用PCA解决overfitting的问题呢?
Ans:NO!应用PCA提取主成分可能会解决一些overfitting的问题,但是呢,不建议用这种方法解决overfitting问题,还是建议用第三章中讲过的加入regularization项(也称为ridge regression)来解决。PCA中主成分分析应用到那部分数据呢?
Ans:Only Training Data!可以用Cross-Validation data 和 test Data进行检验,但是选择主分量的时候只应用training data.不要盲目PCA
Notice:only 当你在原数据上跑到了一个比较好的结果,又嫌它太慢的时候才采取PCA进行降维,不然降了半天白降了~
四、python代码实现
import numpy as np
from numpy.matrixlib.defmatrix import mat
##导入文件
def loadData(filename):
with open(filename, 'r')as fr: ##fr = filename.read()
for line in fr.readlines():
lineArr = line.strip().split('\t')
dataMat.append([float(lineArr[0]), float(lineArr[1]), float(lineArr[2])])##如果维度很大呢,这时候怎么写?下面的尝试都不成功
# for i in range(len(lineArr)-1):
# dataMat.append([float(lineArr[i])]) ###这里用循环感觉比较复杂,不知道有没有简单的写法
# dataMat.append(float(lineArr[i]) for i in range((len(lineArr)-1)))
return np.array(dataMat)
#数据预处理,去均值化和归一化
def dataProcess(dataMat):
dataSize = dataMat.shape[0]
meanData = np.mean(dataMat,axis=0)
newData = dataMat - meanData
dataVariance = ((newData*newData).sum(0))/dataSize
newData = newData/(dataVariance)**0.5
return newData,dataVariance,meanData
##求协方差并用svd求特征值和特征向量
def svd(newData):
covData = np.cov(newData.T,rowvar=True) ##rowvar=True 表示一行代表一个向量,在我们的文本中,一行表示一个样本
U,S,V = np.linalg.svd(covData)
##因为covData是协方差矩阵,则U=V,且U的每一列代表对应特征值的特征向量. S是降序排列的奇异值
return U,S
##根据百分比确定k的值,选择降维到多少
def percentage2n(S,percentage):
eig_sum = sum(S**2) ##特征值的总和
k = 1
eig_k = S[0:k]
eig_k = sum(eig_k**2)
i = eig_k/eig_sum
if i <= percentage:
k += 1.0
else:
return k
##选取k个分量
def lowData(U,k,newData,dataVariance,meanData):
new_U = U[:,0:k] ##列向量表示特征向量,故应该取前k个列向量
lowData = mat(newData) * mat(new_U) ##[m*n]*[n*k]
##重构数据,用于调试 不过这里貌似并不需要。我们已经用了特征值之比来计算error ratio 了
# reconData = (lowData * new_U.T) * (dataVariance**0.5) + meanData
return lowData
if __name__ == '__main__':
dataMat = []
dataMat = loadData('iris_data.txt')
newData,dataVariance,meanData = dataProcess(dataMat)
U,S = svd(newData)
k = percentage2n(S,0.80)
lowData = lowData(U,k,newData,dataVariance,meanData)
print(lowData)