降维实践（PCA，LDA）

17 降维简介

当特征选择完成后，可以直接训练模型了，但是可能由于特征矩阵过大，导致计算量大，训练时间长的问题，因此降低特征矩阵维度也是必不可少的。常见的降维方法除了以上提到的基于L1惩罚项的模型以外，另外还有主成分分析法（PCA）和线性判别分析（LDA），线性判别分析本身也是一个分类模型。PCA和LDA有很多的相似点，其本质是要将原始的样本映射到维度更低的样本空间中，但是PCA和LDA的映射目标不一样：PCA是为了让映射后的样本具有最大的发散性；而LDA是为了让映射后的样本有最好的分类性能。所以说PCA是一种无监督的降维方法，而LDA是一种有监督的降维方法。
PCA、LDA降维一般假设数据集为线性可分，如果用这两种方法，对线性不可分的数据集进行降维，效果往往不理想。本质上PCA和LDA还是一种线性变换。而线性不可分数据应该是很普遍的，对线性不可分数据集该如何进行降维呢？这里我们介绍一种核PCA方法，这样降维方法综合了核技巧及PCA思想，对非线性数据集降维有非常好的效果。
此外，这里我们还介绍SVD方法，这也是一种非常有效的降维方法。

17.1 PCA简介

主成分分析（Principal Components Analysis），简称PCA，是一种数据降维技术，用于数据预处理。一般我们获取的原始数据维度都很高，比如1000个特征，在这1000个特征中可能包含了很多无用的信息或者噪声，真正有用的特征才50个或更少，那么我们可以运用PCA算法将1000个特征降到50个特征。这样不仅可以去除无用的噪声，还能减少很大的计算量。
PCA算法是如何实现的？
简单来说，就是将数据从原特征空间转换到新的特征空间中，例如原始的空间是三维的(x,y,z)，x、y、z分别是原始空间的三个基，我们可以通过某种方法，用新的坐标系(a,b,c)来表示原始的数据，那么a、b、c就是新的基，它们组成新的特征空间。在新的特征空间中，可能所有的数据在c上的投影都接近于0，即可以忽略，那么我们就可以直接用(a,b)来表示数据，这样数据就从三维的(x,y,z)降到了二维的(a,b)。
问题是如何求新的基(a,b,c)?
一般步骤是这样的：
1）对原始数据集做标准化处理。
2）求协方差矩阵。
3）计算协方差矩阵的特征值和特征向量。
4）选择前k个最大的特征向量，k小于原数据集维度。
5）通过前k个特征向量组成了新的特征空间，设为W。
6)通过矩阵W,把原数据转换到新的k维特征子空间。

17.2 PCA算法实现

这里以葡萄酒数据为例，数据集特征如下：

数据来源于：https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data
1）对原数据集做标准化处理
导入需要的库及数据

   

1

2

3

4

5

6

7

8

9

10

import pandas as pd

df_wine = pd.read_csv(‘https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data‘, header=None)

df_wine.columns = [‘Class label‘, ‘Alcohol‘, ‘Malic acid‘, ‘Ash‘,

‘Alcalinity of ash‘, ‘Magnesium‘, ‘Total phenols‘,

‘Flavanoids‘, ‘Nonflavanoid phenols‘, ‘Proanthocyanins‘,

‘Color intensity‘, ‘Hue‘, ‘OD280/OD315 of diluted wines‘, ‘Proline‘]

df_wine.head()

部分内容：

为便于后续处理，把数据集分为训练集和测试集，划分比例为7:3

1

2

3

4

5

from sklearn.cross_validation import train_test_split

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)

对原数据进行标准化处理

1

2

3

4

5

from sklearn.preprocessing import StandardScaler

sc = StandardScaler()

X_train_std = sc.fit_transform(X_train)

X_test_std = sc.fit_transform(X_test)

2) 求协方差矩阵
这里使用numpy.cov函数，求标准化后数据的协方差矩阵

1

2

import numpy as np

cov_mat = np.cov(X_train_std.T)

3）计算协方差矩阵的特征值和特征向量
使用np.linalg.eig函数，求协方差的特征值和特征向量

1

2

eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)

print(‘\nEigenvalues \n%s‘ % eigen_vals)

