PCA与梯度上升法

1、什么是PCA

主成分分析(Principal Component Analysis):

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

假设存在一根直线，将所有的点都映射在该条指直线上，这样的话点的整体分布和原来的点的分布就没有很大的差异（点和点的距离比映射到x轴或者映射到y轴都要大，区分度就更加明显），与此同时所有的点都在一个轴上（理解成一个维度），虽然这个轴是斜着的。用这种方式将二维降到了一维度

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

如何定义样本间间距？ 使用方差(Variance)

方差越大代表样本之间越稀疏，方差越小代表样本之间越紧密。

移动坐标轴，使得样本在每一个维度均值都为0：

demean之后的方差最大其实就是求映射后每个点到(0,0)的距离最大再求和，假设降维后轴的方向是w=(w1, w2)，Xi是映射前的向量，Xi(project)是映射后的向量

目标函数即：

对目标函数求梯度：

转化为：

转化为：

由于最终转换的结果是一个1行m列的矩阵，而我们想要得到一个n行1列的矩阵，所以还要进行一次转置

2、梯度上升法求解PCA

2.1 求数据的第一个主成分

import numpy as np
import matplotlib.pyplot as plt

X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 2 + np.random.normal(0, 10., size=100)


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

# 均值归零函数
def demean(X):
    return X - np.mean(X, axis=0)


X_demean = demean(X)


# print(X_demean)
# print(np.mean(X_demean, 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.001):
    res = np.empty(len(w))
    for i in range(len(w)):
        w_1 = w.copy()
        w_2 = w.copy()
        w_1[i] += epsilon
        w_2[i] -= epsilon
        # 2*epsilon必须用括号包起来，否则就不是除以2倍的epsilon而是除以2乘以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)
    i_iter = 1
    while i_iter < n_iters:
        gradient = df(w=w, X=X)
        last_w = w
        # gradient是对每个维度求偏导得到的列表，如果偏导数为负则w的这个维度加上一个负值，降维后的方差趋于变大
        # 如果偏导数为正，则w的这个维度加上一个正值，降维后的方差趋于变大，因此w加上导数值，降维后的方差趋于变大
        # 在eta合适的情况下，随着循环进行，导数值逐渐趋近0，eta是常数，降维后的方差的变化量会越来越小
        w = w + eta * gradient
        w = direction(w)  # 注意1，每次求一个单位向量
        # abs求绝对值
        if (abs(f(w=w, X=X) - f(w=last_w, X=X)) < epsilon):
            break
        i_iter += 1
    return w


initial_w = np.random.random(X.shape[1])  # 注意2：不能用0向量开始
eta = 0.0001
# print(gradient_ascent(df_debug, X=X_demean, initial_w=initial_w, eta=eta))
print(gradient_ascent(df_math, X=X_demean, initial_w=initial_w, eta=eta, n_iters=100))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

注意：

1、 w为0本来就是一个极值点，只不过这个极值点对应的是最小值的位置，但是我们求的是最大值的位置。在这个最小值的位置上它的梯度值也为0，所以此时把

零向量代进去是不可以的，所以开始的时候随机化一个向量就可以了

2、不能用StandardScaler标准化数据，因为PCA过程本来就是要求一个轴，使得所有的样本映射到那个轴之后样本的方差最大，一旦我们将样本的数据进行标准

化了，样本的方差就为1了，方差的最大值就不存在了，因为在标准化的过程中把样本的方差给打掉了，这样子的话就求不出来PCA真正想最大化的结果了

2.2 求数据的前n个主成分

根据推导过程可以发现，去除第一个主成分剩下的向量X’和第一主成分向量是垂直的关系，而两个垂直的向量相乘等于0，所以可以根据这个特性去验证我们求出的结果是否正确。

代码实现：

import numpy as np
import matplotlib.pyplot as plt

X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 2 + np.random.normal(0, 10., size=100)


# 均值归零函数
def demean(X):
    return X - np.mean(X, axis=0)


X_demean = 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)
    i_iter = 1
    while i_iter < n_iters:
        gradient = df(w=w, X=X)
        last_w = w
        # gradient是对每个维度求偏导得到的列表，如果偏导数为负则w的这个维度加上一个负值，降维后的方差趋于变大
        # 如果偏导数为正，则w的这个维度加上一个正值，降维后的方差趋于变大，因此w加上导数值，降维后的方差趋于变大
        # 在eta合适的情况下，随着循环进行，导数值逐渐趋近0，eta是常数，降维后的方差的变化量会越来越小
        w = w + eta * gradient
        w = direction(w)  # 注意1，每次求一个单位向量
        # abs求绝对值
        if (abs(f(w=w, X=X) - f(w=last_w, X=X)) < epsilon):
            break
        i_iter += 1
    return w


