无需导数的局部极小化算法NEWUOA在Fortran中的使用简例

无需导数的局部极小化算法NEWUOA在Fortran中的使用简例

A simple example of using the local minimization algorithm NEWUOA without derivatives in Fortran

文/Sobereva@北京科音  2020-Mar-1


1 NEWUOA简介

数学上对多维函数做局部极小化的算法非常多,诸如simplex、BFGS、最陡下降法、共轭梯度法等等。NEWUOA(NEW Unconstrained Optimization Algorithm)是已故的剑桥大学教授Powell于2004年提出的一种不需要导数信息的非约束性局部极小化算法(这点类似于simplex),并且给了实现此算法的Fortran子程序,在此文介绍一下使用。

NEWUOA和其它极小化算法一样需要进行迭代,而且结果依赖于初猜。它不需要导数的好处很明显,很多情况函数是个黑箱,本身没有解析导数,或者解析导数代码很难写,如果用有限差分来计算导数,一方面昂贵,一方面还有数值精度层面的问题。根据Experimental Comparisons of Derivative Free Optimization Algorithms一文的对比(Google一搜就有),相对于流行的BFGS方法(假设梯度通过有限差分得到),多数情况下NEWUOA效率更高,即收敛到极小点所需要计算函数值的次数更少。NEWUOA的原理我就不在这里介绍了,大家可以看https://en.wikipedia.org/wiki/NEWUOA