得到13个特征向量：
Eigenvalues
[ 4.8923083 2.46635032 1.42809973 1.01233462 0.84906459 0.60181514 0.52251546 0.08414846 0.33051429 0.29595018 0.16831254 0.21432212 0.2399553 ]
要实现降维，我们可以选择前k个最多信息（或方差最大）特征向量组成新的子集，由于特征值的大小决定了特征向量的重要性，因此，可以通过对特征值的排序，获取前k个特征值。特征值λ_i的方差贡献率是指特征值λ_i与所有特征值和的比例：

我们可以通过numpy.cumsum函数计算累计方差。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

tot = sum(eigen_vals)

var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]

cum_var_exp = np.cumsum(var_exp)

#然后用matplotlib各主成分的方差贡献率图形。

import matplotlib.pyplot as plt

%matplotlib inline

import matplotlib.font_manager as fm

myfont = fm.FontProperties(fname=‘/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf‘)

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(‘方差贡献率‘,fontproperties=myfont,size=12)

plt.xlabel(‘主成分‘,fontproperties=myfont,size=12)

plt.legend(loc=‘best‘)

plt.tight_layout()

plt.show()

从这个图可以看出第一个主成分占了方差总和的40%左右，前两个主成分占了近60%。
4）选择前k个最大的特征向量，k小于原数据集维度
首先，按特征值按降序排序

1

2

3

4

5

# 构成一个元组 (eigenvalue, eigenvector)

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

# Sort the (eigenvalue, eigenvector) tuples from high to low

eigen_pairs.sort(reverse=True)

5）通过前k个特征向量组成了新的特征空间，设为W。
为便于数据可视化，这里我们取k=2，实际上前2个特征值已占了总方差的近60%。

1

2

3

w = np.hstack((eigen_pairs[0][1][:, np.newaxis],

eigen_pairs[1][1][:, np.newaxis]))

print(‘Matrix W:\n‘, w)

这样我们就可得到一个由这两个特征向量构成的13*2矩阵W:
Matrix W:
[[ 0.14669811 0.50417079]
[-0.24224554 0.24216889]
[-0.02993442 0.28698484]
[-0.25519002 -0.06468718]
[ 0.12079772 0.22995385]
[ 0.38934455 0.09363991]
[ 0.42326486 0.01088622]
[-0.30634956 0.01870216]
[ 0.30572219 0.03040352]
[-0.09869191 0.54527081]
[ 0.30032535 -0.27924322]
[ 0.36821154 -0.174365 ]
[ 0.29259713 0.36315461]]

6）通过矩阵W,把原数据转换到新的k维特征子空间
通过这个特征矩阵W,把原样本x转换到PCA的子空间上，得到一个新样本x^，。
x^，=xW
训练集与W点积后，把这个训练集转换到包括两个主成分的子空间上。然后，把子空间的数据可视化。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

X_train_pca = X_train_std.dot(w)

colors = [‘r‘, ‘b‘, ‘g‘]

markers = [‘s‘, ‘x‘, ‘o‘]

for l, c, m in zip(np.unique(y_train), colors, markers):

plt.scatter(X_train_pca[y_train==l, 0],

X_train_pca[y_train==l, 1],

c=c, label=l, marker=m)

plt.xlabel(‘PC 1‘)

plt.ylabel(‘PC 2‘)

plt.legend(loc=‘lower left‘)

plt.tight_layout()

plt.show()

从以上图形可以看出，大部分数据沿PC1方向分布，而且可以线性划分，在可视化图形时，为便于标识点，这里采用了y_train标签信息。

我们用来6步来实现PCA，这个过程还是比较麻烦的，是否有更简单的方法呢？
有的，接下来我们介绍利用Scikit-learn中PCA类进行降维。

17.3 利用Scikit-learn进行主成分分析

我们将使用Scikit-learn中PCA对数据集进行预测处理，然后使用逻辑斯谛回归对转换后的数据进行分类，最后对数据进行可视化。
1）数据预处理

1

2

3

4

5

from sklearn.decomposition import PCA

pca = PCA()

X_train_pca = pca.fit_transform(X_train_std)

pca.explained_variance_ratio_

