学习笔记： 机器学习经典算法-主成分分析PCA与梯度上升法

算法特点
非监督机器学习算法，主要用于数据降维；降维可以提高算法效率，同时帮助可视化，以便于人类理解更好的理解数据；去噪。

1、 PCA的基本原理

样本在其特征空间的分布表现是其各特征轴上记录信息的综合呈现形式。 PCA 分析基于能够捕获原始样本最大信息量为目标在样本的原始特征空间寻找一个新的坐标系，并通过在新坐标系中挑选能够解释样本绝大部分信息的前$n$个轴来描述样本，从而完成数据的 降维。

1.1 二维空间的降维

对于包含二维特征的样本集，如果要将这个数据集从二维平面降到一维空间，通过在样本的二维特征空间中寻找一条直线$w$，将空间上的所有样本点 向直线作投影，使得直线上的 投影点间距离相对较大同时拥有较高的区分度，从而在直线所代表的一维空间上，投影点的分布与原二维空间上样本点的空间分布尽可能相似，这条一维直线上的投影点就是二维空间样本的降维结果。

在寻找空间中最优的样本投影直线过程中，样本间距 由统计学中的描述 数据的离散程度(疏密) 的 方差(Variance) 所度量；方差大样本分布越稀疏，方差小样本分布越紧密。所以最优的投影轴的判断条件即使得样本空间的所有点映射到这条轴后，方差最大。同时，由于这条投影轴上的投影点分布具有最大方差，称为 第一主成分。
$$Var(x) = \frac {1}{m}\sum _{i=1}^{m}(x_i-\overline x)^2$$

2、 PCA 降维操作的基本步骤

• Step.1 将样例 均值归零(demean) ，$x_{i} - \overline x$；该过程不改变样本分布。但样例的方差计算将简化为:

$$\overline x = 0, Var(x) = \frac {1}{m}\sum _{i=1}^{m}x_i^2$$

• Step.2 在 demean 后的样本空间中寻找一个方向为$(w_1,w_2)$的直线$\vec w$（除了方向，投影直线的大小以及偏移都不影响投影点间距，所以直接从中心化后的样本空间坐标原点确定直线的方向即可），使得所有样本映射到$\vec w$所在直线以后，投影点间方差最大:

$$Var(X_{project}) = \frac {1}{m}\sum _{i=1}^{m} {(x^{(i)}_{project} - \overline x_{project})^2}$$
$\because$ 这里由于样例中心为原点，所以$\overline x_{project} = \overline x = 0$
$\therefore Var(X_{project}) = \frac {1}{m}\sum _{i=1}^{m} {\|x^{(i)}_{project} \|^2}$
$\because$ 根据 向量投影 ，结合方向向量 $\vec w$ 的模为 1可有 $\|x^{(i)}_{project}\| = x^{(i)} \cdot \vec w$

$\therefore Var(X_{project}) = \frac {1}{m}\sum _{i=1}^{m} {(x^{(i)}\cdot \vec w )^2}$
推广到$n$维样本点，展开为: $Var(X_{project}) = \frac {1}{m}\sum _{i=1}^{m} {(x^{(i)}_{1}\cdot \vec w_{1} + ... + x^{(i)}_{n}\cdot \vec w_{n})^2}$

• Step.3 求解使得目标函数 $Var(X_{project}) $ 取最大值对应的$\vec w$，从而获得最优的投影坐标轴。其中使用搜索的方式 进行求解的方案为 梯度上升法 。

目标：求$\vec w$，使得$f(x) = \frac {1}{m}\sum _{i=1}^{m} {(x^{(i)}_{1}\cdot \vec w_{1} + ... + x^{(i)}_{n}\cdot \vec w_{n})^2}$ 最大；
注意：PCA的 目标函数形成的曲面 本身没有极大值存在，直接使用梯度法求解的时候并不会收敛，这里前提是基于约束条件进行的梯度法求解。所以这是一个 带约束的求极值问题。这个约束就是，$\vec w$的模要为1。PCA本身具有非常好的 数学解，梯度法搜索并不是最好的求解策略，但可以用于理解PCA过程。
求关于$\vec w$每个维度的梯度(偏导)

