牛顿法

Wolfram 语言所提供的一个明显的优势是它可以符号化地计算导数. 这意味着,当您指定 Method->"Newton" 并且函数明确可微时,符号化求导数计算就会自动执行. 另一方面,如果函数的形式不是明确可微的话,Wolfram 语言将使用有限差分近似计算Hessian,并使用结构信息来最小化所需计算的次数. 或者,您可以指定一个 Wolfram 语言表达式,用以对 Hessian 的变量赋予数值值.     

这里我们加载一个包含一些功用函数的程序包.
In[1]:=
Click for copyable input
在这个例子中,FindMinimum 用符号式计算Hessian并在需要时把 用数值值取代.
In[2]:=
Click for copyable input
Out[2]=
这里定义了一个函数,目的是使该函数只有其变量为数值值时才有值.
In[3]:=
Click for copyable input

我们无法找到该函数的符号式导数,因为该函数被定义为只对其变量的数值值取值.     

这里显示了 FindMinimum 在不得不利用有限差分来计算梯度和Hessian时所采取的步骤.
In[4]:=
Click for copyable input
Out[4]=

当梯度和Hessian都使用有限差分计算的时候,Hessian的误差可能非常大,因而使用不同的方法可能更好. 在里,FindMinimum 很准确地找到极小值,但是因导数信息不准确而不能确保. 此外,函数和梯度计算的次数远大于在自动计算符号式导数的例子中的次数,因为对梯度和Hessian求近似,需要额外的计算.

如果有可能给出梯度 (或者函数的梯度可以自动计算得到的话),这种方法通常会表现得更好. 您可以通过 Gradient 选项提供梯度,它有几种方式来 指定导数.     

这里定义一个函数,返回相对 的数值值的梯度.     
In[5]:=
Click for copyable input
Out[5]=
这里告诉 FindMinimum 使用所提供的梯度. Hessian是用梯度的有限差分计算的.
In[6]:=
Click for copyable input
Out[6]=

如果您能够提供一个产生 Hessian 的程序,您也可以提供这个. 因为该Hessian只用于牛顿法,所以它作为 的一个方法选项被提供.

这里定义一个函数,返回 的数值的Hessian.
In[7]:=
Click for copyable input
Out[7]=
这里告诉 FindMinimum 使用所提供的梯度和Hessian.
In[8]:=
Click for copyable input
Out[8]=

原则上,牛顿方法要么使用通过对符号式导数求值而计算的Hessian,要么使用通过有限差分来计算得到的Hessian. 然而,这种计算方法的收敛决定于我们的函数是凸函数,此时Hessian总是正定的. 我们经常可以见到搜索过程从违反这一条件的位置起始,所以算法必须考虑到这种可能性.     

下面的例子中搜索是从一个局部极大值附近开始的.     
In[9]:=
Click for copyable input
Out[9]=

当足够接近一个局部极大值的时候,Hessian实际上是负定的.     

这里计算局部极大值附近的Hessian的特征值.     
In[10]:=
Click for copyable input
Out[10]=

在 Hessian 不是正定的情况,如果您只是应用牛顿步计算公式,就有可能得到一个不会导致函数值减少的方向.

这里计算通过求解 找到方向的方向导数. 由于它是正的,所以沿着这个方向移动将局部性地增加函数值.
In[11]:=
Click for copyable input
Out[11]=

对于线搜索方法来说,利用正定的二次模型 计算方向是至关重要的,因为由此得到的搜索过程和收敛结果依赖于一个有足够下降度的方向. 参见 "线搜索方法". Wolfram 语言通过元素足够大的对角矩阵 修正 Hessian,以使得 是正定的. 这种方法有时被称为修正牛顿方法. 对于 的修正,无论是对密集的还是稀疏的Hessian,差不多都是在按[GMW81]中所描述的Cholesky分解的计算过程中完成. 而且只有在 不是正定的情况下才做这样的修正. 如果您想独立地使用这种分解方法的话,它可以通过 LinearSolve 访问.     

