机器学习之线性模型

　　　　在求解机器学习算法的模型参数，即无约束优化问题时，梯度下降（Gradient Descent）是最常采用的方法之一，另一种常用的方法是最小二乘法。这里就对梯度下降法做一个完整的总结。

1. 梯度

　　　　在微积分里面，对多元函数的参数求∂偏导数，把求得的各个参数的偏导数以向量的形式写出来，就是梯度。比如函数f(x,y), 分别对x,y求偏导数，求得的梯度向量就是(∂f/∂x, ∂f/∂y)^T,简称grad f(x,y)或者▽f(x,y)。对于在点(x₀,y₀)的具体梯度向量就是(∂f/∂x₀, ∂f/∂y₀)^T.或者▽f(x₀,y₀)，如果是3个参数的向量梯度，就是(∂f/∂x, ∂f/∂y，∂f/∂z)^T,以此类推。

　　　　那么这个梯度向量求出来有什么意义呢？他的意义从几何意义上讲，就是函数变化增加最快的地方。具体来说，对于函数f(x,y),在点(x₀,y₀)，沿着梯度向量的方向就是(∂f/∂x₀, ∂f/∂y₀)^T的方向是f(x,y)增加最快的地方。或者说，沿着梯度向量的方向，更加容易找到函数的最大值。反过来说，沿着梯度向量相反的方向，也就是 -(∂f/∂x₀, ∂f/∂y₀)^T的方向，梯度减少最快，也就是更加容易找到函数的最小值。

 　　　　

2. 梯度下降与梯度上升

　　　　在机器学习算法中，在最小化损失函数时，可以通过梯度下降法来一步步的迭代求解，得到最小化的损失函数，和模型参数值。反过来，如果我们需要求解损失函数的最大值，这时就需要用梯度上升法来迭代了。

　　　　梯度下降法和梯度上升法是可以互相转化的。比如我们需要求解损失函数f(θ)的最小值，这时我们需要用梯度下降法来迭代求解。但是实际上，我们可以反过来求解损失函数 -f(θ)的最大值，这时梯度上升法就派上用场了。

　　　　下面来详细总结下梯度下降法。        

3. 梯度下降法算法详解

3.1 梯度下降的直观解释

　　　　首先来看看梯度下降的一个直观的解释。比如我们在一座大山上的某处位置，由于我们不知道怎么下山，于是决定走一步算一步，也就是在每走到一个位置的时候，求解当前位置的梯度，沿着梯度的负方向，也就是当前最陡峭的位置向下走一步，然后继续求解当前位置梯度，向这一步所在位置沿着最陡峭最易下山的位置走一步。这样一步步的走下去，一直走到觉得我们已经到了山脚。当然这样走下去，有可能我们不能走到山脚，而是到了某一个局部的山峰低处。

　　　　从上面的解释可以看出，梯度下降不一定能够找到全局的最优解，有可能是一个局部最优解。当然，如果损失函数是凸函数，梯度下降法得到的解就一定是全局最优解。

3.2 梯度下降的相关概念

　　　　在详细了解梯度下降的算法之前，我们先看看相关的一些概念。

　　　　1. 步长（Learning rate）：步长决定了在梯度下降迭代的过程中，每一步沿梯度负方向前进的长度。用上面下山的例子，步长就是在当前这一步所在位置沿着最陡峭最易下山的位置走的那一步的长度。

　　　　2.特征（feature）：指的是样本中输入部分，比如样本（x₀,y₀）,（x₁,y₁）,则样本特征为x，样本输出为y。

　　　　3. 假设函数（hypothesis function）：在监督学习中，为了拟合输入样本，而使用的假设函数，记为h_θ(x)。比如对于样本（x_i,y_i）(i=1,2,...n),可以采用拟合函数如下： h_θ(x) = θ₀+θ₁x。

