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維影象，所以最後一步偷了點小懶，直接降維成二維了！
歡迎轉載，轉載請註明出處，謝謝！