得到主成分数据：
array([ 0.37329648, 0.18818926, 0.10896791, 0.07724389, 0.06478595, 0.04592014, 0.03986936, 0.02521914, 0.02258181, 0.01830924, 0.01635336, 0.01284271, 0.00642076])

2）可视化主成分方差贡献率图

1

2

3

4

5

plt.bar(range(1, 14), pca.explained_variance_ratio_, alpha=0.5, align=‘center‘)

plt.step(range(1, 14), np.cumsum(pca.explained_variance_ratio_), where=‘mid‘)

plt.ylabel(‘Explained variance ratio‘)

plt.xlabel(‘Principal components‘)

plt.show()

3）获取前2个主成分

1

2

3

4

5

6

7

8

pca = PCA(n_components=2)

X_train_pca = pca.fit_transform(X_train_std)

X_test_pca = pca.transform(X_test_std)

4）把训练集映射到主成分空间上，并可视化。

plt.scatter(X_train_pca[:,0], X_train_pca[:,1])

plt.xlabel(‘PC 1‘)

plt.ylabel(‘PC 2‘)

plt.show()

5）利用回归模型对数据进行分类。

1

2

3

4

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression()

lr = lr.fit(X_train_pca, y_train)

6）为了更好看到分类后情况，这里我们定义一个函数plot_decision_regions，通过这个函数对决策区域数据可视化。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):

# setup marker generator and color map

markers = (‘s‘, ‘x‘, ‘o‘, ‘^‘, ‘v‘)

colors = (‘red‘, ‘blue‘, ‘lightgreen‘, ‘gray‘, ‘cyan‘)

cmap = ListedColormap(colors[:len(np.unique(y))])

# plot the decision surface

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, xx2, Z, alpha=0.4, cmap=cmap)

plt.xlim(xx1.min(), xx1.max())

plt.ylim(xx2.min(), xx2.max())

# plot class samples

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)

7)把训练数据转换到前两个主成分轴后生成决策区域图形

1

2

3

4

5

6

7

plot_decision_regions(X_train_pca, y_train, classifier=lr)

plt.xlabel(‘PC 1‘)

plt.ylabel(‘PC 2‘)

plt.legend(loc=‘lower left‘)

plt.tight_layout()

# plt.savefig(‘./figures/pca3.png‘, dpi=300)

plt.show()

对高维数据集进行降维除了PCA方法，还有线性判别分析（Linear Discriminant Analysis， LDA）、决策树、核主成分分析、SVD等等。

17.4 LDA 降维

LDA的基本概念与PCA类似，PCA是在数据集中找到方差最大的正交的主成分分量的轴。而LDA的目标是发现可以最优化分类的特征子空间。两者都是可以用于降维的线性转换方法，其中，PCA是无监督算法，LDA是监督算法。与PCA相比，LDA是一种更优越的用于分类的特征提取技术。
LDA的主要步骤：
（1）对d维数据集进行标准化处理（d为特征数量）
（2）对每一类别，计算d维的均值向量
（3）构造类间的散布矩阵S_B以及类内的散布矩阵S_W
（4）计算矩阵〖S_W〗^(-1) S_B的特征值所对应的特征向量，
（5）选取前k个特征值对应的特征向量，构造一个d x k维的转换矩阵W，其中特征向量以列的形式排列
（6）使用转换矩阵W将样本映射到新的特征子空间上.
以下还是以下葡萄酒数据为例，用代码实现以上各步：
（1）对d维数据集进行标准化处理

1

2

3

4

5

6

7

8

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.transform(X_test)

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

1

2

3

4

5

6

7

8

#设置精度

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

运行结果
MV 1: [ 0.9259 -0.3091 0.2592 -0.7989 0.3039 0.9608 1.0515 -0.6306 0.5354 0.2209 0.4855 0.798 1.2017]

MV 2: [-0.8727 -0.3854 -0.4437 0.2481 -0.2409 -0.1059 0.0187 -0.0164 0.1095 -0.8796 0.4392 0.2776 -0.7016]

MV 3: [ 0.1637 0.8929 0.3249 0.5658 -0.01 -0.9499 -1.228 0.7436 -0.7652 0.979 -1.1698 -1.3007 -0.3912]
（3）构造类间的散布矩阵S_B以及类内的散布矩阵S_W
通过均值向量计算类内散布矩阵Sw：