　　　　4. 损失函数（loss function）：为了评估模型拟合的好坏，通常用损失函数来度量拟合的程度。损失函数极小化，意味着拟合程度最好，对应的模型参数即为最优参数。在线性回归中，损失函数通常为样本输出和假设函数的差取平方。比如对于样本（x_i,y_i）(i=1,2,...n),采用线性回归，损失函数为：

              J(θ0,θ1)=∑i=1m(hθ(xi)−yi)2 J(θ0,θ1)=∑i=1m(hθ(xi)−yi)2

 　　　　其中 xi xi表示样本特征x的第i个元素， yi yi表示样本输出y的第i个元素， hθ(xi) hθ(xi)为假设函数。   

3.3 梯度下降的详细算法

　　　　梯度下降法的算法可以有代数法和矩阵法（也称向量法）两种表示，如果对矩阵分析不熟悉，则代数法更加容易理解。不过矩阵法更加的简洁，且由于使用了矩阵，实现逻辑更加的一目了然。这里先介绍代数法，后介绍矩阵法。

 

3.3.1 梯度下降法的代数方式描述

　　　　1. 先决条件： 确认优化模型的假设函数和损失函数。

　　　　比如对于线性回归，假设函数表示为  hθ(x1,x2,...xn)=θ0+θ1x1+...+θnxn hθ(x1,x2,...xn)=θ0+θ1x1+...+θnxn, 其中 θi θi (i = 0,1,2... n)为模型参数， xi xi (i = 0,1,2... n)为每个样本的n个特征值。这个表示可以简化，我们增加一个特征 x0=1 x0=1 ，这样 hθ(x0,x1,...xn)=∑i=0nθixi hθ(x0,x1,...xn)=∑i=0nθixi。

　　　　同样是线性回归，对应于上面的假设函数，损失函数为：

            J(θ0,θ1...,θn)=12m∑i=0m(hθ(x0,x1,...xn)−yi)2 J(θ0,θ1...,θn)=12m∑i=0m(hθ(x0,x1,...xn)−yi)2

 

　　　　2. 算法相关参数初始化：主要是初始化 θ0,θ1...,θn θ0,θ1...,θn,算法终止距离 ε ε以及步长 α α。在没有任何先验知识的时候，我喜欢将所有的 θ θ初始化为0， 将步长初始化为1。在调优的时候再 优化。

　　　　3. 算法过程：

　　　　　　1）确定当前位置的损失函数的梯度，对于 θi θi,其梯度表达式如下：

　　　　　　　　 ∂∂θiJ(θ0,θ1...,θn) ∂∂θiJ(θ0,θ1...,θn)

　　　　　　2）用步长乘以损失函数的梯度，得到当前位置下降的距离，即 α∂∂θiJ(θ0,θ1...,θn) α∂∂θiJ(θ0,θ1...,θn)对应于前面登山例子中的某一步。

　　　　　　3）确定是否所有的 θi θi,梯度下降的距离都小于 ε ε，如果小于 ε ε则算法终止，当前所有的 θi θi(i=0,1,...n)即为最终结果。否则进入步骤4.

　　　　　　4）更新所有的 θ θ，对于 θi θi，其更新表达式如下。更新完毕后继续转入步骤1.

　　　　　　　　 θi=θi−α∂∂θiJ(θ0,θ1...,θn) θi=θi−α∂∂θiJ(θ0,θ1...,θn)

　　　　下面用线性回归的例子来具体描述梯度下降。假设我们的样本是 (x(0)1,x(0)2,...x(0)n,y0),(x(1)1,x(1)2,...x(1)n,y1),...(x(m)1,x(m)2,...x(m)n,yn) (x1(0),x2(0),...xn(0),y0),(x1(1),x2(1),...xn(1),y1),...(x1(m),x2(m),...xn(m),yn),损失函数如前面先决条件所述：

　　　　 J(θ0,θ1...,θn)=12m∑i=0m(hθ(x0,x1,...xn)−yi)2 J(θ0,θ1...,θn)=12m∑i=0m(hθ(x0,x1,...xn)−yi)2。

　　　　则在算法过程步骤1中对于 θi θi 的偏导数计算如下： 　　

 　　　　 ∂∂θiJ(θ0,θ1...,θn)=1m∑j=0m(hθ(xj0,xj1,...xjn)−yj)xji ∂∂θiJ(θ0,θ1...,θn)=1m∑j=0m(hθ(x0j,x1j,...xnj)−yj)xij

