【ML】LinearDiscriminantAnalysis线性判别分析

• 线性判别分析LDA于1936年由Fisher提出&＃xff0c;其基本思想是将给定的训练集设法投影到一条直线 $l$ 上&＃xff0c;并且很自然地&＃xff0c;我们力求投影点集合中同类的点相距越近越好&＃xff0c;异类的点相距越远越好。
• 给定一个新样本&＃xff0c;采用同样的投影方法&＃xff0c;将其投影到 $l$ 上&＃xff0c;再根据投影点的位置来确定其类别。
  

• 给定训练集 D&＃61;{(xi,yi)}i&＃61;1m&＃xff0c;yi∈{0,1,..,K}D&＃61;\{(x_i,y_i)\}_{i&＃61;1}^m&＃xff0c;y_i∈\{0,1,..,K\}D&＃61;{(xi​,yi​)}i&＃61;1m​&＃xff0c;yi​∈{0,1,..,K} &＃xff0c;记 μi,Σi\mu_i,\Sigma_iμi​,Σi​分别为第 $k$ 类样本的均值向量和协方差矩阵。不妨假定将数据集 $D$ 投影到直线 y&＃61;wTxy&＃61;w^Txy&＃61;wTx 上&＃xff0c;截距项只是对直线产生平移效果&＃xff0c;对于第一部分我们所关心的同类点、异类点的距离并无影响。
• 先从二分类问题入手&＃xff0c;即 yi∈{0,1}y_i∈\{0,1\}yi​∈{0,1}&＃xff0c;那么投影后&＃xff0c;我们得到两类样本点的均值向量在直线上的投影分别为 【wTμ0、wTμ1w^T\mu_0、w^T\mu_1wTμ0​、wTμ1​】&＃xff0c;两类样本投影点的协方差分别为 【wTΣ0w、wTΣ1ww^T\Sigma_0w、w^T\Sigma_1wwTΣ0​w、wTΣ1​w】。此时所有的样本投影点点都位于直线上&＃xff0c;易知投影点均值以及投影点协方差皆为实数。
• 考虑第一部分给出的基本思想&＃xff0c;同类点尽可能靠近&＃xff0c;异类点尽可能远离&＃xff0c;即我们希望两类样本投影点的协方差都尽可能小&＃xff0c;而两个均值要差距尽可能大&＃xff0c;形式化的描述如下&＃xff1a;
  
• 其中 $J$ 即为我们希望最大化的目标&＃xff0c;其分母使得类内距离趋向小&＃xff0c;分子使得类间距离趋向大。在原数据维度定义类内散度矩阵以及类间散度矩阵&＃xff0c;改写优化目标如下&＃xff1a;
  
• 上式中 SwS_wSw​ 代表类内散度矩阵with-in&＃xff0c;SbS_bSb​ 代表类间散度矩阵between&＃xff0c;$J$ 称为二者的广义瑞利商Generalized Rayleigh Quotient。目标是最大化 $J$&＃xff0c;无论是直观理解还是从广义瑞利商的定义式出发&＃xff0c;参数 $w$ 的大小并不影响LDA的结果&＃xff0c;只有其方向会影响投影点的分布。
• 不失其一般性&＃xff0c;我们在最大化 $J$ 的过程中加入限制&＃xff0c;令 wTSww&＃61;1w^TS_ww&＃61;1wTSw​w&＃61;1&＃xff0c;那么该优化问题改写为&＃xff1a;
  
• 使用等式约束下的拉格朗日乘数法&＃xff0c;得到如下的推导过程&＃xff1a;
  
• 最终得到了参数 $w$ 的解析解&＃xff0c;实际应用中考虑到数值解的稳定性&＃xff0c;在计算类间散度矩阵 SwS_wSw​ 的逆时会先进行奇异值分解&＃xff0c;而后在求出其逆矩阵。

• 从二分类问题推广到多分类问题&＃xff0c;假定共有 $K$ 个类&＃xff0c;$m$ 个样本&＃xff0c;定义全局散度矩阵total如下所示 &＃xff1a;
  

• 导入数据&＃xff0c;计算均值。

LdaData &＃61; np.loadtxt("data2.txt")
Pos0 &＃61; np.where(LdaData[:,2] &＃61;&＃61; 0)
Pos1 &＃61; np.where(LdaData[:,2] &＃61;&＃61; 1)
X1 &＃61; LdaData[Pos0,0:2]
X1 &＃61; X1[0,:,:]
X2 &＃61; LdaData[Pos1,0:2]
X2 &＃61; X2[0,:,:]MeanX1 &＃61; np.mean(X1,axis&＃61;0)
MeanX2 &＃61; np.mean(X2,axis&＃61;0)
print(X1.shape)
print(X2.shape)
print(MeanX1)


