流形学习之拉普拉斯特征映射

首先，我们说一下流形学习；接着，我们重点介绍拉普拉斯特征映射算法；最后，本文将给出拉普拉斯特征特征映射算法代码。

流形学习

流形学习是一种非线性降维方法，能够从高维数据中发现低维流形结构，得到高维和低维之间的映射关系，从而实现数据的维数约简。**为什么要实现特征约简？**因为高维的数据存在数据量大并且高维数据输入到网络模型中训练难度较大，耗费时间长。
流形学习展示图，其中包括Isomap（等距映射算法），LLP（局部保持投影），LE（拉普拉斯特征映射）等

拉普拉斯特征映射

原理

拉普拉斯特征映射是一种典型的流形学习方法，它是运用图拉普拉斯概念计算高维特征集得到低维流形表示。拉普拉斯特征映射算法降维的实质就是寻找到一个平均意义上保持数据点局部邻近信息，即在原来高维特征空间中是近邻的点在低维表示中也应该是近邻的。拉普拉斯特征映射算法所产生的映射可以看作是对几何流形的一种连续离散逼近的映射，用数据点的邻域图来近似表示流形，并用Laplace-Beltrami 算子近似表示邻域图的权值矩阵，实现高维流形的最优嵌入。这样就将问题简单转化为对拉普拉斯算子的特征函数的求解。

1. 构建近邻图。如果Xi和Xj为近邻点，就在节点i和 j 之间置一条边，的 K 近邻域法构建近邻图，K的取值可以具体参考文献。
2. 确定边的权值 Wij。同样边的权值确定也有两种方法:① 热核法(Heat Kernel)。若第i个和第j个节点之间是连接的，则边的权值为。
   
   若Wij=0。如果第i个和第j个节点之间有边连接，则定义边的权值Wij=1，否则Wij=0。
   3.特征映射。假设前面所建立的近邻图是连通的，那么寻找低维嵌入的问题就归结为对广义特征向量问题的求解:Ly = λDy式中：D是对角权值矩阵，它的每个元素 Dii =∑jWji；其对角元是W矩阵的列或行的和(W 为对称矩阵)，L=D－W是对称的半正定的矩阵，称为拉普拉斯矩阵，可以看作是算子函数。经过特征分解后得到的第2到第d +1 项的特征值对应的特征向量即是算法的d维输出。

拉普拉斯特征映射优点

LE方法将问题转化成对矩阵特征值的求解，不需要进行迭代计算，和其他几种流形学习算法相比，计算量和运算时间大大减少计算量。

拉普拉斯特征映射代码（matlab）

clc
clear

load y12
X=y12;
if ~exist('no_dims', 'var')  %检测no_dim变量是否存在
        no_dims = 1;
    end
    if ~exist('k', 'var') %检测k变量是否存在
        k = 8;
    end
	if ~exist('sigma', 'var') %检测sigma变量是否存在
		sigma = 1;   %sigma为什么设为1
    end
