參考資料

《精通MATLAB最最佳化計算（第二版）》

程式設計工具

Matlab 2019a

目錄

石中居士：最最佳化計算與Matlab實現——目錄

擬牛頓法

原理：

牛頓法的收斂速度雖然較快，但要求海森矩陣可逆，要計算二階導數和逆矩陣，就加大了計算機計算量和儲存量。為克服牛頓法的缺點，同時保持較快收斂速度的優點，就產生了擬牛頓法。

擬牛頓法是牛頓法的直接推廣，透過在試探點附近的二次逼近引進擬牛頓條件來確定線搜尋方向，它主要有DFP和BFGS兩種形式。

一般演算法步驟：

【1】給定初始點

x^{(0)}

，初始對稱正定矩陣

\boldsymbol{H}_{0}

，

g_{0}=g\left(x^{(0)}\right)

及精度

\varepsilon>0

；

【2】計算搜尋方向

p^{(k)}=-\boldsymbol{H}_{k}{g}_{k}

；

【3】作直線搜尋

x^{(k+1)}=F\left(x^{(k)}, p^{(k)}\right)

，計算

f_{k+1}=f\left(x^{(k+1)}\right)

，

g_{k+1}=g\left(x^{(k+1)}\right)

，

S_{k}=x^{(k+1)}-x^{(k)}

，

y_{k}=g_{k+1}-g_{k}

【4】判斷終止準則是否滿足；

【5】令

\boldsymbol{H}_{k+1}=\boldsymbol{H}_{k}+\boldsymbol{E}_{k}

置

k=k+1

，轉【2】。

DFP法

原理

DFP演算法中的校正公式為：

\boldsymbol{H}_{k+1}=\boldsymbol{H}_{k}+\frac{{S}_{k} {S}_{k}^{T}}{{S}_{k}^{T} {y}_{k}}-\frac{\boldsymbol{H}_{k} {y}_{k} {y}_{k}^{T} \boldsymbol{H}_{k}}{{y}_{k}^{T} \boldsymbol{H}_{k} {y}_{k}}

為了保證

\boldsymbol{H}_{k}

的正定性，在下面的演算法步驟中迭代一定次數後，重置初始點和迭代矩陣再進行迭代。

演算法步驟

用DFP法求無約束多維極值問題

\min f(x), x \in \mathbf{R}^{n}

的演算法步驟如下：

【1】給定初始點

x^{(0)}

，初始矩陣

\boldsymbol{H}_{0}={I}_{n}

及精度

\varepsilon>0

；

【2】若

\left\|\nabla f\left(x^{(0)}\right)\right\| \leqslant \varepsilon

，停止，極小點為

x^{(0)}

，否則轉【3】；

【3】取

{p}^{(0)}=-\boldsymbol{H}_{0} \nabla f\left(x^{(0)}\right)

，且令

k=0

；

【4】一維搜尋法求

t_{k}

，使得：

f\left(x^{(k)}+t_{k} p^{(k)}\right)=\min _{r \geq 0} f\left(x^{(k)}+t p^{(k)}\right)

令

x^{(k+1)}=x^{(k)}+t_{k} p^{(k)}

，轉【5】；

【5】若

\left\|\nabla f\left(x^{(k+1)}\right)\right\| \leqslant \varepsilon

，停止，極小點為

x^{(0)}

，否則轉【3】；

【6】若

k+1=n

，令

x^{(0)}=x^{(n)}

，轉【3】，否則轉【7】；

【7】令

\boldsymbol{H}_{k+1}=\boldsymbol{H}_{k}+\frac{\left(x^{(k+1)}-x^{(k)}\right)\left(x^{(k+1)}-x^{(k)}\right)^{T}}{\left(x^{(k+1)}-x^{(k)}\right)^{T}\left(\nabla f\left(x^{(k+1)}\right)-\nabla f\left(x^{(k)}\right)\right)}

-\frac{\boldsymbol{H}_{k}\left(\nabla f\left(x^{(k+1)}\right)-\nabla f\left(x^{(k)}\right)\right)\left(\nabla f\left(x^{(k+1)}\right)-\nabla f\left(x^{(k)}\right)\right)^{T} \boldsymbol{H}_{k}}{\left(\nabla f\left(x^{(k+1)}\right)-\nabla f\left(x^{(k)}\right)\right)^{T} \boldsymbol{H}_{k}\left(\nabla f\left(x^{(k+1)}\right)-\nabla f\left(x^{(k)}\right)\right)}

取

p^{(k)}=-\boldsymbol{H}_{k+1} \nabla f\left(x^{(k+1)}\right)

，置

k=k+1

，轉【4】

Matlab程式碼與試算

用DFP法求函式

f(t, s)=t^{2}-t s+s^{2}+2

的極小值，其中初始點取

x^{0}=(-4,6)

。

test.m

syms

t

s

；

f

=

t

^

2

-

t

*

s

+

s

^

2

+

2

；

［

x_optimization

，

f_optimization

］

=

Quasi_Newton_DFP_Method

（

f

，［

-

4

6

］，［

t

s

］）；

x_optimization

=

double

