無約束最佳化方法及程式碼實現

一、本文內容

相關理論：最佳化方法、數值分析、矩陣論

程式碼：使用Julia語言，在REPL上測試

涉及演算法：一階最速下降法（梯度法）、二階牛頓法、共軛梯度法、變換輪換法、鮑威爾法單形替換法

採用數值分析的理論實現，即採用迭代法，從某一個初始值開始，沿某一方向和長度進行搜尋求解使函式值下降，即

$x^{k+1}=x^k-\alpha^k d^k$

，其中

$\alpha$

是步長，

是第

步的搜尋方向，使

$f(x^{k+1})<f(x^k)$

。問題的核心就在於如何選取搜尋方向和搜尋步長，形成各種演算法。

二、參考

矩陣分析與應用（第二版）－張賢達－清華大學出版社

數值分析（第二版）－Timothy Sauer著裴玉茹馬庚宇譯－機械工業出版社

機械最佳化設計（第4版）－孫靖民梁迎春主編－機械工業出版社

解釋不精確線搜尋中的Armijo-Goldstein準則及Wolfe-Powell準則

小市民：Armijo-Goldstein準則

數值分析6（3共軛梯度法）蘇州大學_嗶哩嗶哩（゜-゜）つロ乾杯~-bilibili

三、無約束最佳化數學模型

3。1、無約束目標函式（二次函式）

$f(x)=\frac{1}{2}x^TGx+b^Tx+c，G為對稱正定矩陣$

一階梯度向量為：

$\bigtriangledown f(x)=Gx+b$

二階海森矩陣為：

3。2、函式1

迭代起始值：

四、無約束最佳化－最速下降法（梯度法）和牛頓法

4。1、搜尋方向選擇

$f(x^{k+1})\approx f(x^k)+(\triangledown f(x^k))^T \Delta x^k \quad \& \quad (\triangledown f(x^k))^T \Delta x^k <0 \Rightarrow f(x^{k+1})<f(x^k)$

因此，下降方向與負梯度方向的夾角成銳角，即能滿足函式值下降（注：

$\Delta x^k = - \alpha^k d^k$

）。

最速下降法：採用負梯度方向

$d^k = -\bigtriangledown f(x^k)$

，僅利用一階泰勒展開求得

牛頓法：

$d^k = G^{-1}\triangledown f(x^k)$

利用二階泰勒展開求得，需要求二階海森矩陣的逆。

4。2、搜尋步長選擇

演算法執行之前，直接給出步長序列

$\{ \alpha^k\}_{k=0}^\infty$

。例如

$\alpha^k =\alpha$

固定步長，或

$\alpha^k=\frac{h}{\sqrt{k+1}}$

全松馳，即在搜尋方向上函式值下降最多，

$\alpha^k=min_\alpha f(x^k - \alpha \triangledown d^k)$

Goldstein-Armijo規則其中，

$0<\alpha < \beta <1$

是兩個固定的引數。詳細解釋見參考4，解釋的很棒。

4。3、使用梯度方向+全松馳步長的最速下降法

4。3。1、搜尋方向的正交性

由上可看出，這種方向相相鄰搜尋方向相互正交。

4。3。2 程式

使用3。2中的函式對4。3。1中的演算法程式設計（julia語言）：

#= program1。jl

使用梯度方向+全松馳演算法迭代法求函式 f（x）=x1^2+25x2^2的極小點。

迭代起始點為 x0=［2；2］。

顯然已知最小點為（0，0），按前向誤差|x*-x|的二範數最小設定終止條件。 =#

using LinearAlgebra， Plots； gr（）

# 將 f（x）轉換為矩陣形式，也可以不用轉，但求梯度方便。

G = ［2 0；0 50］；

gradient_fx（x） = G * x；

f（x） = 0。5 * x‘ * G * x；

f2（x1， x2） = x1^2 + 25x2^2；

x_next = ［2；2］；

x_list = copy（x_next）；

y_list = ［f（x_next）］；

grad_list = ［0。0；0。0］；

alpha_list = ［0。0］；

TOL = 0。0001；

while norm（x_next） > TOL

global x_next， x_list， y_list， grad_list， alpha_list；

# 計算當前點梯度

gradient_k = gradient_fx（x_next）；

# 計算步長

alpha = （x_next’ * G‘ * gradient_k） / （gradient_k’ * G‘ * gradient_k）；

