机器学习系列(十九)——PCA进行特征降维
求数据的前n个主成分
求出第一主成分后,将数据在第一个主成分上的分量去掉,得到新的数据,再在新数据上求第一主成分即可得到原数据的第二主成分。这个过程继续进行可以得到高维数据的第三、第四...主成分。n个维度最多有n个主成分。
以二维特征数据为例,如图所示 :
其第一主成分是w=(w1,w2)对应的方向,要求第二主成分,则在原数据X(i)上减去第一主成分对应分量X(pro),在得到的新数据上进行第一主成分求解。
接下来生成模拟数据集:
'''模拟数据集'''
import numpy as np
import matplotlib.pyplot as plt
X = np.empty(shape=(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)
def demean(X):
'''axis=0,行方向求得每一列均值'''
return X - np.mean(X,axis=0)
X = demean(X)
'''demean之后的X绘图'''
plt.scatter(X[:,0],X[:,1])
plt.show()
demean后数据图示
该数据集两个特征之间的关系大致符合y=0.75x+3,加入了微小噪音。在此数据集上求解主成分函数:
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):
'''将w模变为即单位向量'''
return w/np.linalg.norm(w)
def first_component(X, initial_w, eta, 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)#注意1:每次都要单位化比较好
if(abs(f(w,X)-f(last_w,X)) < epsilon):
break
cur_iter += 1
return w
基本上是复用了上次的代码,不过这里没有df_debug,因为已经测试过知道梯度表达式是正确的了。下面求解第一主成分:
initial_w = np.random.random(X.shape[1])
eta = 0.01
'''第一主成分'''
w = first_component(X,initial_w,eta)
w
求解结果:
第一主成分求解第二主成分之前先给出减去第一主成分分量后新数据的样子:
X2 = np.empty(X.shape)
for i in range(len(X)):
'''计算新数据集'''
X2[i] = X[i] - X[i].dot(w) * w
##向量化 X2 = X - X.dot(w).reshape(-1,1)*w
'''新数据X2的样子'''
plt.scatter(X2[:,0],X2[:,1])
plt.show()
新数据图示
可见符合直线分布,这是因为我们生成的数据性质过于良好,而且可以发现,新数据符合的直线分布的方向向量(实际是第二主成分)和第一主成分方向向量垂直。
下求解第二主成分:
w2 = first_component(X2,initial_w,eta)
w2
求解结果:
第二主成分可以验证,和第一主成分是垂直的。
封装求解前n个主成分的函数:
'''封装求前n主成分函数'''
def first_n_component(n, X, eta=0.01, n_iters = 1e4, epslion=1e-8):
X_pca = X.copy()
X_pca = demean(X_pca)
res = []
for i in range(n):
initial_w = np.random.random(X_pca.shape[1])
w = first_component(X_pca,initial_w,eta)
res.append(w)
X_pca = X_pca - X_pca.dot(w).reshape(-1,1)*w
return res
使用我们封装的函数:
第一第二主成分可见求解和上面是一致的。
高维数据向低维映射——PCA实现降维
PCA的本质就是将数据从一个坐标系转换到另外一个坐标系,即找到一组相互正交的投影方向,使得数据在这些投影方向上的方差最大(找新的一组正交基)。计算原始数据在这些正交基上投影的方差,方差越大,就说明在对应正交基上包含了更多的信息量。
n维特征有n个主成分,往往我们取前k(k<=n)个主成分,因为它们更加重要,这样便达到了降维的目的。
的第k行是第k主成分,于是对原数据X的降维表示为:
接下来模拟sklearn逻辑,封装PCA算法,在play_Ml中添加PCA.py:
import numpy as np
class PCA:
def __init__(self,n_components):
'''初始化PCA'''
self.n_components = n_components
self.components_ = None
def fit(self, X, eta=0.01, n_iters=1e4):
'''获得数据集的前n_component个主成分'''
assert self.n_components <= X.shape[1],"主成分数不能大于总特征数"
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):
'''将w模变为即单位向量'''
return w/np.linalg.norm(w)
def first_component(X, initial_w, eta, 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)#注意1:每次都要单位化比较好
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],"must be the same!"
return X.dot(self.components_.T)
def inverse_transform(self,X):
'''将降维后的数据升维到n维'''
assert X.shape[1]==self.components_.shape[0],"must be the same!否则不能做乘法"
return X.dot(self.components_)
def __repr__(self):
return "PCA(n_components=%d)" % self.n_components
在模拟数据上使用我们的PCA算法:
'''模拟数据集'''
import numpy as np
import matplotlib.pyplot as plt
X = np.empty(shape=(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 play_Ml.PCA import PCA
pca = PCA(n_components=2)
pca.fit(X)
2个主成分求解结果
这里使用的二维特征数据,对二维特征降维只需要求解一个主成分:
降维结果同样,由矩阵的运算过程降维后的数据还可以被升维回去,只是损失了部分数据,升维后的数据不会和原始数据一样:
升维结果下面将升维后的数据和原始数据做图对比:
'''绘制原始数据和降维后升维得到的数据'''
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()
升维数据与原始对比
恢复回去的点会对应主成分的位置,和原数据相比损失了一些信息,因此只是在高维空间表达低维样本而已。