流形学习
一、什么是流形学习?介绍两个常见的非线性流行学习算法的原理和思想;
流形学习:
流行(Manifold)= Many+Fold,字面意思就是很多曲面叠加;流形学习的观点认为我们所观察到的数据实际上是由一个低维流行映射到高维空间的;一个流行就好比是一个d维空间,在一个m维的空间中(m>d)被扭曲之后的结果。需要注意的是流形不是一个形状,而是一个空间,欧式空间是流行的一种特殊情况;
相对的,传统的数据降维很多算法都是用欧氏距离作为评价两个点之间距离的函数;换句话说就是用数据在某个方向上投影的距离最大来表示两点间的距离,显然在二维平面上两点间得距离可以用欧氏距离俩衡量,然而在高维空间上就会产生错误;用地球做例子来说就是,地球上两点之间的距离如果用欧氏距离计算的话,会得到两点间的距离,但是那个距离本身是不可达的,也就是没有意义的;由此可见传统方法在高纬度的问题上有一定的局限性(PCA无法处理非线性的流行);我们可以将高维空间的数据映射到低维空间(MDS)然后再使用欧氏距离来计算两点之间的距离,映射的过程也就Unfold a manifold,也就是在保持流行的几何结构的基础上展开一个流行;
ISOMAP(isometeric feature mapping)
构造邻接图,近点之间的测地距离用欧氏距离逼近,远点之间的测地距离用最短路径逼近离(Dijkstra或者Floyd算法),当数据点无穷多时,估计趋近于真实的测地距离;ISOMAP试图通过保持任意两点之间的测地距离来保持流行的几何结构;
LLE(Locally Linear Embedding) :
本质上是保持邻接点(局部区域)线性关系来保证整体重构的误差最小;
-
计算每一个点 X i X_{i} Xi的临近点,采用k近邻或者ε邻域;
-
计算权值 W ij W_{\text{ij}} Wij使得 X i X_{i} Xi用它的K个邻近点线性表示的误差最小,也就是通过最小化误差求出 W ij W_{\text{ij}} Wij;
min ε ( w ) = ∑ i ∥ X i − ∑ j w ij x ij ∥ 2 \min\varepsilon\left( w \right) = \ \sum_{i}^{}\left\| X_{i} - \sum_{j}^{}{w_{\text{ij}}x_{\text{ij}}} \right\|^{2} minε(w)= i∑∥∥∥∥∥Xi−j∑wijxij∥∥∥∥∥2
min ε ( Y ) = ∑ i ∥ Y i − ∑ j Y ij x ij ∥ 2 \min\varepsilon\left( Y \right) = \ \sum_{i}^{}\left\| Y_{i} - \sum_{j}^{}{Y_{\text{ij}}x_{\text{ij}}} \right\|^{2} minε(Y)= i∑∥∥∥∥∥Yi−j∑Yijxij∥∥∥∥∥2
二、算法实现:
原理:
LLE:
-
对邻域的大小k进行选择,可以选择通过与KNN一样的思想通过距离度量(欧氏距离)来选择k个最近邻;
-
寻找某个样本 x i x_{i} xi的k个最近邻之后找到其余k个最近邻之间的线性关系,也就是找到线性关系的权重系数;假设存在m个n维样本 { x 1 , x 2 , ⋅ ⋅ ⋅ ⋅ ⋅ , x m } \left\{ x_{1},x_{2}, \cdot \cdot \cdot \cdot \cdot {,x}_{m} \right\} {x1,x2,⋅⋅⋅⋅⋅,xm}
J ( w ) = ∑ i ∥ X i − ∑ j ∈ Q w ij x ij ∥ 2 其 中 , ∑ j ∈ Q w ij = 1 J\left( w \right) = \sum_{i}^{}\left\| X_{i} - \sum_{j \in Q}^{}{w_{\text{ij}}x_{\text{ij}}} \right\|^{2}其中,\sum_{j \in Q}^{}w_{\text{ij}} = \ 1 J(w)=i∑∥∥∥∥∥∥Xi−j∈Q∑wijxij∥∥∥∥∥∥2其中,j∈Q∑wij= 1
- 对J(w)进行矩阵化:
J ( w ) = ∑ i = 1 m W i T ( x i − x j ) ( x i − x j ) T W i J\left( w \right) = \ \sum_{i = 1}^{m}{W_{i}^{T}\left( x_{i} - x_{j} \right)\left( x_{i} - x_{j} \right)^{T}W_{i}} J(w)= i=1∑mWiT(xi−xj)(xi−xj)TWi
- 令矩阵 Z i = ( x i − x j ) ( x i − x j ) T Z_{i} = \ \left( x_{i} - x_{j} \right)\left( x_{i} - x_{j} \right)^{T} Zi= (xi−xj)(xi−xj)T,用拉格朗日乘法将其合为一个目标;
L ( w ) = ∑ i = 1 k W i T Z i W i + λ ( W i T 1 k − 1 ) L\left( w \right) = \ \sum_{i = 1}^{k}{W_{i}^{T}Z_{i}W_{i}} + \ \lambda\left( W_{i}^{T}1_{k} - 1 \right) L(w)= i=1∑kWiTZiWi+ λ(WiT1k−1)
- 对w求导并令其值为0,可算得
W i = λ ′ Z i − 1 1 k W_{i} = {\lambda\ }^{'}Z_{i}^{- 1}1_{k} Wi=λ ′Zi−11k
- 带入 W i T 1 k = 1 最 终 W_{i}^{T}1_{k} = 1\ 最终 WiT1k=1 最终化简得到
W i = Z i − 1 1 k 1 k T Z i − 1 1 k W_{i} = \ \frac{Z_{i}^{- 1}1_{k}}{{1_{k}^{T}Z}_{i}^{- 1}1_{k}} Wi= 1kTZi−11kZi−11k
- 这些权重系数对应的线性关系在降维后的低维一样得到保持。假设我们的n维样本集
{
x
1
,
x
2
,
⋅
⋅
⋅
⋅
⋅
,
x
m
}
\left\{ x_{1},x_{2}, \cdot \cdot \cdot \cdot \cdot {,x}_{m} \right\}
{x1,x2,⋅⋅⋅⋅⋅,xm}在低维的d维度对应投影为
{
y
1
,
y
2
,
⋅
⋅
⋅
⋅
⋅
,
y
m
}
\left\{ y_{1},y_{2}, \cdot \cdot \cdot \cdot \cdot {,y}_{m} \right\}
{y1,y2,⋅⋅⋅⋅⋅,ym},
则我们希望保持线性关系,最小化函数J(y):
J ( y ) = ∑ i ∥ y i − ∑ j ∈ Q w ij y ij ∥ 2 , 其 中 Y Y T = m I , ∑ i = 1 m y i = 0 J\left( y \right) = \sum_{i}^{}\left\| y_{i} - \sum_{j \in Q}^{}{w_{\text{ij}}y_{\text{ij}}} \right\|^{2},其中YY^{T} = mI,\sum_{i = 1}^{m}y_{i} = 0 J(y)=i∑∥∥∥∥∥∥yi−j∈Q∑wijyij∥∥∥∥∥∥2,其中YYT=mI,i=1∑myi=0
= ∑ i ∥ YI i − YW i ∥ 2 \ \ \ \ = \sum_{i}^{}\left\| \text{YI}_{i} - \text{YW}_{i} \right\|^{2} =i∑∥YIi−YWi∥2
= tr ( Y ( I − W ) ( I − W ) T Y T ) \text{\ \ \ \ } = \text{tr}\left( Y\left( I - W \right)\left( I - W \right)^{T}Y^{T} \right) =tr(Y(I−W)(I−W)TYT)
- 令矩阵 M = ( I − W ) ( I − W ) T M = \ \left( I - W \right)\left( I - W \right)^{T} M= (I−W)(I−W)T,用拉格朗日乘法将其合为一个目标
L ( Y ) = tr ( YM Y T + λ ( Y Y T − m I ) ) L\left( Y \right) = \ \text{tr}\left( \text{YM}Y^{T} + \lambda\left( YY^{T} - mI \right) \right) L(Y)= tr(YMYT+λ(YYT−mI))
解 得 : M Y T = λ ′ Y T 解得:\text{\ \ \ \ }MY^{T} = \lambda^{'}Y^{T} 解得: MYT=λ′YT
-
求出矩阵M最小的d个特征值所对应的d个特征向量组成的矩阵;
-
利用特征向量做对原始数据做相应的线性变换;
ISOMAP:
(注意:ISOMAP是数据越密集计算映射也就越准确,但由于其中Floyd算法复杂度高,计算机运算速度有限需要对总生成点数和临近点数不断调整以此到达最佳展示效果,而且由于初始点是随机生成的会偶尔出现临近点数量不够导致部分点路径不可达,从而无法计算特征值,程序报错;因此出错只需增大随机点数量或者增加临近点个数即可)
-
对于每一个样本 X i X_{i} Xi寻找其k个临近点,将 X i X_{i} Xi与k个邻近点之间的距离设置为欧氏距离,其他点距离设置为无穷大;
-
调用最短路径算法计算任意两样本之间的距离dist( X i X_{i} Xi, X j X_{j} Xj);
-
将dist( X i X_{i} Xi, X j X_{j} Xj)作为MDS算法的输入,获得到低维空间的映射;
关于MDS算法步骤及推导:
- 假定m个样本在原始空间的距离矩阵为
D
∈
R
m
∗
m
D \in R^{m*m}
D∈Rm∗m
(也就是MDS算法的输入),我们的目标是获得样本d’维空间的表示
Z ∈ R ( d ′ ∗ m ) Z \in R^{\left( d^{'}*m \right)} Z∈R(d′∗m) ,且任意两个样本在 d ′ d^{'} d′
维空间的欧氏距离等于原始空间中的距离。
令 B = Z T Z ∈ R m ∗ m B = \ Z^{T}Z \in R^{m*m} B= ZTZ∈Rm∗m,其中B为降维后的样本的内积矩阵, b ij = z i T z j b_{\text{ij}} = z_{i}^{T}z_{j} bij=ziTzj,则
d ist ij 2 = ∥ z i ∥ 2 + ∥ z j ∥ 2 − x z i T z j {d\text{ist}}_{\text{ij}}^{2}\ = \ \left\| z_{i} \right\|^{2} + \left\| z_{j} \right\|^{2} - x{z_{i}}^{T}z_{j} distij2 = ∥zi∥2+∥zj∥2−xziTzj
= b ii + b jj − 2 b i j ① = b_{\text{ii}} + b_{\text{jj}} - 2b_{ij}① =bii+bjj−2bij①
为了便于推导,领样本降维后的Z去中心化, ∑ i = 1 m z i = 0 \sum_{i = 1}^{m}z_{i} = 0 ∑i=1mzi=0显然B行列之和为0,则
∑ i = 1 m dist ij 2 = tr ( B ) + m b ij ② \sum_{i = 1}^{m}\text{dist}_{\text{ij}}^{2} = \text{tr}\left( B \right) + mb_{\text{ij}}② i=1∑mdistij2=tr(B)+mbij②
∑ i = 1 m dist ij 2 = t r ( B ) + m b ij ③ \sum_{i = 1}^{m}\text{dist}_{\text{ij}}^{2} = tr\left( B \right) + mb_{\text{ij}}③ i=1∑mdistij2=tr(B)+mbij③
∑ i = 1 m ∑ j = 1 m dist ij 2 = 2 m t r ( B ) ④ \sum_{i = 1}^{m}{\sum_{j = 1}^{m}\text{dist}_{\text{ij}}^{2}} = 2m\ tr\left( B \right)④ i=1∑mj=1∑mdistij2=2m tr(B)④
其中 tr ( ⋅ ) \text{tr}\left( \cdot \right) tr(⋅)表示矩阵的迹, tr ( B ) = ∑ i = 1 m ∥ z ∥ i 2 \ \text{tr}\left( B \right) = \sum_{i = 1}^{m}\left\| z \right\|_{i}^{2} tr(B)=∑i=1m∥z∥i2,令
d ist i ⋅ 2 = 1 m ∑ j = 1 m dist ij 2 ⑤ {d\text{ist}}_{i \cdot}^{2}\ = \frac{1}{m}\sum_{j = 1}^{m}\text{dist}_{\text{ij}}^{2}⑤ disti⋅2 =m1j=1∑mdistij2⑤
d ist ⋅ j 2 = 1 m ∑ i = 1 m dist ij 2 ⑥ {d\text{ist}}_{\cdot j}^{2}\ = \frac{1}{m}\sum_{i = 1}^{m}{\text{dist}_{\text{ij}}^{2}⑥} dist⋅j2 =m1i=1∑mdistij2⑥
d ist i j 2 = 1 m 2 ∑ i = 1 m ∑ j m dist ij 2 ⑦ {d\text{ist}}_{ij}^{2}\ = \frac{1}{m^{2}}\sum_{i = 1}^{m}{\sum_{j}^{m}\text{dist}_{\text{ij}}^{2}}⑦ distij2 =m21i=1∑mj∑mdistij2⑦
由式① ~⑦可知
b i j = − 1 2 ( d ist ij 2 − d ist i ⋅ 2 − d ist ⋅ j 2 + d ist ij 2 ) b_{ij} = - \frac{1}{2}({d\text{ist}}_{\text{ij}}^{2} - {d\text{ist}}_{i \cdot}^{2} - {d\text{ist}}_{\cdot j}^{2} + {d\text{ist}}_{\text{ij}}^{2}) bij=−21(distij2−disti⋅2−dist⋅j2+distij2)
以此得到保持D不变的内积矩阵B;
- 对矩阵B做特征值分解,求出对应的特征向量,做出从高维到低维的映射;
源码:
LLE
function [ Y ] = LLE( X,K,d )
%LLE 此处显示有关此函数的摘要
% X为待处理矩阵
% k代表邻接点数
% d代表最终维度数
%LLE代码
[D,N] = size(X);
%计算距离并且找出邻近点
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),:);
%计算权重W
if(K>D)
tol=1e-3;
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
% M=(I-W)'(I-W)
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
% 计算对应低维空间的象
options.disp = 0;
options.isreal = 1;
options.issym = 1;
[Y,eigenvals] = eigs(M,d+1,0,options);
Y = Y(:,2:d+1)'*sqrt(N);
end
%==========================================================
%测试调用代码
N=2000;
K=16;%领域点个数
d=2; %最大嵌入维数
% 绘制原始图像
tt0 = (3*pi/2)*(1+2*[0:0.02:1]);
hh = [0:0.125:1]*30;
xx = (tt0.*cos(tt0))'*ones(size(hh));
yy = ones(size(tt0))'*hh;
zz = (tt0.*sin(tt0))'*ones(size(hh));
cc = tt0'*ones(size(hh));
figure(1)
surf(xx,yy,zz,cc);
title('原始真实图像')
% 生成实例数据
tt = (3*pi/2)*(1+2*rand(1,N));
height = 21*rand(1,N);
X = [tt.*cos(tt); height; tt.*sin(tt)];
% 实例数据可视化
figure(2)
scatter3(X(1,:),X(2,:),X(3,:),12,tt,'+');
title('随机生成的符合原始图像的数据')
% LLE嵌入
Y=LLE(X,K,d);
disp(['邻近点个数为' num2str(K)])
disp(['最大嵌入维度' num2str(d)])
%嵌入的低维数据可视化
figure(3); cla;
scatter(Y(1,:),Y(2,:),12,tt,'+');
title('LLE降维后的结果')
ISOMAP
function [ V1 ] = ISOMAP1( X,k,d )
%UNTITLED 此处显示有关此函数的摘要
% 此处显示详细说明
%计算距离矩阵个
X=X';
[m,n]=size(X);
D=zeros(m,m);
for i=1:m
for j=i:m
D(i,j)=norm(X(i,:)-X(j,:));
D(j,i)=D(i,j);
end
end
%计算矩阵中每行前k个值的位置并赋值(先按大小排列)
W1=zeros(m,m);
for i=1:m
A=D(i,:);
t=sort(A(:));%对每行进行排序后构成一个从小到大有序的列向量
[row,col]=find(A<=t(k),k);%找出每行前K个最小数的位置
for j=1:k
c=col(1,j);
W1(i,c)=D(i,c); %W1(i,c)=1;%给k近邻赋值为距离
end
end
for i=1:m
for j=1:m
if W1(i,j)==0&i~=j
W1(i,j)=inf;
end
end
end
%计算测地距离o就是测地距离
% [dist,mypath,o]=myfloyd(W1,10,20);
% figure
% [col,rol]=size(mypath);
% X1=[];
% for i=1:rol
% ding=mypath(1,i);
% X1=[X1;X(ding,:)];
% end
% plot3(X(:,1),X(:,2),X(:,3),'.')
% hold on
% plot3(X1(:,1),X1(:,2),X1(:,3),'o-r')
N = length(W1);
for k=1:N
for i=1:N
for j=1:N
if W1(i,j)>W1(i,k)+W1(k,j)
W1(i,j)=W1(i,k)+W1(k,j);
end
end
end
end
%mds算法降维
d=W1;
n=size(d,1);
B=zeros(n,n);
for i=1:n
for j=1:n
B(i,j)=-0.5*(d(i,j)^2 -1/n*d(i,:)*d(i,:)' -1/n*d(:,j)'*d(:,j) +1/n^2*sum(sum(d.^2)));
end
end
[V,D]=eig(B);
A=sort(diag(D),'descend');%把对角矩阵变成列向量
% V1=[V(:,1),V(:,2)]*[A(1,1),0;0,A(2,1)];
V1=V(:,1:2)*D(1:2,1:2).^(1/2);
V1 = V1';
end
%==================================================================================================
%测试调用代码
N=500;
K=8;%领域点个数
d=2; %最大嵌入维数
% 绘制原始图像
tt0 = (3*pi/2)*(1+2*[0:0.02:1]);
hh = [0:0.125:1]*30;
xx = (tt0.*cos(tt0))'*ones(size(hh));
yy = ones(size(tt0))'*hh;
zz = (tt0.*sin(tt0))'*ones(size(hh));
cc = tt0'*ones(size(hh));
figure(1)
surf(xx,yy,zz,cc);
title('原始真实图像')
% 生成实例数据
tt = (3*pi/2)*(1+2*rand(1,N));
height = 21*rand(1,N);
X = [tt.*cos(tt); height; tt.*sin(tt)];
% 实例数据可视化
figure(2)
scatter3(X(1,:),X(2,:),X(3,:),12,tt,'+');
title('随机生成的符合原始图像的数据')
% ISOMAP嵌入
Y=ISOMAP(X,K,d);
disp(['邻近点个数为' num2str(K)])
disp(['最大嵌入维度' num2str(d)])
%嵌入的低维数据可视化
figure(3); cla;
scatter(Y(1,:),Y(2,:),12,tt,'+');
title('ISOMAP降维后的结果')
实验结果:
LLE
ISOMAP
(注意:ISOMAP是数据越密集计算映射也就越准确,但由于其中Floyd算法复杂度高,计算机运算速度有限需要对总生成点数和临近点数不断调整以此到达最佳展示效果,而且由于初始点是随机生成的会偶尔出现临近点数量不够导致部分点路径不可达,从而无法计算特征值,程序报错;因此出错只需增大随机点数量或者增加临近点个数即可)