机器学习中最优化问题总结

记录一下自己学习机器学习中最优化问题的过程。不断更新中。

一、基础知识

1.1 梯度向量、Jacobian矩阵和Hessian矩阵

这里讨论的三个概念：梯度向量、Jacobian矩阵和Hessian矩阵

它的自变量：$\mathbf{x}=(x_1,x_2,\cdots ,x_n{)}^T$

因变量有两种情况：

• 一维$f(\mathbf{x})$：
  • 此时的一阶导数构成的向量为梯度向量$g(\mathbf{x})$
  • 二阶导数构成的矩阵为Hessian矩阵
• 多维$\mathbf{f}(\mathbf{x})={\left(f_1(\mathbf{x}),f_2(\mathbf{x}),\cdots ,f_m(\mathbf{x})\right)}^T$
  • 此时的一阶导数构成的矩阵为Jacobian矩阵

1.1.1 梯度向量

即目标函数$f$为单变量，它是关于自变量向量$\mathbf{x}=(x_1,x_2,\cdots ,x_n{)}^T$的函数，此时，单变量函数$f$对向量$\mathbf{x}$求梯度，得到的结果为一个与向量$\mathbf{x}$同维度的向量，称之为梯度向量
$$g(\mathbf{x})=▽f(\mathbf{x})={\left(\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2},\cdots ,\frac{\partial f}{\partial x_n}\right)}^T$$

1.1.2 Jacobian矩阵

Jacobian矩阵中文称为雅可比矩阵。即目标函数$\mathbf{f}$为一个函数向量，它的每一个元素都是关于自变量向量$\mathbf{x}=(x_1,x_2,\cdots ,x_n{)}^T$的函数，此时，函数向量$\mathbf{f}$对向量$\mathbf{x}$求梯度，得到的结果为一个矩阵，行数为函数向量$\mathbf{f}$的维度，列数为自变量$\mathbf{x}$的维度，称之为Jacobian矩阵
$$▽\mathbf{f}(\mathbf{x})={\left[\begin{array}{cccc}\frac{\partial f_1(\mathbf{x})}{\partial x_1}& \frac{\partial f_1(\mathbf{x})}{\partial x_2}& \cdots & \frac{\partial f_1(\mathbf{x})}{\partial x_n}\\ \frac{\partial f_2(\mathbf{x})}{\partial x_1}& \frac{\partial f_2(\mathbf{x})}{\partial x_2}& \cdots & \frac{\partial f_2(\mathbf{x})}{\partial x_n}\\ \cdots & \cdots & \cdots & \cdots \\ \frac{\partial f_m(\mathbf{x})}{\partial x_1}& \frac{\partial f_m(\mathbf{x})}{\partial x_2}& \cdots & \frac{\partial f_m(\mathbf{x})}{\partial x_n}\end{array}\right]}_{m\times n}=\left[\begin{array}{c}{g}_1(\mathbf{x}{)}^T\\ g_2(\mathbf{x}{)}^T\\ \cdots \\ g_m(\mathbf{x}{)}^T\end{array}\right]$$

它的每一行都是由相应的函数的梯度向量的转置构成的。

梯度向量和Jacobi矩阵的关系：实际上，梯度向量是Jacobi矩阵的一个特例！当目标函数为标量函数时，Jacobi矩阵就是梯度向量。

1.1.3 Hessian矩阵

Hessian矩阵中文称为海赛矩阵。实际上，Hessian矩阵就是梯度向量$g(\mathbf{x})$对自变量$\mathbf{x}$的Jacobian矩阵，是一个方阵：
$$\begin{gathered} G(\mathbf{x})={▽}_{\mathbf{x}}[g(\mathbf{x})]={▽}_{\mathbf{x}}\left[{\left(\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2},\cdots ,\frac{\partial f}{\partial x_n}\right)}^T\right] \\ =\left[\begin{array}{cccc}\frac{\partial f}{\partial x_1{x}_1}& \frac{\partial f}{\partial x_1{x}_2}& \cdots & \frac{\partial f}{\partial x_1{x}_n}\\ \frac{\partial f}{\partial x_2{x}_1}& \frac{\partial f}{\partial x_2{x}_2}& \cdots & \frac{\partial f}{\partial x_2{x}_n}\\ \cdots & \cdots & \cdots & \cdots \\ \frac{\partial f}{\partial x_n{x}_1}& \frac{\partial f}{\partial x_n{x}_2}& \cdots & \frac{\partial f}{\partial x_n{x}_n}\end{array}\right] \end{gathered}$$

二、无约束优化

2.1 梯度下降法

梯度下降法(Gradient descent)或最速下降法(Steepest descent)是求解无约束最优化问题的一种最常见的方法，有实现简单的优点。梯度下降法是迭代算法，每一步需要求解目标函数的梯度向量。

