费舍尔线性分辩分析(Fisher‘s Linear Discriminant Analysis, FLDA)

费舍尔线性分辩分析(Fisher’s Linear Discriminant Analysis, FLDA)

1. 问题描述

为解决两个或多个类别的分类问题，大多数机器学习（ML）算法的工作方式相同。

通常，它们采用某种形式的转换来对输入数据进行处理，以降低原始输入维度到一个新的（更小）维度。其目的是将数据投影到新的空间中。然后，在投影后，它们尝试通过找到线性分离来对数据点进行分类。例如，我们有如下数据，

对数据直接进行线性分类显然不是最佳的方法，但是如果我们将数据投影到一维空间，我们可以找到一个线性分类器，将数据分为两个类别。这就是费舍尔线性判别分析（FLDA）的基本思想。我们将数据做如下操作：

$$y=x_0^2+x_1^2$$

其中，$x_0$和$x_1$是原始数据的两个特征。我们可以看到，通过这种方式，我们将数据投影到了一维空间，然后我们可以找到一个线性分类器，将数据分为两个类别。投影后的数据如下图所示：

通常，我们要探寻一种将数据从高维向低维度转换的方式，这被称为表征学习（Representation Learning）。深度学习也是表征学习的 一种，但在深度学习中，我们不需要猜测哪种转换会导致数据的最佳表示，算法会自行解决。

但是，请记住，无论是表示学习还是手工特征，模式都是相同的。我们需要以某种方式改变数据，使其更加适用于分类任务。

2. 二分类情况

假设我们有$C_1$，$C_2$两个类别的样本，每个样本维度为$D$，样本数为$n_1$和$n_2$。我们的目标是找到一个投影矩阵$W$，将数据投影到一维空间：

$$\hat{x}=W^{\top }x$$

设新样本$\hat{x}$的维度为1，则$W$的维度为$D\times 1$。

那么，我们该如何寻找$W$呢？换句话说，我们寻找的$W$应该符合什么条件呢？

这就是费舍尔线性判别（Fisher’s Linear Discriminant, FLD）发挥作用的地方。

费舍尔提出的想法是最大化一个函数，该函数将在投影后的类均值之间产生大的分离，同时在每个类内部给出小的方差，从而最小化类之间的重叠。

换句话说，FLD选择最大化类别间分离的投影方法。为此，它最大化类别间方差与类别内方差之比。

简而言之，为了将数据投影到更小的维度并避免类别重叠，FLD保持了两个属性：

• 数据集类别间具有很大的方差。

• 数据集每个类别内部具有较小的方差。

请注意，较大的类别间方差意味着投影后的类别平均值应该尽可能远离彼此。相反，较小的类别内方差会使投影后的数据点更加接近。

计算每个类别的均值，我们可以得到：

$${\mu }_1=\frac{1}{n_1}\sum _{x\in C_1}x,{\mu }_2=\frac{1}{n_2}\sum _{x\in C_2}x$$

其中，$m_1$和$m_2$分别是$C_1$和$C_2$类的均值。经过投影后，

$${\hat{\mu }}_1=W^{\top }{\mu }_1,{\hat{\mu }}_2=W^{\top }{\mu }_2$$

其中，${\hat{\mu }}_1$和${\hat{\mu }}_2$分别是$\widehat{C_1}$和$\widehat{C_2}$的均值。

我们计算类间方差（between-class variance）得到：

$${\hat{\mu }}_1-{\hat{\mu }}_2=W^{\top }\left({\mu }_1-{\mu }_2\right)$$

类内方差（within-class variance）为：

$${\hat{s}}_1=\sum _{i\in C_1}{\left({\hat{x}}_i-{\hat{\mu }}_1\right)}^2,{\hat{s}}_2=\sum _{i\in C_2}{\left({\hat{x}}_i-{\hat{\mu }}_2\right)}^2$$

其中，${\hat{s}}_1$和${\hat{s}}_2$分别是$\widehat{C_1}$和$\widehat{C_2}$的方差。

我们的目标是最大化类间方差与类内方差之比：

$$J(W)=\frac{{\left({\hat{\mu }}_1-{\hat{\mu }}_2\right)}^2}{{\hat{s}}_1^2+{\hat{s}}_2^2}$$

在此基础上，我们可以对$J(W)$进行进一步的变换处理。

我们定义一些散度(Scatter)的度量如下：

