LLE解析
LLE Algorithm and Implement in Matlab
「数学モデリング事例分析」の大作業はLLEアルゴリズムを用いるが、原作者のウェブサイトで提供されたソースコードにはいくつかの問題がある.主に異なるバージョンのMatlabのため、内蔵関数eigsが返す特徴ベクトルの順序が異なる.古いバージョンに対応する特徴値は昇順であり、新しいバージョンは降順である.
この問題では、0は常にフィーチャー値であり、対応するフィーチャーベクトルは(1,1,...,1)であり、これは私たちが望んでいるものではありません(これを入れると、swiss_rollデータセットに次元を下げて太い直線を得ることができます.)
Matlab 6を使っているなら5,ライン65の
Y = Y(:,2:d+1)'*sqrt(N);
に改心
Y = Y(:,n-d:n-1)'*sqrt(N);
このパッチを打った後、lle.mは正常に動作します!
LLEアルゴリズムはとてもきれいだと思います.次はそのデザインと実現~原作者のpaperがlleをキーワードにGoogleで検索できます.
問題の提案:
N個のD次元列ベクトルX 1,...,XNを与え、マッピングによりN個のd(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がSiに属さないという制約を満たす.
3、求像:Y 1,…,YN,極小化用W再構成Y 1,…,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に注目し、ここで<,>は内積である.
X 2=sum(X.^2,1)は1行N列の行列であり、j番目の要素はX番目j列の二乗和、すなわち|Xj|^2である.
repmatはタイリング関数であり,repmat(X 2,N,1)はNを得る×N方程式は,各行X 2,同理repmat(X 2',1,N)各列がX 2','は転置の意味である.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が最小なので除外します.
任意の2つの点は少なくとも1回内積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\ones(K,1);%解方程式Gw=e,Gauss消去複雑度O(K^3)
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アルゴリズムを用いるが、原作者のウェブサイトで提供されたソースコードにはいくつかの問題がある.主に異なるバージョンのMatlabのため、内蔵関数eigsが返す特徴ベクトルの順序が異なる.古いバージョンに対応する特徴値は昇順であり、新しいバージョンは降順である.
この問題では、0は常にフィーチャー値であり、対応するフィーチャーベクトルは(1,1,...,1)であり、これは私たちが望んでいるものではありません(これを入れると、swiss_rollデータセットに次元を下げて太い直線を得ることができます.)
Matlab 6を使っているなら5,ライン65の
Y = Y(:,2:d+1)'*sqrt(N);
に改心
Y = Y(:,n-d:n-1)'*sqrt(N);
このパッチを打った後、lle.mは正常に動作します!
LLEアルゴリズムはとてもきれいだと思います.次はそのデザインと実現~原作者のpaperがlleをキーワードにGoogleで検索できます.
問題の提案:
N個のD次元列ベクトルX 1,...,XNを与え、マッピングによりN個のd(
アルゴリズムスケルトン
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がSiに属さないという制約を満たす.
3、求像:Y 1,…,YN,極小化用W再構成Y 1,…,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
X 2=sum(X.^2,1)は1行N列の行列であり、j番目の要素はX番目j列の二乗和、すなわち|Xj|^2である.
repmatはタイリング関数であり,repmat(X 2,N,1)はNを得る×N方程式は,各行X 2,同理repmat(X 2',1,N)各列がX 2','は転置の意味である.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が最小なので除外します.
任意の2つの点は少なくとも1回内積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\ones(K,1);%解方程式Gw=e,Gauss消去複雑度O(K^3)
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,D);
% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS
fprintf(1,'-->Finding %d nearest neighbours.
',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 WEIGHTS
fprintf(1,'-->Solving for reconstruction weights.
');
if(K>D)
fprintf(1,' [note: K>D; regularization will be used]
');
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)
fprintf(1,'-->Computing embedding.
');
% 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 = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0
fprintf(1,'Done.
');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% other possible regularizers for K>D
% C = C + tol*diag(diag(C)); % regularlization
% C = C + eye(K,K)*tol*trace(C)*K; % regularlization