机器学习初学者——主成分分析与梯度上升法学习笔记

主成分分析法（Princiipal Component Analysis）简称PCA，是一个非监督学习算法，用处：降维、去噪、实现数据的可视化
本质：讲一个坐标系转化为另一个同样维度的坐标系
原理：寻找一条轴，使得所有样本点垂直映射到该轴上，轴上映射的点之间间距最大，区分度最明显，以轴上样本点的方差表示来衡量样本间间距
数学算法：

• 1.样本均值归零(demean)
  目的：便于化简方差公式，
  

相当于每个方向上的坐标减去该方向上的均值，坐标轴向正的方向平移了均值的单位长度

• 2、求出轴的方向向量w=（w1,w2,……wn）（n表示样本的维度）

该方向向量满足条件：所有样本映射到w方向后，轴上的样本点之间的方差最大，X_project表示映射到最佳轴上的点的坐标（n维的向量），被减数表示各个特征的均值所组成的向量（0,0,0，……0），那么括号中减法得到的就是一个向量，向量的平方可以看成同样的两个相同向量相乘，由于夹角为0，因此相当于向量的长度相乘，也就是范数相乘。

2.1用样本点的原坐标表示映射后的坐标X_project

** 以二维为例：w=(w₁,w₂)表示最佳轴方向的的单位向量，||w||=1， Xⁱ表示样本点，X_pr表示映射后的点，那么图中蓝色的向量的长度就是||X_projectⁱ||，即所求值即为单位方向向量与样本点代表的向量相乘，**

因此，我们的目标可以转化为：

因此，主成分分析的目标函数为：

样本本身限定了目标函数的最大值，因此F(X)的最大值一定存在，这里用梯度上升法，通过求导，逐渐达到最大值。

f（X）关于w求导：

经过向量化处理结果如下

以上可知，目标函数和梯度值的表达式分别为

（一）梯度上升法实现主成分分析

根据以上数学原理：写出以下代码：
功能包括：均值归零、求梯度、、求函数值、将方向向量化为单位方向向量、梯度上升函数
注意w不能为零初始化的w值不能为0，因为w为0时，w与任何向量点乘之后都是0，程序会一直止步在最小点
均值归零

def demean(X):
    return X - np.mean(X, axis=0)

求梯度

def df(w, X):
    return X.T.dot(X.dot(w)) * 2. / len(X)

求目标函数

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

对方向向量的处理

def direction(w):
    return w / np.linalg.norm(w) //np.linalg.norm(w)表示 求向量的模    

上面这行代码不仅仅使向量的模变成了一，同时因为梯度上升法中，主要原因是计算目标函数（效用函数）和其导数时，认为w为单位向量，如果w的比较大，会造成搜索的过程不顺利，为了使搜索的结果合理，那么就要保证eta的值小，相应的会增加循环次数

梯度上升法

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(w) # 注意1：每次求一个单位方向
        if(abs(f(w, X) - f(last_w, X)) < epsilon):
            break            
        cur_iter += 1
    return w

另外为了验证求导的算法正确性，需要根据导数为割线的斜率的极限来求导，用低维小批量的样本用于验证两种求导结果是否相同。代码如下

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

(二)求前N个主成分

样本点经过第一次主成分分析之后，所得到的方向向量w代表第一主成分，由于样本有n个维度，因此，样本最多可以求出n个主成分，每次的轴的重要程度一次递减
如何求得第二个主成分：
求出第一主成分之后，将样本z在第一主成分上的分量去掉，由下图可知，数据在第一主成分上的分量为X_project，其模长为Xⁱw，那么Xⁱw*w就表示一个向量，也就是X_project的坐标。
那么X’ⁱ =Xⁱ-Xⁱ_project 就表示样本取出第一主成分的分量之后的结果，如果要求第二主成分，只需对X’ⁱ进行第一主成分的操作，之后的主成分同理可求。

代码实现：

方法一：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

def n_components(n, X, eta=0.01):
    X_pca = X.copy()
    res = []
    for i in range(n):
        initial_w = np.random.random(X_pca.shape[1])
        w = ascent_gradient(X_pca, initial_w, eta)
        res.append(w)
        
        X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
        
    return res

（三）高维数据向低维数据的映射和数据的降噪

高维数据向低维数据映射：
X表示m个样本，每个样本n个特征，W表示前k个主成分，每个主成分对应一个单位向量。

那么X与W_k的转置点乘之后得到一个mk的矩阵，表示样本在前k个主成分上映射所对应的坐标，这一过程就完成了样本从n维降为k维的过程

上图中X_k代表降维之后的数据，X_k点乘W_k会重新得到一个mn维的矩阵，但，由于在降维的过程中一部分数据被丢掉，因此新的矩阵不是原来的样本，这就是降噪的原理

代码封装：

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

（四）scikit-learn中实现主成分分析

from sklearn.decomposition import PCA
pca=PCA(n_components=???)
pca.fit(x)
之后便可以查看结果
pca.components_
x_reduction=pca.transform(x)//降维
x_restore=pca.inverse_transform(x_reduction)//得到和原来样本同样维度的数据，但是不在等于原来的数据

• 怎样判定降噪时应降至多少维既能减少噪音又能不损失过多信息？

pca.explained_variance_ratio_ //该指标表示经过样本经过降维之后，每个主成分方向向量上的样本所能解释的原样本的方差相应的比例

主成分从前到后重要程度依次递减，因此，所能解释的方差也会递减，因此，如果原样本有n个维度，那么求出n个主方向，然后查看每个主方向对样本的方差解释的程度，它们从大到小排列,加和为1，为了判断应该舍弃多少主成分，这里我们绘制一下累加曲线

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

上图分别表示一个64维的数据的64个主成分对样本的方差的解释程度和频率累加曲线，上图中，我们可以看到，当降维至大约30时，已经可以解释95%以上原样本的方差，
因此，在创建对象时，可以给PCA类输入一个0~1之间的数值，代表要解释的样本方差比例，之后便可进行其他操作

pca = PCA(0.95)
pca.fit(X_train)