PCA算法和SVD算法简介

为什么需要降维？

我们习惯于三维空间的生活，当维度高于三维时，很难想象和描述。而且在高维空间中，许多事物的行为都迥然不同。例如，如果在一个1*1的正方形中随机选择一个点，那么这个点距离边界小于0.001的概率只有约0.4%，但是，在一个10000维的单位超立方体中，这个概率大于99.99999%，在高维超立方体中，大多数点都非常接近边界。

还有一个更麻烦的区别：如果你在单位平面中随机挑两个点，这两个点之间的平均距离大约为0.52。如果在三维的单位立方体中随机挑两个点，两点之间的平均距离大约为0.66。但是，如果在一个100万维的超立方体中随机挑两个点呢？平均距离大约为408.25。这就说明，高维数据集有很大可能是非常稀疏的，也就是说大数据训练数据实例彼此之间距离很远，简单来说就是根据训练集训练出来的模型，预测非常不可靠，训练集的维度越高，过拟合的风险越大。

基于此，我们需要对训练数据集进行降维。

不同于聚类分析，降维是另一类无监督学习算法，其中主成分分析(PCA)和奇异值分解(SVD)是两类比较经典的算法，我们这分别看看这两类算法。PCA是通过正交变换将一组由线性相关变量表示的数据转换为几个由线性无关变量表示的数据，这几个线性无关的变量就是所谓的主成分。奇异值分解是将矩阵进行分解，揭示了矩阵分析的本质变换。奇异值分解不光可以用在降维算法中的特征分解，还可以用于推荐系统和自然语言处理等多个领域。其实PCA本身就有两种实现方法，一种是基于特征值分解的，另一种就是基于奇异值分解的。sklearn就是实用SVD来实现的PCA。

特征分解

假设存在n*n矩阵A和n维向量x，λ为任意实数，矩阵特征值与特征向量定义如下：

Ax = λx，其中λ就叫作矩阵A的一个特征值，n维向量x就叫作矩阵A的特征值λ对应的特征向量。计算出矩阵的特征值和特征向量就可以对矩阵进行分解。假设矩阵A有n个特征值，λ1<=λ2<=...<=λn，每个特征值对应的特征向量为w1,w2,...,wn，由特征向量构成特征矩阵，那么矩阵A可分解为：

《机器学习》中使用的是基于特征值分解的PCA算法，我先看一下这种算法：

(1) 对m行n列的数据X按照列均值为0、方差为1进行标准化处理；

(2) 计算标准化后的X的协方差矩阵C=1/m*X*(X的转置)

(3) 计算协方差矩阵C的特征值和对应的特征向量；

(4) 将特征向量按照对应特征值大小排列成矩阵，取前k行构成的矩阵P；

(5) 计算Y=PX即可得到经过PCA降维后的k维数据。

下面我就根据这个算法来用python实现一下。

import numpy as np

# 计算协方差矩阵
def cov(X):
    m = X.shape[0]
    # 数据标准化
    X = (X-np.mean(X,axis=0))/np.var(X,axis=0)
    return 1/m*np.matmul(X.T, X)

# PCA算法实现
# X是需要进行主成分分析的矩阵，k是最后保留的主成分个数
def pca(X, k):
    cov_m = cov(X)
    # 使用Numpy进行特征值分解
    eigen_value,eigen_vector = np.linalg.eig(cov_m)
    # argsort对数组排序，返回序号数组，-1代表降序排列
    idx = eigen_value.argsort(-1) 
    print(idx)
    idx_arr = idx[:k] # 取前k个特征值的序号
    print(idx_arr)
    vectors = []
    for i in idx_arr:
        vectors.append(eigen_vector[i]) # 根据前k个特征值的序号取出特征向量
    print(X.shape)
    vectors = np.array(vectors) # 这k个特征向量组成新的数组
    print(vectors.shape)
    return np.matmul(X, vectors.T) # 两个矩阵的乘积就是降维后的结果

下面我们用这个函数对鸢尾花数据进行降维：

from sklearn import datasets
import matplotlib.pyplot as plt

iris = datasets.load_iris()
X,y = iris.data, iris.target
print(X.shape) # 输出：(150, 4)
print(y) 
# 输出：
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2]
print(iris.target_names) # 输出：['setosa' 'versicolor' 'virginica']
# 把数据降到3个主成分
X_trans = pca(X,3)
print(X_trans.shape) # 输出：(150, 3)，可以看到数据维度从4降到了3

colors = ['r','g','b'] # 三种鸢尾花数据的颜色
for x_,y_ in zip(X_trans, y):
    plt.scatter(x_[0],x_[1],color=colors[y_],lw=2)

plt.show()

可以看到三种鸢尾花还是很明显的被分开了，说明降维后的数据一样可以达到降维之前的效果。

如果用sklearn实现PCA将会更加简洁：

from sklearn import decomposition
from sklearn import datasets
import matplotlib.pyplot as plt

