5 主成分分析PCA

最新推荐文章于 2025-06-03 20:18:29 发布

weixin_33725270

最新推荐文章于 2025-06-03 20:18:29 发布

阅读量121

点赞数

CC 4.0 BY-SA版权

原文链接：https://juejin.im/post/5cf5471f51882576be275b9d

主成分分析(PCA)是最常见的降维算法。

PCA是非监督的机器学习算法
主要用于数据的降维
其他应用：可视化、去噪

1 什么是PCA

1.1 分析

上图为数据在2个特征维度的分布情况，分别垂直到X轴、Y轴，这是两个维度，从上图中来看，无法说明这两个维度哪个更好。那么能否找到一个轴，使得投射到此轴的上的点，间距要比投射到X轴、Y轴间距更大（区分度更高）。如下图：

那么如何找到这个让样本间间距最大的轴？
如何定义样本间间距？

使用方差， $Var(x) = \frac{1}{m}\sum_1^m(x_i - \overline{x})^2$

1.2步骤

将样本的均值归0（demean），转化后如下图，此时 $Var(x) = \frac{1}{m}\sum_1^mx_i ^2$

求一个轴的方向,使得所有样本映射到w以后，下面值更大

Var(x) = \frac{1}{m}\sum_1^m(X_{project}^{(i)} - \overline{X}_{project})^2

归0化之后， $\overline{X}_{project} = 0$ ，此时转化为求下面式子最大值

Var(x) = \frac{1}{m}\sum_1^m||X_{project}^{(i)}||^2

令 w = (w1,w2) 为一个单位向量

1.3 目标函数

目标：求w，使得下值最大

Var(x) = \frac{1}{m}\sum_1^m(X^{(i)}.w)^2 = \frac{1}{m}\sum_1^m(X_1^{(i)}.w_1 + ... +X_n^{(i)}.w_n)^2

一个目标函数的最优化问题，使用梯度上升法解决

1.4 PCA与线性回归的区别

PCA是找到一个轴，使得样本空间点映射到这个轴后方差最大。
线性回归尝试的是最小化预测误差。
线性回归的目的是预测结果，而主成分分析不作任何预测。

上图中，左边的是线性回归的误差（垂直于横轴投影），右边则是主要成分分析（垂直于红线投影）。

2 梯度上升法求PCA问题

2.1 推导过程

目标：求w，使得下值最大

f(X) = \frac{1}{m}\sum_1^m(X_1^{(i)}.w_1 + ... +X_n^{(i)}.w_n)^2

求梯度▽f，转化为向量点乘的过程：w是列向量，X是m*n维矩阵

2.2 代码实现

def demean(X):
    """均值归0化"""
    # 减去每一列的均值
    return X - np.mean(X, axis=0)

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):
    """梯度调试，因为w是单位向量，比较小，因此epsilon也应该小一些"""
    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):
    """让w为单位向量"""
    return w / np.linalg.norm(w)

def first_component(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):
    """梯度上升法求w"""
    # 几个注意点：
    # 1、w初始向量不能为0向量   
    # 2、不能使用StandardScalar 标准化数据 ，因为我们要求这个轴的最大方差，如果使用标准化，方差就为1了，就不存在最大方差了
    w = direction(initial_w)  # 注意，w初始向量不能为0向量
    cur_iter = 0

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

    return w
复制代码

2.3 前n个主成分

求出第一主成分之后，如何求下一个主成分？答：将数据在第一个主成分上的分量去掉

如下图： $X^{`(i)}$ 即为下一个主成分

def first_n_components(n, X, eta=0.01, n_iters = 1e4, epsilon=1e-8):
    """求前n个主成分"""
    X_pca = X.copy()
    X_pca = demean(X_pca) # 归0化
    res = [] # 记录前n个主成分
    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
复制代码

3 高维数据向低维数据映射

3.1 映射分析

W_k 为前k个主成分， X.W_k^T = X_k ,就完成了n个维度向k个维度的映射

当然，也可以反向映射，由低维数据映射到高维数据，但是此时会有信息丢失

3.2 PCA的封装

import numpy as np

class PCA:

    def __init__(self, n_components):
        """初始化PCA"""
        assert n_components >= 1
        self.n_components = n_components
        self.components_ = None # 前n个主成分

    def fit(self, X, eta=0.01, n_iters=1e4):
        """获得数据集X的前n个主成分"""
        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):
            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

复制代码

3.3 降到多少维才好

我们先看下问题的产生过程：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.decomposition import PCA

digits = datasets.load_digits()
X = digits.data
y = digits.target
# X_train.shape = (1347, 64)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)

# %%time 使用kNN算法，时间64.5 ms
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
knn_clf.score(X_test, y_test) #0.98666666666666669

"""使用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  时间为2.93 ms，时间大大节省
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)
knn_clf.score(X_test_reduction, y_test)  #0.60666666666666669
复制代码

原来有64维度，如果降到2维，时间虽然提高了，但是精度却降低的太多了，那么应该降维到多少为好呢？

sklearn 中提供了一个方法

pca.explained_variance_ratio_   # array([ 0.14566817,  0.13735469])
# 表示2个主成分 分别能代表14.5%的方差和13.7%的方差
复制代码

下面打印所有维度能代表方差的维度

from sklearn.decomposition import PCA

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()
复制代码

sklearn 中又提供了一个方法，输入参数为保留多少方差

pca = PCA(0.95)
pca.fit(X_train)
pca.n_components_  # 28

X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
# %%time  19.7 ms
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)
knn_clf.score(X_test_reduction, y_test) #0.97999999999999998
复制代码