各类降维算法

一、降维的概念

降维就是指采用某种映射方法，将原高维空间中的数据点映射到低维度的空间中.

1、降维的必要性

• 多重共线性
  预测变量之间相互关联。多重共线性会导致解空间的不稳定，从而可能导致结果的不连贯
• 高维空间本身具有稀疏性
  一维正态分布有68%的值落于正负标准差之间，而在十维空间上只有0.02%。
• 过多的变量会妨碍查找规律的建立
• 仅在变量层面上分析可能会忽略变量之间的潜在联系。
  例如几个预测变量可能落入仅反映数据某一方面特征的一个组内。

2、降维的目的

减少变量的个数，缓解"维度灾难"
提供一个框架来解释结果，对数据进行可视化

3、降维算法分类

4、降维效果的评估

降维的好坏没有一个直接的标准。通常通过对数据进行降维，然后用降维后的数据进行学习，再根据学习的效果选择一个恰当的降维方式和一个合适的降维模型参数。

二、投影与线性降维

在大多数的真实问题，有些特征是强相关的，训练样例实际上可以投影到高维空间的低维子空间中。这类技术统称线性降维技术，它们侧重让不相似的点在低维表示中分开。

1、SVD算法（奇异值分解）

为什么先介绍SVD算法，因为在后面的PCA算法的实现用到了SVD算法。SVD算法不光可以用于降维算法中的特征分解，还可以用于推荐系统，以及自然语言处理等领域，是很多机器学习算法的基石。
许多数学对象可以通过将它们分解成多个组成部分或者找到它们地 一些属性来更好的理解。例如：我们可以用十进制或二进制等方式表示12，但12=2X2X3永远是对的。
在介绍奇异值分解前需要先介绍一些基础概念

（1）基本概念

特征值与特征向量
设 A 是n阶方阵，如果存在数λ和非零n维列向量 x，使得 Ax=λx 成立，则：

• 称λ是矩阵A的一个特征值
• 非零向量x称为A的对应于特征值λ的特征向量

特征分解
特征分解，又称谱分解是将矩阵（方阵）分解为由其特征值和特征向量表示的矩阵之积的方法。
先把上式Ax=λx写成(λE-A)x=0的形式（E是单位矩阵），它有非零解（x为非零解）的充分必要条件是系数行列式|λE-A|=|A-λE|=0

特征分解例题
求解矩阵A的特征值与特征向量

奇异值分解
在现实的世界中，遇到的大部分矩阵都不是方阵，比如说有N个学生，每个学生有M科成绩；有N个用户，每个用户购买了M件商品，这样形成的一个N * M的长方形矩阵。
怎样描述这样普通的矩阵的重要特征？
可以使用奇异值分解来解决。

r（或者k）是一个远小于m、n的数，这样矩阵的乘法看起来像是下面的样子：

A是一个M * N的矩阵，那么得到的U是一个M * k的矩阵，Σ是一个k * k的矩阵，V^T是一个k * N的矩阵。
U的列 和 V^T 的行（也就是V的列）分别叫做 A 的 左奇异向量（left-singular vectors）和 右奇异向量（right-singular vectors），它们自身叫左奇异矩阵和右奇异矩阵，Σ 的对角线上的值叫做 A 的奇异值（singular values）它自身叫奇异值矩阵。（Σ 的对角元素来源于A^TA或AA^T的特征值的平方根，并且是按从大到小的顺序排列的）

（2）执行步骤

step 1:对AA^T进行特征值分解，得到左奇异矩阵和奇异值矩阵。
step 2:对奇异值矩阵进行降序，更新左奇异矩阵排序
step 3:间接求部分右奇异矩阵

（3）代码实现

python实现奇异值分解

# coding = utf-8
import os
from PIL import Image
import numpy as np
FILE_DIR = os.path.dirname(os.path.abspath(__file__))

if __name__ == '__main__':
    A = Image.open(FILE_DIR+'/0.jpg')
    # 转换成矩阵
    a = np.array(A)[:, :, 0]/255.0
    
    # 特征分解
    sigma, u = np.linalg.eigh(a.dot(a.T))

    # 对特征值和左奇异矩阵进行降序排列
    sort_idx = np.argsort(sigma)[::-1]
    sigma = np.sort(sigma)[::-1]
    u = u[:,sort_idx]

    # 对特征值开平方根得到奇异值
    sigma = np.sqrt(sigma)

    # 间接求右奇异矩阵
    sigma_inv = np.linalg.inv(np.diag(sigma))
    v = sigma_inv.dot((u.T).dot(a))

    print(u,sigma,v)

用SVD做图像压缩（注意：使用numpy.linalg.svd分解得到的右奇异矩阵是一个方阵，需要取其前k行才能参与矩阵重建）

# coding = utf-8
from PIL import Image
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams['font.sans-serif'] = [u'simHei']

#========================================================
#  图像重构
#========================================================