　　　　由于样本中没有 x0 x0上式中令所有的 xj0 x0j为1.

　　　　步骤4中 θi θi的更新表达式如下：

            θi=θi−α1m∑j=0m(hθ(xj0,xj1,...xjn)−yj)xji θi=θi−α1m∑j=0m(hθ(x0j,x1j,...xnj)−yj)xij

　　　　从这个例子可以看出当前点的梯度方向是由所有的样本决定的，加 1m 1m 是为了好理解。由于步长也为常数，他们的乘机也为常数，所以这里 α1m α1m可以用一个常数表示。

　　　　在下面第4节会详细讲到的梯度下降法的变种，他们主要的区别就是对样本的采用方法不同。这里我们采用的是用所有样本。

3.3.2 梯度下降法的矩阵方式描述

　　　　这一部分主要讲解梯度下降法的矩阵方式表述，相对于3.3.1的代数法，要求有一定的矩阵分析的基础知识，尤其是矩阵求导的知识。

　　　　1. 先决条件： 和3.3.1类似， 需要确认优化模型的假设函数和损失函数。对于线性回归，假设函数 hθ(x1,x2,...xn)=θ0+θ1x1+...+θnxn hθ(x1,x2,...xn)=θ0+θ1x1+...+θnxn的矩阵表达方式为：

　　　　　 hθ(x)=Xθ hθ(x)=Xθ ，其中， 假设函数 hθ(X) hθ(X)为mx1的向量, θ θ为nx1的向量，里面有n个代数法的模型参数。 X X为mxn维的矩阵。m代表样本的个数，n代表样本的特征数。

             损失函数的表达式为： J(θ)=12(Xθ−Y)T(Xθ−Y) J(θ)=12(Xθ−Y)T(Xθ−Y), 其中 Y Y是样本的输出向量，维度为mx1.

　　　　2. 算法相关参数初始化:  θ θ向量可以初始化为默认值，或者调优后的值。算法终止距离 ε ε，步长 α α和3.3.1比没有变化。

　　　　3. 算法过程：

　　　　　　1）确定当前位置的损失函数的梯度，对于 θ θ向量,其梯度表达式如下：

　　　　　　　　 ∂∂θJ(θ) ∂∂θJ(θ)

　　　　　　2）用步长乘以损失函数的梯度，得到当前位置下降的距离，即 α∂∂θJ(θ) α∂∂θJ(θ)对应于前面登山例子中的某一步。

　　　　　　3）确定 θ θ向量里面的每个值,梯度下降的距离都小于 ε ε，如果小于 ε ε则算法终止，当前 θ θ向量即为最终结果。否则进入步骤4.

　　　　　　4）更新 θ θ向量，其更新表达式如下。更新完毕后继续转入步骤1.

　　　　　　　　 θ=θ−α∂∂θJ(θ) θ=θ−α∂∂θJ(θ)

　　　

　　　　还是用线性回归的例子来描述具体的算法过程。

　　　　损失函数对于 θ θ向量的偏导数计算如下：

　　　　　　 ∂∂θJ(θ)=XT(Xθ−Y) ∂∂θJ(θ)=XT(Xθ−Y)

　　　　步骤4中 θ θ向量的更新表达式如下： θ=θ−αXT(Xθ−Y) θ=θ−αXT(Xθ−Y)

　　　　对于3.3.1的代数法，可以看到矩阵法要简洁很多。这里面用到了矩阵求导链式法则，和两个矩阵求导的公式。

　　　　　　公式1： ∂∂X(XXT)=2X ∂∂X(XXT)=2X

　　　　　　公式2： ∂∂θ(Xθ)=XT ∂∂θ(Xθ)=XT

　　　　如果需要熟悉矩阵求导建议参考张贤达的《矩阵分析与应用》一书。

 

3.4 梯度下降的算法调优

　　　　在使用梯度下降时，需要进行调优。哪些地方需要调优呢？