# 計算下一點

x_next = x_next - alpha * gradient_k；

# 將運算過程中的值儲存

x_list = hcat（x_list， x_next）；

y_list = push！（y_list， f（x_next））；

grad_list = hcat（grad_list， gradient_k）；

alpha_list = push！（alpha_list， alpha）；

end

data = hcat（hcat（hcat（x_list’， y_list）， alpha_list）， grad_list‘）；

data = @。 round（data， digits = 4）；

println（“k \t x1 \t x2 \t f（x） \t alpha \t Gradient_x1 \t Gradient_x2”）；

row_len， col_len = size（data）；

for row in 1：row_len

print（row， “\t”）

for col in 1：col_len

print（data［row， col］， “\t”）；

end

println（“”）；

end

lims_x1 = -3：0。01：3；

lims_x2 = -3：0。01：3；

p = contour（lims_x1， lims_x2， f2， fill = true）；

display（p）

png（“contour”）；

結果

從上圖的結果中計算可得出，相鄰的梯度向量之間是正交的。

等高線圖

4。4 使用梯度方向+Armijo-Goldstein準則（以下簡寫AG準則）

4。4。1 解釋

在參考4和5中已經完美的解釋了大部分內容，但是編碼的時候發現還有有些問題，想了一個晚上加除錯發現了端倪，充分證明了實踐是檢驗真理的唯一標準（@-@）。以下是對參考4和5的補充說明。

首先對參考4進行修正說明，

注意：上圖中默認了搜尋方向

與負梯度方向

的夾角是銳角，因為只有這樣才能滿足函式值下降。參考4中只解釋了將式（1）泰勒級數展開滿足要求，但卻沒有說將（2）式展開泰勒展開。如果按上圖上把兩式都一階泰勒展開，因為

$\alpha > 0$

，顯然（1）式永遠符合要求，（2）式不符合要求，產生了悖論AG原則多寫了式（1），而且式（2）永遠不能滿足。出現問題，還是先從自身的解釋找原因比較好。

Q：只要搜尋方向與負梯度方向成銳角，函式值就下降？

A：這句話要加上一個限定條件，在當前點的某個小範圍內。因為只有在這個小的範圍內，才能滿足泰勒展開的條件，推匯出函式值下降。從另一個角度考慮，也可以說只有在上述搜尋方向上的小範圍內的步長才能滿足函式值下降。

如上圖所示，假設取搜尋方向為負梯度方向，顯示當搜尋步長在KA兩點之間的時候才能滿足函式值下降，否則函式值就會上升，因此AG準則中要求了函式值必須下降的要求。參考5中的公式：

解釋的相對比較清楚，明確了函式值的下降，也說明了不能用一階泰勒展開解釋，因為如果用一階泰勒展開，上式可等價於：

$\mu_1 \alpha_k \leq \alpha_k \leq \mu_2 \alpha_k$

顯然中間不能成立。下面先給出程式及結果，然後繼續解釋。

4。4。2 程式

#= program2。jl

使用梯度方向+Goldstein-Armijo準則求函式 f（x）=x1^2+25x2^2的極小點。

迭代起始點為 x0=［2；2］，Armijo準則引數 u1=0。2， u2=0。8， p1=0。5， p2=1。5。

顯然已知最小點為（0，0），按前向誤差|x*-x|的二範數最小設定終止條件，

也可按照|x［k+1］-x［k］|< TOL來設定終止準則。 =#

using

LinearAlgebra

# 將 f（x）轉換為矩陣形式，也可以不用轉，但求梯度方便。

［

；

］；

gradient_fx

（

）

；

（

）

0。5

’

；

x_next

［

；

］；

x_list

copy

（

x_next

）；

y_list

［

（

x_next

）］；

grad_list

［

0。0

；

0。0

］；

alpha_list

［

0。0

］；

# 設定Armijo準則引數

0。2

；

0。8

；

0。5

；

1。5

；

TOL

0。0001

；

while

norm

（

x_next

）

TOL

global

x_next

，

x_list

，

y_list

，

grad_list

，

alpha_list

；

# 計算當前點梯度

G_k

gradient_fx

（

x_next

）；

y_k

（

x_next

）；

# 設定搜尋方向

G_k

；

# 計算步長

alpha

1。

；

