sklearn与机器学习系列专题之降维（一）一文弄懂PCA特征筛选&降维

伴随着通信与互联网技术的不断发展，人们收集和获取数据的能力越来越强，而这些数据已呈现出维数高、规模大和结构复杂等特点，当数据量非常大时，会面临维度灾难，即:1.在高维情况下，数据样本稀疏；2.涉及距离、内积的计算变得困难。缓解灾难的一个重要途径就是降维。

本篇博客主要以PCA为例，讲解sklearn实现降维。

1.标准PCA

如下图所示，我们用两个特征——一个是attack、一个是defense来描述样本。作两条垂直的参考线，方差相差很大。在方差较大的方向上，所包含的信息就较多。降维要尽可能保持原始数据的信息，故要取使方差最大的一个维度（本例中从二维降到一维）。

主成分分析（Principal Component Analysis，PCA）是一种最为常见的无监督降维方法，通过降维技术把多个变量化为少数几个主成分，这些主成分能够反映原始变量的大部分信息。

PCA降维有两个主要准则:
1.最近重构性——样本重构后的点距离原来的点误差之和最小；
2.最大可分性——样本点在低维空间的投影尽可能分开。

以鸢尾花（iris）数据集为例，采用PCA降维。

#加载matplotlib用于数据的可视化
import matplotlib.pyplot as plt   
#加载PCA算法包        
from sklearn.decomposition import PCA           
from sklearn.datasets import load_iris

data=load_iris()
print("鸢尾花特征：",data.feature_names)
y=data.target
x=data.data
x.shape,y.shape

#加载PCA算法，设置降维后主成分数目为2
pca=PCA(n_components=2)   

#对样本进行降维
reduced_x=pca.fit_transform(x)
reduced_x.shape

red_x,red_y=[],[]
blue_x,blue_y=[],[]
green_x,green_y=[],[]

for i in range(len(reduced_x)):
    if y[i] ==0:
        red_x.append(reduced_x[i][0])
        red_y.append(reduced_x[i][1])
    elif y[i]==1:
        blue_x.append(reduced_x[i][0])
        blue_y.append(reduced_x[i][1])
    else:
        green_x.append(reduced_x[i][0])
        green_y.append(reduced_x[i][1])
#可视化
plt.scatter(red_x,red_y,c='r',marker='x')
plt.scatter(blue_x,blue_y,c='b',marker='D')
plt.scatter(green_x,green_y,c='g',marker='.')
plt.show()

运行程序，可知，特征在降维前的维度为(150,4)，特征分别为’sepal length (cm)’, ‘sepal width (cm)’, ‘petal length (cm)’, ‘petal width (cm)’，降维后的维度为(150,2)。

降维后，得到如下图所示的效果。

2.非负矩阵分解

当数据集由非负元素组成时，可以使用非负矩阵分解（Non-Negative Matrix Factorization，NNMF）代替标准PCA，该方法通过添加正则化项提高精度。非负矩阵分解法及其它因式分解法对于推荐系统和主体建模等更先进的技术十分实用。矩阵的因式分解允许将其表达为共享一个实体的乘积。值得一提的是，非负矩阵分解对其参数（特别是初始化和正则化）非常敏感。

依然以鸢尾花为例，示意sklearn实现NNMF。初始化参数可以选择不同的值，它们决定了数据矩阵最初处理的方式，并随机选择用于缩放的矩阵。

#加载matplotlib用于数据的可视化
import matplotlib.pyplot as plt     
#加载PCA算法包     
from sklearn.decomposition import NMF           
from sklearn.datasets import load_iris

data=load_iris()
y=data.target
x=data.data
x.shape

#加载PCA算法，设置降维后主成分数目为2
nmf=NMF(n_components=2,init='random',l1_ratio=0.1)     
#对样本进行降维
reduced_x=nmf.fit_transform(x)
reduced_x.shape

red_x,red_y=[],[]
blue_x,blue_y=[],[]
green_x,green_y=[],[]

for i in range(len(reduced_x)):
    if y[i] ==0:
        red_x.append(reduced_x[i][0])
        red_y.append(reduced_x[i][1])
    elif y[i]==1:
        blue_x.append(reduced_x[i][0])
        blue_y.append(reduced_x[i][1])
    else:
        green_x.append(reduced_x[i][0])
        green_y.append(reduced_x[i][1])
#可视化
plt.scatter(red_x,red_y,c='r',marker='x')
plt.scatter(blue_x,blue_y,c='b',marker='D')
plt.scatter(green_x,green_y,c='g',marker='.')
plt.show()

3.稀疏PCA

sklearn提供了可以解决特定问题的不同PCA变体。当使用稀疏PCA（Sparse PCA）时，可以在提取主元的同时利用数据的自然稀疏性降维。例如，在对手写数字分类时，数据的初始维度可能非常高，而应用标准PCA仅会选择平均最重要的特征。当矩阵中的非零元素非常少时，可以利用稀疏PCA，重建每个元素，在大多数情况下，这些特定的主元将始终是最重要的，可以弥补标准PCA丢弃部分主元的缺点。

以下代码显示了具有60个主元的稀疏PCA，可以通过L1范数正则化控制稀疏度，正则化参数越大，降维结果越稀疏。

import numpy as np
import matplotlib.pyplot as plt

from sklearn.datasets import load_digits
from sklearn.decomposition import SparsePCA

# 加载数据集并打印数据形状
digits = load_digits()
digits.data.shape

# 随机显示部分图片
selection = np.random.randint(0, 1797, size=100)

fig, ax = plt.subplots(10, 10, figsize=(10, 10))