　　　　1. 算法的步长选择。在前面的算法描述中，我提到取步长为1，但是实际上取值取决于数据样本，可以多取一些值，从大到小，分别运行算法，看看迭代效果，如果损失函数在变小，说明取值有效，否则要增大步长。前面说了。步长太大，会导致迭代过快，甚至有可能错过最优解。步长太小，迭代速度太慢，很长时间算法都不能结束。所以算法的步长需要多次运行后才能得到一个较为优的值。

　　　　2. 算法参数的初始值选择。 初始值不同，获得的最小值也有可能不同，因此梯度下降求得的只是局部最小值；当然如果损失函数是凸函数则一定是最优解。由于有局部最优解的风险，需要多次用不同初始值运行算法，关键损失函数的最小值，选择损失函数最小化的初值。

　　　　3.归一化。由于样本不同特征的取值范围不一样，可能导致迭代很慢，为了减少特征取值的影响，可以对特征数据归一化，也就是对于每个特征x，求出它的期望 x⎯⎯⎯ x¯和标准差std(x)，然后转化为：

　　　　　　 x−x⎯⎯⎯std(x) x−x¯std(x)

　　　　这样特征的新期望为0，新方差为1，迭代次数可以大大加快。

4. 梯度下降法大家族（BGD，SGD，MBGD）

4.1 批量梯度下降法（Batch Gradient Descent）

　　　　批量梯度下降法，是梯度下降法最常用的形式，具体做法也就是在更新参数时使用所有的样本来进行更新，这个方法对应于前面3.3.1的线性回归的梯度下降算法，也就是说3.3.1的梯度下降算法就是批量梯度下降法。　　

　　　　 θi=θi−α∑j=0m(hθ(xj0,xj1,...xjn)−yj)xji θi=θi−α∑j=0m(hθ(x0j,x1j,...xnj)−yj)xij

　　　　由于我们有m个样本，这里求梯度的时候就用了所有m个样本的梯度数据。

4.2 随机梯度下降法（Stochastic Gradient Descent）

　　　　随机梯度下降法，其实和批量梯度下降法原理类似，区别在与求梯度时没有用所有的m个样本的数据，而是仅仅选取一个样本j来求梯度。对应的更新公式是：

　　　　 θi=θi−α(hθ(xj0,xj1,...xjn)−yj)xji θi=θi−α(hθ(x0j,x1j,...xnj)−yj)xij

　　　　随机梯度下降法，和4.1的批量梯度下降法是两个极端，一个采用所有数据来梯度下降，一个用一个样本来梯度下降。自然各自的优缺点都非常突出。对于训练速度来说，随机梯度下降法由于每次仅仅采用一个样本来迭代，训练速度很快，而批量梯度下降法在样本量很大的时候，训练速度不能让人满意。对于准确度来说，随机梯度下降法用于仅仅用一个样本决定梯度方向，导致解很有可能不是最优。对于收敛速度来说，由于随机梯度下降法一次迭代一个样本，导致迭代方向变化很大，不能很快的收敛到局部最优解。

　　　　那么，有没有一个中庸的办法能够结合两种方法的优点呢？有！这就是4.3的小批量梯度下降法。

4.3 小批量梯度下降法（Mini-batch Gradient Descent）

　　小批量梯度下降法是批量梯度下降法和随机梯度下降法的折衷，也就是对于m个样本，我们采用x个样子来迭代，1<x<m。一般可以取x=10，当然根据样本的数据，可以调整这个x的值。对应的更新公式是：

　　　　 θi=θi−α∑j=tt+x−1(hθ(xj0,xj1,...xjn)−yj)xji θi=θi−α∑j=tt+x−1(hθ(x0j,x1j,...xnj)−yj)xij

5. 梯度下降法和其他无约束优化算法的比较

　　　　在机器学习中的无约束优化算法，除了梯度下降以外，还有前面提到的最小二乘法，此外还有牛顿法和拟牛顿法。