# 設定一個布林量控制迴圈跳出

gate_alpha

true

；

while

gate_alpha

（

alpha

G_k

‘

）

（

y_k

（

x_next

alpha

））

alpha

；

continue

；

end

（

y_k

（

x_next

alpha

））

（

alpha

G_k

’

）

alpha

；

continue

；

else

gate_alpha

false

；

end

# 計算下一點

x_next

alpha

G_k

；

# 將運算過程中的值儲存

x_list

hcat

（

x_list

，

x_next

）；

y_list

push！

（

y_list

，

（

x_next

））；

grad_list

hcat

（

grad_list

，

G_k

）；

alpha_list

push！

（

alpha_list

，

alpha

）；

end

data

hcat

（

hcat

（

hcat

（

x_list

‘

，

y_list

），

alpha_list

），

grad_list

’

）；

data

。

round

（

data

，

digits

）；

println

（

“k

f（x）

alpha

Gradient_x1

Gradient_x2”

）；

row_len

，

col_len

size

（

data

）；

for

row

：

row_len

（

row

，

“

”

）

for

col

：

col_len

（

data

［

row

，

col

］，

“

”

）；

end

println

（

“”

）；

end

結果一部分結果：

4。4。3 解釋一階泰勒的問題

首先使用下列程式碼繪製圖形：

#= image_of_f。jl

繪製等高線圖解釋Armijo-Goldstein準則，函式f（x，y）=x^2+25y^2 =#

using LinearAlgebra， Plots；gr（）

f（x，y） = x^2 + 25y^2

f1（x， z） = sqrt（z - x^2） / 5；

f2（x， z） = -sqrt（z - x^2） / 5；

# 繪製（2，2）點處的等高線

z = f（2， 2）；

x = -3：0。001：3；

plot（x， hcat（f1。（x， z）， f2。（x， z）），c = ［：blue ：blue］，legend = false， framestyle = ：zerolines）

# 繪製（2，2）到下一迭代點（1。875， -1。125）的長度

plot！（［2；1。875］，［2；-1。125］， c = ［：red］， marker = ：dot， arrow = 2）

# 繪製（1。875， -1。125）處的等高線

z = f（1。875， -1。125）；

plot！（x， hcat（f1。（x， z）， f2。（x， z）），c = ［：yellow ：yellow］，legend = false， framestyle = ：zerolines）

# 繪製函式值下降最多處（1。9199， -0。0031）點的等高線

z = f（1。9199， -0。0031）；

x_max = floor（sqrt（z）， digits = 6）；

x2 = -x_max：0。01：x_max；

plot！（x2， hcat（f1。（x2， z）， f2。（x2， z）），c = ［：black ：black］，legend = false， framestyle = ：zerolines）

# 繪製第二個點的一階泰勒展開的等高線

alpha = 0。0312；

z = f（1。875， -1。125） - alpha * dot（［4；100］，［4；100］）；

gui（）

程式結果

上面的程式和結果是針對初始點和下一點之間的幾何關係影象，然後將函式在第二點處二階泰勒展開：

$f(x_{k+1})=f(x_k - \alpha d_k) \approx f(x_k) - \alpha \triangledown^T f(x_k) d_k+\frac{1}{2}\alpha^2 d_k^TGd_k$

，取

$d_k = \triangledown f(x_k)$

，代入第1點和第2點資料得：

$f(1.875, -1.125) = \\ f(2,2) - 0.0312*[4,100]*[4;100] +0.5*0.0312^2*[4,100]*[2,0;0,50]*[4;100]\\ \Downarrow \\ 35.1562 \approx 104 -312.4992+243.376$

這裡就很完美的解釋了4。4。1中的悖論，即使用一階泰勒近似時結論有問題。從上面的資料可以看出問題出在一階泰勒的誤差是非常大的，也說明了

$\alpha d_k$

的值不足以小到可以使用一階泰勒近似。

4。4。4 解釋收斂性

透過AG和全松馳的收斂性比較，可以看出本例中AG的收斂速度是比全松馳慢的，猜測這個可能跟終止準則的設定不合理有關。

4。5 牛頓法及其變種

4。5。1 理論推導

牛頓法：在

處將

二階泰勒展開，如下所示：

$f(x) \approx f(x_k) + \triangledown^Tf(x-x_k)(x-x_k)+\frac{1}{2}(x-x_k)^T H (x-x_k)$

