7-2 使用梯度上升法解决主成分分析问题

主成分分析

转化为：

转化为：
得到的是一个1xN的向量，而我们需要的是一个Nx1的向量。所以需要将得到的结果转置一下，最终结果为

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

第一主成分

import numpy as np
import matplotlib.pyplot as plt

X = np.empty((100,2))#100个样本，每个样本有两个特征
X[:,0] = np.random.uniform(0.,100.,size=100)#特征一
X[:,1] = 0.75 * X[:,0] + 3. + np.random.normal(0,10.,size=100)#特征二

X = np.empty((100,2))#100个样本，每个样本有两个特征

X[:,0] = np.random.uniform(0.,100.,size=100)#特征一
X[:,1] = 0.75 * X[:,0] + 3. + np.random.normal(0,10.,size=100)#特征二
plt.scatter(X[:,0],X[:,1])
plt.show()

输出图片：

demean

def demean(X):
	#在行的方向上求平均值，即求每一列的平均值，得到一个向量，X中的每一个样本都减去这个向量
    return X - np.mean(X,axis=0)
    
X_demean = demean(X)
plt.scatter(X_demean[:,0],X_demean[:,1])
plt.show()

输出图片：

np.mean(X_demean[:,0])#验证在第0个维度上的均值
输出：-1.5774048733874226e-14

np.mean(X_demean[:,1])#验证在第1个维度上的均值
输出：-1.6626700016786343e-14

梯度上升法求解数据的主成分（这里求的是第一主成分）

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 + eta * gradient
        w = direction(initial_w)

        if(abs(f(w,X) - f(last_w,X)) < epsilon):
            break

        cur_iter += 1

    return w

注意：
1、 w为0本来就是一个极值点，只不过这个极值点对应的是最小值的位置，但是我们求的是最大值的位置。在这个最小值的位置上它的梯度值也为0，所以此时把零向量代进去是不可以的，所以开始的时候随机化一个向量就可以了
2、不能用StandardScaler标准化数据，因为PCA过程本来就是要求一个轴，使得所有的样本映射到那个轴之后样本的方差最大，一旦我们将样本的数据进行标准化了，样本的方差就为1了，方差的最大值就不存在了，因为在标准化的过程中把样本的方差给打掉了，这样子的话就求不出来PCA真正想最大化的结果了

initial_w = np.random.random(X.shape[1])
initial_w
输出：array([0.72295935, 0.94989331])

eta = 0.001

gradient_ascent(df_debug,X_demean,initial_w,eta)
输出：array([0.6056356 , 0.79574212])

gradient_ascent(df_math,X_demean,initial_w,eta)
输出：array([0.6056356 , 0.79574212])

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')#因为w是个单位向量，我们可以放大30倍观察
plt.show()

输出图片：

 

求数据的前n个主成分

问：求出第一主成分以后如何求出下一个主成分？
答：将数据进行改变，将数据在第一个主成分上的分量去掉，然后在新的数据上求第一主成分

即：将$X^{(i)}$在$X_{project}^{(i)}$上的分量去掉，即：将$X^{(i)}$在映射在和$X_{project}^{(i)}$这个轴相垂直的一个轴上

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