$\nabla f= \begin{bmatrix} \partial f/ \partial w_{1} \\ \partial f/ \partial w_{2} \\ ... \\ \partial f/ \partial w_{n} \end{bmatrix} =\frac{2}{m} \begin{bmatrix} \sum_{i=1}^{m}{(x^{(i)}_{1}\cdot \vec w_{1} + ... + x^{(i)}_{n}\cdot \vec w_{n})\cdot (x^{(i)}_1)} \\ \sum_{i=1}^{m}{(x^{(i)}_{1}\cdot \vec w_{1} + ... + x^{(i)}_{n}\cdot \vec w_{n})\cdot (x^{(i)}_2)}\\... \\ \sum_{i=1}^{m}{(x^{(i)}_{1}\cdot \vec w_{1} + ... + x^{(i)}_{n}\cdot \vec w_{n})\cdot (x^{(i)}_n)} \end{bmatrix} = \frac{2}{m} \begin{bmatrix} \sum_{i=1}^{m}{(x^{(i)}\cdot \vec w)\cdot x^{(i)}_1} \\ \sum_{i=1}^{m}{(x^{(i)}\cdot \vec w)\cdot x^{(i)}_2} \\... \\ \sum_{i=1}^{m}{(x^{(i)}\cdot \vec w)\cdot x^{(i)}_n} \end{bmatrix} = (\frac{2}{m} \cdot (X \cdot \vec w)^{T} \cdot X) ^{T} = \frac{2}{m} \cdot X^{T} \cdot (X\vec w)$

2.2 PCA与数据归一化

• 归一化 (Standar Scaler) 处理样本的特征数据，会非线性地改变数据的特征量纲，从而进行主成分分析时需要考虑实际需求选择性进行数据的归一化。比如只是利用PCA对原始数据进行降噪不降维的处理，则不能进行归一化的数据预处理，只需要进行数据平移(如 demean 操作)即可。
• 对于需要无偏差考量不同特征对样本输出标记贡献率的分析任务，则需要注意样本特征数量级带来的方差影响影响（如一个特征的量纲从km更改为cm就会放大它的方差），其中具有较大方差的特征将主导第一主成分的方差贡献。因此对每个特征进行归一化后再进行PCA降维处理才变得有意义（比如在聚类分析上，样本的相似性度量就需要在同一尺度上来处理特征对样本输出标记的影响）。

2.3 梯度上升的简单pyhton过程

import numpy as np
### Prepare dataset
X = np.empty((100,2))
X[:,0] = np.random.uniform(0.,100.,size = 100)
X[:,1] =0.75 * X[:,0] + np.random.normal(0.,10.,size = 100)
### 特征均值归零
X = np.hstack([np.array(X[:,0] - np.mean(X[:,0])).reshape((-1,1)),np.array(X[:,1] - np.mean(X[:,1])).reshape((-1,1))])


### 目标函数
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))

### 梯度上升求解过程
w = np.random.randint(10,size = X.shape[1])+ np.random.normal(size = X.shape[1]) ### 初始化参数集
w =  w / np.linalg.norm(w) ### 约束w为单位向量
for cur_iter in range(int(1e4)):
    last_w = w
    w = w + .001 * df(w,X)
    w = w / np.linalg.norm(w) ### 约束w为单位向量
    if(abs(f(w,X) - f(last_w,X)) < 1e-8):break

print(w)

2.4 数据的前$n$个主成分

PCA 工作的本质是 对特征空间的坐标轴重新排列，在样本的原始特征空间内寻找一套新的坐标轴替换掉样本的原始坐标轴；在这套新的坐标轴里，第一主成分轴捕获了样本最大的方差，第二主成分轴轴次之，第三主成分轴再次之，以此类推；

所以，对属于$n$维特征空间的样本，如果仅仅找出包含它们投影最大方差的 第一主成分 轴($\vec w$)，在该轴上的投影向量结果 $\vec x^{(i)}_{project} = \| x^{(i)}_{project}\|\cdot \vec w = x^{(i)} \vec w \cdot \vec w$ 并不是一维数据， 依旧还是$n$维向量，也就并没有完成数据的降维。为了完成 降维 ,有必要继续寻找其特征空间的其它$n-1$个主成分。

