通过降维压缩数据

在前面已经介绍了几种不同的特征选择技术对数据集进行降维的方法。另一种常用于降维的特征选择方法就是特征抽取。数据压缩也是机器学习领域中的一个重要内容。数据压缩技术可以帮助我们对数据及逆行存储和分析。

这里主要涵盖三个主题：

1.无监督数据压缩——主成分分析（PCA）

2.基于类别可分最大化的监督降维技术——线性判别分析（LDA）

3.通过核主成分分析进行非线性降维

无监督数据降维技术——主成分分析

特征抽取可以理解为：在尽可能多地保持相关信息的情况下，对数据进行压缩的一种方法。特征抽取通常用于提高计算效率，同样也可以帮助我们降低“维度灾难”——尤其当模型不适于正则化处理时。

主成分分析（PCA）是一种广泛应用于不同领域的无监督线性数据转换技术，其突出作用是降维。简而言之，PCA的目标是在高维数据中找到最大最大方差的方向，并将数据映射到一个维度不大于原始数据的新的子空间上。

如果使用PCA降维，我们将构建一个d×k维的转换矩阵W，这样就可以将一个样本向量x映射到一个新的k维特征子空间上去，此空间的维度小于原始的d维特征空间：

X=[x1,x2,….,xd]  

Z=[z1,z2,….,zk]

完成从原始d维数据到新的k维子空间（一般情况下k<<d）的转换后，第一主成分的方差应该是最大的，由于各主成分之间是不相关的，后续各主成分也具备尽可能大的方差。需注意的是，主成分的方向对数据值的范围高度敏感，如果特征的值使用不同的度量标准，我们需要先对特征进行标准化处理，以让各特征具有相同的重要性。

在详细讨论使用PCA算法降维之前，我们先通过以下几个步骤来概括一下算法的流程：

1）对原始d维数据集做标准化处理

2）构造样本的协方差矩阵

3）计算协方差矩阵的特征值和特征向量

4）选择与前k个最大特征值对应的特征向量，其中k为新特征空间的维度

5）通过前k个特征向量构建映射矩阵W

6)通过映射矩阵W将d维的输入数据集X转换到新的k维特征子空间

总体方差与贡献方差

主成分分析算法的前四个步骤：数据标准化、构造协方差矩阵、获得协方差矩阵的特征值和特征向量，以及按降序排列特征值所对应的特征向量：

import pandas as pd
df_wine=pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data',header=None)
from sklearn.cross_validation import train_test_split
from sklearn.preprocessing import StandardScaler
X,y=df_wine.iloc[:,1:].values,df_wine.iloc[:,0].values
X_train,X_test,y_train,y_test=\
                train_test_split(X,y,test_size=0.3,random_state=0)
sc=StandardScaler()
X_train_std=sc.fit_transform(X_train)
X_test_std=sc.fit_transform(X_test)

通过上述代码完成数据的预处理后，进入第二步：构造协方差矩阵。此d×d维协方差矩阵是沿主对角线对称的，其中d为数据集的维度，此矩阵成对地存储了不同特征之间的协方差。

import numpy as np
cov_mat=np.cov(X_train_std.T)
eigen_vals,eigen_vecs=np.linalg.eig(cov_mat)
print('\n Eigenvalues \n%s'% eigen_vals)

应用numpy.cov函数，我们计算得到了经标准化处理的训练数据集的协方差矩阵。使用linalg.eig函数，通过特征分解。

因为要将数据集压缩到一个新的特征子空间上来实现数据降维，所以我们只选择那些包含最多信息（方差最大）的特征向量（主成分）组成子集。由于特征值的大小决定了特征向量的重要性，因此需要将特征值按降序排列，我们感兴趣的是排序在前k个的特征值所对应的特征向量。在整理包含信息量最大的前k个特征向量前，我们先绘制特征值的方差贡献率图像。

特征值λj的方差贡献率是指，特征值λj与所有特征值和的比值

使用Numpy的cumsum函数，我们可以计算出累计方差，其图像可通过matplotlib的step函数绘制

