一、LDA原理

1.LDA的基本思想

LDA是一种监督学习的降维技术，也就是说它的数据集的每个样本是有类别输出的。这点和PCA不同。PCA是不考虑样本类别输出的无监督降维技术。LDA的思想可以用一句话概括，就是“投影后类内方差最小，类间方差最大”。什么意思呢？ 我们要将数据在低维度上进行投影，投影后希望每一种类别数据的投影点尽可能的接近，而不同类别的数据的类别中心之间的距离尽可能的大。
假设我们有两类数据 分别为红色和蓝色，如下图所示，这些数据特征是二维的，我们希望将这些数据投影到一维的一条直线，让每一种类别数据的投影点尽可能的接近，而红色和蓝色数据中心之间的距离尽可能的大。

上图中国提供了两种投影方式，哪一种能更好的满足我们的标准呢？从直观上可以看出，右图要比左图的投影效果好，因为右图的黑色数据和蓝色数据各个较为集中，且类别之间的距离明显。左图则在边界处数据混杂。以上就是LDA的主要思想了，当然在实际应用中，我们的数据是多个类别的，我们的原始数据一般也是超过二维的，投影后的也一般不是直线，而是一个低维的超平面。

2.瑞利商和广义瑞利商

3.二分类LDA原理

LDA的思想是设法将样本投影到一条直线上，使得：

• 同类样本的投影点尽可能近
• 异类样本的投影点尽可能远

首先从比较简单的二类LDA入手，分析LDA原理。

现在我们回到LDA的原理上，我们在第一节说讲到了LDA希望投影后希望同一种类别数据的投影点尽可能的接近，而不同类别的数据的类别中心之间的距离尽可能的大，但是这只是一个感官的度量。现在我们首先从比较简单的二类LDA入手，严谨的分析LDA的原理。
假设我们的数据集D={(x1,y1),(x2,y2),…,((xm,ym))},其中任意样本xi为n维向量，yi∈{0,1}。我们定义Nj(j=0,1)为第j类样本的个数，Xj(j=0,1)为第j类样本的集合，而μj(j=0,1)为第j类样本的均值向量，定义Σj(j=0,1)为第j类样本的协方差矩阵（严格说是缺少分母部分的协方差矩阵）。
μj的表达式为：

Σj的表达式为：

由于是两类数据，因此我们只需要将数据投影到一条直线上即可。假设我们的投影直线是向量w,则对任意一个样本本xi,它在直线w的投影为wTxi,对于我们的两个类别的中心点μ0,μ1,在在直线w的投影为wTμ0和wTμ1。由于LDA需要让不同类别的数据的类别中心之间的距离尽可能的大，也就是我们要最大化||wTμ0−wTμ1||22,同时我们希望同一种类别数据的投影点尽可能的接近，也就是要同类样本投影点的协方差wTΣ0w和wTΣ1w尽可能的小，即最小化wTΣ0w+wTΣ1w。综上所述，我们的优化目标为：
我们一般定义类内散度矩阵Sw为：
同时定义类间散度矩阵Sb为：

这样我们的优化目标重写为：
仔细一看上式，这不就是我们的广义瑞利商嘛！这就简单了，利用我们第二节讲到的广义瑞利商的性质，我们知道我们的J(w′)最大值为矩阵S−12wSbS−12w的最大特征值，而对应的w′为S−12wSbS−12w的最大特征值对应的特征向量! 而S−1wSb的特征值和S−12wSbS−12w的特征值相同,S−1wSb的特征向量w和S−12wSbS−12w的特征向量w′满足w=S−12ww′
的关系!
注意到对于二类的时候，Sbw的方向恒平行于μ0−μ1,不妨令Sbw=λ(μ0−μ1)，将其带入：(S−1wSb)w=λw，可以得到w=S−1w(μ0−μ1)， 也就是说我们只要求出原始二类样本的均值和方差就可以确定最佳的投影方向w了。

4.多分类LDA原理

