机器学习(三)：一文读懂线性判别分析（LDA）

一.什么是线性判别分析？

线性判别分析（Linear Discriminant Analysis，LDA）的一种经典的线性学习方法（属于监督学习），这里先借用周志华教授的《机器学习》中的图片来做一个直观的展示：

正如该图中展示的那样，LDA需要寻找一条合适的直线$y=w^{T}x$，使得数据集中的样例投影到该直线时同类样例的投影点尽可能接近，不同类样例的投影点尽可能远离，这样一来就可以对新样例进行分类了，具体做法同样是先将新样例投影到该直线上，然后根据投影点的位置来判断样本的类别。但是实际上LDA用于分类的情况并不是特别常见，其更常见的用法是用来进行高维数据的降维，但是本文还是对LDA分类进行了尝试。

二.推导过程(以二分类为例)

假设数据集$D$中包含$m$个样本，每个样本中$x$代表样例，$y$代表对应样例的标签（target），即 D = { ( x i , y i ) } i = 1 m ,    y i ∈ { 0 , 1 } D=\{(x_i,y_i)\}_{i=1}^m,\;yi \in \{0,1\} D={(xi​,yi​)}i=1m​,yi∈{0,1}。后面的讨论都是基于该假设，为了方便理解，这里先给出一些数学符号介绍：

数学符号	含义
$X_i$	表示第$i$类样例的集合，其中 i ∈ { 0 , 1 } i \in \{0,1\} i∈{0,1}
${\mu }_i$	表示第$i$类样例的均值向量，其中 i ∈ { 0 , 1 } i \in \{0,1\} i∈{0,1}
$\vec{w}$	投影直线的方向向量

2.1 投影

若已知向量$\vec{x}$和向量$\vec{w}$，则向量$\vec{x}$与向量$\vec{w}$上内积可表示为$\vec{x}\cdot \vec{w}=|\vec{x}||\vec{w}|cos\theta$，当$\vec{w}$的模长为1时，两个向量的内积值就变成了$\vec{x}$在$\vec{w}$上的投影，即
$$\begin{array}{rl}\vec{x}\cdot \vec{w}& =|\vec{x}||\vec{w}|cos\theta \\ & =|\vec{x}|cos\theta \end{array}$$

注意：后面内容中的向量都没有采用数学中向量的书写方式，而是表示为矩阵的形式！

2.2 类间散度矩阵

由于需要考虑如何使得不同类样例的投影点尽可能远离，因此需要采用一个标准来衡量类间距离，LDA中采用的是两个类别的均值向量 μ i ,    i ∈ { 0 , 1 } \mu_i,\;i\in\{0,1\} μi​,i∈{0,1}在直线上的投影间的距离，即
$$\begin{array}{rl}(w^T{\mu }_0-w^T{\mu }_1{)}^2& =(w^T({\mu }_0-{\mu }_1){)}^2\\ & =(w^T({\mu }_0-{\mu }_1))(({\mu }_0-{\mu }_1{)}^{T}w)\\ & =w^T(\mu 0-\mu 1)(\mu 0-\mu 1{)}^{T}w\end{array}$$
在《机器学习》一书中用$S_b=(\mu 0-\mu 1)(\mu 0-\mu 1{)}^T$表示类间散度矩阵，则上述公式中的类间距离被改写为了$w^T{S}_{b}w$。

2.3 类内散度矩阵

同样的，也需要提出一个指标来衡量同类样例投影点的接近程度，很容易的我们可以想到方差，方差可以用来度量一组数据间的离散程度，方差越大则数据就越分散。基于此我们可以分别计算出$X_0$和$X_1$的方差（不是真正标准的方差，没有除以样本量）
$$\begin{gathered} var(X_0)=\sum _{x\in X_0}(w^{T}x-w^T{\mu }_0{)}^2 \\ var(X_1)=\sum _{x\in X_1}(w^{T}x-w^T{\mu }_1{)}^2 \end{gathered}$$
两个类别的方差和为
$$\begin{array}{rl}var(X_0)+var(X_1)& =\sum _{x\in X_0}(w^{T}x-w^T{\mu }_0{)}^2+\sum _{x\in X_1}(w^{T}x-w^T{\mu }_1{)}^2\\ & =\sum _{x\in X_0}{w}^T(x-{\mu }_0)(x-{\mu }_0{)}^{T}w+\sum _{x\in X_1}{w}^T(x-{\mu }_1)(x-{\mu }_1{)}^{T}w\\ & =w^T(\sum _{x\in X_0}(x-{\mu }_0)(x-{\mu }_0{)}^T+\sum _{x\in X_1}(x-{\mu }_1)(x-{\mu }_1{)}^T)w\end{array}$$
在《机器学习》一书中，用$\sum _0$和$\sum _1$分别表示$\sum _{x\in X_0}(x-{\mu }_0)(x-{\mu }_0{)}^T$和$\sum _{x\in X_1}(x-{\mu }_1)(x-{\mu }_1{)}^T$，并令$S_w=\sum _0+\sum _1$，则两个类别间的方差可以改写为
$$var(X_0)+var(X_1)=w^T{S}_{w}w$$

