最小化

下降法

我们考虑多变量函数$f$∶$R^N→$以及为求最小值的向量$\bar{x}$的问题。即使是光滑函数，我们也会立即面对存在局部极小值和全局极小值的事实;现在我们满足于找到一个局部最小值。

我们不是在一次大型计算中寻找最小值，而是使用迭代策略，从一个点$\bar{x}$开始，将其更新为$\bar{x}+\bar{h}$，并重复这个过程。更新向量$\bar{h}$是通过直线搜索策略找到的：我们确定搜索方向$\bar{h}$，沿着该方向找到步长$\tau$，然后更新 $$ \bar{x}\leftarrow \bar{x}+\tau \bar{h}. $$

最陡下降

利用泰勒展开公式的多维形式（第16节），我们可以将任何足够光滑的函数写成 $$ f(\bar{x}+\bar{h})=f(\bar{x})+\bar{h}^{t} \nabla f+\cdots $$ 不难看出，选择$\bar{h}=−\nabla f$会得到最大的最小化，所以我们设置 $$ x_{\text {new }} \equiv x-\tau \nabla f $$ 这种方法被称为梯度下降或最陡下降。

对于新函数的值 $$ f(\bar{x}-\tau \nabla f) \approx f(\bar{x})-\tau|\nabla f|^{2} $$ 因此，对于$\tau$足够小的情况，这使得新函数的值既为正，又小于$f(\bar{x})$。

对于二次函数$f$，步长$\tau$可以计算出来，否则近似为: $$ f(\bar{x}+\tau \bar{h})=f(\bar{x})+\tau \bar{h}^{t} \nabla f+\frac{\tau^{2}}{2} h^{t}(\nabla \cdot \nabla f) h+\cdots $$ 忽略高阶项并设置$\delta f/\delta \tau= 0$给出 $$ \tau=-\bar{h}^{t} \nabla f / h^{t}(\nabla \cdot \nabla f) h $$ 另一种寻找合适步长的策略叫做回溯（backtracking）：

  • 从默认的步长开始，例如$\tau \equiv 1$。

  • 然后直到$f(\bar{x}+\tau \bar{h})<f(\bar{x})$，执行 $$ \tau \leftarrow \tau/2 $$

随机梯度下降法

虽然梯度在任何给定的步骤中都是最优的搜索方向，但总的来说它不一定是最好的。例如，二次问题，其中包括许多问题的偏微分方程，一个更好的选择是正交搜索方向。另一方面，在ML应用中，随机梯度下降(SGD)是一个很好的选择，其中坐标方向被用作搜索方向。在这种情况下 $$ f\left(x+\tau e{i}\right)=f(x)+\tau \frac{d f}{d e{i}}+\frac{\tau^{2}}{2} \frac{\delta^{2} f}{\delta e{i}^{2}} $$ 那么 $$ \tau=-(\nabla f){i} / \frac{\delta^{2} f}{\delta e_{i}^{2}} $$

代码

预知识

我们声明了一个vector类，它是一个标准vector类，并在其上定义了诸如加法和旋转等操作。有一个派生类unit_vector做了显而易见的工作。

框架

我们首先将函数定义为纯虚类，这意味着任何函数都需要支持这里提到的方法：

  1. // minimlib.h
  2. class function { 
  3.   public:
  4.   virtual double eval( const valuevector& coordinate ) const = 0;
  5.   virtual valuevector grad( const valuevector& coordinate ) const = 0;
  6.   virtual std::shared_ptr<matrix> delta( const valuevector& coordinate ) const = 0; 
  7.   virtual int dimension() const = 0;
  8. };

使用这样的函数可以定义各种更新步骤。例如，最陡的下降阶梯使用梯度：

  1. // minimlib.cxx
  2. valuevector steepest_descent_step
  3.   ( const function& objective,const valuevector& point ) {
  4.   auto grad = objective.grad(point); auto delta = objective.delta(point);
  5.   auto search_direction( grad );
  6.   auto hf = grad.inprod(search_direction);
  7.   auto hfh = search_direction.inprod( delta->mvp(search_direction) ); 
  8.   auto tau = - hf / hfh;
  9.   auto new_point = point + ( search_direction * tau );
  10.   return new_point;
  11. };

另一方面，随机下降是基于单位向量:

  1. valuevector stochastic_descent_step
  2. ( const function& objective,const valuevector& point, int idim,double damp) {
  3.   int dim = objective.dimension(); 
  4.   auto grad = objective.grad(point); 
  5.   auto delta = objective.delta(point);
  6.   auto search_direction( unit_valuevector(dim,idim) ); 
  7.   auto gradfi = grad.at(idim);
  8.   auto hf = gradfi;
  9.   auto hfh = delta->d2fdxi2(idim);
  10.   auto tau = - hf / hfh;
  11.   auto new_point = point + ( search_direction * (tau * damp) );
  12.   return new_point; 
  13. };