有了二类LDA的基础，我们再来看看多类别LDA的原理。
假设我们的数据集D={(x1,y1),(x2,y2),…,((xm,ym))},其中任意样本xi为n维向量，yi∈{C1,C2,…,Ck}。我们定义Nj(j=1,2…k)为第j类样本的个数Xj(j=1,2…k)为第j类样本的集合，而μj(j=1,2…k)为第j类样本的均值向量，定义Σj(j=1,2…k)为第j类样本的协方差矩阵。在二类LDA里面定义的公式可以很容易的类推到多类LDA。
由于我们是多类向低维投影，则此时投影到的低维空间就不是一条直线，而是一个超平面了。假设我们投影到的低维空间的维度为d，对应的基向量为(w1,w2,…wd)，基向量组成的矩阵为W, 它是一个n×d的矩阵。
此时我们的优化目标应该可以变成为:

其中Sb=∑j=1kNj(μj−μ)(μj−μ)T,μ为所有样本均值向量。Sw=∑j=1kSwj=∑j=1k∑x∈Xj(x−μj)(x−μj)T但是有一个问题，就是WTSbW
和WTSwW都是矩阵，不是标量，无法作为一个标量函数来优化！也就是说，我们无法直接用二类LDA的优化方法，怎么办呢？一般来说，我们可以用其他的一些替代优化目标来实现。
常见的一个LDA多类优化目标函数定义为：

其中∏diagA为A的主对角线元素的乘积，W为n×d的矩阵。
J(W)的优化过程可以转化为：
仔细观察上式最右边，这不就是广义瑞利商嘛！最大值是矩阵S−1wSb
的最大特征值,最大的d个值的乘积就是矩阵S−1wSb的最大的d个特征值的乘积,此时对应的矩阵W为这最大的d个特征值对应的特征向量张成的矩阵。
由于W是一个利用了样本的类别得到的投影矩阵，因此它的降维到的维度d最大值为k-1。为什么最大维度不是类别数k呢？因为Sb中每个μj−μ的秩为1，因此协方差矩阵相加后最大的秩为k(矩阵的秩小于等于各个相加矩阵的秩的和)，但是由于如果我们知道前k-1个μj后，最后一个μk可以由前k-1个μj线性表示，因此Sb的秩最大为k-1，即特征向量最多有k-1个。

5.LDA算法流程

LDA降维的流程:

• 输入：数据集D={(x1,y1),(x2,y2),…,((xm,ym))},其中任意样本xi为n维向量，yi∈{C1,C2,…,Ck}，降维到的维度d。
• 输出：降维后的样本集 D ′ D′ D′

计算类内散度矩阵Sw
计算类间散度矩阵Sb
计算矩阵S−1wSb
计算S−1wSb的最大的d个特征值和对应的d个特征向量(w1,w2,…wd),得到投影矩阵WW
对样本集中的每一个样本特征xi,转化为新的样本zi=WTxi
得到输出样本集D′={(z1,y1),(z2,y2),…,((zm,ym))}

6.小结

LDA降维和PCA降维有很多相似之处：
　　(1) 两者在降维时都使用了特征分解的思想
　　(2) 两者都假设数据符合高斯分布，因此LDA和PCA都不适合对非高斯分布的样本进行降维
　　相对于PCA，LDA又有所不同：
　　(1) LDA是有监督的降维方法，降维过程中可以使用类别的先验知识经验，而PCA不行
　　(2) LDA选择分类性能最好的投影方向，而PCA选择最大方差的投影方向，因此LDA有过拟合的风险
　　(3) LDA最多能降到 N − 1 N-1 N−1的维数，如果降维维度大于 N − 1 N-1 N−1，则不能使用LDA，而PCA没有这个限制
　　(4) 当样本分类信息依赖均值时LDA效果较好；依赖方差的时候PCA效果较好

二、线性分类算法

支持向量机（SVM）

支持向量机（support vector machines）是一种二分类模型，它的目的是寻找一个超平面来对样本进行分割，分割的原则是间隔最大化，最终转化为一个凸二次规划问题来求解。由简至繁的模型包括：

当训练样本线性可分时，通过硬间隔最大化，学习一个线性可分支持向量机；
当训练样本近似线性可分时，通过软间隔最大化，学习一个线性支持向量机；
当训练样本线性不可分时，通过核技巧和软间隔最大化，学习一个非线性支持向量机；

三、算法代码实现

1.python编程实现

LDA算法