通过累加各类别i的散布矩阵Si来计算：

1

2

3

4

5

6

7

8

9

10

d = 13 # number of features

S_W = np.zeros((d, d))

for label,mv in zip(range(1, 4), mean_vecs):

class_scatter = np.zeros((d, d)) # scatter matrix for each class

for row in X_train_std[y_train == label]:

row, mv = row.reshape(d, 1), mv.reshape(d, 1) # make column vectors

class_scatter += (row-mv).dot((row-mv).T)

S_W += class_scatter                             # sum class scatter matrices

print(‘Within-class scatter matrix: %sx%s‘ % (S_W.shape[0], S_W.shape[1]))

运行结果
Within-class scatter matrix: 13x13

计算各类标样本数

1

print(‘Class label distribution: %s‘ % np.bincount(y_train)[1:])

运行结果为：
Class label distribution: [40 49 35]
由此看出，各类记录数不很均匀，为此，需要对SB进行归一化处理：

1

2

3

4

5

6

d = 13 # number of features

S_W = np.zeros((d, d))

for label,mv in zip(range(1, 4), mean_vecs):

class_scatter = np.cov(X_train_std[y_train==label].T)

S_W += class_scatter

print(‘Scaled within-class scatter matrix: %sx%s‘ % (S_W.shape[0], S_W.shape[1]))

运行结果
Scaled within-class scatter matrix: 13x13

计算类间散布矩阵：

1

2

3

4

5

6

7

8

9

10

mean_overall = np.mean(X_train_std, axis=0)

d = 13 # number of features

S_B = np.zeros((d, d))

for i,mean_vec in enumerate(mean_vecs):

n = X_train[y_train==i+1, :].shape[0]

mean_vec = mean_vec.reshape(d, 1) # make column vector

mean_overall = mean_overall.reshape(d, 1) # make column vector

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

运行结果
Between-class scatter matrix: 13x13

1

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

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

求得广义特征值之后，按照降序对特征值排序

1

2

3

4

5

6

7

8

9

10

11

# 生成特征值与特征向量构成的元组

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

# Sort the (eigenvalue, eigenvector) tuples from high to low

eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)

# Visually confirm that the list is correctly sorted by decreasing eigenvalues

print(‘Eigenvalues in decreasing order:\n‘)

for eigen_val in eigen_pairs:

print(eigen_val[0])

运行结果
Eigenvalues in decreasing order:

452.721581245
156.43636122
7.05575044266e-14
5.68434188608e-14
3.41129233161e-14
3.40797229523e-14
3.40797229523e-14
1.16775565372e-14
1.16775565372e-14
8.59477909861e-15
8.59477909861e-15
4.24523361436e-15
2.6858909629e-15
d x d维协方差矩阵的秩最大为d-1，得到两个非0的特征值。
与PCA一样，我们可视化各特征贡献率

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

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.ylabel(‘"discriminability" ratio‘)

plt.xlabel(‘Linear Discriminants‘)

plt.ylim([-0.1, 1.1])

plt.legend(loc=‘best‘)

plt.tight_layout()

# plt.savefig(‘./figures/lda1.png‘, dpi=300)

plt.show()

运行结果

（6）使用转换矩阵W将样本映射到新的特征子空间上.
由上面两个新得到两个特征构成一个新矩阵

1

2

3

w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,

eigen_pairs[1][1][:, np.newaxis].real))

print(‘Matrix W:\n‘, w)

d x d维协方差矩阵的秩最大为d-1，得到两个非0的特征值。Matrix W:
[[-0.0662 -0.3797]
[ 0.0386 -0.2206]
[-0.0217 -0.3816]
[ 0.184 0.3018]
[-0.0034 0.0141]
[ 0.2326 0.0234]
[-0.7747 0.1869]
[-0.0811 0.0696]
[ 0.0875 0.1796]
[ 0.185 -0.284 ]
[-0.066 0.2349]
[-0.3805 0.073 ]
[-0.3285 -0.5971]]
将样本映射到新的特征空间

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

X_train_lda = X_train_std.dot(w)

