1. 程式人生 > >【轉載】黃金比例搜尋演算法(Golden Section Search)的實現

【轉載】黃金比例搜尋演算法(Golden Section Search)的實現

出處:

https://www.codelast.com/%E5%8E%9F%E5%88%9B-%E9%BB%84%E9%87%91%E6%AF%94%E4%BE%8B%E6%90%9C%E7%B4%A2%E7%AE%97%E6%B3%95%EF%BC%88gold

 

黃金比例搜尋演算法/黃金分割法/0.618法/Golden Section Search/Golden Ratio Search 表達的都是同一個東西,它是一種精確的一維搜尋(line search)演算法。黃金比例搜尋演算法只需要計算目標函式值,而無需計算導數之類的值,因此用起來非常容易,應用廣泛。
這裡要瞎扯一下:我們都知道黃金比例是0.618,但是你可能沒注意到黃金比例還有一些神奇的特性,例如1÷0.618=1+0.618。在建築學中,也經常使用黃金比例來達到美學效果,例如,代表了全希臘建築藝術的最高水平的帕特農神廟的立面高與寬的比例為19比31,接近黃金比例。

要使用黃金比例搜尋演算法,函式 f(x) 需要在某一區間內滿足單峰unimodal)條件。那麼什麼是單峰呢?看這裡

在這種情況下,就具備使用該演算法的條件了。

有人說,我想用golden section search,但是我怎麼能保證在我的搜尋區間之內函式是單峰的呢?這就涉及到另一種技術了——劃界。試想一下,類似於f(x)=sinx這樣的函式圖形,在一個區間之內可能有多個波峰波谷,你如何能確定它在一個區間內的極值呢?靠的就是劃界。但是,這不在本文的討論範圍內,我們假設你的搜尋區間已經使得目標函式在其上是單峰的了。

怎麼那麼麻煩呢?每次只要一涉及到一個演算法,就會用到另一個演算法/技術,沒辦法,數學就是這樣層層相套。像最優化中的“單純形法”那樣對其他演算法依賴很少的演算法,應該算比較少的吧(相比之下,Powell演算法真是個麻煩的角色)。

若已知f(x)在[a,b]上單峰,則可用一個含有f(x)最小值的子區間來代替區間I。有一個方法是:選取兩個內點c,d(其中c<d),即有a<c<d<b。f(x)的單峰保證了f(c) < max {f(a), f(b)} 且 f(d) < max {f(a), f(b)}。

 

這裡又提到了另一個概念:內點。什麼是內點呢?

<定義>數學上,集合 S 的內部(又稱開核)含有所有直觀上“不在 S 的邊界上”的 S 的點。S 的內部中的點稱為 S 的內點。

現在有兩種情況需要考慮:

(1) f(c)≤f(d) ,則極小值出現在子區間 [a,d] 中,所以我們用d代替b,繼續在子區間 [a,d] 中搜索;

(2) f(d)<f(c) ,則極小值出現在子區間 [c,b] 中,所以我們用c代替a,繼續在子區間 [a,c] 中搜索。

這裡要注意了:我們選取的兩個內點不是隨意選定的,而是特殊選定的,我們人為地使[a, c]與[d, b]對稱,即b-d = c-a,從而:

(3)c=a+(1-r)(b-a)=ra+(1-r)b

(4)d=b-(1-r)(b-a)=(1-r)a+rb

此處r∈(0.5, 1)可以保證c < d

其中,r是未知的,我們要求出來。

我們想讓r值在每一個子區間內均保持不變,另外,在每一次搜尋的過程中,其中一箇舊的內點將被用作新的子區間的一個內點,而一箇舊的內點則會成為新的子區間的一個端點。這樣,每一次迭代過程中只需要找到一個新點,並且只需要計算一次新函式值。如果f(c) <= f(d)且只需計算一次新函式值的話,那麼我們就有:

也可以寫成:

 

將(3)、(4)式代入此式,得:

由於r∈(0.5, 1),故:

此r值即所謂的黃金比例。平常我們老是說“黃金分割點”,就是這個0.618了。

同理,當 f(d)<f(c) 時,也可求得同樣的 r 值。
下面我用一幅圖,從另一個角度來說明黃金比例搜尋演算法的尋優過程。

