机器学习入门算法之主成分分析法(PCA)

一 、概述

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

(一)分析并得出目标函数

以2个特征值为例,曲线是二维平面。

要找到一个直线,使得特征值1和特征值2都拟合到直线上,这样就从二维降维到了一维, 这样就可以用拟合的这个新的拟合值同时反映原来2个特征值。同时拟合后的值尽可能的使得特征值排列稀疏,这样使得特征更明显

问题提出:如何找到这个让样本间间距最大(特征值稀疏)的直线 → 如何定义样本件间距。

答:使用方差 $\quad Var(x) = \frac{1}{m} \displaystyle\sum_{i=1}^m (x_i - \bar x)^2$

第一步:将样本数据的均值归为0。$\bar x = 0$

  $ \Rightarrow\quad Var(x) = \frac{1}{m} \displaystyle\sum_{i=1}^m (x_i - \bar x)^2 = \frac{1}{m} \displaystyle\sum_{i=1}^m x_i^2$

第二步:设想要求的一个直线的方向为$ w = (w_1, w_2)$

  使得样本数据映射到 w 以后,有:

  $Var(X_{project}) = \frac{1}{m} \displaystyle\sum_{i=1}^m ||X_{project}^{(i)} - \bar X_{project}||^2$取得最大值。($X_{project}$为样本数据映射到w轴后的点)

  因为样本数据均值为0,所以:$Var(X_{project}) = \frac{1}{m} \displaystyle\sum_{i=1}^m ||X_{project}^{(i)} ||^2$

第三步:求$Var(X_{project}) = \frac{1}{m} \displaystyle\sum_{i=1}^m ||X_{project}^{(i)} ||^2$最大值

$||X_{project}^{(i)} || = ||X^{(i)} || \cdot cos\theta = ||X^{(i)} || \cdot ||w|| \cdot cos\theta = X^{(i)} \cdot w $

$\Rightarrow Var(X_{project}) = \frac{1}{m} \displaystyle\sum_{i=1}^m (X^{(i)} \cdot w)^2 = \frac{1}{m} \displaystyle\sum_{i=1}^m (X^{(i)}_1 w_1 +X^{(i)}_2 w_2 +...+ X^{(n)}_n w_n )^2 $

(w为单位方向向量)

即得出目标函数为$f(x) = \frac{1}{m} \displaystyle\sum_{i=1}^m(X^{(i)} \cdot w)^2 $

(二)计算目标函数的梯度

目标:w,使得上述式子取得最大值,对目标函数的最优化问题,使用梯度上升法

公式化简:

$$ \nabla f = \left(\begin{matrix} \frac{\partial f}{\partial w_1} \\ \frac{\partial f}{\partial w_2} \\ \cdots \\ \frac{\partial f}{\partial w_n} \\ \end{matrix}\right) = \frac{2}{m} \left(\begin{matrix} \displaystyle\sum_{i=1}^m (X^{(i)}_1 w_1 + X^{(i)}_2 w_2 +...+ X^{(i)}_n w_n)X^{(i)}_1 \\ \displaystyle\sum_{i=1}^m (X^{(i)}_1 w_1 + X^{(i)}_2 w_2 +...+ X^{(i)}_n w_n)X^{(i)}_2 \\ \cdots \\ \displaystyle\sum_{i=1}^m (X^{(i)}_1 w_1 + X^{(i)}_2 w_2 +...+ X^{(i)}_n w_n)X^{(i)}_n \end{matrix}\right) = \frac{2}{m} \left(\begin{matrix} \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_1 \\ \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_2 \\ \cdots \\ \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_1 \\ \end{matrix}\right) $$

公式向量化:

$$ \begin{cases} \frac{2}{m} \cdot (X^{(1)} w, X^{(2)} w, X^{(3)} w,...,X^{(m)} w ) \\ \\ \left(\begin{matrix} X_1^{(1)} & X_2^{(1)} & X_3^{(1)} &\cdots& X_n^{(1)} \\ X_1^{(2)} & X_2^{(2)} & X_3^{(2)} &\cdots& X_n^{(2)} \\ X_1^{(3)} & X_2^{(3)} & X_3^{(3)} &\cdots& X_n^{(3)} \\ \cdots & \cdots & \cdots & \ddots & \cdots & \\ X_1^{(m)} & X_2^{(m)} & X_3^{(m)} &\cdots& X_n^{(m)} \\ \end{matrix}\right) \\ \end{cases} $$

