（六）局部线性嵌入（LLE）及python实现

一、简述：

试图保持邻域内样本之间的线性关系。
如上图，$X_{n\times p}$中样本点$x_{(i)}$的坐标可通过邻域样本$x_{(j)}$，$x_{(k)}$，$x_{(l)}$线性组合重构出来，计算出取值之后再通过算出的权重矩阵重构$X_{n\times d}$.

二、算法

（1）计算线性重构的系数矩阵

原数据$X_{n\times p}=\left(\begin{array}{c}{x}_{(1)}^T\\ x_{(2)}^T\\ .\\ .\\ x_{(n)}^T\end{array}\right)$,$x_{(i)}=(x_{i1},x_{i2},...,x_{ip}{)}^T$.假设$Q_{(i)}$是样本$x_{(i)}$的$k$近邻，用$x_{(1)},x_{(2)},...x_{(k)}$来表示，则$$x_{(i)}=\sum _{j\in Q_{(i)}}{w}_{ij}{x}_{(j)}=\sum _{j=1}^k{w}_{ij}{x}_{(j)}$$其中${\sum }_{j\in Q_{(i)}}{w}_{ij}=1$。使用均方差作为回归问题的损失函数，即求：$$\arg \underset{w_i}{\min }\sum _{i=1}^n||x_{(i)}-\sum _{j\in Q_{(i)}}{w}_{ij}{x}_{(j)}|{|}_2^2\text{\hspace{1em}(1)}$$$$s.t.\sum _{j\in Q_{(i)}}{w}_{ij}=1$$将$(1)$式矩阵化，$$\begin{array}{rl}(1)& =\arg \underset{w_i}{\min }\sum _{i=1}^n||\sum _{j\in Q_{(i)}}{w}_{ij}(x_{(i)}-x_{(j)})|{|}_2^2\\ & =\arg \underset{w_i}{\min }\sum _{i=1}^n||\left(\begin{array}{c}(x_{(i)}-x_{(1)}),(x_{(i)}-x_{(1)}),...(x_{(i)}-x_{(k)})\end{array}\right)\left(\begin{array}{c}{w}_{i1}\\ w_{i2}\\ .\\ .\\ w_{ik}\end{array}\right)|{|}_2^2\\ & =\arg \underset{w_i}{\min }\sum _{i=1}^n||(X_i-N_i)w_i|{|}_2^2\\ & =\arg \underset{w_i}{\min }\sum _{i=1}^n[w_i^T(X_i-N_i{)}^T(X_i-N_i)w_i]\\ & =\arg \underset{w_i}{\min }\sum _{i=1}^n[w_i^T{S}_i{w}_i]\end{array}$$其中，$X_i=(x_{(i)},x_{(i)},...x_{(i)}),N_i=(x_{(1)},x_{(2)},...x_{(k)}),w_i=\left(\begin{array}{c}{w}_{i1}\\ w_{i2}\\ .\\ .\\ w_{ik}\end{array}\right),S_i=(X_i-N_i{)}^T(X_i-N_i)$,使用矩阵和拉格朗日乘子法求解该问题，令$$\begin{gathered} Y(w_i)=\sum _{i=1}^n[w_i^T{S}_i{w}_i]+\lambda (1_k^T{w}_i-1) \\ \frac{\partial Y(w_i)}{\partial w_i}=2w_i^T{S}_i+\lambda 1_k^T=0分子布局求导⟹2S_i{w}_i+\lambda 1_k=0⟹w_i=-\frac{1}{2}\lambda S_i^{-1}{1}_k \end{gathered}$$由$1_k^T{w}_i=1⟹-\frac{1}{2}\lambda =1_k{S}_i{1}_k^T⟹w_i=\frac{S_i^{-1}{1}_k}{1_k^T{S}_i{1}_k}$
注：计算得到的$w_i$为$k\times 1$向量,把非近邻的值设为0，得到$n\times 1$维向量。

（2）确定低维空间的坐标

$$X_{n\times p}=\left(\begin{array}{c}{x}_{(1)}^T\\ x_{(2)}^T\\ .\\ .\\ x_{(n)}^T\end{array}\right)⟹Y_{n\times d}=\left(\begin{array}{c}{x}_{(1)}^T\\ x_{(2)}^T\\ .\\ .\\ x_{(n)}^T\end{array}\right)$$求解：$$\arg \underset{Y}{\min }\sum _{i=1}^n||y_{(i)}-\sum _{j=1}^n{w}_{ij}{y}_{(j)}|{|}_2^2\text{\hspace{1em}(2)}$$$$s.t.Y^{T}Y=mI$$ $$\begin{array}{rl}(2)& =\arg \underset{Y}{\min }\sum _{i=1}^n||Y^T{I}_i-Y^T{w}_i|{|}_2^2\\ & =\arg \underset{Y}{\min }\sum _{i=1}^n||Y^T(I_i-w_i)|{|}_2^2\\ & =\arg \underset{Y}{\min }tr(Y^T(I-W)(I-W{)}^{T}Y)\\ & =\arg \underset{Y}{\min }tr(Y^{T}MY)\end{array}$$补充：$若A=(a_1,a_2,...,a_N),则{\sum }_i(a_i{)}^2={\sum }_i{a}_i^T{a}_i=tr(A{A}^T)$.使用矩阵和拉格朗日乘子法求解该问题，令$$\begin{gathered} J(y)=tr(Y^{T}MY)+\lambda (Y^{T}Y-mI) \\ \frac{\partial J(Y)}{\partial Y}=2MY+2\lambda Y=0 \end{gathered}$$即$MY={\lambda }^{{}^{\prime }}Y$,$Y$为$M$对应的特征向量矩阵。

