线性判别分析 LDA（Linear Discriminant Analysis）

最新推荐文章于 2024-04-20 22:13:34 发布

姜萌芽

最新推荐文章于 2024-04-20 22:13:34 发布

阅读量1.8k

点赞数

分类专栏：机器学习文章标签：机器学习算法人工智能

本文链接：https://blog.csdn.net/jiangmengya1/article/details/127010384

版权

机器学习专栏收录该内容

7 篇文章 0 订阅

订阅专栏

算法目标

线性判别分析，主要用于分类（Classification）、降维（Dimensionality Reduction）领域。

基本思想是，找到一个单位向量 $w$ ，将数据集中的样例映射到 $w$ 上，使得同类别的样例的投影点尽可能接近，不同类别的样例的投影点尽可能远离。

如上图所示，比起左侧的映射，右侧的映射更符合我们对 LDA 的期待。

正交投影

正交即垂直，我们来看如何将向量 $x$ 正交投影到向量 $w$ 上，及其具体公式是什么。

向量 $p$ 和 $d$ 承载了由 $x$ 向 $w$ 的这种正交投影关系，所以我们对这种关系，进行计算具象化：

由于向量平行，则 $p=cw$ ，

由于向量加减，则 $d=x-p =x-cw$

由于向量垂直，则 $0 = p^{T}d = (cw)^{T}(x-cw) = cw^{T}x-c^{2}w^{T}w$ ，即

$cw^{T}x = c^{2}w^{T}w$

$c=\frac{w^{T}x}{w^{T}w}$

那么， $p = cw$ 可以重写为：

$p= (\frac{w^{T}x}{w^{T}w})w$

即为，正交投影公式

投影简化

$w$ 代表的是被投影的方向。为了简化正交投影公式，我们不妨假设 $w$ 为单位向量。

单位向量具有性质： $w^{T}w = 1$

于是正交投影公式被简化为： $\hat{x_i} = (\frac{w^{T}x_i}{w^{T}w})w = (w^{T}x_i)w = a_iw$

这里， $a_i$ 可以看成是 $\hat{x_i}$ 在 $w$ 方向上的坐标。

对包含 $n$ 个样例的数据集执行正交投影，可以从 $\{x_1, x_2,...,x_n\}$ 得到 $\{a_1, a_2,...,a_n\}$ 。

这就代表了从 $R^d\rightarrow R$ 的降维映射。

数据集简化

记数据集为 $\{x_i, y_i\}\ for\ i=0,...n$ ，为了简化处理：

假设样例的特征维数为 2，即 $x_i \in R^2$

假设样例类别数也为 2，即 $y_i = class1\ or\ class2$

投影点均值

针对数据集中属于 $class1$ 的样例，

记其样例集合为 $D_1$ ，样例个数为 $n1$ ，投影均值为 $m_1$ ，样本均值为 $\mu_1$ ，则：

$m_1 = \frac{1}{n_1}\sum a_i$

$=\frac{1}{n_1}\sum_{x_i\in D_1}^{}w^Tx_i$

$= w^T(\frac{1}{n_1}\sum_{x_i \in D_1}^{}x_i)$

$= w^T\mu_1$

相应的，针对数据集中属于 $class2$ 的样例，有：

$m_2 = w^T\mu_2$

抽象点

将针对 “不同类别的样例的投影点尽可能远离” 的数学刻画，描述为最大化两个类别的样例的投影点的均值差，即最大化下式：

$\left | m_1 - m_2 \right |$

将针对 "同类别的样例的投影点尽可能接近" 的数学刻画，描述为最小化类别内部的样例的投影点的方差，即最小化下式子：

$s_i^2 = \sum_{j=1}^{n_i}(a_j-m_i)^2$

Fisher's LDA

我们将上述的这两个思想：

1. 最大化类间投影点均值之间的差值。

2. 最小化类内投影点的方差。

融入到一个定义中，便得到了 Fisher's LDA :

$\underset{w}{max}J(w) = \frac{(m_1-m_2)^2}{s_1+s_2}$

现在，我们的目标转换为：找到一个向量 $w$ ，使得目标函数 $J(w)$ 取得最大值。

所找到的向量 $w$ 被称为最优线性判别 ( Optimal Linear Discriminant )

将样本数据代入目标函数

$(m_1-m_2)^2 = (w^T\mu_1 - w^T\mu_2)$

$= [w^T(\mu_1-\mu_2)]^2$

$= w^T(\mu_1-\mu_2)(\mu_1-\mu_2)^Tw$

$= w^TBw$

这里的 $B$ 被称为类间散度矩阵（Between class scatter matrix），维度为 $d*d$

$s_1^2 = \sum_{j=1}^{n_1}(a_j-m_1)^2$

$= \sum_{x_i \in D_1}^{}(w^Tx_i-w^T\mu_1)^2$

$= \sum_{x_i \in D_1}^{} (w^T(x_i - \mu_1))^2$

$= w^T(\sum_{x_i \in D_i}^{}(x_i - \mu_1)(x_i - \mu_1)^T)w$

$= w^TS_1w$

这里的 $S1$ 被称为 $D_1$ 的散度矩阵，

同样由 $s_2^2 = w^TS_2w$ ，可得：

$s_1^2 + s_2^2 = w^TS_1w + w^TS_2w = w^T(S_1 + S_2)w = w^TSw$

那么，LDA 的目标函数可以重写为：

$\underset{w}{max}J(w) = \frac{(m_1-m_2)^2}{s_1+s_2} = \frac{w^TBw}{w^TSw}$

