学习日记_241110_局部线性嵌入（Locally Linear Embedding, LLE）

前言

提醒：
文章内容为方便作者自己后日复习与查阅而进行的书写与发布，其中引用内容都会使用链接表明出处（如有侵权问题，请及时联系）。
其中内容多为一次书写，缺少检查与订正，如有问题或其他拓展及意见建议，欢迎评论区讨论交流。

相关链接：
LLE原理及推导过程
局部线性嵌入原理总结
代码抄录源

代码抄录

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

def make_swiss_roll(n_samples=100,noise=0.0,random_state=None):
    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 cal_pairwise_dist(x):
    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 get_n_neighbors(data,n_neighbors=10):
    dist=cal_pairwise_dist(data)
    dist[dist<0]=0
    dist=dist**0.5
    n=dist.shape[0]
    N=np.zeros((n,n_neighbors))
    for i in range(n):
        index_=np.argsort(dist[i])[1:n_neighbors+1]
        N[i]=N[i]+index_
    return N.astype(np.int32)
def lle(data,n_dims=2,n_neighbors=10):
    N=get_n_neighbors(data,n_neighbors)
    n,D=data.shape
    if n_neighbors>D:
        tol=1e-3
    else:
        tol=0
    W=np.zeros((n_neighbors,n))
    I=np.ones((n_neighbors,1))
    for i in range(n):
        Xi=np.tile(data[i],(n_neighbors,1)).T
        Ni=data[N[i]].T
        Si=np.dot((Xi-Ni).T,(Xi-Ni))
        Si=Si+np.eye(n_neighbors)*tol*np.trace(Si)
        Si_inv=np.linalg.pinv(Si)
        wi=(np.dot(Si_inv,I))/(np.dot(np.dot(I.T,Si_inv),I)[0,0])
        W[:,i]=wi[:,0]
    print("Xi.shape:",Xi.shape)
    print("Ni.shape:",Ni.shape)
    print("Si.shape:",Si.shape)
    W_y=np.zeros((n,n))
    for i in range(n):
        index=N[i]
        for j in range(n_neighbors):
            W_y[index[j],i]=W[j,i]
    I_y=np.eye(n)
    M=np.dot((I_y-W_y),(I_y-W_y).T)

    eig_val,eig_vector=np.linalg.eig(M)
    index_=np.argsort(np.abs(eig_val))[1:n_dims+1]
    print("index_:",index_)
    Y=eig_vector[:,index_]
    return Y
    
if __name__=="__main__":
    X,Y=make_swiss_roll(n_samples=500,noise=0.1,random_state=42)
    data_1=lle(X,n_neighbors=30)
    print(data_1.shape)
    data_2=LocallyLinearEmbedding(n_components=2,n_neighbors=30).fit_transform(X)

    plt.figure(figsize=(8,4))
    plt.subplot(121)
    plt.title("LLE")
    plt.scatter(data_1[:,0],data_1[:,1],c=Y)

    plt.subplot(122)
    plt.title("sklearn LLE")
    plt.scatter(data_2[:,0],data_2[:,1],c=Y)
    plt.show()

代码分析

LLE算法介绍

局部线性嵌入（Locally Linear Embedding, LLE）是一种非线性降维算法，致力于保持高维数据的局部几何结构。它通常用于揭示数据的内在流形结构。下面是对 LLE 算法的详细介绍，包括其数学原理和步骤。

LLE 算法步骤

1. 构建邻域

   对于每个数据点 ${\mathbf{x}}_i$，确定其 $k$ 个最近邻。最近邻的选择通常基于欧氏距离或其他距离度量。

2. 线性重建权重

   以每个数据点及其邻居为基础，构建重建权重矩阵 $\mathbf{W}$。重建权重使得数据点可以由其邻居的线性组合表示，即最小化以下重构误差：

   $$ϵ(W)=\sum _{i=1}^N{\Vert {\mathbf{x}}_i-\sum _{j\in \mathcal{N}(i)}{W}_{ij}{\mathbf{x}}_j\Vert }^2$$

   其中 $\mathcal{N}(i)$ 是点 ${\mathbf{x}}_i$ 的邻居集合，权重 $W_{ij}$ 满足归一化条件：

   $$\sum _{j\in \mathcal{N}(i)}{W}_{ij}=1$$

   上述优化问题可以通过局部线性系统求解，常用的方法是通过拉格朗日乘子或约束二次规划解决。

