04机器学习--主成分分析(PCA)与梯度上升法及python实现_pca.inverse_transform

目录

①PCA概述

②梯度上升法解决PCA问题 

③第二主成分

④高维数据映射为低维数据 

⑤scikit-learn中的PCA  

---

①PCA概述

• 一个非监督的机器学习算法
• 主要用于数据的降维
• 通过降维，可以发现更便于人类理解的特征
• 其他应用：可视化、去噪

例如：一个二维的数据

经过降维（保留其中一个维度，去除另一个维度，这里保留了特征）

另一种降维方式

比较容易看出第一种降维方案比较好，因为映射在x轴上的特征之间的距离相对较大，容易区分。

第一种方案就是最好的方案吗？

再看一种情况

都映射到了一个轴上，这种方案所有的点更加趋近于原来的分布情况，这个时候的映射的区分度是最明显的，我们的目标就变成了：如何找到使样本间的间距最大的轴？

首先定义一下样本间的间距：使用方差（描述样本疏密的指标）

其实就是找一个轴，样本映射上去之后，方差达到最大值。

第一步：将样例的均值归零（所有样本都减去整体的均值）

原来：

 变成了：

样本的分布没有改变，只是变动了一下坐标轴，使得样本在每一个维度的均值都是0

这个时候方差的公式就便形成：

小结一下：主成分分析思路

对所有的样本进行均值归零处理

我们想要求一个轴（方向向量w=（w1，w2））

使得我们所有的样本，映射到w以后，有：

X是向量，所以更准确的表示：

最终等价于：

看一下数学映射关系：

 所以上面的式子还可以变形为：

 最终目标为：

既然是求最大值，那我们可以想到，后面进行优化的时候可以使用梯度上升法解决 

小区分：PCA与线性回归的样子比较像，都有点和一条直线，但是线性回归求得是回归值与样本值的差值，他与PCA的映射关系是不一样的，他们坐标轴的含义也是不一样的，下图为线性回归

②梯度上升法解决PCA问题 

求最值，常规思路：求偏导

向量化处理（可以简化编程过程，前面几章都提到过） 

最后

 代码实现：

# coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
 
# 编一组样例
X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)  # 特征1
X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 10., size=100)  # 特征2
# 绘制一下散点
# plt.scatter(X[:, 0], X[:, 1])
# plt.show()
 
 
# 均值归零
def demean(X):
    return X - np.mean(X, axis=0)  # 向量X中表示每一个特征的均值为0，这里就是对矩阵每一列求均值
 
 
X_demean = demean(X)
# 绘制一下 现在每一维度的均值基本都是0（无限趋近于0）
# plt.scatter(X_demean[:, 0], X_demean[:, 1])
# plt.show()
 
 
# 梯度上升法求解主成分
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, 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 = w + eta * gradient  # 新的w
        w = direction(w)  # 注意1：每次求一个单位方向
        if abs(f(w, X) - f(last_w, X)) < epsilon:
            break
 
        cur_iter += 1
 
    return w
 
 
# 调用一下                                  # 线性回归是可以的，因为线性回归的目标函数求导带入0时，初始梯度不为0
initial_w = np.random.random(X.shape[1])  # 注意2：不能用0向量开始，观察损失函数的导数会发现，如果初始为0那么，本身就是一个极值点，梯度为0
eta = 0.001
# 注意3： 不能使用StandardScaler标准化数据，一旦归一化，那么样本间的方差就一样了，不存在最大方差
print(gradient_ascent(df_debug, X_demean, initial_w, eta))  # 本办法观察一下结果
print(gradient_ascent(df_math, X_demean, initial_w, eta)) # 公式推导的方法
 
 
w = gradient_ascent(df_math, X_demean, initial_w, eta)
# 绘制一下
# plt.scatter(X_demean[:, 0], X_demean[:, 1])
# plt.plot([0, w[0]*30], [0, w[1]*30], color='r')
# plt.show()

结果输出：

原始数据散点图：

均值归零后：

梯度下降后的结果：(逼近法和公式法的结果是一样的，说明我们的数学推导是正确的)

最终结果可视化输出：（红色直线就是最优那条直线的方向向量，这条轴就是第一主成分）

③第二主成分

我们上面把样本点映射到了一个方向向量上，但是其实这些数据还是在二维空间，刚才的方向向量上方差最大，所以是第一主成分，但是还有一个维度我们没找出来，他就是第二主成分，对于n维数据来说也是类似的，换句话说，主成分分析法其实就是把样本点从一组数坐标系转移到另一组坐标系，上面我们只求出来对于新的坐标系的第一条坐标轴所在的方向，那么怎么求另一条轴坐在的方向呢？

思路：对数据进行改变，将数据在第一个主成分上的分量去掉 

我的理解，就是X向量减去X在第一主成分方向上映射的向量得到的就是和第一主成分方向垂直的另一个分量，如图所示，绿色的就是第二主成分，其实就是进行了一个向量的减法，也就是说我们要在绿色方向上求X的第一主成分，这样求出来的方向其实就是X的第二主成分，以此类推就可以求出高维度数据的第三主成分、第四主成分......循环往复的过程

# coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
 