$$ \Downarrow \\ = \frac{2}{m} \left(\begin{matrix} \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_1 \\ \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_2 \\ \cdots \\ \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_1 \\ \end{matrix}\right) = \frac{2}{m} \cdot (Xw)^T \cdot X \\ \Rightarrow 转置获得n×1的矩阵: \frac{2}{m} \cdot X^T \cdot(Xw) $$

最后得出结果:

$$ \nabla f = = \frac{2}{m} \left(\begin{matrix} \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_1 \\ \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_2 \\ \cdots \\ \displaystyle\sum_{i=1}^m (X^{(i)} w)X^{(i)}_1 \\ \end{matrix}\right) = \frac{2}{m} \cdot X^T \cdot(Xw) $$

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

(一)模拟实现

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

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

2. 进行demean操作(把样本数据均值归为0)
#demean
def demean(X):
    return X - np.mean(X, axis=0)

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

3. 梯度上升法模拟
def J(w, X):
    return np.sum((X.dot(w) ** 2)) / len(X)

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

def dJ_debug(w, X_b, 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] = (J(w_1, X_b) - J(w_2, X_b)) / (2 * epsilon)
    return res

def gradient_ascent(dJ, X_b, initial_w, eta, n_iters=1e4, epsilon=1e-8):
    
    def direction(w):
        return w / np.linalg.norm(w)
    
    w = direction(initial_w)
    cur_iter = 0
    while cur_iter < n_iters:
        gradient = dJ( w, X_b)
        last_w = w
        w = w + eta * gradient # -改为+
        w = direction(w) # 每次求一个单位方向,即转换为单位向量
            
        if abs(J(w, X_b) - J(last_w, X_b)) < epsilon:
            break

        cur_iter += 1
    return w

initial_w = np.random.random(X.shape[1]) # w不能用0向量开始
initial_w

eta = 0.001
# 注意不能标准化数据,标准化后消掉了方差(方差为1)
w = gradient_ascent(dJ, X_demean, initial_w, eta)
w

w2 = gradient_ascent(dJ_debug, X_demean, initial_w, eta)
w2
#使用debug看看结果是否一致

plt.scatter(X_demean[:,0], X_demean[:,1])
plt.plot([w[0]*(-30), w[0]*30], [w[1]*(-30), w[1]*30], color='r') # w是方向向量,比较小,所以乘系数放大
plt.show()

三、以二维样本分析——求数据的前n个主成分

(一 )分析

问:求出第一个主成分以后,如何求出下一个主成分?

答:数据进行改变,将数据在第一个主成分上的分量去掉。在新的数据上求第一主成分。

如图:去掉第一主成分,实质上是将向量映射在与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()

#向量化
X3 = X - X.dot(w).reshape(-1,1) * w
plt.scatter(X3[:,0], X3[:,1])
plt.show()

(三)求前n个主成分

def first_n_components(dJ, n, X, eta=0.01, n_iters = 1e4, epsilon=1e-8):
    
    X_pca = X.copy()
    X_pce = demean(X_pca)
    res = []
    for i in range(n):
        initial_w = np.random.random(X_pca.shape[1])
        w = gradient_ascent(dJ, X_pca, initial_w, eta)
        res.append(w)
        
        X_pca = X_pca - X_pca.dot(w).reshape(-1,1) * w
    return res
    

ws = first_n_components(dJ, 2, X)
ws

ws[0].dot(ws[1]) # 2个w之积为0

四、高维数据向低维数据映射——PCA代码封装

(一 )分析

对于m × n的矩阵X,与k × n的矩阵W

$$ X = \left(\begin{matrix} X_1^{(1)} & X_2^{(1)} & X_3^{(1)} &\cdots& X_n^{(1)} \\ X_1^{(2)} & X_2^{(2)} & X_3^{(2)} &\cdots& X_n^{(2)} \\ X_1^{(3)} & X_2^{(3)} & X_3^{(3)} &\cdots& X_n^{(3)} \\ \cdots & \cdots & \cdots & \ddots & \cdots & \\ X_1^{(m)} & X_2^{(m)} & X_3^{(m)} &\cdots& X_n^{(m)} \\ \end{matrix}\right) $$

