机器学习算法之PCA

最新推荐文章于 2024-06-18 14:28:36 发布

cmvaxx

最新推荐文章于 2024-06-18 14:28:36 发布

阅读量1.6k

点赞数 1

分类专栏：机器学习文章标签：算法人工智能数据分析

本文链接：https://blog.csdn.net/m0_52004821/article/details/127728950

版权

机器学习专栏收录该内容

5 篇文章 0 订阅

订阅专栏

1.PCA简介

PCA（Principal Component Analysis）：也是一个梯度分析的应用，不仅是机器学习的算法，也是统计学的经典算法
1-1

1.1 举个栗子

例如下面一个两个特征的一个训练集，我们可以选择一个特征，扔掉一个特征
1.1-1
下图分别是扔掉了特征一和特征二的两种方案，很明显右边这种的效果会更好一些，因为访问二扔掉特征二以后，点之间的分布情况更接近与原图，但是这不是更好的
1.1-2
我们希望有一根直线，是斜着的，我们希望将所有的点都映射到这条直线上，那么这个时候我们就成功的将二维降到了一维，与此同时，这些点更加趋近与原来的点的分布情况，换句话说，点和点之间的距离比无论是映射到x还是映射到y周，他们之间的区分度都更加的大，也就更加容易区分
1.1-3

1.2 总结

那么如何找到这个让样本间间距最大的轴?
如何定义样本间间距?
事实上有一个指标可以之间定义样本间的距离，就是方差（Variance）（方差：描述样本整体之间的疏密的一个指标，方差越大，代表样本之间越稀疏，方差越小，代表样本之间越紧密）
1.2-1：方差
1-6

第一步: 将样例的均值归为0 (demean)（归0：所有样本都减去他们的均值），使得均值为0，这样可以简化方差的公式

1.2-2

1.3 推导

1.3-1
进行均值归0操作以后，就是下面的式子
1.3-2

注：|X_project|的平均值也是一个向量

X⁽ⁱ⁾映射到w的距离实际上就是X⁽ⁱ⁾与w的点乘（蓝色的线），根据定义推导，其值实际上就是X_project

1.3-3
此时我们的目标函数就可以化简成

1.3-4

这是一个目标函数的最优化问题，使用梯度上升法解决。
当然我们也可以之间使用数学原理推导出结果，这里我们主要关注使用搜索的策略来求解主成分分析法，这样我们对梯度上升发和梯度下降法也可以有一个更深刻的认识

1.4 与线性回归的区别

1.4-1
1.4-2
1.主成分分析法的两个轴都是特征，线性回归y轴是目标结果值
2.主成分分析法的点是垂直于方差轴直线的，线性回归的点事垂直于x轴的

1.5 对PCA和线性代数中的特征向量的对比

一个矩阵可以把一个向量拉伸或者缩短λ倍，这个向量就是特征向量，λ是特征值；
那么能够拉伸最长的特征向量，实际上就是这个矩阵的第一主成分；
一组线性不相关的特征向量可以组成一个特征向量空间；
所以主成分分析实际上就是把矩阵转换到了一个特殊的特征向量空间；
这个特征向量空间是由能够被拉伸最长的前k个特征向量组成的；而且这k个特征向量相互正交；

2.使用梯度上升法解决PCA问题

2-1
1.注意上面式子里的每一个**(X₁⁽ⁱ⁾·w₁+X₂⁽ⁱ⁾·w₂+…X_n⁽ⁱ⁾·w_n)都是一个X⁽ⁱ⁾和w的点乘，所以式子可以进一步化解，
2.化简过后可以进行向量化，即每一个∑(X⁽ⁱ⁾·w₁)·X₁⁽ⁱ⁾** 可以看成是**(X·w)这个向量的转置（本来是个行向量，转置后是1行m列的列向量）与X**这个矩阵（m行n列）做点乘等到的其中一项的相乘相加的结果
3.最后根据转置法则 ((A*B)^T=B^T*A^T)转换成最后的结果
2-2
2-3

3.主成分PCA的实现

使用梯度上升法求解主成分

import numpy as np
import matplotlib.pyplot as plt

构造一个两个样本之间有基本线性关系的数据集，可以使得我们降维的效果更加明显

