最速下降法的python实现——基于数值计算

最速下降法的python实现——基于数值计算

最速下降法的python实现——基于数值计算

• 1. 进退法——获取单峰区间
• 2. 在单峰区间进行一维搜索——黄金分割法(0.618法)
• 3. 最速下降法——求解最优值
• 4. 基于数值方法的实现

1. 进退法——获取单峰区间

∙ bullet ∙　进退法通过“试探”的方式来确定单峰区间，目的是为了使函数值呈现“高-低-高”的三点。

quad 计算步骤：
quad 步骤 ( 1 ) : (1): (1): 给定初始点 x ( 0 ) ∈ R x^{(0)}in R x(0)∈R，初始步长 h 0 > 0 h_0>0 h0​>0，令 h = h 0 h=h_0 h=h0​， x ( 1 ) = x ( 0 ) x^{(1)}=x^{(0)} x(1)=x(0)，计算 f ( x ( 1 ) ) f(x^{(1)}) f(x(1))，置 k = 0 k=0 k=0；
quad 步骤 ( 2 ) : (2): (2): 令 x ( 4 ) = x ( 1 ) + h x^{(4)}=x^{(1)}+h x(4)=x(1)+h，计算 f ( x ( 4 ) ) f(x^{(4)}) f(x(4))，置 k = k + 1 k=k+1 k=k+1；
quad 步骤 ( 3 ) : (3): (3): 如果 f ( x ( 4 ) ) < f ( x ( 1 ) ) f(x^{(4)}) quadqquadqquad 如果 f ( x ( 4 ) ) > f ( x ( 1 ) ) f(x^{(4)})>f(x^{(1)}) f(x(4))>f(x(1))，则转步骤 ( 5 ) (5) (5)；
quad 步骤 ( 4 ) : (4): (4): 令 x ( 2 ) = x ( 1 ) x^{(2)}=x^{(1)} x(2)=x(1)， x ( 1 ) = x ( 4 ) x^{(1)}=x^{(4)} x(1)=x(4)， f ( x ( 2 ) ) = f ( x ( 1 ) ) f(x^{(2)})=f(x^{(1)}) f(x(2))=f(x(1))， f ( x ( 1 ) ) = f ( x ( 4 ) ) f(x^{(1)})=f(x^{(4)}) f(x(1))=f(x(4))，置 h = 2 h h=2h h=2h，转步骤 ( 2 ) (2) (2)；
quad 步骤 ( 5 ) : (5): (5): 如果 k = 1 k=1 k=1，则转步骤 ( 6 ) (6) (6)；否则，转到步骤 ( 7 ) (7) (7)；
quad 步骤 ( 6 ) : (6): (6): 置 h = − h h=-h h=−h， x ( 2 ) = x ( 4 ) x^{(2)}=x^{(4)} x(2)=x(4)， f ( x ( 2 ) ) = f ( x ( 4 ) ) f(x^{(2)})=f(x^{(4)}) f(x(2))=f(x(4))，转步骤 ( 2 ) (2) (2)；
quad 步骤 ( 7 ) : (7): (7): 令 x ( 3 ) = x ( 2 ) x^{(3)}=x^{(2)} x(3)=x(2)， x ( 2 ) = x ( 1 ) x^{(2)}=x^{(1)} x(2)=x(1)， x ( 1 ) = x ( 4 ) x^{(1)}=x^{(4)} x(1)=x(4)，停止计算。

quad 可以得到包含最小值的区间 [ x ( 1 ) , x ( 3 ) ] [x^{(1)},x^{(3)}] [x(1),x(3)] 或者 [ x ( 3 ) , x ( 1 ) ] [x^{(3)},x^{(1)}] [x(3),x(1)]。
qquad
∙ bullet ∙　进退法两种典型的搜索示意图如下图所示：
qquad