（

x_optimization

）；

f_optimization

=

double

（

f_optimization

）；

x_optimization

f_optimization

Quasi_Newton_DFP_Method.m

function

［x_optimization，f_optimization］

=

Quasi_Newton_DFP_Method

（

f，x0，var_x，epsilon

）

format

long

；

% f：目標函式

% x0：初始點

% var_x：自變數向量

% epsilon：精度

% x_optimization：目標函式取最小值時的自變數值

% f_optimization：目標函式的最小值

if

nargin

==

3

epsilon

=

1。0e-6

；

end

x0

=

transpose

（

x0

）；

var_x

=

transpose

（

var_x

）；

n

=

length

（

var_x

）；

syms

t

；

H0

=

eye

（

n

，

n

）；

grad_fx

=

jacobian

（

f

，

var_x

）；

grad_fx0

=

subs

（

grad_fx

，

var_x

，

x0

）；

p0

=

-

H0

*

transpose

（

grad_fx0

）；

k

=

0

；

xk

=

x0

；

pk

=

p0

；

Hk

=

H0

；

while

1

grad_fxk

=

subs

（

grad_fx

，

var_x

，

xk

）；

if

norm

（

grad_fxk

）

<

=

epsilon

x_optimization

=

x0

；

break

；

end

yk

=

xk

+

t

*

pk

；

fyk

=

subs

（

f

，

var_x

，

yk

）；

［

xmin

，

xmax

］

=

Advance_and_Retreat_Method

（

fyk

，

0

，

0。00001

）；

tk

=

Golden_Selection_Method

（

fyk

，

xmin

，

xmax

）；

xk_next

=

xk

+

tk

*

pk

；

grad_fxk_next

=

subs

（

grad_fx

，

var_x

，

xk_next

）；

if

norm

（

grad_fxk_next

）

<

=

epsilon

x_optimization

=

xk_next

；

break

；

end

if

k

+

1

==

n

xk

=

xk_next

；

continue

；

else

dx

=

xk_next

-

xk

；

dgrad_fx

=

grad_fxk_next

-

grad_fxk

；

dgrad_fx

=

transpose

（

dgrad_fx

）；

dx_T

=

transpose

（

dx

）；

dgrad_fx_T

=

transpose

（

dgrad_fx

）；

m_dx

=

dx

*

dx_T

；

m_dgrad_fx

=

dgrad_fx

*

dgrad_fx_T

；

Hk_next

=

Hk

+

m_dx

/

（

dx_T

*

dgrad_fx

）

-

（

Hk

*

m_dgrad_fx

*

Hk

）

/

（

dgrad_fx_T

*

Hk

*

dgrad_fx

）；

pk

=

-

Hk_next

*

transpose

（

grad_fxk_next

）；

k

=

k

+

1

；

xk

=

xk_next

；

Hk

=

Hk_next

；

end

end

f_optimization

=

subs

（

f

，

var_x

，

x_optimization

）；

format

short

；

命令列視窗

x_optimization =

1。0e-07 *

0。9065

0。6844

f_optimization =

2。0000

成功求得函式

f(t, s)=t^{2}-t s+s^{2}+2

的極小點。

BFGS法

原理

BFGS演算法中的校正公式為

\boldsymbol{H}_{k+1}=\boldsymbol{H}_{k}+\frac{S^{(k)}\left(S^{(k)}\right)^{T}}{\left(S^{(k)}\right)^{T} y^{(k)}}\left[1+\frac{\left(y^{(k)}\right)^{T} \boldsymbol{H}_{k} y^{(k)}}{\left(S^{(k)}\right)^{T} y^{(k)}}\right]

-\frac{1}{\left(S^{(k)}\right)^{T} y^{(k)}}\left[S^{(k)}\left(y^{(k)}\right)^{T} \boldsymbol{H}_{k}+\boldsymbol{H}_{k} y^{(k)}\left(S^{(k)}\right)^{T}\right]

為了保證

\boldsymbol{H}_{k}

的正定性，在下面的演算法步驟中迭代一定次數後，重置初始點和迭代矩陣再進行迭代。

演算法步驟

用BFGS法求無約束多維極值問題

\min f(x), x \in \mathrm{R}^{n}

的演算法步驟如下：

【1】給定初始點

x^{(0)}

，初始矩陣

\boldsymbol{H}_{0}={I}_{n}

及精度

\varepsilon>0

；

【2】若

\left\|\nabla f\left(x^{(0)}\right)\right\| \leqslant \varepsilon

，停止，極小點為

x^{(0)}

，否則轉【3】；

【3】取

{p}^{(0)}=-\boldsymbol{H}_{0} \nabla f\left(x^{(0)}\right)

，且令

k=0

；

【4】一維搜尋法求

t_{k}

，使得：

f\left(x^{(k)}+t_{k} p^{(k)}\right)=\min _{r \geq 0} f\left(x^{(k)}+t p^{(k)}\right)

令

x^{(k+1)}=x^{(k)}+t_{k} p^{(k)}

，轉【5】；

【5】若

\left\|\nabla f\left(x^{(k+1)}\right)\right\| \leqslant \varepsilon