目标函数求解

目标函数对 $w$ 求导，令求导后的式子为零，便可求得最优的 $w$ 值。

1. 回顾一下分数的微分公式：

$\frac{d}{dx}(\frac{f(x)}{g(x)}) = \frac{f'(x)g(x) - g'(x)f(x)}{g(x)^2}$

2. 回顾一下矩阵的微分公式：

$\frac{\partial a^Tx}{\partial x} = \frac{\partial x^Ta}{\partial x} = a$

$\frac{\partial x^TAx}{\partial x} = (A+A^T)x$

3. 同时我们发现

$B^T = B$ ， $S^T = S$

所以，令：

$\frac{d}{dw}J(w) = \frac{d}{dw}(\frac{w^TBw}{w^TSw})$

$= \frac{((B+B^T)w)(w^TSw) - ((S+S^T)w)(w^TBw)}{(w^TSw)^2}$

$= \frac{(2Bw)(w^TSw) - (2Sw)(w^TBw)}{(w^TSw)^2} = 0$

即得到：

$Bw(w^TSw) = Sw(w^TBw)$

$Bw = Sw(\frac{w^TBw}{w^TSw})$

$Bw = J(w)Sw$

注意到，这里 $J(w)$ 是个 1*1 的值，而不是一个矩阵，得到：

$Bw = \lambda Sw$

如果 $S$ 是一个可逆矩阵，那么得到：

$S^{-1}Bw = \lambda S^{-1}Sw$

$(S^{-1}B)w = \lambda w$

观察该等式的形式，我们已经将问题转化成了特征值、特征向量的分解问题。

特征向量为 $w$ ，对应的特征值为 $J(w)$ 。

我们要做的，就是对矩阵 $S^{-1}B$ 进行特征值分解（Eigen Decomposition），然后找到最大特征值所对应的那个特征向量，即为我们要找的最优向量 $w$

针对多分类

相对于二分类而言，多分类在类内散度矩阵 S 上没有变化，只在类间散度矩阵（Between class scatter matrix） $B$ 上有如下变化：

$B = \sum_{i=1}^{C}n_i(\mu_i-\mu)(\mu_i-\mu)^T$

这里：

$C$ 为多分类的类别数

$n_i$ 为单个类别的样本数

$\mu_i$ 为单个类别的样本均值

$\mu$ 为所有类别的总体样本均值

代码实现

import numpy as np
from sklearn import datasets
from sklearn import preprocessing
import matplotlib.pyplot as plt

class LDA:
    
    """
    入口函数
    @param X: The type of X is np.array 
    @param y: The type of y is np.array
    """
    def fit(self, X, y):
        self.X = X
        self.y = y
        self.target_classes = np.unique(y)
        # 每个类别的样本均值
        self.sub_mean = []
        for cls in self.target_classes:
            self.sub_mean.append(np.mean(X[list(y == cls)], axis=0))    
        # 所有类别的样本均值
        self.all_mean = np.mean(X, axis=0).reshape(1, X.shape[1])
        B = self.calculateB()
        S = self.calculateS()
        # 计算 S 的逆矩阵，乘以 B
        S_inv_B = np.linalg.inv(S).dot(B)
        # 计算对应的特征向量
        eig_vecs = self.sortEigenVectors(S_inv_B)
        return eig_vecs
    
    
    """
    计算类间散度矩阵 B
    """
    def calculateB(self):
        B = np.zeros((X.shape[1], X.shape[1]))
        for cls, mean in enumerate(self.sub_mean):
            # 每个类别的样本数
            n = self.X[self.y == cls].shape[0]
            mui_mu = mean.reshape(1, self.X.shape[1]) - self.all_mean
            B += n * np.dot(mui_mu.T, mui_mu)
        return B
        
    
    """
    计算类内散度矩阵 S
    """
    def calculateS(self):
        si_list = []
        for cls, mean in enumerate(self.sub_mean):
            si = np.zeros((X.shape[1], X.shape[1]))
            for xi in X[self.y == cls]:
                xi_mu = (xi - mean).reshape(1, self.X.shape[1])
                si += np.dot(xi_mu.T, xi_mu)
            si_list.append(si)
        S = np.zeros((X.shape[1], X.shape[1]))
        for si in si_list:
            S += si
        return S
    
    
    """
    计算最大特征向量
    """
    def sortEigenVectors(self, M):
        # 计算特征值，特征向量
        eig_values, eig_vectors = np.linalg.eig(M)
        # 排序找到最大特征值对应的特征向量
        idx = eig_values.argsort()[::-1]
        eig_values = eig_values[idx]
        eig_vectors = eig_vectors[:,idx]
        return eig_vectors
    
    """
    初始化函数
    """
    def __init__(self):
        pass

"""
加载 Iris 数据集
"""
def load_iris_data(cols):
    iris = datasets.load_iris()
    X = iris['data']
    y = iris['target']
    return X, y

# 加载数据
X, y = load_data(None)
X = X[:,(2,3)]
lda = LDA()
eig_vecs = lda.fit(X, y)
W = eig_vecs[:, :1]

# 绘制数据
colors = ['red', 'green', 'blue']
fig, ax = plt.subplots(figsize=(10, 8))
for point, pred in zip(X, y):
    ax.scatter(point[0], point[1], color=colors[pred], alpha=0.3)
    proj = (np.dot(point, W) * W) / np.dot(W.T, W)
    ax.scatter(proj[0], proj[1], color=colors[pred], alpha=0.3)
plt.show()

结果展示：