费舍尔线性分辩分析(Fisher’s Linear Discriminant Analysis, FLDA)
目录
- 费舍尔线性分辩分析(Fisher's Linear Discriminant Analysis, FLDA)
- 1. 问题描述
- 2. 二分类情况
- 3. 多分类情况
- 4. 代码实现
- 4.1 二分类情况
- 4.2 多分类情况
- 5. 参考资料
1. 问题描述
为解决两个或多个类别的分类问题,大多数机器学习(ML)算法的工作方式相同。
通常,它们采用某种形式的转换来对输入数据进行处理,以降低原始输入维度到一个新的(更小)维度。其目的是将数据投影到新的空间中。然后,在投影后,它们尝试通过找到线性分离来对数据点进行分类。例如,我们有如下数据,
对数据直接进行线性分类显然不是最佳的方法,但是如果我们将数据投影到一维空间,我们可以找到一个线性分类器,将数据分为两个类别。这就是费舍尔线性判别分析(FLDA)的基本思想。我们将数据做如下操作:
y = x 0 2 + x 1 2 y=x_{0}^2+x_{1}^2 y=x02+x12
其中, x 0 x_{0} x0和 x 1 x_{1} x1是原始数据的两个特征。我们可以看到,通过这种方式,我们将数据投影到了一维空间,然后我们可以找到一个线性分类器,将数据分为两个类别。投影后的数据如下图所示:
通常,我们要探寻一种将数据从高维向低维度转换的方式,这被称为表征学习(Representation Learning)。深度学习也是表征学习的 一种,但在深度学习中,我们不需要猜测哪种转换会导致数据的最佳表示,算法会自行解决。
但是,请记住,无论是表示学习还是手工特征,模式都是相同的。我们需要以某种方式改变数据,使其更加适用于分类任务。
2. 二分类情况
假设我们有 C 1 C_1 C1, C 2 C_2 C2两个类别的样本,每个样本维度为 D D D,样本数为 n 1 n_1 n1和 n 2 n_2 n2。我们的目标是找到一个投影矩阵 W W W,将数据投影到一维空间:
x ^ = W ⊤ x \hat{x}=W^{\top}x x^=W⊤x
设新样本 x ^ \hat{x} x^的维度为1,则 W W W的维度为 D × 1 D \times 1 D×1。
那么,我们该如何寻找 W W W呢?换句话说,我们寻找的 W W W应该符合什么条件呢?
这就是费舍尔线性判别(Fisher’s Linear Discriminant, FLD)发挥作用的地方。
费舍尔提出的想法是最大化一个函数,该函数将在投影后的类均值之间产生大的分离,同时在每个类内部给出小的方差,从而最小化类之间的重叠。
换句话说,FLD选择最大化类别间分离的投影方法。为此,它最大化类别间方差与类别内方差之比。
简而言之,为了将数据投影到更小的维度并避免类别重叠,FLD保持了两个属性:
-
数据集类别间具有很大的方差。
-
数据集每个类别内部具有较小的方差。
请注意,较大的类别间方差意味着投影后的类别平均值应该尽可能远离彼此。相反,较小的类别内方差会使投影后的数据点更加接近。
计算每个类别的均值,我们可以得到:
μ 1 = 1 n 1 ∑ x ∈ C 1 x , μ 2 = 1 n 2 ∑ x ∈ C 2 x \mu_{1}=\frac{1}{n_{1}} \sum_{x \in C_{1}} x, \quad \mu_{2}=\frac{1}{n_{2}} \sum_{x \in C_{2}} x μ1=n11x∈C1∑x,μ2=n21x∈C2∑x
其中, m 1 m_1 m1和 m 2 m_2 m2分别是 C 1 C_1 C1和 C 2 C_2 C2类的均值。经过投影后,
μ ^ 1 = W ⊤ μ 1 , μ ^ 2 = W ⊤ μ 2 \hat{\mu}_{1}=W^{\top} \mu_{1}, \quad \hat{\mu}_{2}=W^{\top} \mu_{2} μ^1=W⊤μ1,μ^2=W⊤μ2
其中, μ ^ 1 \hat{\mu}_{1} μ^1和 μ ^ 2 \hat{\mu}_{2} μ^2分别是 C 1 ^ \hat{C_1} C1^和 C 2 ^ \hat{C_2} C2^的均值。
我们计算类间方差(between-class variance)得到:
μ ^ 1 − μ ^ 2 = W ⊤ ( μ 1 − μ 2 ) \hat{\mu}_{1} - \hat{\mu}_{2}=W^{\top}\left(\mu_{1}-\mu_{2}\right) μ^1−μ^2=W⊤(μ1−μ2)
类内方差(within-class variance)为:
s ^ 1 = ∑ i ∈ C 1 ( x ^ i − μ ^ 1 ) 2 , s ^ 2 = ∑ i ∈ C 2 ( x ^ i − μ ^ 2 ) 2 \hat{s}_{1}=\sum_{i\in C_{1}}\left(\hat{x}_{i}-\hat{\mu}_{1}\right)^{2}, \quad \hat{s}_{2}=\sum_{i\in C_{2}}\left(\hat{x}_{i}-\hat{\mu}_{2}\right)^{2} s^1=i∈C1∑(x^i−μ^1)2,s^2=i∈C2∑(x^i−μ^2)2
其中, s ^ 1 \hat{s}_{1} s^1和 s ^ 2 \hat{s}_{2} s^2分别是 C 1 ^ \hat{C_1} C1^和 C 2 ^ \hat{C_2} C2^的方差。
我们的目标是最大化类间方差与类内方差之比:
J ( W ) = ( μ ^ 1 − μ ^ 2 ) 2 s ^ 1 2 + s ^ 2 2 J(W)=\frac{\left(\hat{\mu}_{1}-\hat{\mu}_{2}\right)^{2}}{\hat{s}_{1}^{2}+\hat{s}_{2}^{2}} J(W)=s^12+s^22(μ^1−μ^2)2
在此基础上,我们可以对 J ( W ) J(W) J(W)进行进一步的变换处理。
我们定义一些散度(Scatter)的度量如下:
S B = ( μ 1 − μ 2 ) ( μ 1 − μ 2 ) ⊤ S k = ∑ x ∈ C k ( x − μ k ) ( x − μ k ) ⊤ S W = S 1 + S 2 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} SB=(μ1−μ2)(μ1−μ2)⊤ Sk=x∈Ck∑(x−μk)(x−μk)⊤ SW=S1+S2
经过一些变换,我们得到:
J ( W ) = W ⊤ S B W W ⊤ S W W J(W)=\frac{W^{\top} S_{B} W}{W^{\top} S_{W} W} \\ \ \\ J(W)=W⊤SWWW⊤SBW
我们的目标是最大化 J ( W ) J(W) J(W),我们现在对 J ( W ) J(W) J(W)求导,得到:
∂ J ( W ) ∂ W = ( W ⊤ S W W ) ∂ ( W ⊤ S B W ) ∂ W − ( W ⊤ S B W ) ∂ ( W ⊤ S W W ) ∂ 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} ∂W∂J(W)=(W⊤SWW)∂W∂(W⊤SBW)−(W⊤SBW)∂W∂(W⊤SWW)
令上式为0,我们得到:
( W ⊤ S W W ) 2 S B W − ( W ⊤ S B W ) 2 S W W = 0 ( W ⊤ S W W ) S B W = ( W ⊤ S B W ) S W W (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 (W⊤SWW)2SBW−(W⊤SBW)2SWW=0 (W⊤SWW)SBW=(W⊤SBW)SWW
由于投影操作,我们只关心 W W W的方向,上面的式子,可以去掉 ( W ⊤ S B W ) , ( W ⊤ S W W ) (W^{\top}S_{B}W),(W^{\top}S_{W}W) (W⊤SBW),(W⊤SWW),根据 S B S_{B} SB的定义, S B W S_BW SBW的方向与 ( μ 1 − μ 2 ) (\mu_{1}−\mu_{2}) (μ1−μ2)一致,我们可以得到:
W ∝ S W − 1 ( μ 2 − μ 1 ) W∝S_{W}^{-1}(\mu_{2}−\mu_{1}) W∝SW−1(μ2−μ1)
3. 多分类情况
由于投影不再是一个标量,我们这里假设维度为 D ′ D^{\prime} D′,我们使用Scatter矩阵的行列式来获得一个标量目标函数:
有
J ( W ) = ∣ W ⊤ S B W ∣ ∣ W ⊤ S W W ∣ J(W)=\frac{|W^{\top} S_{B} W|}{|W^{\top} S_{W} W|} J(W)=∣W⊤SWW∣∣W⊤SBW∣
对于上式的目标函数,最优投影矩阵 W W W的列向量是最大的 D ′ D^{\prime} D′个特征向量,对应于 S W − 1 S B S_{W}^{-1}S_{B} SW−1SB的最大的 D ′ D^{\prime} D′个特征值。
对此,我们有求解公式:
W = max D ′ ( eig ( S W − 1 S B ) ) W=\text{max}_{D^{\prime}}(\text{eig}(S_{W}^{-1}S_{B})) W=maxD′(eig(SW−1SB))
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 dataif 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 W 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_classesself.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_dictdef 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, labelsdataset = 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_ / NSB = []
for class_i in three_class_data.keys():tmp = mk[class_i] - mSB.append(np.multiply(Nk[class_i], np.outer(tmp, tmp.T)))
SB = np.sum(SB,axis=0) # sum of K (# of classes) matrices
计算投影矩阵 W W 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 yyk = []
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