机器学习算法笔记(十三):求数据的前n个主成分

上一篇文章中我们探讨了主成分分析法的基本思路以及第一主成分的求法。本文将着重探讨如何求数据的前n个主成分以及把我们的PCA算法封装起来,以便后续方便的调用。

一、求数据的前n个主成分的思路

以上一篇文章的数据为例,我们将这些数据映射到了这根红色的轴上。同时,它又是个二维的数据,它还应该有另外的一个轴(当然对一个n维数据来说它应该还有另外n-1个轴)。在主成分分析中,我们将一个n维数据的n个轴重新进行排列,使得第一个轴对于这些样本的方差是最大的,第二个轴次之,第三个轴再次之……换句话说,主成分分析法本质上是从一组坐标系转到另一组坐标系的过程。之前我们只求出了对于这个新坐标系来说第一个轴所在的方向,那么我们如何求出下一个轴所在的方向呢(求出第一主成分后,如何求出下一个主成分)?由于我们已经求出第一主成分,所以我们可以将数据进行改变,将数据在第一个主成分上的分量去掉。

所谓的“去掉”其实就是向量运算中的“减掉”,也就是上图中X(i) ,记作X’(i)(图示绿色向量)。现在我们要求第二主成分就非常简单了,就是求新数据(X’(i))上的第一主成分。这个过程以此类推,我们就可以求第三主成分,第四主成分了。

二、编码实现求数据的前n个主成分

新建一个工程,创建main.py文件,实现以下代码:

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] + 3. + np.random.normal(0, 10., size=100)

#样本均值归零
def demean(X):
    return X - np.mean(X, axis=0) #axis=0,沿着行的方向计算均值,实则求每一列均值

X_demean = demean(X)

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

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

#定义f的导数
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 + eta * gradient
        w = direction(w)
        if (abs(f(w, X) - f(last_w, X)) < epsilon):
            break

        cur_iter += 1

    return w

initial_w = np.random.random(X.shape[1]) #随机生成初始搜索点
w = first_component(X, initial_w, 0.01)
print(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, 0.01)
print(w2)

打印的两个主成分方向如下:

在计算X2的过程中,我们同样可以对其进行向量化处理从而加快运行速度,第51至53行可修改如下:

X2 = X - X.dot(w).reshape(-1, 1) * w

上面仅仅是实现了求数据前两个主成分的情形,现在我们综合上述知识,封装一个函数来实现求数据的前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):
        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

调用示例如下:

result = first_n_components(2, X) #X只有两个维度,最多有两个主成分

三、将PCA算法封装成一个类

在MyML包中新建一个PCA.py文件,将前面的代码封装起来,实现如下:

import numpy as np

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 __repr__(self):
        return "PCA(n_components=%d)" % self.n_components

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注