2.1.1 梯度下降与梯度上升

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

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

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

2.1.2 梯度下降的直观解释

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

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

2.1.3 梯度下降的相关概念

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

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

2. 特征（feature）：指的是样本中输入部分，比如2个单特征的样本$(x^{(0)},y^{(0)}),(x^{(1)},y^{(1)})$，则第一个样本特征为$x^{(0)}$，第一个样本输出为$y^{(0)}$。

3. 假设函数（hypothesis function）：在监督学习中，为了拟合输入样本，而使用的假设函数，记为$h_{\theta }(x)$。比如对于单个特征的$m$个样本$(x^{(i)},y^{(i)})(i=1,2,...m)$，可以采用拟合函数如下：$h_{\theta }(x)={\theta }_0+{\theta }_{1}x$。

4. 损失函数（loss function）：为了评估模型拟合的好坏，通常用损失函数来度量拟合的程度。损失函数极小化，意味着拟合程度最好，对应的模型参数即为最优参数。在线性回归中，损失函数通常为样本输出和假设函数的差取平方。比如对于$m$个样本$（x_i,y_i）(i=1,2,...m)$，采用线性回归，损失函数为：
   $$J({\theta }_0,{\theta }_1)=\sum _{i=1}^m(h_{\theta }(x_i)-y_i{)}^2$$
   其中$x_i$表示第$i$个样本特征，$y_i$表示第$i$个样本对应的输出，$h_{\theta }(x_i)$为假设函数。

2.1.4 梯度下降的详细算法

梯度下降法的算法可以有代数法和矩阵法（也称向量法）两种表示。这里先介绍代数法，后介绍矩阵法。

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

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

比如对于线性回归，假设函数表示为$h_{\theta }(x_1,x_2,...x_n)={\theta }_0+{\theta }_1{x}_1+...+{\theta }_n{x}_n$，其中${\theta }_i(i=0,1,...,n)$为模型参数，$x_i(i=0,1,...,n)$为每个样本的$n$个特征值。这个表示可以简化，我们增加一个特征$x_0=1$ ，即：$h_{\theta }(x_0,x_1,...x_n)=\sum _{i=0}^n{\theta }_i{x}_i$

线性回归的损失函数为：
$$J({\theta }_0,{\theta }_1...,{\theta }_n)=\frac{1}{2m}\sum _{j=0}^m(h_{\theta }(x_0^{(j)},x_1^{(j)},...x_n^{(j)})-y_j{)}^2$$

2、 算法相关参数初始化：主要是初始化${\theta }_0,{\theta }_1...,{\theta }_n$，算法终止距离$\epsilon$以及步长$\alpha$。$\epsilon$可以初始化为0， 步长可以初始化为1。在调优的时候再优化。

3、算法过程：
1）确定当前位置的损失函数的梯度，对于${\theta }_i$，其梯度表达式如下：
$$\frac{\partial }{\partial {\theta }_i}J({\theta }_0,{\theta }_1...,{\theta }_n)$$

2）用步长乘以损失函数的梯度，得到当前位置下降的距离，即$\alpha \frac{\partial }{\partial {\theta }_i}J({\theta }_0,{\theta }_1...,{\theta }_n)$对应于前面登山例子中的某一步。

3）确定是否所有的${\theta }_i$，梯度下降的距离都小于$\epsilon$，如果小于ε则算法终止，当前所有的${\theta }_i$即为最终结果。否则进入步骤4.

4）更新所有的$\theta$，对于${\theta }_i$，其更新表达式如下。更新完毕后继续转入步骤1.
$${\theta }_i={\theta }_i-\alpha \frac{\partial }{\partial {\theta }_i}J({\theta }_0,{\theta }_1...,{\theta }_n)$$

下面用线性回归的例子来具体描述梯度下降。假设我们的样本是$(x_1^{(0)},x_2^{(0)},...x_n^{(0)},y_0),(x_1^{(1)},x_2^{(1)},...x_n^{(1)},y_1),...(x_1^{(m)},x_2^{(m)},...x_n^{(m)},y_m)$，损失函数如前面先决条件所述：$J({\theta }_0,{\theta }_1...,{\theta }_n)=\frac{1}{2m}\sum _{j=0}^m(h_{\theta }(x_0^{(j)},x_1^{(j)},...x_n^{(j)})-y_j{)}^2$

则在算法过程步骤1中对于${\theta }_i$的偏导数计算如下： 　　
$$\frac{\partial }{\partial {\theta }_i}J({\theta }_0,{\theta }_1...,{\theta }_n)=\frac{1}{m}\sum _{j=0}^m(h_{\theta }(x_0^{(j)},x_1^{(j)},...x_n^{(j)})-y_j)x_i^{(j)}$$

由于样本中没有$x_0$上式中令所有的$x_0^j$为1.