　　　　梯度下降法和最小二乘法相比，梯度下降法需要选择步长，而最小二乘法不需要。梯度下降法是迭代求解，最小二乘法是计算解析解。如果样本量不算很大，且存在解析解，最小二乘法比起梯度下降法要有优势，计算速度很快。但是如果样本量很大，用最小二乘法由于需要求一个超级大的逆矩阵，这时就很难或者很慢才能求解解析解了，使用迭代的梯度下降法比较有优势。

　　　　梯度下降法和牛顿法/拟牛顿法相比，两者都是迭代求解，不过梯度下降法是梯度求解，而牛顿法/拟牛顿法是用二阶的海森矩阵的逆矩阵或伪逆矩阵求解。相对而言，使用牛顿法/拟牛顿法收敛更快。但是每次迭代的时间比梯度下降法长。

线性模型基本形式

3.1基本形式

　　对于给定d个属性描述的示例x=（x1，x2，......，xd）,通过属性的线性组合来进行预测。一般的写法如下：

　　因此，线性模型具有很好的解释性（understandability，comprehensibility），参数w代表每个属性在回归过程中的重要程度。

3.2 线性回归

　　对于线性回归，我们先考虑简单的问题，输入的属性数目只有一个。

　　　　对于线性回归而言，

　　均方误差有非常好的几何意义，它对应了常用的欧几里得距离（欧式距离），

　　基于均方误差最小化来进行模型求解的方法称“最小二乘法”

　　在求解时，我们考虑XTX可能不满秩，因此将对应多个接都能使得均方误差最小化，选择哪个解作为输出，将由学习算法的偏好决定，最常见的方法是引入正则化。

　　广义线性回归，其中函数g（.）称为“联系函数”

线性判别分析LDA原理总结

　　　　在主成分分析（PCA）原理总结中，我们对降维算法PCA做了总结。这里我们就对另外一种经典的降维方法线性判别分析（Linear Discriminant Analysis, 以下简称LDA）做一个总结。LDA在模式识别领域（比如人脸识别，舰艇识别等图形图像识别领域）中有非常广泛的应用，因此我们有必要了解下它的算法原理。

　　　　在学习LDA之前，有必要将其自然语言处理领域的LDA区别开来，在自然语言处理领域， LDA是隐含狄利克雷分布（Latent Dirichlet Allocation，简称LDA），他是一种处理文档的主题模型。我们本文只讨论线性判别分析，因此后面所有的LDA均指线性判别分析。

1. LDA的思想

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

　　　　可能还是有点抽象，我们先看看最简单的情况。假设我们有两类数据 分别为红色和蓝色，如下图所示，这些数据特征是二维的，我们希望将这些数据投影到一维的一条直线，让每一种类别数据的投影点尽可能的接近，而红色和蓝色数据中心之间的距离尽可能的大。

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

　　　　在我们将上面直观的内容转化为可以度量的问题之前，我们先了解些必要的数学基础知识，这些在后面讲解具体LDA原理时会用到。

2. 瑞利商（Rayleigh quotient）与广义瑞利商（genralized Rayleigh quotient） 

　　　　我们首先来看看瑞利商的定义。瑞利商是指这样的函数 R(A,x) R(A,x):

R(A,x)=xHAxxHx R(A,x)=xHAxxHx

　　　　其中 x x为非零向量，而 A A为 n×n n×n的Hermitan矩阵。所谓的Hermitan矩阵就是满足共轭转置矩阵和自己相等的矩阵，即 AH=A AH=A。如果我们的矩阵A是实矩阵，则满足 AT=A AT=A的矩阵即为Hermitan矩阵。

　　　　瑞利商 R(A,x) R(A,x)有一个非常重要的性质，即它的最大值等于矩阵 A A最大的特征值，而最小值等于矩阵 A A的最小的特征值，也就是满足

λmin≤xHAxxHx≤λmax λmin≤xHAxxHx≤λmax

　　　　具体的证明这里就不给出了。当向量 x x是标准正交基时，即满足 xHx=1 xHx=1时，瑞利商退化为： R(A,x)=xHAx R(A,x)=xHAx，这个形式在谱聚类和PCA中都有出现。