在 第一主成分 的前提下，计算 第二组主成分。
以 向量投影 一章的内容作为铺垫。当 样本向量$\vec c$ 向第一主成分轴$\vec w$ 进行投影后，在第一主成分轴上的 投影向量$\vec b$ 仅包含了原样本向量的部分信息，其中未被捕获的信息由 原样本向量 在 投影向量垂直轴上 的 投影向量$\vec a$ 所捕获。

因此，在求算 第二主成分 时，原始向量 需要先除去被第一主成分捕获的信息，即将数据在第一主成分上的分量减掉; 然后对去掉第一主成分分量的原始样本向量继续找使得它们方差最大的投影直线，这条直线所在的轴就是这组数据的 第二主成分， 该过程就是 求PCA第一主成分过程的重复。

• 总结
• Step.1 去掉 第一主成分 分量 $\vec x^{'(i)} = \vec x^{(i)} - \vec x^{(i)}_{project}$ ；
• Step.2 求新样本向量的第一主成分；
• 求解第$i$个主成分也是类似该过程的重复.

def cur_component(X,w,eta,epsilon):
    w =  w / np.linalg.norm(w) ### 约束w为单位向量
    for cur_iter in range(int(1e4)):
        last_w = w
        gradient = X.T.dot(X.dot(w)) *(2./len(X)) ### 求梯度
        w = w + eta * gradient
        w = w / np.linalg.norm(w) 
        target_f_0 = np.sum(X.dot(last_w)**2) / len(X) ### 目标函数值，求方差
        target_f_1 = np.sum(X.dot(w)**2) / len(X) ### 目标函数值，求方差
        if(abs(target_f_0 - target_f_1) < epsilon):break
            
    return w

    
def n_components(n,X,eta = .01,n_iters = 1e4 , epsilon = 1e-8):
    X_use = X.copy()
    X_use = X_use - np.mean(X_use,axis=0) ### demean
    res = [] ### 存储每个主成分结果
    for i in range(n):
        initial_w = np.random.random(X_use.shape[1])
        w = cur_component(X_use,initial_w,eta,epsilon) ### 计算当前主成分
        X_use = X_use - X_use.dot(w).reshape((-1,1))*w   ### 除掉上一主成分上的分量
        res.append(w)
    return res



#### Prepare data
import numpy as np
X = np.empty((100,2))
X[:,0] = np.random.uniform(0.,100.,size = 100)
X[:,1] =0.75 * X[:,0] + np.random.normal(0.,10.,size = 100)
#### Test method
n_components(2,X)### 求前两个主成分轴方向

2.5 高维数据向低维数据的映射

当在样本的$n$维特征空间中，按 方差最大化 搜索了前$k(k \lt n)$个方向轴，使用这$k$个轴来描述样本分布完成数据的降维。

根据向量投影，它在投影轴上的 投影模长 ($\| x^{(i)}_{project}\| = x^{(i)} \vec w $) 就是它在投影轴上的分布，也就说是这是 向量对象 脱离原始空间在 投影轴生成空间 内的描述坐标。原始数据到$k$个方向轴的映射（原始坐标系到新坐标系的转换）为:
$$X = \begin{bmatrix} x_{1}^{(1)}&x_{2}^{(1)}&...&x_{n}^{(1)} \\ x_{1}^{(2)}&x_{2}^{(2)}&...&x_{n}^{(2)}\\ ...&...&...&... \\ x_{1}^{(m)}&x_{2}^{(m)}&...&x_{n}^{(m)}\end{bmatrix},W_{k} = \begin{bmatrix} w_{1}^{(1)}&w_{2}^{(1)}&...&w_{n}^{(1)} \\ w_{1}^{(2)}&w_{2}^{(2)}&...&w_{n}^{(2)}\\ ...&...&...&... \\ w_{1}^{(k)}&w_{2}^{(k)}&...&w_{n}^{(k)}\end{bmatrix} \rightarrow X_{k} = X\cdot W^{T}$$
数据的维度恢复
$$X_{k\_(m*k)} \cdot W_{k\_(k*n)} = X_{recover\_(m*n)} \neq X_{\_(m*n)}$$
数据维度恢复在数学运算上虽然成立。但是降维后的数据$X_k$，意味着原始数据$X$发生了 信息丢失 ，所以通过前$k$个主成分来恢复到$n$维的话，恢复出来的数据对象并不是降维前的对象.

PCA降维与数据噪音
从噪音角度出发，对于降维造成的 信息丢失 ，这部分丢失的信息存在于解释方差比非常小的成分轴，这些丢失的信息很大一部分对于描述对象来说属于 噪音。原始数据在测量过程中，由于仪器误差、测量方法缺陷、记录错误等因素，都会造成现实世界采集的数据存在噪音。通过 PCA降维去除这些噪音，能有效提升样本对象进行识别分类回归等分析时的准确率。

利用 维度恢复 方法将降维去噪后的数据对象恢复到原始维度之后，相当于在不损失数据特征数量的前提下对原始数据完成了 去噪 。这在图像识别处理上应用较多。

3、scikit-learn 框架内封装的PCA

from sklearn.decomposition import PCA
pca = PCA(n_components= 2) ### sklearn内的pca算法为pca的正规数学解
pca.fit(X)
pca.components_  ### 主成分方向
pca.transform(X) ### 根据n_components降维

PCA的主成分轴方向向量（pca.components_ ）内各分量的大小代表了对应特征对该主成分的 方差贡献，或者说该方向轴捕获了方向轴向量内最大分量对应特征的最多信息。

3.2 获取$n$个主成分

方法属性 ： pca.explained_variance_ratio_ 返回每个主成分所解释 样本总方差 的比例（特征空间内样本的总方差 为样本各特征上的方差之和）。解释方差比也可理解为捕获的信息量。

初始化PCA方法对象时如果不指定参数输入一个浮点数如 pca = PCA(0.95) 来表示降维后希望能够保留下的 信息量，sklearn-pca 方法将挑选相应主成分数目进行降维。

pca = PCA(0.95) 
pca.fit(X)
pca.transform(X)
pca.n_components_ ### 保留用于降维的主成分个数

3.3 数据集测试

### Prepare data
import numpy as np
from sklearn import datasets
digits = datasets.load_digits() ### 载入手写数字数据集
X = digits.data
y = digits.target
### test_train_split
from sklearn.model_selection import train_test_split
Train_X,Test_X,Train_Y,Test_Y = train_test_split(X,y,test_size= 0.2,random_state=666)

###### scikit-learn PCA reducion #####
from sklearn.decomposition import PCA
pca = PCA(0.8)  ### 降维到仅保留原始样本80%信息
pca.fit(Train_X)
Train_X_reduc = pca.transform(Train_X) ### 根据n_components降维训练集
Test_X_reduc  = pca.transform(Test_X ) ### 根据n_components降维测试集
print(pca.n_components_) ### 保留用于降维的主成分个数

###### scikit-learn kNN classifier #####
from sklearn.neighbors import KNeighborsClassifier
kNN_clf = KNeighborsClassifier(n_neighbors=3) 
kNN_clf.fit(Train_X_reduc,Train_Y)
kNN_clf.score(Test_X_reduc,Test_Y)  ### 预测准确率

3.4 PCA 去噪

降噪不降维数学过程：
$X_{n} \rightarrow PCA_{reduction} \rightarrow X_{k} \rightarrow X_{k}\cdot W_{k} = X_{n}^{'}$

### Prepare dataset
import numpy as np
X = np.empty((100,2))
X[:,0] = np.random.uniform(0.,100.,size = 100)
X[:,1] = 0.75 * X[:,0] + np.random.normal(0.,4.,size = 100) ### 给y轴添加一个方差为4的正态噪音

###### scikit-learn PCA Noise reduction #####
from sklearn.decomposition import PCA
pca = PCA(.5)  ### 降维到仅保留原始样本50%信息
pca.fit(X) 
X_reduc = pca.transform(X) ### 对原始矩阵进行降维，该过程伴随维度降低和消除噪音
X_restore = pca.inverse_transform(X_reduc) ### 恢复降噪降维对象到原始对象的维度，从而达到仅降噪的效果