iris = datasets.load_iris()
X, y = iris.data, iris.target
pca = decomposition.PCA(n_components=3)
pca.fit(X)
X_trans = pca.transform(X)
colors = ['r','g','b']
for c, i, target_name in zip(colors, [0,1,2], iris.target_names):
    # print(y==i)
    # 这里X_trans[y==i,0]代表了在y==i时候的数据的第0维，
    # y==i返回的其实一个大小为150的boolean数组，y等于i的地方返回True，其余返回False
    plt.scatter(X_trans[y == i, 0], X_trans[y == i, 1],
                color=c, lw=2, label=target_name)
# 添加图例
plt.legend()
plt.show();

可以看到结果一样很明显的把三种鸢尾花分隔出来了，和我们自己实现的结果略有不同是因为实现的算法不同，sklearn是基于奇异值分解来做的PCA，不过两种方法都可以对数据进行降维。

奇异值分解

要计算特征值和特征向量有一个问题，就是矩阵A必须要是一个n*n的方阵，当遇到A是一个m*n矩阵（m不等于n时），就无法使用特征值进行分解了，当遇到非方阵矩阵的时候就需要使用SVD算法了。

假设有m*n的非方阵矩阵A，对其进行矩阵分解：

其中U为m*m矩阵，λ为m*n对角矩阵，V为n*n矩阵。同时，

也就是说，U和V都是酉矩阵，也就是矩阵的转置等于矩阵的逆矩阵。奇异值分解SVD，可以让任何一个实数矩阵进行分解，鸢尾花书《矩阵力量》中把奇异值分解用一副很生动的图来展现出来：

下面我们用numpy来实现一下奇异值分解：

import numpy as np

A = np.array([[0,1],[1,1],[1,0]])
# SVD分解
U,S,Vt = np.linalg.svd(A)
print(U)
print(S) # 结果只给出了奇异值向量，省略了奇异值矩阵中为0的部分
print(Vt)
# 输出：
[[-4.08248290e-01  7.07106781e-01  5.77350269e-01]
 [-8.16496581e-01  7.45552182e-17 -5.77350269e-01]
 [-4.08248290e-01 -7.07106781e-01  5.77350269e-01]]
[1.73205081 1.        ]
[[-0.70710678 -0.70710678]
 [-0.70710678  0.70710678]]

U是一个3X3的矩阵，Vt是一个2X2的矩阵，S是一个3X2的对角矩阵，因为 𝐴=𝑈S𝑉t，而A是一个3X2矩阵。

我们可以从这些分解后的矩阵还原出原始的矩阵：

# 先用U和奇异值向量相乘，由于返回的S是一个1行2列的向量，因此U也只取前两列，否则无法相乘，最后再和Vt做矩阵相乘
np.dot(U[:,:2]*S,Vt) 
# 输出：
array([[ 1.02120423e-16,  1.00000000e+00],
       [ 1.00000000e+00,  1.00000000e+00],
       [ 1.00000000e+00, -2.11898069e-16]])

可以看到除了精度上的小差异，结果基本还原了原始矩阵A。

我们可以使用SVD算法对图像进行压缩，对图像进行奇异值分解后，可以取出前K个奇异值对图像进行重组，实现对图像的压缩。

def get_image(u,s,v,k): # u,s,v就是奇异值分解的结果，k就是取前k个奇异值
    m = len(u)
    n = len(v[0])
    # 因为原始矩阵的第一维度和分解后的U的第一维度一样，原始矩阵的第二维度和分解后的V的第二维度一样
    a = np.zeros((m,n)) # 构建一个大小和原始矩阵一样的矩阵
    for i in range(k):
        uk = u[:,i].reshape(m,1) # 取出U第i列的数据，转成m行1列的向量
        vk = v[i].reshape(1,n) # 取出V第i行的数据，转成1行n列的向量
        a += s[i]*np.dot(uk,vk)
    # numpy的clip函数代表把数组a的元素控制在0到255之间，小于0的置为0，大于255的置为255
    a = a.clip(0,255) # 因为图像的像素值在0到255之间，所以需要将结果控制在0到255之间
    return np.array(a).astype('uint8')

下面我们实现一下压缩和重组的过程：

from PIL import Image
import matplotlib.pyplot as plt
# 展示原图
img = np.array(Image.open('lena.jpg','r'))
plt.imshow(img)
plt.show()
# 对RGB图像进行奇异值分解
u_r,s_r,v_r = np.linalg.svd(img[:,:,0]) # R分量
u_g,s_g,v_g = np.linalg.svd(img[:,:,1]) # G分量
u_b,s_b,v_b = np.linalg.svd(img[:,:,2]) # B分量
k = 30 # 取前30个奇异值
# 对图像进行重组
r = get_image(u_r,s_r,v_r,k)
g = get_image(u_g,s_g,v_g,k)
b = get_image(u_b,s_b,v_b,k)
img2 = np.stack((r,g,b),axis=2)
plt.imshow(img2)
plt.show()

奇异值k取到30时，还略有点模糊，当k取到100时，和原图就非常接近了。