《lle 解析》由会员分享，可在线阅读，更多相关《lle 解析（7页珍藏版）》请在金锄头文库上搜索。
1、【转】 LLE 解析【转】 LLE 解析2011年07月08日数学建模案例分析的大作业用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
2、.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为NN矩阵，Ei=Xi-jSiWi
3、jXj为D维残差向量，极小化E(W)=i|Ei|2，满足如下约束：归一化W行和为1、局部性对于任意i，Wij=0若j不属于Si 3、求像：改变Y1,YN，极小化用W重构Y1,YN的方差 Ei=Yi-jSiWijYj为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
4、,N)-2*X*X; sorted,index = sort(distance); neighborhood = index(2:(1+K),:); 这里充分利用了矩阵运算技巧避免了编程中使用循环，不仅使代码紧凑，而且Matlab矩阵运算底层优化了性能。下面逐一解释 a)距离方阵 首先注意到这样一个事实|Xi-Xj|2=|Xi|2+|Xj|2-2，其中为内积。 X2 = sum(X.2,1)是1行N列的矩阵，第j个元素为X第j列的平方和，也就是|Xj|2。 repmat是平铺函数，repmat(X2,N,1)得到NN方阵，每行都是X2，同理repmat(X2,1,N)每列都是X2，是转置的意思
5、。X*X得到NN方阵，(i,j)位置元素为。因此distance就是我们要的距离方阵的平方。 b)省去开平凡 因为我们只需要K近邻，而不需要具体的距离值，而开平方是严格单调函数，保序，因此我们可以省去开平方，节省了计算量。 c)反向索引 sort给矩阵每列排序，并返回置换前后下标的映射关系index：index(i,j)是distance第j列第i小的元素的下标。我们要每列前K小的下标，也就是每个Xi的K近邻。取2:(1+K)是因为Xi到自己的距离总是0最小，排除掉。 任意两个点至少做一次内积O(D)，共O(N2)个点对，故复杂度：O(DN2) 2、求权 首先，因为对于任意i，Wij=0若j不
6、属于Si，故W只需要存K行N列即可。 其次，易见E(W)极小当且仅当每一求和项极小，因此我们依次计算W的每一列。固定列i，记x=Xi，w=W第i列，j=Xj，极小化|x-j=1.Kwjj|2，满足归一化约束 j=1.k wj=1。用矩阵语言描述：记B=( 1-x, k-x)为DK矩阵，G=BB为KK方阵（讲义中称之为Gram方阵，半正定，在摄动意义下总可以假设它非奇异），e=(1,1)为K维单位列向量，则问题化为 min |Bw|2也就是min wGw（二次型） s.t. ew=1 用拉格朗日乘数法求此条件极值：做辅助函数F(w,)= wGw-(ew -1) 对每个wj求偏导数令为0得Gw=e
7、，反解出w=G-1e，代入到归一化约束得 =(eG-1e)-1，即最优解w=(eG-1e)-1 G-1e 实际操作时，我们先解线性方程组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(DK2) C = C + eye(K,K)*tol*trace(C); %必要时摄动 W(:,ii) = Cones(K,1); %解方程Gw=e，高斯消去复杂度O(K3) W
8、(:,ii) = W(:,ii)/sum(W(:,ii); %解向量w归一化 end; 故总复杂度O(DNK3) 3、求像 将上一步得到的W视为NN方阵，Y为dN矩阵，Y第j列Yj为降维映射下Xi的像。mini|Yi-jWijYj|2 注意到Yi-jWijYj=j Wij (Yi-Yj)，因此若Y为最优解，则所有Yi平移任一相同的向量也是最优解，为了定解，我们不妨假设jYj=0。事实上，若jYj=Z，则有jYj-Z/N=0。 此外，Y=0为平凡最优解，为了避免这种退化情形，我们不妨假设jYjYj/N=I，即jYajYbj=N(a,b)，即Y的d个行向量，都在半径为sqrt(N)的N维单位球上，
9、且两两正交。 验证：i|Yi-jWijYj|2=i,jMijYiYj，其中M=(Mij)=(I-W)(I-W) 左 =iYiYi-2jWijYiYj+(jWijYj)(kWikYk) =iYiYi -2i,jWijYiYj +i,j,kWijWikYjYk =i,jij YiYi -i,jWijYiYj -i,jWjiYiYj +j,kiWijWikYjYk =i,j(ij-Wij-Wji+kWkiWkj)YiYj =右 而在单位球上极小化二次型i,jMijYiYj当且仅当Y每行是M最小的d个特征值对应的特征向量，排除0对应的特征向量。（d=1的情形为Rayleigh-Ritz定理。） 这步复
10、杂度取决于计算矩阵（部分）特征值/向量算法的复杂度。 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
11、= size(X); fprintf(1,LLE running on %d points in %d dimensionsn,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),:); %
12、 STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS fprintf(1,-Solving for reconstruction weights.n); if(KD) fprintf(1, note: KD; regularization will be usedn); 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(:
13、,ii),1,K); % shift ith pt to origin C = z*z; % local covariance C = C + eye(K,K)*tol*trace(C); % regularlization (KD) W(:,ii) = Cones(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