步骤4中${\theta }_i$的更新表达式如下：
$${\theta }_i={\theta }_i-\alpha \frac{1}{m}\sum _{j=0}^m(h_{\theta }(x_0^{(j)},x_1^{(j)},...x_n^j)-y_j)x_i^{(j)}$$

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

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

1、 先决条件： 和上一节类似， 需要确认优化模型的假设函数和损失函数。对于线性回归，假设函数$h_{\theta }(x_1,x_2,...x_n)={\theta }_0+{\theta }_1{x}_1+...+{\theta }_n{x}_n$的矩阵表达方式为：$h_{\theta }(\mathbf{X})=\mathbf{X}\theta$，其中， 假设函数$h_{\theta }(\mathbf{X})$为$m\ast 1$的向量，$\theta$为$(n+1)\ast 1$的向量，里面有$n+1$个代数法的模型参数。$\mathbf{X}$为$m\ast (n+1)$维的矩阵。$m$代表样本的个数，$n+1$代表样本的特征数。

损失函数的表达式为：$J(\theta )=\frac{1}{2}(\mathbf{X}\theta -\mathbf{Y}{)}^T(\mathbf{X}\theta -\mathbf{Y})$, 其中$\mathbf{Y}$是样本的输出向量，维度为$m\ast 1$.

2、算法相关参数初始化：与上一节一样。

3、算法过程：

1）确定当前位置的损失函数的梯度，对于$\theta$向量，其梯度表达式如下：$\frac{\partial }{\partial \theta }J(\theta )$

2）用步长乘以损失函数的梯度，得到当前位置下降的距离，即$\alpha \frac{\partial }{\partial \theta }J(\theta )$对应于前面登山例子中的某一步。

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

4）更新$\theta$向量，其更新表达式如下。更新完毕后继续转入步骤1。
$$\theta =\theta -\alpha \frac{\partial }{\partial \theta }J(\theta )$$

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

损失函数对于$\theta$向量的偏导数计算如下：
$$\frac{\partial }{\partial \theta }J(\theta )={\mathbf{X}}^T(\mathbf{X}\theta -\mathbf{Y})$$

步骤4中$\theta$向量的更新表达式如下:
$$\theta =\theta -\alpha {\mathbf{X}}^T(\mathbf{X}\theta -\mathbf{Y})$$

相比于代数法，可以看到矩阵法要简洁很多。

注意，$J(\theta )$求导用到了矩阵求导链式法则，和两个矩阵求导的公式。

公式1：$\frac{\partial }{\partial \mathbf{x}}({\mathbf{x}}^{\mathbf{T}}\mathbf{x})=2\mathbf{x}$，$\mathbf{x}$为向量

公式2：${\nabla }_{X}f(AX+B)=A^T{\nabla }_{Y}f,Y=AX+B,f(Y)$为标量

2.1.5 梯度下降的算法调优

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

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

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

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

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

2.1.6 BGD，SGD，MBGD

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

批量梯度下降法，是梯度下降法最常用的形式，具体做法也就是在更新参数时使用所有的样本来进行更新，我们在上一小节介绍的线性回归梯度下降算法是批量梯度下降法。　
$${\theta }_i={\theta }_i-\alpha \frac{1}{m}\sum _{j=0}^m(h_{\theta }(x_0^{(j)},x_1^{(j)},...x_n^{(j)})-y_j)x_i^{(j)}$$

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

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

随机梯度下降法，其实和批量梯度下降法原理类似，区别在与求梯度时没有用所有的$m$个样本的数据，而是仅仅选取一个样本$j$来求梯度。对应的更新公式是：
$${\theta }_i={\theta }_i-\alpha (h_{\theta }(x_0^{(j)},x_1^{(j)},...x_n^{(j)})-y_j)x_i^{(j)}$$

随机梯度下降是通过每个样本来迭代更新一次，对比上面的批量梯度下降，迭代一次需要用到所有训练样本（往往如今真实问题训练数据都是非常巨大），一次迭代不可能最优，如果迭代10次的话就需要遍历训练样本10次。但是，SGD伴随的一个问题是噪音较BGD要多，使得SGD并不是每次迭代都向着整体最优化方向。但是大体上是往着最优值方向移动。

③MBGD

小批量梯度下降法是批量梯度下降法和随机梯度下降法的折中，也就是对于$m$个样本，我们采用$x$个样本来迭代，$1<x<m$。一般可以取$x=10$，当然根据样本的数据，可以调整这个$x$的值。对应的更新公式是：
$${\theta }_i={\theta }_i-\alpha \frac{1}{x}\sum _{j=t}^{t+x-1}(h_{\theta }(x_0^{(j)},x_1^{(j)},...x_n^{(j)})-y_j)x_i^{(j)}$$

④三种梯度下降方法的总结

