最优化
随着大数据的到来，并行计算的流行，实际上机器学习领域的很多研究者会把重点放在最优化方法的研究上，如large scale computation。那么为什么要研究最优化呢？我们先从机器学习研究的目的说起。机器学习理论主要是设计和分析一些让计算机可以自动“学习”的算法，这些算法可以从数据中自动分析获得规律，并利用规律对未知数据进行预测，并可用于发现数据之间隐藏的关系，解释某些现象的发生。至于为什么要让机器去做这些事情，原因很简单，数据量和维度过于庞大，无法通过人脑简单的处理或分析这些数据。比如，我们无法通过百万级的DNA序列分析各序列和疾病发生的关系，也无法在一定有限的时间内标定出1万张图像上人脸的位置。所以 研究者倾向于建立一些机器学习模型，通过输入一个样本（一个人的DNA序列或者是一副图片），输出样本的标签（是否患病、头像的位置）。这些学习模型里有大量可以调整的参数，它们通过调整参数组合，拟合高维空间训练样本数据的分布，识别出数据和标签之间复杂的关系。目前已有的神经网络、支持向量机、AdaBoost、卷积神经网络等算法，它们的共同特点是通过调整参数让模型的目标函数尽可能大或者小（如logistic回归是使得分类错误率尽量小）。为了达到该目的，不同的机器学习算法首先会设定不同的目标函数，然后给参数赋上随机初值，最后用各种下降法更快更好地寻找能让分类错误率更小的参数组合。
所以，从广义上讲，机器学习算法的本质上是优化问题求解，例如，梯度下降、牛顿法、共轭梯度法都是常见的机器学习算法优化方法。那么有人肯定会有疑问，这不还是调参数、选择参数么？这个参数优化与之前的调参的概念是不同的，之前说的调参是针对算法中的自由参数（即通过经验指定的参数，用过weka或者R的人应该知道，如SVM中的gamma或者logistic回归中的惩罚因子lamda），这些参数主要是控制学习模型的泛化能力，防止过拟合。而这里通过优化算法迭代出的参数则是目标函数中待求解的参数，例如，神经网络中各隐含层节点的权值、logistic回归中各自变量的系数。对于不同目标函数和不同优化算法，产生的问题也各不相同，比如某些目标函数不是凸函数也不是无约束的优化问题，无法直接利用梯度下降算法求解，又例如梯度下降法往往只能保证找到目标函数的局部最小值，找不到全局最小值，那么该怎么解决这些问题呢？答案是不一味的强行采用梯度下降，而是引入其他变量（拉格朗日乘子）将有约束问题转化为无约束问题，也可以适当爬爬山，说不定能跳出小水沟（局部极小值）找到真正的深井（全局极小值），这种算法叫模拟退火。也可以增大搜索范围，让一群蚂蚁（蚁群算法）或者鸟儿（粒子群算法）一齐搜索，或者让参数巧妙地随机改变（遗传算法）。所以，如何更快的找到目标函数的全局最小值或者全局最大值，如何避免陷入局部最优解是优化算法研究的重点。
讲了这么多，主要是为了说明机器学习与最优化问题的联系，也为大家更好的理解后续机器学习算法提供基础。接下来，我们会把讲解重点放在放在最优化及凸优化理论上。
1. 最优化问题
数值优化问题或者简称为优化问题主要是求问题的解，优化问题具有以下一般形式：
minimize f0(x)subject to fi(x)<bi, i=1,…,m（1.1）minimize f0(x)subject to fi(x)<bi, i=1,…,m（1.1）
其中，函数f0f0为目标函数，函数fi：Rn→Rfi：Rn→R为不等式，或约束函数。x=(x1,…,xn)x=(x1,…,xn)为向量，是目标函数的待优化参数（optimization variables），也称为优化问题的可行解，常量b1,…,bmb1,…,bm称为边界或约束。当存在向量x∗x∗，使得满足上式约束的任意向量zz，存在f0(x∗)<f0(z)f0(x∗)<f0(z)，则向量x∗x∗称为上式（1.1）的最优解。当不存在约束时，优化问题称为无约束的优化问题，反之，称为有约束的优化问题。对于有约束的优化问题，当目标函数以及约束函数fi(x), i=0,…,mfi(x), i=0,…,m满为线性函数，则称该优化问题为线性规划（linear programming）。例如，fi(x)fi(x)满足fi(αx+βy)=αfi(x)+βfi(y)fi(αx+βy)=αfi(x)+βfi(y)；如果目标函数或者约束函数不满足线性关系，则称优化问题为非线性规划（nonlinear programming）。如果目标函数和约束函数都满足凸函数不等式的性质，即：fi(αx+βy)≦αfi(x)+βfi(y)fi(αx+βy)≦αfi(x)+βfi(y)，其中α+β=1, α≧0,β≧0α+β=1, α≧0,β≧0则该优化问题可称为凸优化问题。很明显， 线性规划问题也服从凸优化问题的条件，只有αα和ββ具有一定的特殊取值时，不等式才能变成等式成立，因此，凸优化问题是线性规划问题的一般形式。
对于机器学习而言，其一般表现为从一系列潜在的模型集合中寻找一个模型能最好的满足先验知识以及拟合观测数据。模型训练的目的在于找到分类错误或者预测误差最小的模型参数值。相对于优化问题，这里所说的先验知识可以对应成优化问题的约束函数，而目标函数则是一个评价观测值和估计值差别的函数（均方误差函数），或者是评价观测数据服从某模型下某组参数值的可能性的函数（似然函数）。因此，还是像之前所说过的，机器学习其实就是一个优化问题，优化问题的相关研究进展会极大影响机器学习领域的发展。
最优化算法即用于求解优化问题的方法，啊，简直是废话。。。最优化算法的计算效率（时间复杂度和空间复杂度）会严重限制算法的应用，例如，对于非光滑函数，传统的优化算法一般不会有较好的效果；优化参数维度较高时，有些优化算法也不适用。优化问题中，最常见的是最小二乘问题和线性规划，接下来会对这两个基本的优化问题做简要介绍。
a. 最小二乘问题
最小二乘问题（Least－Square problems）是无约束优化问题，同时，目标函数是具有aTix−biaiTx−bi形式的线性表达式的平方和，其一般形式可记为：
minimize f0(x)=∥Ax−b∥2=k∑i=1(aTix−bi)2（1.2）minimize f0(x)=∥Ax−b∥2=∑i=1k(aiTx−bi)2（1.2）
其中，A∋Rk×n, k≧nA∋Rk×n, k≧n为观测样本集合，向量xx为待优化参数。
最小二乘问题是回归分析的根本，最小二乘问题很容易辨认，当目标函数为二次函数时，即可认为该优化问题为最小二乘问题。我们学过的解决该问题的最简单的方法是最小二乘法，我们可以将式(1.2)简化为求解线性等式，(ATA)x=ATb(ATA)x=ATb。因此，我们可以获得求解该问题的解析式x=(ATA)−1ATbx=(ATA)−1ATb。该方法的时间复杂度为n2kn2k，当样本数量以及特征维度较低时（百维、千维），一般计算机会在几秒内求的最优解，当使用更好的计算资源时，对于同样的样本量计算时间会呈指数衰减（摩尔定律）。但是，对于需要满足实时计算要求时，如果样本特征维度高于百万级别，采用最小二乘法求解就会变的不可行。所以，我们一般会采用梯度下降等迭代优化方法求解上述目标函数的可行解，当然为了防止过拟合，可选择惩罚（regularization）或者偏最小二乘（weighted least－squares）解决该问题。
b. 线性规划
线性规划是优化问题的另一个重要分支，其一般形式为：
miniminze cTxsubject to aTix≦bi, i=1,…,m(1.3)miniminze cTxsubject to aiTx≦bi, i=1,…,m(1.3)
对于线性规划问题，不存在类似最小二乘问题的一步求解方法，即最小二乘法，但是可用于解决线性规划问题的方法很多，如simplex方法，内插点法。虽然无法直接一步求解，但是我们可以通过控制迭代次数或设置停止迭代条件来减少运算的时间复杂度，例如，内插点法的时间复杂度为n2mn2m，其中m≧nm≧n。另外，采用迭代法求解优化问题不一定能像最小二乘法一样求得全局最优解，但目前的迭代算法大多场合下可以达到最小二乘法一样的准确度，而且，可满足实时计算的需求。
同时，很多优化问题都可以转换成线性规划问题，如Chebyshev approximation problem：
minimize maxi=1,…,k∣aTix−bi∣(1.4)minimize maxi=1,…,k∣aiTx−bi∣(1.4)
其中，x为待优化参数。Chebyshev优化问题与最小二乘问题类似，但最小二乘问题可微（矩阵半正定），而Chebyshev目标函数不可微，所以无法采用最小二乘法求解。我们需要将目标函数进行转化，即可转化成线性规划：
minimize tsubject to aTix−t≦bi, −aTix−t≦−bi, wherei=1,…,k(1.5)minimize tsubject to aiTx−t≦bi, −aiTx−t≦−bi, wherei=1,…,k(1.5)
这样，我们就可采用simplex等方法求解该对偶线性规划问题。
c. 凸优化
凸函数的定义在上面已经介绍过了，即：
minimize f0(x)subject to fi(x)<bi, i=1,…,m（1.6）minimize f0(x)subject to fi(x)<bi, i=1,…,m（1.6）
其中，函数f0,…,fm:Rn→Rf0,…,fm:Rn→R为凸函数。
凸函数定义为：
fi(tx+(1−t)y)≦tfi(x)+(1−t)fi(y)fi(tx+(1−t)y)≦tfi(x)+(1−t)fi(y)，其中t≧0（1.7）t≧0（1.7）
也就是说，凸函数其函数图像上方的点集是凸集，函数图像上的任意两点确定的弦在其图像上方，这里需要点明的是国内某些教材上关于凸函数和凹函数的定义与这里写的正好相反，这里的凸函数是直观视觉上的凹函数，即碗形函数。凸函数的定义在定义域C上的凸函数可表现为
凸函数的判定方法一般是求解其二阶导数（如果存在），如果其二阶导数在区间上非负，则该函数为凸函数。当且仅当，二阶导数在区间上恒大于0，则函数称为严格凸函数。凸函数具有如下性质：
(1) 凸函数的一阶导数在区间内单调不减；
(2) 凸函数具有仿射不变性，即f(x)f(x)是凸函数，则g(y)=f(Ax+b)g(y)=f(Ax+b)也是凸函数；
(3) 凸函数的任何 极小值都是最小值，严格凸函数最多有一个最小值；
(4) 凸函数在区间内（凸集内部）是正定的。最小二乘问题和线性规划问题都属于凸优化问题。
因为凸函数具有局部最优解就是全局最优解的优良性质，我们可以在求解过程不用过多考虑局部最优解和全局最优解的问题，因此，现有优化问题研究更多放在将一般形式的目标函数转化为凸函数后求解。而对于凸优化问题，我们可以采用熟知的内插法、梯度下降法、牛顿拉斐逊算法以及BFGS算法等。这些算法以及如何将优化目标函数转换为凸函数是本系列博客的主要阐述的问题。
数学优化入门：梯度下降法、牛顿法、共轭梯度法
2016-10-15
数学优化入门：梯度下降法、牛顿法、共轭梯度法。前面讲的都是从n维空间到一维空间的映射函数，对于从n维欧式空间变换到m维欧式空间的函数f，也可以表示成由m个实函数组成y=f(x)=[y1(x1,…xn),…ym(x1,…,xn)]T。
1、基本概念
1.1 方向导数
1.2 梯度的概念
如果考虑z=f(x,y)描绘的是一座在点(x,y)的高度为f(x,y)的山。那么，某一点的梯度方向是在该点坡度最陡的方向，而梯度的大小告诉我们坡度到底有多陡。
对于含有n个变量的标量函数，其梯度表示为
1.3 梯度与方向导数
函数在某点的梯度是这样一个向量，它的方向与取得最大方向导数的方向一致，而它的模为方向导数的最大值。
1.4 梯度与等高线
函数z=f(x)在点P(x,y)的梯度的方向与过点的等高线f(x,y)=c在这点的法线的一个方向相同，且从数值较低的等高线指向数值较高的等高线，而梯度的模等于函数在这个法线方向的方向导数。这个法线方向就是方向导数取得最大值的方向。
即负梯度方向为最速下降方向。
1.5 等高面
对于二次函数
其中A为n*n的对称正定矩阵，x-为一定点，则函数f(x)的等值面f(x,y)=c是一个以x-为中心的椭球面。
此椭球面的形状受 Hesse 矩阵的条件数影响，长轴与短轴对应矩阵的最小特征值和最大特征值的方向，其大小与特征值的平方根成反比，最大特征值与最小特征值相差越大，椭球面越扁。
矩阵的条件数：矩阵A的条件数等于A的范数与A的逆的范数的乘积，即cond(A)=‖A‖·‖A^(-1)‖，是用来判断矩阵病态与否的一种度量，条件数越大矩阵越病态。
1.6 Hesse 矩阵
在牛顿法中应用广泛，定义为
1.7 Jacobi矩阵
前面讲的都是从n维空间到一维空间的映射函数，对于从n维欧式空间变换到m维欧式空间的函数f，也可以表示成由m个实函数组成y=f(x)=[y1(x1,…xn),…ym(x1,…,xn)]T。对于函数f，如果其偏导数都存在，可以组成一个m行n列的矩阵，即所谓的Jacobi矩阵：
显然， 当f(x) 是一个标量函数时，Jacobi矩阵是一个向量，即f(x)的梯度，此时Hesse 矩阵是一个二维矩阵;当f(x)是一个向量值函数时，Jacobi 矩阵是一个二维矩阵，Hesse 矩阵是一个三维矩阵。
1.8 共轭方向
先给出共轭方向的定义：
当A为单位阵时，这两个方向关于A共轭将等价于两个方向正交，可见共轭是正交概念的推广。
我们在来看共轭方向的几何意义。
前面提到过二次函数
的等值面f(x,y)=c是一个以x-为中心的椭球面。设x^(1)为此椭球面边缘的一点，则x^(1)处的法向量为
将其中后面一项记作
即由x(1)指向椭圆面中心x-的向量。
下面，设d^(1)为此椭球面在x(1)处的切向量，由于切向量d^(1)与法向量delta f(x(1))正交，即
可见，等值面上一点处的切向量与由此点指向极小点的向量是关于A共轭的。
因此，极小化上述二次函数，若依次沿着d^(1)和d^(2)进行一维搜索，则经过两次迭代必达到极小点。
1.9 一维搜索
在许多迭代下降算法中，具有一个共同点，就是得到x(k)后，按某种规则确定一个方向d(k)，再从x(k)除法，沿方向d(k)在直线上求目标函数f(x(k)+lambda*d(k))的的极小点，从而得到x(k)的后继点x(k+1)，这里求目标函数在直线上的极小点，称为一维搜索，或者线搜索，可以归结为单变量lambda的极小化问题。确定的lambda可以成为步长。
一维搜索方法很多，大致可以分为试探法和函数逼近法(插值法)。当然，这两种方法都是求得即小的的近似值。
试探法包括0.618法、Fibonacci法、进退法、平分法等。
函数逼近法包括牛顿法、割线法、抛物线法、三次插值法、有理插值法等。
2、梯度下降法(最速下降法)
即利用一阶的梯度信息找到函数局部最优解的一种方法。核心迭代公式为
其中，pk是第k次迭代时选取的移动方向，在梯度下降法中，移动的方向设定为梯度的负方向。
ak是第k次迭代是移动的步长，每次移动的步长可以固定也可以改变。在数学上，步长可以通过line search令导数为零找到该方向上的最小值，但是在实际编程时，这样计算的代价太大，我们一般可以将它设定位一个常量。
因此，梯度下降法计算步骤可以概括为
如果目标函数是一个凸优化问题，那么梯度下降法获得的局部最优解就是全局最优解，理想的优化效果如下图，值得注意一点的是，每一次迭代的移动方向都与出发点的等高线垂直：
实际上，梯度还可以提供不在最快变化方向的其他方向上坡度的变化速度，即在二维情况下，按照梯度方向倾斜的圆在平面上投影成一个椭圆。椭球面的形状受 Hesse 矩阵的条件数影响，椭球面越扁，那么优化路径需要走很大的弯路，计算效率很低。这就是常说的锯齿现象( zig-zagging)，将会导致收算法敛速度变慢。
3、牛顿法
前面提到的梯度下降法，主要利用的是目标函数的局部性质，具有一定的“盲目性”。
牛顿法则是利用局部的一阶和二阶偏导信息，去推测整个目标函数的形状，进而可以求得近似函数的全局最小值，然后将当前的最小值设定为近似函数的最小值。也就是说，牛顿法在二阶导数的作用下，从函数的凸性出发，直接搜索怎样到达极值点，即在选择方向时，不仅考虑当前坡度是否够大，还会考虑走了一步之后，坡度是否会变得更大。
相比最速下降法，牛顿法带有一定对全局的预测性，收敛性质也更优良。当然由于牛顿法是二阶收敛的，比梯度下降法收敛的更快。
假设我们要求f(x)的最小值，首先用泰勒级数求得其二阶近似
显然这里x是自变量，x^(k)是常量，求解近似函数phi(x)的极值，即令其倒数为0，很容易得到
从而得到牛顿法的迭代公式
显然除了f(x)二次可微外，还要求f(x)的Hesse矩阵可逆。此外，由于矩阵取逆的计算复杂为 n 的立方，当问题规模比较大时，计算量很大，解决的办法是采用拟牛顿法，如 BFGS, L-BFGS, DFP, Broyden’s Algorithm 进行近似。
此外，如果初始值离局部极小值太远，泰勒展开并不能对原函数进行良好的近似，导致牛顿法可能不收敛。
我们给出阻尼牛顿法的计算步骤，其实阻尼牛顿法相较原始牛顿法只是增加了沿牛顿方向的一维搜索：
4、共轭梯度法
共轭梯度法是介于最速下降法与牛顿法之间的一个方法，它仅需一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点。
共轭梯度法的基本思想是把共轭性与最速下降法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜索，求出目标函数的极小点。根据共轭方向的基本性质，这种方法具有二次终止性。
共轭梯度法中的核心迭代过程可以采取不同的方案，一种是直接延续，即总是用d^(k+1)=-g(k+1)+beta_k*d^(k)构造搜索方向;一种是把n步作为一轮，每搜索一轮之后，取一次最速下降方向，开始下一轮，此种策略称为“重置”。
下面我们介绍一种传统的共轭梯度法
注意，上述算法中均假设采用的精确一维搜索，但实际计算中，精确一维搜索会带来一定困难，代价较大，通常还是采用不精确的一维搜索。但此时(4)中构造的搜索方向可能就不是下降方向了，解决这个问题有两个方法。
其一，当d^(k+1)不是下降方向时，以最速下降方向重新开始。事实上，这也存在问题，但一维搜索比较粗糙时，这样重新开始可能是大量的，会降低计算效率。
其二，在计算过程中增加附加的检验，具体细节可以参考陈宝林老师的“最优化理论与方法”的P301。