LLE 解析

LLE Algorithm and Implement in Matlab

《数学建模案例分析》的大作业用LLE算法，但是原作者网站上提供的源代码有些问题，主要是因为不同版本的Matlab，内置函数eigs返回的特征向量的顺序不同：老版本对应的特征值是升序，而新版本的是降序。

在这个问题中，0总是特征值，对应的特征向量为(1,1,…,1)，这不是我们要的（如果把它放进来，则用它给swiss_roll数据集降维总得到一条粗直线。）

如果你用的是Matlab6.5，则应该把line 65的

Y = Y(:,2:d+1)'*sqrt(N);

改为

Y = Y(:,n-d:n-1)'*sqrt(N);

打了这个补丁之后，lle.m就可以正常工作了！

我认为LLE算法很漂亮，下面过一下它的设计与实现~原作者的paper用lle做关键词即可在Google上搜到。

问题的提法：

给N个D维列向量X1,…,XN，希望通过映射得到N个d(<D)维列向量Y1,…,YN，要求保持邻域关系：原来离的近的点，映射过来也近。

LLE算法的思想：如果原像Xi能够表示成他邻域内点的线性组合，则像Yi也应该用相同的组合系数表示成对应像点的线性组合。即：局部线性嵌入。

算法骨架

1、 求近邻：计算每个Xi的邻居集合Si

直接用2范数下的K近邻，其中K作为算法的参数

2、 求权：改变组合系数W，极小化局部线性表出原像的方差

W为N×N矩阵，Ei=Xi-∑{j∈Si}{WijXj}为D维残差向量，极小化E(W)=∑{i}{|Ei|^2}，满足如下约束：归一化——W行和为1、局部性——对于任意i，Wij=0若j不属于Si

3、 求像：改变Y1,…,YN，极小化用W重构Y1,…,YN的方差

Ei=Yi-∑{j∈Si}{WijYj}为d维残差向量，极小化E(Y)=∑{i}{|Ei|^2}

算法实现细节

1、 求近邻

X为D行N列的矩阵，X的第j列为输入向量Xj

向量Xi到Xj的距离为sqrt(|Xi-Xj|^2)，得到距离方阵distance，然后每行排序，取前k个，对应原下标就是近邻点的编号。在Matlab中实现如下

X2 = sum(X.^2,1);

distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;

[sorted,index] = sort(distance);

neighborhood = index(2:(1+K),:);

这里充分利用了矩阵运算技巧避免了编程中使用循环，不仅使代码紧凑，而且Matlab矩阵运算底层优化了性能。下面逐一解释——

a)       距离方阵

首先注意到这样一个事实|Xi-Xj|^2=|Xi|^2+|Xj|^2-2<Xi, Xj>，其中< , >为内积。

X2 = sum(X.^2,1)是1行N列的矩阵，第j个元素为X第j列的平方和，也就是|Xj|^2。

repmat是平铺函数，repmat(X2,N,1)得到N×N方阵，每行都是X2，同理repmat(X2',1,N)每列都是X2'，'是转置的意思。X'*X得到N×N方阵，(i,j)位置元素为<Xi, Xj>。因此distance就是我们要的距离方阵的平方。

b)      省去开平凡

因为我们只需要K近邻，而不需要具体的距离值，而开平方是严格单调函数，保序，因此我们可以省去开平方，节省了计算量。

c)      反向索引

sort给矩阵每列排序，并返回置换前后下标的映射关系index：index(i,j)是distance第j列第i小的元素的下标。我们要每列前K小的下标，也就是每个Xi的K近邻。取2:(1+K)是因为Xi到自己的距离总是0最小，排除掉。

       任意两个点至少做一次内积O(D)，共O(N^2)个点对，故复杂度：O(DN^2)

2、 求权

首先，因为对于任意i，Wij=0若j不属于Si，故W只需要存K行N列即可。

其次，易见E(W)极小当且仅当每一求和项极小，因此我们依次计算W的每一列。固定列i，记x=Xi，w=W第i列，ŋj=Xj，极小化|x-∑{j=1..K}{wjηj}|^2，满足归一化约束∑{ j=1..k }{wj}=1。用矩阵语言描述：记B=( ŋ1-x,…, ηk-x)为D×K矩阵，G=B'B为K×K方阵（讲义中称之为Gram方阵，半正定，在摄动意义下总可以假设它非奇异），e=(1,…,1)'为K维单位列向量，则问题化为——

min |Bw|^2也就是min w'Gw（二次型）

s.t. e'w=1

用拉格朗日乘数法求此条件极值：做辅助函数F(w,λ)= w'Gw-λ(e'w -1)

对每个wj求偏导数令为0得Gw=λe，反解出w=G^{-1}λe，代入到归一化约束得

λ=(e'G^{-1}e)^{-1}，即最优解w=(e'G^{-1}e)^{-1} G^{-1}e