%确定三个参数，no_dims,k,sigma这里全部是输入
    
    % Construct neighborhood graph
    disp('Constructing neighborhood graph...');   %disp是显示每个变量的值
    if size(X, 1) < 2000                          %返回矩阵A所对应的行数
        G = squareform(pdist(X, 'euclidean'));  %pdist（x）是计算x中各对行向量的相互距离，计算欧式距离
        % Compute neighbourhood graph
        [tmp, ind] = sort(G);    %sort函数是将参量中得元素按照升序和降序排列，ind是一个大小等于size(G)得数组
        for i=1:size(G, 1)
            G(i, ind((2 + k):end, i)) = 0;  
        end
        G = sparse(double(G));   %首先先将G转化为浮点型矩阵，再将举证转化为稀疏矩阵。sparse是稀疏矩阵%稀疏矩阵的作用是便于计算量大。
        G = max(G, G');          %取矩阵最大的    % Make sure distance matrix is symmetric确保距离矩阵是对称的
    else
        G = find_nn(X, k);    %K近邻求点
    end
    G = G .^ 2;        %.^2是矩阵中的每个元素都求平方
	G = G ./ max(max(G));
    
    % Compute weights (W = G)
    disp('Computing weight matrices...');    %可以直接将要显示的字符放入括号中，进行表达。
    
    % Compute Gaussian kernel (heat kernel-based weights)
    G(G ~= 0) = exp(-G(G ~= 0) / (2 * sigma ^ 2));     %高斯核公式
        
    % Construct diagonal weight matrix
    D = diag(sum(G, 2)); %构建对角权重矩阵
    
    % Compute Laplacian
    L = D - G;%G是邻接矩阵
    L(isnan(L)) = 0; D(isnan(D)) = 0;    %确定哪些数组元素为NAN
	L(isinf(L)) = 0; D(isinf(D)) = 0;    %%返回维数与A相同的数组A中元素为正无穷、负无穷时，返回1，否则返回0
    %isinf()测试某个浮点数是否是无限大，其中INF表示无穷大
    % Construct eigenmaps (solve Ly = labda*Dy)
    
    disp('Constructing Eigenmaps...');  
	tol = 0;
   
        options.disp = 0;
        options.isreal = 1;  %确定数组是否为实数数组
        options.issym = 1;   %对称
        [mappedX, lambda] = eigs(L, D, no_dims + 1, tol, options);	%为什么是这样的  %eigs函数主要是通过迭代法来求解矩阵特征值和特征向量 		% only need bottom (no_dims + 1) eigenvectors
     
    % Sort eigenvectors in ascending order
    lambda = diag(lambda); %构造对角函数
    [lambda, ind] = sort(lambda, 'ascend');%'ascend' 升序排列(默认)。ind表示最小的位置
    lambda = lambda(1:no_dims + 1);%求值
    
    % Final embedding
	mappedX3 = mappedX(:,ind(2:no_dims +1));

拉普拉斯特征映射代码（python）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from mpl_toolkits.mplot3d import Axes3D

'''
author: heucoder
email: 812860165@qq.com
date: 2019.6.13
'''

def make_swiss_roll(n_samples=100, noise=0.0, random_state=None):
    #Generate a swiss roll dataset.
    t = 1.5 * np.pi * (1 + 2 * np.random.rand(1, n_samples))
    x = t * np.cos(t)
    y = 83 * np.random.rand(1, n_samples)
    z = t * np.sin(t)
    X = np.concatenate((x, y, z))
    X += noise * np.random.randn(3, n_samples)
    X = X.T
    t = np.squeeze(t)
    return X, t

def rbf(dist, t = 1.0):
    '''
    rbf kernel function
    '''
    return np.exp(-(dist/t))

def cal_pairwise_dist(x):

    '''计算pairwise 距离, x是matrix
    (a-b)^2 = a^2 + b^2 - 2*a*b
    '''
    sum_x = np.sum(np.square(x), 1)
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    #返回任意两个点之间距离的平方
    return dist

def cal_rbf_dist(data, n_neighbors = 10, t = 1):

    dist = cal_pairwise_dist(data)
    dist[dist < 0] = 0
    n = dist.shape[0]
    rbf_dist = rbf(dist, t)

    W = np.zeros((n, n))
    for i in range(n):
        index_ = np.argsort(dist[i])[1:1+n_neighbors]
        W[i, index_] = rbf_dist[i, index_]
        W[index_, i] = rbf_dist[index_, i]

    return W

def le(data,
          n_dims = 2,
          n_neighbors = 5, t = 1.0):
    '''
    :param data: (n_samples, n_features)
    :param n_dims: target dim
    :param n_neighbors: k nearest neighbors
    :param t: a param for rbf
    :return:
    '''
    N = data.shape[0]
    W = cal_rbf_dist(data, n_neighbors, t)
    D = np.zeros_like(W)
    for i in range(N):
        D[i,i] = np.sum(W[i])

    D_inv = np.linalg.inv(D)
    L = D - W
    eig_val, eig_vec = np.linalg.eig(np.dot(D_inv, L))

    sort_index_ = np.argsort(eig_val)

    eig_val = eig_val[sort_index_]
    print("eig_val[:10]: ", eig_val[:10])

    j = 0
    while eig_val[j] < 1e-6:
        j+=1

    print("j: ", j)

    sort_index_ = sort_index_[j:j+n_dims]
    eig_val_picked = eig_val[j:j+n_dims]
    print(eig_val_picked)
    eig_vec_picked = eig_vec[:, sort_index_]

    # print("L: ")
    # print(np.dot(np.dot(eig_vec_picked.T, L), eig_vec_picked))
    # print("D: ")
    # D not equal I ???
    print(np.dot(np.dot(eig_vec_picked.T, D), eig_vec_picked))

    X_ndim = eig_vec_picked
    return X_ndim

if __name__ == '__main__':
    # X, Y = make_swiss_roll(n_samples = 2000)
    # X_ndim = le(X, n_neighbors = 5, t = 20)
    #
    # fig = plt.figure(figsize=(12,6))
    # ax1 = fig.add_subplot(121, projection='3d')
    # ax1.scatter(X[:, 0], X[:, 1], X[:, 2], c = Y)
    #
    # ax2 = fig.add_subplot(122)
    # ax2.scatter(X_ndim[:, 0], X_ndim[:, 1], c = Y)
    # plt.show()

    X = load_digits().data
    y = load_digits().target

    dist = cal_pairwise_dist(X)
    max_dist = np.max(dist)
    print("max_dist", max_dist)
    X_ndim = le(X, n_neighbors = 20, t = max_dist*0.1)
    plt.scatter(X_ndim[:, 0], X_ndim[:, 1], c = y)
    plt.savefig("LE2.png")
    plt.show()

总结

流形学习的每一种降维方法都有其各自的特点(对数据的适用性)，不存在哪种方法一定比哪种方法好。

参考文献：

1.10种常用降维算法源代码(python) https://zhuanlan.zhihu.com/p/68754729
2.黄宏臣, 韩振南, 张倩倩,等. 基于拉普拉斯特征映射的滚动轴承故障识别[J]. 振动与冲击, 2015, 34(5):128-134.