# 编一组样例
X = np.empty((100, 2))
np.random.seed(0)  #随机种子
X[:, 0] = np.random.uniform(0, 100, size=100)  # 特征1
np.random.seed(0)
X[:, 1] = 0.75 * X[:, 0] + 3 + np.random.normal(0, 10, size=100)  # 特征2
 
# 均值归零
def demean(X):
    return X - np.mean(X, axis=0)  # 向量X中表示每一个特征的均值为0，这里就是对矩阵每一列求均值
X = demean(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_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 = w + eta * gradient  # 新的w
        w = direction(w)  # 注意1：每次求一个单位方向
        if abs(f(w, X) - f(last_w, X)) < epsilon:
            break
 
        cur_iter += 1
 
    return w
 
 
# 调用一下                                  # 线性回归是可以的，因为线性回归的目标函数求导带入0时，初始梯度不为0
np.random.seed(0)
initial_w = np.random.random(X.shape[1])  # 注意2：不能用0向量开始，观察损失函数的导数会发现，如果初始为0那么，本身就是一个极值点，梯度为0
eta = 0.01
 
w = first_component(X, initial_w, eta) # 第一主成分对应的方向向量
print('w的值', w)
 
# 去掉第一主成分方向的向量
X2 = X - X.dot(w).reshape(-1, 1) * w # 向量化球阀
# 也可以用循环
# X2 = np.empty(X.shape)
# 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_component(X2, initial_w, eta) # 求另一个维度对应的方向向量
print('w2的值', w2)
 
# 绘制一下
plt.scatter(X2[:, 0], X2[:, 1])
plt.plot([0, w2[0]*10], [0, w2[1]*10], color='r')
plt.show()
 
# 验证一下w其实和w2是垂直关系
print('验证是否垂直', w.dot(w2))  # 结果肯定趋近于0
 
 
# 最后总和一下
# 求前n个主成分
def first_n_components(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):
        np.random.seed(0)
        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
print('直接求前连个主成分',first_n_components(2, X))

输出结果：

第二主成分

数据结果

④高维数据映射为低维数据 

如何将样本从n维转换成k维（k个比较重要的主成分）？

Xm和原来的X是不一样的 ，我们用主成分分析法求得前k个主成分后，就可以实现从高维度向低维度的映射，反过来也可以对低维度数据映射回高维度

 
 
import numpy as np
 
# 自己封装一下PCA类
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) # Wk的转置
 
    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
 
 
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)
 
pca = PCA(n_components=1)
pca.fit(X)
 
# 降维
X_reduction = pca.transform(X)
print(X_reduction.shape) # 2特征变成2特征
 
# 恢复维度 这过程是丢失信息的  因为将为过程中丢失了恢复不了
X_restore = pca.inverse_transform(X_reduction)
print(X_restore.shape)

结果输出：

降维的思路：可以理解为，我们为数据找了一个新的坐标系，这个坐标系每一个轴依次可以表达原来样本的主成分的重要程度，我们取出前k个最重要的主成分，然后把数据映射上去，得到了一个低维的数据。

⑤scikit-learn中的PCA  

简单实现

import numpy as np
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 sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(X)
print(pca.components_)  # 这里主成分方向和上面我们自己求出来方向是相反的，不影响最后的降维结果
                        # 因为他这里用的不是梯度上升法，使用的其他的数学方法
X_reduction = pca.transform(X)
X_restore = pca.inverse_transform(X_reduction)
print(X_reduction.shape)
print(X_restore.shape)

结果输出：

手写识别数据集演示

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)
 
print('数据集维度', X_train.shape)
 
# kNN训练
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
 
print('KNN', knn_clf.score(X_test, y_test))
 
# PCA降维
from sklearn.decomposition import PCA
 
pca = PCA(n_components=2)  # 直接降维到2
pca.fit(X_train)
X_train_reduction = pca.transform(X_train) # 训练集降维
X_test_reduction = pca.transform(X_test)  # 测试集降维
 
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)
 
print('降维后的KNN', knn_clf.score(X_test_reduction, y_test))  # 降维后识别时间减少，但是识别精度也相应降低了
# 所以要找到合适的降维数
 
from sklearn.decomposition import PCA
pca = PCA(n_components=X_train.shape[1])
pca.fit(X_train)
# print(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()
 
pca = PCA(0.95) # 找到95%的还原度
pca.fit(X_train)
print('还原95%需要的维度', pca.n_components_)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)
print(knn_clf.score(X_test_reduction, y_test)) #虽然比全样本要低，但是时间性能很棒这里快了10倍

输出结果：

使用PCA对数据进行降维可视化 

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
 
# 手写识别数据集
digits = datasets.load_digits()
X = digits.data
y = digits.target
 
# PCA降维
from sklearn.decomposition import PCA
 
pca = PCA(n_components=2)  # 直接降维到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()

 结果输出：

这里想说明，上面把64维数据直接降维到2，并不没有意义的，这样做的的好处是可以进行数据可视化，会发现就算降维到2，在图中还是有许多数据有明显的差异，这样我们就可以进行一些分析，比如说我们只想分析蓝色和紫色的数据，那么我们就直接降维到2维也是可以的

注：有时候降维能去除噪音，这样准确率可能会更高！