colors = [‘r‘, ‘b‘, ‘g‘]

markers = [‘s‘, ‘x‘, ‘o‘]

for l, c, m in zip(np.unique(y_train), colors, markers):

plt.scatter(X_train_lda[y_train==l, 0] * (-1),

X_train_lda[y_train==l, 1] * (-1),

c=c, label=l, marker=m)

plt.xlabel(‘LD 1‘)

plt.ylabel(‘LD 2‘)

plt.legend(loc=‘lower right‘)

plt.tight_layout()

# plt.savefig(‘./figures/lda2.png‘, dpi=300)

plt.show()

运行结果

17.5 利用Scikit-learn进行LDA分析

下面我们利用scikit-learn中对LDA类的实现
这里先定义一个函数，后面需要用到

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):

# setup marker generator and color map

markers = (‘s‘, ‘x‘, ‘o‘, ‘^‘, ‘v‘)

colors = (‘red‘, ‘blue‘, ‘lightgreen‘, ‘gray‘, ‘cyan‘)

cmap = ListedColormap(colors[:len(np.unique(y))])

# plot the decision surface

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, xx2, Z, alpha=0.4, cmap=cmap)

plt.xlim(xx1.min(), xx1.max())

plt.ylim(xx2.min(), xx2.max())

# plot class samples

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)

对数据先LDA处理，然后用逻辑回归进行分类。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

lda = LDA(n_components=2)

X_train_lda = lda.fit_transform(X_train_std, y_train)

# 逻辑回归在相对低维数据上的表现

from sklearn.linear_model import LogisticRegression

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.tight_layout()

# plt.savefig(‘./images/lda3.png‘, dpi=300)

plt.show()

运行结果

还有几个点划分错误，下面通过正则化，效果将更好

1

2

3

4

5

6

7

8

9

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.tight_layout()

# plt.savefig(‘./images/lda4.png‘, dpi=300)

plt.show()

运行结果

17.6使用核PCA降维

前面我们介绍了两种降维方法，PCA及LDA.这两种方法，如果用于线性不可分数据集上进行分类，效果往往不很理想，原因是通过他们无法把线性不可分数据集变为线性可分数据集。如果遇到线性不可分数据集（这样的数据集往往比较普遍），有什么好方法，既降维，又可把线性不可分数据集变为线性可分数据集？
在SVM中，我们了解到核函数的神奇，把可以通过把线性不可分的数据集映射到一个高维空间，变得线性可分。基于这点，如果我们在降维时也采用核技术是否也可以呢？可以的，这就是接下来我们要介绍的内容---核PCA.
核PCA=核技术+PCA,具体步骤如下：
（1）计算核矩阵，也就是计算任意两个训练样本。这里以向基核函数（RBF）为例
经向基函数核（又称高斯核）为：

得到以下矩阵：

（2）对核矩阵K进行中心化处理

其中,是n*n的矩阵，n=训练集样本数，中每个元素都等于.l_n中的每个元素都是1/n
（3）求核矩阵的特征向量，并按降序排列，提取前k个特征向量。
不同于标准PCA，这里的特征向量并不是主成分轴。
下面我们根据以上三个步骤，实现一个核PCA。借助SciPy和NumPy，其实实现核PCA很简单。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

from scipy.spatial.distance import pdist, squareform

from scipy import exp

from scipy.linalg import eigh

import numpy as np

def rbf_kernel_pca(X, gamma, n_components):

"""

RBF kernel PCA implementation.

Parameters

------------

X: {NumPy ndarray}, shape = [n_samples, n_features]

gamma: float

Tuning parameter of the RBF kernel

n_components: int

Number of principal components to return

Returns

------------

X_pc: {NumPy ndarray}, shape = [n_samples, k_features]

Projected dataset

"""

# Calculate pairwise squared Euclidean distances

# in the MxN dimensional dataset.

sq_dists = pdist(X, ‘sqeuclidean‘)

# Convert pairwise distances into a square matrix.

mat_sq_dists = squareform(sq_dists)

# Compute the symmetric kernel matrix.

K = exp(-gamma * mat_sq_dists)

# Center the kernel matrix.

N = K.shape[0]

one_n = np.ones((N,N)) / N