左图： ( 2 ) → ( 3 ) → ( 4 ) ⏟ → ( 2 ) → ( 3 ) → ( 4 ) ⏟ → ( 2 ) → ( 3 ) → ( 5 ) ⏟ → ( 7 ) underbrace{(2)to(3)to(4)}tounderbrace{(2)to(3)to(4)}tounderbrace{(2)to(3)to(5)}totextcolor{red}{(7)} (2)→(3)→(4)​→ (2)→(3)→(4)​→ (2)→(3)→(5)​→(7)，单峰区间为 [ x 3 , x 1 ] = [ 0.1 , 2.5 ] [x_3,x_1]=[0.1,2.5] [x3​,x1​]=[0.1,2.5]
　　　　 f ( x 4 ) < f ( x 1 ) small f(x_4) f ( x 1 ) small f(x_4)>f(x_1)   f(x4​)>f(x1​)
　
右图： ( 2 ) → ( 3 ) → ( 5 ) → ( 6 ) ⏟ → ( 2 ) → ( 3 ) → ( 4 ) ⏟ → ( 2 ) → ( 3 ) → ( 5 ) ⏟ → ( 7 ) underbrace{(2)to(3)to(5)to(6)}tounderbrace{(2)to(3)to(4)}tounderbrace{(2)to(3)to(5)}totextcolor{red}{(7)} (2)→(3)→(5)→(6)​→ (2)→(3)→(4)​→ (2)→(3)→(5)​→(7)，单峰区间为 [ x 1 , x 3 ] = [ 0.1 , 1.3 ] [x_1,x_3]=[0.1,1.3] [x1​,x3​]=[0.1,1.3]
　　　　　　 f ( x 4 ) > f ( x 1 ) small f(x_4)>f(x_1) f(x4​)>f(x1​)　　　　　 f ( x 4 ) < f ( x 1 ) small f(x_4) f ( x 1 ) small f(x_4)>f(x_1) f(x4​)>f(x1​)
　　　　　掉头搜索　 h = − h h=-h h=−h

qquad
∙ bullet ∙　进退法两种典型的搜索示意图的代码

import numpy as np
def getInterval(f,x0,h):
    k = 0
    x1 = x0 
    f1 = f(x1)     
    while (1):
        x4 = x1 + h
        f4 = f(x4)
        k += 1
        if (k == 1) and (f4 > f1):
            x2 = x4
            h = -h 
            continue        
        if f4 < f1:
            x2 = x1
            x1 = x4
            h = 2*h
        if f4 > f1:
            x3 = x2
            x2 = x1
            x1 = x4            
            break           
    return np.fmin(x1,x3), np.fmax(x1,x3)    
if __name__ == "__main__": 
    f = lambda x: 2*(x-1)**2+1    
    #a,b = getInterval(f,-0.3,0.4)    	# 左图 [0.1,2.5]
    a,b = getInterval(f,1.3,0.4)		# 右图 [0.1,1.3]
    print('[a,b] = [%.1f,%.1f]'%(np.round(a,2),np.round(b,2)))      

2. 在单峰区间进行一维搜索——黄金分割法(0.618法)

∙ bullet ∙　 0.618 0.618 0.618 法的计算步骤

qquad 假设初始搜索区间为 [ a 1 , b 1 ] [a_1,b_1] [a1​,b1​]，第 1 1 1 次迭代需要取两个试探点：

