Fisher线性判别

分类问题

现实世界中，我们常常需要对事物进行分类，而分类的所依据的标准往往是多样的，这尤其体现在使用电脑解决分类问题之时。

我们如果要使用电脑进行分类，首先需要的便是数据。比如：我们需要将两个不同品种的鱼进行分类，我们首先需要根据找不同的方式找出两种鱼在某些特征上的不一致之处（如：长度，宽度，鱼鳍数目、长度等等）。我们将找出来的这些特征构建为一个行向量，同时对每一类都选许多条鱼进行测量，每一条鱼都能得到这样个行向量，将所有行向量拼接起来，构成一个$N\times M$ 的矩阵（$N$ 为鱼的数目，$M$ 为依据的特征数目），这样就构成了一个简单的有标签的数据集。

有了数据集，接下来所要做的便是分类。什么是分类呢？现实生活中此概念很简单，比如我们买鱼时便有许多选择：鲤鱼，草鱼，带鱼等等等等，而用以区分这些鱼类的，大都是凭经验根据某些特征来区分。其实电脑也一样，分类要做的，就是将两类数据通过一个分类面区分开来，如下粗糙的图所示，分类所做的就是根据已有的数据，找出中间那条红线，将二者分开。

而所谓分类的方法，其归根到底就是利用各种方法，来找到这样一个最佳分类面。而这篇博客所介绍的，仅是所有方法中非常简单的一种——Fisher线性判别。

Fisher线性判别

Fisher 线性判别是一种简单的二分类方法，其主要思路就是将高维空间的数据进行降维投影至一维空间，从而降低分类的难度。其一般步骤如下：

将高维数据投影到一维平面上
通过运算使得每一类之间更加紧凑（使同一类数据之间的距离变小），两类之间的距离尽可能远，从而方便将两类数据分开。
一般以两类数据中心点的中点最为分类面。

数学原理

显然，Fisher线性判别的重难点在于前两步，即如何找到一个投影面，使得类内距离尽可能小而类间距离尽可能大。

基础概念：

在 d 维的 X 空间中（样本空间）：

各类样本的均值向量$m_i$ ：$m_i = \frac{1}{n_i}\sum_\limits{X\in D_i}{X},\ i=1,2$
样本类内离散矩阵$S_i$和总样本类内离散度矩阵$S_w$ （其实可以看出来，所谓类内离散度矩阵便是每一个样本与样本均值向量之间的欧氏距离所构成的矩阵） \[ \begin{align*} & S_i = \sum_{X\in D_i}(x-m_i)(x-m_i)^T,\ i = 1, 2\\ &S_w = S_1+S_2\\ \end{align*} \]
样本类间离散度矩阵 $ S_b $ （即两类样本均值之间的欧氏距离） \[ S_b = (m_1-m_2)(m_1-m_2)^T \]

在 1 维的 Y 空间中（投影空间）：

各类样本的均值 $\tilde{m_i}$ ：$\tilde{m_i} = \frac{1}{n_i}\sum_\limits{y\in D_i}{y},\ i=1,2$
样本类内离散度 $\tilde{S_i}$ 和总样本类内离散度 $\tilde{S_w}$ \[ \begin{align*} \tilde{S_i} &= \sum_{y\in D_i}(y- \tilde{m_i})(y- \tilde{m_i})^T,\ i = 1, 2\\ &=\sum_{y\in D_i}(y- \tilde{m_i})^2,\ i = 1, 2\\ \tilde{S_w} &= \tilde{S_1}+ \tilde{S_2}\\ \end{align*} \]
样本类间离散度： \[ \begin{align*} \tilde{S_b} &= (\tilde{m_1}-\tilde{m_2})(\tilde{m_1}-\tilde{m_2})^T\\ &=(\tilde{m_1} - \tilde{m_2})^2 \end{align*} \]

接下来要让投影后的类内距离尽可能小而类间距离尽可能大，也就是要使得 $\tilde{S_w}$ 尽可能小， $\tilde{S_b}$ 尽可能大，换言之也就是令 $\frac{\tilde{S_b}}{\tilde{S_w}}$ 尽可能取得最大值，我们令其比值为 $J(w)$，如下式子： \[ \begin{align*} J(w) &= \frac{\tilde{S_b}}{\tilde{S_w}}= \frac{(\tilde{m_1} - \tilde{m_2})^2}{\tilde{S_1}+ \tilde{S_2}} \end{align*} \] 又因为，由各类样本均值可推出： \[ \begin{align*} \tilde{m_i} &= \frac{1}{n_i}\sum_\limits{y\in D_i}{y}\\ &= \frac{1}{n_i}\sum_\limits{X\in D_i}{w^Tx}\\ &= {w^T}(\frac{1}{n_i}\sum_\limits{X\in D_i}{x})\\ &= w^Tm_i \end{align*} \]

