大数据,机器学习,人工智能机器学习、深度学习与人工智能人工智能/模式识别/机器学习精华专题

机器学习系列(十九)——PCA进行特征降维

2019-07-13  本文已影响1人  Ice_spring

求数据的前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)个主成分,因为它们更加重要,这样便达到了降维的目的。
X= \begin{pmatrix} X_{1}^{(1)} , X_{2}^{(1)},...,X_{n}^{(1)}\\X_{1}^{(2)} , X_{2}^{(2)},...,X_{n}^{(2)}\\...\\X_{1}^{(m)} , X_{2}^{(m)},...,X_{n}^{(m)} \end{pmatrix} , W_{k}= \begin{pmatrix} W_{1}^{(1)} , W_{2}^{(1)},...,W_{n}^{(1)}\\W_{1}^{(2)} , W_{2}^{(2)},...,W_{n}^{(2)}\\...\\W_{1}^{(k)} , W_{2}^{(k)},...,W_{n}^{(k)} \end{pmatrix}
W_{k}的第k行是第k主成分,于是对原数据X的降维表示为:
X_{k}=X.W_{k}^{T}=\begin{pmatrix} X_{1}^{(1)} , X_{2}^{(1)},...,X_{k}^{(1)}\\X_{1}^{(2)} , X_{2}^{(2)},...,X_{k}^{(2)}\\...\\X_{1}^{(m)} , X_{2}^{(m)},...,X_{k}^{(m)} \end{pmatrix}
接下来模拟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()
升维数据与原始对比

恢复回去的点会对应主成分的位置,和原数据相比损失了一些信息,因此只是在高维空间表达低维样本而已。

上一篇下一篇

猜你喜欢

热点阅读