X = np.empty((100,2))
X[:,0] = np.random.uniform(0.,100.,size=100)
# 0.75倍的X[:,0]加上3加上一个噪音
X[:,1] = 0.75*X[:,0]+3.+np.random.normal(0.,10.,size=100)

plt.scatter(X[:,0],X[:,1])

<matplotlib.collections.PathCollection at 0x111f2e4e0>

3-1

均值归0

def demean(X):
    return X-np.mean(X,axis=0)

X_demean = deamean(X)
plt.scatter(X_demean[:,0],X_demean[:,1])

<matplotlib.collections.PathCollection at 0x111d8e438>

3-2

梯度上升法

def f(w,X):
    return np.sum((X.dot(w)**2))/len(X)

def df_math(w,X):
    return X.T.dot(X.dot(w)) * 2. /len(X)

def df_debug(w,X,epsilon=0.0001):
    res = np.empty(len(w))
    for i in range(len(w)):
        w_1 = w.copy()
        w_1[i] += epsilon
        w_2 = w.copy()
        w_2[i] -= epsilon
        res[i] = (f(w_1,X)-f(w_2,X)) / (2*epsilon)
    return res

def direction(w):
    """计算单位向量"""
    return w / np.linalg.norm(w)

def gradient_ascent(df,X,inital_w,eta,n_iters = 1e4,epsilon=1e-8):
    
    w = direction(inital_w)
    cur_iter = 0
    
    while cur_iter < n_iters:
        gradient = df(w,X)
        last_w = w
        w = w + eta * gradient
        # 注意1：每次求单位向量
        w = direction(w) 
        if abs(f(w,X)-f(last_w,X)) < epsilon:
            break
            
        cur_iter = cur_iter+1
    return w

# 初始值不能为0，因为将0带入求导公式，会发现得0，没有任何方向
# 因为对于我们的目标函数来说，w=0本身就是一个最小值点

## 注意2：不能从0向量开始
inital_w = np.random.random(X.shape[1])
eta = 0.01

# 注意3：不能使用StandardScaler标准化数据
# 因为我们本来就是要使得方差最大，而标准化的目的是使得方差为1

# 使用debug模式
gradient_ascent(df_debug,X_demean,inital_w,eta)
# 输出  array([0.75934073, 0.65069321])

# 使用math数学解
gradient_ascent(df_math,X_demean,inital_w,eta)
# 输出     array([0.75934073, 0.65069321])

w = gradient_ascent(df_math,X_demean,inital_w,eta)
plt.scatter(X_demean[:,0],X_demean[:,1])
# 这个轴就是我们求出的第一个主成分
plt.plot([0,w[0]*30],[0,w[1]*30],color='r')

3-3

4.求数据的前N个主成分

求出第一主成分以后，如何求出下一个主成分?
1.数据进行改变，将数据在第一个主成分上的分量去掉

X⁽ⁱ⁾·w = ||X_project⁽ⁱ⁾|| 即X⁽ⁱ⁾映射到w上的值，那么||X_project⁽ⁱ⁾||（大小） ·w（方向）就是X⁽ⁱ⁾在w上的分向量记为X_project⁽ⁱ⁾= ||X_project⁽ⁱ⁾|| ·w
X⁽ⁱ⁾-X_project⁽ⁱ⁾就可以实现将X样本在X_project相应上的分量去掉，相减之后的集合意义就是讲X样本映射到了X_project向量相垂直的一个轴上，记为X^`(i) = X_project⁽ⁱ⁾
4-1

2.在新的数据上求第一主成分
得到的X^` 是X中的所有样本都去除了第一主成分上的分量得到的结果，要求第二主成分，只要在新的数据上，重新求一下第一主成分

4.1 获得前n个主成分实现

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_componet( X,inital_w,eta,n_iters = 1e4,epsilon=1e-8):
    w = direction(inital_w)
    cur_iter = 0
    
    while cur_iter < n_iters:
        gradient = df(w,X)
        last_w = w
        w = w + eta * gradient
        # 注意1：每次求单位向量
        w = direction(w) 
        if abs(f(w,X)-f(last_w,X)) < epsilon:
            break
            
        cur_iter = cur_iter+1
    return w

inital_w = np.random.random(X.shape[1])
eta = 0.01

first_componet(X_demean,inital_w,eta)

array([0.75934077, 0.65069317])