1、批梯度下降每次更新使用了所有的训练数据，最小化损失函数，如果只有一个极小值，那么批梯度下降是考虑了训练集所有数据，是朝着最小值迭代运动的，但是缺点是如果样本值很大的话，更新速度会很慢。

2、随机梯度下降在每次更新的时候，只考虑了一个样本点，这样会大大加快训练数据，也恰好是批梯度下降的缺点，但是有可能由于训练数据的噪声点较多，那么每一次利用噪声点进行更新的过程中，就不一定是朝着极小值方向更新，但是由于更新多轮，整体方向还是大致朝着极小值方向更新，又提高了速度。

3、小批量梯度下降法是为了解决批梯度下降法的训练速度慢，以及随机梯度下降法的准确性综合而来，但是这里注意，不同问题的batch是不一样的，只能通过实验结果来进行超参数的调整。

2.1.7 梯度下降法与最小二乘法的比较

我在这里没有写最小二乘法，实际上最小二乘法与我们上面介绍的很像，在梯度下降法的矩阵方式描述中，介绍了$\theta$的更新方式，我们稍微整理下就得到了最小二乘法求参数$\theta$的公式：
$$\theta =({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{Y}$$

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

2.2 牛顿法和拟牛顿法

牛顿法（Newton method）和拟牛顿法（quasi Newton method）也是求解无约束最优化问题的常用方法，有收敛速度快的优点。牛顿法是迭代算法，每一步都要求解目标函数的Hessian矩阵的逆矩阵，计算比较复杂。拟牛顿法通过正定矩阵近似Hessian矩阵的逆矩阵或者Hessian矩阵，简化了这一计算过程。

2.2.1 牛顿法

一般来说, 牛顿法主要应用在两个方面。1, 求方程的根; 2, 非线性最优化问题。

1) 求解方程

并不是所有的方程都有求根公式，或者求根公式很复杂，导致求解困难。 利用牛顿法，可以迭代求解。

原理是利用泰勒公式, 在$x_0$处展开, 且展开到一阶, 即
$$f(x)=f(x_0)+(x–x_0)f^{\prime }(x_0)\text{(1)}$$

求解方程$f(x)=0$, 即$f(x_0)+(x–x_0)f^{\prime }(x_0)=0$, 求解$x=x_1=x_0–f(x_0)/f^{\prime }(x_0)$, 因为这是利用泰勒公式的一阶展开, $f(x)=f(x_0)+(x–x_0)f^{\prime }(x_0)$处并不是完全相等, 而是近似相等, 这里求得的$x_1$并不能让$f(x)=0$, 只能说$f(x_1)$的值比$f(x_0)$更接近$f(x)=0$, 于是乎, 迭代求解的想法就很自然了, 可以进而推出$x_{n+1}=x_n–f(x_n)/f^{\prime }(x_n)$, 通过迭代, 这个式子必然在$f(x^{\ast })=0$的时候收敛. 整个过程如下图：

上图中的Funktion是德语，函数的意思；Tangente是切线的意思。

2) 非线性最优化问题

在非线性优化最优化的问题中, 牛顿法提供了一种求解的办法. 假设任务是优化一个目标函数$f$, 求函数$f$的极大极小问题, 可以转化为求解函数$f$的导数$f^{\prime }=0$的问题, 这样求可以把优化问题看成方程求解问题($f^{\prime }=0$). 剩下的问题就和第一部分提到的牛顿法求解很相似了。

为了求解$f^{\prime }=0$的根, 首先把$f(x)$在探索点$x_n$处二阶泰勒展开：

$$f(x)=f(x_n)+f^{\prime }(x_n)(x–x_n)+\frac{f^{\prime \prime }(x_n)}{2}(x–x_n{)}^2$$

类比公式（1），即有:
$$f^{\prime }(x)=f^{\prime }(x_n)+f^{\prime \prime }(x_n)(x–x_n)=0$$

求得出迭代公式：
$$x=x_{n+1}=x_n–\frac{f^{\prime }(x_n)}{f^{\prime \prime }(x_n)},n=0,1,\dots \text{(2)}$$

一般认为牛顿法可以利用到曲线本身的信息, 比梯度下降法更容易收敛（迭代更少次数）, 如下图是一个最小化一个目标方程的例子, 红色曲线是利用牛顿法迭代求解, 绿色曲线是利用梯度下降法求解.

上图讨论的是2维情况, 高维情况的牛顿迭代公式是：
$$x_{n+1}=x_n–H_n^{–1}{g}_n,n\ge 0\text{(3)}$$

其中，$g_n=g(x_n)=\nabla f(x_n)$是$f(x)$的梯度向量在点$x_n$的值。$H_n=H(x_n)$是$f(x)$的海赛矩阵
$$H(x)=[\frac{{\partial }^{2}f}{\partial x_i\partial x_j}{]}_{n\times n}$$

在点$x_n$的值。

或者，把公式（3）重写为：