则，投影样本均值之差可以展开为： \[ \begin{align*} (\tilde{m_1} - \tilde{m_2})^2 &= (w^Tm_1 - w^Tm_2)^2\\ &= w^T(m_1-m_2)(m_1-m_2)^Tw\\ &= w^TS_bw \end{align*} \] 同理，可将 $\tilde{S_i}$ 以相同的方式进行变化： \[ \begin{align*} \tilde{S_i} &= \sum_{y\in D_i}(y- \tilde{m_i})^2\\ &= \sum_{X\in D_i}{(w^Tx-w^Tm_i)^2}\\ &= \sum_{X\in D_i}{w^T(x-m_2)(x-m_2)^Tw}\\ &= w^TS_iw \end{align*} \] 则 $J(w)$ 可化为： \[ \begin{align*} J(w) &= \frac{(\tilde{m_1} - \tilde{m_2})^2}{\tilde{S_1}+ \tilde{S_2}}\\ &= \frac{w^TS_bw}{w^TS_ww} \end{align*} \] 之后令分母为非零常数，然后采用拉格朗日乘子法确定最佳变换向量，定义拉格朗日函数如下： \[ \begin{align*} L(w, \lambda) = w^TS_bw-\lambda (w^TS_ww-c) \end{align*} \] 求解过程如下： \[ \begin{align*} &\frac{\partial L(w,\lambda)}{\partial w} = 2S_bw-2\lambda S_ww =0\\ \\ 即：\qquad &S_bw^* = \lambda S_ww^* \qquad S^{-1}_wS_bw^* = \lambda w^*\\ \\ \end{align*} \]

\[ \begin{align*} w^* &= \frac{1}{\lambda}S^{-1}_wS_bw^*\\ &= \frac{1}{\lambda}S^{-1}(m_1-m_2)(m_1-m_2)^Tw^*\\ &= \frac{R}{\lambda}S^{-1}(m_1-m_2) \end{align*} \]

省去常数，则最后可得 $w^* = S^{-1}(m_1-m_2)$

数据便可通过向量 $w^*$ 投影至一维平面上，变成一个个点，而要将两类分开，便直接可以找出阈值点 $w_0$ ，将之分开。常用确定方法如下： \[ \begin{align*} &w^0 = \frac{\tilde{m_1} - \tilde{m_2}}{2}\\ &w_0 = \frac{n_1\tilde{m_1} + n_2\tilde{m_2}}{n_1+n_2} = \tilde{m}\\ &w^0 = \frac{\tilde{m_1} - \tilde{m_2}}{2}+\frac{ln[P(w_1)/P(w_2)]}{n_1+n_2-2} \end{align*} \]

实例验证

要求：在 UCI 数据集上的 Iris和 sonar 数据上验证算法的有效性。【Iris 数据3类，4维，150个数据；Sonar 数据2类，60维，208个样本】

说明：训练和测试样本有三种方式进行划分：（三选一）

将数据随机分训练和测试，多次平均求结果
k折交叉验证
留1法

下载数据集

到 UCI 官网上即可直接下载两个数据集

零碎说明

我选择的是第一种划分方式
为了简便，我使用了 MATLAB 进行编程

开始验证

首先将数据读入，并划分矩阵，此为基础问题，不再赘述
然后，利用 mean() 函数求每一类样本的均值向量

根据公式，求取训练数据的类内散度矩阵 $S_i$ 和 $S_w$ ，代码实现函数如下：

function [s] = getScatter(sample, size, m, n1, n2)
	% getScatter
    %   计算数据的类内离散度矩阵
    s = {zeros(size), zeros(size), zeros(size)};
    for i = 1:n1
        tmp = sample{i:i};
        mi = m{i:i};
        for j = 1:n2
            xj = tmp(j:j, :);
            A = (xj-mi);
            si = A' * A;
            s{i:i} = s{i:i} + si;
        end
    end
end

根据公式，求取最佳变换向量 $w^*$ 和阈值 $w_0$ ，代码实现函数如下：

function [w, w0] = bestVector(s1, s2, m1, m2)
    % bestVector 
    %   求取最佳变换向量
    sw = s1 + s2;
    w = sw \ ((m1 - m2)');
    w0 = (w'*m1' + w'*m2') / 2;
end

将样本矩阵在投影方向上投影，以方便进行检验和画图，代码实现如下：

function [D] = rearProjection(sample, size, i1, i2, w1, n)
    % rearProjection
    %   计算投影后的点
    D = {zeros(size, 1), zeros(size, 1)};
    k = 1;
    for i = [i1 i2]
        tmp = sample{i:i};
        for j = 1:n(k)
            xj = tmp(j:j, :);
            Di = w1' * xj';
            D{k:k}(j,1) = Di;
        end
        k = k+1;
    end
end

进行检验，我写了两个函数，分别检验两类分的是否正确，代码实现如下：

function [rate] = testLeft(class, i, n, w1, w0)
    % test
    %   检验剩余测试用例
    count = 0;
    for j = i:n
        xj = class(j:j, :);
        Di = w1' * xj';
        if Di < w0
            count = count + 1;
        end
    end
    rate = (count*100)/(n+1-i);
    fprintf("%d/%d = %.2f%%\n", count, (n-i+1), rate);
end

function [rate] = testRight(class, i, n, w1, w0)
    % test
    %   检验剩余测试用例
    count = 0;
    for j = i:n
        xj = class(j:j, :);
        Di = w1' * xj';
        if Di > w0
            count = count + 1;
        end
    end
    rate = (count*100)/(n+1-i);
    fprintf("%d/%d = %.2f%%\n", count, (n-i+1), rate);
end

将投影后的数据画出散点图，画图函数如下：

function [] = drawPicture(x1, x2, n, w0)
    % drawPicture
    %   画图函数，分类点为五角星
    figure
    scatter(x1, zeros(n(1), 1), 'o');
    hold on
    scatter(x2, zeros(n(2), 1), 'x');
    hold on
    scatter(w0, 0, 400, 'p');
end