线性判别分析(LDA)二分类的原理及python实现

问题引入

已知一个样本存在两个属性，根据属性的不同可以分为三个类别，现需要依据下面的训练集，给出某样本的属性值判断该样本属于下面的哪一个类别。

属性1	属性2	类别
4	2	1
2	4	1
2	3	1
3	6	1
4	4	1
9	10	0
6	8	0
9	5	0
8	7	0
10	8	0

---

符号说明

$n$表示列属性的个数。
${\text{x}}_i$为nx1的列向量，在表格中代表一个不包括类别属性的行。
$c=0,1,\cdots ,k,\cdots$表示类别的集合，对应上述表格中的0、1。
$m^k$表示类别为c的样本总个数。
${\text{X}}^k$表示类别k组成的矩阵，大小为$m^k\times n$。
${\text{x}}_i^k,k\in c$为nx1列向量，在表格中表示类别为k的属性行。
$m={\sum }_{k\in c}{m}^k$表示所有样本数。
$\theta$表示位于直线$l$上的nx1的单位列向量。

---

基本思想

之前我们讲过logistic回归进行分类，现在我们要以一种新的方式LDA进行分类，LDA常常用于二分类问题以及多分类的降维处理，本篇文章先讲解关于二分类的问题。
由于两类数据过于繁多，导致无法直接对数据属于哪一类进行判断，LDA的基本思想是对于两类数据，想办法将两类数据映射到一条合适的直线上，此时只需要找到直线上一个可以将两类数据分开到两边的点，就可以直观且简便地判断某个样本数据属于哪一类了。
但并不是所有直线都可以找到这个分界点，如下图所示，紫色与红色分别代表两类数据，左侧将两类数据的一部分映射到一起了显然增加了判断难度，右侧可以找到一个蓝色的点将两类数据区分开来,相对于左侧图像而言，右侧映射到直线的数据，具有类内尽可能聚集，不同类数据尽可能分散的特点，如何利用数学公式去描述这种特点，就是Linear Discriminant Analysis(LDA)的任务了。

---

推导过程

首先需要理解映射的概念，具体内容在PCA降维的投影部分中讲解过，简单来说，若单位向量$\theta$是直线$l$上的单位向量，那么样本数据${\text{x}}_i$在直线$l$上的投影距离为${\theta }^T{\text{x}}_i$。
类别k的样本平均向量为
$${\text{u}}^k=\frac{1}{m^k}\sum _{i=1}^{m^k}{\text{x}}_i,k\in c\text{\hspace{1em}(1)}$$
基本思想中说类内距离尽可能小，描述类内距离的方式千千万，让所有样本向量向平均向量聚集却不失为一种好的策略，那么类别k中的样本向量${\text{x}}_i^k$在映射到直线$l$之后，与平均向量之间的距离(注意该距离不是欧式距离)可以表示为
$$\sum _{i=1}^{m^k}({\theta }^T{\text{x}}_i^k-\theta {}^T{\text{u}}^k{)}^2$$
所有类别的距离和为
$$\sum _{k\in c}\sum _{i=1}^{m^k}({\theta }^T{\text{x}}_i^k-\theta {}^T{\text{u}}^k{)}^2$$
这就表示我们在基本思想中提到过的类内距离和。
全局最小并不代表每一类中都是最小的，即便存在某条直线$l1$，使其中某一类的类内距离比在$l$中还要小，但是求全局最小也不失为一种比较好的策略。
将该公式改写成下面的形式
$$\begin{array}{c}\begin{array}{rl}& \sum _{k\in c}\sum _{i=1}^{m^k}({\theta }^T{\text{x}}_i^k-\theta {}^T{\text{u}}^k{)}^2\\ =& \sum _{k\in c}\sum _{i=1}^{m^k}[{\theta }^T({\text{x}}_i^k-{\text{u}}^k)][{\theta }^T({\text{x}}_i^k-{\text{u}}^k)]\\ =& \sum _{k\in c}\sum _{i=1}^{m^k}[{\theta }^T({\text{x}}_i^k-{\text{u}}^k)][({\text{x}}_i^k-{\text{u}}^k{)}^T\theta ]\\ =& \sum _{k\in c}\sum _{i=1}^{m^k}{\theta }^T({\text{x}}_i^k-{\text{u}}^k)({\text{x}}_i^k-{\text{u}}^k{)}^T\theta \\ =& {\theta }^T[\sum _{k\in c}\sum _{i=1}^{m^k}({\text{x}}_i^k-{\text{u}}^k)({\text{x}}_i^k-{\text{u}}^k{)}^T]\theta \\ =& {\theta }^T{\text{S}}_w\theta \end{array}\end{array}$$
数学感觉较好的网友肯定可以意识到,${\text{S}}_w$(w表示within class)就是所有类别的散度矩阵和了，本人虽然没有拜读过先人们的论文，但能隐约地感觉到散度矩阵就是像这样慢慢地抽象并产生其数学意义的。
我们已经了解到类内距离了，接下来开始表述类间距离。
对于二分类问题，类间距离可以定义为两类平均向量之间的距离，如下。
$$\begin{array}{c}\begin{array}{rl}& ({\theta }^T{\text{u}}^0-{\theta }^T{\text{u}}^1{)}^2\\ =& {\theta }^T({\text{u}}^0-{\text{u}}^1)({\text{u}}^0-{\text{u}}^1{)}^T\theta \\ =& {\theta }^T{\text{S}}_b\theta \end{array}\end{array}\text{\hspace{1em}(2)}$$
与${\text{S}}_w$相似，${\text{S}}_b$(b表示between class)表示类间的散度矩阵。
在基本思想中我们提到，需要使类间距离最大，类内距离最小，我们的求解目标变为
$$\arg \underset{\theta }{\max }\frac{{\theta }^T{\text{S}}_b\theta }{{\theta }^T{\text{S}}_w\theta }\text{\hspace{1em}(3)}$$
即存在分子又存在分母，利用求导的方式求导求极值比较复杂，分母求最小值，总会存在一个最小值t使${\theta }^T{\text{S}}_w\theta =t$,t为一个常数，一般我们默认求极小值，分子部分可以加个负号，变为$-{\theta }^T{\text{S}}_b\theta$，此时求解目标变为
$$\begin{array}{c}\arg {\min }_{\theta }-{\theta }^T{\text{S}}_b\theta \\ {\theta }^T{\text{S}}_w\theta =t\end{array}\text{\hspace{1em}(4)}$$
教材中的求解方式都是以分数形式（3）求解的，个人认为直接列出公式（4）求解也并无不妥，此时就可以利用拉格朗日求极值。
$$J(\theta )=-{\theta }^T{\text{S}}_b\theta +\lambda ({\theta }^T{\text{S}}_w\theta -t)$$
求导数,令其为0。
$$\frac{dJ(\theta )}{d\theta }=-2{\text{S}}_b\theta +\lambda 2{\text{S}}_w\theta =0$$
$⟺$
$${\text{S}}_b\theta =\lambda {\text{S}}_w\theta \text{\hspace{1em}(5)}$$
由（2）知
$${\text{S}}_b\theta =({\text{u}}^0-{\text{u}}^1)({\text{u}}^0-{\text{u}}^1{)}^T\theta$$
由于$({\text{u}}^0-{\text{u}}^1{)}^T\theta$是一个常量，不妨令${\lambda }_{\theta }=({\text{u}}^0-{\text{u}}^1{)}^T\theta$
那么
$${\text{S}}_b\theta =({\text{u}}^0-{\text{u}}^1){\lambda }_{\theta }=\lambda {\text{S}}_w\theta$$
最终可以求得
$$\theta ={\text{S}}_w^{-1}({\text{u}}^0-{\text{u}}^1)\frac{{\lambda }_{\theta }}{\lambda }$$
由于$\theta$是与直线$l$同向的向量，所以常数项可以省略，最终目标解${\theta }^{\ast }$为
$${\theta }^{\ast }={\text{S}}_w^{-1}({\text{u}}^0-{\text{u}}^1)\text{\hspace{1em}(6)}$$
同时，这一结论也表明LDA在二分类中，只能求得一类平行的向量。