$$x_{n+1}=x_n+p_n\text{(4)}$$

其中，
$$H_n{p}_n=-g_n\text{(5)}$$

具体算法如下：

输入：目标函数 $f(x)$ ，梯度 $g(x)=\nabla f(x)$ ，海塞矩阵 $H(x)$ ，精度要求 $\epsilon$

输出： $f(x)$ 的极小点 $x_{\ast }$

(1) 取初始值 $x^{(0)}$ ，置 $n=0$
(2) 计算 $g_n=g(x_n)$
(3) 若 $||g_n||<\epsilon$ ，则停止计算，得近似解 $x_{\ast }=x_n$
(4) 计算 $H_n=H(x_n)$ ，并求 $p_n$ ： $H_n{p}_n=-g_n$
(5) 置 $x_{n+1}=x_n+p_n$
(6) 置 $n=n+1$ ，转(2)

2.2.2 拟牛顿法

(一) 拟牛顿法的思路

在牛顿法的迭代中，需要计算海赛矩阵的逆矩阵，这一计算比较复杂，考虑用一个 $n$ 阶矩阵 $G_n=G(x_n)$ 来近似替代 $H_n^{-1}=H^{-1}(x_n)$ 。这就是拟牛顿法的基本想法。

注意拟牛顿法的 关键词：

• 不用算二阶导数
• 不用求逆

先看牛顿法迭代中海赛矩阵 $H_n$ 满足的条件。首先， $H_n$ 满足以下关系。根据泰勒展开，有 $\nabla f(x)=g_n+H_n(x-x_n)$，在其中取 $x=x_{n+1}$ ，得
$$g_{n+1}-g_n=H_n(x_{n+1}-x_n)$$

记 $y_n=g_{n+1}-g_n$ ， ${\delta }_n=x_{n+1}-x_n$ ，则
$$y_n=H_n{\delta }_n\text{(6)}$$

或
$$H_n^{-1}{y}_n={\delta }_n\text{(7)}$$

公式(6)和(7)即称为拟牛顿条件。如果 $H_n$ 是正定的( $H_n^{-1}$也是正定的)，那么可以保证牛顿法搜索方向 $p_n$ 是下降方向。这是因为搜索方向是 $p_n=-H_n^{-1}{g}_n$ ，结合牛顿法迭代式 $x_{n+1}=x_n-H_n^{-1}{g}_n$ 有
$$x=x_n+\lambda p_n=x_n-\lambda H_n^{-1}{g}_n$$

结合上式，有$x-x_n=-\lambda H_n^{-1}{g}_n$，所以 $f(x)$ 在 $x_n$ 的泰勒展开式可以近似写成
$$f(x)=f(x_n)-\lambda g_n^T{H}_n^{-1}{g}_n$$

因 $H_n^{-1}$ 正定，故有 $g_n^T{H}_n^{-1}{g}_n>0$ 当 $\lambda$ 为一个充分小的正数时，总有 $f(x)<f(x_n)$ ，也就是说 $p_n$ 是下降方向。

拟牛顿法将 $G_n$ 作为 $H_n^{-1}$ 的近似，要求矩阵 $G_n$ 满足同样的条件。首先，每次迭代矩阵 $G_n$ 是正定的。同时， $G_n$ 满足下面的拟牛顿条件：
$$G_{n+1}{y}_n={\delta }_n\text{(8)}$$

注意公式（8）中$G$的下标为$n+1$。

按照拟牛顿条件，选择 $G_n$ 作为 $H_n^{-1}$ 的近似或选择 $B_n$ 作为 $H_n$ 的近似的算法称为拟牛顿法（注：$B_n$可以看作是 $G_n$的逆，分别对应着接下来的BFGS算法和DFP算法）。

按照拟牛顿条件，在每次迭代中可以选择更新矩阵 $G_{n+1}$ ：
$$G_{n+1}=G_n+\Delta G_n$$

这种选择有一定的灵活性，因此有多种实现方法。下面DFP算法和BFGS算法。

(二) DFP(Davidon-Fletcher-Powell)算法

DFP算法选择 $G_{n+1}$ 的方法是，假设每一步迭代中矩阵 $G_{n+1}$ 是由 $G_n$ 加上两个附加项构成的，即
$$G_{n+1}=G_n+P_n+Q_n$$

其中 $P_n,Q_n$ 是待定矩阵。这时，
$$G_{n+1}{y}_n=G_n{y}_n+P_n{y}_n+Q_n{y}_n$$

为使 $G_{n+1}$ 满足拟牛顿条件，可使 $P_n$ 和 $Q_n$ 满足：
$$P_n{y}_n={\delta }_n$$

$$Q_n{y}_n=-G_n{y}_n$$

这样即可满足公式（8）。