# X2 = np.empty(X.shape)
# for i in range(len(X)):
#    X2[i] = X[i] - X[i].dot(w)*w
# 向量化，X.dot(w)为m*1的向量，reshape后变成了1*m的列向量，再乘以w（方向）就是X的每一个值在w上 的分量矩阵
X2 = X - X.dot(w).reshape(-1,1)*w

# 相减得到的样本分布几乎垂直于原来的样本分布
plt.scatter(X2[:,0],X2[:,1])

4.1-1

# 求出第二主成分
w2 = first_componet(X2,inital_w,eta)
w2
# 输出     array([-0.65068927,  0.75934411])

# 因为w和w2都是单位向量，所以他们两个点乘得到的结果就是他们夹角的cos值，
# 又因为w和w2应该是互相垂直的，所以他们夹角的cos值等于0
w.dot(w2)
# 输出5.1820671385094386e-06 几乎为0

def first_n_componet(n,X,eta = 0.01,n_iters = 1e4,epsilon = 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_componet(X_pca,initial_w,eta)
        res.append(w)
        X_pca = X_pca - X_pca.dot(w).reshape(-1,1) *w
    return res

first_n_componet(2,X)
# 输出  [array([0.75934077, 0.65069316]), array([-0.65068972,  0.75934372])]

5. 高维数据向低维数据进行映射

对于一个数据集X来说，这个X有m行n列，代表有m个样本n个特征，通过我们前面学习的主成分分析法，假设我们已经求出了针对这个数据来说的前k个主成分，每一个主成分对应一个单位方向，用W来表示,W也是一个矩阵，他有k行，代表我们求出的前K个主成分，每一行有n列，代表每一个主成分的坐标轴应该是有n个元素的。这是因为我们的主成分分析法主要就是将数据从一个坐标系转化成了另外一个坐标系，原来这个坐标系有n个维度，现在这个坐标系也应该有n个维度，只不过对于转化的坐标系来说，我们取出来前k个，这k个方向更加重要。

如何将我们的样本X从n维转化成k维呢，回忆们之前学到的，对于一个X样本，与一个W进行点乘，其实就是讲一个样本映射到了w这个坐标轴，得到的模，如果讲这一个样本和这k个w分别做点乘，得到的就是这一个样本，在这k个方向上做映射后每一个方向上的大小，这k个元素合在一起，就代表这一个样本映射到新的k个轴所代表的坐标系上相应的这个样本的大小

X1分别乘以W1到Wn，得到的k个数组成的向量，就是样本1映射到Wk这个坐标系上得到的k维的向量，由于k<n，所以我们就完成了一个样本从n维到k维的映射，这个过程依次类推从样本1到样本m都这么做，我们就将m个样本都从N维映射到了k维-----其实我们就是做了一个乘法X·W^T(为什么是转置呢，因为我们是拿X的每一行去和W的每一行做点乘的，但是矩阵乘法规定是拿X的每一行和W的每一列做乘法)

5-1

我们得到新的降维后的矩阵X_k以后，是可以通过和W_k想乘回复回来的，但是由于我们在降维的过程中丢失了一部分信息，这时及时回复回来也和原来的矩阵不一样了，但是这个从数据角度成立的

5-2

从高维数据向地维数据的映射

class PCA:
    
    def __init__(self,n_components):
        """初始化PCA"""
        assert n_components>=1, "n_components must be vaild"
        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 be greater then the feature number of X"
        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_componet( X,inital_w,eta,n_iters = 1e4,epsilon=1e-8):
            w = direction(inital_w)
            cur_iter = 0

            while cur_iter < n_iters:
                gradient = df(w,X)
                last_w = w
                w = w + eta * gradient
                # 注意1：每次求单位向量
                w = direction(w) 
                if abs(f(w,X)-f(last_w,X)) < epsilon:
                    break

                cur_iter = 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_componet(X_pca,initial_w,eta)
            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

PCA 降维的基本原理:找到另外一个坐标系，这个坐标系每一个轴依次可以表达原来的样本他们的重要程度，也就是称为所有的主成分，我们取得前k个最重要的主成分，就可以将所有的样本映射到这k个轴上，获得一个低维的数据信息

6.sklearn中的PCA

import matplotlib.pyplot as plt
plt.scatter(X[:,0],X[:,1],color='b',alpha=0.5)
plt.scatter(X_restore[:,0],X_restore[:,1],color='r',alpha=0.5)

<matplotlib.collections.PathCollection at 0x1a1e701198>

7-1

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
X_train,X_test,y_train,y_test = train_test_split(X,y,random_state = 666)

使用64个维度的数据集，训练knn算法

%%time
from sklearn.neighbors import KNeighborsClassifier

knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train,y_train)