　　　　以上就是瑞利商的内容，现在我们再看看广义瑞利商。广义瑞利商是指这样的函数 R(A,B,x) R(A,B,x):

R(A,x)=xHAxxHBx R(A,x)=xHAxxHBx

　　　　其中 x x为非零向量，而 A,B A,B为 n×n n×n的Hermitan矩阵。 B B为正定矩阵。它的最大值和最小值是什么呢？其实我们只要通过将其通过标准化就可以转化为瑞利商的格式。我们令 x′=B−1/2x x′=B−1/2x,则分母转化为：

xHBx=x′H(B−1/2)HBB−1/2x′=x′HB−1/2BB−1/2x′=x′Hx′ xHBx=x′H(B−1/2)HBB−1/2x′=x′HB−1/2BB−1/2x′=x′Hx′

　　　　而分子转化为：

xHAx=x′HB−1/2AB−1/2x′ xHAx=x′HB−1/2AB−1/2x′

　　　　此时我们的 R(A,B,x) R(A,B,x)转化为 R(A,B,x′) R(A,B,x′):

R(A,B,x′)=x′HB−1/2AB−1/2x′x′Hx′ R(A,B,x′)=x′HB−1/2AB−1/2x′x′Hx′

　　　　利用前面的瑞利商的性质，我们可以很快的知道， R(A,B,x) R(A,B,x)的最大值为矩阵 B−1/2AB−1/2 B−1/2AB−1/2的最大特征值，或者说矩阵 B−1A B−1A的最大特征值，而最小值为矩阵 B−1A B−1A的最小特征值。如果你看过我写的谱聚类（spectral clustering）原理总结第6.2节的话，就会发现这里使用了一样的技巧，即对矩阵进行标准化。

3. 二类LDA原理

　　　　现在我们回到LDA的原理上，我们在第一节说讲到了LDA希望投影后希望同一种类别数据的投影点尽可能的接近，而不同类别的数据的类别中心之间的距离尽可能的大，但是这只是一个感官的度量。现在我们首先从比较简单的二类LDA入手，严谨的分析LDA的原理。

　　　　假设我们的数据集 D={(x1,y1),(x2,y2),...,((xm,ym))} D={(x1,y1),(x2,y2),...,((xm,ym))},其中任意样本 xi xi为n维向量， yi∈{0,1} yi∈{0,1}。我们定义 Nj(j=0,1) Nj(j=0,1)为第j类样本的个数， Xj(j=0,1) Xj(j=0,1)为第j类样本的集合，而 μj(j=0,1) μj(j=0,1)为第j类样本的均值向量，定义 Σj(j=0,1) Σj(j=0,1)为第j类样本的协方差矩阵（严格说是缺少分母部分的协方差矩阵）。

　　　　 μj μj的表达式为：

μj=1Nj∑x∈Xjx(j=0,1) μj=1Nj∑x∈Xjx(j=0,1)

　　　　 Σj Σj的表达式为：

Σj=∑x∈Xj(x−μj)(x−μj)T(j=0,1) Σj=∑x∈Xj(x−μj)(x−μj)T(j=0,1)

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

argmaxwJ(w)=||wTμ0−wTμ1||22wTΣ0w+wTΣ1w=wT(μ0−μ1)(μ0−μ1)TwwT(Σ0+Σ1)w argmax⏟wJ(w)=||wTμ0−wTμ1||22wTΣ0w+wTΣ1w=wT(μ0−μ1)(μ0−μ1)TwwT(Σ0+Σ1)w

　　　　我们一般定义类内散度矩阵 Sw Sw为：

Sw=Σ0+Σ1=∑x∈X0(x−μ0)(x−μ0)T+∑x∈X1(x−μ1)(x−μ1)T Sw=Σ0+Σ1=∑x∈X0(x−μ0)(x−μ0)T+∑x∈X1(x−μ1)(x−μ1)T

　　　　同时定义类间散度矩阵 Sb Sb为：

Sb=(μ0−μ1)(μ0−μ1)T Sb=(μ0−μ1)(μ0−μ1)T

　　　　这样我们的优化目标重写为：

