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

前言

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

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

文章目录

前言
代码抄录
代码分析
- LLE算法介绍
- - LLE 算法步骤
  - 数学性质
  - 优点
  - 缺点
- def make_swiss_roll(n_samples=100,noise=0.0,random_state=None)
- def cal_pairwise_dist(x)
- def get_n_neighbors(data,n_neighbors=10)
- def lle(data,n_dims=2,n_neighbors=10)
- - 输入
  - 步骤
  - 输出
- python函数库
- 代码运行结果
遗留问题

代码抄录

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 算法步骤

构建邻域

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

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

$\epsilon(W) = \sum_{i=1}^{N} \left\| \mathbf{x}_i - \sum_{j \in \mathcal{N}(i)} W_{ij} \mathbf{x}_j \right\|^2$

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

$\sum_{j \in \mathcal{N}(i)} W_{ij} = 1$

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

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

$\Phi(\mathbf{Y}) = \sum_{i=1}^{N} \left\| \mathbf{y}_i - \sum_{j \in \mathcal{N}(i)} W_{ij} \mathbf{y}_j \right\|^2$

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

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

其中矩阵 $\mathbf{M} = (\mathbf{I} - \mathbf{W})^T (\mathbf{I} - \mathbf{W})$ ， $\mathbf{I}$ 是单位矩阵。
求解特征值问题

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

数学性质

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

参数:

$n\_samples$ : 样本数。
$n o i se$ : 噪声水平。
$random\_state$ : 随机种子。

生成参数 $t$ :

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

生成数据:

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

合并数据:

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

添加噪声:

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

输出:

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

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

$\begin{bmatrix} 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{bmatrix}$

其中 $\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

得到一个 $\times n$ 的矩阵，其中每个元素是对应点对之间的欧几里得距离的平方：
$sum_x [ i ] + sum_x [ j ] − 2 ( x ⋅ x T ) i j \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)

输出可以表达为：

$\text{argsort}(D[i])[1:n_{\text{neighbors}} + 1] \quad \text{对于所有 } i \in \{1, 2, \ldots, 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：

data[i]：
data 是一个二维数组，形状为 $(n, D)$ ，表示 $n$ 个样本，每个样本有 $D$ 个特征。
data[i] 选择第 $i$ 个样本，结果是一个一维数组，形状为 $(D,)$ 。

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}}$ 。

步骤

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

计算权重矩阵 $\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} \mathbf{1}}{\mathbf{1}^T \mathbf{S}_i^{-1} \mathbf{1}}$
其中 $\mathbf{1}$ 是一个全为1的列向量。
将权重 $\mathbf{w}_i$ 存储在权重矩阵 $\mathbf{W}$ 中。

构建嵌入矩阵 $\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$

特征值分解：

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

生成低维嵌入：

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

输出

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