LLE算法与实现
来源:互联网 发布:nba比尔数据 编辑:程序博客网 时间:2024/05/21 09:51
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\on
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 WEIGHTSfprintf(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 conditionedelse tol=0;endW = 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\on es(K,1); % solve Cw=1 W(:,ii) = W(:,ii)/sum(W(:,ii)); % enforce sum(w)=1end;% 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 elementsM = sparse(1:N,1:N,on es(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 EMBEDDINGoptions.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 ev al 0fprintf(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
网友评论:"; }else{ if(cmtname=="" || cmtname=="匿名网友"){ if(cmturl==""){ html1="匿名网友"; }else{ html1=""; } }else{ if(cmturl==""){ html1="网友:"+cmtname+""; }else{ html1="网友:"; } } } document.write(html1); } function filterCmtContent(n){ if(!BdBrowser.isIE){ var defaultfilter1=''; var defaultfilter2=''; var commentDiv=document.getElementById(n); var divs=commentDiv.getElementsByTagName('div'); var d,tmp; for( var i=0,len=divs.length;i/ig,defaultfilter1); tmp=tmp.replace('',defaultfilter2); d.innerHTML=tmp; } } } } 12009/07/27 10:52 | 回复22009/10/09 16:49 | 回复
- LLE算法与实现
- LLE算法与实现
- LLE算法及实现
- LLE算法
- LLE算法
- LLE算法
- LLE算法
- LLE与Sam Roweis
- LLE算法及一个例子
- LLE局部线性嵌入算法
- LLE算法及一个例子
- LLE及其改进算法介绍
- LLE(线性局部嵌)入算法
- 局部线性嵌入(LLE)算法解释
- Locally linear embedding (LLE)算法简介
- LLE流行嵌入式降维算法
- 流形学习——局部线性嵌入算法LLE
- LLE 解析
- 黑马程序员——IO流
- Introduce srvctl command
- 如何在PPT中插入压缩包
- LCA向RMQ转化
- 覆盖
- LLE算法与实现
- VirtualBox使用教程图解
- make介绍与使用
- 跟工作环境相关的一些问题
- Google Code Jam 2012 Qualification Round
- 匿名类
- 黑马程序员--多线程
- 局域网上网配置流程
- hadoop wordcount 相关资料