{ λ 1 = a 1 + 0.382 ( b 1 − a 1 ) μ 1 = a 1 + 0.618 ( b 1 − a 1 ) qquadqquadbegin{cases}lambda_1=a_1+0.382(b_1-a_1) \ mu_1=a_1+0.618(b_1-a_1)end{cases} {λ1​=a1​+0.382(b1​−a1​)μ1​=a1​+0.618(b1​−a1​)​

( 1 ) qquad(1) (1) 置 k = 1 k=1 k=1，精度要求为 ε varepsilon ε

( 2 ) qquad(2) (2) 若 b k − a k < ε b_k-a_k

qquadqquad 若 f ( λ k ) > f ( μ k ) f(lambda_k)>f(mu_k) f(λk​)>f(μk​)，那么 { a k + 1 = λ k , b k + 1 = b k λ k + 1 = μ k μ k + 1 = a k + 1 + 0.618 ( b k + 1 − a k + 1 ) begin{cases}a_{k+1}=lambda_k,quad b_{k+1}=b_k\ lambda_{k+1}=mu_k \ mu_{k+1}=a_{k+1}+0.618(b_{k+1}-a_{k+1})end{cases} ⎩⎪⎨⎪⎧​ak+1​=λk​,bk+1​=bk​λk+1​=μk​μk+1​=ak+1​+0.618(bk+1​−ak+1​)​

qquadqquad 若 f ( λ k ) < f ( μ k ) f(lambda_k)

( 3 ) qquad(3) (3) 置 k = k + 1 k=k+1 k=k+1，返回步骤 ( 2 ) (2) (2)

import numpy as np
def getInterval(f,x0,h):
(略)
def GoldenSection(f,a,b,eps=1e-4):
    ak = a
    bk = b
    Lambda = ak + 0.382*(bk-ak)
    Mu = ak + 0.618*(bk-ak)
    while True:
        if np.abs(bk-ak) < eps:
            break
        if f(Lambda) > f(Mu):
            ak = Lambda
            Lambda = Mu
            Mu = ak + 0.618*(bk-ak)
        else:
            bk = Mu
            Mu = Lambda
            Lambda = ak + 0.382*(bk-ak)
        print(ak,bk)
    return (bk+ak)/2.
if __name__ == "__main__":    
    f = lambda x: (x-1)**2 + 1        
    a,b = getInterval(f,1.3,0.4)
    print('[a,b] = [%.1f,%.1f]'%(np.round(a,2),np.round(b,2)))
    minval = GoldenSection(f,a,b)
    print(minval)

3. 最速下降法——求解最优值

qquad 考虑无约束优化问题

min ⁡   f ( x ) , x ∈ R n qquadqquadqquadmin f(boldsymbol x),qquadboldsymbol xin R^n min f(x),x∈Rn

qquad
qquad 假设在第 k k k 步时的迭代点为 x ( k ) boldsymbol x^{(k)} x(k)，采用最速下降法求迭代后继点 x ( k + 1 ) boldsymbol x^{(k+1)} x(k+1) 的过程为：

( 1 ) qquad(1) (1) 一维搜索的方向为该点的负梯度方向 d ( k ) = − ∇ f ( x ( k ) ) boldsymbol d^{(k)}=-nabla f(boldsymbol x^{(k)}) d(k)=−∇f(x(k))

( 2 ) qquad(2) (2) 计算一维搜索的步长 λ k = arg min ⁡ λ > 0 φ ( λ ) = arg min ⁡ λ > 0 f ( x ( k ) + λ d ( k ) ) lambda_k=argmin_{lambda>0}varphi(lambda)=argmin_{lambda>0} f(boldsymbol x^{(k)}+lambdaboldsymbol d^{(k)}) λk​=λ>0argmin​φ(λ)=λ>0argmin​f(x(k)+λd(k))

( 3 ) qquad(3) (3) 求得第 k + 1 k+1 k+1 步时的迭代点为 x ( k + 1 ) = x ( k ) + λ k d ( k ) boldsymbol x^{(k+1)}=boldsymbol x^{(k)}+lambda_kboldsymbol d^{(k)} x(k+1)=x(k)+λk​d(k)

qquad
算法步骤

qquad 步骤 ( 1 ) : (1): (1): 给定初始点 x ( 1 ) ∈ R n boldsymbol x^{(1)}in R^n x(1)∈Rn，设置允许误差 ε > 0 varepsilon>0 ε>0，令 k = 1 k=1 k=1

qquad 步骤 ( 2 ) : (2): (2): 计算搜索方向 d ( k ) = − ∇ f ( x ( k ) ) boldsymbol d^{(k)}=-nabla f(boldsymbol x^{(k)}) d(k)=−∇f(x(k))

qquad 步骤 ( 3 ) : (3): (3): 若 ∥ d ( k ) ∥ ≤ ε Vertboldsymbol d^{(k)}Vertlevarepsilon ∥d(k)∥≤ε，则停止计算；

qquadqquadqquad 否则，从 x ( k ) boldsymbol x^{(k)} x(k) 出发，沿着 d ( k ) boldsymbol d^{(k)} d(k) 进行一维搜索求出 λ k lambda_k λk​，使得：

f ( x ( k ) + λ k d ( k ) ) = min ⁡ λ > 0 f ( x ( k ) + λ d ( k ) ) qquadqquadqquadqquad f(boldsymbol x^{(k)}+lambda_kboldsymbol d^{(k)})=displaystylemin_{lambda>0} f(boldsymbol x^{(k)}+lambdaboldsymbol d^{(k)}) f(x(k)+λk​d(k))=λ>0min​f(x(k)+λd(k))

qquad 步骤 ( 4 ) : (4): (4): 令 x ( k + 1 ) = x ( k ) + λ k d ( k ) boldsymbol x^{(k+1)}=boldsymbol x^{(k)}+lambda_kboldsymbol d^{(k)} x(k+1)=x(k)+λk​d(k)，且 k = k + 1 k=k+1 k=k+1，转步骤 ( 2 ) (2) (2)

qquad

4. 基于数值方法的实现

qquad 与“符号计算”的实现相比，使用数值方法的主要问题：
qquad
( 1 ) qquad(1) (1) 在符号计算中， x ( k ) + λ d ( k ) boldsymbol x^{(k)}+lambdaboldsymbol d^{(k)} x(k)+λd(k) 可以直接代入到目标函数 f ( x ) f(boldsymbol x) f(x) 中，从而得到了一个单变量 λ lambda λ 的一维搜索函数 φ ( λ ) = f ( x ( k ) + λ d ( k ) ) ,   λ > 0 varphi(lambda)=f(boldsymbol x^{(k)}+lambdaboldsymbol d^{(k)}), lambda>0 φ(λ)=f(x(k)+λd(k)), λ>0

★ quadtextcolor{red}{bigstar} ★ 在《最速下降法的实现及可视化——基于matlab符号计算》中，“使用 0.618 0.618 0.618 法求单峰区间最优值”的过程是基于单变量 λ lambda λ 的一维搜索函数 φ ( λ ) = f ( x ( k ) + λ d ( k ) ) ,   λ > 0 varphi(lambda)=f(boldsymbol x^{(k)}+lambdaboldsymbol d^{(k)}), lambda>0 φ(λ)=f(x(k)+λd(k)), λ>0。
　
★ quadtextcolor{red}{bigstar} ★ 此时， “ 0.618 textbf{0.618} 0.618法”和“进退法”的求解过程都是一维的情形，相当于求一维函数 φ ′ ( λ ) = 0 varphi^{prime}(lambda)=0 φ′(λ)=0 的“解析解”
　　　 (a) text{(a)} (a) 先求一维搜索的步长 λ k lambda_k λk​
　　　 (b) text{(b)} (b) 再计算 x ( k + 1 ) = x ( k ) + λ k d ( k ) boldsymbol x^{(k+1)}=boldsymbol x^{(k)}+lambda_kboldsymbol d^{(k)} x(k+1)=x(k)+λk​d(k)

quad
( 2 ) qquad(2) (2) 因此，在使用数值方法计算的时候，“使用 0.618 0.618 0.618 法求单峰区间最优值”不再是“求解析解 ”，而是直接求有限维空间中单峰区间的最优值。

本文前两部分的代码都是一维的情形 ( x ∈ R ) (xin R) (x∈R)，而优化问题一般考虑的是有限维的情况 ( x ∈ R n ) (boldsymbol xin R^n) (x∈Rn)

qquad
∙ bullet ∙　（有限维情况的）进退法

def getUnimodalInterval(f,x0,d,h):		# 参数多了n维方向 d
    k = 0
    x1 = x0 
    f1 = f(x1)      
    while (1):
        x4 = x1 + h*d					# 在 d 方向上进退
        f4 = f(x4)
        k += 1
        if (k == 1) and (f4 > f1):
            x2 = x4
            h = -h            
            continue        
        if f4 < f1:
            x2 = x1
            x1 = x4
            h = 2*h
        if f4 > f1:
            x3 = x2
            x2 = x1
            x1 = x4
            break           
    return np.fmin(x1,x3), np.fmax(x1,x3)

∙ bullet ∙　（有限维情况的）0.618法

def GoldenSection(f,a,b,eps=1e-6):
    ak = a
    bk = b
    Lambda = ak + 0.382*(bk-ak)
    Mu = ak + 0.618*(bk-ak)
    while True:
        if np.linalg.norm(bk-ak) < eps:		# "区间长度"改用n维空间的欧氏距离
            break
        if f(Lambda) > f(Mu):
            ak = Lambda
            Lambda = Mu
            Mu = ak + 0.618*(bk-ak)
        else:
            bk = Mu
            Mu = Lambda
            Lambda = ak + 0.382*(bk-ak)        
    return (bk+ak)/2.

∙ bullet ∙　（有限维情况的）最速下降法

# 微分对应了 f = lambda x: 2*(x[0]-x0)**2 + (x[1]-x1)**2 的函数形式
def diff1(f,x,h=0.001):							# 定义一阶微分
    dx = np.array([h,0])
    dy = np.array([0,h])
    gx = (f(x+dx)-f(x-dx))/(2*h)
    gy = (f(x+dy)-f(x-dy))/(2*h)
    return np.array([gx,gy])

def steepestdescent(f,x0):
    k=1
    xk = x0
    h = 0.1
    while True:
        print('#',k)
        k = k+1
        d = -diff1(f,xk)						# 最速下降方向 d
        if np.linalg.norm(d)<0.01:
            break
        d = d/np.linalg.norm(d)
        print('d:',d)
        a,b = getUnimodalInterval(f,xk,d,h)		# 进退法求单峰区间
        print('interval:[',np.round(a,2),np.round(b,2),'], len:',np.linalg.norm(b-a))
        if np.linalg.norm(b-a)<0.5:				# 当区间长度小于0.5，进退法步长设置为0.001
            h = 0.001
        xk = GoldenSection(f,a,b)				# 求单峰区间的最小值
        print('min:',xk)
    return xk  

进退法默认采用 0.1 0.1 0.1 的初始步长，但是当区间长度非常小的时候（假设小于 0.5 0.5 0.5 ），进退法的步长也要减小到一定程度，否则进退法会失效。

qquad
∙ bullet ∙　实现代码

import numpy as np
import matplotlib.pyplot as plt
import time
def getUnimodalInterval(f,x0,d,h):
(略)
def GoldenSection(f,a,b,eps=1e-6):
(略)
def diff1(f,x,h=0.001):
(略)
def steepestdescent(f,x0): 
(略)
if __name__ == "__main__":       
    f = lambda x: 2*(x[0]-0)**2 + (x[1]-0)**2
    x0 = np.array([1,1],dtype='float')
    time0 = time.process_time()
    minval = steepestdescent(f,x0)
    time1 = time.process_time()
    print(minval,np.round(f(minval),4))
    print('time: %fs'%(time1-time0))

运行结果：

# 1
d: [-0.8944 -0.4472]
interval:[ [-1.77 -0.39] [0.37 0.69] ], len: 2.4
min: [-0.1111  0.4444]
# 2
d: [ 0.4472 -0.8944]
interval:[ [ 0.02 -0.9 ] [0.56 0.18] ], len: 1.2000000000000002
min: [ 0.3144 -0.3144]
# 3
d: [-0.8944  0.4472]
interval:[ [-1.03 -0.18] [0.05 0.36] ], len: 1.2000000000000002
min: [-0.0741  0.2963]
# 4
d: [ 0.4472 -0.8944]
interval:[ [-0.03 -0.33] [0.24 0.21] ], len: 0.6000000000000001
min: [ 0.0904 -0.0904]
# 5
d: [-0.8944  0.4472]
interval:[ [-0.18 -0.09] [0.09 0.04] ], len: 0.30000000000000004
min: [ 0.0003 -0.0012]
# 6
[ 0.0003 -0.0012] 0.0
time: 0.000000s