LDA在多分类、降维中的应用以及python实现

问题引入

现有三种鸢尾花分别标为0,1,2类，存在4种不同的属性，即萼片长(setosa)，萼片宽(versicolor)，花瓣长与花瓣的宽度，现在需要对这三种鸢尾花的分布在二维平面中进行可视化分析。数据来自于sklearn中的dataset部分数据。

高维数据进行可视化分析，常用思路是对数据进行降维。

符号说明

$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的单位列向量。

前情提要

线性判别分析在二分类的应用中，我们有类内距离为
$${\theta }^T{\text{S}}_w\theta ={\theta }^T\sum _{k\in c}\sum _{i=1}^{m^k}({\text{x}}_i^k-u^k)({\text{x}}_i^k-u^k{)}^T\theta$$
在看很多资料的时候，会发现人们习惯于直接定义类内散度为
$${\text{S}}_w=\sum _{k\in c}\sum _{i=1}^{m^k}({\text{x}}_i^k-u^k)({\text{x}}_i^k-u^k{)}^T$$
然后再定义类内散度
$${\text{S}}_b=({\text{u}}^0-{\text{u}}^1)({\text{u}}^0-{\text{u}}^1{)}^T$$
最后再在散度两边乘以直线向量$\theta$得到映射后的类内距离与类间距离，这种定义初次看来毫无实际意义，但是熟悉PCA以及LDA的二分类之后，就发现这种定义方法所代表数学意义很有用，具体描述如下。
如果存在$\Theta =({\theta }_1,{\theta }_2,...,{\theta }_i,...{\theta }_n)$组成一个新的标准正交基，那么${\Theta }^T{\text{S}}_w\Theta$就是新的标准正交基下的类内散度矩阵。对于对角线上的元素代表位于直线向量${\theta }_i$上的类内距离即${\theta }_i{\text{S}}_w{\theta }_i$。
既然这种定义方式存在其意义，那么我们在本文中就采用这种直接定义散度矩阵的方式进行介绍。

推导过程

首先定义类内散度矩阵，与LDA的二分类一样
$${\text{S}}_w=\sum _{k\in c}\sum _{i=1}^{m^k}({\text{x}}_i^k-u^k)({\text{x}}_i^k-u^k{)}^T\text{\hspace{1em}(1)}$$
再定义类间散度矩阵，与LDA稍微有些不同
$$\begin{array}{c}\begin{array}{rl}{\text{S}}_b& =\sum _{k_1\in c}\sum _{k_2\in c}{m}^{k1}{m}^{k2}({\text{u}}^{k_1}-{\text{u}}^{k_2})({\text{u}}^{k_1}-{\text{u}}^{k_2}{)}^T\\ & =\sum _{k\in c}{m}^k({\text{u}}^k-\text{u})({\text{u}}^k-\text{u}{)}^T\end{array}\end{array}\text{\hspace{1em}(2)}$$
其中
$$\text{u}=\frac{1}{m}\sum _{i=1}^m{\text{x}}_i$$
每个平均向量可以代表该类，将每个类之间的散度矩阵相加，可以表示总体的散度矩阵，这些倒还可以理解，肯定会有人质疑前面为什么加上系数，为什么不定义成${\text{S}}_b={\sum }_{k_1\in c}{\sum }_{k_2\in c}({\text{u}}^{k_1}-{\text{u}}^{k_2})({\text{u}}^{k_1}-{\text{u}}^{k_2}{)}^T$形式，简单地阐述一下这种定义方式的合理性，我们知道二分类中定义成${\text{S}}_b=({\text{u}}^0-{\text{u}}^1)({\text{u}}^0-{\text{u}}^1{)}^T$，前面没有系数，${\text{S}}_b$就算前面加系数也只不过是对矩阵等比例放大而已，对结果并不会产生任何影响，但是对于类别较多的数据就不一样了，比如说0类与1类的训练数据有100个，而2类的训练数据只有1个，那么再求解的过程中将权重更加偏向于那些数据集本身较多的类别会变得更加合理，所以${\text{S}}_b$如此定义。
那么对于某个直线向量$\theta$而言，类内距离为
$${\theta }^T{\text{S}}_w\theta \text{\hspace{1em}(3)}$$
类间距离为
$${\theta }^T{\text{S}}_b\theta \text{\hspace{1em}(4)}$$
求解目标仍然可以写为
$$\arg \underset{\theta }{\max }\frac{{\theta }^T{\text{S}}_b\theta }{{\theta }^T{\text{S}}_w\theta }\text{\hspace{1em}(5)}$$
这种形式仍然可以用于分类，只不过是多分类罢了。
至于降维的应用，实际上还是需要数学理论的支持，首先明确一点，当数据从高维降到低维之后，各个维度之间一定是正交的，而且数据的特征会尽可能地被保留，LDA既然可以分类肯定可以保留特征，那么下面的推导就是证明LDA求解出来的多个极值解${\theta }_i$之间是正交的直线罢了。
具体内容见LDA二分类,构造拉格朗日求极值，其中
$$\text{J}(\theta )={\theta }^T{\text{S}}_w\theta +\lambda ({\theta }^T{\text{S}}_b\theta -t)$$
对$\theta$求导=0之后可以得到
$${\text{S}}_b\theta =\lambda {\text{S}}_w\theta$$
$⟺$
$${\text{S}}_w^{-1}{\text{S}}_b\theta =\lambda \theta \text{\hspace{1em}(6)}$$