K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

# Obtaining eigenpairs from the centered kernel matrix

# numpy.eigh returns them in sorted order

eigvals, eigvecs = eigh(K)

# Collect the top k eigenvectors (projected samples)

X_pc = np.column_stack((eigvecs[:, -i]

for i in range(1, n_components + 1)))

return X_pc

下面以一分离同心数据集为例，分别用PCA和核PCA对数据集进行处理，然后处理后的结果，具体请看以下代码及生成的图形：

1

2

3

4

5

6

7

8

9

from sklearn.datasets import make_circles

X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)

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.tight_layout()

plt.show()

这是典型线性不可数据集，现在我们分别用PCA及核PCA进行处理。
（1）用PCA处理，然后进行分类

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

from sklearn.decomposition import PCA

from sklearn.preprocessing import StandardScaler

import matplotlib.pyplot as plt

scikit_pca = PCA(n_components=2)

X_spca = scikit_pca.fit_transform(X)

plt.figure( figsize=(5,3))

plt.scatter(X_spca[y==0, 0], X_spca[y==0, 1],

color=‘red‘, marker=‘^‘, alpha=0.5)

plt.scatter(X_spca[y==1, 0], X_spca[y==1, 1],

color=‘blue‘, marker=‘o‘, alpha=0.5)

plt.xlabel(‘PC1‘)

plt.ylabel(‘PC2‘)

plt.show()

（2）用核PCA处理，然后进行分类

1

2

3

4

5

6

7

8

9

10

11

12

X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)

plt.figure( figsize=(5,3))

plt.scatter(X_kpca[y==0, 0], X_kpca[y==0, 1],

color=‘red‘, marker=‘^‘, alpha=0.5)

plt.scatter(X_kpca[y==1, 0], X_kpca[y==1, 1],

color=‘blue‘, marker=‘o‘, alpha=0.5)

plt.xlabel(‘PC1‘)

plt.ylabel(‘PC2‘)

plt.show()

（3）使用sklearn实现核PCA
源数据的图形为

这里通过核PCA把该数据变为线性可分数据集，实现代码如下：

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

from sklearn.decomposition import KernelPCA

from sklearn.datasets import make_moons

X, y = make_moons(n_samples=100, random_state=123)

scikit_kpca = KernelPCA(n_components=2, kernel=‘rbf‘, gamma=15)

X_skernpca = scikit_kpca.fit_transform(X)

plt.scatter(X_skernpca[y==0, 0], X_skernpca[y==0, 1],

color=‘red‘, marker=‘^‘, alpha=0.5)

plt.scatter(X_skernpca[y==1, 0], X_skernpca[y==1, 1],

color=‘blue‘, marker=‘o‘, alpha=0.5)

plt.xlabel(‘PC1‘)

plt.ylabel(‘PC2‘)

plt.tight_layout()

plt.show()

17.7 SVD矩阵分解

（1）SVD奇异值分解的定义
假设有一个mxn矩阵，如果存在一个分解

其中U为的mxm酉矩阵，∑为mxn的半正定对角矩阵，除了对角元素不为0，其他元素都为0，并且对角元素是从大到小排列的，前面的元素比较大，后面的很多元素接近0。这些对角元素就是奇异值。V^T为V的共轭转置矩阵，且为nxn的酉矩阵。这样的分解称为的奇异值分解，对角线上的元素称为奇异值，U称为左奇异矩阵，V^T称为右奇异矩阵。
SVD在信息检索（隐性语义索引）、图像压缩、推荐系统、金融等领域都有应用。
（2）SVD奇异值分解与特征值分解的关系
特征值分解与SVD奇异值分解的目的都是提取一个矩阵最重要的特征。然而，特征值分解只适用于方阵，而SVD奇异值分解适用于任意的矩阵，不一定是方阵。

这里M^T M和MM^T都是方阵，UU^T和VV^T都是单位矩阵，V是M^T M的特征向量，U是MM^T的特征向量。
（3）SVD奇异值分解的作用和意义
奇异值分解最大的作用就是数据的降维，当然，还有其他很多的作用，这里主要讨论数据的降维，对于mxn的M矩阵，进行奇异值分解