$$\begin{gathered} S_B=\left({\mu }_1-{\mu }_2\right){\left({\mu }_1-{\mu }_2\right)}^{\top } \\ {S}_k=\sum _{x\in C_k}\left(x-{\mu }_k\right){\left(x-{\mu }_k\right)}^{\top } \\ {S}_W=S_1+S_2 \end{gathered}$$

经过一些变换，我们得到：

$$J(W)=\frac{W^{\top }{S}_{B}W}{W^{\top }{S}_{W}W}$$

我们的目标是最大化$J(W)$，我们现在对$J(W)$求导，得到：

$$\frac{\partial J(W)}{\partial W}=(W^{\top }{S}_{W}W)\frac{\partial (W^{\top }{S}_{B}W)}{\partial W}-(W^{\top }{S}_{B}W)\frac{\partial (W^{\top }{S}_{W}W)}{\partial W}$$

令上式为0，我们得到：

$$\begin{gathered} (W^{\top }{S}_{W}W)2S_{B}W-(W^{\top }{S}_{B}W)2S_{W}W=0 \\ (W^{\top }{S}_{W}W)S_{B}W=(W^{\top }{S}_{B}W)S_{W}W \end{gathered}$$

由于投影操作，我们只关心$W$的方向，上面的式子，可以去掉$(W^{\top }{S}_{B}W),(W^{\top }{S}_{W}W)$,根据$S_B$的定义，$S_{B}W$的方向与$({\mu }_1-{\mu }_2)$一致，我们可以得到：
$$W\propto S_W^{-1}({\mu }_2-{\mu }_1)$$

3. 多分类情况

由于投影不再是一个标量，我们这里假设维度为$D^{\prime }$，我们使用Scatter矩阵的行列式来获得一个标量目标函数：

有
$$J(W)=\frac{|W^{\top }{S}_{B}W|}{|W^{\top }{S}_{W}W|}$$

对于上式的目标函数，最优投影矩阵$W$的列向量是最大的$D^{\prime }$个特征向量，对应于$S_W^{-1}{S}_B$的最大的$D^{\prime }$个特征值。

对此，我们有求解公式：

$$W={\text{max}}_{D^{\prime }}(\text{eig}(S_W^{-1}{S}_B))$$

4. 代码实现

4.1 二分类情况

rom sklearn.datasets import load_iris
import numpy as np
import pandas as pd
from torchvision import datasets
import matplotlib.pyplot as plt
from collections import Counter
from numpy.linalg import pinv
import matplotlib.lines as mlines

在MNIST上进行二分类，我们选择数字0和1，代码如下：

mnist_tr = datasets.MNIST('data', train=True, download=True)
mnist_te = datasets.MNIST('data', train=False, download=True)
x_train = mnist_tr.data.numpy()
y_train = mnist_tr.targets.numpy()
x_test = mnist_te.data.numpy()
y_test = mnist_te.targets.numpy()
two_class_data = []
two_class_target = []
for x, y in zip(x_train, y_train):
  # two class data
  if y == 0 or y == 1:
    two_class_data.append(x.flatten())
    two_class_target.append(y.squeeze())
two_class_data = np.asarray(two_class_data)
two_class_target = np.asarray(two_class_target)

划分数据为两类

C1_input = []
C1_target = []

C2_input = []
C2_target = []

for i in range(len(two_class_target)):
  y = two_class_target[i]
  x = two_class_data[i]
  if y == 0:
    C1_input.append(x.flatten())
    C1_target.append(y.squeeze())
  elif y == 1:
    C2_input.append(x.flatten())
    C2_target.append(y.squeeze())
    
C1_input = np.asarray(C1_input)
C1_target = np.asarray(C1_target)

C2_input = np.asarray(C2_input)
C2_target = np.asarray(C2_target)

计算类内均值

m1 = np.mean(C1_input,axis=0)
m2 = np.mean(C2_input,axis=0)

计算类内方差

tmp = np.subtract(C1_input, m1)
a = np.dot(tmp.T, tmp)

tmp = np.subtract(C2_input, m2)
b = np.dot(tmp.T, tmp)
SW = np.add(a,b)

计算变换矩阵$W$

inv_SW = pinv(SW)
s = m2 - m1
W = np.dot(inv_SW, np.expand_dims(s,1))

投影后的数据

y = np.dot(two_class_data,W)