事实上，不难找出这样的 $P_n$ 和 $Q_n$ ，例如取
$$P_n=\frac{{\delta }_n{\delta }_n^T}{{\delta }_n^T{y}_n}$$

$$Q_n=-\frac{G_n{y}_n{y}_n^T{G}_n}{y_n^T{G}_n{y}_n}$$

这样就可得到矩阵 $G_{n+1}$ 的迭代公式：
$$G_{n+1}=G_n+\frac{{\delta }_n{\delta }_n^T}{{\delta }_n^T{y}_n}-\frac{G_n{y}_n{y}_n^T{G}_n}{y_n^T{G}_n{y}_n}\text{(9)}$$

称为DFP算法。

可以证明，如果初始矩阵 $G_0$ 是正定的，则迭代过程中的每个矩阵 $G_n$ 都是正定的。

DFP具体算法如下：

输入：目标函数 $f(x)$ ，梯度 $g(x)=\nabla f(x)$ ，精度要求 $\epsilon$ ；

输出： $f(x)$ 的极小点 $x_{\ast }$

(1) 选定初始点 $x_0$ ，取 $G_0$ 为正定对称矩阵，置 $n=0$
(2) 计算 $g_n=g(x_n)$。 若 $||g_n||<\epsilon$ ，则停止计算，得近似解 $x_{\ast }=x_n$
(3) 置 $p_n=-G_n{g}_n$（注：参考公式(5)）
(4) 一维搜索：求 ${\lambda }_n$ 使得 $$f(x_n+{\lambda }_n{p}_n)=\underset{\lambda \ge 0}{\min }f(x_n+\lambda p_n)$$

（注：找到${\lambda }_n$使得$f(x_n+\lambda p_n)$最小）
(5) 置 $x_{n+1}=x_n+{\lambda }_n{p}_n$
(6) 计算 $g_{n+1}=g(x_{n+1})$。若 $||g_{n+1}||<\epsilon$ ，则停止计算，得近似解 $x_{\ast }=x_{n+1}$ ；否则，根据公式（9）计算出 $G_{n+1}$
(7) 置 $n=n+1$ ，转(3)

(三) BFGS(Broyden-Fletcher-Goldfard-Shanno)算法

BFGS算法是最流行的拟牛顿算法。可以考虑用 $G_n$ 逼近海赛矩阵的逆矩阵 $H^{-1}$ ，也可以考虑用 $B_n$ 逼近海赛矩阵 $H$ 。这时，相应的拟牛顿条件是
$$B_{n+1}{\delta }_n=y_n\text{(10)}$$

可以用同样的方法得到另一迭代公式。首先令
$$B_{n+1}=B_n+P_n+Q_n$$

$$B_{n+1}{\delta }_n=B_n{\delta }_n+P_n{\delta }_n+Q_n{\delta }_n$$

考虑使 $P_n$ 和 $Q_n$ 满足：
$$P_n{\delta }_n=y_n$$

$$Q_n{\delta }_n=-B_n{\delta }_n$$

找出适合条件的 $P_n$ 和 $Q_n$，得到BFGS算法矩阵 $B_{n+1}$ 的迭代公式：
$$B_{n+1}=B_n+\frac{y_n{y}_n^T}{y_n^T{\delta }_n}-\frac{B_n{\delta }_n{\delta }_n^T{B}_n}{{\delta }_n^T{B}_n{\delta }_n}\text{(11)}$$

我们可以看到，BFGS与DFP是对偶的，只是$G\leftrightarrow B$，$y\leftrightarrow \delta$互调。

可以证明，如果初始矩阵 $B_0$ 是正定的，则迭代过程中的每个矩阵 $B_n$ 都是正定的。

但是，我们需要注意的是：我们在计算搜索方向$p_n$的时候，还需要计算$B_n^{-1}$，即由 $B_n{p}_n=-g_n$ 求出 $p_n$。在这里，不也一样要求逆吗？这就要引入谢尔曼莫里森公式了。

---

Sherman-Morrison 公式

对于任意非奇异方阵$A$，$u,v\in R^n$是$n$维向量，若$1+v^T{A}^{-1}u\ne 0$，则
$$(A+u{v}^T{)}^{-1}=A^{-1}-\frac{(A^{-1}u)(v^T{A}^{-1})}{1+v^T{A}^{-1}u}$$

该公式描述了在矩阵$A$发生某种变化时，如何利用之前求好的逆，求新的逆。

对迭代公式(11)引入两次 Sherman-Morrison 公式就能得到
$$B_{n+1}^{-1}=\left(I_n-\frac{\Delta {\delta }_n\Delta y_n^T}{\Delta {\delta }_n^T\Delta g_n}\right)B_n^{-1}\left(I_n-\frac{\Delta y_n\Delta {\delta }_n^T}{\Delta {\delta }_n^T\Delta y_n}\right)+\frac{\Delta {\delta }_n\Delta {\delta }_n^T}{\Delta {\delta }_n^T\Delta y_n}\text{(12)}$$