上圖來自Wiki
x1,x2,x3 是已知的三個點(函式值“高→低→高”),它們劃界了一個單峰區間,函式值分別為 f1,f2,f3 。現在的問題是:如何通過一個新點 x4 ,不斷縮小該區間的範圍,從而逼近極小值點?

首先,在最大的區間 (x2,x3) 內安插 x4 是最高效的,因為我們更傾向於相信極小值更有可能出現在大的區間內。
當 x4 的函式值為圖中 f4a 所示時(即 f4a>f2 ),此時,區間縮小為 x1,x2,x4 (函式值“高→低→高”的3個點);當 x4 的函式值為圖中 f4b 所示時(即 f4b<f2 ),此時,區間縮小為 x2,x4,x3(函式值“高→低→高”的3個點),每次迭代均按此規則進行選取。
其次, x4 應該放在什麼位置上呢?隨意放置的話,若運氣不好,有可能會導致收斂很慢。而事實證明,黃金比例是極好的,即 cb=1−0.618 。
黃金比例搜尋演算法是線性收斂的,那麼其收斂準則是什麼呢?也就是說,我們逐步逼近了極小值所在的區間,到什麼程度的時候,我們就可以“收手”,不用再搜尋下去了呢?各種參考資料上給出了不同的收斂準則,例如:
|x3−x1|<τ⋅|x2|+|x4| 
或:
|x3−x1|<τ⋅INIT ( INIT 為初始區間長度)
需要一提的是,在給定一個精度 δ 的情況下,我們可以計算出試驗的次數(即需要多少次能找到極小值點)為:
n≥lgδlg0.618+1 取整數。
例如:求函式 f(x)=x2−x+2 在 [−1,3] 上的極小值點,精度要求 δ=0.08 (即搜尋區間縮小為初始區間長的8%時停止搜尋,就認為找到了極小值)。
則:估計試驗次數 n≥lg0.08lg0.618+1=6.24 ,實際實驗次數 n=7 (按第二個收斂準則)。
最後來看一下簡單的黃金比例搜尋演算法程式碼(我根據Wiki中的程式碼修改的,更易於理解):

public static double PHI = (1 + Math.sqrt(5)) / 2;    // 0.618...
public static double RESPHI = 2 - PHI;  // 0.382...
/**
* Golden section search.
*
* @author Darran Zhang @ codelast.com
* @param x1  The left bound of current bounds.
* @param x2  An interior point between (x1, x3).
* @param x3  The right bound of current bounds.
* @param tau The error tolerance to judge the convergence, recommended value is e^0.5, where e is is the required absolute precision of f(x).
*/
public double goldenSectionSearch(double x1, double x2, double x3, double tau) {
double x4;
if (x3 - x2 > x2 - x1) {    // b > a,右半區間長度較大
x4 = x2 + RESPHI * (x3 - x2);   // x4安插在右半區間
} else {    // b <= a,左半區間長度較大
x4 = x2 - RESPHI * (x2 - x1);   // x4安插在左半區間
}
if (Math.abs(x3 - x1) < tau * (Math.abs(x2) + Math.abs(x4))) {    // 判斷是否達到收斂條件
return (x3 + x1) / 2;   // 達到收斂條件,取區間長度的一半作為極小值點輸出
}
assert (function.func(x4) != function.func(x2));
/* 未達到收斂條件,繼續搜尋 */
if (function.func(x4) < function.func(x2)) {    // x4點的函式值如圖中f4b所示時
if (x3 - x2 > x2 - x1) {    // b > a,右半區間長度較大
return goldenSectionSearch(x2, x4, x3, tau);  // 以x2,x4,x3作為新三點,繼續搜尋
} else {    // b <= a,左半區間長度較大
return goldenSectionSearch(x1, x4, x2, tau);  // 以x1,x4,x2作為新三點,繼續搜尋
}
} else {    // x4點的函式值如圖中f4a所示時
if (x3 - x2 > x2 - x1) {    // b > a,右半區間長度較大
return goldenSectionSearch(x1, x2, x4, tau);  // 以x1,x2,x4作為新三點,繼續搜尋
} else {  // b <= a,左半區間長度較大
return goldenSectionSearch(x4, x2, x3, tau);  // 以x4,x2,x3作為新三點,繼續搜尋
}
}
}

這段使用遞迴方式(並不好)實現的程式碼要配合上面的圖來看。