此时，$\lambda$为${\text{S}}_w^{-1}{\text{S}}_b$的特征值，$\theta$为对应的特征向量，根据特征向量的性质，对于不同的${\lambda }_i$之间${\theta }_i$肯定是正交的，那么对于相同的$\lambda$可以通过斯密特正交化方法正交化单位化处理成正交的形式，现在LDA可以进行降维才得到了数学上的保证。
对于多个$\theta$而言，选取哪一个$\theta$作为降维后的轴与PCA一样需要看$\lambda$取值，证明如下：
${\text{S}}_w$为实对称矩阵，那么存在
$${\text{S}}_w\theta =t_1\theta$$
带入(6)中得到
$${\text{S}}_b\theta =\lambda t_1\theta \text{\hspace{1em}(7)}$$
与PCA选取较大的$\lambda$一样，此处不再累述，在公式(6)中有，$\lambda t_1$越大，${\theta }^T{\text{S}}_b\theta$越大，类间会越分散。
根据拉格朗日的约定条件${\theta }^T{\text{S}}_w\theta =t$可以知道$t_1$为一个定值，且$t_1=t$。
所以，$\lambda$越大，${\theta }^T{\text{S}}_b\theta$越大，类间会越分散。
最终得到结论是选取$\lambda$较大的矩阵。
那么对于降维而言，公式(5)就变成
$$\arg \text{imax}\frac{{\theta }^T{\text{S}}_b\theta }{{\theta }^T{\text{S}}_w\theta }\text{\hspace{1em}(8)}$$
其中imax代表求极值
求解方式与公式(6)一样，如下
$${\text{S}}_w^{-1}{\text{S}}_b\theta =\lambda \theta$$
如果需要降维到d维，只需要求出前d个较大的$\lambda$所对应的$\theta$即可。

补充知识