然後將上式對

求導並取0，得到下一個近似的極小值點

$x_{k+1}$

。

注意：上式也可以理解為搜尋方向

$d_k=-H^{-1} \triangledown f(x_k)$

步長

$\alpha_k=1$

的搜尋，但是上式並不未包含沿下降方向搜尋的含義，只是求了個極值點。除此之外，搜尋方向還需要計算二階海森矩陣的逆，可能無法實現。

阻尼牛頓法：定義

$d_k=-H^{-1} \triangledown f(x_k)$

為牛頓方向，也是此演算法的搜尋方向，改變的是步長，定義阻尼因子

$\alpha_k$

為沿牛頓方向一維搜尋的最佳步長，即：

$f(x_{k+1})=f(x_k+\alpha_k d_k)=min_\alpha f(x_k+\alpha d_k)$

擬牛頓法：針對牛頓法中需要計算二階海森矩陣及其逆的問題，使用一個對稱正定的矩陣來代替H進行計算，並在計算中更新此近似矩陣。

4。5。2 程式（阻尼牛頓法）

#= program3。jl

使用阻尼牛頓法求函式 f（x）=x1^2+25x2^2的極小點。 =#

using

LinearAlgebra

；

# 將 f（x）轉換為矩陣形式，也可以不用轉，但求梯度方便。

［

；

］；

gradient_fx

（

）

；

（

）

0。5

‘

；

（

，

）

；

x_next

［

；

］；

x_list

copy

（

x_next

）；

y_list

［

（

x_next

）］；

grad_list

［

0。0

；

0。0

］；

alpha_list

［

0。0

］；

TOL

0。0001

；

while

norm

（

x_next

）

TOL

global

x_next

，

x_list

，

y_list

，

grad_list

，

alpha_list

；

# 計算當前點梯度

gradient_k

gradient_fx

（

x_next

）；

# 設定搜尋方向為牛頓方向

d_k

inv

（

）

gradient_k

；

# 計算步長為最佳步長

alpha

（

x_next

’

［

；

］）

（

gradient_k

‘

［

；

］）；

# 計算下一點

x_next

alpha

d_k

；

# 將運算過程中的值儲存

x_list

hcat

（

x_list

，

x_next

）；

y_list

push！

（

y_list

，

（

x_next

））；

grad_list

hcat

（

grad_list

，

gradient_k

）；

alpha_list

push！

（

alpha_list

，

alpha

）；

end

data

hcat

（

hcat

（

hcat

（

x_list

’

，

y_list

），

alpha_list

），

grad_list

‘

）；

data

。

round

（

data

，

digits

）；

println

（

“k

f（x）

alpha

Gradient_x1

Gradient_x2”

）；

row_len

，

col_len

size

（

data

）；

for

row

：

row_len

（

row

，

“

”

）

for

col

：

col_len

（

data

［

row

，

col

］，

“

”

）；

end

println

（

“”

）；

end

從結果上看，阻尼牛頓法是真的666，精確搜尋，一步到位。從結果中也能看出這個例子比較特別，正好阻尼牛頓和牛頓法是一致的。猜測應該是對於二次函式，這兩種方法是一致的。

4。5。3 擬牛頓法（使用變尺度矩陣實現）

尺度矩陣－以二次型函式為例

使用

$x \leftarrow Qx$

代入標準二次型函式中，二次項變為：

$\frac{1}{2}x^TGx \leftarrow \frac{1}{2}x^TQ^TGQx$

，目的是改變二次項的偏心程度，使函式偏心率減小，函式的等值面變為球面，這樣任意點的梯度都指定球心（即最小點）。通俗的說，將一個橢圓形的等高線轉變為同心圓，這樣每點的負梯度方向指定球心。

若G是對稱正定矩陣，則透過LDLt分解（或相似對角化或相似正交對角化）可使

，這樣變換後的二次項為單位矩陣，即偏心率為0。再對上式進行變換可得：

$Q^TG=Q^{-1} \Rightarrow QQ^TG=I \Rightarrow H=QQ^T=G^{-1}$

，這樣可能透過Q求G的逆矩陣。當用H來代替代替牛頓法中的二階海森矩陣的逆裡，即為擬牛頓法的一種實現方法。

變尺度矩陣的建立－對於一般函式