tot=sum(eigen_vals)
var_exp=[(i/tot) for i in sorted(eigen_vals,reverse=True)]
cum_var_exp=np.cumsum(var_exp)
import matplotlib.pyplot as plt
plt.bar(range(1,14),var_exp,alpha=0.5,align='center',label='individual explained variance')
plt.step(range(1,14),cum_var_exp,where='mid',label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.show()

特征转换

在将协方差矩阵分解为特征对后，我们继续执行PCA方法的最后三个步骤。

我们将对特征值按降序进行排列，并通过挑选出对应的特征向量构造出映射矩阵，然后使用映射矩阵将数据转换到低维的子空间上。

首先，按特征值的降序排列特征对：

eigen_pairs=[(np.abs(eigen_vals[i]),eigen_vecs[:,i])for i in range(len(eigen_vals))]
eigen_pairs.sort(reverse=True)

使用scikit-learn进行主成分分析

from matplotlib.colors import ListedColormap
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
df_wine=pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data',header=None)
from sklearn.cross_validation import train_test_split
from sklearn.preprocessing import StandardScaler
X,y=df_wine.iloc[:,1:].values,df_wine.iloc[:,0].values
X_train,X_test,y_train,y_test=\
                train_test_split(X,y,test_size=0.3,random_state=0)
sc=StandardScaler()
X_train_std=sc.fit_transform(X_train)
X_test_std=sc.fit_transform(X_test)
def plot_decision_regions(X,y,classifier,resolution=0.02):
    markers=('s','x','o','^','v')
    colors=('red','blue','lightgreen','gray','cyan')
    cmap=ListedColormap(colors[:len(np.unique(y))])

    x1_min,x1_max=X[:,0].min()-1,X[:.0].max()+1
    x2_min, x2_max = X[:, 1].min() - 1, X[:.1].max() + 1
    xx1,xx2=np.meshgrid(np.arange(x1_min,x1_max,resolution),np.arange(x2_min,x2_max,resolution))
    Z=classifier.predict(np.array([xx1.ravel(),xx2.ravel()]).T)
    Z=Z.reshape(xx1.shape)
    plt.contourf(xx1,xx1,Z,alpha=0.4,cmap=cmap)
    plt.xlim(xx1.min(),xx1.max())
    plt.ylim(xx2.min(),xx2.max())

    for idx,cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y==cl,0],y=X[y==cl,1],alpha=0.8,c=cmap(idx),marker=markers[idx],label=cl)

from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
pca=PCA(n_components=2)
lr=LogisticRegression()
X_train_pca=pca.fit_transform(X_train_std)
X_test_pca=pca.transform(X_test_std)
lr.fit(X_train_pca,y_train)
plot_decision_regions(X_train_pca,y_train,classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()

通过线性判别分析压缩无监督数据

线性判别分析（LDA）是一种可作为特征抽取的技术，它可以提高数据分析过程中的计算效率，同时，对于不适用于正则化的模型，它可以降低因维度灾难带来的过拟合。

LDA的基本概念与PCA非常相似，PCA试图在数据集中找到方差最大的正交的主成分分量的轴，而LDA的目标是发现可以最优化分类的特征子空间。LDA与PCA都是可用于降低数据集维度的线性转换技巧。其中，PCA是无监督算法，而LDA是监督算法。因此我们可以这样直观地认为：与PCA相比，LDA是一种更优越的用于分类的特征提取技术。

总结一下LDA方法的关键步骤：

1）对d维数据集进行标准化处理（d为特征的数量）

2）对于每一类别，计算d维的均值向量

3）构造类间的散布矩阵S_B以及类内的散布矩阵S_W

4)计算矩阵S_w^-1 S_B的特征值及对应的特征向量

5）选取前k个特征值所对应的特征向量，构造一个d×k维的转换矩阵W,其中特征向量以列的形式排列

6）使用转换矩阵W将样本映射到新的特征子空间上

计算散布矩阵

通过均值向量，计算类内散布矩阵S_w：

import numpy as np
np.set_printoptions(precision=4)
mean_vecs=[]
for label in range(1,4):
    mean_vecs.append(np.mean(X_train_std[y_train==label],axis=0))
    print('MV %s:%s \n'%(label,mean_vecs[label-1]))