def img_reconst(u, sigma, v, k):
    '''
    根据奇异值进行图像重建
    INPUT -> u矩阵, sigma矩阵, v矩阵, 奇异值个数
    '''
    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')

#========================================================
#  主程序
#========================================================

if __name__ == '__main__':
    A = Image.open('day10.png')
    # 转换成矩阵
    a = np.array(A)  
    # 对彩色图像的3个通道进行SVD分解
    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])
    plt.figure(facecolor = 'w', figsize = (10, 10))
    # 奇异值个数依次取：1,2,...,12。来看看一下效果
    K = 12
    for k in range(1, K + 1):
        print(k)
        R = img_reconst(u_r, sigma_r, v_r, k)
        G = img_reconst(u_g, sigma_g, v_g, k)
        B = img_reconst(u_b, sigma_b, v_b, k)
        I = np.stack((R, G, B), axis = 2)
        # 将图片重构后的显示出来
        plt.subplot(3, 4, k)
        plt.imshow(I)
        plt.axis('off')
        plt.title('奇异值个数：%d' %  k)
 
    plt.suptitle('SVD与图像分解', fontsize = 20)
    plt.tight_layout(0.1, rect = (0, 0, 1, 0.92))
    plt.show()

2、PCA算法

PCA（principal Component Analysis），主成分分析法。顾名思义，就是提取出数据中主要的成分，是一种数据压缩方法，常用于去除噪声、数据预处理，也是机器学习中常见的降维方法。
PCA把原先的n个特征用数目更少的m个特征取代，新特征是旧特征的线性组合，这些线性组合最大化样本方差。从旧特征到新特征的映射会捕获数据中的固有变异性。

PCA的本质就是找一些投影方向，使得数据在这些投影方向上的方差最大

限制：旧特征之间存在线性关联。

（1）执行步骤

step 1:对所有样本进行中心化处理。
即每一维的数据都减去该维的均值。
step 2:计算协方差矩阵及特征值、特征向量
step 3:对特征值进行排序，选取大的特征值对应的特征向量，得到新的数据集

（2）代码实现

# coding: utf-8
import numpy as np

def PCA(data, n_components):
    '''
    PCA算法
    INPUT  -> 原数据, 目标维数
    '''
    # 数据标准化
    m = np.mean(data, axis=0)
    data -= m

    # 协方差矩阵
    C = np.cov(data.T)

    # 计算矩阵特征向量特征值、特征向量，按降序排序
    eigvals, eigvecs = np.linalg.eigh(C)
    indices = np.argsort(eigvals) # 返回从小到大的索引值
    indices = indices[::-1] # 反转

    eigvals = eigvals[indices] # 特征值从大到小排列
    eigvecs = eigvecs[:, indices] # 排列对应特征向量
    eigvecs_n_max = eigvecs[:, :n_components] # 取最大的前n个特征值对应的特征向量

    # 产生新的数据矩阵
    finaldata = np.dot(data, eigvecs_n_max)
    return finaldata

#========================================================
#  主程序
#========================================================

if __name__ == '__main__':
    data = np.array([[-1,2,66,-1], 
                    [-2,6,58,-1], 
                    [-3,8,45,-2], 
                    [1,9,36,1], 
                    [2.1,10,62,1], 
                    [3.2,5,83,2]])
    newdata = PCA(data, 2)
    print(newdata)

3、LDA算法

线性判别法（Linear Discriminant Analysis，LDA）与PCA很像，但不同于PCA寻找使全部数据方差最大化的主坐标轴成分，LDA关心的是能够最大化类间区分度的坐标轴成分，即在特征空间投影降维的同时还能能保持区分类别的信息。 所以在实际应用中LDA降维不仅降低分类任务的计算量，还能减小参数估计的误差，避免过拟合。

（1）执行步骤

step 1:计算中不同类别数据的 n 维均值向量
step 2:计算散布矩阵，包括类间、类内散布矩阵
step 3:计算散布矩阵的本征向量 ,,…, 和对应的本征值 ,,…,
step 4:将本征向量按本征值大小降序排列，然后选择前 k 个最大本征值对应的本征向量，组建一个 n×k 维矩阵——即每一列就是一个本征向量 （丢掉一些并没有那么重要的向量）
step 5:用这个n×k维本征向量矩阵将样本变换到新的子空间，即输出降维的结果。这一步可以写作矩阵乘法 Y=X×W 。 X 是 N×d 维矩阵，表示 N 个样本； y 是变换到子空间后的 N×k 维样本。

（2）代码实现

# coding = utf-8
from sklearn import datasets
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pylab import *
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

def compute_withinclass_scatter(samples, mean):
    # 获取样本维数, 样本个数  
    dimens,nums = samples.shape[:2]
    # 将所有样本向量减去均值向量
    samples_mean = samples - mean
    return np.cov(samples.T)