採用牛頓法，更新方程為

$x_{k+1} = x_k -\alpha_k G_k^{-1} \triangledown f(x_k)$

，其中

$\alpha_k$

採用一維搜尋的最佳步長（即函式值下降最大）。

構造一個矩陣序列

$\{ H_k \}$

來逼近海森矩陣的逆序列

$\{ G^{-1} \}$

，每迭代一次，尺度矩陣就改變一次，稱之為變尺度。為了使

滿足近似且容易計算，需要加上以下附加條件：

$\{ H_k \}$

中每一個元素都是對稱正定，這樣才能滿足搜尋方向

$d_k=-H_k \triangledown f(x_k)$

為下降方向

迭代簡單。一個簡單的校正公式為

$H_{k+1} = H_K +E_k$

，其中

為校正矩陣

$\{ H_k \}$

滿足擬牛頓條件。直接給結論，

$H_{k+1}(g_{k+1}-g_k) = x_{k+1} - x_k \\ \Downarrow \\ y_k=(g_{k+1}-g_k) \quad s_k=x_{k+1} - x_k \\ \Downarrow \\ H_{k+1} y_k = s_k$

，其中

表示第 k 次迭代處的梯度。推導方法跟共軛方向法有相通的地方，對於二次函式，能產生對G共軛的搜尋方向，也可以看作一種共軛方向法。

條件2和3看似是衝突的，因為在每次的迭代中

是可能求出的，然後僅按照這兩個向量是無法辨識出唯一的

$H_{k+1}$

，所以需要使用條件2 進行迭代。

條件2的實現需要

的確定，

$(H_k +E_k)y_k = s_k \Rightarrow E_k y_k = s_k -H_k y_k$

即是

需要滿足的條件。按照這一條件構造不同的

，就形成了不同的變尺度法，DFP就是其中的一種。直接給出DFP的校正公式如下：

$H_{k+1} = H_k + \frac{s_k s_k^T}{s_k^T y_k} - \frac{H_k y_k y_k^T H_k}{y_k^T H_k y_k}, k=0, 1, 2, ...$

4。5。4 程式框圖及程式碼－變尺度法

變尺度法的程式框圖

程式（使用DFP的變尺度法）：

#= program4。jl

使用擬牛頓法（基於DFP的變尺度法）求函式 f（X）=0。5*X^T*G*X+b^TX的極小點，

其中G=［2 -2； -2 4］， b=［-4；0］。 =#

using

LinearAlgebra

，

Plots

；

pyplot

（）

# 定義函式

［

；

］；

［

；

］；

（

）

0。5

’

‘

；

（

，

）

；

# 求梯度

g_fx

（

）

；

# 設定初始引數

x_k

［

1。0

；

］；

H_k

［

1。0

；

1。0

］；

# 設定誤差

TOL

0。0001

；

# 設定最大迴圈次數

N_MAX

；

# 設定布林量控制迴圈

gate

true

# 設定陣列來儲存運算過程中的結果

x_list

［］；

y_list

［］；

g_list

［］；

alpha_list

［］；

H_list

［］；

while

gate

global

x_k

，

H_k

，

x_list

，

y_list

，

g_list

，

alpha_list

，

H_list

，

gate

# 儲存當前的 x 向量

push！

（

x_list

，

x_k

）；

# 儲存當前的 y 函式值

push！

（

y_list

，

（

x_k

））；

# 計算當前點梯度

g_k

g_fx

（

x_k

）；

push！

（

g_list

，

g_k

）；

# 如果梯度為0，則退出

norm

（

g_k

）

TOL

gate

false

；

end

# 設定變尺度矩陣來近似海森矩陣的逆

length

（

H_list

）

y_k

g_k

g_list

［

end

］；

s_k

x_k

x_list

［

end

］；

# 透過DFP演算法求Ek

E_k

（

s_k

’

）

（

s_k

‘*

y_k

）

（

H_k

y_k

’*

H_k

）

（

y_k

‘*

H_k

y_k

）；

# 更新Hk

H_k

E_k

；

push！

（

H_list

，

H_k

）；

else

push！

（

H_list

，

H_k

）；

end

# 設定搜尋方向

d_k

H_k

g_k

；

# 計算步長為最佳步長

alpha

（

d_k

’*

（

x_k

））

（

d_k

‘*

d_k

）；

push！

（

