一、简述:
试图保持邻域内样本之间的线性关系。
如上图,
X
n
×
p
X_{n\times p}
Xn×p中样本点
x
(
i
)
x_{(i)}
x(i)的坐标可通过邻域样本
x
(
j
)
x_{(j)}
x(j),
x
(
k
)
x_{(k)}
x(k),
x
(
l
)
x_{(l)}
x(l)线性组合重构出来,计算出取值之后再通过算出的权重矩阵重构
X
n
×
d
X_{n\times d}
Xn×d.
二、算法
(1)计算线性重构的系数矩阵
原数据
X
n
×
p
=
(
x
(
1
)
T
x
(
2
)
T
.
.
x
(
n
)
T
)
X_{n\times p}=\begin{pmatrix} x_{(1)}^T \\x_{(2)}^T\\.\\.\\ x_{(n)}^T \end{pmatrix}
Xn×p=⎝⎜⎜⎜⎜⎛x(1)Tx(2)T..x(n)T⎠⎟⎟⎟⎟⎞,
x
(
i
)
=
(
x
i
1
,
x
i
2
,
.
.
.
,
x
i
p
)
T
x_{(i)}=(x_{i1},x_{i2},...,x_{ip})^T
x(i)=(xi1,xi2,...,xip)T.假设
Q
(
i
)
Q_{(i)}
Q(i)是样本
x
(
i
)
x_{(i)}
x(i)的
k
k
k近邻,用
x
(
1
)
,
x
(
2
)
,
.
.
.
x
(
k
)
x_{(1)},x_{(2)},...x_{(k)}
x(1),x(2),...x(k)来表示,则
x
(
i
)
=
∑
j
∈
Q
(
i
)
w
i
j
x
(
j
)
=
∑
j
=
1
k
w
i
j
x
(
j
)
x_{(i)}=\sum_{j \in Q_{(i)}} w_{ij} x_{(j)}=\sum_{j =1}^k w_{ij} x_{(j)}
x(i)=j∈Q(i)∑wijx(j)=j=1∑kwijx(j)其中
∑
j
∈
Q
(
i
)
w
i
j
=
1
\sum_{j \in Q_{(i)}} w_{ij}=1
∑j∈Q(i)wij=1。使用均方差作为回归问题的损失函数,即求:
arg
min
w
i
∑
i
=
1
n
∣
∣
x
(
i
)
−
∑
j
∈
Q
(
i
)
w
i
j
x
(
j
)
∣
∣
2
2
(1)
\arg \min \limits_{w_i}\sum_{i=1}^n ||x_{(i)}-\sum_{j \in Q_{(i)}} w_{ij} x_{(j)}||_2^2 \tag{1}
argwimini=1∑n∣∣x(i)−j∈Q(i)∑wijx(j)∣∣22(1)
s
.
t
.
∑
j
∈
Q
(
i
)
w
i
j
=
1
s.t.\sum_{j \in Q_{(i)}} w_{ij}=1
s.t.j∈Q(i)∑wij=1将
(
1
)
(1)
(1)式矩阵化,
(
1
)
=
arg
min
w
i
∑
i
=
1
n
∣
∣
∑
j
∈
Q
(
i
)
w
i
j
(
x
(
i
)
−
x
(
j
)
)
∣
∣
2
2
=
arg
min
w
i
∑
i
=
1
n
∣
∣
(
(
x
(
i
)
−
x
(
1
)
)
,
(
x
(
i
)
−
x
(
1
)
)
,
.
.
.
(
x
(
i
)
−
x
(
k
)
)
)
(
w
i
1
w
i
2
.
.
w
i
k
)
∣
∣
2
2
=
arg
min
w
i
∑
i
=
1
n
∣
∣
(
X
i
−
N
i
)
w
i
∣
∣
2
2
=
arg
min
w
i
∑
i
=
1
n
[
w
i
T
(
X
i
−
N
i
)
T
(
X
i
−
N
i
)
w
i
]
=
arg
min
w
i
∑
i
=
1
n
[
w
i
T
S
i
w
i
]
\begin{aligned}(1)&=\arg \min \limits_{w_i}\sum_{i=1}^n ||\sum_{j \in Q_{(i)}} w_{ij} (x_{(i)}-x_{(j)})||_2^2 \\ &=\arg \min \limits_{w_i}\sum_{i=1}^n ||\begin{pmatrix}(x_{(i)}-x_{(1)}),(x_{(i)}-x_{(1)}),...(x_{(i)}-x_{(k)})\end{pmatrix} \begin{pmatrix} w_{i1} \\w_{i2}\\.\\.\\ w_{ik} \end{pmatrix}||_2^2 \\ &=\arg \min \limits_{w_i}\sum_{i=1}^n ||(X_i-N_i)w_i||_2^2 \\ &=\arg \min \limits_{w_i}\sum_{i=1}^n [w_i^T(X_i-N_i)^T(X_i-N_i)w_i] \\ &=\arg \min \limits_{w_i}\sum_{i=1}^n [w_i^TS_iw_i] \end{aligned}
(1)=argwimini=1∑n∣∣j∈Q(i)∑wij(x(i)−x(j))∣∣22=argwimini=1∑n∣∣((x(i)−x(1)),(x(i)−x(1)),...(x(i)−x(k)))⎝⎜⎜⎜⎜⎛wi1wi2..wik⎠⎟⎟⎟⎟⎞∣∣22=argwimini=1∑n∣∣(Xi−Ni)wi∣∣22=argwimini=1∑n[wiT(Xi−Ni)T(Xi−Ni)wi]=argwimini=1∑n[wiTSiwi]其中,
X
i
=
(
x
(
i
)
,
x
(
i
)
,
.
.
.
x
(
i
)
)
,
N
i
=
(
x
(
1
)
,
x
(
2
)
,
.
.
.
x
(
k
)
)
,
w
i
=
(
w
i
1
w
i
2
.
.
w
i
k
)
,
S
i
=
(
X
i
−
N
i
)
T
(
X
i
−
N
i
)
X_i=(x_{(i)},x_{(i)},...x_{(i)}),N_i=(x_{(1)},x_{(2)},...x_{(k)}),w_i=\begin{pmatrix} w_{i1} \\w_{i2}\\.\\.\\ w_{ik}\end{pmatrix},S_i=(X_i-N_i)^T(X_i-N_i)
Xi=(x(i),x(i),...x(i)),Ni=(x(1),x(2),...x(k)),wi=⎝⎜⎜⎜⎜⎛wi1wi2..wik⎠⎟⎟⎟⎟⎞,Si=(Xi−Ni)T(Xi−Ni),使用矩阵和拉格朗日乘子法求解该问题,令
Y
(
w
i
)
=
∑
i
=
1
n
[
w
i
T
S
i
w
i
]
+
λ
(
1
k
T
w
i
−
1
)
∂
Y
(
w
i
)
∂
w
i
=
2
w
i
T
S
i
+
λ
1
k
T
=
0
分
子
布
局
求
导
⟹
2
S
i
w
i
+
λ
1
k
=
0
⟹
w
i
=
−
1
2
λ
S
i
−
1
1
k
Y(w_i)=\sum_{i=1}^n [w_i^TS_iw_i]+\lambda(1_k^Tw_i-1) \\ \frac {\partial Y(w_i)}{\partial{w_i}} =2w_i^TS_i+\lambda 1_k^T=0\,\,\,\,\,\,分子布局求导\Longrightarrow 2S_iw_i+\lambda 1_k=0\Longrightarrow w_i=-\frac{1}{2}\lambda S_i^{-1}1_k
Y(wi)=i=1∑n[wiTSiwi]+λ(1kTwi−1)∂wi∂Y(wi)=2wiTSi+λ1kT=0分子布局求导⟹2Siwi+λ1k=0⟹wi=−21λSi−11k由
1
k
T
w
i
=
1
⟹
−
1
2
λ
=
1
k
S
i
1
k
T
⟹
w
i
=
S
i
−
1
1
k
1
k
T
S
i
1
k
1_k^Tw_i=1 \Longrightarrow -\frac{1}{2}\lambda=1_kS_i1_k^T \Longrightarrow w_i=\frac{S_i^{-1}1_k}{1_k^TS_i1_k}
1kTwi=1⟹−21λ=1kSi1kT⟹wi=1kTSi1kSi−11k
注:计算得到的
w
i
w_i
wi为
k
×
1
k\times 1
k×1向量,把非近邻的值设为0,得到
n
×
1
n\times 1
n×1维向量。
(2)确定低维空间的坐标
X n × p = ( x ( 1 ) T x ( 2 ) T . . x ( n ) T ) ⟹ Y n × d = ( x ( 1 ) T x ( 2 ) T . . x ( n ) T ) X_{n\times p}=\begin{pmatrix} x_{(1)}^T \\x_{(2)}^T\\.\\.\\ x_{(n)}^T \end{pmatrix} \Longrightarrow Y_{n\times d}=\begin{pmatrix} x_{(1)}^T \\x_{(2)}^T\\.\\.\\ x_{(n)}^T \end{pmatrix} Xn×p=⎝⎜⎜⎜⎜⎛x(1)Tx(2)T..x(n)T⎠⎟⎟⎟⎟⎞⟹Yn×d=⎝⎜⎜⎜⎜⎛x(1)Tx(2)T..x(n)T⎠⎟⎟⎟⎟⎞求解: arg min Y ∑ i = 1 n ∣ ∣ y ( i ) − ∑ j = 1 n w i j y ( j ) ∣ ∣ 2 2 (2) \arg \min \limits_{Y} \sum_{i=1}^n ||y_{(i)}-\sum_{j=1}^n w_{ij} y_{(j)}||_2^2 \tag{2} argYmini=1∑n∣∣y(i)−j=1∑nwijy(j)∣∣22(2) s . t . Y T Y = m I s.t.\,\,\,\,Y^TY=mI s.t.YTY=mI ( 2 ) = arg min Y ∑ i = 1 n ∣ ∣ Y T I i − Y T w i ∣ ∣ 2 2 = arg min Y ∑ i = 1 n ∣ ∣ Y T ( I i − w i ) ∣ ∣ 2 2 = arg min Y t r ( Y T ( I − W ) ( I − W ) T Y ) = arg min Y t r ( Y T M Y ) \begin{aligned}(2)&=\arg \min \limits_{Y} \sum_{i=1}^n ||Y^TI_i-Y^Tw_i||_2^2 \\ &=\arg \min \limits_{Y} \sum_{i=1}^n ||Y^T(I_i-w_i)||_2^2 \\ &=\arg \min \limits_{Y}tr(Y^T(I-W)(I-W)^TY) \\&=\arg \min \limits_{Y}tr(Y^TMY) \end{aligned} (2)=argYmini=1∑n∣∣YTIi−YTwi∣∣22=argYmini=1∑n∣∣YT(Ii−wi)∣∣22=argYmintr(YT(I−W)(I−W)TY)=argYmintr(YTMY)补充: 若 A = ( a 1 , a 2 , . . . , a N ) , 则 ∑ i ( a i ) 2 = ∑ i a i T a i = t r ( A A T ) 若A=(a_1,a_2,...,a_N),则\sum_i(a_i)^2=\sum_ia_i^Ta_i=tr(AA^T) 若A=(a1,a2,...,aN),则∑i(ai)2=∑iaiTai=tr(AAT).使用矩阵和拉格朗日乘子法求解该问题,令 J ( y ) = t r ( Y T M Y ) + λ ( Y T Y − m I ) ∂ J ( Y ) ∂ Y = 2 M Y + 2 λ Y = 0 J(y)=tr(Y^TMY)+\lambda (Y^TY-mI) \\ \frac {\partial J(Y)}{\partial{Y}}=2MY+2\lambda Y=0 J(y)=tr(YTMY)+λ(YTY−mI)∂Y∂J(Y)=2MY+2λY=0即 M Y = λ ′ Y MY=\lambda^{'}Y MY=λ′Y, Y Y Y为 M M M对应的特征向量矩阵。
(3)注意事项
取M的除0之外的最小特征值所对应的特征向量,第 2 , 3 , . . . d + 1 2,3,...d+1 2,3,...d+1个特征值对应 Y = ( y 2 , y 3 , . . . , y d + 1 ) Y=(y_2,y_3,...,y_d+1) Y=(y2,y3,...,yd+1).
- 为什么不能为0?M的特征值为0不能反应数据特征,对应的特征向量全为1. W T e = e → ( W T − I ) e = 0 → W T − I = 0 → W − I = 0 → ( I − W ) ( I − W ) T e = 0 e W^Te=e\rightarrow (W^T-I)e=0\rightarrow W^T-I=0 \rightarrow W-I=0\rightarrow (I-W)(I-W)^Te=0e WTe=e→(WT−I)e=0→WT−I=0→W−I=0→(I−W)(I−W)Te=0e,M的最小特征值为0,需要排除。
- 为什么取最小?Y是实对称矩阵的特征向量,故两两正交, Y T Y = m I ( m 看 作 1 ) Y^TY=mI(m看作1) YTY=mI(m看作1),则Y为正交矩阵, Y T = Y − 1 Y^T=Y^{-1} YT=Y−1。求 arg min Y t r ( Y T M Y ) \arg \min \limits_{Y}tr(Y^TMY) argYmintr(YTMY)相当于求特征值的和最小。
三、python实现
(1)自己写的
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from mpl_toolkits.mplot3d import Axes3D
def load_data():
swiss_roll =datasets.make_swiss_roll(n_samples=1000)
return swiss_roll[0],np.floor(swiss_roll[1])
#my_LLE任意计算两点之间距离
def cal_pairwise_dist(x): #任意两点之间距离的平方
'''计算pairwise 距离, x是matrix
(a-b)^2 = a^2 + b^2 - 2*a*b
'''
sum_x = np.sum(np.square(x), 1) #square计算各元素平方 1 按行相加 0按列相加
dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
#返回任意两个点之间距离的平方
return dist
#my_LLe近邻
def floyd(D,n_neighbors): #k近邻得到距离最近的指标
D_arg = np.argsort(D,axis=1)#返回从小到大的索引值,每一列进行排序
print(type(D_arg))
return D_arg[:,1:n_neighbors+1]
def linear_weight(X,X_floyd):
n1,n2=X_floyd.shape
n3,n4=X.shape
D1 = np.ones((n4,n2))
w=np.zeros((n1,n1))
one_k=np.ones((n2,1))
for i in range(n1):
for j in range(n2):
D1[:,j]=np.transpose(X[i,:]-X[X_floyd[i,j],:])
a=np.dot(np.transpose(one_k),np.linalg.inv(np.dot(np.transpose(D1),D1))).dot(one_k)
W_I=np.linalg.inv(np.dot(np.transpose(D1),D1)).dot(one_k)/a
for k in range(n2):
w[X_floyd[i,k],i]=W_I[k]
return w
def reconstruction(w,d):
n1,n2=w.shape
danwei=np.eye(n1)
M=(danwei-w).dot(np.transpose(danwei-w))
eig_val, eig_vec = np.linalg.eig(M) #计算矩阵特征值和特征向量
eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(n2)] #abs绝对值 把特征值对应的特征向量放一块
# sort eig_vec based on eig_val from highest to lowest
eig_pairs.sort(key=lambda x:x[0]) #升序排列
# select the top d eig_vec
feature=np.array([ele[1] for ele in eig_pairs[1:d+1]]) #选出d个特征值对应的向量,此时为特征向量为行向量
#get new data
new_data=np.transpose(feature) #向量显示为列向量
return new_data
def my_LLE(X,Y,Neighbors):
D = cal_pairwise_dist(X)
D[D < 0] = 0
D = D**0.5 #距离矩阵
X_floyd=floyd(D,Neighbors)
w=linear_weight(X,X_floyd)
new_data=reconstruction(w,2)
return new_data
X,Y=load_data()
fig = plt.figure('data')
ax = Axes3D(fig)
ax.scatter(X[:, 0], X[:, 1], X[:, 2],marker='o',c=Y)
fig=plt.figure("my_LLE",figsize=(9, 9))
Neighbors=[1,2,3,4,5,15,30,100,Y.size-1]
for i,k in enumerate(Neighbors): #i位置,k代表的值
lle=my_LLE(X,Y,k)
ax=fig.add_subplot(3,3,i+1)
ax.scatter(lle[:,0],lle[:,1],marker='o',c=Y,alpha=0.5) #c代表颜色
ax.set_title("k = %d"%k)
plt.xticks(fontsize=10, color="darkorange")
plt.yticks(fontsize=10, color="darkorange")
#plt.suptitle("LLE") #总图标题
plt.show()
(2)使用manifold.LocallyLinearEmbedding
import numpy as np
import operator
import matplotlib.pyplot as plt
from sklearn import datasets,decomposition,manifold
from itertools import cycle
from mpl_toolkits.mplot3d import Axes3D
def load_data():
swiss_roll =datasets.make_swiss_roll(n_samples=1000)
return swiss_roll[0],np.floor(swiss_roll[1])
def LLE_components(*data):
X,Y=data
for n in [3,2,1]:
lle=manifold.LocallyLinearEmbedding(n_components=n) #LLE降维,n_components降的维数 #n_neighbors近邻数
lle.fit(X)
print("n = %d 重建误差:"%n,lle.reconstruction_error_)
def LLE_neighbors(*data):
X,Y=data
Neighbors=[1,2,3,4,5,15,30,100,Y.size-1]
fig=plt.figure("LLE",figsize=(9, 9))
for i,k in enumerate(Neighbors): #i位置,k代表的值
lle=manifold.LocallyLinearEmbedding(n_components=2,n_neighbors=k,eigen_solver='dense')
X_r=lle.fit_transform(X)
ax=fig.add_subplot(3,3,i+1)
ax.scatter(X_r[:,0],X_r[:,1],marker='o',c=Y,alpha=0.5) #c代表颜色
ax.set_title("k = %d"%k)
plt.xticks(fontsize=10, color="darkorange")
plt.yticks(fontsize=10, color="darkorange")
plt.suptitle("LLE")
plt.show()
X,Y=load_data()
fig = plt.figure('data')
ax = Axes3D(fig)
ax.scatter(X[:, 0], X[:, 1], X[:, 2],marker='o',c=Y)
LLE_components(X,Y)
LLE_neighbors(X,Y)
四、附录
- 1)n_neighbors:即我们搜索样本的近邻的个数,默认是5。 n_neighbors个数越大,则建立样本局部关系的时间会越大,也就意味着算法的复杂度会增加。当然n_neighbors个数越大,则降维后样本的局部关系会保持的更好。在下一节我们可以通过具体的例子看出这一点。一般来说,如果算法运行时间可以接受,我们可以尽量选择一个比较大一些的n_neighbors
- 2)n_components:即我们降维到的维数。如果我们降维的目的是可视化,则一般可以选择2-5维。
- 3)reg :正则化系数,在n_neighbors大于n_components时,即近邻数大于降维的维数时,由于我们的样本权重矩阵不是满秩的,LLE通过正则化来解决这个问题。默认是0.001。一般不用管这个参数。当近邻数远远的大于降维到的维数时可以考虑适当增大这个参数。
- 4)eigen_solver:特征分解的方法。有‘arpack’和‘dense’两者算法选择。当然也可以选择’auto’让scikit-learn自己选择一个合适的算法。‘arpack’和‘dense’的主要区别是‘dense’一般适合于非稀疏的矩阵分解。而‘arpack’虽然可以适应稀疏和非稀疏的矩阵分解,但在稀疏矩阵分解时会有更好算法速度。当然由于它使用一些随机思想,所以它的解可能不稳定,一般需要多选几组随机种子来尝试。
- 5)method: 即LLE的具体算法。LocallyLinearEmbedding支持4种LLE算法,分别是’standard’对应我们标准的LLE算法,'hessian’对应原理篇讲到的HLLE算法,‘modified’对应原理篇讲到的MLLE算法,‘ltsa’对应原理篇讲到的LTSA算法。默认是’standard’。一般来说HLLE/MLLE/LTSA算法在同样的近邻数n_neighbors情况下,运行时间会比标准的LLE长,当然降维的效果会稍微好一些。如果你对降维后的数据局部效果很在意,那么可以考虑使用HLLE/MLLE/LTSA或者增大n_neighbors,否则标准的LLE就可以了。需要注意的是使用MLLE要求n_neighbors > n_components,而使用HLLE要求n_neighbors > n_components * (n_components + 3) / 2
- 6)neighbors_algorithm:这个是k近邻的搜索方法,和KNN算法的使用的搜索方法一样。算法一共有三种,第一种是蛮力实现,第二种是KD树实现,第三种是球树实现。这三种方法在K近邻法(KNN)原理小结中都有讲述,如果不熟悉可以去复习下。对于这个参数,一共有4种可选输入,‘brute’对应第一种蛮力实现,‘kd_tree’对应第二种KD树实现,‘ball_tree’对应第三种的球树实现, ‘auto’则会在上面三种算法中做权衡,选择一个拟合最好的最优算法。需要注意的是,如果输入样本特征是稀疏的时候,无论我们选择哪种算法,最后scikit-learn都会去用蛮力实现‘brute’。个人的经验,如果样本少特征也少,使用默认的 ‘auto’就够了。 如果数据量很大或者特征也很多,用"auto"建树时间会很长,效率不高,建议选择KD树实现‘kd_tree’,此时如果发现‘kd_tree’速度比较慢或者已经知道样本分布不是很均匀时,可以尝试用‘ball_tree’。而如果输入样本是稀疏的,无论你选择哪个算法最后实际运行的都是‘brute’。