实际操作时，我们先解线性方程组Gw=e，然后再将解向量w归一化，易见得到的就是上述最优解。

在Matlab中如下实现：

W = zeros(K,N);

for ii=1:N

   z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % 计算B

   C = z'*z;                               % 计算G，复杂度O(DK^2)

   C = C + eye(K,K)*tol*trace(C);             % 必要时摄动

   W(:,ii) = C\ones(K,1);                  % 解方程Gw=e，高斯消去复杂度O(K^3)

   W(:,ii) = W(:,ii)/sum(W(:,ii));               % 解向量w归一化

end;

故总复杂度O(DNK^3)

3、 求像

将上一步得到的W视为N×N方阵，Y为d×N矩阵，Y第j列Yj为降维映射下Xi的像。min ∑{i}{|Yi-∑{j}{WijYj}|^2}

注意到Yi-∑{j}{WijYj}=∑{j}{ Wij (Yi-Yj)}，因此若Y为最优解，则所有Yi平移任一相同的向量也是最优解，为了定解，我们不妨假设∑{j}{Yj}=0。事实上，若∑{j}{Yj}=Z，则有∑{j}{Yj-Z/N}=0。

此外，Y=0为平凡最优解，为了避免这种退化情形，我们不妨假设∑{j}{YjYj'}/N=I，即∑{j}{YajYbj}=Nδ(a,b)，即Y的d个行向量，都在半径为sqrt(N)的N维单位球上，且两两正交。

验证：∑{i}{|Yi-∑{j}{WijYj}|^2}=∑{i,j}{MijYi'Yj}，其中M=(Mij)=(I-W)'(I-W)

左

=∑{i}{Yi'Yi-2∑{j}{WijYi'Yj}+(∑{j}{WijYj'})(∑{k}{WikYk})}

=∑{i}{Yi'Yi}

-2∑{i,j}{WijYi'Yj}

+∑{i,j,k}{WijWikYj'Yk}

=∑{i,j}{δij Yi'Yi}

-∑{i,j}{WijYi'Yj}

-∑{i,j}{WjiYi'Yj}

+∑{j,k}{∑{i}{WijWik}Yj'Yk}

=∑{i,j}{(δij-Wij-Wji+∑{k}{WkiWkj})Yi'Yj}

=右

而在单位球上极小化二次型∑{i,j}{MijYi'Yj}当且仅当Y每行是M最小的d个特征值对应的特征向量，排除0对应的特征向量。（d=1的情形为Rayleigh-Ritz定理。）

这步复杂度取决于计算矩阵（部分）特征值/向量算法的复杂度。

LLE matlab代码

% LLE ALGORITHM (using K nearest neighbors)
%
% [Y] = lle(X,K,dmax)
%
% X = data as D x N matrix (D = dimensionality, N = #points)
% K = number of neighbors
% dmax = max embedding dimensionality
% Y = embedding as dmax x N matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Y] = lle(X,K,d)

[D,N] = size(X);
fprintf(1,'LLE running on %d points in %d dimensions\n',N,D);


% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS 
fprintf(1,'-->Finding %d nearest neighbours.\n',K);

X2 = sum(X.^2,1);
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;

[sorted,index] = sort(distance);
neighborhood = index(2:(1+K),:);



% STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS
fprintf(1,'-->Solving for reconstruction weights.\n');

if(K>D) 
  fprintf(1,'   [note: K>D; regularization will be used]\n'); 
  tol=1e-3; % regularlizer in case constrained fits are ill conditioned
else
  tol=0;
end

W = zeros(K,N);
for ii=1:N
   z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % shift ith pt to origin
   C = z'*z;                                        % local covariance
   C = C + eye(K,K)*tol*trace(C);                   % regularlization (K>D)
   W(:,ii) = C\ones(K,1);                           % solve Cw=1
   W(:,ii) = W(:,ii)/sum(W(:,ii));                  % enforce sum(w)=1
end;


% STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF COST MATRIX M=(I-W)'(I-W)
fprintf(1,'-->Computing embedding.\n');

% M=eye(N,N); % use a sparse matrix with storage for 4KN nonzero elements
M = sparse(1:N,1:N,ones(1,N),N,N,4*K*N); 
for ii=1:N
   w = W(:,ii);
   jj = neighborhood(:,ii);
   M(ii,jj) = M(ii,jj) - w';
   M(jj,ii) = M(jj,ii) - w;
   M(jj,jj) = M(jj,jj) + w*w';
end;

% CALCULATION OF EMBEDDING
options.disp = 0; options.isreal = 1; options.issym = 1; 
[Y,eigenvals] = eigs(M,d+1,0,options);
Y = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0


fprintf(1,'Done.\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% other possible regularizers for K>D
%   C = C + tol*diag(diag(C));                       % regularlization
%   C = C + eye(K,K)*tol*trace(C)*K;                 % regularlization