非约束最优化:用 Wolfram 语言设置最优化问题

指定导数

函数 FindRoot 有一个 Jacobian 选项;函数 FindMinimumFindMaximumFindFit 有一个 Gradient 选项;而"牛顿" 方法有一个方法选项 Hessian. 所有这些导数都用同样的基本结构来指定. 下面是用来指定导数计算方法的方式总结.    

Automatic寻找函数的一个符号化导数,当符号化导数无法找到时,则用有限差分近似
SymbolicAutomatic 相同,但在使用有限差分时,给出一个警告信息
FiniteDifference使用有限差分对导数进行近似
expression使用给定的 expression 及变量的局部数值值来计算导数

计算梯度、雅可比和 Hessian 导数的方法.

对于一个导数的基本指定就是用来计算它的方法. 然而,所有导数也都采用选项. 这些可以通过使用列表 {method,opts} 来指明. 下面是一个导数选项的总结.

选项名
默认值
"EvaluationMonitor"None当每次进行导数计算的时候,利用变量的局部值来计算的表达式. 经常用 :> 而不是 -> 指定来防止符号化的计算
"Sparse"Automatic导数的稀疏结构;可以是 AutomaticTrueFalse,或者是给出非零结构的一种模式 SparseArray
"DifferenceOrder"1当使用有限差分来计算导数时,所使用的差分阶数

计算梯度、雅可比和 Hessian 导数的选项.

我们将举例说明以上这些方法和选项如何配合使用.     

这里加载一个包含一些功用函数的程序包:
这里定义了一个函数,唯一的目的是只对数值变量值进行计算:

如果只使用 Method->"Newton"FindMinimum 会释放一条 lstol 消,这是因为,由于缺乏良好的导数信息,它不能够足够精确地解出极小值.

这里显示了 FindMinimum 所采取的步骤,这里它不得不使用有限差分来计算梯度和 Hessian 的:

下面介绍如何使用梯度选项来指定导数.

这里使用梯度的符号表达式来计算 f[x,y] 的极小值:

符号化的导数不会总是有的. 如果在有限差分上,您需要更高的准确度,您可以从默认的值 1 增加差分阶数,但代价是额外的函数计算.

这里通过使用二阶有限差分计算梯度的方法来计算 f[x,y] 的极小值:

请注意,函数计算的次数要高得多,因为函数计算被用来计算梯度,进而近似计算 Hessian. (Hessian 是使用有限差分计算的,因为从提供的信息看,我们没有符号化的表达式可以用来计算. )

FindMinimumPlot 所提供的一些关于函数,梯度,和 Hessian 计算次数的信息非常有用,EvaluationMonitor 的选项使之成为可能. 下面这个例子,只是合计各种类型计算的次数. (绘图时用 ReapSow 收集那些执行过计算的值. )

这里在计算极小值的同时用计数器来记录步骤数目和函数,梯度,以及 Hessian 的计算次数:

在决定什么方法和/或方法参数对于具有相似特点的一类问题可能是最成功的时候,使用这种诊断方法会是非常有益的.

当 Wolfram 语言能够访问函数的符号化结构时,它会自动对函数及其导数做结构分析并在适当的情况下,使用 SparseArray 对象来代表导数. 由于随后的数值线性代数可以使用稀疏结构,这可以对搜索的整体效率产生深刻的影响. 当 Wolfram 语言不能做结构分析时,一般来说,它必须假定该结构是稠密的. 但是,如果您知道该导数的稀疏结构是什么,您可以使用 "Sparse" 方法选项对它进行指定,从而在导数计算(使用有限差分,计算的次数可以显著地降低)和随后的线性代数计算中,获得巨大的效率优势. 这个问题对于使用向量值的变量工作时尤为重要. 一个说明这方面问题的好例子是扩展的 Rosenbrock 问题,它有一个非常简单的稀疏结构.