(注：$I_n$中的$n$是维度的含义，不是迭代的次数)

就得到了逆矩阵之间的推导。可能有人会问，第一个矩阵不也要求逆吗？其实这是一个迭代算法，初始矩阵设为单位矩阵（对角阵也可以）就不用求逆了。

这个公式的详细推导可以参考这里。

---

BFGS具体算法如下：

输入：目标函数 $f(x)$ ，梯度 $g(x)=\nabla f(x)$ ，精度要求 $\epsilon$ ；

输出： $f(x)$ 的极小点 $x_{\ast }$

(1) 选定初始点 $x_0$ ，取 $B_0^{-1}$ 为正定对称矩阵，置 $n=0$
(2) 计算 $g_n=g(x_n)$ 若 $||g_n||<\epsilon$ ，则停止计算，得近似解 $x_{\ast }=x_n$
(3) 由 $B_n{p}_n=-g_n$ 求出 迭代方向$p_n=-B_n^{-1}{g}_n$
(4) 一维搜索：求 ${\lambda }_n$ 使得
$$f(x_n+{\lambda }_n{p}_n)=\underset{\lambda \ge 0}{\min }f(x_n+\lambda p_n)$$

(5) 置 $x_n=x_n+{\lambda }_n{p}_n$
(6) 计算 $g_{n+1}=g(x_{n+1})$。若 $||g_{n+1}||<\epsilon$ ，则停止计算，得近似解 $x_{\ast }=x_{n+1}$；否则，按公式(12)计算出 $B_{n+1}$
(7) 置 $n=n+1$ ，转(3)

该算法的过程与DFP算法的流程几乎一样。

(四)L-BFGS（Limited-memory BFGS）算法

对于$n$维参数，BFGS算法需要保存一个$O(d^2)$大小的$B^{-1}$矩阵，实际上只需要每一轮的$\Delta \delta$和$\Delta y$，也可以递归计算出当前迭代的$B^{-1}$矩阵，L-BFGS就是基于这种思想，实现了节省内存的BFGS。

L-BFGS推导

BFGS的递推公式(12)如下：
$$B_{n+1}^{-1}=\left(I_n-\frac{\Delta {\delta }_n\Delta y_n^T}{\Delta {\delta }_n^T\Delta g_n}\right)B_n^{-1}\left(I_n-\frac{\Delta y_n\Delta {\delta }_n^T}{\Delta {\delta }_n^T\Delta y_n}\right)+\frac{\Delta {\delta }_n\Delta {\delta }_n^T}{\Delta {\delta }_n^T\Delta y_n}\text{(12)}$$

现在假设${\rho }_n=\frac{1}{\Delta {\delta }_n^T\Delta y_n}$，$V_n=I_n-{\rho }_n\Delta y_n\Delta {\delta }_n^T$，则递推公式可以写成
$$B_{n+1}^{-1}=V_n^T{B}_n^{-1}{V}_n+{\rho }_n\Delta {\delta }_n\Delta {\delta }_n^T\text{(13)}$$

给定的初始矩阵$B_0^{-1}$后，之后的每一轮都可以递推计算
$$B_1^{-1}=V_0^T{B}_0^{-1}{V}_0+{\rho }_0\Delta {\delta }_0\Delta {\delta }_0^T$$

$$B_2^{-1}=V_1^T{B}_0^{-1}{V}_1+{\rho }_1\Delta {\delta }_1\Delta {\delta }_1^T=(V_1^T{V}_0^T)B_0^{-1}(V_0{V}_1)+V_1^T{\rho }_0\Delta {\delta }_0\Delta {\delta }_0^T{V}_1+{\rho }_1\Delta {\delta }_1\Delta {\delta }_1^T$$

一直到最后$B_{k+1}^{-1}$可以由$n=0$到$n=k$的$\Delta {\delta }_n$和$\Delta y_n$表示：
$$\begin{array}{cccc}{B}_{n+1}^{-1}& =& & (V_n^T{V}_{n-1}^T\cdots V_1^T{V}_0^T)B_0^{-1}(V_0\cdots V_{n-1}{V}_n)\\ & & +& (V_n^T{V}_{n-1}^T\cdots V_2^T{V}_1^T)({\rho }_0\Delta x_0\Delta x_0^T)(V_1\cdots V_{n-1}{V}_n)\\ & & +& \cdots \\ & & +& V_n^T({\rho }_{n-1}\Delta x_{n-1}\Delta x_{n-1}^T)V_n\\ & & +& {\rho }_n\Delta x_n\Delta x_n^T\end{array}$$