线性判别式分析(Linear Discriminant Analysis, LDA)，也叫做Fisher线性判别(Fisher Linear
Discriminant
,FLD)，是模式识别的经典算法，它是在1996年由Belhumeur引入模式识别和人工智能领域的。性鉴别分析的基本思想是将高维的模式样本投影到最佳鉴别矢量空间，以达到抽取分类信息和压缩特征空间维数的效果，投影后保证模式样本在新的子空间有最大的类间距离和最小的类内距离，即模式在该空间中有最佳的可分离性。因此，它是一种有效的特征抽取方法。使用这种方法能够使投影后模式样本的类间散布矩阵最大，并且同时类内散布矩阵最小。就是说，它能够保证投影后模式样本在新的空间中有最小的类内距离和最大的类间距离，即模式在该空间中有最佳的可分离性。

处理鸢尾花数据集：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobsclass LDA():def Train(self, X, y):"""X为训练数据集，y为训练label"""X1 = np.array([X[i] for i in range(len(X)) if y[i] == 0])#array()函数：创建数组X2 = np.array([X[i] for i in range(len(X)) if y[i] == 1])# 求中心点mju1 = np.mean(X1, axis=0)  # mju1是ndrray类型mju2 = np.mean(X2, axis=0)#mean()函数：计算每一列的均值# dot(a, b, out=None) 计算矩阵乘法cov1 = np.dot((X1 - mju1).T, (X1 - mju1))cov2 = np.dot((X2 - mju2).T, (X2 - mju2))Sw = cov1 + cov2# 计算ww = np.dot(np.mat(Sw).I, (mju1 - mju2).reshape((len(mju1), 1)))# 记录训练结果self.mju1 = mju1  # 第1类的分类中心self.cov1 = cov1self.mju2 = mju2  # 第2类的分类中心self.cov2 = cov2self.Sw = Sw  # 类内散度矩阵self.w = w  # 判别权重矩阵def Test(self, X, y):"""X为测试数据集，y为测试label"""# 分类结果y_new = np.dot((X), self.w)# 计算fisher线性判别式nums = len(y)c1 = np.dot((self.mju1 - self.mju2).reshape(1, (len(self.mju1))), np.mat(self.Sw).I)c2 = np.dot(c1, (self.mju1 + self.mju2).reshape((len(self.mju1), 1)))c = 1/2 * c2  # 2个分类的中心h = y_new - c# 判别y_hat = []for i in range(nums):if h[i] >= 0:y_hat.append(0)else:y_hat.append(1)# 计算分类精度count = 0for i in range(nums):if y_hat[i] == y[i]:count += 1precise = count / nums# 显示信息print("测试样本数量:", nums)print("预测正确样本的数量:", count)print("测试准确度:", precise)return precise
if '__main__' == __name__:# 产生分类数据n_samples = 500X, y = datasets.make_classification(n_samples=n_samples, n_features=2, n_redundant=0, n_classes=2,n_informative=1, n_clusters_per_class=1, class_sep=0.5, random_state=10)# LDA线性判别分析(二分类)lda = LDA()# 60% 用作训练，40%用作测试Xtrain = X[:299, :]Ytrain = y[:299]Xtest = X[300:, :]Ytest = y[300:]lda.Train(Xtrain, Ytrain)precise = lda.Test(Xtest, Ytest)# 原始数据plt.scatter(X[:, 0], X[:, 1], marker='o', c=y)plt.xlabel("x1")plt.ylabel("x2")plt.title("Test precise:" + str(precise))plt.show()