initial_w = np.random.random(X.shape[1])  # 注意2：不能用0向量开始
eta = 0.01
w = first_component(X=X, initial_w=initial_w, eta=eta)
print(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

# plt.scatter(X2[:,0],X2[:,1])
# plt.show()

w2 = first_component(X2, initial_w=initial_w, eta=eta)
print(w2)

print(w.dot(w2))

# 求前n个主成分
def first_n_component(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_component(X=X_pca, initial_w=initial_w, eta=eta)
        res.append(w)

        X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
    return res

print(first_n_component(2, X))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

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

编写PCA.py文件

import numpy as np

class PCA:

    def __init__(self, n_components):
        self.n_components = n_components
        # 数据的前n个主成分，公式中的Wk
        self.components_ = None

    # 取得数据的前n个主成分
    def fit(self, X, eta=0.01, n_iters=1e4):
        # 均值归零函数
        def demean(X):
            return X - np.mean(X, axis=0)

        X_demean = 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)
            i_iter = 1
            while i_iter < n_iters:
                gradient = df(w=w, X=X)
                last_w = w
                # gradient是对每个维度求偏导得到的列表，如果偏导数为负则w的这个维度加上一个负值，降维后的方差趋于变大
                # 如果偏导数为正，则w的这个维度加上一个正值，降维后的方差趋于变大，因此w加上导数值，降维后的方差趋于变大
                # 在eta合适的情况下，随着循环进行，导数值逐渐趋近0，eta是常数，降维后的方差的变化量会越来越小
                w = w + eta * gradient
                w = direction(w)  # 注意1，每次求一个单位向量
                # abs求绝对值
                if (abs(f(w=w, X=X) - f(w=last_w, X=X)) < epsilon):
                    break
                i_iter += 1
            return w

        X_pca = demean(X)
        # 数据的前n个主成分，公式中的Wk
        self.components_ = np.empty(shape=(self.n_components, X_pca.shape[1]))
        for i in range(self.n_components):
            initial_w = np.random.random(X_pca.shape[1])
            w = first_component(X=X_pca, initial_w=initial_w, eta=eta, n_iters=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映射到各个主成分上去，结果就是公式中的Xk
        return X.dot(self.components_.T)

    def inverse_transform(self,X):
        # 将X反向映射回原来的特征空间，参数是公式中的Xk，结果是公式中的Xm
        return X.dot(self.components_)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

测试：

import numpy as np
import matplotlib.pyplot as plt

from common.PCA import PCA

X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 2 + np.random.normal(0, 10., size=100)

pca = PCA(n_components=1)
pca.fit(X=X)
print(pca.components_)

pca_reduction = pca.transform(X)
# print(pca_reduction)

X_restore = pca.inverse_transform(pca_reduction)
print(X_restore)
plt.scatter(X[:,0], X[:,1], color='r')
plt.scatter(X_restore[:,0], X_restore[:,1])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

3、scikitlearn中的PCA

scikitlearn中的PCA使用的不是梯度上升法求解出来的，而是数学方法。

3.1 scikitlearn中PCA基本使用

import numpy as np
import matplotlib.pyplot as plt

from sklearn.decomposition import PCA

X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 2 + np.random.normal(0, 10., size=100)

pca = PCA(n_components=1)
pca.fit(X=X)
X_reduction = pca.transform(X=X)
print(X_reduction.shape)
X_restore = pca.inverse_transform(X=X_reduction)

plt.scatter(X[:,0], X[:,1], color='r')
plt.scatter(X_restore[:,0], X_restore[:,1])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

3.2 用scikitlearn中PCA分析digits数据集

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
import time

digits = datasets.load_digits()
X = digits.data
y = digits.target

X_train,X_test,y_train,y_test = train_test_split(X, y, random_state=666)
# print(X_train.shape)
start_time = time.time()
knn_clf = KNeighborsClassifier(n_neighbors=5)
knn_clf.fit(X_train, y_train)
print(knn_clf.score(X_test,y_test))
end_time = time.time()
print(end_time - start_time)

start_time = time.time()
pca = PCA(n_components=41)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)

knn_clf = KNeighborsClassifier(n_neighbors=5)
knn_clf.fit(X_train_reduction, y_train)
print(knn_clf.score(X_test_reduction, y_test))

end_time = time.time()
print(end_time - start_time)
# explained_variance_ratio_涵盖了原有数据的方差的比例
# print(pca.explained_variance_ratio_)

# 查看每个特征在原数据方差的占比
# pca = PCA(n_components=X_train.shape[1])
# pca.fit(X_train)
# print(pca.explained_variance_ratio_)

# 保留占比超过0.99的前n个主成分
pca = PCA(n_components=0.99)
pca.fit(X_train)
print(pca.n_components_)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

3.3 用scikitlearn中PCA分析mnist数据集

import numpy as np
from sklearn.datasets import fetch_openml
from sklearn.neighbors import KNeighborsClassifier
from sklearn.decomposition import PCA
import time

# fetch_mldata()不再能够使用是因为其所依赖的资源不再适用。将之替换为fetch_openml()
# 需要注意的是这个替换并不是一个无缝替换。例如mnist数据集需要改为mnist_784，具体数据集是在https://www.openml.org/
mnist = fetch_openml("mnist_784")
# print(mnist)
X,y = mnist['data'],mnist['target']
# print(X.shape)

X_train = np.array(X[:60000], dtype=float)
y_train = np.array(y[:60000], dtype=float)
X_test = np.array(X[60000:], dtype=float)
y_test = np.array(y[60000:], dtype=float)
print(X_train.shape, y_train.shape)
print(X_test.shape, y_test.shape)

start_time = time.time()
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
end_time = time.time()
print("模型训练时间：", end_time - start_time)

start_time = time.time()
print(knn_clf.score(X_test, y_test))
end_time = time.time()
print("模型预测时间：", end_time - start_time)

# 使用PCA降维
pca = PCA(0.9)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)

start_time = time.time()
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)
end_time = time.time()
print("模型训练时间：", end_time - start_time)

start_time = time.time()
print(knn_clf.score(X_test_reduction, y_test))
end_time = time.time()
print("模型预测时间：", end_time - start_time)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

特征脸方法就是将PCA方法应用到人脸识别中，将人脸图像看成是原始数据集，使用PCA方法对其进行处理和降维，得到“主成分”——即特征脸，然后每个人脸都可以用特征脸的组合进行表示。这种方法的核心思路是认为同一类事物必然存在相同特性（主成分），通过将同一目标（人脸图像）的特性寻在出来，就可以用来区分不同的事物了。人脸识别嘛，就是一个分类的问题，将不同的人脸区分开来。