（3）注意事项

取M的除0之外的最小特征值所对应的特征向量，第$2,3,...d+1$个特征值对应$Y=(y_2,y_3,...,y_d+1)$.

• 为什么不能为0？M的特征值为0不能反应数据特征，对应的特征向量全为1.$W^{T}e=e\to (W^T-I)e=0\to W^T-I=0\to W-I=0\to (I-W)(I-W{)}^{T}e=0e$，M的最小特征值为0，需要排除。
• 为什么取最小？Y是实对称矩阵的特征向量，故两两正交，$Y^{T}Y=mI(m看作1)$，则Y为正交矩阵,$Y^T=Y^{-1}$。求$\arg \underset{Y}{\min }tr(Y^{T}MY)$相当于求特征值的和最小。

三、python实现

（1）自己写的

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

def load_data():
    swiss_roll =datasets.make_swiss_roll(n_samples=1000)
    return swiss_roll[0],np.floor(swiss_roll[1])

#my_LLE任意计算两点之间距离
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)   #square计算各元素平方      1 按行相加  0按列相加
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    #返回任意两个点之间距离的平方
    return dist

#my_LLe近邻        
def floyd(D,n_neighbors):    #k近邻得到距离最近的指标
    D_arg = np.argsort(D,axis=1)#返回从小到大的索引值，每一列进行排序
    print(type(D_arg))
    return D_arg[:,1:n_neighbors+1] 

def linear_weight(X,X_floyd):
    n1,n2=X_floyd.shape
    n3,n4=X.shape
    D1 = np.ones((n4,n2))
    w=np.zeros((n1,n1))
    one_k=np.ones((n2,1))
    for i in range(n1):
        for j in range(n2):
            D1[:,j]=np.transpose(X[i,:]-X[X_floyd[i,j],:])
            
        a=np.dot(np.transpose(one_k),np.linalg.inv(np.dot(np.transpose(D1),D1))).dot(one_k)
        W_I=np.linalg.inv(np.dot(np.transpose(D1),D1)).dot(one_k)/a

        for k in range(n2):
            w[X_floyd[i,k],i]=W_I[k]
        
    return w

def reconstruction(w,d):
    n1,n2=w.shape
    danwei=np.eye(n1)
    M=(danwei-w).dot(np.transpose(danwei-w))
    eig_val, eig_vec = np.linalg.eig(M)     #计算矩阵特征值和特征向量
    eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(n2)]    #abs绝对值   把特征值对应的特征向量放一块
    # sort eig_vec based on eig_val from highest to lowest
    eig_pairs.sort(key=lambda x:x[0])  #升序排列
    # select the top d eig_vec
    feature=np.array([ele[1] for ele in eig_pairs[1:d+1]])   #选出d个特征值对应的向量,此时为特征向量为行向量
    #get new data
    new_data=np.transpose(feature)   #向量显示为列向量
    return new_data

def my_LLE(X,Y,Neighbors):
    
    D = cal_pairwise_dist(X)
    D[D < 0] = 0
    D = D**0.5       #距离矩阵
    X_floyd=floyd(D,Neighbors)
    w=linear_weight(X,X_floyd)
    new_data=reconstruction(w,2)
    return new_data
    
    
X,Y=load_data()
fig = plt.figure('data')
ax = Axes3D(fig)
ax.scatter(X[:, 0], X[:, 1], X[:, 2],marker='o',c=Y)


fig=plt.figure("my_LLE",figsize=(9, 9))
Neighbors=[1,2,3,4,5,15,30,100,Y.size-1]
for i,k in enumerate(Neighbors): #i位置，k代表的值
        lle=my_LLE(X,Y,k)
        ax=fig.add_subplot(3,3,i+1)
        ax.scatter(lle[:,0],lle[:,1],marker='o',c=Y,alpha=0.5) #c代表颜色
        ax.set_title("k = %d"%k)   
        plt.xticks(fontsize=10, color="darkorange")  
        plt.yticks(fontsize=10, color="darkorange") 
#plt.suptitle("LLE")   #总图标题
plt.show()
    

（2）使用manifold.LocallyLinearEmbedding