CPU times: user 2.65 ms, sys: 1.46 ms, total: 4.11 ms
Wall time: 2.65 ms

knn_clf.score(X_test,y_test)

0.9866666666666667

使用PCA进行降维，然后再训练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)

%%time
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction,y_train)

CPU times: user 1.67 ms, sys: 935 µs, total: 2.61 ms
Wall time: 1.67 ms

# 从64个维度降到两个维度以后，虽然运行速度提高了，但是识别精度大大降低了
knn_clf.score(X_test_reduction,y_test)

0.6066666666666667

explained_variance_ratio_-解释方差的比例

0.14566817 代表第一个主成分可以解释14%的原数据
0.13735469 代表第二个主成分可以解释13%的原数据
两个主成分加起来可以解释百分之27的原数据，而其他的信息丢失了

可以使用explained_variance_ratio_这个参数来查看每个主成分所解释的原数据，来判断要取多少个主成分

pca.explained_variance_ratio_

array([0.14566817, 0.13735469])

pca = PCA(n_components=64)
pca.fit(X_train)
# 这个数据可以近乎表示每个主成分轴的重要程度
pca.explained_variance_ratio_

array([1.45668166e-01, 1.37354688e-01, 1.17777287e-01, 8.49968861e-02,
5.86018996e-02, 5.11542945e-02, 4.26605279e-02, 3.60119663e-02,
3.41105814e-02, 3.05407804e-02, 2.42337671e-02, 2.28700570e-02,
1.80304649e-02, 1.79346003e-02, 1.45798298e-02, 1.42044841e-02,
1.29961033e-02, 1.26617002e-02, 1.01728635e-02, 9.09314698e-03,
8.85220461e-03, 7.73828332e-03, 7.60516219e-03, 7.11864860e-03,
6.85977267e-03, 5.76411920e-03, 5.71688020e-03, 5.08255707e-03,
4.89020776e-03, 4.34888085e-03, 3.72917505e-03, 3.57755036e-03,
3.26989470e-03, 3.14917937e-03, 3.09269839e-03, 2.87619649e-03,
2.50362666e-03, 2.25417403e-03, 2.20030857e-03, 1.98028746e-03,
1.88195578e-03, 1.52769283e-03, 1.42823692e-03, 1.38003340e-03,
1.17572392e-03, 1.07377463e-03, 9.55152460e-04, 9.00017642e-04,
5.79162563e-04, 3.82793717e-04, 2.38328586e-04, 8.40132221e-05,
5.60545588e-05, 5.48538930e-05, 1.08077650e-05, 4.01354717e-06,
1.23186515e-06, 1.05783059e-06, 6.06659094e-07, 5.86686040e-07,
1.71368535e-33, 7.44075955e-34, 7.44075955e-34, 7.15189459e-34])

绘制曲线观察取前i个主成分的时候，所能解释的原数据比例

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])])

[<matplotlib.lines.Line2D at 0x1a1fc6eef0>]

7-2

sklearn中的PCA算法支持传入一个小于1的数来表示我们希望能解释多少比例的主成分

pca = PCA(0.95)
pca.fit(X_train)
# 说明前28个主成分表示了百分之95的信息
pca.n_components_

X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)

%%time
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction,y_train)

CPU times: user 2.43 ms, sys: 1.21 ms, total: 3.63 ms
Wall time: 2.55 ms

# 虽然训练出来的精度丢失了一些，但是效率却大大提高了
knn_clf.score(X_test_reduction,y_test)

0.98

PCA降维到两维的意义-可以方便可视化展示，帮助人们理解

下图每个颜色代表一个数字在降维到二维空间中的分布情况
仔细观察后可以发现，很多数字的区分还是比较明细的
比如如果只是区分蓝色的数字和紫色的数字，那么使用二个维度就足够了

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)