处理月亮数据集：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_moons
class LDA():def Train(self, X, y):"""X为训练数据集，y为训练label"""X1 = np.array([X[i] for i in range(len(X)) if y[i] == 0])X2 = np.array([X[i] for i in range(len(X)) if y[i] == 1])# 求中心点mju1 = np.mean(X1, axis=0)  # mju1是ndrray类型mju2 = np.mean(X2, axis=0)# dot(a, b, out=None) 计算矩阵乘法cov1 = np.dot((X1 - mju1).T, (X1 - mju1))cov2 = np.dot((X2 - mju2).T, (X2 - mju2))Sw = cov1 + cov2# 计算ww = np.dot(np.mat(Sw).I, (mju1 - mju2).reshape((len(mju1), 1)))# 记录训练结果self.mju1 = mju1  # 第1类的分类中心self.cov1 = cov1self.mju2 = mju2  # 第1类的分类中心self.cov2 = cov2self.Sw = Sw  # 类内散度矩阵self.w = w  # 判别权重矩阵def Test(self, X, y):"""X为测试数据集，y为测试label"""# 分类结果y_new = np.dot((X), self.w)# 计算fisher线性判别式nums = len(y)c1 = np.dot((self.mju1 - self.mju2).reshape(1, (len(self.mju1))), np.mat(self.Sw).I)c2 = np.dot(c1, (self.mju1 + self.mju2).reshape((len(self.mju1), 1)))c = 1/2 * c2  # 2个分类的中心h = y_new - c# 判别y_hat = []for i in range(nums):if h[i] >= 0:y_hat.append(0)else:y_hat.append(1)# 计算分类精度count = 0for i in range(nums):if y_hat[i] == y[i]:count += 1precise = count / (nums+0.000001)# 显示信息print("测试样本数量:", nums)print("预测正确样本的数量:", count)print("测试准确度:", precise)return precise
if '__main__' == __name__:# 产生分类数据X, y = make_moons(n_samples=100, noise=0.15, random_state=42)# LDA线性判别分析(二分类)lda = LDA()# 60% 用作训练，40%用作测试Xtrain = X[:60, :]Ytrain = y[:60]Xtest = X[40:, :]Ytest = y[40:]lda.Train(Xtrain, Ytrain)precise = lda.Test(Xtest, Ytest)# 原始数据plt.scatter(X[:, 0], X[:, 1], marker='o', c=y)plt.xlabel("x1")plt.ylabel("x2")plt.title("Test precise:" + str(precise))plt.show()

2.对月亮数据集进行SVM分类

import matplotlib.pyplot as plt
from sklearn.pipeline import Pipeline
import numpy as np
import matplotlib as mpl
from sklearn.datasets import make_moons
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
# 为了显示中文
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False#rc配置或rc参数，通过rc参数可以修改默认的属性，包括窗体大小、每英寸的点数、线条宽度、颜色、样式、坐标轴、坐标和网络属性、文本、字体等。
X, y = make_moons(n_samples=100, noise=0.15, random_state=42)#生成月亮数据集
def plot_dataset(X, y, axes):#绘制图形plt.plot(X[:, 0][y==0], X[:, 1][y==0], "bs")plt.plot(X[:, 0][y==1], X[:, 1][y==1], "g^")plt.axis(axes)plt.grid(True, which='both')plt.xlabel(r"$x_1$", fontsize=20)plt.ylabel(r"$x_2$", fontsize=20, rotation=0)plt.title("月亮数据",fontsize=20)
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
plt.show()

使用SVM对月亮数据集聚类：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm
from sklearn.datasets import make_moons# 导入数据集
X,y = make_moons(n_samples=200,random_state=0,noise=0.05)h = .02  # 网格中的步长# 创建支持向量机实例，并拟合出数据
C = 1.0  # SVM正则化参数
svc = svm.SVC(kernel='linear', C=C).fit(X, y) # 线性核
rbf_svc = svm.SVC(kernel='rbf', gamma=0.7, C=C).fit(X, y) # 径向基核
poly_svc = svm.SVC(kernel='poly', degree=3, C=C).fit(X, y) # 多项式核
lin_svc = svm.LinearSVC(C=C).fit(X, y) #线性核# 创建网格，以绘制图像
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h),np.arange(y_min, y_max, h))# 图的标题
titles = ['SVC with linear kernel','LinearSVC (linear kernel)','SVC with RBF kernel','SVC with polynomial (degree 3) kernel']for i, clf in enumerate((svc, lin_svc, rbf_svc, poly_svc)):# 绘出决策边界，不同的区域分配不同的颜色plt.subplot(2, 2, i + 1) # 创建一个2行2列的图，并以第i个图为当前图plt.subplots_adjust(wspace=0.4, hspace=0.4) # 设置子图间隔Z = clf.predict(np.c_[xx.ravel(), yy.ravel()]) #将xx和yy中的元素组成一对对坐标，作为支持向量机的输入，返回一个array# 把分类结果绘制出来Z = Z.reshape(xx.shape) #(220, 280)plt.contourf(xx, yy, Z, cmap=plt.cm.Paired, alpha=0.8) #使用等高线的函数将不同的区域绘制出来# 将训练数据以离散点的形式绘制出来plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)plt.xlabel('Sepal length')plt.ylabel('Sepal width')plt.xlim(xx.min(), xx.max())plt.ylim(yy.min(), yy.max())plt.xticks(())plt.yticks(())plt.title(titles[i])plt.show()

四、参考文献

1.线性判别准则与线性分类编程实践