3. 低维嵌入

   找到低维坐标 ${\mathbf{y}}_i$ ，以最小化以下嵌入误差：

   $$\Phi (\mathbf{Y})=\sum _{i=1}^N{\Vert {\mathbf{y}}_i-\sum _{j\in \mathcal{N}(i)}{W}_{ij}{\mathbf{y}}_j\Vert }^2$$

   该问题转化为求解特征值问题：

   $$\mathbf{MY}=\lambda \mathbf{Y}$$

   其中矩阵 $\mathbf{M}=(\mathbf{I}-\mathbf{W}{)}^T(\mathbf{I}-\mathbf{W})$，$\mathbf{I}$ 是单位矩阵。

4. 求解特征值问题

   通过求解上述特征值问题，选择对应于最小非零特征值的特征向量作为低维嵌入坐标。通常，去掉最小的零特征值对应的特征向量（即常数特征向量）。

数学性质

• LLE 本质上是保持局部线性关系的降维方法，假设高维数据位于低维流形上，且该流形在局部近似为线性。
• 通过保持数据点在其邻居上的重建关系，LLE 能够揭示数据的内在几何结构。
• LLE 不需要明确的参数化模型，也不需要迭代优化及非凸目标的求解。

优点

• 能够有效处理非线性流形结构。
• 保持数据的局部几何结构。
• 不依赖于参数化模型，减少了对模型选择的依赖。

缺点

• 对噪声和参数选择（如邻居数量）敏感。
• 算法复杂度较高，特别是在处理大规模数据时。

LLE 算法通过保持数据的局部线性关系，能够有效地实现高维到低维的非线性降维，被广泛应用于模式识别、图像处理等领域。

def make_swiss_roll(n_samples=100,noise=0.0,random_state=None)

def make_swiss_roll(n_samples=100,noise=0.0,random_state=None):
    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

make_swiss_roll函数生成一个三维的瑞士卷数据集。这个数据集是一个形似瑞士卷蛋糕的螺旋形状，常用于测试机器学习算法。函数的数学表达可以分为以下几步：

1. 参数:

• $n\_samples$: 样本数。
• $noise$: 噪声水平。
• $random\_state$: 随机种子。

2. 生成参数 $t$:

• $t=1.5\pi (1+2U)$, 其中 $U$ 是一个均匀分布在 $[0,1]$ 上的随机数。

3. 生成数据:

• 第一维 (x轴): $x=t\cos (t)$
• 第二维 (y轴): $y=83V$, 其中 $V$ 是一个均匀分布在 $[0,1]$ 上的随机数。
• 第三维 (z轴): $z=t\sin (t)$

4. 合并数据:

• 数据点 $X=(x,y,z)$ 通过将三个维度的值合并形成一个三维点。

5. 添加噪声:

• $X$ 的每个坐标加上正态分布的随机噪声 $\mathcal{N}(0,\text{noise})$。

6. 输出:

• 返回值是一个由样本点组成的矩阵 $X$，以及对应的参数 $t$，其中 $X$ 是形状为 $(n\_samples,3)$ 的矩阵。

因此，这个函数生成的瑞士卷的数学公式为：

$$X=\left[\begin{array}{c}t\cos (t)+\text{noise}\cdot N(0,1)\\ 83V+\text{noise}\cdot N(0,1)\\ t\sin (t)+\text{noise}\cdot N(0,1)\end{array}\right]$$

其中 $V\sim U(0,1)$ 是在 $[0,1]$ 上的均匀随机数，$N(0,1)$ 是正态分布的随机噪声。

使以下代码生成图像：

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

def make_swiss_roll(n_samples=100, noise=0.0, random_state=None):
    # 设置随机种子以保证结果的可复现性
    if random_state is not None:
        np.random.seed(random_state)

    # 生成瑞士卷数据
    t = 1.5 * np.pi * (1 + 2 * np.random.rand(1, n_samples))  # 角度
    x = t * np.cos(t)  # 瑞士卷的X坐标
    y = 83 * np.random.rand(1, n_samples)  # 瑞士卷的Y坐标
    z = t * np.sin(t)  # 瑞士卷的Z坐标

    # 将坐标组合成一个矩阵，并添加噪声
    X = np.concatenate((x, y, z))
    X += noise * np.random.randn(3, n_samples)
    X = X.T  # 转置以匹配样本数和特征数
    t = np.squeeze(t)  # 压缩维度
    return X, t