看起来很长，其实可以写成一个求和项
$$B_{n+1}^{-1}=\left(\prod _{i=n}^0{V}_i^T\right)B_0^{-1}\left(\prod _{i=0}^n{V}_i\right)+\sum _{j=0}^n\left(\prod _{i=n}^{j+1}{V}_i^T\right)\left({\rho }_j\Delta x_j\Delta x_j^T\right)\left(\prod _{i=j+1}^n{V}_i\right)$$

这个求和项包含了从0到$n$的所有$\Delta \delta$和$\Delta y$，而根据实际需要，可以只取最近的$m$个，也就是：
$$B_n^{-1}=\left(\prod _{i=n-1}^{n-m}{V}_i^T\right)B_0^{-1}\left(\prod _{i=n-m}^{n-1}{V}_i\right)+\sum _{j=n-1}^{n-m}\left(\prod _{i=n}^{j+1}{V}_i^T\right)\left({\rho }_j\Delta x_j\Delta x_j^T\right)\left(\prod _{i=j+1}^n{V}_i\right)$$

工程上的L-BFGS

我们关心的其实不是$B_t^{-1}$本身如何，算$B_t^{-1}$的根本目的是要算本轮搜索方向$p_n=-B_n^{-1}{g}_n$

以下算法摘自《Numerical Optimization》，它可以高效地计算出拟牛顿法每一轮的搜索方向。仔细观察一下，你会发现它实际上就是复现上面推导的那一堆很长的递推公式，你所需要的是最近$m$轮的$\Delta \delta$和$\Delta y$，后向和前向算完得到最终的 $r$就是搜索方向$B_n^{-1}{g}_n$，之后要做一维搜索或者什么的都可以。

解释一下算法的符号和本文符号之间的对应关系，$s_i=\Delta {\delta }_i$，$y_i=\Delta y_i$，$H_k=B_k^{-1}$

L-BFGS具体算法如下：

（1）给定初始点 $x_0$，允许误差$\epsilon$，预定保留最近$m$个向量，设置 $B_0^{-1}$，$n=0$
（2）用Algorithm 9.1计算搜索方向 $p_n=-B_n^{-1}{g}_n$
（3）从点 $x_n$ 出发，沿着 $p_n$ 做一维搜索，获得最优步长并更新参数：
$$\begin{gathered} f(x_n+{\lambda }_n{p}_n)=\underset{\lambda \ge 0}{\min }f(x_n+\lambda p_n) \\ {x}_{n+1}=x_n+{\lambda }_n\cdot p_n \end{gathered}$$

（4）判断精度，若$||g_{n+1}||<\epsilon$则停止迭代，否则转（5）
（5）判断 $n>m$，删掉存储的$\Delta {\delta }_{n-m}$和 $\Delta y_{n-m}$
（5）计算$\Delta \delta ={\delta }_{n+1}-{\delta }_n$，$\Delta y=y_{n+1}-y_n$，令$n=n+1$，转（2）

BFGS与DFP算法的比较

通常BFGS效果要好于DFP，为什么呢？

BFGS优于DFP的原因在于，BFGS有自校正的性质(self-correcting property)。通俗来说，如果某一步BFGS对Hessian阵的估计偏了，导致优化变慢，那么BFGS会在较少的数轮迭代内（取决于线搜索的质量），校正估计的Hessian阵。

参考自这里。

2.3 坐标下降法

坐标下降法是一种非梯度的优化方法，和梯度下降方法不同，它每次沿着单个维度方向进行搜索，交替搜索，当得到一个当前维度最小值之后再循环使用不同的维度方向，最终收敛得到最优解。

对于函数$f(x)$ ， $x=(x_1,x_2,...,x_n)$ ，

坐标下降方法可以表示为：其中$k$是迭代次数。

可以知道，这一序列与最速下降具有类似的收敛性质。如果在某次迭代中，函数得不到优化，说明一个驻点已经达到。

这一过程可以用下图表示。

但是，对于非平滑函数，坐标下降法可能会遇到问题。下图展示了当函数等高线非平滑，即不可微时，算法可能在非驻点中断执行。

如上图所示，最后会在红线交汇处就停止迭代了。这是因为，在红线交汇处，无论往哪个轴走，红线交汇处的值都是最小的，算法自然会在这里停止。

总结一下：坐标下降法与梯度下降法的不同之处在于不需要计算目标函数的梯度，每次迭代只用在单一维度上进行线性搜索，但是坐标下降法只适用于光滑函数，如果是非光滑函数可能会陷入非驻点从而无法更新。

参考文献

【1】Jacobian矩阵，Hessian矩阵和牛顿法
【2】机器学习中的最优化算法总结
【3】无约束优化
【4】最优化理论·光滑函数·Hessian矩阵·Jacobian矩阵·方向导数
【5】梯度下降法、牛顿法和拟牛顿法
【6】拟牛顿法（DFP、BFGS、L-BFGS）（这篇文章写的很好）
【7】梯度下降（Gradient Descent）小结