argmaxwJ(w)=wTSbwwTSww argmax⏟wJ(w)=wTSbwwTSww

　　　　仔细一看上式，这不就是我们的广义瑞利商嘛！这就简单了，利用我们第二节讲到的广义瑞利商的性质，我们知道我们的 J(w) J(w)最大值为矩阵 S−1wSb Sw−1Sb的最大特征值，而对应的 w w为 S−1wSb Sw−1Sb的最大特征值对应的特征向量!

　　　　注意到对于二类的时候， Sbw Sbw的方向恒为 μ0−μ1 μ0−μ1,不妨令 Sbw=λ(μ0−μ1) Sbw=λ(μ0−μ1)，将其带入： (S−1wSb)w=λw (Sw−1Sb)w=λw，可以得到 w=S−1w(μ0−μ1) w=Sw−1(μ0−μ1)， 也就是说我们只要求出原始二类样本的均值和方差就可以确定最佳的投影方向 w w了。

4. 多类LDA原理

　　　　有了二类LDA的基础，我们再来看看多类别LDA的原理。

　　　　假设我们的数据集 D={(x1,y1),(x2,y2),...,((xm,ym))} D={(x1,y1),(x2,y2),...,((xm,ym))},其中任意样本 xi xi为n维向量， yi∈{C1,C2,...,Ck} yi∈{C1,C2,...,Ck}。我们定义 Nj(j=1,2...k) Nj(j=1,2...k)为第j类样本的个数， Xj(j=1,2...k) Xj(j=1,2...k)为第j类样本的集合，而 μj(j=1,2...k) μj(j=1,2...k)为第j类样本的均值向量，定义 Σj(j=1,2...k) Σj(j=1,2...k)为第j类样本的协方差矩阵。在二类LDA里面定义的公式可以很容易的类推到多类LDA。

　　　　由于我们是多类向低维投影，则此时投影到的低维空间就不是一条直线，而是一个超平面了。假设我们投影到的低维空间的维度为d，对应的基向量为 (w1,w2,...wd) (w1,w2,...wd)，基向量组成的矩阵为 W W, 它是一个 m×d m×d的矩阵。

　　　　此时我们的优化目标应该可以变成为:

WTSbWWTSwW WTSbWWTSwW

　　　　其中 Sb=∑j=1kNj(μj−μ)(μj−μ)T Sb=∑j=1kNj(μj−μ)(μj−μ)T, μ μ为所有样本均值向量。 Sw=∑j=1kSwj=∑j=1k∑x∈Xj(x−μj)(x−μj)T Sw=∑j=1kSwj=∑j=1k∑x∈Xj(x−μj)(x−μj)T

　　　　但是有一个问题，就是 WTSbW WTSbW和 WTSwW WTSwW都是矩阵，不是标量，无法作为一个标量函数来优化！也就是说，我们无法直接用二类LDA的优化方法，怎么办呢？一般来说，我们可以用其他的一些替代优化目标来实现。

　　　　常见的一个LDA多类优化目标函数定义为：

argmaxWJ(W)=∏diagWTSbW∏diagWTSwW argmax⏟WJ(W)=∏diagWTSbW∏diagWTSwW

　　　　其中 ∏diagA ∏diagA为 A A的主对角线元素的乘积， W W为 m×d m×d的矩阵。

 　　　　 J(W) J(W)的优化过程可以转化为：

J(W)=∏i=1dwTiSbwi∏i=1dwTiSwwi=∏i=1dwTiSbwiwTiSwwi J(W)=∏i=1dwiTSbwi∏i=1dwiTSwwi=∏i=1dwiTSbwiwiTSwwi

　　　　仔细观察上式最右边，这不就是广义瑞利商嘛！最大值是矩阵 S−1wSb Sw−1Sb的最大特征值,最大的d个值的乘积就是矩阵 S−1wSb Sw−1Sb的最大的d个特征值的乘积,此时对应的矩阵 W W为这最大的d个特征值对应的特征向量张成的矩阵。