Powell分享了他的Fortran77写的NEWUOA程序,后来有人写了Fortran95的wrapper(https://github.com/ralna/MJDP_software/tree/master/newuoa),我又进一步做了轻微修改使之用着更舒服,可以在这里下载:http://sobereva.com/attach/536/newuoa_module.f90。下面就基于这个文件,通过两个很简单的问题示例怎么使用NEWUOA方法解决实际问题,你会发现超级容易。


2 例:二元函数极小化

本例我们求下面这个函数的极小点位置和函数值

我们可以先自行求一下解析解。让函数对x1和x2求导分别得0,联立求解方程组,可知精确极小点位置是x1=1/3、x2=10/3。

极小化以上函数的最简Fortran代码如下

PROGRAM test_newuoa
use newuoa_module
implicit real*8 (a-h,o-z)
real*8,allocatable :: X(:)
external CALFUN

nvar=2
allocate(X(nvar))
X(:)=0
call newuoa_min(CALFUN, X, RHOBEG=0.1D0, RHOEND=1D-6, IPRINT=2, MAXFUN=50000)
END PROGRAM
        
subroutine CALFUN(X,F)
real*8 :: X(:),F
F=(X(1)-2)**2+(X(2)-3.5D0)**2+X(1)*X(2)
end subroutine

在编译的时候要将上面的源文件和newuoa_module.f90放在一起编译。newuoa_module.f90中定义了newuoa_module的module,并提供了包含newuoa_min在内的各种相关函数。运行期间newuoa_min会反复调用计算函数值的子程序CALFUN直到达到收敛。CALFUN这个函数的名字以及里面的参数名也可以用其它名字。

此例中,我们定义了含两个元素的数组X,一开始我们给它赋的值是初猜值,在调用newuoa_min子程序做NEWUOA极小化后,X里记录的就是我们要求的终值。CALFUN是计算被极小化的函数的数值的子程序,传入X数组,返回函数值F,这个函数有且只能有这两个参数。NEWUOA算法迭代过程牵扯到置信半径(rho),初值和终值分别由RHOBEG和RHOEND定义,它们直接影响收敛成功几率以及总耗时。二者都为正,且显然RHOEND<=RHOBEG。Powell建议把RHOBEG设为估计的参数最大变化量的十分之一。RHOEND控制最终X收敛时能达到的精度,显然要求精度越高就应当被设得越小。IPRINT决定newuoa_min运行过程中在屏幕上输出的信息量,0代表什么也不输出,1代表只输出最终结果,2代表运行过程中每当更新RHO的时候输出到现在为止最佳的X及相应的函数值,3代表每次调用CALFUN的时候都输出当前的X及对应的函数值。MAXFUN设的是最多调用CALFUN的次数,如果达到了这个值还没正常结束就算失败。

此例我们以X(1)=0、X(2)=0作为初猜,运行过程输出的信息如下
  New RHO = 1.0000D-02     Current number of function evaluations =    17
  Least value of F =  3.918956916835460D+00
  The corresponding X array is:
   3.772896D-01   3.340358D+00

  New RHO = 1.0000D-03     Current number of function evaluations =    22
  Least value of F =  3.916670957635140D+00
  The corresponding X array is:
   3.356529D-01   3.331668D+00

  New RHO = 1.0000D-04     Current number of function evaluations =    26
  Least value of F =  3.916666676331397D+00
  The corresponding X array is:
   3.334374D-01   3.333242D+00

  New RHO = 1.0000D-05     Current number of function evaluations =    30
  Least value of F =  3.916666666692425D+00
  The corresponding X array is:
   3.333392D-01   3.333330D+00

  New RHO = 1.0000D-06     Current number of function evaluations =    33
  Least value of F =  3.916666666666715D+00
  The corresponding X array is:
   3.333336D-01   3.333333D+00

  At the return from NEWUOA     Total times of function evaluations =    37
  Least value of F =  3.916666666666666D+00
  The corresponding X array is:
   3.333333D-01   3.333333D+00

可见每次RHO更新时两个X值都被输出,且相应的函数值F也被输出,迄今总共调用了多少次CALFUN也显示了。到最后达到收敛时,X(1)=0.3333333、X(2)=3.333333,和我们前面手动求的解析解1/3和10/3精确一致。把这两个值代入公式得到的值也正是最终显示的3.916666666。

如果你想节约耗时可以把RHOEND设大,这样正常结束时调用CALFUN的总次数会明显变少,但会发现优化出的参数精度也有所打折扣。一般来说需要较精确结果时用1E-6是比较适合的。

当前被极小化的函数很简单,所以用(0,0)这个初猜就得到了想要的结果,然而有时候被极小化的函数可能有不止一个极小点,如果想找全,或者找全局最低的,显然就需要考虑不止一个初猜了。可以结合比如Differential Evolution (DE)等全局最小化算法使用,诸如J Cheminform, 8, 57 (2016)这篇文章就把DE和NEWUOA相结合来优化电负性均衡方法(EEM)电荷的参数。


3 例:单变量求解

本例我们要对下面的方程求解x

如果以极小化的思想来考虑,实际上这就等于对|x+2cos(x)-0.5|进行极小化。因此此例的代码应当这样写:

PROGRAM test_newuoa
use newuoa_module
implicit real*8 (a-h,o-z)
real*8,allocatable :: X(:)
external :: CALFUN
  
nvar=1
allocate(X(nvar))
X(1)=0
call newuoa_min(CALFUN, X, RHOBEG=0.1D0, RHOEND=1D-6, IPRINT=1, MAXFUN=50000)

END PROGRAM

subroutine CALFUN(X,F)
real*8 :: X(:),F
F=abs(X(1)+2*cos(X(1))-0.5D0)
end subroutine

由于这次用了IPRINT=1,所以就只有最终结束时的信息输出了:
At the return from NEWUOA     Total times of function evaluations =    24
Least value of F =  2.144288913097370D-08
The corresponding X array is:
-8.379604D-01

结果为x=-0.83796,此时方程的求解精度可见已经很好了,x+2cos(x)与0.5之间仅相差2.1442889E-8。


4 总结

本文简单介绍了实现NEWUOA局部极小化算法的Fortran子程序的使用,通过实例,充分体现出使用这种算法相当便利,值得大家在平时的研究中尝试使用。此算法不仅可以用于参数较少的情况,用于甚至含有几百个参数的情况也同样可以。但是参数越多时,达到收敛所需计算函数值的次数就会相应地明显越多。另外,如果被优化的函数已经有解析梯度计算代码了,那么NEWOUA的价值就不大了,因为此时用BFGS达到收敛所需的代价一般会更少。总的来说,笔者感觉此算法目前受到的关注程度相对于其价值来说还偏低,以后估计会受到越来越多的重视。顺带一提,Powell还开发过用于带约束条件的无需导数的局部极小化子程序,见https://zhangzk.net/software.html