DIRECT算法及其實現

葛仁磊 王亞男 王毅 張明浩

關鍵詞： 為了有效地解決傳統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：didjI3=i∈I：di=dj]

則[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.