$$ W = \left(\begin{matrix} W_1^{(1)} & W_2^{(1)} & W_3^{(1)} &\cdots& W_n^{(1)} \\ W_1^{(2)} & W_2^{(2)} & W_3^{(2)} &\cdots& W_n^{(2)} \\ W_1^{(3)} & W_2^{(3)} & W_3^{(3)} &\cdots& W_n^{(3)} \\ \cdots & \cdots & \cdots & \ddots & \cdots & \\ W_1^{(k)} & W_2^{(k)} & W_3^{(k)} &\cdots& W_n^{(k)} \\ \end{matrix}\right) $$

高维降低维实质上是:$X \cdot W_k^{T} = X_k$(m×n · n×k = m×k),即n维降到了k维。

相反,恢复维度(k维升到n维),对于m × k的矩阵X,与k × n的矩阵W

$X_k \cdot W_k^ = X_m$(m×k · k×n = m×n)。

(二)编码

将完整的PCA代码封装成类,在文件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 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

(三)测试

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)

from playML.PCA import PCA

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

X_reduction = pca.transform(X)
X_restore = pca.inverse_transform(X_reduction)

plt.scatter(X[:,0], X[:,1], color="b", alpha=0.5 )
plt.scatter(X_restore[:,0], X_restore[:,1], color="r", alpha=0.5 )
plt.show()

(四)scikit-learn中的PCA

from sklearn.decomposition import PCA

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

pca.components_
X_reduction = pca.transform(X)
X_restore = pca.inverse_transform(X_reduction)

plt.scatter(X[:,0], X[:,1], color="b", alpha=0.5 )
plt.scatter(X_restore[:,0], X_restore[:,1], color="r", alpha=0.5 )
plt.show()

五、真实数据模拟PCA降维

(一 )手写识别数据

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

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

#数据集分离
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
X_train.shape


#使用kNN算法进行分类
from sklearn.neighbors import KNeighborsClassifier
%%time
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
knn_clf.score(X_test, y_test)

# 使用PCA
from sklearn.decomposition import PCA

## 降到2维
pca = PCA(n_components=2)
pca.fit(X_train)

X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)

%%time
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)

knn_clf.score(X_test_reduction, y_test)

##降到2维以后效果比较差,
##调用pca提供的函数,解释每一个指标能够解释的方差有多少

pca.explained_variance_ratio_

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])]
)
plt.show()

##构造函数传 精度 参数,这个精度即上面的方差覆盖率
pca = PCA(0.95)
pca.fit(X_train)

pca.n_components_

X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)

%%time
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)

knn_clf.score(X_test_reduction, y_test)

#降到两维进行可视化
pca = PCA(n_components=2)
pca.fit(X_train)
X_reduction = pca.transform(X)

X_reduction.shape
y.shape

for i in range(10):
    plt.scatter(X_reduction[y==i,0], X_reduction[y==i,1], alpha=0.8)
plt.show()

(二)更大更正规更经典的手写识别数据——MNIST

MNIST数据集,是由美国人口普查局的高中学生和雇员手写的一组70,000个小数字图像。每个图像都标有它所代表的数字。这个集已经被研究了很多,它通常被称为机器学习的“Hello World”

1. 加载数据集

方法一:在程序中调用代码直接下载。

import numpy as np
from sklearn.datasets import fetch_mldata
#下载MNIST数据集,耐心等待一段时间
mnist = fetch_mldata("MNIST original")

方法二:如果上述下载失败,使用此方法。

(1)下载数据集文件,地址:mnist-original文件下载

(2)文件解压出来获得mnist-original.mat文件,放入自己指定的目录,我把他和python文件放一起了。

(3)在程序中执行代码加载即可,其中loadmat函数参数里写mnist-original.mat文件的地址,因为我放在了同一个目录下,所以用./即可。