取其前k个非零奇异值，可以还原原来的矩阵，即前k个非零奇异值对应的奇异向量代表了矩阵的主要特征。可以表示为

17.8 用Python实现SVD，并用于图像压缩
（1）首先读取一张图片（128*128*3）：

1

2

3

4

5

6

7

8

9

10

11

12

13

#!python

#  -*- coding:utf-8 -*-

from PIL import Image

import os

import numpy as np

import matplotlib as mpl

import matplotlib.pyplot as plt

if __name__ == ‘__main__‘:

mpl.rcParams[‘font.sans-serif‘] = [u‘simHei‘]

mpl.rcParams[‘axes.unicode_minus‘] = False

A = Image.open(‘02.jpg‘)

a = np.array(A)  #转换成矩阵

（2）然后可以利用python的numpy库对彩色图像的3个通道进行SVD分解

1

2

3

4

5

numpy库中有SVD分解函数：np.linalg.svd

#由于是彩色图像，所以3通道。a的最内层数组为三个数，分别表示RGB，用来表示一个像素

u_r, sigma_r, v_r = np.linalg.svd(a[:, :, 0])

u_g, sigma_g, v_g = np.linalg.svd(a[:, :, 1])

u_b, sigma_b, v_b = np.linalg.svd(a[:, :, 2])

（3）然后便可以根据需要压缩图像（丢弃分解出来的三个矩阵中的数据），利用的奇异值个数越少，则压缩的越厉害。下面来看一下不同程度压缩后，重构图像的清晰度：

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

plt.figure(facecolor = ‘w‘, figsize = (10, 10))

# 奇异值个数依次取：1,2,...,16。来看看一下效果

K = 16

for k in range(1, K + 1):

R = restore(u_r, sigma_r, v_r, k)

G = restore(u_g, sigma_g, v_g, k)

B = restore(u_b, sigma_b, v_b, k)

I = np.stack((R, G, B), axis = 2)

# 将图片重构后的显示出来

plt.subplot(4, 4, k)

plt.imshow(I)

plt.axis(‘off‘)

plt.title(u‘奇异值个数：%d‘ %  k)

plt.suptitle(u‘SVD与图像分解‘, fontsize = 20)

plt.tight_layout(0.1, rect = (0, 0, 1, 0.92))

plt.show()

（4）其中restore函数定义为

1

2

3

4

5

6

7

8

9

10

11

12

13

def restore(u, sigma, v, k):

m = len(u)

n = len(v)

a = np.zeros((m, n))

# 重构图像

a = np.dot(u[:, :k], np.diag(sigma[:k])).dot(v[:k, :])

# 上述语句等价于：

# for i in range(k):

#     ui = u[:, i].reshape(m, 1)

#     vi = v[i].reshape(1, n)

#     a += sigma[i] * np.dot(ui, vi)

a[a < 0] = 0 a[a > 255] = 255

return np.rint(a).astype(‘uint8‘)

原文地址：https://www.cnblogs.com/Knight66666/p/12635772.html

时间： 2024-10-13 16:37:43

降维实践（PCA，LDA）的相关文章

sklearn中的降维算法PCA和SVD

sklearn中的降维算法PCA和SVD 1 概述 1.1 从什么叫“维度”说开来 1.2 sklearn中的降维算法 2 PCA与SVD 2.1 降维究竟是怎样实现? 2.2 重要参数n_components 2.2.1 迷你案例:高维数据的可视化 2.2.2 最大似然估计自选超参数 2.2.3 按信息量占比选超参数 2.3 PCA中的SVD 2.3.1 PCA中的SVD哪里来? 2.3.2 重要参数svd_solver 与 random_state 2.3.3 重要属性components_

机器学习基础与实践（三）----数据降维之PCA

写在前面:本来这篇应该是上周四更新,但是上周四写了一篇深度学习的反向传播法的过程,就推迟更新了.本来想参考PRML来写,但是发现里面涉及到比较多的数学知识,写出来可能不好理解,我决定还是用最通俗的方法解释PCA,并举一个实例一步步计算,然后再进行数学推导,最后再介绍一些变种以及相应的程序.(数学推导及变种下次再写好了) 正文: 在数据处理中,经常会遇到特征维度比样本数量多得多的情况,如果拿到实际工程中去跑,效果不一定好.一是因为冗余的特征会带来一些噪音,影响计算的结果:二是因为无关的特征会加大计