2.4 目标函数

根据LDA的目标，我们需要使得同类样例间的投影点尽可能接近，即使得$w^T{S}_{w}w$尽可能小，异类样例间的投影点尽可能远离，即使得$w^T{S}_{b}w$尽可能大，综合可以得到目标函数$J$
$$J=\frac{w^T{S}_{b}w}{w^T{S}_{w}w}$$
可知$J$的值越大，就越符合LDA的目标。观察目标函数$J$可得，若$w^{\prime }$是其解时，对于任意常数$\alpha$来说，$\alpha w$也是目标函数的解，也就是说目标函数的目标解与$w$的长度无关，而与其方向有关。鉴于此可令$w^T{S}_{w}w=1$，则目标函数等价于在约束条件为$w^T{S}_{w}w=1$的情况下，求解$w^T{S}_{b}w$的最大值，也等价于在同样约束条件下，求解$-w^T{S}_{b}w$的最小值，利用拉格朗日乘子法可得
$$L(w,\lambda )=-w^T{S}_{b}w+\lambda (w^T{S}_{w}w-1)$$
其中$\lambda$为拉格朗日乘子，对上式求导并令导数为零可得
$$\frac{\partial L(w,\lambda )}{\partial w}=-2S_{b}w+2\lambda S_{w}w=0$$
即$S_{b}w=\lambda S_{w}w$，又由于$S_{b}w=({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}w$中$({\mu }_0-{\mu }_1{)}^{T}w$为标量，设该标量的值为$k$，因此$S_{b}w=\lambda S_{w}w$可以转换为：
$$({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}w=({\mu }_0-{\mu }_1)k=\lambda S_{w}w$$
化简可得$w$的解为
$$w=S_w^{-1}({\mu }_0-{\mu }_1)\frac{k}{\lambda }$$
而$\frac{k}{\lambda }$省略也对最后的求解不影响，因为只关心其方向，因此可得
$$w=S_w^{-1}({\mu }_0-{\mu }_1)$$

三.多分类LDA

对于多分类任务，假设存在$N$个类，样本总数为$m$，样本中第$i$类示例集合为$X_i$，其所对应的类别样本数为$m_i$，定义全局散度矩阵$S_t$如下：
$$\begin{array}{rl}{S}_t& =S_b+S_w\\ & =\sum _{i=1}^m(x_i-\mu )(x_i-\mu {)}^T\end{array}$$
上式中$\mu$表示所以样本的均值向量，全局散度矩阵实际上就是类内散度矩阵和类间散度矩阵之和。对于类内散度矩阵，重定义为了每个类别的散度矩阵之和，即
$$\begin{array}{rl}{S}_w& =\sum _{i=1}^N{S}_{w_i}\\ & =\sum _{i=1}^N\sum _{x\in X_i}(x-{\mu }_i)(x-{\mu }_i{)}^T\end{array}$$
式子中${\mu }_i$表示类别$i$的均值向量，由上可得$S_b$的表达式，即
$$\begin{array}{rl}{S}_b& =S_t-S_w\\ & =\sum _{i=1}^N\sum _{j=1}^{m_i}(x_j-\mu )(x_j-\mu {)}^T-\sum _{i=1}^N\sum _{x\in X_i}(x-{\mu }_i)(x-{\mu }_i{)}^T\\ & =\sum _{i=1}^N\sum _{j=1}^{m_i}(x_j-\mu )(x_j-\mu {)}^T-\sum _{i=1}^N\sum _{j=1}^{m_i}(x_j-{\mu }_i)(x_j-{\mu }_i{)}^T\\ & =\sum _{i=1}^N\sum _{j=1}^{m_i}[(x_j-\mu )(x_j-\mu {)}^T-(x_j-{\mu }_i)(x_j-{\mu }_i{)}^T]\\ & =\sum _{i=1}^N\sum _{j=1}^{m_i}[(x_j{x}_j^T-x_j{\mu }^T-\mu x_j^T+\mu {\mu }^T)-(x_j{x}_j^T-x_j{\mu }_i^T-{\mu }_i{x}_j^T+{\mu }_i{\mu }_i^T)]\\ & =\sum _{i=1}^N[\sum _{j=1}^{m_i}{x}_j(-{\mu }^T+{\mu }_i^T)+(-\mu +{\mu }_i)\sum _{j=1}^{m_i}{x}_j^T+\sum _{j=1}^{m_i}(\mu {\mu }^T-{\mu }_i{\mu }_i^T)]\\ & =\sum _{i=1}^N[m_i{\mu }_i(-{\mu }^T+{\mu }_i^T)+(-\mu +{\mu }_i)m_i{\mu }_i^T+m_i(\mu {\mu }^T-{\mu }_i{\mu }_i^T)]\\ & =\sum _{i=1}^N{m}_i[-{\mu }_i{\mu }^T+{\mu }_i{\mu }_i^T-\mu {\mu }_i^T+\mu {\mu }^T]\\ & =\sum _{i=1}^N{m}_i({\mu }_i-\mu )({\mu }_i-\mu {)}^T\end{array}$$
多分类LDA可以有多种实现方法，使用$S_b$，$S_w$，$S_t$三者中的任何两个即可，常见的一种时采用优化目标
$$\underset{W}{\max }\frac{tr(W^T{S}_{b}W)}{tr(W^T{S}_{w}W)}$$
其中$w\in R^{d\times (N-1)}$，$tr(.)$表示矩阵的迹，则上式求解同样和二分类一样可以令$tr(W^T{S}_{w}W)=1$，然后使用拉格朗日乘子法（具体过程就不展开了），最后通过求导可以得到下式
$$S_{b}W=\lambda S_{w}W$$
而$W$的闭式解就是**$S_w^{-1}{S}_b$的$d^{\prime }$个最大非零广义特征值所对应的特征向量组成的矩阵**，$d^{\prime }\le N-1$。

四.降维和分类

前面两节中介绍了如何求取二分类和多分类任务中的$w$，那么在获取到$w$后就可以进行降维及建立在降维之上的分类操作了。假设数据集 D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . ( x m , y m ) } D=\{(x_1,y_1),(x_2,y_2),...(x_m,y_m)\} D={(x1​,y1​),(x2​,y2​),...(xm​,ym​)}，其中$x_k=(x_{k1},x_{k2},...,x_{kd}),k\in [1,m]$，基于该数据集求取的投影矩阵为$W$，那么降维操作

就是对数据集$D$中的每个样本特征$x_i$进行投影操作，即
$$x_i^{\prime }=w^T{x}_i$$
至于分类，在线性判别分析LDA原理总结一文中我发现作者提到了一种常用的思想：假设各个类别的样本数据符合高斯分布，在LDA进行降维操作后，先计算各个类别投影数据的均值和方差，进而得到该类别高斯分布的概率密度函数。然后对新样本先进行同样的降维操作，然后分别代入各个类别的高斯分布概率密度函数中去计算概率，概率最高的类别即为新样本的分类类别。

五.实现

5.1 数据集说明

本次实验采用的数据集是UCI上的鸢尾花卉数据集一共包含150个样本，分为三类(Setosa，Versicolour，Virginica)，每类50个数据，该数据集的部分数据展示如下：

其中各个属性的描述具体见下表：

属性	属性描述
sepal.length	花萼长度
sepal.width	花萼宽度
petal.length	花瓣长度
petal.width	花瓣宽度
variety	花类型(0表示Setosa, 1表示Versicolour 2表示Virginica)

5.2 示例代码

首先，由于iris数据集较小，因此训练集和测试集的划分比例为$7:3$。

from os import startfile
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split

def load_iris(file_path='./data/iris.csv',test_size=0.3):
    """
    功能：加载数据集并划分训练集和测试集
    """
    df = pd.read_csv(file_path)
    x,y = df.iloc[:,0:-1].values,df.iloc[:,-1].values
    #打乱并划分测试集和训练集
    train_x,test_x,train_y,test_y = train_test_split(x,y,test_size=test_size,random_state=5)
    return train_x,test_x,train_y,test_y

然后，我将LDA模型其封装为一个类以供调用，该类需要传入三个参数：数据特征，数据标签以及降维后的维度。

import numpy as np

class LDA(object):
    """
    线性判别分析
    """
    def __init__(self,data,target,d) -> None:
        self.data = data
        self.target = target
        self.d = d
        self.labels = set(target)
        self.mu = self.data.mean(axis=0)
    
    def divide(self):
        """
        功能：将传入的数据集按target分成不同的类别集合并求出对应集合的均值向量
        """
        self.classify,self.classmu = {},{}
        for label in self.labels:
            self.classify[label] = self.data[self.target==label]
            self.classmu[label] = self.classify[label].mean(axis=0)
    
    def getSt(self):
        """
        功能：计算全局散度矩阵
        """
        self.St = np.dot((self.data-self.mu).T,(self.data-self.mu))

    def getSb(self):
        """
        功能：计算类内散度矩阵
        """
        self.Sb = np.zeros((self.data.shape[1],self.data.shape[1]))
        for i in self.labels:
            #获取类别i样例的集合
            classi = self.classify[i]
            #获取类别i的均值向量
            mui = self.classmu[i]
            self.Sb += len(classi) * np.dot((mui - self.mu).reshape(-1,1),(mui - self.mu).reshape(1,-1))

    def getW(self):
        """
        功能：计算w
        """
        self.divide()
        self.getSt()
        self.getSb()
        #St = Sw + Sb
        self.Sw = self.St - self.Sb 
        #计算Sw-1*Sb的特征值和特征向量
        #eig_vectors[:i]与 eig_values相对应
        eig_values, eig_vectors = np.linalg.eig(np.linalg.inv(self.Sw).dot(self.Sb))
        #寻找d个最大非零广义特征值
        topd = (np.argsort(eig_values)[::-1])[:self.d]
        #用d个最大非零广义特征值组成的向量组成w
        self.w = eig_vectors[:,topd]
        



if __name__ == "__main__":
    x = np.array([[1,2,3],[2,1,3],[2,4,1],[1,3,2],[3,6,4],[3,1,1]])
    y = np.array([0,1,2,0,1,2])
    lda = LDA(x,y,2)
    lda.getW()

最后，对于模型的训练与验证，我是先利用训练集获取投影矩阵$W$，然后对训练集数据进行降维操作，按照第四节中的分类思想需要将对降维后各个类别的数据类别求各自的均值和标准差（确定高斯分布的概率密度函数），然后再对测试集利用投影矩阵进行降维，将样本依次代入各个类别的高斯分布函数，选取概率最大的类别作为样本的预测类别。

from LDAModel import LDA
from load_data import load_iris
import numpy as np
from scipy.stats import norm
from sklearn.metrics import accuracy_score

def judge(gauss_dist,x):
    """
    功能：判断样本x属于哪个类别
    """
    #将样本带入各个类别的高斯分布概率密度函数进行计算
    outcome = [[k,norm.pdf(x,loc=v['loc'],scale=v['scale'])] for k,v in gauss_dist.items()]
    #寻找计算结果最大的类别
    outcome.sort(key=lambda s:s[1],reverse=True)
    return outcome[0][0]

def Test():
    """
    功能：对测试集进行分类并返回准确率
    """
    #加载数据集
    train_x,test_x,train_y,test_y = load_iris()
    #创建模型
    lda = LDA(train_x,train_y,1)
    #获取投影矩阵w
    lda.getW()
    #对训练集进行降维
    train_x_new = np.dot((train_x), lda.w)
    #获取训练集各个类别对应的高斯分布的均值和方差
    gauss_dist = {}
    for i in lda.labels:
        classi = train_x_new[train_y==i]
        loc = classi.mean()
        scale = classi.std()
        gauss_dist[i] = {'loc':loc,'scale':scale}
    test_x_new = np.dot(test_x,lda.w)
    pred_y = np.array([judge(gauss_dist,x) for x in test_x_new])
    
    return accuracy_score(test_y,pred_y)
    

if __name__ == "__main__":
    acc = Test()
    print(acc)

运行结果表明LDA模型在测试集上的预测准确率约为$0.93$，对于训练集的降维操作最后也作一个可视化展示，由于最后降维至一维（直线上）为了方便绘制散点图我给这些点的纵坐标统一置为1，只需看横轴即可，从结果可以看出鸢尾花投影后的确将几个类别区分了出来。

参考资料

《机器学习》周志华

LDA深入补充推导和python实现（二分类问题）

Dimensionality Reduction——LDA线性判别分析原理篇

线性判别分析（LDA）实现二分类的思路

以上便是本文的全部内容，要是觉得不错的话，就点个赞或关注一下吧，你们的支持是博主创作的不竭动力！