alpha_list

，

alpha

）；

# 計算下一點

x_k

alpha

d_k

；

# 判斷是否達到精度要求

norm

（

x_k

x_list

［

end

］）

TOL

length

（

x_list

）

N_MAX

gate

false

；

end

# 將資料整合在一起

row_len

length

（

x_list

）；

col_len

sum

（

map

（

length

，

［

x_list

［

］，

g_list

［

］，

H_list

［

］］））

；

data

zeros

（

Float64

，

（

row_len

，

col_len

））；

for

：

length

（

x_list

）

row_x

，

row_g

，

row_H

map

（

vec

，

［

x_list

［

］，

g_list

［

］，

H_list

［

］］）；

row

hcat

（

hcat

（

hcat

（

hcat

（

row_x

’

，

y_list

［

］），

alpha_list

［

］），

row_g

‘

），

row_H

’

）；

data

［

，

：

］

row

；

end

# 為了顯示美觀，將資料保留4位小數

data

。

round

（

data

，

digits

）；

println

（

“k

f（x）

alpha

Gradient_x1

Gradient_x2

H-Matrix（2x2）”

）；

for

row

：

row_len

（

row

，

“

”

）

for

col

：

col_len

（

data

［

row

，

col

］，

“

”

）；

end

println

（

“”

）；

end

lims_x1

：

0。01

：

；

lims_x2

：

0。01

：

；

contour

（

lims_x1

，

lims_x2

，

fill

true

）；

display

（

）

png

（

“contour”

）；

結果：

原函式的等高線圖

五、無約束最佳化－共軛梯度法

5。1、共軛

方向

法原理（不是共軛梯度）

以二次函式

$f(x)=\frac{1}{2}x^TGx+b^Tx+c$

為例。從初始點

，沿搜尋方向

行走步長

$\alpha_0$

到達

點，其中

$\alpha_0$

是此方向上的最佳步長，即在

點處的負梯度方向與

垂直，如下圖所示：

負梯度方向與共軛方向

如上圖所示，在選擇

方向時選擇一種比較貪婪的方法，希望直接指定最優點

，那麼該怎麼確定

？最優點

處的的導數為0（也等價於梯度為0），所以：

$\triangledown f(x_*) = Gx_* + b = 0 \Rightarrow \\ \triangledown f(x_1+\alpha_1 d^1) = G(x_1+\alpha_1 d^1) +b \Rightarrow \\ (Gx_1+b) + \alpha_1Gd_1 = \triangledown f(x_1)+ \alpha_1Gd_1=0 \Rightarrow \\ d_0^T \triangledown f(x_1) + \alpha_1 d_0^TGd_1=0 \Rightarrow \\ \alpha_1 d_0^TGd_1=0$

注意：

$d_0^T \triangledown f(x_1)=0$

是由上步中搜索步長最佳得到的，

$\alpha_1$

不能為0，否則當前點即為最優點。因此由最後一個等式可得：

當相鄰兩個搜尋方向關於G共軛的時候，那麼每步的搜尋方向均指向最優點

。稱這樣的向量為關於G共軛或共軛方向。當G是單位矩陣歸，共軛退化為正交。

共軛向量組有一些性質：

非m個非零向量對一個對稱正定矩陣G共軛，則這m個向量是線性無關的。反之，可推匯出若m個向量線性相關，則必不相對於G共軛。線性無關是一個必要條件。

在n維空間中互相共軛的非零向量的個數不超過n。

從任意初始點出發，順次沿n個G的共軛方向進行一維搜尋，最多經過n次迭代就可以找到二次函式的極小點。

注：共軛方向與牛頓方向推導區別，牛頓方向是在迭代點附近二階泰勒展開求導找最優值，共軛方向是整個函式求導尋最優解。

5。2、共軛方向的構造

由上述可知，演算法的難點在於共軛向量組的構建。對應於不同的構造方法，形成不同的共軛方向法，如共軛梯度法，鮑威爾法等。施密特正交法的推廣也可以從n個線性無關的得到共軛向量組。由於篇幅，不再細述，參考6中有解釋。

5。3、共軛梯度法

共軛梯度法是共軛方向法中的一種，因為每一個共軛向量依賴於迭代點處的負梯度而構造出來。推導如下：還以上文中提到的二次函式為例

