DIRECT算法及其实现
2021-06-11葛仁磊王亚男王毅张明浩
葛仁磊 王亚男 王毅 张明浩
关键词： 为了有效地解决传统Lipschitz函数极值求解方法向高维推广时遇到的复杂度及计算量急剧增加的问题，DIRECT算法使用超矩形的中心点替代顶点，降低了算法的复杂度及运算量。在进行初始化后，迭代进行“查找潜在超矩形”和“细分超矩形”，直至运算结果满足要求。本文从实用角度出发，对该算法的原理及实现的细节进行了全面介绍，并对算法进行了实现，验证了算法的可行性。
关键词： DIRECT算法; 全局优化; 细分超矩形; 潜在超矩形
中图分类号：TP312 文獻标识码：A 文章编号：1006-8228（2021）05-01-05
DIRECT algorithm and its Python implementation
Ge Renlei， Wang Yanan， Wang Yi， Zhang Minghao
（Offshore Oil Engineering Co.， Ltd.， Qingdao， Shandong 266500， China）
Abstract： In order to effectively solve the problem that the complexity and calculation increase sharply when the traditional Lipschitz function extremum solving method is extended to high dimension， DIRECT algorithm uses the center point of hyper-rectangles to replace the vertices， which reduces the complexity and calculation of the algorithm. After initialization， iteratively "finding potentially optimal hyper-rectangles" and "subdividing the hyper-rectangles" until the calculation results meet the requirements. In this paper， from the practical point of view， the principle and implementation details of the algorithm are comprehensively introduced， and the algorithm is implemented to verify its feasibility.
Key words： DIRECT algorithm; global optimization; subdividing the hyper-rectangle; potentially optimal hyper-rectangle
0 引言
DIRECT（DIvidingRECTangles，细分矩形法）算法是属于Lipschitz函数在有限区间内极值问题的一种求解方法，其本质是一种采样算法。DIRECT算法不需要求解目标函数的参数，非常适合于黑箱函数的优化仿真[1]。本文抛开DIRECT算法的技术与理论细节，展示其具体实现，使得读者能够对照本文快速实现应用，并在文章最后给出实际案例以供参考。
1 Lipschitz函数
函数[fx]为闭区间[a，b]上的函数，如果存在正常数[K]使[fx]满足：
[fx-fx'≤Kx-x' x，x'∈a，b] ⑴
则称[fx]为区间[a，b]上的Lipschitz函数，[K]为Lipschitz常数，所有满足条件中的[K]中最小的为最佳Lipschitz常数。一阶可导且有界的函数都是Lipschitz连续的[2]。根据定义⑴可得，对任意[x∈a，b]都有：
[ fx≥fa-Kx-a fx≥fb+Kx-b] ⑵
则直线[y=fa-Kx-a]与[y=fb+Kx-b]都位于函数曲线下方。两条直线的交点[x1， gK， fa， fb]为函数[fx]的下界值。
使用相同的方法可以求得在[[a， x1]]和[[x1， b]]两个区间内下[fx]的下界值。依此类推，区间的划分和最小值的求解可以一直迭代下去，直到下界值无限逼近[fx]的最小值。如图1所示。
2 DIRECT算法
传统的算法对于一维的Lipschitz函数最小值的求解问题非常有效，但不便于向n维空间推广。对一个n维的Lipschitz函数，进行一次最小值的迭代就需要对其[2n]个边界点的函数值进行计算，效率较低。
DIRECT算法改进了此问题。令[x1=m+n/2]并带入（1）式可得[6]：
[ fx≥fx1+Kx-x1，x≤x1 fx≥fx1-Kx-x1，x≥x1] ⑶
式⑶的几何意义如图2步骤（a）所示。后续对每一个半区间使用相同的方法进行迭代（图2步骤（b）），就可以使下界值逐渐逼近[fx]的最小值。
为了便于后文的讨论，此处引入以下与DIRECT算法相关的概念。
超矩形：n维空间内的闭区间。
潜在超矩形：可能存在函数最小值的超矩形。
细分超矩形：将超矩形划分成更小的超矩形的过程。
在不考虑效率的前提下，我们对所有超矩形都进行细分，再计算细分区间的下界值就可以无限逼近函数的最小值。实际上只有满足一定条件的超矩形才是潜在超矩形，在此前提下，我们只需要对潜在超矩形进行细分就能更快的逼近最小值。所以“寻找潜在超矩形”及“细分超矩形”是DIRECT算法的两个主要任务。
2.1 寻找潜在超矩形
假设[ε]为一较小正常数，[fmin]为当前分割下的最优解，[ I]为所有超矩形索引的集合，[di]为超矩形中心到超矩形顶点的欧氏距离。如果[j∈I]为潜在矩形，首先根据[I]内元素与[j]的关系将[I]划分成三个集合：
[I1=i∈I：di
则[j]满足以下三个条件[3]：
条件一：
[fcj≤fci， i∈I3] ⑷
条件二：
[max i∈I1fcj-fcidj-di≤min i∈I2fcj-fcidj-di] ⑸
条件三：当[fmin≠0]时：
[ε≤fmin-fcifmin+djfminmin i∈I2fci-fcjdi-dj] ⑹
当[fmin=0]时：
[fcj≤djmin i∈I2fci-fcjdi-dj] ⑺
[ε]被称作琼斯因子（Jones Factor）[8]，这个值为经验值。实验证明[1×10-2≤ε≤1×10-7]时对计算影响不大，[ε]的推荐值为[1×10-4]。以上引理的正确性本文不再论述，如读者对引理的证明感兴趣，可到引文[3]第2.4.3节中查看详细的证明过程。
2.2 细分超矩形
为了使超矩形在各个维度上的边长都能不断减小，保证分割的效率，超矩形的细分只在一条或多条最长边上进行[4]。图3、图4和图5展示了二维超矩形和三维超矩形的一次分割方法。
2.3 迭代终止条件
在引文[1，2，5，6]中，都使用了当迭代次数达到预设值T时，结束该算法。这种方法在规模特定的问题中得到了很好的应用[5-6]。在引文[3]中，对停止迭代条件进行了更深入的探讨，同时也建议将目标函数的评估次数及进行超矩形分割的深度当作终止迭代的条件，这种做法在采样算法中经常用到。
在实际应用中，特别是工程应用中函数的全局区间最小值[fglobal]往往是未知的。而且正如前文所述，DIRECT算法无法对函数的参数及Lipschitz常数进行估计，所以在[fglobal]未知时也无法对[fglobal]进行估计。这导致在实际应用中更常用、更有效的为“精度”或“百分精度”很难作为DIRECT算法的终止条件。
在寻找Lipschitz函数的最小值过程中，除了获得函数最小值外，有时函数在何处取得最小值对用户来说也非常重要。用户自然的会以[fmin]所在超矩形尺寸达到某一特定要求作为终止迭代的条件[8]。在后面的DIRECT算法实现中，我们将以[fmin]所在的超矩形最长被细分的次数定义为深度，并以深度小于等于给定值作为终止条件之一。
表1总结了DIRECT算法常用终止条件。
3 性能优化
第2节给定的算法逻辑清晰有效，但效率有待优化。本节从超矩形细分维度次序、边长的指数式存储和距离的存储与计算三个方面讨论如何提高DIRECT算法的效率。
3.1 超矩形的细分维度次序
2.2节中进行超矩形的细分讨论时对哪个维度优先进行细分并没有讨论。但实际情况不同维度优先会造成不同的效果，如图6所示。
图6中（a）与（b）显示出相同的中心点进入了不同尺寸的超矩形中。如在细分（a）中，左侧中心点最终进入了边长为[13×13]的超矩形中。而在细分（b）中，该中心点进入了边长为[13×1]超矩形中。这种差别会影响后续潜在超矩形的查找及下一步分割的工作量。
那么（a）与（b）两种方式，哪个更有利于后续运算呢？我们使用的策略是将最优的函數值放入到更大的超矩形中，这是由于较大的超矩形总是能够更快的被再次分割。经过实验也发现，这样的分割策略的确能够提高算法速度[6]。
3.2 边长的指数式存储
第2.2节细分超矩形的过程（图2、图3、图4）中，每次超矩形的最长边都会变成原来边长的[13]。即一条边如果经历[i]次分割，那么边的长度会变成原边长的[3-i]。
为了便于边长的存储，在算法开始前将边长不规则的初始超矩形映射到所有边长都为1的超正方形内。此后就可以使用指数[i]来替代原有边长[3-i]来进行存储、排序等操作，可以节省空间并提高运算效率。
3.3 距离的存储与计算
在寻找潜在超矩形时，公式⑸⑹⑺都用到了中心到顶点的距离。DIRECT算法在提出时[6]，作者很自然的使用了欧氏距离。欧氏距离是二范数距离，其计算涉及平方和开方运算，计算量大，且产生的结果为实数。为了保证计算精度，往往需要较大的存储空间。这些原因促使研究者转而研究其他类型的距离是否也能产生相同的效果，以及其他类型的距离能不能加快运算速度，节省存储空间。
在引文[7]中，作者提出并证明了使用无穷范数距离可以达到与欧氏距离相同的效果。超矩形中心到超矩形顶点的无穷范数距离可以表示为：
[dist∞=3-l/2]
其中[l]表示超矩形最长边被分割的次数，无穷范数距离本质上就是最长边的一半。使用无穷范数距离有如下好处。
⑴ 距离的表示与存储可以直接使用自然数[l]，减少存储空间。
⑵ 距离的计算转化为边长排序，减少计算量。
⑶ 在第2.1节的运算中，进行条件一运算时，能够得到更少的潜在超矩形，减少后续查找潜在超矩形的运算量。而且会进一步减少中心点处函数值的运算次数。
为了比较二范数与无穷范数的效率，我们使用GP函数[8，13]这个例子，分别使用欧氏距离和无穷范数距离两种方法进行运算，得到统计数据如表2所示。
從表2数据可以看出，使用无穷范数作为距离度量时，虽然有时迭代次数会超过使用二范数作为距离度量，但在目标函数的计算次数上有很大优势。所以本文给出的示例代码中以无穷范数作为距离度量。
4 DIRECT算法实现
4.1 算法流程图
DIRECT算法总体流程图7所示。
图7中步骤a即为第2.2节“查找潜在超矩形”，其细节请参照图8。步骤b和步骤c合在一起为2.1节“细分超矩形”，其细节请参照图9。
4.2 DIRECT算法算例
按照4.1节提供的流程图，笔者用Python编写了DIRECT.py。我们使用引文[4，9]中采用的GoldsteinPrice（GP）测试函数：
[fx=1+x1+x2+1219-14x1+3x21-14x2+6x1x2+3x22]
[×][30+2x1-3x2218-32x1+12x21+48x2-36x1x2+27x22]
同样，采用与引文[4]相同的范围[-2≤xi≤2]，
[i∈1，2]，图10、图11表示了在此区间范围内的图像：
5 结束语
DIRECT算法稳定性强，适用性广泛，适合在对性能、精度要求不高的情况下使用。本文从实用角度，全面而详细的介绍了DIRECT算法，提供了一种开箱可用的算法，并给出了详细的算法流程图及其实现，读者可以按照自身需求直接使用。
参考文献（References）：
[1] 安华萍，李文静.DIRECT优化算法计算机仿真研究[J].惠州学院学报，2019.3：16
[2] Sohrab H H. Basic real analysis[M]. Boston， Basel， Berlin：Birkh?user，2003.
[3] Jones D R， Perttunen C D， Stuckman B E. Lipschitzianoptimization without the Lipschitz constant[J].Journal of optimization Theory and Applications，1993.79（1）：157-181
[4] Jorg Maximilian XaverGablonsky. Modifications of thedirect algorithm[M]，2001.
[5] Finkel D. DIRECT optimization algorithm user guide[R].North Carolina State University. Center for Research in Scientific Computation，2003.
[6] 宋科康，冯文涛.采用DIRECT算法的外辐射源雷达高效直接定位方法[J].信号处理，2020.
[7] 王云宏.基于DIRECT算法的微震震源快速网格搜索定位方法研究[J].地球物理学进展，2016.31（4）：1700-1708
[8] 武汉大学测绘学院测量平差学科组. 误差理论与测量平差基础-第2版[M].武汉大学出版社，2009.
[9] Cox S E， Haftka R T， Baker C， et al. Globalmultidisciplinary optimization of a high speed civil transport[C]//Proc.Aerospace Numerical Simulation Symposium，1999.99：23-28
[10] Dixon L，Szego G.Towards Global Optimisation 2 North-Holland[J]，1978.