---

计算步骤以及python代码

1.计算平均值$u^k$
2.计算散度矩阵${\text{S}}_w$
3.计算逆${\text{S}}_w^{-1}$
4.计算${\theta }^{\ast }=textbf{S}_w^{-1}({\text{u}}^0-{\text{u}}^1)$
5.计算样本在直线上的位置${\theta }^T{\text{x}}_i$
将上述问题计算之后，并绘图

计算结果

python代码

import matplotlib.pyplot as plt
import numpy as np


def draw1(data,theta,num,colors):
	x=[[] for i in range(2)]
	for ele in data:
		x[ele[-1]].append(ele[:-1])
	#绘制原始数据
	for i in range(num):
		x[i]=np.array(x[i])
		plt.scatter(x[i][:,0],x[i][:,1],color=colors[i])
	#绘制被的映射直线
	plt.plot([0,theta[0]*15],[0,theta[1]*15])
	#绘制映射到直线上的点
	for i in range(num):
		for ele in x[i]:
			ta=theta*np.dot(ele,theta)
			plt.plot([ele[0],ta[0]],[ele[1],ta[1]],color=colors[i],linestyle="--")
			plt.scatter(ta[0],ta[1],color=colors[i])
	plt.show()

def c2(data,num):
	n=data.shape[1]-1
	#存储x值
	x=[[] for i in range(num)]
	#存储u值
	u=[[]for i in range(num)]
	#存储Sw
	sw=np.zeros([n,n])
	for ele in data:
		x[ele[-1]].append(ele[:-1])
	for i in range(num):
		x[i]=np.array(x[i])
		u[i]=np.mean(x[i],axis=0)
	print("1计算平均值：\n",u)
	for i in range(num):
		x[i]=x[i]-u[i]
		sw=sw+np.dot(x[i].T,x[i])
	print("2x_i去中心化:\n",x)
	print("3计算散度矩阵S_w:\n",sw)
	#计算theta
	theta=np.dot(np.linalg.inv(sw),(u[0]-u[1]).T)
	#单位化
	fm=0
	for i in range(n):
		fm=fm+theta[i]**2
	return theta/np.sqrt(fm)


data=np.array([[4,2,1],[2,4,1],[2,3,1],[3,6,1],[4,4,1],[9,10,0],[6,8,0],[9,5,0],[8,7,0],[10,8,0]])
colors=['red','green']
theta=c2(data,2)
print("4计算直线向量theta:\n",theta)
draw1(data,theta,2,colors)