# 生成瑞士卷数据
X, t = make_swiss_roll()

# 绘制三维瑞士卷
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=t, cmap='plasma')  # 使用颜色映射
ax.set_title("Swiss Roll in 3D")
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.show()

图片如下所示：

def cal_pairwise_dist(x)

def cal_pairwise_dist(x):
    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

得到一个$n\times n$的矩阵，其中每个元素是对应点对之间的欧几里得距离的平方：
$$\text{dist}(i,j)=\text{sum\_x}[i]+\text{sum\_x}[j]-2(x\cdot x^T{)}_{ij}$$

def get_n_neighbors(data,n_neighbors=10)

def get_n_neighbors(data,n_neighbors=10):
    dist=cal_pairwise_dist(data)
    dist[dist<0]=0
    dist=dist**0.5
    n=dist.shape[0]
    N=np.zeros((n,n_neighbors))
    for i in range(n):
        index_=np.argsort(dist[i])[1:n_neighbors+1]
        N[i]=N[i]+index_
    return N.astype(np.int32)

输出可以表达为：

N [ i ] = argsort ( D [ i ] ) [ 1 : n neighbors + 1 ] 对于所有  i ∈ { 1 , 2 , … , n } N[i] = \text{argsort}(D[i])[1:n_{\text{neighbors}} + 1] \quad \text{对于所有 } i \in \{1, 2, \ldots, n\} N[i]=argsort(D[i])[1:nneighbors​+1]对于所有 i∈{1,2,…,n}
描述了如何为每个数据点找到 $n_{\text{neighbors}}$ 个最近邻的索引，并将结果存储在矩阵$N$中。

关于代码：
for i in range(n):
index_=np.argsort(dist[i])[1:n_neighbors+1]
N[i]=N[i]+index_

• N 是一个二维数组，形状为 (n, n_neighbors)，用于存储每个数据点的最近邻索引。
• N[i] 是矩阵 $N$ 的第$i$行，用于存储第$i$ 个数据点的最近邻的索引。
• index_ 是一个包含$n_{\text{neighbors}}$ 个最近邻索引的一维数组。
• N[i] = N[i] + index_ 实际上是将 index_ 中的值赋给 N[i]。这里 N[i] 的初始值是一个零向量（因为 N 是通过 np.zeros 初始化的），加上 index_ 之后，结果就是 index_。

def lle(data,n_dims=2,n_neighbors=10)

def lle(data,n_dims=2,n_neighbors=10):
    N=get_n_neighbors(data,n_neighbors)
    n,D=data.shape
    if n_neighbors>D:
        tol=1e-3
    else:
        tol=0
    W=np.zeros((n_neighbors,n))
    I=np.ones((n_neighbors,1))
    for i in range(n):
        Xi=np.tile(data[i],(n_neighbors,1)).T
        Ni=data[N[i]].T
        Si=np.dot((Xi-Ni).T,(Xi-Ni))
        Si=Si+np.eye(n_neighbors)*tol*np.trace(Si)
        Si_inv=np.linalg.pinv(Si)
        wi=(np.dot(Si_inv,I))/(np.dot(np.dot(I.T,Si_inv),I)[0,0])
        W[:,i]=wi[:,0]
    print("Xi.shape:",Xi.shape)
    print("Ni.shape:",Ni.shape)
    print("Si.shape:",Si.shape)
    W_y=np.zeros((n,n))
    for i in range(n):
        index=N[i]
        for j in range(n_neighbors):
            W_y[index[j],i]=W[j,i]
    I_y=np.eye(n)
    M=np.dot((I_y-W_y),(I_y-W_y).T)

    eig_val,eig_vector=np.linalg.eig(M)
    index_=np.argsort(np.abs(eig_val))[1:n_dims+1]
    print("index_:",index_)
    Y=eig_vector[:,index_]
    return Y

关于代码 Xi=np.tile(data[i],(n_neighbors,1)).T：