这里使用 来计算步骤,其中 是在计算 Hessian 的 Cholesky 因子时决定的.
In[12]:=
Click for copyable input
Out[12]=
所得步骤是沿着一个下降的方向的.
In[13]:=
Click for copyable input
Out[13]=

除了(修正的)牛顿法的稳健性,另一个重要方面是它的收敛速度. 一旦搜索足够接近一个局部极小值,收敛速度就可以说是 二次的,这意味着如果 是局部极小值点,那么对于某常数

在机器精度上,这并不总是有显著的差别,因为通常绝大部分步骤花在接近局部极小值上. 但是,如果您想要一个极高精度的根,牛顿方法通常因其快速收敛性而成为最好的选择.     

这里使用牛顿法计算具有非常高精度的解. 精度自适应地从机器精度(初始点的精度)提高到具有100000位数的最大工作精度. 这里 ReapSow 一起使用以保存所采取的步骤数. 并用计数器来跟踪并且打印函数计算的次数和使用的步骤数目.
In[14]:=
Click for copyable input
Out[14]=

当使用选项 WorkingPrecision->prec 时,AccuracyGoalPrecisionGoal默认值prec/2. 因此,这个例子应该找到至少精确到50000位数的极小值.

这里通过计算求解我们的搜索不断接近的极小值位置,这里采用符号化的表示方法.
In[15]:=
Click for copyable input
Out[15]=
这里计算每个步骤结束时从搜索点到极小值精确点的距离的范数.     
In[16]:=
Click for copyable input
Out[16]=

这里,我们需要比步骤数目更多的函数计算次数,这是因为 Wolfram 语言自适应地把精度从初始值的精度增加到所要求的最大 WorkingPrecision. 精确度序列的选择,是要在假定点收敛于极小值的情况下,使得在最昂贵的最后精度上做尽量少的计算. 有时候,当 Wolfram 语言改变精度时,我们有必要在更高的精度上重新计算函数.

这里是一个关于每个点的精度的表格,同时表格中还提供了误差的范数.
In[17]:=
Click for copyable input
Out[17]=

请注意,通常精度基本上是误差尺度 的两倍. 对于牛顿法,这是恰当的,因为当我们计算步骤时,根据二次收敛,误差尺度实际上将翻倍.

FindMinimum 总是起始于您所提供的初始值的精度. 因此,如果您不想让它使用自适应精度控制,您可以选择精确的或者至少有最大WorkingPrecision 的数值作为初始值.

这里在整个计算过程中只使用 100000 位数的精度求解. (警告: 这需要很长的时间才能完成.)
In[18]:=
Click for copyable input
Out[18]=

虽然这可能使用较少的函数计算,但它们都是在最高精度上进行的,因此通常情况下,自适应精度节省了很多时间. 例如,前面没有使用自适应精度的命令比起从机器精度开始计算花了超过50倍的时间.

在牛顿方法下,线搜索方法信赖域方法的步控制都得以实现. 默认值是线搜索,前面的例子中使用的都是这个方法. 但是,其中任何一个例子也都可以采用信赖域方法完成. 该方法通常每步都需要更多的数值线性代数计算,但是由于这些步骤被更好地控制,所以它可能在较少的迭代次数里达到收敛.

这里使用非约束问题的程序包来建立经典Rosenbrock函数,它有一个狭窄弯曲的谷部.
In[19]:=
Click for copyable input
Out[19]=
这里显示了对于一个Rosenbrock函数,使用信赖域牛顿方法的 FindMinimum 所采取的步骤.
In[20]:=
Click for copyable input
Out[20]=
这里显示了对于同样的函数,使用线搜索牛顿方法的 FindMinimum 所采取的步骤.
In[21]:=
Click for copyable input
Out[21]=

从这两个图示的比较当中,您可以看到由于搜索沿着谷的方向进行,因而信赖域方法更好地控制了搜索步骤,进而用较少的函数计算达到收敛.

下面的表格总结了使用牛顿方法时您可以使用的选项.

选项名
默认值
"Hessian"Automatic用于计算Hessian矩阵的表达式
"StepControl""LineSearch"如何控制步骤; 选项包括 ,或者 None

Method->"Newton" 的方法选项.