samples = [digits.data[x].reshape((8, 8)) for x in selection]

for i in range(10):
    for j in range(10):
        ax[i, j].set_axis_off()
        ax[i, j].imshow(samples[(i * 8) + j], cmap='gray')

plt.show()

# 对数据集降维
spca = SparsePCA(n_components=36, alpha=0.1)
X_spca = spca.fit_transform(digits.data / 255)

print('SPCA components shape:',spca.components_.shape)

#展示降维后的效果
fig, ax = plt.subplots(10, 10, figsize=(10, 10))
samples = [X_spca[x].reshape((6, 6)) for x in selection]
for i in range(10):
    for j in range(10):
        ax[i, j].set_axis_off()
        ax[i, j].imshow(samples[(i * 6) + j], cmap='gray')

plt.show()

降维前，原数据集的维度为(1797,64)，输出的示意图如下：

利用稀疏矩阵拟合后的维度为(1797,36)，输出的示意图如下：

4.核PCA

标准PCA方法假设从高维空间到低维空间的函数映射是线性的，但在很多任务中，可能需要非线性映射才能找到合适的降维空间来降维。非线性降维的一种常用方法是基于核技巧对线性降维方法进行核化，核PCA（kernelized PCA）是对PCA的一种推广。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist,squareform
from scipy import exp
from scipy.linalg import eigh
from sklearn.datasets import make_moons
from sklearn.datasets import make_circles
from sklearn.decomposition import PCA
from matplotlib.ticker import FormatStrFormatter

def rbf_kernel_pca(X,gama,n_components):

    #1：计算样本对欧几里得距离，并生成核矩阵k(x,y)=exp(-gama *||x-y||^2)，
#x和y表示样本，构建一个NXN的核矩阵，矩阵值是样本间的欧氏距离值。
    #计算两两样本间欧几里得距离
    sq_dists = pdist (X, 'sqeuclidean')     
    ##距离平方
    mat_sq_dists=squareform(sq_dists)     
    ##计算对称核矩阵
    K=exp(-gama * mat_sq_dists)     
    
    #2:聚集核矩阵K'=K-L*K-K*L + L*K*L，其中L是一个nXn的矩阵(和核矩阵K
#的维数相同，所有的值都是1/n。聚集核矩阵的必要性是：样本经过标准化处理
#后，当在生成协方差矩阵并以非线性特征的组合替代点积时，所有特征的均值为
#0；但用低维点积计算时并没有精确计算新的高维特征空间，也无法确定新特征
#空间的中心在零点。
    N=K.shape[0]
    one_n = np.ones((N,N))/N #NXN单位矩阵
    K=K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n) 
       
    #3：对聚集后的核矩阵求取特征值和特征向量 
    eigvals,eigvecs = eigh(K)
    
    #4：选择前k个特征值所对应的特征向量，和PCA不同，KPCA得到的K个特
#征，不是主成分轴，而是高维映射到低维后的低维特征数量核化过程是低维映射
#到高维，pca是降维，经过核化后的维度已经不是原来的特征空间。核化是低维
#映射到高维，但并不是在高维空间计算(非线性特征组合)而是在低维空间计算(点
#积)，做到这点关键是核函数，核函数通过两个向量点积来度量向量间相似度，
#能在低维空间内近似计算出高维空间的非线性特征空间。
    X_pc = np.column_stack((eigvecs[:,-i] for i in range(1,n_components+1))) 
    return X_pc 

###分离半月形数据
##生成二维线性不可分数据
X,y=make_moons(n_samples=100,random_state=123)
plt.scatter(X[y==0,0],X[y==0,1],color='red',marker='^',alpha=0.5)
plt.scatter(X[y==1,0],X[y==1,1],color='blue',marker='o',alpha=0.5)
plt.show()

##PCA降维，映射到主成分，仍不能很好线性分类
sk_pca = PCA(n_components=2)
X_spca=sk_pca.fit_transform(X)
fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(7,3))
ax[0].scatter(X_spca[y==0,0],X_spca[y==0,1],color='red',marker='^',alpha=0.5)
ax[0].scatter(X_spca[y==1,0],X_spca[y==1,1],color='blue',marker='o',alpha=0.5)
ax[1].scatter(X_spca[y==0,0],np.zeros((50,1))+0.02,color='red',marker='^',alpha=0.5)
ax[1].scatter(X_spca[y==1,0],np.zeros((50,1))-0.02,color='blue',marker='^',alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1,1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()

##利用基于RBF核的KPCA来实现线性可分
X_kpca=rbf_kernel_pca(X, gama=15, n_components=2)
fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(7,3))
ax[0].scatter(X_kpca[y==0,0],X_kpca[y==0,1],color='red',marker='^',alpha=0.5)
ax[0].scatter(X_kpca[y==1,0],X_kpca[y==1,1],color='blue',marker='o',alpha=0.5)
ax[1].scatter(X_kpca[y==0,0],np.zeros((50,1))+0.02,color='red',marker='^',alpha=0.5)
ax[1].scatter(X_kpca[y==1,0],np.zeros((50,1))-0.02,color='blue',marker='^',alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1,1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
ax[0].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))
plt.show()

所生成的初始图像为：

用标准PCA降维时，仍然不能很好地线性可分。

使用核PCA降维，具有很好的区分度。

下篇预告

sklearn与机器学习系列专题之降维（二）一文弄懂LDA特征筛选&降维。
敬请期待！

欢迎关注公众号“码点联盟”，不定期分享机器学习、深度学习干货，及各种资料大礼包！