1. data[i]：
   • data 是一个二维数组，形状为 $(n,D)$，表示 $n$ 个样本，每个样本有 $D$ 个特征。
   • data[i] 选择第 $i$ 个样本，结果是一个一维数组，形状为 $(D,)$。
2. np.tile(data[i], (n_neighbors, 1))：
   • np.tile 函数用于重复数组。
   • (n_neighbors, 1) 是重复参数，表示将 data[i] 沿着第一个维度重复 $n_{\text{neighbors}}$ 次，沿着第二个维度重复 1 次。
   • 结果是一个二维数组，形状为 $(n_{\text{neighbors}},D)$，其中每一行都是 data[i]。

输入

• 数据集 $\mathbf{X}\in {\mathbb{R}}^{n\times D}$，其中 $n$ 是样本数量，$D$ 是特征维度。
• 目标维度 $n_{\text{dims}}$。
• 最近邻数量 $n_{\text{neighbors}}$。

步骤

1. 找到最近邻： 使用 get_n_neighbors 函数找到每个数据点的最近邻索引： $$\mathbf{N}=\text{get\_n\_neighbors}(\mathbf{X},n_{\text{neighbors}})$$

2. 计算权重矩阵 $\mathbf{W}$：

   对于每个数据点 ${\mathbf{x}}_i$：

   • 取出 ${\mathbf{x}}_i$ 的最近邻 ${\mathbf{N}}_i$，用 ${\mathbf{X}}_{{\mathbf{N}}_i}$ 表示这些邻居。
   • 构造局部协方差矩阵：
     $${\mathbf{S}}_i=({\mathbf{X}}_i-{\mathbf{X}}_{{\mathbf{N}}_i}{)}^T({\mathbf{X}}_i-{\mathbf{X}}_{{\mathbf{N}}_i})$$
     其中 ${\mathbf{X}}_i$ 是将数据点 ${\mathbf{x}}_i$ 复制为适当形状的矩阵。
   • 添加正则化项以确保 ${\mathbf{S}}_i$ 可逆（根据 ${\mathbf{S}}_i$ 的迹）：
     $${\mathbf{S}}_i={\mathbf{S}}_i+\text{tol}\cdot \text{tr}({\mathbf{S}}_i)\cdot \mathbf{I}$$
   • 计算伪逆：
     $${\mathbf{S}}_i^{-1}=\text{pinv}({\mathbf{S}}_i)$$
   • 计算加权系数：
     $${\mathbf{w}}_i=\frac{{\mathbf{S}}_i^{-1}1}{1^T{\mathbf{S}}_i^{-1}1}$$
     其中 $1$ 是一个全为1的列向量。
   • 将权重 ${\mathbf{w}}_i$ 存储在权重矩阵 $\mathbf{W}$ 中。

3. 构建嵌入矩阵 $\mathbf{M}$：

   构造稀疏矩阵 ${\mathbf{W}}_y$： $${\mathbf{W}}_y=\sum _{i=1}^n\sum _{j\in {\mathbf{N}}_i}{\mathbf{w}}_{ij}{\mathbf{e}}_j{\mathbf{e}}_i^T$$ 其中 ${\mathbf{e}}_j$ 和 ${\mathbf{e}}_i$ 为标准基向量。

   构造矩阵 $\mathbf{M}$： $$\mathbf{M}=(\mathbf{I}-{\mathbf{W}}_y)(\mathbf{I}-{\mathbf{W}}_y{)}^T$$

4. 特征值分解：

   求解特征值问题： $$\mathbf{Mv}=\lambda \mathbf{v}$$
   选择对应于最小的 $n_{\text{dims}}+1$ 个特征值的特征向量（忽略最小的那个特征值，对应于平移不变性）。

5. 生成低维嵌入：

   最终的低维嵌入是特征向量矩阵的前 $n_{\text{dims}}$ 列： $$\mathbf{Y}=\mathbf{V}[:,\text{index\_}]$$

输出

• 低维嵌入 $\mathbf{Y}\in {\mathbb{R}}^{n\times n_{\text{dims}}}$。

python函数库

data_2=LocallyLinearEmbedding(n_components=2,n_neighbors=30).fit_transform(X)

代码运行结果

遗留问题

1. 关于权重矩阵与嵌入矩阵，，