🌞

最优化:Python3 梯度下降法、共轭梯度下降法、牛顿法

> [最优化:Python3 模拟退火算法](https://hxd.red/2019/12/31/python3-la-methode-du-recuit-simule/)

无约束最优化

什么是无约束最优化?看下面这张图。图中有两个点:一个为全局最小值点,另一个为局部最小值点

目前大部分无约束最优化算法只能保证求取局部最小值点。在实际应用中,许多情形被抽象为函数形式后均为凸函数,对于凸函数来说局部最小值点即为全局最小值点,因此只要能求得这类函数的一个最小值点,该点就为全局最小值点。

在初始点选择好之后,就可以按照各种不同的无约束最优化求解算法,求解最小值点。求解过程中需要从初始点开始沿“哪个方向”(方向)以及“走多远”(步长)到达下一个点。

导数与梯度

考虑一座高度在((x, y) ) 点是(H(x, y) ) 的山。 ( H ) 这一点的梯度是在该点坡度(或者说斜度)最陡的方向。梯度的大小告诉我们坡度到底有多陡。标量函数 ( f\colon \mathbb {R} ^{n}\mapsto \mathbb {R} ) 的梯度表示为: ( \nabla f ) 或 ( {grad} f ) ,其中 ( \nabla ) (nabla) 表示向量微分算子。

\nabla f在三维直角坐标系中表示为:

( \nabla f={\begin{pmatrix}{\frac {\partial f}{\partial x}},{\frac {\partial f}{\partial y}},{\frac {\partial f}{\partial z}}\end{pmatrix}}={\frac {\partial f}{\partial x}}\mathbf {i} +{\frac {\partial f}{\partial y}}\mathbf {j} +{\frac {\partial f}{\partial z}}\mathbf {k} )

i, j, k 为标准的单位向量,分别指向 x, y 跟 z 坐标的方向。举例来讲,函数 ( f(x,y,z)=2x+3y^{2}-\sin(z) ) 的梯度为:

( \nabla f={\begin{pmatrix}{2},{6y},{-\cos(z)}\end{pmatrix}}=2\mathbf {i} +6y\mathbf {j} -\cos(z)\mathbf {k} )

梯度定义如下:函数在某一点的梯度是这样一个向量,它的方向与取得最大方向导数的方向一致,而它的模为方向导数的最大值。这里注意三点:

  • 梯度是一个向量,即有方向有大小;
  • 梯度的方向是最大方向导数的方向;
  • 梯度的值是最大方向导数的值。

梯度下降法

梯度下降法是一个一阶最优化算法,通常也称为最速下降法。 要使用梯度下降法找到一个函数的局部极小值,必须向函数上当前点对应梯度(或者是近似梯度)的反方向的规定步长距离点进行迭代搜索。

假设f(x)是 ( R∗ ) 上具有一阶连续偏导数的函数,需要求解的无约束最优化问题是

( \min_{x\varepsilon R^n}f(x) )

(x^* ) 是目标函数的极小值点。

梯度下降法是一种迭代算法,选取适当的初值 ( x^(0) ) ,反复迭代,更新x的值,进行目标函数的极小化,直至收敛。由于梯度方向是函数增长最快的方向,那么负梯度方向就是函数值下降最快的方向。 因此,以负梯度方向作为极小化的下降方向,在迭代的每一步,以负梯度方向来更新x的值,从而达到减小函数值目的,这种方法就是梯度下降法(也叫最速下降法)。

由于 ( f(x) ) 具有一阶连续偏导数,若第k次迭代值为 ( x^(k) ) , ( \nabla{f(x^{(k)})} ) 为 ( f(x) ) 在 ( x^(k) ) 的梯度。

第k+1次迭代值 ( x^(k+1) ) :

( x^{(k+1)} = x^{(k)}+\rho ^{(k)}d^{(k)} )

其中, ( d^{(k)} ) 是搜索方向,取负梯度方向 ( d^{(k)}= \nabla{f(x^{(k)})} ) , ( \rho ^{(k)} ) 是步长,这个值优先选择常量。

梯度下降法的算法

Même si ces m´méthodes sont conceptuellement très simples et qu’elles peuvent être programmées directement, elles sont souvent lentes dans la pratique. Elles convergent mais sous des conditions de convergence souvent complexes. 在实际求解过程中,它的收敛速度是比较慢的。所以就出现了下面两种方法。

这里有几幅很生动的动图可以直观地查看它是如何寻找最小值的。

共轭梯度下降法

共轭梯度法是介于最速下降法与牛顿法之间的一个方法, 它仅需一阶导数信息,但克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点。

共轭梯度法的基本思想是把共轭性与最速下降法相结合,利用已知点处的梯度构造一组共轭方向,并沿这组方向进行搜索,求出目标函数的极小点。根据共轭方向的基本性质,这种方法具有二次终止性。

共轭

设G为对称正定矩阵,若 ( a_m^TGa_n ) ,则称 ( a_m ) 和 ( a_n ) 为“G共轭”。当G为单位向量时,有 ( a_m a_n =0 ) ,所以“共轭”是“正交”的推广。

那么沿着一系列的共轭方向做迭代,这些共轭方向组成的集合叫做共轭方向集,这就是共轭方向法。即对于二次正定目标函数,从任意点 ( x_0 ) 出发,沿任意下降方向 ( a_m ) 做直线搜索得到 ( x_1 ) ,沿与 ( a_m ) 共轭的方向 ( a_n ) 做直线搜索,反复迭代,即可得到目标函数的极小值点。

Python3 代码参考

首先说一下,这方面的东西我是一个外行,也没有特别深入的研究过,这些例子是学校“数学、信息“课程(其实就是”最优化“)的内容。其理论部分详见这篇博文。内容和解释一定有不详尽、甚至错误之处,可以比对多方资料加深理解。本文中做列出的代码在Python3.7环境下可用并能够正常执行。

# -*- coding: utf-8 -*-
"""
Methodes de gradients.
"""
###################
##### 导入模块 #####
###################
import numpy as np
import matplotlib.pyplot as plt
import itertools
import pickle  #python2中为cPickle,python3中为pickle
import scipy.linalg as spl

################################
##### 定义函数 #####
################################

def function(A,b,c,xx):
    shape=xx[0].shape #shape查看矩阵或者数组的维数
    ZZ=np.zeros_like(xx[0]) #其维度与矩阵xx[0]一致,并为其初始化为全0
    for item in itertools.product(*map(range,shape)): #map()会根据提供的函数对指定序列做映射。笛卡尔积,相当于嵌套的for循环.
        vect=list()
        for elem in xx:
            vect.append(elem[item])
        vec=np.array(vect,ndmin=1)
        ZZ[item]=0.5*np.dot(vec.T,A.dot(vec))-np.dot(vec.T,b) + c
    return ZZ

def gradient(A,b,x):
    return A.dot(x)-b

def pas_fixe(x0,fonction,pas=1.0e-2,tol=1.0e-10,itermax=10000):
    A=fonction['A']
    b=fonction['b']
    c=fonction['c']

    #***** 初始化 Initialisation *****
    xx=[x0]
    dir = - gradient(A,b,xx[0])
    residu = [ np.linalg.norm(gradient(A,b,xx[0])) ]
    k=0
    
    while residu[-1] >= tol and k <= itermax:
        #----- Calcul de x(k+1) -----
        xx.append(xx[-1]+pas*dir)
    
        #----- Calcul de la nouvelle direction de descente d(x+1) -----
        dir = - gradient(A,b,xx[-1])
    
        #----- Calcul du residu r(k+1) -----
        residu.append(np.linalg.norm(gradient(A,b,xx[-1])))
        k = k + 1
    
    return {'xx': np.array(xx), 'residu': np.array(residu),'k': k}

def conjugate(x0,fonction,tol=1.0e-10,itermax=10000):
    A=fonction['A']
    b=fonction['b']
    c=fonction['c']

    #***** 初始化 Initialisation *****
    xx = [x0]
    dir = - gradient(A,b,xx[0])
    residu =  [ np.linalg.norm(gradient(A,b,xx[0])) ]
    k = 0
    
    while residu[-1] >= tol and k <= itermax:
        #----- Calcul de rho(k) -----
        rho = - gradient(A,b,xx[-1]).dot(dir)/ A.dot(dir).dot(dir)
    
        #----- Calcul de x(k+1) -----
        xx.append(xx[-1]+rho*dir)
    
        #----- Calcul de la nouvelle direction de descente d(x+1) -----
        beta = gradient(A,b,xx[-1]).dot(dir) / A.dot(dir).dot(dir)
        dir = - gradient(A,b,xx[-1]) + beta * dir
    
        #----- Calcul du residu r(k+1) -----
        residu.append(np.linalg.norm(gradient(A,b,xx[-1])))
        k = k+1
    
    return  {'xx': np.array(xx), 'residu': np.array(residu),'k': k}

def plot(xx,res):
    plt.figure()

    X1, X2 = np.meshgrid(np.linspace(-5.0,5.0,101),np.linspace(-5.0,5.0,101))
    Z=function(A,b,c,[X1,X2])
    
    plt.contour(X1,X2,Z)
    
    plt.plot(xx.T[0],xx.T[1],'k-x')
    
    plt.axes().set_aspect('equal')
    
    plt.figure()
    plt.plot(res)
    plt.yscale('log')
    plt.grid()
    plt.xlabel('Iterations')
    plt.ylabel(r'$||\nabla f(x) ||_2$')
    plt.title('Convergence')
    
    plt.show()

#######################
##### 执行部分 #####
#######################

if __name__=="__main__":
    """
    在本练习中,最小值是[-0.1429,-0.4286]
    Pour cet exercice le minimum est [-0.1429,-0.4286]
    """
    x0=np.array([1.,1.])
    A=np.array([[4.,1.],[1.,2.]])
    b=np.array([-1.,-1.])
    c=np.zeros(1)

    fonction={'A':A, 'b':b, 'c':c}
    
    #----- 固定步长 Pas fixe -----
    cas_01=pas_fixe(x0,fonction,1.0e-1,1.0e-6,10000)
    print (cas_01['xx'][-1], len(cas_01['xx'])-1)
    
    plot(cas_01['xx'],cas_01['residu'])
    
    #----- 共轭梯度 Gradient conjugue -----
    cas_02=conjugate(x0,fonction,1.0e-6,10000)
    print (cas_02['xx'][-1], len(cas_02['xx'])-1)
    
    plot(cas_02['xx'],cas_02['residu'])

牛顿法

牛顿法又称为牛顿-拉弗森方法,它是一种在实数域和复数域上近似求解方程的方法。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(y)=0的根。

  1. 先找到一个接近函数f(x)零点的x0,计算对应的f(x0)f'(x0)。f’表示f的导数。
  2. 计算( 0=(x-x_{0})\cdot f'(x_{0})+f(x_{0}))
  3. 从而得到新的x坐标x1x1会比x0更接近方程f(x)=0的解。
  4. 再次迭代计算( x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}} )

参考资料

本文的主要参考资料:梯度下降法,共轭梯度下降法,随机梯度下降法、Introduction à l’optimisation numérique

更新记录

  • 2019.10.07 v0.1.0
updatedupdated2020-01-252020-01-25