1. 为什么求解公式定义成$$\arg \underset{\theta }{\max }\frac{{\theta }^T{\text{S}}_b\theta }{{\theta }^T{\text{S}}_w\theta }$$
却不直接定义成$$\arg \underset{\theta }{\min }\frac{{\theta }^T{\text{S}}_w\theta }{{\theta }^T{\text{S}}_b\theta }$$
实际上，求最小值的方式也对，利用拉格朗日求解之后就会发现，$\lambda$越小，${\theta }^T{\text{S}}_w\theta$越小，我们目标是求解类内距离最小，最终降维后选取的是前d个较小的$\lambda$，与PCA结论恰好相反，可能为了方便记忆，故意写成第一个公式进行求解。
2. 公式(8)是本人根据自己理解写出的公式，实际上教材中的公式常常写作
$$\arg \max \frac{tr({\Theta }^T{\text{S}}_b\Theta )}{tr({\Theta }^T{\text{S}}_w\Theta )}$$
$$\Theta =({\theta }_1,{\theta }_2,...,{\theta }_i,...{\theta }_n)$$
这种写法比较简洁，但是对于刚刚入门的人而言可能有点难以理解，求迹就是求对角线和，等价于
$$\arg \max \frac{{\sum }_{i=1}^n({\theta }_i^T{\text{S}}_b{\theta }_i)}{{\sum }_{i=1}^n({\theta }_i^T{\text{S}}_w{\theta }_i)}$$
写成拉格朗日的形式就会发现，对${\theta }_i$求导，那么${\theta }_j,i\ne j$相当于是个常数，求导之后含有${\theta }_j$的项变成0，结果不变。
3. 构造拉格朗日的时候，本人将分母写作常数t，很多教材都写作1，求导结果都一样，等于1的写法，是数学中的一个求解定理，需要拥有相应的数学知识才能立即反应过来，实际上$\theta$是个直线，也就是说，不管是$\theta$还是$\sqrt{t}\theta$都是正确解，所以分母等于常数t还是等于1，都可以约掉，本人希望尽可能将数学知识限定在非数学系大学本科数学教学范围之内，不愿做过多的拓展，故意写作成常数t。

python代码

利用sklearn降维之后的图如下，代码来自于sklearn官方文档

import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

iris = datasets.load_iris()

X = iris.data
y = iris.target
target_names = iris.target_names

lda = LDA(n_components=2)
X_r2 = lda.fit(X, y).transform(X)

plt.figure()
for c, i, target_name in zip("rgb", [0, 1, 2], target_names):
    plt.scatter(X_r2[y == i, 0], X_r2[y == i, 1], c=c, label=target_name)
plt.legend()
plt.title('LDA of IRIS dataset')

plt.show()

下面是自己的python代码实现结果，虽然形式上有区别，但两者都可以区分三个不同种类

源码如下

from sklearn import datasets
import numpy as np
import matplotlib.pyplot as plt
#有num个种类
#降到k维
def desc(num,k):
	iris = datasets.load_iris()
	xa = iris.data
	y = iris.target
	target_names = iris.target_names
	#所有元素的平均值
	ua=np.array([np.mean(xa,axis=0)])
	n=xa.shape[1]
	#存储x值
	x=[[]for i in range(num)]
	u=[[]for i in range(num)]
	sw=np.zeros([n,n])
	sb=np.zeros([n,n])
	for i in range(y.shape[0]):
		x[y[i]].append(xa[i])
	for i in range(num):
		x[i]=np.array(x[i])
		u[i]=np.array([np.mean(x[i],axis=0)])
	#去中心化,计算Sw
	for i in range(num):
		x[i]=x[i]-u[i]
		sw=sw+np.dot(x[i].T,x[i])
		sb=sb+x[i].shape[0]*np.dot((u[i]-ua).T,(u[i]-ua))
	#求Sw_1Sb的特征值与特征向量
	eig_value,eig_vec=np.linalg.eig(np.dot(np.linalg.inv(sw),sb))
	print(eig_vec)
	#对特征值进行排序，返回标号
	index=np.argsort(eig_value)
	#获取取最后k个特征向量
	eig_vec_k=eig_vec[:,index[:-k-1:-1]]
	colors=["r","g","b"]
	#绘图
	for i in range(xa.shape[0]):
		plt.scatter(np.dot(xa[i],eig_vec_k[:,0]),np.dot(xa[i],eig_vec_k[:,1]),color=colors[y[i]])
	plt.title("my LDA")
	plt.show()

desc(3,2)