这里用 UnconstrainedProblems` 程序包,获得随时可以用 FindRoot 求解的符号化表示的具有 1000 个变量的扩展的 Rosenbrock 函数:
这里使用函数的符号化表示的形式求解问题:

对于一个形式像这样简单的函数来说,很容易写出函数的向量形式,这可以比符号化表示的形式更快地计算出来,即使是使用自动编译.

这里定义了一个扩展的 Rosenbrock 函数的向量形式,它可以非常快地进行计算:
这里从问题结构中把初始点作为向量提取出来:
下面在求解问题的计算中用了一个向量变量和向量函数:

通过函数来求解,在函数计算上很快,但总体上却较慢。原因是,由于 x_List 模式使得符号化分析不透明,雅可比必须使用有限差分进行计算. 这倒不是说有限差分多么慢,而是它需要做 100 次函数计算,以获得雅可比的所有列. 如果对结构有了解,这可以减少到两次计算就可得到雅可比. 对于这个函数,雅可比的结构很简单.

这里定义了一个模式 SparseArray,它具有扩展的 Rosenbrock 函数的雅可比的非零部分的结构. (通过指定规则中的值为 _SparseArray 被当作一个 Pattern 类型的模版,正如输出形式中所显示的那样.)
这里在了解实际雅可比结构的情况下,求解该问题,并呈现出对成本的显著节约:

当一个稀疏结构给定时,我们也可以由符号化表达式计算得到值,这些值的计算对应于稀疏结构模板中给出的位置. 请注意,该值必须直接对应于 SparseArray 中排序给出的位置 (使用 ArrayRules 可以看到排序). 获得一个一致的指标排序的一个方法是对矩阵进行两次转置,这可以产生一个指标按字典顺序排序的 SparseArray.

这里对非零结构矩阵进行两次转置,以使指标排序:
这里定义了一个函数,它将返回在雅可比中对应于非零结构矩阵中指标位置的非零值:
这里用产生的稀疏符号化雅可比求解问题:

在这种情况下,使用稀疏雅可比没有显著地更快,因为该雅可比如此稀疏,以至于只要两次函数计算就可以找到它的有限差分近似;也因为该问题在接近极小值处定义得足够好,以至于雅可比的更高的准确度并没有带来任何显著的差异.

变量和初始条件

所有函数 FindMinimumFindMaximum,和 FindRoot 采取相同形式的变量指定形式. 函数 FindFit 使用相同形式的参数指定.

FindMinimum[f,vars]找出 f 关于 vars 中所给变量的一个局部极小值
FindMaximum[f,vars]找出 f 关于 vars 中所给变量的一个局部极大值,
FindRoot[f,vars]找出 关于 vars 中所给变量的一个根
FindRoot[eqns,vars]找出方程组 eqns 关于 vars 中所给变量的一个根
FindFit[data,expr,pars,vars]找到参数 pars 的值,使 expr 对作为 vars 的函数的 data 给出一个最好的拟合

"Find" 函数中的变量和参数.

列表 vars (对于 FindFitpars)应该是由单个变量指定组成的列表. 每个变量指定应该采取如下形式.

{var,st}变量 var 有初始值 st
{var,st1,st2}变量 var 有两个初始值 st1st2;第二个初始条件仅在主轴法和割线法中使用
{var,st,rl,ru}变量 var 有初始值 st; 当 var 的值超过区间 [rl,ru] 时,搜索将被终止
{var,st1,st2,rl,ru}变量 var 有两个初始值 st1st2;当 var 的值超过区间 [rl,ru] 时,搜索将被终止

"Find" 函数中单个变量的指定形式.

vars 中的变量指定都需要有相同数量的起始值. 当区域边界没有具体指定的时候,它们将被认为是无界的,即 rl=-ru=.

向量和矩阵值变量

变量最常见的用途是代表数字. 然而,变量输入语法对可视为向量,矩阵,或高阶张量处理的变量提供支持. 一般情况下,"Find" 命令,除了目前只适用于标量型变量的 FindFit,都认为变量将采用具有与其初始条件相同的矩形结构的值.

下面是一个矩阵:
这里用 FindRoot 寻找 A 一个特征值和相应的归一化特征向量:

当然,这不是计算特征值的最好方式,但是它显示了变量维度是如何从初始值中选定的. 由于 具有初始值 1,它被视为标量. 另一方面,给予 的一个初始值是一个长度为 3 的向量,所以, 它始终被作为一个长度为 3 的向量来使用.

如果您对变量使用多个初始值,则这些值必须有一致的维度,并且初始值的每个分量是不同的.

这里对每个变量使用两个初始条件寻找一个不同的特征值:

能够采用向量和矩阵值的变量的一个优势是允许您写函数,而这些函数对于较大的问题和/或自动处理不同大小的问题非常有效.

这里定义了一个函数,给出一个相当于 UnconstrainedProblems 程序包中 ExtendedRosenbrock 问题的目标函数. 该函数所期待的 值,是一个由两个行组成的矩阵.

请注意,因为除非 具有正确的结构,函数的值将会毫无意义,该定义仅限于具有该结构的变量. 例如,如果您对于任何模式 x_ 定义该函数,那么使用一个未定义的符号 x 计算(这就是 FindMinimum 所做的)会产生毫无意义的意外结果. 通常的情况是,在用到具有向量值的变量的函数时,您将不得不对这些定义进行限制. 请注意,上述的定义不排除具有正确结构的符号化的值. 例如,ExtendedRosenbrockObjective[{{x11,x12},{x21,x22}}] 对于标量 x11 给出一个函数的符号化表示方法.

这里使用 FindMinimum 求解问题,在该问题中,对于问题的大小,我们提供了一个通用值. 您可以更改 的值,而不改变其他任何东西,以求解不同大小的问题:

这里我们求得的解没有达到默认的容差,原因是 Wolfram 语言无法对函数取得符号化的导数,因此它必须求助于并不准确的有限差分法.

使用向量-和矩阵-值变量的一个缺点是 Wolfram 语言目前无法计算它们的符号化表示的导数. 有时不难开发一个产生正确导数的函数. (如果做不到这一点,如果您确实需要更高的准确性,您可以使用高阶有限差分.)

这里定义一个返回 ExtendedRosenbrockObjective 函数的梯度的函数. 请注意,该梯度是根据变量位置压平矩阵所得到的一个向量:
这里用梯度的符号化表示的值求解该问题:

雅可比和 Hessian 导数往往是稀疏的. 如果适当,您也可以对这些导数指定结构稀疏度,它可以相当程度上降低整体求解复杂度.

终止条件

从数学上来说,平滑函数的局部极小值的充分条件是很简单的: 是一个局部极小值,如果 并且 Hessian 是正定的.(Hessian 半正定是一个必要条件.)对于一个根来说,条件更简单. 然而,当函数 在一台计算机上计算时,它的值至多只在一定精度上已知,并且实际操作上只可能有有限次数的函数计算,这时我们必须使用误差估计来决定何时搜索已经足够接近一个极小值或是一个根,同时把解的计算只近似到有限容差范围内. 在大多数情况下,这些估计已足够好,但在某些情况下,通常因为该函数在细小尺度内没有被分辨的行为,它们可以是错误的.

容差影响在搜索结束前搜索可以接近一个根或局部极小值的程度. 假设函数本身有一些误差(这在使用数值值进行计算时是典型的),通常不可能比所用的数字精确度的一半更好地定位一个极小值. 这是因为局部极小值的二次性质. 在接近抛物线的底部位置,当您从极小值移动,高度变化得相当缓慢. 因此,如果在函数中有任何误差噪音,它通常会在大致相当于噪音平方根的宽度范围内掩盖抛物线的实际增长. 这可以通过一个例子很好地说明.    

这里加载一个包含一些功用函数的程序包:
下面的命令在不断递减的范围内,画出显示函数 极小值的一系列图形. 使用机器精度计算的曲线用黑色显示;实际曲线(用 100 位数的精度算得)用蓝色显示:

从这一系列的图示,很显然看出对于 量级(大约是机器精度的一半或者更小)的改变,函数中的误差掩盖了极小值附近的曲线的实际形状. 如果只在这样的精度上对函数进行采样,是无法确定对于相近的容差,一个给定的点是否给出该函数的最小局部值.

如果采用符号化表示的方法计算导数值,是更加可靠的,但在一般情况下,仅仅依靠导数值是不够的;搜索需要找到函数的一个局部极小值,而在该位置导数必须小到可以在一般情况下满足容差. 另外,请注意,如果您的函数的符号化表示的导数不能被计算,而是使用有限差分或者是一个无导数方法,解的准确度可能会进一步降低.    

根的求解可能会遭受在函数中同样的不准确度. 虽然通常情况下问题不会那么严重,一些误差估计是基于一个优值函数的,而优值函数确实有一个二次型的形状.

由于这种限制的原因,Find 函数的默认容差都将被设置为最后工作精度的一半大小. 根据函数具有多大误差,这可能可以实现也可能无法实现,但是在多数情况下,这是一个合理的目标. 您可以使用 AccuracyGoalPrecisionGoal 选项调整容差. 当 AccuracyGoal->agPrecisionGoal->pg,这就定义了容差 .

给定后,FindMinimum 尝试找到一个满足 . 当然,由于极小值的确切位置 是未知的,因而许要对 进行估计. 对此一般的做法是根据过去的步骤和导数值. 为了在极小值可以符合导数条件,我们对它强加额外的要求 . 对于 FindRoot,相应的条件只是在根的位置,残差要小:.

这里找到了 ,它满足至少 12 位的准确度,或在容差 之内. 的精度目标意味着 ,因此它没有对公式产生任何影响. (注:您不能同样地把精度目标设为 ,因为它总是用于残差的大小.)
下面表明,结果满足所要求误差容差:
这里尝试在 8 位数的准确度内找到函数 的极小值. FindMinimum 给出一个警告信息,这是由于从图中我们可以看到的函数的误差.    
这里显示了,虽然在极小值处找到的值基本上是在机器精度(machine epsilon)内,位置却发现只达到了 左右的量级:

在多维空间里,情况就更为复杂,因为在一些方向可能有比其他方向更多的误差,比如像在 FreudensteinRoth 问题 中,当极小值在沿着相对狭窄的谷处找到时. 对于这样的搜索,搜索参数经常被缩放,这进而又影响了误差估计. 尽管如此,典型的情况仍然是,极小值的二次形状影响实际可以实现的容差.

当您需要在超出默认容差的范围内寻找一个根或极小值时,您可能需要增加最后工作精度. 为此,您可以需要使用WorkingPrecision 选项. 当您使用 WorkingPrecision->prec 时,搜索从初始值的精度开始,并随着搜索的收敛,自适应地增加到 prec. 默认情况下,WorkingPrecision->MachinePrecision,所以机器精度数字被使用,这通常要快得多. 达到更高的精度花费的时间可能会显著增加,但是可能带来更准确的结果,如果您的函数定义适当的话. 对于非常高精度的解,推荐用牛顿方法,因为它的二次收敛速度大大减少了最终需要的步骤数目.

值得注意的是,如果函数用较低精度的数字定义,WorkingPrecision 选项则不会带来任何好处. 一般而言,为了使 WorkingPrecision->prec 是有效的,用来定义函数的数字应该是精确的或者至少达到精度 prec. 在可能的时候,函数中的数字精度用 SetPrecision 被人为地增加到 prec,因此收敛仍然可以达到,但是这并非总是可能的. 在任何情况下,当函数和导数使用数值计算的时候,如果必要的话,结果的精度都提高到 prec,这样使得内部算法可以使用 prec 位数的精度来完成. 即使如此,根或者极小值的实际精度或准确度受函数中的准确度限制. 记住这一点在使用 FindFit 时尤为重要,这种情况下数据通常只在一定精度上已知.

这是使用机器精度的数字定义的函数:
即使有较高的工作精度,极小值也不能被更好地分辨,因为实际函数仍然有如图中显示的同样的误差. 这里对导数做了具体指定,以保持其他信息与前面所示的机器精度下的计算的一致性:
这是当函数没有机器精度数字时,使用 20 位数精度的计算:

如果您指定 WorkingPrecision->prec,但是没有明确指定 AccuracyGoalPrecisionGoal 选项,那么它们的 Automatic 的默认设置自动被视为 AccuracyGoal->prec/2PrecisionGoal->prec/2. 正如前文所述,这将导致一般现实情况下可以预期的最小容差.

以下是使用 50 位数精度进行的计算,这里没有明确指定 AccuracyGoal 或者 PrecisionGoal 选项的设置:
这里显示了在极小值处的数值实际上被发现比默认的 25 位容差更好:

下表总结了影响精度和容差的选项.

选项名
默认值
WorkingPrecisionMachinePrecision使用的最终工作精度 prec; 精度自适应地从 prec 和初始条件的精度中较小的那个增加到 prec
AccuracyGoalAutomatic设置 ag 通过tola=10-ag 决定一个绝对容差; 当 Automatic 时,ag=prec/2
PrecisionGoalAutomatic设置 pg 通过 tolr=10-pg 决定一个绝对容差; 当 Automatic 时,pg=prec/2

"Find" 函数中的精度和容差选项.

一个搜索过程有时候会缓慢地收敛. 为了防止缓慢的搜索无限制地持续下去,所有 Find 命令都有一个在终止前所允许的最大迭代(步骤)数目. 这可以通过具有默认值 MaxIterations->100 的选项 MaxIterations 控制. 当在此条件下,搜索终止的时候,命令将会发出 cvmit 讯息.

这里从 Optimization`UnconstrainedProblems` 程序包取得 Brown-Dennis 问题.     
这里试图用默认方法,即 Levenberg-Marquardt 方法,求解问题,因为这里函数是一个平方和的形式:

Levenberg-Marquardt 方法在这个问题收敛缓慢,因为在极小值附近,残差不为零,需要 Hessian 的二阶部分. 虽然该方法最终仅在 400 个步骤下收敛,也许更好的选择是使用一个可能更快收敛的方法.

在更大规模的计算中,当达到迭代极限时,一种可能的解决方案是使用返回的最后搜索点,作为继续搜索的一个初始条件,最好是使用另一种方法来搜索.

符号计算

函数 FindMinimumFindMaximumFindRoot 具有 HoldAll 属性,因此对它们参数的计算具有特殊的语义. 首先,变量由第二个参数确定,然后被局部化. 其次,函数先进行符号化计算,然后处理成用以进行数值计算的有效形式. 最后,在命令的执行过程中,该函数使用不同的数值值反复计算. 下表列出了上述步骤以及进一步的说明.

决定变量处理第二个参数;如果第二个参数不是正确的形式(变量和初始值的列表),将对其进行计算,以得到正确的形式
局部化变量用与 BlockTable 类似的方法,对变量增加规则使得对它们的任何赋值将不会影响 "Find" 命令范围以外的 Wolfram 系统进程,也使得以往的赋值不影响该值(变量会在这个阶段对计算为自身)
计算函数使用变量的局部未定义的(符号)值,计算第一个参数(函数或者方程). 注:这是在 Wolfram 系统 5 中引入的一个变动,所以可能需要对运行在以前的版本上的代码进行一些调整. 如果您的函数在符号计算时不能保持您想要的函数形式或者慢得无法接受,那么您在定义函数时应该使它只对变量的数值值进行计算. 对此最简单的办法是使用 PatternTest (?),如f[x_?NumberQ]:=definition 中那样.
预处理函数对函数进行分析,以帮助确定应该使用的算法(例如,平方和 -> LevenbergMarquardt);如果可能的话,对函数进行优化和编译以使数值计算更快:对于FindRoot 这首先包括从方程到函数的过程
计算导数如果可能的话,计算所有需要的符号导数;否则,进行使用有限差分计算导数所需要的预处理
数值计算反复使用不同的数值计算函数(和在必要时的导数)

"Find" 命令中处理函数的步骤.

FindFit 不具有 HoldAll 属性,因此它的参数在命令开始之前都被计算. 但是,它使用如上所述所有步骤,除了一个例外,即它不进行计算函数,而是根据模型函数,变量,和所提供的数据构建一个函数来进行极小化.

有时候,您要防止进行符号计算,这往往是当您的函数不是一个明确的公式,而是通过运行一个程序所得的一个值. 我们通过一个例子来显示会发生什么情况,以及如何防止符号计算.

这里试图使用打靶法数值化地求解一个简单的边值问题:

由于函数的符号计算,命令失败. 当您在 Block 中计算它的时候,您可以看到是怎么回事.

这里用 xp 的一个局部(未定义)的值,对赋给 FindRoot 的函数进行计算:

当然,这根本不是该函数的目的所在;它甚至不依赖于 xp. 问题出现在,没有 xp 的一个数值值时,NDSolve 失败,所以ReplaceAll (/.) 失败,因为没有规则. First 只返回它的第一个参数,即 x[1]. 由于函数在 xp 没有数值值时是毫无意的,我们应该适当地定义它.    

这里定义一个函数,它以 x'[t]t=-1 的一个数值值的函数的形式返回值 x[1]

FindRoot 之外有一个简单的函数定义的一个优点是,它可以独立地被测试,以确保它就是您的目的所在.

这里画 fx1 的图形:

从图中,您推断出两个包围根的值,因此有可能利用布伦特方法的优势,来快速而准确地求解问题.

这里求解打靶问题:

看起来似乎符号计算只是创造了麻烦,因为您必须专门地定义函数来避免它. 然而,如果没有符号计算,Wolfram 语言很难利用它的综合了数值和符号功能的独特优势. 符号计算意味着命令可以持续利用来自符号分析带来的好处,比如算法决定,导数的自动计算,自动优化和编,以及结构分析.