$x_{k+1}-x_k=\alpha_k d_k \Rightarrow\\ g_k = Gx^k+b, \quad g_{k+1} = Gx^{k+1}+b \Rightarrow \\ g_{k+1}-g_k=G(x_{k+1}-x_k)=\alpha_kGd_k$

若已知某一方向向量

，即成共軛方向，則：

$d_j^T(g_{k+1}-g_k)=0$

，說明

的一個共軛方向與相鄰點的梯度差正交，不用計算G即可算出共軛方向。

共軛梯度法的幾何說明

由於推導複雜，直接給出結論：共軛方向的遞推公式為：

其中，

是第

步的梯度。

5。4、共軛梯度法程式框圖

5。5、程式（共軛梯度法）

#= program5。jl

使用共軛梯度法求函式 f（X）=0。5*X^T*G*X+b^TX的極小點，

其中G=［2 -2； -2 4］， b=［-4；0］。 =#

using LinearAlgebra；

# 定義函式

G = ［2 -2； -2 4］；

b = ［-4；0］；

f（X） = 0。5 * X‘ * G * X + b’ * X；

# 求梯度

g_fx（X） = G * X + b；

# 設定初始引數

x_k = ［1。0；1］；

# 設定誤差

TOL = 0。0001；

# 設定最大迴圈次數

N_MAX = 30；

# 設定布林量控制迴圈

gate = true

# 設定陣列來儲存運算過程中的結果

x_list = ［］；

y_list = ［］；

g_list = ［］；

alpha_list = ［］；

conj_list = ［］；

while gate

global x_k， conj_k， x_list， y_list， g_list， alpha_list， conj_list， gate

# 儲存當前的 x 向量

push！（x_list， x_k）；

# 儲存當前的 y 函式值

push！（y_list， f（x_k））；

# 計算當前點梯度

g_k = g_fx（x_k）；

push！（g_list， g_k）；

# 如果梯度為0，則退出

if norm（g_k） < TOL

gate = false；

end

# 求共軛方向

if length（conj_list） > 0

beta_k = norm（g_k）^2 / norm（g_list［end - 1］）^2；

conj_k = -g_k + beta_k * conj_k；

push！（conj_list， conj_k）；

else

conj_k = -g_k；

push！（conj_list， conj_k）；

end

# 設定搜尋方向

d_k = conj_k；

# 計算步長為最佳步長

alpha = -（d_k‘ * （G * x_k + b）） / （d_k’ * G * d_k）；

push！（alpha_list， alpha）；

# 計算下一點

x_k = x_k + alpha * d_k；

# 判斷是否達到精度要求

if norm（x_k - x_list［end］） < TOL || length（x_list） > N_MAX

gate = false；

end

# 將資料整合在一起

row_len = length（x_list）；

col_len = sum（map（length，［x_list［1］， g_list［1］， conj_list［1］］）） + 2；

data = zeros（Float64，（row_len， col_len））；

for i in 1：length（x_list）

row_x， row_g， row_H = map（vec，［x_list［i］， g_list［i］， conj_list［i］］）；

row = hcat（hcat（hcat（hcat（row_x‘， y_list［i］）， alpha_list［i］）， row_g’）， row_H‘）；

data［i，：］ = row；

end

# 為了顯示美觀，將資料保留4位小數

data = @。 round（data， digits = 4）；

println（“k \t x1 \t x2 \t f（x） \t alpha \t Gradient_x1 \t Gradient_x2 \t Conj-Vector（2）”）；

for row in 1：row_len

print（row， “\t”）

for col in 1：col_len

print（data［row， col］， “\t”）；

end

println（“”）；

end

結果：

六、座標輪換法及鮑威爾方法

6。1、座標輪換法

每次搜尋只允許一個變數變化，其餘變數保持不變，即沿座標方向輪流進行一維搜尋最佳步長的尋優方法。把多變數的最佳化問題輪流地轉化為單變數的最佳化問題。計算過程中不需要目標函式的導數，相對簡單。缺點是要經過多次曲折迂迴的路徑才能達到極值點，尤其在極值點附近步長很小，收斂很慢。

座標輪換法的搜尋過程

6。2、鮑威爾方法

基本思想是在不用導數的前提下，直接利用函式值來迭代構造共軛方向的一種共軛方向法

6。2。1 共軛方向的產生

透過一維方向確定共軛方向

上圖中，