import numpy as np
import operator
import matplotlib.pyplot as plt
from sklearn import datasets,decomposition,manifold
from itertools import cycle
from mpl_toolkits.mplot3d import Axes3D
def load_data():
    swiss_roll =datasets.make_swiss_roll(n_samples=1000)
    return swiss_roll[0],np.floor(swiss_roll[1])
 
def LLE_components(*data):
    X,Y=data
    for n in [3,2,1]:
        lle=manifold.LocallyLinearEmbedding(n_components=n)  #LLE降维，n_components降的维数 #n_neighbors近邻数
        lle.fit(X)
        print("n = %d 重建误差："%n,lle.reconstruction_error_)
 
def LLE_neighbors(*data):
    X,Y=data
    Neighbors=[1,2,3,4,5,15,30,100,Y.size-1]
 
    fig=plt.figure("LLE",figsize=(9, 9))
  
    for i,k in enumerate(Neighbors): #i位置，k代表的值
        lle=manifold.LocallyLinearEmbedding(n_components=2,n_neighbors=k,eigen_solver='dense')
        X_r=lle.fit_transform(X)
        ax=fig.add_subplot(3,3,i+1)
        ax.scatter(X_r[:,0],X_r[:,1],marker='o',c=Y,alpha=0.5) #c代表颜色
        ax.set_title("k = %d"%k)   
        plt.xticks(fontsize=10, color="darkorange")  
        plt.yticks(fontsize=10, color="darkorange") 
    plt.suptitle("LLE")
    plt.show()
 
X,Y=load_data()
fig = plt.figure('data')
ax = Axes3D(fig)
ax.scatter(X[:, 0], X[:, 1], X[:, 2],marker='o',c=Y)
LLE_components(X,Y)
LLE_neighbors(X,Y)

四、附录

• 1）n_neighbors：即我们搜索样本的近邻的个数，默认是5。 n_neighbors个数越大，则建立样本局部关系的时间会越大，也就意味着算法的复杂度会增加。当然n_neighbors个数越大，则降维后样本的局部关系会保持的更好。在下一节我们可以通过具体的例子看出这一点。一般来说，如果算法运行时间可以接受，我们可以尽量选择一个比较大一些的n_neighbors
• 2）n_components：即我们降维到的维数。如果我们降维的目的是可视化，则一般可以选择2-5维。
• 3)reg ：正则化系数，在n_neighbors大于n_components时，即近邻数大于降维的维数时，由于我们的样本权重矩阵不是满秩的，LLE通过正则化来解决这个问题。默认是0.001。一般不用管这个参数。当近邻数远远的大于降维到的维数时可以考虑适当增大这个参数。
• 4）eigen_solver：特征分解的方法。有‘arpack’和‘dense’两者算法选择。当然也可以选择’auto’让scikit-learn自己选择一个合适的算法。‘arpack’和‘dense’的主要区别是‘dense’一般适合于非稀疏的矩阵分解。而‘arpack’虽然可以适应稀疏和非稀疏的矩阵分解，但在稀疏矩阵分解时会有更好算法速度。当然由于它使用一些随机思想，所以它的解可能不稳定，一般需要多选几组随机种子来尝试。
• 5）method： 即LLE的具体算法。LocallyLinearEmbedding支持4种LLE算法，分别是’standard’对应我们标准的LLE算法，'hessian’对应原理篇讲到的HLLE算法，‘modified’对应原理篇讲到的MLLE算法，‘ltsa’对应原理篇讲到的LTSA算法。默认是’standard’。一般来说HLLE/MLLE/LTSA算法在同样的近邻数n_neighbors情况下，运行时间会比标准的LLE长，当然降维的效果会稍微好一些。如果你对降维后的数据局部效果很在意，那么可以考虑使用HLLE/MLLE/LTSA或者增大n_neighbors，否则标准的LLE就可以了。需要注意的是使用MLLE要求n_neighbors > n_components，而使用HLLE要求n_neighbors > n_components * (n_components + 3) / 2
• 6）neighbors_algorithm：这个是k近邻的搜索方法，和KNN算法的使用的搜索方法一样。算法一共有三种，第一种是蛮力实现，第二种是KD树实现，第三种是球树实现。这三种方法在K近邻法(KNN)原理小结中都有讲述，如果不熟悉可以去复习下。对于这个参数，一共有4种可选输入，‘brute’对应第一种蛮力实现，‘kd_tree’对应第二种KD树实现，‘ball_tree’对应第三种的球树实现， ‘auto’则会在上面三种算法中做权衡，选择一个拟合最好的最优算法。需要注意的是，如果输入样本特征是稀疏的时候，无论我们选择哪种算法，最后scikit-learn都会去用蛮力实现‘brute’。个人的经验，如果样本少特征也少，使用默认的 ‘auto’就够了。 如果数据量很大或者特征也很多，用"auto"建树时间会很长，效率不高，建议选择KD树实现‘kd_tree’，此时如果发现‘kd_tree’速度比较慢或者已经知道样本分布不是很均匀时，可以尝试用‘ball_tree’。而如果输入样本是稀疏的，无论你选择哪个算法最后实际运行的都是‘brute’。