健全性测试

在圆周上使用最陡的下降可以在一步内收敛：

  1. Code:
  2. // ellipse.cxx
  3. ellipse circle
  4. ( valuevector( {1.,1.} ),valuevector( {0.,0.} ) );
  5. valuevector search_point( { 1.,1. } ); auto value = circle.eval(search_point);
  6. auto new_point = steepest_descent_step(circle, search_point);
  7. auto new_value = circle.eval(new_point); cout << "Starting point "
  8. << "(" << search_point.at(0)
  9. << "," << search_point.at(1) << "); " << "with value: " << value << "\n"
  10. << " becomes "
  11. << "(" << new_point.at(0)
  12. << "," << new_point.at(1) << "); "
  13. << "with value: " << new_value << endl;

另一方面，椭圆上最陡的下降需要多次迭代：

  1. Code:
  2. ellipse circle( valuevector( {1.,.1} ),valuevector( {0.,0.} ) );
  3. valuevector search_point( { 1.,1. } );
  4. auto value = circle.eval(search_point);
  5. for (int step=0; step<5 and value>.0001; step++) {
  6. auto new_point = steepest_descent_step(circle, search_point);
  7. auto new_value = circle.eval(new_point); cout << "Starting point " << fixed
  8. << "(" << setprecision(4) << search_point.at(0)
  9. << "," << setprecision(4) << search_point.at(1) << "); "
  10. << "with value: " << value << "\n"
  11. << " becomes "
  12. << "(" << setprecision(4) << new_point.at(0)
  13. << "," << setprecision(4) << new_point.at(1) << ")
  14. ;"
  15. << "with value: " << new_value << endl;
  16. search_point = new_point; value = new_value; }

牛顿法

牛顿法（或牛顿-拉弗森法）是一个求函数零点的迭代过程，即求出$f(x)= 0$的值$x$。它需要知道函数的导数$f’$，可以从17.1等数字中证明。

另一个理由来自于最小化:如果函数$f$是两次可微的，我们可以这样写 $$ f(x+h)=f(x)+h^{t} \nabla f+\frac{1}{2} h^{t}\left(\nabla^{2} f\right) h $$ 最小值为 $$ x_{\text {new }}=x-\left(\nabla^{2} f\right)^{-1} \nabla f $$ 这相当于求出梯度的零点： $$ 0=\nabla f(x+h)=\nabla f(x)+h^{t}\left(\nabla^{2} f(x)\right) h $$ 练习 17.1 让$f(x_1,x_2)=(10x_1^2 +x_2^2)/2 + 5\log(1+e^{−x_1−x_2})$。梯度下降法和牛顿法的收敛速度有多快？为了了解它们的不同行为，可以根据函数的关卡曲线绘制若干次迭代。

这个练习告诉我们梯度下降是错误的：它总是垂直于当前的水平曲线。然而，最小值并不一定在那个方向上。

Newton

迭代： $$ x_{n+1}=x_n-f(x_n)/f’(x_n) $$ 没有证明，我们声明：

  • 如果迭代的起始点足够接近于0，并且函数在0处是可微的，那么牛顿法将收敛到0。

  • 对于许多函数，牛顿方法不会收敛，但可以通过引入阻尼或做不精确的更新来获得收敛： $$ x_{n+1}=x_n-\alpha f(x_N)/f’(x_n) $$ 其中$\alpha<1$。

练习 17.2 牛顿的方法有可能是循环的。假设这是一个长度为2的循环： $$ x_0\rightarrow x_1\rightarrow x_2 = x_0 $$ 如果写出这个循环的方程，就会得到$f$的微分方程。解决方案是什么？为什么牛顿法在这个函数中不收敛?

在多维情况下，即以函数$f∶R^N\rightarrow R$牛顿法成为一个迭代线性系统解： $$ \bar{x}{n+1}=\bar{x}{n}-F\left(\bar{x}{n}\right)^{-1} f\left(\bar{x}{n}\right) $$ 其中$F$是$f$的雅可比矩阵。

由于牛顿法是一个迭代过程，我们不需要对这个线性系统的解达到完全精度，所以不精确的解，例如通过迭代求解，往往是可行的。

不精确牛顿法

牛顿方法不收敛的原因有很多（事实上，对于合适的函数，从它的收敛点绘制区域可以得到很好的分形）。由于这个原因，在实践中使用的是不精确的牛顿法。这种不精确表现在两个方面:

  • 代替“最佳”步长，我们使用一个分数。这通常是通过“回溯”来选择的，并采用了至少观察到一些函数值下降的选择。
  • 我们不用雅可比矩阵的逆而是用这个算子的近似值。例如，我们可以计算雅可比矩阵（通常是一个昂贵的过程），然后将其用于多个牛顿步骤。

可以证明，这些技术的组合可以保证收敛[209]。