局部线性嵌入(Local Linear Embedding, LLE)算法是一种经典的流形学习算法，不同于等度量映射(Isometric Mapping, Isomap)算法考虑结点的全局连接信息，LLE算法只考虑每个结点的k近邻信息，所以速度比Isomap快很多。LLE的思想如下图所示，LLE在低维空间保持了原始高维空间样本邻域内的线性关系，这也是局部线性的一个最直接的解释：高维数据在局部上具有欧式空间的性质，满足线性相关性。
给定LLE算法的输入χ={x1,x2,...,xN}∈RD×N，其中xi∈RD×1,i=1,...,N，即有N个样本，每个样本都是维度为D的列向量。算法输出为Y={y1,y2,...,yN}∈Rd×N，其中d<D，yi∈Rd×1,i=1,...,N。假定原始输入空间的向量xi满足xi=∑j=1kwjixji，xji为xi的第j个近邻点，wji∈(0,1)为权重，那么LLE算法认为在低维空间同样应该保持上述线性关系：yi=∑j=1kwjiyji，所以LLE算法的核心就是如何根据高维空间的邻域线性关系求解具有同邻域线性关系的低维流形嵌入。
首先看高维空间。LLE算法第一步是求出原始输入空间的每个样本与其邻域的线性关系，即argminWi=1∑N∥xi−j=1∑kwjixji∥2 s.t. j=1∑kwji=1记Φ(W)=i=1∑N∥xi−j=1∑kwjixji∥2=i=1∑N∥j=1∑k(xi−xji)wji∥2=i=1∑N∥(Xi−Ni)Wi∥2=i=1∑NWiT(Xi−Ni)T(Xi−Ni)Wi s.t. WiT1k=1 其中Xi,Ni∈RD×k,Wi∈Rk×1。记S=(Xi−Ni)T(Xi−Ni)∈Rk×k，则上述优化目标变成了：Φ(W)=i=1∑NWiTSWi s.t. WiT1k=1 写出该优化目标的拉格朗日方程： L=WiTSWi+λ(WiT1k−1) 求L对Wi的偏导数，并令其等于零：∂Wi∂L=2SWi+λ1k=0 则可以求出Wi： Wi=1kTSi−11kSi−11k 上式中的分母1kTSi−11k是为了对Wi∈Rk×1进行归一化，使得Wi各维度之和为1，所以这里可以忽略系数λ，总体的权重系数W矩阵为W={W1,W2,...,WN}∈Rk×N。
然后看低维空间Y={y1,y2,...,yN}∈Rd×N，Y的每一列都是一个样本，样本长度为d，共N个样本，因为低维空间要保证局部线性与原始输入空间相同，所以可以表示为：argminYi=1∑N∥yi−j=1∑kwjiyji∥2 s.t. i=1∑Nyi=0i=1∑NyiyiT=NId×d 这里有两个约束条件，第一个是说Y被中心化，所有的样本在同一个维度上的数值之和为0，第二个约束条件指的是Y中的低维样本yi都是单位向量。现建立稀疏矩阵W′={wij′}∈RN×N，且wij′={wij, xj∈kNN(xi)0, otherwise 因此，j=1∑Nwji′yji=j=1∑kwjiyji=YWi′ 上式中，Y∈Rd×N每一列代表一个样本，Wi′∈RN×1为矩阵W′的列向量，得到的结果刚好是yi∈Rd×1。记Φ(Y)=i=1∑N∥yi−j=1∑kwjiyji∥2=i=1∑N∥Y(Ii−Wi′)∥2=i=1∑N[Y(Ii−Wi′)(Ii−Wi′)TYT]=tr(Y(Ii−Wi′)(Ii−Wi′)TYT) 这里用到了矩阵的一个技巧，假定A=[a1,a2,,...,an]，an∈Rd×1为列向量，则∑i=1Nai2=∑i=1NaiTai=tr(AAT)（结果最后为一个数值，即∑i=1d∑j=1NAij2），可以使用三元数组举一个toy problem，篇幅所限就不列举出来了。所以同样，∑i=1NyiTyi=tr(YYT)。令M=(Ii−Wi′)(Ii−Wi′)T∈RN×N，则优化问题可写为：min Φ(Y)=tr(YMYT) s.t. tr(YYT−NId×d)=0 写出拉格朗日方程：L=tr[YMYT+λ(YYT−NId×d)] 对Y求导并令其等于零，则可以得到MYT+2λYT=0 MYT=λ′YT 然后，借鉴PCA的思想，要得到降维后的d维数据集，只需要取矩阵M最小的d个特征值对应的特征向量组成的矩阵Y∈Rd×N，由于M最小的特征值一般为0，所以去掉最小的特征值，选择M的第2到第d+1特征值对应的特征向量。
然后把代码贴出来，详细代码见我的GitHub，首先是求距离矩阵然后计算K近邻的索引：
def cal_pairwise_dist(data):
expand_ = data[:, np.newaxis, :]
repeat1 = np.repeat(expand_, data.shape[0], axis=1)
repeat2 = np.swapaxes(repeat1, 0, 1)
D = np.linalg.norm(repeat1 - repeat2, ord=2, axis=-1, keepdims=True).squeeze(-1)
return D
def get_n_neighbors(data, n_neighbors=10):
dist = cal_pairwise_dist(data)
dist[dist < 0] = 0
n = dist.shape[0]
N = np.zeros((n, n_neighbors))
for i in range(n):
# np.argsort 列表从小到大的索引
index_ = np.argsort(dist[i])[1:n_neighbors+1]
N[i] = N[i] + index_
return N.astype(np.int32) # [n_features, n_neighbors]
然后是LLE算法的核心：
def lle(data, n_dims=2, n_neighbors=10):
N = get_n_neighbors(data, n_neighbors) # k近邻索引
n, D = data.shape # n_samples, n_features
# prevent Si to small
if n_neighbors > D:
tol = 1e-3
else:
tol = 0
# calculate W
W = np.zeros((n_neighbors, n))
I = np.ones((n_neighbors, 1))
for i in range(n): # data[i] => [1, n_features]
Xi = np.tile(data[i], (n_neighbors, 1)).T # [n_features, n_neighbors]
# N[i] => [1, n_neighbors]
Ni = data[N[i]].T # [n_features, n_neighbors]
Si = np.dot((Xi-Ni).T, (Xi-Ni)) # [n_neighbors, n_neighbors]
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]
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]
Y = eig_vector[:, index_]
return Y
最后贴一张实验效果图，瑞士卷降维的：