Locally Linear Embedding原始碼剖析
LLE是一個很經典的降維方法，關於它的理論我就不贅述了。最近我在學習LLE的時候，查看了很多LLE相關的論文，最初偷懶想在CNKI上找兩篇中文文獻作為入門，發現那些論文都是純扯蛋，竟然連原文的意思都沒有搞懂就敢發論文，對於這些混吃等S的教授，只能豎著中指說無數遍“SHIT“才能洩我心頭之恨啊。
先推薦一篇LLE原論文作者對LLE的詳細闡述的文章，http://www.cs.nyu.edu/~roweis/lle/publications.html， 在這個網址可以找到“Nonlinear dimensionality reduction by locally linear embedding”這篇論文，還可以看到一篇對LLE做二次闡述的長文“An Introduction to Locally Linear Embedding
很多人可能在看了論文後依然心存各種疑惑，畢竟論文很多地方都是略寫了，這時候看看Matlab原始碼可能會對你帶來很大幫助。這篇文章就是為了引導大家更好地理解LLE的演算法，對原始碼做了一些分析，有些地方我理解的不一定正確，歡迎大家拍磚。
% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
先對演算法的輸入和輸出做簡單介紹，輸入D*N的一個矩陣，D是每個樣本的維數，N代表有N個樣本。注意一下資料是怎麼放的，列是樣本，行是每一維度的資料。假如你有1000幅關於人臉的大小為128*128的灰度影象，那麼D=128*128=16384，N=1000。每一列是一個16384的向量，代表一張人臉影象，每一行代表不同人臉中某個位置的畫素。dmax是你希望降維之後獲得的維數，K是你採用K近鄰演算法設定的K。輸出Y是dmax*N的矩陣，如果你把人臉影象降成2維，最後就是2*N的矩陣，如果你的人臉是由向左轉或向右轉的資料組成，最後在左邊的可能就是往左轉，在右邊就是往右轉。
function [Y] = lle(X,K,d) [D,N] = size(X); fprintf(1,'LLE running on %d points in %d dimensions\n',N,D);
X是D*N的樣本矩陣，把維數D放在D，樣本數N放在N。接下來就是構建K近鄰圖。
% 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),:);%取出每個點的K個鄰居（即距離最小的前K個鄰居），因為第一個是自身到自身的距離0，所以從第二個開始。獲取一個K*N的矩陣Matlab程式碼都是非常精簡的，當然，這也跟使用他的人有很大關係，上面這麼幾行程式碼就完成了距離矩陣構建、距離排序、構建鄰接圖與索引這麼多事，Amazing！
在Matlab中X.^2是對X矩陣中每個元素做平方運算，對於Matlab不熟悉的童鞋可以自己嘗試以下程式碼：
獲得的結果就是：
sum(X.^2,1)是將X各元素平方後按列求和，擺成一個1*N的行向量X2。將各個列向量組成的樣本看成一個整體，實際上就是：
接下來我們看下一行：
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;
首先說明一下repmat(matrix,m,n)函式，它是將matrix矩陣填充到一個m*n的矩陣，這個矩陣的每個元素就是matrix。所以repmat(X2,N,1)是擺成一個N*1的矩陣，每個元素是X2，
最後就是如下矩陣：
repmat(X2',1,N)得到的矩陣如下：
X‘*X得到的矩陣是：
我們知道兩個向量的距離公式為：
所以distance就在一行程式碼上完成了距離矩陣的計算。繼續來看這行程式碼：
[sorted,index] = sort(distance);matlab的sort函式預設是對列排序，即對每一列的元素從小到大排序，可以看下面的例子：
sorted就是按列排序後的矩陣，index對應排序之前的列索引。
neighborhood = index(2:(1+K),:);因為index儲存是相應的索引，所以對於第i列來說（實際上就是考察第i個樣本），前面的K+1個數就是到樣本i距離最小的K+1個樣本，這裡包括了自身，因為自己到自己的距離為0，所以肯定會出現在前面K個，這就是取K+1個的原因，從第2到K+1取出最近K鄰居放在neighborhood。neighborhood就是截取了index矩陣的第2行起到第K+1行，形成一個K*N的矩陣。
第二步開始，求權值矩陣W。
% 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;因為K>D的時候，C是奇異矩陣，不能求逆，所以對C進行了調整，但這種調整不能過大，以免獲得的結果偏離真值過遠。這裡採取的方法是加上這麼一個矩陣：
這是K*K的矩陣，對角線的元素值非常小，也就是保證在求解W過程中對原來的W影響非常小。
我們先回到論文，看看怎麼求得權值矩陣W。在“An Introduction to locally linear embedding”這篇論文中 是這樣表述的：
上式是目標函式，在下述限制條件下求得最小值：
其中C是由下式獲得的：
上面的式子都是顯而易見的，接下來的問題是如何通過上述條件怎麼求得W。這裡要用到拉格朗日乘子法，關於什麼是拉格朗日乘子法請自行查閱相關資料。總之，我們會得到一個如下的函式，稱為朗格朗日函式：
對這個函式求極值獲得的W就是使epsilon獲得最小值的W。對各變數分別求偏導，得到如下方程組：
將上述方程組寫成矩陣的形式，得到下式：
這裡求得的C矩陣和論文的C有點出入，作者的C矩陣對角線元素是沒有2倍因子的，因為作者並沒有寫出推導過程，直說用拉格朗日乘子法就可以得到答案，也沒說這個是近似解還是精確解，這是我很不解的地方。如果在求K近鄰的時候，把自己算上的話，C矩陣的對角線元素是為0的，那樣的話我這個C對角元素的兩倍因子也是可以去掉的。可以進一步推出：
這就求出了跟第ii個樣本相關的W(ii)，這是一個1*K的向量。現在回過頭來看看Matlab原始碼。
W = zeros(K,N);這是獲得一個所有元素為0的K*N矩陣W。
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;
上面的Z是這樣獲得的，先取出X的所有鄰接節點，即K個鄰居，然後每個鄰居減去X(ii)，求解過程如下：
C=z'*z是什麼也就很清楚了。得到的C矩陣的每一項就是論文所述的C矩陣：
求出C矩陣後，還進行了調整，這是為了避免C是奇異矩陣的情況。
C = C + eye(K,K)*tol*trace(C); % regularlization (K>D)這一行程式碼就是調整，已在前面說過。
W(:,ii) = C\ones(K,1);上面這一行就是一個所有元素為1的K*1的行向量乘以C的逆。
W(:,ii) = W(:,ii)/sum(W(:,ii));上面這一行就是使W(ii)各元素相加的和為1。
經過上述過程，第二步求解W權值矩陣就算完結了。現在看第三步，怎麼求解Y。
% 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;M是一個N*N的稀疏單位矩陣矩陣，即對角線元素為1，其餘元素均為0。它的構建方式不同，儲存方式不同，但作用還是跟普通矩陣一樣，不必深究稀疏矩陣是什麼，當普通矩陣用就是了。
w = W(:,ii);這是取出W(ii)，沒什麼好解釋道的。
jj = neighborhood(:,ii);取出ii的所有鄰居放在jj這個向量裡，它當然是K*1的向量。
根據論文的推導，M這個矩陣元素如下：
因為M已經是一個單位矩陣，所有前面的delta就可以去掉了。反過來思考，假設我們現在取出了W(ii)，然後去修改跟W(ii)相關的M元素也是等價的。可以由一下推導過程給出：
M(ii,jj) = M(ii,jj) - w';
M(jj,ii) = M(jj,ii) - w;
M(jj,jj) = M(jj,jj) + w*w';其實這幾行程式碼實現的功能就是上述過程。M本來是N*N的矩陣，但是在處理的時候，只用到了跟第ii個節點相鄰的K個分量。獲得了M，再對M求特徵值、特徵向量就比較簡單了。
% 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 0Matlab的eigs是求特徵值與特徵向量的函式，詳細資料請自行查閱相關資料。
最後取出對應Y的最小的d個特徵值的特徵向量即為Y。M是對稱矩陣，並且行列式為0，所以它一定有特徵值0，但是有特徵值0就一定有特徵向量[1 1 ...1]嗎？這一直搞不明白。