$x_k 和 x_{k+1}$

為從不同點出發，沿同一方向

（平行向量也是同一方向）進行一維搜尋最佳步長得到的點，顯然由之前的推導可得：

$(d^j)^Tg_k = (d^j)^Tg_{k+1}=0$

，又兩點的梯度表示為：

$g_k = Gx_k+b \quad g_{k+1}=Gx_{k+1}+b \Rightarrow g_{k+1}-g_k = G(x_{k+1}-x_k)$

，因此有：

$(d^j)^T(g_{k+1}-g_k)=(d^j)^TGd_k=0,\quad d^k=x_{k+1}-x_k$

，說明

關於G共軛。

6。2。2 鮑威爾演算法

二維鮑威爾法的幾何表示

按照座標輪換的思路，首先選擇各座標軸上的單位向量輪流搜尋最佳步長得到

點，然後使用起始點

與此點的連續方向取代第一個座標方向

組成新的線性無關向量組（即

）作為下一輪的搜尋方向得到點

，由6。2。1中的結論可知

的連線方向

與

是共軛方向，然後進行新的搜尋。

6。2。3 改進鮑威爾

由於鮑威爾演算法中的新生成的向量組可能會線性相關，因上，在改進演算法中需要判斷向量組中的哪個向量最壞，然後用新產生的向量替換這個最換的向量，以保證逐次生成共軛方向。具體不再細述。

七、單形替換法

上述的方法都使用了函式的導數，在不計算導數的情況下，先演算法若干個點的函式值，從它們的大小關係中看出函式的變化趨勢，尋找下降方向，這就是單形替換法的思路。

單純形指的是n維空間中具有n+1個頂點的多面體。利用單純形的頂點，將函式值進行大小排序，從而確定有利的搜尋方向和步長，找到一個好的點替換單純形中較差的點，使單純形不斷地向目標函式的極小值靠近，直到搜尋到極小點為止，這就是單形替換法的基本思路。

單形替換法

八、總結本文中約束最佳化方法搜尋方向的關係

無約束最佳化方法搜尋方向之間的相互聯絡

後記：

最後的鮑威法和單形變換法沒寫太多，一方面是因為要轉向約束最佳化的內容，另一方面是貌似知乎出了一個bug，文章圖片太多了瀏覽器總崩潰。

上述最佳化方法各有各的優點和缺點，應該是針對不同的問題使用不同的方法（這個中國人都知道@-@），具體來說，比如當高維最佳化中目標函式的二維海森矩陣可不是好求的，更不用說其逆的運算，因此牛頓法就不是很適用了。當一階梯度求解起來麻煩或難以求解時，就要使用零階方法了。綜合來講，覺得共軛梯度法比較好些。這些都是搜尋方向的選擇。對於步長來說，看起來用全松馳（即最佳步長）的求解比較好，應該共軛方向的推導過程中嚴重依賴最佳步長時的梯度與搜尋方向垂直這個條件，Armijio-Goldstein應該也有自己擅長的地方。歡迎讀者在這個地方展開激烈討論。

除此之外，還有一個比較大的問題就是，上述演算法求解的都是區域性最優解，如何保證區域性最佳是全域性最優，這就需要凸最佳化上場了。

特此宣告：我在知乎上的所有文章、言論、見解等只代表我個人的觀念，與任何研究機構、商業公司等無關。所有文章僅供個人學習交流之用，側耳傾聽任何與文章相關的評論。在文章開頭會註明相關引用內容，所引用內容版權歸原作者所有，在此一併表示感謝。若未註明之處，由此可能造成的侵權影響請聯絡本人進行刪除或者修改。若個別讀者不遵守此宣告進行轉載，並因此產生商業行為，並造成對文章內引用內容的侵權行為與本人無關。如果文章能對個人學習、工作帶來一些思路，那麼文章的目的已經達到。如果引用文章相關內容，請註明出處。

小蜜蜂問答

小蜜蜂問答

無約束最佳化方法及程式碼實現

推薦文章

小蜜蜂問答

小蜜蜂問答

無約束最佳化方法及程式碼實現

相關文章

p!=NULL與p？

C#能不能直接把一個LIST1直接複製給另外一個LIST2？

人工神經網路學習筆記2——MATLAB神經網路工具箱

Linux核心中常用的資料結構和演算法

推薦文章