• 计算类内、类间散度矩阵&＃xff0c;本例中基于式 Sbw&＃61;λSwwS_bw&＃61;\lambda S_wwSb​w&＃61;λSw​w&＃xff0c;从而有&＃xff1a;Sw−1Sbw&＃61;λwS_w^{-1}S_bw&＃61;\lambda wSw−1​Sb​w&＃61;λw 对矩阵 Sw−1SbS_w^{-1}S_bSw−1​Sb​ 进行特征分解后选择垂直于 μ0−μ1\mu_0-\mu_1μ0​−μ1​ 的特征向量即可。(下面代码中直接选择主特征值对应的特征向量&＃xff0c;现在看应该是诈胡了&＃xff0c;我也记不清自己当时怎么想的了&＃xff0c;修改后的部分在最后).

def ComputeS_w(X,MeanX):### START THE CODE ###S &＃61; np.dot((X-MeanX).T,X-MeanX) # np.dot 计算类内散度### END THE CODE ###return STestS_i &＃61; ComputeS_w(X1,MeanX1)
print("TestS_i &＃61; ",TestS_i)def ComputeS_b(X1Mean,X2Mean):### START THE CODE ###S &＃61; np.dot(X1Mean-X2Mean,(X1Mean-X2Mean).T) # np.dot 计算类间散度### END THE CODE ###return SSb &＃61; ComputeS_b(MeanX1,MeanX2)
print("Sb &＃61; ", Sb)Sw_1 &＃61; ComputeS_w(X1,MeanX1)
Sw_2 &＃61; ComputeS_w(X2,MeanX2)
Sw &＃61; Sw_1 &＃43; Sw_2
[V,L] &＃61; np.linalg.eig(np.dot(np.linalg.inv(Sw),Sb))
print(V)
print(L)
index &＃61; np.argsort(-V)
W &＃61; L[:,index[0]]
print("W &＃61; ", W)def TranLine(X,K):X_i &＃61; []Y_i &＃61; []for i in range(np.size(X,0)):y_0 &＃61; X[i,1]x_0 &＃61; X[i,0]x1 &＃61; (k*(y_0-b)&＃43;x_0)/(k**2&＃43;1)y1 &＃61; K*x1&＃43;bX_i.append(x1)Y_i.append(y1)return X_i,Y_ib &＃61; 0
x &＃61; np.arange(2,10)
k &＃61; W[1]/W[0]
y&＃61;k*x&＃43;b
plt.plot(x,y)
plt.scatter(X1[:,0],X1[:,1],marker&＃61;&＃39;&＃43;&＃39;,color&＃61;&＃39;r&＃39;,s&＃61;40)
plt.scatter(X2[:,0],X2[:,1],marker&＃61;&＃39;*&＃39;,color&＃61;&＃39;b&＃39;,s&＃61;40)
plt.grid()
plt.show()
X_1,Y_1 &＃61; TranLine(X1,k)
X_2,Y_2 &＃61; TranLine(X2,k)
plt.plot(x,y)
plt.scatter(X1[:,0],X1[:,1],marker&＃61;&＃39;&＃43;&＃39;,color&＃61;&＃39;r&＃39;,s&＃61;60)
plt.scatter(X2[:,0],X2[:,1],marker&＃61;&＃39;*&＃39;,color&＃61;&＃39;b&＃39;,s&＃61;60)
plt.grid()
plt.plot(X_1,Y_1,&＃39;r&＃43;&＃39;)
plt.plot(X_2,Y_2,&＃39;b>&＃39;)


• 【勘误】

def ComputeS_w(X,MeanX):### START THE CODE ###S &＃61; np.dot((X-MeanX).T,X-MeanX) # np.dot 计算类内散度### END THE CODE ###return STestS_i &＃61; ComputeS_w(X1,MeanX1)
print("TestS_i &＃61; ",TestS_i)def ComputeS_b(X1Mean,X2Mean):### START THE CODE ###S &＃61; np.dot(X1Mean-X2Mean,(X1Mean-X2Mean).T) # np.dot 计算类间散度### END THE CODE ###return SSb &＃61; ComputeS_b(MeanX1,MeanX2)
print("Sb &＃61; ", Sb)Sw_1 &＃61; ComputeS_w(X1,MeanX1)
Sw_2 &＃61; ComputeS_w(X2,MeanX2)
Sw &＃61; Sw_1 &＃43; Sw_2
# [V,L] &＃61; np.linalg.eig(np.dot(np.linalg.inv(Sw),Sb))
# print(V)
# print(L)
# index &＃61; np.argsort(-V)
# W &＃61; L[:,index[0]]
W &＃61; np.dot(np.linalg.inv(Sw),MeanX1-MeanX2)
print("W &＃61; ", W)