d=13#特征数
S_W=np.zeros((d,d))
for label,mv in zip(range(1,4),mean_vecs):
    class_scatter=np.zeros((d,d))
    for row  in X[y==label]:
        row,mv=row.reshape(d,1),mv.reshape(d,1)
        class_scatter+=(row-mv).dot((row-mv).T)
    S_W+=class_scatter
print('Within-class scatter matrix:%sx%s'%(S_W.shape[0],S_W.shape[1]))

计算类间散布矩阵S_B:

mean_overall=np.mean(X_train_std,axis=0)
d=13
S_B=np.zeros((d,d))
for i,mean_vec in enumerate(mean_vecs):
    n=X[y==i+1,:].shape[0]
    mean_vec=mean_vec.reshape(d,1)
    mean_overall=mean_overall.reshape(d,1)
    S_B+=n*(mean_vec-mean_overall).dot((mean_vec-mean_overall).T)
print('Between-class scatter matrix:%sx%s'%(S_B.shape[0],S_B.shape[1]))

在新特征子空间上选取线性判别算法

LDA余下的步骤与PCA的步骤相似。不过，这里我们不对协方差矩阵做特征分解，而是求解矩阵S_w^-1S_B  的广义特征值

eigen_vals,eigen_vecs=np.linalg.eig(np.linalg.inv(S_W).dot(S_B))

在求得了特征对之后，我们按照降序对特征值进行排序：

eigen_pairs=[(np.abs(eigen_vals[i]),eigen_vecs[:,i]) for i in range(len(eigen_vals))]
eigen_pairs=sorted(eigen_pairs,key=lambda k: k[0],reverse=True)

d×d维协方差矩阵的秩最大为d-1，而且确实可以发现，我们只得到了两个非零特征值（实际得到的第3~13个特征值并非完全为0，而是趋近于0的实数）请注意，在极少的情况下可达到完美的共线性（所有样本的点落在一条直线上），这时协方差矩阵的秩为1，将导致矩阵只有一个含非零特征值的特征向量。

为了度量线性判别（特征向量）可以获取多少可区分类别的信息，与前面PCA中对累计方差的绘制类似，我们按照特征值降序绘制出特征对线性判别信息保持程度的图像。

tot=sum(eigen_vals.real)
discr=[(i/tot)for i in sorted(eigen_vals.real,reverse=True)]
cum_discr=np.cumsum(discr)
plt.bar(range(1,14),discr,alpha=0.5,align='center',label='individual "discriminability"')
plt.step(range(1,14),cum_discr,where='mid',label='cumulative "discriminability"')
plt.xlabel('Linear Discriminants')
plt.ylabel('"discriminability" ratio')
plt.ylim([-0.1,1.1])
plt.legend(loc='best')
plt.show()

下面将这两个判别能力最强的特征向量列来构建转换矩阵W:

w=np.hstack((eigen_pairs[0][1][:,np.newaxis].real,eigen_pairs[1][1][:,np.newaxis].real))
print("Matrix W:\n",w)

使用scikit-learn进行LDA分析

from sklearn.lda import  LDA
lda=LDA(n_components=2)
X_train_lda=lda.fit_transform(X_train_std,y_train)
lr=LogisticRegression()
lr=lr.fit(X_train_lda,y_train)
plot_decision_regions(X_train_lda,y_train,classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.show()


x_test_lda=lda.transform(X_test_std)
plot_decision_regions(x_test_lda,y_test,classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.show()

使用核主成分分析进行非线性映射

使用核PCA，我们将非线性可分的数据转换到一个适合对其进行线性分类的新的低维子空间上。

核函数与核技巧

通过将非线性可分问题映射到维度更高的特征空间，使其在新的特征空间上线性可分。为了将样本x转换到维度更高的k维子空间，我们定义如下非线性映射函数使R^d——>R^k(k>>d)

换句话说，利用核PCA，我们可以通过非线性映射将数据转换到一个高维空间，然后在此高维空间中使用标准PCA，将其映射到另外一个低维空间中。这种方法的缺点是会带来高昂的计算成本。