from scipy.io import loadmat
mnist = loadmat("./mnist-original.mat")
mnist_data = mnist["data"].T
mnist_label = mnist["label"][0]
2. 编码实现
# MNIST
from scipy.io import loadmat
mnist = loadmat("./mnist-original.mat")
mnist_data = mnist["data"].T
mnist_label = mnist["label"][0]

X , y = mnist_data, mnist_label

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)

X_train.shape
y_test.shape

# 使用kNN
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_train, y_train)

%time knn_clf.score(X_test, y_test)

# PCA降维
from sklearn.decomposition import PCA
pca = PCA(0.9)
pca.fit(X_train)

X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)

## 降维后的维度
X_train_reduction.shape

knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_train_reduction, y_train)

%time knn_clf.score(X_test_reduction, y_test)

1564571545390

观察可以发现使用PCA降维本身丢失了一部分数据,为什么准确率会比kNN的准确还要高呢?

答:因为PCA具有降噪的特性。

六、降噪示例

(一 )手写识别的例子

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

#准备数据
digits = datasets.load_digits()
X = digits.data
y = digits.target
#添加噪音
noisy_digits = X + np.random.normal(0,4,size=X.shape)
#每个数字选出前10行样本
example_digits = noisy_digits[y==0,:][:10]
for num in range(1,10):
    X_num = noisy_digits[y==num,:][:10]
    example_digits = np.vstack([example_digits, X_num])

example_digits.shape    

#绘图函数
def plot_digits(data):
    fig, axes = plt.subplots(10,10, figsize=(10, 10),
                             subplot_kw={'xticks':[],'yticks':[]},
                             gridspec_kw=dict(hspace=0.1,wspace=0.1)) 
    for i, ax in enumerate(axes.flat):
        ax.imshow(data[i].reshape(8, 8),
                 cmap='binary', interpolation='nearest',
                 clim=(0,16))
        
    plt.show()
    
plot_digits(example_digits)

#PCA降噪处理
from sklearn.decomposition import PCA
# 精确度为50%
pca = PCA(0.5)
pca.fit(noisy_digits)
#降维后的维度
pca.n_components_
#降维、再恢复维度,降维过程中丢失了信息也就是降噪了
components = pca.transform(example_digits)
filtered_digits = pca.inverse_transform(components)
plot_digits(filtered_digits)

(二)人脸识别之特征脸

1. 下载人脸识别数据集要注意。

方法一:直接执行faces = fetch_lfw_people()

方法二:本地下载:链接:https://pan.baidu.com/s/1QD3IFl4IzK3CCa_si6mGeQ 提取码:r3tk

把下载好的.tgz格式文件直接放入C:\Users\xxx\scikit_learn_data\lfw_home,其中xxx是用户名,再次执行faces = fetch_lfw_people(),就好了。如果还是无法加载,Restart & Run All一下即可。

2. 下载/加载数据
import numpy as np
import matplotlib.pyplot as plt

#下载/加载人脸识别数据库
from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people()

# faces数据相关属性,数据的维度
faces.keys()
faces.data.shape
faces.images.shape

# 对数据进行随机化处理(这里是随机化下标)
random_indexes = np.random.permutation(len(faces.data))
X = faces.data[random_indexes]
# 随机取36个样本
example_faces = X[:36,:]
example_faces.shape

# 画图函数
def plot_faces(faces):
    fig, axes = plt.subplots(6, 6, figsize=(10, 10),
                            subplot_kw={'xticks':[],'yticks':[]},
                            gridspec_kw=dict(hspace=0.1,wspace=0.1))
    for i, ax in enumerate(axes.flat):
        ax.imshow(faces[i].reshape(62, 47), cmap='bone')
        
    plt.show()
 
# 绘图
plot_faces(example_faces)

3. 提取特征脸——PCA降噪处理
%%time
from sklearn.decomposition import PCA
pca = PCA(svd_solver='randomized')
pca.fit(X)

# 查看降维后的维度
pca.components_.shape

# 绘图
plot_faces(pca.components_[:36,:])

每一个人脸都是这些特征脸的线性组合,特征脸按照重要程度顺次的排列,越是排在前面的图像就越笼统,

Last modification:February 7th, 2020 at 12:50 pm
如果觉得我的文章对你有用,请随意赞赏~