计算分类阈值，这里我们选择两类数据投影后的均值作为阈值

m1 = np.mean(y[C1_target==0])
m2 = np.mean(y[C2_target==1])
threshold = (m1+m2)/2

计算分类准确率

y[y<threshold] = 0
y[y>=threshold] = 1
acc = np.sum(y.squeeze()==two_class_target)/len(two_class_target)
print('acc:',acc)

4.2 多分类情况

three_class_data = {}
for x, y in zip(x_train, y_train):
  if y == 0 or y == 1 or y == 2:
    if y not in three_class_data:
      three_class_data[y] = [x.flatten()]
    else:
      three_class_data[y].append(x.flatten())
      
three_class_data[0] = np.asarray(three_class_data[0])
three_class_data[1] = np.asarray(three_class_data[1])
three_class_data[2] = np.asarray(three_class_data[2])

class DataSet:
  def __init__(self, data, targets, valid_classes=None):
    if valid_classes is None:
      self.valid_classes = np.unique(targets)
    else:
      self.valid_classes = valid_classes
     
    self.number_of_classes = len(self.valid_classes)
    self.data = self.to_dict(data,targets)
    
  def to_dict(self,data,targets):
    data_dict = {}
    for x, y in zip(data, targets):
      if y in self.valid_classes:
        if y not in data_dict:
          data_dict[y] = [x.flatten()]
        else:
          data_dict[y].append(x.flatten())
     
    for i in self.valid_classes:
      data_dict[i] = np.asarray(data_dict[i])

    return data_dict
  
  def get_data_by_class(self,class_id):
    if class_id in self.valid_classes:
      return self.data[class_id]
    else:
      raise("Class not found.")
  
  def get_all_data(self):
    data = []
    labels = []
    for label, class_i_data in self.data.items():
      data.extend(class_i_data)
      labels.extend(class_i_data.shape[0] * [label])
    data = np.asarray(data)
    labels = np.asarray(labels)
    return data, labels

dataset = DataSet(x_train, y_train, valid_classes=[0, 1, 2])

inputs, targets = dataset.get_all_data()

定义类别数量和目标维度

number_of_classes = three_class_data.keys()
D_prime = 2

计算类内均值

mk = []
for class_i, input_vectors in three_class_data.items():
  mk.append(np.mean(input_vectors,axis=0))
  mk[class_i] = np.asarray(mk[class_i])

计算类内方差

Sks = []
for (class_i, input_vectors), m in zip(three_class_data.items(),mk):
  tmp = np.subtract(input_vectors, m)
  Sks.append(np.dot(np.transpose(tmp), tmp))
  
Sks = np.asarray(Sks)

计算类间方差

N = 0
Nk = []
sum_ = 0
for class_i, data in three_class_data.items():
  Nk.append(data.shape[0])
  sum_ += np.sum(data,axis=0)

N = sum(Nk)
# m is the mean of the total data set
m = sum_ / N

SB = []
for class_i in three_class_data.keys():
  tmp = mk[class_i] - m
  SB.append(np.multiply(Nk[class_i], np.outer(tmp, tmp.T)))
SB = np.sum(SB,axis=0) # sum of K (# of classes) matrices

计算投影矩阵$W$

from numpy.linalg import eig
matrix = np.dot(pinv(Sw),SB)
print("Out:",matrix.shape)

# find eigen values and eigen-vectors pairs for np.dot(pinv(SW),SB)
eigen_values, eigen_vectors = eig(matrix)

# sort eigen values and eigen-vectors pairs
idx = eigen_values.argsort()[::-1]
eigen_values = eigen_values[idx]
eigen_vectors = eigen_vectors[:,idx]
# find the projection matrix W
W = eigen_vectors[:,:D_prime]

投影后的数据

def inference(x,W):
  y = np.dot(x,W)
  return y

yk = []
for class_i, data in three_class_data.items():
  yk.extend(inference(data,W))
y = np.asarray(yk)
print(y.shape)

5. 参考资料

• https://sthalles.github.io/fisher-linear-discriminant/

• https://www.csd.uwo.ca/~oveksler/Courses/CS434a_541a/Lecture8.pdf

• https://www.ccs.neu.edu/home/vip/teach/MLcourse/5_features_dimensions/lecture_notes/LDA/LDA.pdf

• http://webdancer.is-programmer.com/posts/37867.html