一、PCA
1.PCA思想
PCA(principle component analysis),即主成分分析法,是一个非监督的机器学习算法,主要用于对数据的降维,通过降维可以发现更便于人理解的特征,加快对样本有价值信息的处理速度,此外还可以应用于可视化(降到二维)和去噪。
例如下图所示,样本有2个特征,现在对该样本进行降维处理。首先考虑直接选择特征1或者特征2降维,经过降维后的样本由2维降到1维,如图所示。
选择特征1降维。
选择特征2降维。
可以发现,相对于选择特征2,选择特征1降维后,样本间的距离更大,即样本的可区分度更高。那是否还有其他的映射方式,使得映射后样本的间距更大,事实上还可以选择某个轴线,例如下图所示,样本映射到该轴线之后,有更大的间距。
PCA降维的思想就是寻找某个轴线,使得样本映射到该轴线上后,能够有最大的可区分度。衡量可区分度的指标即方差(事实上,方差背后的意义就是数据的稀疏程度),现在的问题就是如何求得该轴线,使得方差的值最大。
2.求解思路
用方差来定义样本的间距,方差越大表示样本分布越稀疏,方差越小表示样本分布越密集,方差的公式如下。
在求解最大方差前,为了方便计算,可以先对样本进行demean(去均值)处理,即减去每个特征的均值,这种处理方式不会改变样本的相对分布(效果就像坐标轴进行了移动)。去均值后,样本x每个特征维度上的均值都是0,方差的公式转换成图示的公式。在这里,xi代表已经经过映射后的某样本。
对于只有2个维度的样本,现在的目标就是:求一个轴的方向w=(w1,w2),使得映射到w方向后,方差最大。目标函数表示如下。表示经过映射后的样本。
如下图所示,原始样本和映射后的样本有如下的空间几何关系。w为单位向量,所以w的模为1。
经过什么的推导,最终求解的公式可以转换成如下形式。
样本x为已知,现在的目标转换成如下。
这就转换成了求目标函数的最优化问题,可以用梯度上升法进行求解。
二、梯度上升法求解
1.公式求解
最终公式可以化简为如下形式。
梯度上升法的思路和梯度下降法的思路完全一样,只是前者求极大值,后者求极小值。可以参考梯度下降法博客内容。
2.代码实现
import numpy as np
import matplotlib.pyplot as plt
x = np.empty((100,2)) #x为100个样本,2个特征值的矩阵
x[:,0]= np.random.uniform(0.,100.,size = 100) #样本x的第1个特征值为均值为0,方差为100的100个值
x[:,1]= 0.75*x[:,0]+3.+np.random.normal(0,10.,size=100) #样本x的第2个特征值和第1个特征值保持一定的线性关系,这里引入噪声。这里之所以使特征1和特征2有一定的线性关系,仅仅是为了降维后,能够有较好的展示。
plt.scatter(x[:,0],x[:,1]) #将特征1和特征2绘制出来
plt.show()
初始样本准备完毕,现在就开始对该样本进行降维处理。
def demean(x):
return x-np.mean(x,axis=0) #首先对每个特征值进行去均值,axis=0代表按列求均值
def f (w,x):
return np.sum((x.dot(w)**2))/len(x) #即方差公式,len(x)即样本的个数
def df_math(w,x):
return x.T.dot(x.dot(w))*2./len(x) #即通过数学推导得出的求梯度公式
def direstion(w):
return w /np.linalg.norm(w) #将w进行归一化,因为w仅代表方向,而每次求解后,w模不一定为1
def gradietn_ascent(df,x,initial_w,eta,n_iters = 1e4,epsilon=1e-8): #梯度上升法,传入初始w,迭代次数和停止迭代的最小差值
w = direstion(initial_w) #先将传入的w进行归一化
cur_iter =0
while cur_iter <n_iters: #这里的思想和梯度下降法完全一致,只是此时,eta为正,即朝着函数值增大方向
gradient = df(w,x)
last_w = w
w = w+eta*gradient
w = direstion(w)
if(abs(f(w,x)-f(last_w,x))<epsilon): #如果目标函数的值差距小于某个值,则表示差不多已经到达极值点,此时停止迭代
break
return w
编写完核心代码后,现在就开始对样本x进行降维。
x_mean = demean(x) #首先去均值处理
initial_w = np.random.random(x.shape[1]) #初始化不能为0向量,因为0方向的向量并不代表某个具体方向,相当于没初始化w,此外w的维度和x特征维度保持一致,由上面推导的公式可以得出该结论。
eta = 0.01 #初始化步长
w=gradietn_ascent(df_math,x_mean,initial_w,eta) #降维并得出第一个方向,即第一个主成分,使得在该方向上进行映射可以有最大的间距(方差)
plt.scatter(x_mean[:,0],x_mean[:,1]) #首先将原始数据绘制
plt.plot([-w[0]*60,w[0]*60],[-w[1]*60,w[1]*60],color = 'r') #将w表示的方向进行绘制,为了更清楚地展示,这里适当扩大w的表示范围,但是w0和w1比例要保持一致。
plt.show()
这里求得的仅仅是1个主成分,在更高维度的空间,可能还需要继续求出其他的主成分,即降更多的维度。
三、求前n个主成分
1.求解思路
在上面,对于只有2维特征的样本x,我们求解出了第一个主成分,即w=(w1,w2),为了求解第二个主成分,这里需要将样本x在第一个主成分上的分量去除掉,这里使用的方法即空间几何的向量减法,得到的结果即下图中的绿线部分。
在得到绿线部分的分量后,再对该分量重新求第一主成分,以此类推,可以得到前n个主成分。
2.代码实现
首先给出求第一主成分的核心代码,这个和上面的核心代码完全一样,只是函数名字稍微有点变化。
def demean(x):
return x-np.mean(x,axis=0)#对每个特征值进行求均值
def f (w,x):
return np.sum((x.dot(w)**2))/len(x)
def df(w,x):
return x.T.dot(x.dot(w))*2./len(x)
def direstion(w):
return w /np.linalg.norm(w) #将w进行归一化,因为w仅代表方向
def first_conponent(df,x,initial_w,eta,n_iters = 1e4,epsilon=1e-8):
w = direstion(initial_w)
cur_iter =0
while cur_iter <n_iters:
gradient = df(w,x)
last_w = w
w = w+eta*gradient
w = direstion(w)
if(abs(f(w,x)-f(last_w,x))<epsilon):
break
return w
接下来求解第2个主成分,首先求去除掉第一个主成分后样本的值。
x_mean = demean(x)
initial_w = np.random.random(x.shape[1])
eta = 0.01
w = first_conponent(df,x_mean,initial_w,eta) #这里求解出第一个主成分
x2 = np.empty(x.shape) #初始1个矩阵,和x的维度保持一致,用来存储去除掉第一个主成分分量后的样本
for i in range(len(x)):
x2[i] =x[i]-x[i].dot(w)*w #得到每个样本在去除掉第一个主成分分量后的值
plt.scatter(x2[:,0],x2[:,1]) #将该分量绘制出来
plt.show()
接下来,对剩下的分量求第一主成分。
w2 = first_conponet(df,x2,initial_w,eta)
w2.dot(w) #验证,结果相乘接近0,因为第一主成分和第二主成分垂直,所以矩阵相乘为0
![](https://i-blog.csdnimg.cn/blog_migrate/0efb89d51c5549bd6d9c5579fd9ee2c7.png)
3.封装代码
将上面的代码全部进行封装。
def demean(x):
return x-np.mean(x,axis=0)
def f (w,x):
return np.sum((x.dot(w)**2))/len(x)
def df(w,x):
return x.T.dot(x.dot(w))*2./len(x)
def direstion(w):
return w /np.linalg.norm(w)
def first_conponent(x,initial_w,eta,n_iters = 1e4,epsilon=1e-8):
w = direstion(initial_w)
cur_iter =0
while cur_iter <n_iters:
gradient = df(w,x)
last_w = w
w = w+eta*gradient
w = direstion(w)
if(abs(f(w,x)-f(last_w,x))<epsilon):
break
return w
def first_n_componets(n,x,eta=0.01,n_iters =1e4,epsilon=1e-8):#这里将求前n个主成分进行封装
x_pca = x.copy() #首先将原始样本拷贝给x_pca
x_pca = demean(x_pca)
res =[] #用来存n个主成分
for i in range(n):
initial_w = np.random.random(x_pca.shape[1])
w = first_conponent(x_pca,initial_w,eta)
res.append(w)
x_pca = x_pca-x_pca.dot(w).reshape(-1,1)*w #之前用的是for循环,实际上可以进行向量化,向量化后可以优化计算速度
return res
使用代码,这个时候可以传入参数,指明需要求几个主成分,即保留多少个特征。
first_n_componets(2,x)
上面仅仅使用2个特征值的样本,接下来使用真实的数据集进行降维,并求原始样本向n个主成分降维的结果。
四、高维数据向低维数据转化
1.求解思路
现在有原始样本x,并求得了前k个主成分,形成wk矩阵,如下图所示。现在需要将x映射到新的轴上,形成只有k个特征的新样本xk,用该样本来描述原来的样本x,求解公式如下。
需要注意的是,映射后数据已经丢失了一部分信息,所以做逆运算时,并不能回到最初x的数据。
2.代码实现
在这里创建PCA.py文件(放在ml文件夹下),写入如下代码。
import numpy as np
class PCA:
def __init__(self, n_components):
"""初始化PCA"""
assert n_components >= 1, "n_components must be valid"
self.n_components = n_components #传入需要保留多少
self.components_ = None
def fit(self, X, eta=0.01, n_iters=1e4):
"""获得数据集X的前n个主成分"""
assert self.n_components <= X.shape[1], \
"n_components must not be greater than the feature number of X"
def demean(X):
return X - np.mean(X, axis=0)
def f(w, X):
return np.sum((X.dot(w) ** 2)) / len(X)
def df(w, X):
return X.T.dot(X.dot(w)) * 2. / len(X)
def direction(w):
return w / np.linalg.norm(w)
def first_component(X, initial_w, eta=0.01, n_iters=1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w)
if (abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
X_pca = demean(X)
self.components_ = np.empty(shape=(self.n_components, X.shape[1]))
for i in range(self.n_components):
initial_w = np.random.random(X_pca.shape[1])
w = first_component(X_pca, initial_w, eta, n_iters)
self.components_[i,:] = w
X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
return self
def transform(self, X):
"""将给定的X,映射到各个主成分分量中"""
assert X.shape[1] == self.components_.shape[1]
return X.dot(self.components_.T)
def inverse_transform(self, X):
"""将给定的X,反向映射回原来的特征空间"""
assert X.shape[1] == self.components_.shape[0]
return X.dot(self.components_)
def __repr__(self):
return "PCA(n_components=%d)" % self.n_components
使用代码:
import numpy as np
import matplotlib.pyplot as plt
x = np.empty((100,2))
x[:,0]= np.random.uniform(0.,100.,size = 100)
x[:,1]= 0.75*x[:,0]+3.+np.random.normal(0,10.,size=100) #这个数据集和上面使用的是一样的数据
from ml.PCA import PCA
pca = PCA(n_components=1)
pca.fit(x)
x_reduction = pca.transform(x)
x_restore = pca.inverse_transform(x_reduction) #即将一维的数据重新映射回二维,仅仅是展示用
plt.scatter(x[:,0],x[:,1],color='b',alpha=0.5)
plt.scatter(x_restore[:,0],x_restore[:,1],color='r',alpha=0.5)
plt.show() #可以发现重新映射回原来矩阵形式时,只剩下主成分方向上的数据,restore的功能只是在高维的空间表达低维的数据
![](https://i-blog.csdnimg.cn/blog_migrate/9302756e2f87ef2dff8200f4758f3b19.png)
3.sklearn中的PCA
from sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(x)
pca.components_
可以发现求得的方向和之前的自己编写的相反。需要注意的是,sklearn降维并不是用梯度上升的计算方法,且其内部有一定的优化。自己编写的代码只是为了更好地理解PCA原理。
4.对手写识别数据进行降维
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
digits = datasets.load_digits()
x = digits.data
y = digits.target
from sklearn.model_selection import train_test_split #这里依然需要区分测试数据集和测试数据集,不理解的可以参考knn博客
x_train,x_test,y_train,y_test = train_test_split(x,y,random_state = 666)
from sklearn.neighbors import KNeighborsClassifier #因此目前为止,仅仅使用过knn算法进行分类,所以使用knn进行预测。
knn = KNeighborsClassifier()
knn.fit(x_train,y_train)
knn.score(x_test,y_test)
接下来使用降维的方式,再次用knn算法进行分类,并求分值。
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(x_train)
x_train_reduction = pca.transform(x_train)
x_test_reduction = pca.transform(x_test)#测试数据集也必须降维
knn = KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
knn.score(x_test_reduction,y_test)
可以发现,经过降维后,因为维度就剩下2个,得出的分值很低,为了提高分值,应该保留更多的特征值。但是具体保留多少特征合适?有什么衡量指标可以作为参考?
sklearn-PCA中有个变量为pca.explained_variance_ratio_,可以用来衡量映射后的数据保留了原始数据的多少信息(方差)。对于上面求得的结果,有如下的值。
这表示:经过映射后的数据,只保留原始数据大约百分之(14+13=27)的方差,丢掉了原来约73的方差,这个损失是比较大的,所以得出来的分数很低。
现在的问题是,应该保留多少的维度使得信息丢失在允许法范围内。假设选取的主成分的个数和原始样本特征维度一样,由于explained_variance_ratio_是按照大小顺序排列,这时候观察explained_variance_ratio_发现,越到后面,值越小,甚至几乎可以忽略。
pca = PCA(n_components=x_train.shape[1]) #保留所有的维度
pca.fit(x_train)
pca.explained_variance_ratio_
现在将这些值从前往后依次进行累加,并在坐标轴上展现。
plt.plot([i for i in range(x_train.shape[1])],[np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(x_train.shape[1])])
plt.show()
观察该图可以发现,当横轴为0时,纵轴也为0(丢失了所有信息),当横轴为特征值大小时,纵轴为1(保留了全部信息)。因此,比如想要保留95%的信息,就可以通过该图知道应该保留多少个特征维度。
pca = PCA(0.95) #传入的参数为0.95,代表保留的百分之95的信息
pca.fit(x_train)
pca.n_components_ #结果为28
x_train_reduction = pca.transform(x_train)
x_test_reduction = pca .transform(x_test)
knn = KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
knn.score(x_test_reduction,y_test)
值明显得到了提高。而且如果进行时间测试,会发现时间缩短了很多。这就是pca的作用。
之前已经说过,pca经过降维可以用于可视化。将高维度难以想象的数据,降为2维可以在平面上展示。手写识别数据降为2维的图像如下。
pca = PCA(n_components=2)
pca.fit(x)
x_reduction = pca.transform(x)
for i in range(10):
plt.scatter(x_reduction[y==i,0],x_reduction[y==i,1],alpha=0.8) #将每种数字进行统一绘制,并用不同颜色展示
plt.show() #可以发现即使是2维,也有较好的区分度
![](https://i-blog.csdnimg.cn/blog_migrate/156ab0251ce2df81b3bd221377cbfa24.png)
5.MNIST数据集降维
import numpy as np
from sklearn.datasets import fetch_mldata
mnist = fetch_mldata("MNIST original") #第一次加载会需要一些时间
x,y = mnist['data'],mnist['target']
x_train = np.array(x[:60000],dtype=float)
y_train = np.array(y[:60000],dtype=float)
x_test = np.array(x[:60000],dtype=float)
y_test = np.array(y[60000:],dtype=float)
不使用pca时,直接经过knn会需要很长的时间。
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier()
%time knn.fit(x_train,y_train)
%time knn.score(x_test,y_test)
使用pca后发现,时间提高很多。奇怪的是,虽然特征减低了,但是分值却提高了,这就涉及到pca的另一个重要的作用-降噪,即在降维的过程中也将噪声减小了。
from sklearn.decomposition import PCA
pca = PCA(0.9)
pca.fit(x_train)
x_train_reduction = pca.transform(x_train)
x_test_reduction = pca.transform(x_test) #测试数据集也必须降维
knme knn.fit(x_train_reduction,y_train)
%time knn.score(x_test_reduction,y_test)
6.降噪
首先看之前的一个数据。
import numpy as np
import matplotlib.pyplot as plt
x = np.empty((100,2))
x[:,0]= np.random.uniform(0.,100.,size = 100)
x[:,1]= 0.75*x[:,0]+3.+np.random.normal(0,5.,size=100)
from sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(x)
x_reduction = pca.transform(x)
x_restore = pca.inverse_transform(x_reduction)
plt.scatter(x[:,0],x[:,1],color='r')
plt.scatter(x_restore[:,0],x_restore[:,1])
plt.show()
原来的数据可能是因为测量的过程中出现了问题,所以发生了一定的抖动(即产生了噪声),经过降维后可以发现最终的数据在一条直线上,也就是将噪声进行了去除,因此说pca有一定的降噪能力。下面使用手写识别例子继续探讨降维中的去燥功能。
from sklearn import datasets
digits = datasets.load_digits()
x = digits.data
y = digits.target
noisy_digits = x + np.random.normal(0,4,size=x.shape) #这里加入随机噪声
example_digits = noisy_digits[y==0,:][:10]
for num in range(1,10):
x_num = noisy_digits[y==num,:][:10]
example_digits = np.vstack([example_digits,x_num])#将全部样本按照0-10的顺序叠起来,共100个样本
def plot_digits(data): #这个函数当做模块使用即可,不必深究
fig,axes = plt.subplots(10,10,figsize=(10,10),subplot_kw={'xticks':[],'yticks':[]},gridspec_kw=dict(hspace=0.1,wspace=0.1))
for i,ax in enumerate(axes.flat):
ax.imshow(data[i].reshape(8,8),cmap='binary',interpolation='nearest',clim=(0,16))
plt.show()
plot_digits(example_digits) #调用函数
可以看到图像十分模糊,现在通过pca进行降噪处理。
![](https://i-blog.csdnimg.cn/blog_migrate/4ea3691745f71d04459eeb447f29d996.png)
pca = PCA(0.5) #因为本身噪声比较大,所以这里仅保留原始数据0.5的信息
pca.fit(noisy_digits)
pca.n_components_ #使用了12个特征
components = pca.transform(example_digits)
filtered_digits = pca.inverse_transform(components)
plot_digits(filtered_digits)
![](https://i-blog.csdnimg.cn/blog_migrate/023fae602e37213f0098806660f640bd.png)
发现图像的清晰度有所提高。