LLE局部线性嵌入算法

LLE局部线性嵌入算法

• 导语：很久没发博文了，今日抽个小空，整理下上个学期做过的东西，写成博文，供给初学者参考：

一、LLE算法概念
LLE算法，即局部线性嵌入算法，是一种非线性降维算法，它利用线性重构的局部对称性找出高维数据空间中的非线性结构，并在保持各数据点临近位置关系情况下，把高维空间数据点映射为低维空间对应的数据点，它能够使降维后的数据保持原有的流行结构。
二、LLE算法步骤
1、求每个样本点的k个最近邻点；
2、从每个样本点的近邻点计算该样本点的局部权值矩阵；
3、由样本点的近邻点和样本点重建权值矩阵来计算样本点的输出值。
理解：对于图像X中每一个块X1，第一步是计算出每个样本点Xi的k个近邻点，即找出样本点的近邻域集合，此处计算可以通过欧几里德距离法求；第二步是计算出样本点的局部重建权值矩阵W来使得重构块X1的误差最小化；第三步将所有样本点映射到低维空间中达到降维的目的。
三、LLE算法MATLAB实现-By自己

%% LLEByczc
% a little test of image Embedding by LLE
% by xiaochou
% last modified 2016.05.4
%% 清空环境变量
close all;
clear all;
clc;
%% 读取数据
pic = imread('1.bmp');
subplot(1,2,1);
imshow(pic);
%% 处理数据
[m,n,d0] = size(pic);
%F= zeros(m*n,3,'double');
m=floor(m/4)*4; n=floor(n/4)*4;
X1 = pic(1:m,1:n,:);
h = fspecial('gaussian',7,2);
Xt = imfilter(X1,h); % 高斯滤波
Xt = Xt(1:4:end,1:4:end,:);% 降质（降采样）
subplot(1,2,2);
imshow(Xt);
[r,c,d2] = size(Xt);
k = 1;
for i = 1:r
    for j = 1:c
        F(k,1) = Xt(i,j,1);
        F(k,2) = Xt(i,j,2);
        F(k,3) = Xt(i,j,3);
        k = k+1;
    end
end
X = F(1:2:end,:)';
%X = rgb2gray(X);
[D,N] = size(X); % X是原数据 N是原数据个数，即矩阵的列 D是维数，即矩阵的行
tol = 0; % tol是绝对误差限
k = 5; % 邻域点数
d = 2; % 最大嵌入维数
if(k>D) 
  tol=1e-3; % 因为k大于D时，C是奇异矩阵，行列式为0，不能求逆
else
  tol=0;
end
%% LLE algorithm
Z = sum(X.^2,1); % 矩阵X中每个元素以2为指数求幂值，并且竖向相加，即把每列加起来
X = double(X);
distance = repmat(Z,N,1)+repmat(Z',1,N)-2*X'*X; % repmat就是在行方向把X2复制成N份，列方向1份
[distance,ID] = sort(distance); % sort是对矩阵进行排序 distance是返回对每列排序的结果 ID是返回排序后矩阵中每个数据在矩阵未排序前每列中的位置，这个ID是与低分辨率块对应的5个高分辨率块的位置索引
neighborhood = ID(2:(1+k),:);
W = zeros(k,N);
for i = 1:N
    A = X(:,neighborhood(:,i))-repmat(X(:,i),1,k); % 求A是把找到的5个相似块减去与其相似的低分辨率快的差值
    C = A'*A; % 计算局部协方差 求的是相似块之间的协方差 计算五个快的相关性
    C = C+eye(k,k)*tol*trace(C); % C+的部分迭代终止误差限 防止C过小 对数据正则化 为了数据的统一
    W(:,i) = C\ones(k,1); % 解方程Gw=e，高斯消去复杂度O(K^3)
    W(:,i) = W(:,i)/sum(W(:,i)); % 解向量w归一化 计算权重矩阵
end
M=sparse(1:N,1:N,ones(1,N),N,N,4*k*N); % M=(I-W)'(I-W); 稀疏矩阵
% 计算矩阵M=(I-W)'(I-W)的最小d个非零特征值对应的特征向量
for i=1:N
    w=W(:,i);
    j=neighborhood(:,i);
    M(i,j)=M(i,j)-w';
    M(j,i)=M(j,i)-w;
    M(j,j) = M(j,j) + w*w';
end
opts.disp=0;
opts.isreal=1; % M为实数
opts.issym=1; % M对称
[Y,eigenvals]=eigs(M,d+1,0,opts); % 返回 d+1 个最小特征值 指定选项结构 由大到小降序 0相当于sm 即取最小的特征值
% Y=Y(:,2:d+1)'*sqrt(N);
Y = Y(:,3:-1:2)'*sqrt(N);
Y = Y';
%% scatterplot of embedding
%Yt(:,:,1) = Y;
%Yt(:,:,2) = Y;
%Yt(:,:,3) = Y;
%Yt = ntsc2rgb(Yt);
%figure,imshow(Yt,[]);

四、LLE算法MATLAB实现-网络版

%% lle
% a litter test of image Embedding
% by xiaochou
% last modified 2016.03.15
% 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);
%D是矩阵的行数，N是矩阵的列数
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);
%矩阵X中的每个元素以2为指数求幂值，并且竖向相加
%if two point X=(x1,x2),Y=(y1,y2)
%than the distance between X and Y is sqtr(x1-y1)+sqtr(x2-y2)
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X; 
%repmat就是在行方向把X2复制成N份,列方向为1份
[sorted,index] = sort(distance);
%sort是对矩阵排序,sorted是返回对每列排序的结果,index是返回排序后矩阵中每个数在矩阵未排序前每列中的位置
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)
%计算矩阵M=(I-W)'(I-W)的最小d个非零特征值对应的特征向量
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,eigenvals] = jdqr(M,d+1);%change in using JQDR func
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

五、结束语
利用网络版，也就是作者提供的代码并没有能实现降维，出现了一些未知的错误，经过一些处理，自己好好修改了下代码，总算是把LLE算法理解了，顺便把降维操作实现，由于自己使用的是3维图像，所以最后一步偷了点小懒，直接降维成二维了！
欢迎转载，转载请注明出处，谢谢！