if __name__=='__main__':
    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 计算样本均值
    names = locals()
    for label in np.unique(y):
        names['mean_'+str(label)] = np.mean(X[y==label], axis=0)  # 每一类样本的均值
        names['sw_'+str(label)] = compute_withinclass_scatter(X[y==label], names['mean_'+str(label)])
    # 求总样本均值mean_overall
    mean_overall = np.mean(X, axis=0) # 所有样本的总均值
    # 求总类内离散度矩阵s_w
    s_w = 0
    for label in np.unique(y):
        s_w += names['sw_'+str(label)]
    # 求类间离散度矩阵s_b
    s_b = 0
    for label in np.unique(y):
        s_b += (names['sw_'+str(label)]-mean_overall).dot((names['sw_'+str(label)] - mean_overall).T)
    # 计算特征矩阵和特征值
    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)

    w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real, eigen_pairs[1][1][:, np.newaxis].real))
    X_lda = X.dot(w)
    print(X_lda)

拓展阅读

PCA和LDA的适用场合
当不同类的均值有明显差异时，LDA降维方法较优；当不同类的协方差有明显差异时，PCA降维方法较优。在实际应用中也常结合LDA和PCA一起使用，先用PCA降维去消除噪声，再用LDA降维。

三、非线性降维与流形学习

所谓流形（manifold）就是一般的几何对象的总称，包括各种维数的曲线曲面等。和线性降维方法一样，流形学习把一组在高维空间中的数据在低维空间中重新表示。不同的是，在流形学习中有一个假设，就是所处理的数据采样于一个潜在的流形上，或是说对于这组数据存在一个潜在的流形。
关于流形学习(Manifold Learning)最形象的解释莫过于下面的瑞士卷（Swiss roll）数据：

如果简单的使用投影（project）降维（例如通过压平第3维），那么会变成如下左图的样子，不同类别的样例都混在了一起，而我们的预想是变成右下图的形式。

瑞士卷（Swiss roll）是二维流形的例子，它在高维（三维）空间中弯曲。更一般地，一个d维流形在n维空间弯曲（其中d＜n）。

1、基于核的——KPCA算法

核主成分分析是对PCA的一种推广，是一种非线性映射的方法。在Kernel PCA分析之中，我们认为原有数据有更高的维数，我们可以在更高维的空间（Hilbert Space）中做PCA分析（即在更高维空间里，把原始数据向不同的方向投影）。
这样做的优点有：对于在通常线性空间难于线性分类的数据点，我们有可能再更高维度上找到合适的高维线性分类平面。
下图的数据就不适合用PCA进行降维

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

data, label = make_moons(n_samples=100, random_state=123)

plt.scatter(data[label==0, 0], data[label==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(data[label==1, 0], data[label==1, 1], color='blue', marker='o', alpha=0.5)
plt.show()

用PCA先试下

接下来，我们用基于核函数（RBF核）的主成分分析来试下

（1）执行步骤

step 1:计算核矩阵
step 2:中心化核矩阵
step 3:计算中心化核矩阵的特征值、特征向量
step 4:对特征值进行排序，选取大的特征值对应的特征向量，得到新的数据集

（2）代码实现

这里使用的是rbf核，还可以使用多项式核或sigmoid核，要想得到好的结果就要不断改变核参数，rbf效果好调一些，sigmoid不易找到好的参数，多项式核在一些范围效果不错，但两个参数调整需要好好找找

# coding: utf-8
import numpy as np
from scipy.spatial.distance import pdist, squareform
from sklearn.datasets import make_moons
import matplotlib.pyplot as plt

def RBF_kernel_PCA(X, gamma, n_components):
    '''
    基于RBF核的KPCA降维
    INPUT  -> 原数据, 目标维数
    '''
    # 计算核矩阵
    sq_dists = pdist(X, 'sqeuclidean')  # 平方欧式距离
    mat_sq_dists = squareform(sq_dists) # 将冗余的距离矩阵转换为简洁的距离矩阵
    K = np.exp(-gamma * mat_sq_dists) # 计算对称核矩阵

    # 中心化核矩阵
    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)

    # 从中心化核矩阵中获取特征值、特征向量
    eigvals, eigvecs = np.linalg.eigh(K)

    # 产生新的数据矩阵
    finaldata = np.column_stack((eigvecs[:, -i] for i in range(1, n_components + 1)))
    return finaldata

#========================================================
#  主程序
#========================================================

if __name__ == '__main__':
    data, label = make_moons(n_samples=100, random_state=123)
    newdata = RBF_kernel_PCA(data, 15, 2)
    print(newdata)

    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(14,6))
    ax[0].scatter(newdata[label==0, 0], newdata[label==0, 1], color='red', marker='^', alpha=0.5)
    ax[0].scatter(newdata[label==1, 0], newdata[label==1, 1], color='blue', marker='o', alpha=0.5)
    ax[1].scatter(newdata[label==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
    ax[1].scatter(newdata[label==1, 0], np.zeros((50,1))+0.02, color='blue', marker='o', alpha=0.5)
    plt.show()

KPCA对上面数据的降维效果