线性判别分析LDA算法与python实现

  降维指的是通过某种数学变换将高维原始空间的属性转变为低维子空间,根据变换形式可将该数学变换分为线性变换和非线性变换,对应的降维算法也被称为线性降维和非线性降维.其中,线性降维算法主要有线性判别分析(linear discriminant analysis，LDA)和主成分分析(Principal Component Analysis，PCA),非线性降维主要有核化思想(如Kernelized PCA)和流形学习(Isomap,LLE,LE等)两类.
  假定有原始空间高维数据$X\in R^{n\times m}$,其中$n$为样本数,$m$为样本长度,现我们要求$X$的低维嵌入$Y\in R^{n\times d}$,其中$d<<m$.线性降维的思想就是求一个权重矩阵$W\in R^{m\times d}$,使用$W$对$X$进行线性变换$Y=XW$,使得变换前后的数据分布一致.如上所述,LDA和PCA都是线性降维算法,不同的是LDA是监督学习算法,而PCA是面向无标签的数据样本.本文介绍LDA算法.

  首先盗用西瓜书里的一张图,这张图很清晰地阐释了LDA的核心思想:图中的数据简化为二维降维到一维,降维过程中,LDA算法使得低维空间中,同一类的数据尽可能接近,使得不同类数据尽可能远离.如上文所述,LDA是一种监督学习算法,即数据具有label,这里与西瓜书保持一致,使用了二分类问题的数据,分别记为$X_0\in R^{n_0\times m}$与$X_1\in R^{n_1\times m}$.我们的目标是找到一个变换矩阵$W\in R^{m\times d}$对原始数据$X$进行线性变换$Y=XW$,且变换后的$Y$满足上述性质.记${\mu }_0\in R^{m\times 1},{\mu }_1\in R^{m\times 1},{\Sigma }_0\in R^{m\times m},{\Sigma }_1\in R^{m\times m}$分别为$X_0$的均值,$X_1$的均值,$X_0$的协方差,$X_1$的协方差,则: $${\mu }_i=\frac{1}{n_i}\sum _{x\in X_i}x$$ $${\Sigma }_i=\sum _{x\in X_i}(x-{\mu }_i)(x-{\mu }_i{)}^T$$ 首先我们希望变换后的类间距离越大越好,我们定义类间距离为类中心的$l_2$距离,所以该步骤我们的目标是: $${\mathrm{m}\mathrm{a}\mathrm{x}}_W||W^T{\mu }_0-W^T{\mu }_1|{|}_2^2$$ 即:$${\mathrm{m}\mathrm{a}\mathrm{x}}_W{W}^T({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}W$$ 其次我们希望变换后的类内协方差越小越好,即: $${\mathrm{m}\mathrm{i}\mathrm{n}}_W{W}_T({\Sigma }_0+{\Sigma }_1)W$$   现定义两个矩阵,类内散度矩阵(intra-class scatter matrix) $S_a\in R^{m\times m}$与类间散度矩阵(inter-class scatter matrix) $S_r\in R^{m\times m}$:$$S_a={\Sigma }_0+{\Sigma }_1$$ $$S_r=({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^T$$ 我们约束$W^T{S}_{a}W=1$,所以最后的优化问题可以写成: $${\mathrm{m}\mathrm{i}\mathrm{n}}_W-W^T{S}_{r}W$$ $$st.W^T{S}_{a}W=1$$ 定义拉格朗日函数为: $$L(W)=-W^T{S}_{r}W+\lambda (W^T{S}_{a}W-1)$$ 对上述方程求$W$的偏导,得到: $$S_{r}W=\lambda S_{a}W$$ 由上式可知,$W\in R^{m\times d}$的闭解为矩阵$S_a^{-1}{S}_r$最大的$d$个特征值对应的$m$维特征向量.这里公布一下代码和实验结果,代码略简略,只考虑了三维降到二维.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets.samples_generator import make_classification
from mpl_toolkits.mplot3d import Axes3D
np.random.seed(0)

def createDataSet(n=100):
	X0 = np.array(np.random.random((n, 3)) + 1)
	X1 = np.array(np.random.random((n, 3)) - 1)
	return X0, X1

def compute_W(X0, X1):
	mean0 = np.mean(X0, axis=0)
	mean1 = np.mean(X1, axis=0)
	diff = np.expand_dims(mean0 - mean1, axis = 1)
	Sr = np.matmul(diff, diff.T)
	cov0 = np.matmul((X0 - mean0).T, (X0 - mean0))
	cov1 = np.matmul((X1 - mean1).T, (X1 - mean1))
	Sa = cov0 + cov1
	# print(Sa)
	tmp = np.matmul(np.linalg.pinv(Sa), Sr)
	# tmp = np.matmul(np.linalg.pinv(-Sr), Sa)
	eig_val, eig_vector = np.linalg.eig(tmp)
	return eig_vector[:, [0,1]]

def LDA(X0, X1):
	W = compute_W(X0, X1)
	return np.matmul(X0, W), np.matmul(X1, W)

if __name__ == '__main__':
	X0, X1 = createDataSet()
	Y0, Y1 = LDA(X0, X1)

	fig = plt.figure()
	ax = fig.add_subplot(121, projection='3d')
	ax.scatter(X0[:, 0], X0[:, 1], X0[:, 2], c=[1,0,0], cmap=plt.cm.hot)
	ax.scatter(X1[:, 0], X1[:, 1], X1[:, 2], c=[0,0,1], cmap=plt.cm.hot)
	ax2 = fig.add_subplot(122)
	ax2.scatter(Y0[:, 0], Y0[:, 1], c=[1,0,0], cmap=plt.cm.hot)
	ax2.scatter(Y1[:, 0], Y1[:, 1], c=[0,0,1], cmap=plt.cm.hot)
	plt.show()

  在LDA中,我们约束$W^T{S}_{a}W=1$,可能是提出算法的学者觉得类内相似对比类间差异不那么重要吧,现在我们探索一下另一种情况,我们约束$W^T{S}_{r}W=1$,那么优化问题变成了:$${\mathrm{m}\mathrm{i}\mathrm{n}}_W{W}^T{S}_{a}W$$ $$st.W^T{S}_{r}W=1$$ 定义拉格朗日函数为: $$L(W)=W^T{S}_{a}W+\lambda (W^T{S}_{r}W-1)$$ 对上述方程求$W$的偏导,得到: $$-S_{a}W=\lambda S_{r}W$$ 由上式可知,$W\in R^{m\times d}$的闭解为矩阵$-S_r^{-1}{S}_a$最大的$d$个特征值对应的$m$维特征向量.代码中只需将21行注释,并恢复22行即可,下图展示了用这种约束得到的实验结果,可以看出两种约束并没有什么很大的差异,当然可能在高阶上第一种方法表现更优异,这里就不往下探索了.