　　　　由于 W W是一个利用了样本的类别得到的投影矩阵，因此它的降维到的维度d最大值为k-1。为什么最大维度不是类别数k呢？因为 Sb Sb中每个 μj−μ μj−μ的秩为1，因此协方差矩阵相加后最大的秩为k(矩阵的秩小于等于各个相加矩阵的秩的和)，但是由于如果我们知道前k-1个 μj μj后，最后一个 μk μk可以由前k-1个 μj μj线性表示，因此 Sb Sb的秩最大为k-1，即特征向量最多有k-1个。

5. LDA算法流程

　　　　在第三节和第四节我们讲述了LDA的原理，现在我们对LDA降维的流程做一个总结。

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

　　　　输出：降维后的样本集$D′$

　　　　1) 计算类内散度矩阵 Sw Sw

　　　　2) 计算类间散度矩阵 Sb Sb

　　　　3) 计算矩阵 S−1wSb Sw−1Sb

　　　　4）计算 S−1wSb Sw−1Sb的最大的d个特征值和对应的d个特征向量 (w1,w2,...wd) (w1,w2,...wd),得到投影矩阵[Math Processing Error]W

　　　　5) 对样本集中的每一个样本特征 xi xi,转化为新的样本 zi=WTxi zi=WTxi

　　　　6) 得到输出样本集 D′={(z1,y1),(z2,y2),...,((zm,ym))} D′={(z1,y1),(z2,y2),...,((zm,ym))}

 

　　　　以上就是使用LDA进行降维的算法流程。实际上LDA除了可以用于降维以外，还可以用于分类。一个常见的LDA分类基本思想是假设各个类别的样本数据符合高斯分布，这样利用LDA进行投影后，可以利用极大似然估计计算各个类别投影数据的均值和方差，进而得到该类别高斯分布的概率密度函数。当一个新的样本到来后，我们可以将它投影，然后将投影后的样本特征分别带入各个类别的高斯分布概率密度函数，计算它属于这个类别的概率，最大的概率对应的类别即为预测类别。

　　　　由于LDA应用于分类现在似乎也不是那么流行，至少我们公司里没有用过，这里我就不多讲了。

6. LDA vs PCA

　　　　LDA用于降维，和PCA有很多相同，也有很多不同的地方，因此值得好好的比较一下两者的降维异同点。

　　　　首先我们看看相同点：

　　　　1）两者均可以对数据进行降维。

　　　　2）两者在降维时均使用了矩阵特征分解的思想。

　　　　3）两者都假设数据符合高斯分布。

　　　　我们接着看看不同点：

　　　　1）LDA是有监督的降维方法，而PCA是无监督的降维方法

　　　　2）LDA降维最多降到类别数k-1的维数，而PCA没有这个限制。

　　　　3）LDA除了可以用于降维，还可以用于分类。

　　　　4）LDA选择分类性能最好的投影方向，而PCA选择样本点投影具有最大方差的方向。

　　　　这点可以从下图形象的看出，在某些数据分布下LDA比PCA降维较优。

　　　　当然，某些某些数据分布下PCA比LDA降维较优，如下图所示：

7. LDA算法小结

　　　　LDA算法既可以用来降维，又可以用来分类，但是目前来说，主要还是用于降维。在我们进行图像识别图像识别相关的数据分析时，LDA是一个有力的工具。下面总结下LDA算法的优缺点。

　　　　LDA算法的主要优点有：

　　　　1）在降维过程中可以使用类别的先验知识经验，而像PCA这样的无监督学习则无法使用类别先验知识。

　　　　2）LDA在样本分类信息依赖均值而不是方差的时候，比PCA之类的算法较优。

　　　　LDA算法的主要缺点有：

　　　　1）LDA不适合对非高斯分布样本进行降维，PCA也有这个问题。

　　　　2）LDA降维最多降到类别数k-1的维数，如果我们降维的维度大于k-1，则不能使用LDA。当然目前有一些LDA的进化版算法可以绕过这个问题。

　　　　3）LDA在样本分类信息依赖方差而不是均值的时候，降维效果不好。

　　　　4）LDA可能过度拟合数据。