KNN PCA LDA

http://blog.csdn.net/scyscyao/article/details/5987581 这学期选了门模式识别的课.发现最常见的一种情况就是,书上写的老师ppt上写的都看不懂,然后绕了一大圈去自己查资料理解,回头看看发现,Ah-ha,原来本质的原理那么简单,自己一开始只不过被那些看似formidable的细节吓到了.所以在这里把自己所学的一些点记录下来,供备忘,也供参考. 1. K-Nearest Neighbor K-NN可以说是一种最直接的用来分类未知数据的方法.基本通过下

降维技术---PCA

数据计算和结果展示一直是数据挖掘领域的难点,一般情况下,数据都拥有超过三维,维数越多,处理上就越吃力.所以,采用降维技术对数据今夕简化一直是数据挖掘工作者感兴趣的方向. 对数据进行简化的好处:使得数据集更易于使用,降低很多算法的计算开销,去除噪声,使得结果易懂. 主成分分析法(PCA)是一种常用的降维技术.在PCA中,数据从原来的坐标系转换到了新的坐标系,新坐标系的选择是由数据本身决定的.第一个新坐标轴选择的是原始数据中方差最大的方向,第二个新坐标轴的选择和第一个坐标轴正交且具有最大方差的方向.

ML: 降维算法-PCA

        PCA (Principal Component Analysis) 主成份分析 也称为卡尔胡宁-勒夫变换(Karhunen-Loeve Transform),是一种用于探索高维数据结构的技术.PCA通常用于高维数据集的探索与可视化.还可以用于数据压缩,数据预处理等.PCA可以把可能具有相关性的高维变量合成线性无关的低维变量,称为主成分( principal components).新的低维数据集会尽可能的保留原始数据的变量.PCA将数据投射到一个低维子空间实现降维.例如,二维数

降维之pca算法

pca算法: 算法原理: pca利用的两个维度之间的关系和协方差成正比,协方差为0时,表示这两个维度无关,如果协方差越大这表明两个维度之间相关性越大,因而降维的时候, 都是找协方差最大的. 1 将XX中的数据进行零均值化,即每一列都减去其均值. 2 计算协方差矩阵C=1mXTXC=1mXTX 3 求出CC的特征值和特征向量 4 将特征向量按对应特征值大小从上到下按行排列成矩阵,取前k行组成矩阵P 5 Y=XPY=XP就是降维到k维后的数据. 代码: # coding=utf-8 import m

非监督的降维算法--PCA

PCA是一种非监督学习算法,它能够在保留大多数有用信息的情况下,有效降低数据纬度. 它主要应用在以下三个方面: 1. 提升算法速度 2. 压缩数据,减小内存.硬盘空间的消耗 3. 图示化数据,将高纬数据映射到2维或3维 总而言之,PCA干的事情就是完成一个将原始的n维数据转化到k维的映射.其中,k<n 它的核心算法如下: 1. 将数据均一化 x' = [x-mean(x)] / range(x) 2. 计算它的协方差矩阵 即:Sigma = 1/m * x' * x 3. 进行svd分解,计算特

机器学习之路：python 特征降维 主成分分析 PCA

python3 学习api使用 主成分分析方法实现降低维度 使用了网络上的数据集,我已经下载到了本地,可以去我的git上参考 git:https://github.com/linyi0604/MachineLearning 代码: 1 from sklearn.svm import LinearSVC 2 from sklearn.metrics import classification_report 3 from sklearn.decomposition import PCA 4 impo

sklearn 特征降维利器 —— PCA TSNE

同为降维工具,二者的主要区别在于, 所在的包不同(也即机制和原理不同) from sklearn.decomposition import PCA from sklearn.manifold import TSNE 因为原理不同,导致,tsne 保留下的属性信息,更具代表性,也即最能体现样本间的差异: TSNE 运行极慢,PCA 则相对较快: 因此更为一般的处理,尤其在展示(可视化)高维数据时,常常先用 PCA 进行降维,再使用 tsne: data_pca = PCA(n_components