，停止，極小點為

x^{(0)}

，否則轉【3】；

【6】若

k+1=n

，令

x^{(0)}=x^{(n)}

，轉【3】，否則轉【7】；

【7】令

\boldsymbol{H}_{k+1}=\boldsymbol{H}_{k}+\frac{S^{(k)}\left(S^{(k)}\right)^{T}}{\left(S^{(k)}\right)^{T} y^{(k)}}\left[1+\frac{\left(y^{(k)}\right)^{T} \boldsymbol{H}_{k} y^{(k)}}{\left(S^{(k)}\right)^{T} y^{(k)}}\right]

-\frac{1}{\left(S^{(k)}\right)^{T} y^{(k)}}\left[S^{(k)}\left(y^{(k)}\right)^{T} \boldsymbol{H}_{k}+\boldsymbol{H}_{k} y^{(k)}\left(S^{(k)}\right)^{T}\right]

其中

\begin{aligned} &S^{(k)}=x^{(k+1)}-x^{(k)}\\ &y^{(k)}=\nabla f\left(x^{(k+1)}\right)-\nabla f\left(x^{(k)}\right) \end{aligned}

，取

p^{(k)}=-\boldsymbol{H}_{k+1} \nabla f\left(x^{(k+1)}\right)

，置

k=k+1

，轉【4】。

Matlab程式碼與試算

用BFGS法求函式

f(t, s)=t^{2}-t s+s^{2}+2

的極小值，初始點取

x^{0}=(-1,1)

。

test.m

syms

t

s

；

f

=

t

^

2

-

t

*

s

+

s

^

2

+

2

；

［

x_optimization

，

f_optimization

］

=

Quasi_Newton_BFGS_Method

（

f

，［

-

1

1

］，［

t

s

］）；

x_optimization

=

double

（

x_optimization

）；

f_optimization

=

double

（

f_optimization

）；

x_optimization

f_optimization

Quasi_Newton_BFGS_Method.m

function

［x_optimization，f_optimization］

=

Quasi_Newton_BFGS_Method

（

f，x0，var_x，epsilon

）

format

long

；

% f：目標函式

% x0：初始點

% var_x：自變數向量

% epsilon：精度

% x_optimization：目標函式取最小值時的自變數值

% f_optimization：目標函式的最小值

if

nargin

==

3

epsilon

=

1。0e-6

；

end

x0

=

transpose

（

x0

）；

var_x

=

transpose

（

var_x

）；

n

=

length

（

var_x

）；

syms

t

；

H0

=

eye

（

n

，

n

）；

grad_fx

=

jacobian

（

f

，

var_x

）；

grad_fx0

=

subs

（

grad_fx

，

var_x

，

x0

）；

p0

=

-

H0

*

transpose

（

grad_fx0

）；

k

=

0

；

xk

=

x0

；

pk

=

p0

；

Hk

=

H0

；

while

1

grad_fxk

=

subs

（

grad_fx

，

var_x

，

xk

）；

if

norm

（

grad_fxk

）

<

=

epsilon

x_optimization

=

x0

；

break

；

end

f_xk_tkpk

=

subs

（

f

，

var_x

，（

xk

+

t

*

pk

））；

［

xmin

，

xmax

］

=

Advance_and_Retreat_Method

（

f_xk_tkpk

，

0

，

0。00001

）；

tk

=

Golden_Selection_Method

（

f_xk_tkpk

，

xmin

，

xmax

）；

xk_next

=

xk

+

tk

*

pk

；

grad_fxk_next

=

subs

（

grad_fx

，

var_x

，

xk_next

）；

if

norm

（

grad_fxk_next

）

<

=

epsilon

x_optimization

=

xk_next

；

break

；

end

if

k

+

1

==

n

xk

=

xk_next

；

continue

；

else

Sk

=

xk_next

-

xk

；

yk

=

grad_fxk_next

-

grad_fxk

；

yk

=

transpose

（

yk

）；

Sk_T

=

transpose

（

Sk

）；

yk_T

=

transpose

（

yk

）；

m_Sk

=

Sk

*

Sk_T

；

Hk_next

=

Hk

+

（

m_Sk

/

（

Sk_T

*

yk

））

*

（

1

+

（

yk_T

*

Hk

*

yk

）

/

（

Sk_T

*

yk

））

-

（

Sk

*

yk_T

*

Hk

+

Hk

*

yk

*

Sk_T

）

/

（

Sk_T

*

yk

）；

pk

=

-

Hk_next

*

transpose

（

grad_fxk_next

）；

k

=

k

+

1

；

xk

=

xk_next

；

Hk

=

Hk_next

；

end

end

f_optimization

=

subs

（

f

，

var_x

，

x_optimization

）；

format

short

；

命令列視窗

x_optimization =

1。0e-11 *

0。7324

-0。7324

f_optimization =

2

擬牛頓法的好處在於不用計算海森矩陣，大大減少了計算量。

最最佳化計算與matlab實現（10）——無約束多維極值問題——擬牛頓法

目錄

石中居士：最最佳化計算與Matlab實現——目錄