非约束最优化:用 Wolfram 语言设置最优化问题
指定导数
函数 FindRoot 有一个 Jacobian 选项;函数 FindMinimum、FindMaximum 和 FindFit 有一个 Gradient 选项;而"牛顿" 方法有一个方法选项 Hessian. 所有这些导数都用同样的基本结构来指定. 下面是用来指定导数计算方法的方式总结.
Automatic | 寻找函数的一个符号化导数,当符号化导数无法找到时,则用有限差分近似 |
Symbolic | 与 Automatic 相同,但在使用有限差分时,给出一个警告信息 |
FiniteDifference | 使用有限差分对导数进行近似 |
expression | 使用给定的 expression 及变量的局部数值值来计算导数 |
对于一个导数的基本指定就是用来计算它的方法. 然而,所有导数也都采用选项. 这些可以通过使用列表 {method,opts} 来指明. 下面是一个导数选项的总结.
选项名
|
默认值
| |
"EvaluationMonitor" | None | 当每次进行导数计算的时候,利用变量的局部值来计算的表达式. 经常用 :> 而不是 -> 指定来防止符号化的计算 |
"Sparse" | Automatic | 导数的稀疏结构;可以是 Automatic,True,False,或者是给出非零结构的一种模式 SparseArray |
"DifferenceOrder" | 1 | 当使用有限差分来计算导数时,所使用的差分阶数 |
如果只使用 Method->"Newton",FindMinimum 会释放一条 lstol 消息,这是因为,由于缺乏良好的导数信息,它不能够足够精确地解出极小值.
符号化的导数不会总是有的. 如果在有限差分上,您需要更高的准确度,您可以从默认的值 1 增加差分阶数,但代价是额外的函数计算.
请注意,函数计算的次数要高得多,因为函数计算被用来计算梯度,进而近似计算 Hessian. (Hessian 是使用有限差分计算的,因为从提供的信息看,我们没有符号化的表达式可以用来计算. )
FindMinimumPlot 所提供的一些关于函数,梯度,和 Hessian 计算次数的信息非常有用,EvaluationMonitor 的选项使之成为可能. 下面这个例子,只是合计各种类型计算的次数. (绘图时用 Reap 和 Sow 收集那些执行过计算的值. )
在决定什么方法和/或方法参数对于具有相似特点的一类问题可能是最成功的时候,使用这种诊断方法会是非常有益的.
当 Wolfram 语言能够访问函数的符号化结构时,它会自动对函数及其导数做结构分析并在适当的情况下,使用 SparseArray 对象来代表导数. 由于随后的数值线性代数可以使用稀疏结构,这可以对搜索的整体效率产生深刻的影响. 当 Wolfram 语言不能做结构分析时,一般来说,它必须假定该结构是稠密的. 但是,如果您知道该导数的稀疏结构是什么,您可以使用 "Sparse" 方法选项对它进行指定,从而在导数计算(使用有限差分,计算的次数可以显著地降低)和随后的线性代数计算中,获得巨大的效率优势. 这个问题对于使用向量值的变量工作时尤为重要. 一个说明这方面问题的好例子是扩展的 Rosenbrock 问题,它有一个非常简单的稀疏结构.
对于一个形式像这样简单的函数来说,很容易写出函数的向量形式,这可以比符号化表示的形式更快地计算出来,即使是使用自动编译.
通过函数来求解,在函数计算上很快,但总体上却较慢。原因是,由于 x_List 模式使得符号化分析不透明,雅可比必须使用有限差分进行计算. 这倒不是说有限差分多么慢,而是它需要做 100 次函数计算,以获得雅可比的所有列. 如果对结构有了解,这可以减少到两次计算就可得到雅可比. 对于这个函数,雅可比的结构很简单.
当一个稀疏结构给定时,我们也可以由符号化表达式计算得到值,这些值的计算对应于稀疏结构模板中给出的位置. 请注意,该值必须直接对应于 SparseArray 中排序给出的位置 (使用 ArrayRules 可以看到排序). 获得一个一致的指标排序的一个方法是对矩阵进行两次转置,这可以产生一个指标按字典顺序排序的 SparseArray.
在这种情况下,使用稀疏雅可比没有显著地更快,因为该雅可比如此稀疏,以至于只要两次函数计算就可以找到它的有限差分近似;也因为该问题在接近极小值处定义得足够好,以至于雅可比的更高的准确度并没有带来任何显著的差异.
变量和初始条件
所有函数 FindMinimum,FindMaximum,和 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 给出一个最好的拟合 |
列表 vars (对于 FindFit 为 pars)应该是由单个变量指定组成的列表. 每个变量指定应该采取如下形式.
{var,st} | 变量 var 有初始值 st |
{var,st1,st2} | 变量 var 有两个初始值 st1 和 st2;第二个初始条件仅在主轴法和割线法中使用 |
{var,st,rl,ru} | 变量 var 有初始值 st; 当 var 的值超过区间 [rl,ru] 时,搜索将被终止 |
{var,st1,st2,rl,ru} | 变量 var 有两个初始值 st1 和 st2;当 var 的值超过区间 [rl,ru] 时,搜索将被终止 |
vars 中的变量指定都需要有相同数量的起始值. 当区域边界没有具体指定的时候,它们将被认为是无界的,即 rl=-∞,ru=∞.
向量和矩阵值变量
变量最常见的用途是代表数字. 然而,变量输入语法对可视为向量,矩阵,或高阶张量处理的变量提供支持. 一般情况下,"Find" 命令,除了目前只适用于标量型变量的 FindFit,都认为变量将采用具有与其初始条件相同的矩形结构的值.
当然,这不是计算特征值的最好方式,但是它显示了变量维度是如何从初始值中选定的. 由于 具有初始值 1,它被视为标量. 另一方面,给予 的一个初始值是一个长度为 3 的向量,所以, 它始终被作为一个长度为 3 的向量来使用.
如果您对变量使用多个初始值,则这些值必须有一致的维度,并且初始值的每个分量是不同的.
能够采用向量和矩阵值的变量的一个优势是允许您写函数,而这些函数对于较大的问题和/或自动处理不同大小的问题非常有效.
请注意,因为除非 具有正确的结构,函数的值将会毫无意义,该定义仅限于具有该结构的变量. 例如,如果您对于任何模式 x_ 定义该函数,那么使用一个未定义的符号 x 计算(这就是 FindMinimum 所做的)会产生毫无意义的意外结果. 通常的情况是,在用到具有向量值的变量的函数时,您将不得不对这些定义进行限制. 请注意,上述的定义不排除具有正确结构的符号化的值. 例如,ExtendedRosenbrockObjective[{{x11,x12},{x21,x22}}] 对于标量 x11,… 给出一个函数的符号化表示方法.
这里我们求得的解没有达到默认的容差,原因是 Wolfram 语言无法对函数取得符号化的导数,因此它必须求助于并不准确的有限差分法.
使用向量-和矩阵-值变量的一个缺点是 Wolfram 语言目前无法计算它们的符号化表示的导数. 有时不难开发一个产生正确导数的函数. (如果做不到这一点,如果您确实需要更高的准确性,您可以使用高阶有限差分.)
雅可比和 Hessian 导数往往是稀疏的. 如果适当,您也可以对这些导数指定结构稀疏度,它可以相当程度上降低整体求解复杂度.
终止条件
从数学上来说,平滑函数的局部极小值的充分条件是很简单的: 是一个局部极小值,如果 并且 Hessian 是正定的.(Hessian 半正定是一个必要条件.)对于一个根来说,条件更简单. 然而,当函数 在一台计算机上计算时,它的值至多只在一定精度上已知,并且实际操作上只可能有有限次数的函数计算,这时我们必须使用误差估计来决定何时搜索已经足够接近一个极小值或是一个根,同时把解的计算只近似到有限容差范围内. 在大多数情况下,这些估计已足够好,但在某些情况下,通常因为该函数在细小尺度内没有被分辨的行为,它们可以是错误的.
容差影响在搜索结束前搜索可以接近一个根或局部极小值的程度. 假设函数本身有一些误差(这在使用数值值进行计算时是典型的),通常不可能比所用的数字精确度的一半更好地定位一个极小值. 这是因为局部极小值的二次性质. 在接近抛物线的底部位置,当您从极小值移动,高度变化得相当缓慢. 因此,如果在函数中有任何误差噪音,它通常会在大致相当于噪音平方根的宽度范围内掩盖抛物线的实际增长. 这可以通过一个例子很好地说明.
从这一系列的图示,很显然看出对于 量级(大约是机器精度的一半或者更小)的改变,函数中的误差掩盖了极小值附近的曲线的实际形状. 如果只在这样的精度上对函数进行采样,是无法确定对于相近的容差,一个给定的点是否给出该函数的最小局部值.
如果采用符号化表示的方法计算导数值,是更加可靠的,但在一般情况下,仅仅依靠导数值是不够的;搜索需要找到函数的一个局部极小值,而在该位置导数必须小到可以在一般情况下满足容差. 另外,请注意,如果您的函数的符号化表示的导数不能被计算,而是使用有限差分或者是一个无导数方法,解的准确度可能会进一步降低.
根的求解可能会遭受在函数中同样的不准确度. 虽然通常情况下问题不会那么严重,一些误差估计是基于一个优值函数的,而优值函数确实有一个二次型的形状.
由于这种限制的原因,Find 函数的默认容差都将被设置为最后工作精度的一半大小. 根据函数具有多大误差,这可能可以实现也可能无法实现,但是在多数情况下,这是一个合理的目标. 您可以使用 AccuracyGoal 和 PrecisionGoal 选项调整容差. 当 AccuracyGoal->ag 和 PrecisionGoal->pg,这就定义了容差 和 .
当 和 给定后,FindMinimum 尝试找到一个满足 值 . 当然,由于极小值的确切位置 是未知的,因而许要对 进行估计. 对此一般的做法是根据过去的步骤和导数值. 为了在极小值可以符合导数条件,我们对它强加额外的要求 . 对于 FindRoot,相应的条件只是在根的位置,残差要小:.
在多维空间里,情况就更为复杂,因为在一些方向可能有比其他方向更多的误差,比如像在 Freudenstein–Roth 问题 中,当极小值在沿着相对狭窄的谷处找到时. 对于这样的搜索,搜索参数经常被缩放,这进而又影响了误差估计. 尽管如此,典型的情况仍然是,极小值的二次形状影响实际可以实现的容差.
当您需要在超出默认容差的范围内寻找一个根或极小值时,您可能需要增加最后工作精度. 为此,您可以需要使用WorkingPrecision 选项. 当您使用 WorkingPrecision->prec 时,搜索从初始值的精度开始,并随着搜索的收敛,自适应地增加到 prec. 默认情况下,WorkingPrecision->MachinePrecision,所以机器精度数字被使用,这通常要快得多. 达到更高的精度花费的时间可能会显著增加,但是可能带来更准确的结果,如果您的函数定义适当的话. 对于非常高精度的解,推荐用牛顿方法,因为它的二次收敛速度大大减少了最终需要的步骤数目.
值得注意的是,如果函数用较低精度的数字定义,WorkingPrecision 选项则不会带来任何好处. 一般而言,为了使 WorkingPrecision->prec 是有效的,用来定义函数的数字应该是精确的或者至少达到精度 prec. 在可能的时候,函数中的数字精度用 SetPrecision 被人为地增加到 prec,因此收敛仍然可以达到,但是这并非总是可能的. 在任何情况下,当函数和导数使用数值计算的时候,如果必要的话,结果的精度都提高到 prec,这样使得内部算法可以使用 prec 位数的精度来完成. 即使如此,根或者极小值的实际精度或准确度受函数中的准确度限制. 记住这一点在使用 FindFit 时尤为重要,这种情况下数据通常只在一定精度上已知.
如果您指定 WorkingPrecision->prec,但是没有明确指定 AccuracyGoal 和 PrecisionGoal 选项,那么它们的 Automatic 的默认设置自动被视为 AccuracyGoal->prec/2 和 PrecisionGoal->prec/2. 正如前文所述,这将导致一般现实情况下可以预期的最小容差.
选项名 | 默认值 | |
WorkingPrecision | MachinePrecision | 使用的最终工作精度 prec; 精度自适应地从 prec 和初始条件的精度中较小的那个增加到 prec |
AccuracyGoal | Automatic | 设置 ag 通过tola=10-ag 决定一个绝对容差; 当 Automatic 时,ag=prec/2 |
PrecisionGoal | Automatic | 设置 pg 通过 tolr=10-pg 决定一个绝对容差; 当 Automatic 时,pg=prec/2 |
一个搜索过程有时候会缓慢地收敛. 为了防止缓慢的搜索无限制地持续下去,所有 Find 命令都有一个在终止前所允许的最大迭代(步骤)数目. 这可以通过具有默认值 MaxIterations->100 的选项 MaxIterations 控制. 当在此条件下,搜索终止的时候,命令将会发出 cvmit 讯息.
Levenberg-Marquardt 方法在这个问题收敛缓慢,因为在极小值附近,残差不为零,需要 Hessian 的二阶部分. 虽然该方法最终仅在 400 个步骤下收敛,也许更好的选择是使用一个可能更快收敛的方法.
在更大规模的计算中,当达到迭代极限时,一种可能的解决方案是使用返回的最后搜索点,作为继续搜索的一个初始条件,最好是使用另一种方法来搜索.
符号计算
函数 FindMinimum、FindMaximum 和 FindRoot 具有 HoldAll 属性,因此对它们参数的计算具有特殊的语义. 首先,变量由第二个参数确定,然后被局部化. 其次,函数先进行符号化计算,然后处理成用以进行数值计算的有效形式. 最后,在命令的执行过程中,该函数使用不同的数值值反复计算. 下表列出了上述步骤以及进一步的说明.
决定变量 | 处理第二个参数;如果第二个参数不是正确的形式(变量和初始值的列表),将对其进行计算,以得到正确的形式 |
局部化变量 | 用与 Block 和 Table 类似的方法,对变量增加规则使得对它们的任何赋值将不会影响 "Find" 命令范围以外的 Wolfram 系统进程,也使得以往的赋值不影响该值(变量会在这个阶段对计算为自身) |
计算函数 | 使用变量的局部未定义的(符号)值,计算第一个参数(函数或者方程). 注:这是在 Wolfram 系统 5 中引入的一个变动,所以可能需要对运行在以前的版本上的代码进行一些调整. 如果您的函数在符号计算时不能保持您想要的函数形式或者慢得无法接受,那么您在定义函数时应该使它只对变量的数值值进行计算. 对此最简单的办法是使用 PatternTest (?),如f[x_?NumberQ]:=definition 中那样. |
预处理函数 | 对函数进行分析,以帮助确定应该使用的算法(例如,平方和 -> Levenberg–Marquardt);如果可能的话,对函数进行优化和编译以使数值计算更快:对于FindRoot 这首先包括从方程到函数的过程 |
计算导数 | 如果可能的话,计算所有需要的符号导数;否则,进行使用有限差分计算导数所需要的预处理 |
数值计算 | 反复使用不同的数值计算函数(和在必要时的导数) |
FindFit 不具有 HoldAll 属性,因此它的参数在命令开始之前都被计算. 但是,它使用如上所述所有步骤,除了一个例外,即它不进行计算函数,而是根据模型函数,变量,和所提供的数据构建一个函数来进行极小化.
有时候,您要防止进行符号计算,这往往是当您的函数不是一个明确的公式,而是通过运行一个程序所得的一个值. 我们通过一个例子来显示会发生什么情况,以及如何防止符号计算.
由于函数的符号计算,命令失败. 当您在 Block 中计算它的时候,您可以看到是怎么回事.
当然,这根本不是该函数的目的所在;它甚至不依赖于 xp. 问题出现在,没有 xp 的一个数值值时,NDSolve 失败,所以ReplaceAll (/.) 失败,因为没有规则. First 只返回它的第一个参数,即 x[1]. 由于函数在 xp 没有数值值时是毫无意的,我们应该适当地定义它.
在 FindRoot 之外有一个简单的函数定义的一个优点是,它可以独立地被测试,以确保它就是您的目的所在.
从图中,您推断出两个包围根的值,因此有可能利用布伦特方法的优势,来快速而准确地求解问题.
看起来似乎符号计算只是创造了麻烦,因为您必须专门地定义函数来避免它. 然而,如果没有符号计算,Wolfram 语言很难利用它的综合了数值和符号功能的独特优势. 符号计算意味着命令可以持续利用来自符号分析带来的好处,比如算法决定,导数的自动计算,自动优化和编,以及结构分析.