數(shù)學(xué)參考

有限差方法求導(dǎo)，Finite Difference Approximations of Derivatives，是數(shù)值計算中常用的求導(dǎo)方法。數(shù)學(xué)上也比較簡單易用。本文主要針對的是向量值函數(shù)，也就是$$f(x):{\mathbb{R}}^n\to \mathbb{R}$$當(dāng)然，普通的標(biāo)量值函數(shù)是向量值函數(shù)的一種特例。

本文采用的數(shù)學(xué)參考是：有限差方法
參考的主要是Central Difference Approximations小節(jié)中的Second-order derivatives based on gradient calls的那個公式。

代碼

用法

將下面代碼中的Hessian矩陣一節(jié)中的Hessian函數(shù)直接復(fù)制到你的代碼中，然后就可以按照用法示例使用。

特別要注意，eps的選擇比較關(guān)鍵，直接決定了有限差方法的精度。建議大家根據(jù)函數(shù)參數(shù)的數(shù)量級動態(tài)的設(shè)置，例如某參數(shù)變化范圍1-10，就可以設(shè)置為0.001；而某參數(shù)變化范圍為0-0.0001，則可設(shè)置為0.000001，之類的。

用法示例

def func(x):
    x_0 = x[0]
    x_1 = x[1]
    return x_0**2 + x_1**2
hessian(func, [0,0], esp = [0.01, 0.01])

得到結(jié)果：

array([[2., 0.],
       [0., 2.]], dtype=float32)

函數(shù)主體

準(zhǔn)備

本文的方法只需要numpy包，幾乎可以說不需要任何包，而且不受到什么限制，只要滿足輸入格式就能求取，比所謂autograd，numdifftools好用的多。

梯度函數(shù)

為了求Hessian矩陣，本文采用的方法需要首先求取梯度。首先需要有一個函數(shù)func，示例的func如下：

def func(x, **args):
    x_0 = x[0]
    x_1 = x[1]
    return x_0**2 + x_1**2

該函數(shù)是一個${\mathbb{R}}^2\to \mathbb{R}$的函數(shù)。將該函數(shù)輸入進下面的函數(shù)grad_func_generator中之后，就可以返回梯度函數(shù)，支持在任何一點求取梯度。這里輸入x應(yīng)該是一個列表，是各個維度的輸入。例如x = [0,0].

def grad_func_generator(func, eps = 0.00001):
    def gradient_func(point):
        n_var = len(point)
        gradient = np.zeros(n_var, np.float32)
        # nth gradient
        for i in range(n_var):
            # 初始化左點和右點，同時不改變原來的展開點
            left_point = point.copy()
            right_point = point.copy()
            left_point[i] = point[i] - eps
            right_point[i] = point[i] + eps
            gradient[i] = (func(right_point) - func(left_point))/(2*eps)
        return gradient
    return gradient_func

求取梯度：

grad_f = grad_func_generator(func) # 生成梯度函數(shù)
grad_f([1,1])

可以得到結(jié)果：

array([2., 2.], dtype=float32)

Hessian矩陣

利用已經(jīng)實現(xiàn)的梯度函數(shù)，可以實現(xiàn)Hessian矩陣。

# -*- coding: utf-8 -*-
# @author: Dasheng Fan
# @email: fandasheng1999@163.com

def hessian(func, point = [0, 0], eps = [0.001, 0.001]):
    """
    Hessian matrix of func at expendung point.
    """
    n_var = len(point)
    def grad_func_generator(func):
        def gradient_func(point):
            gradient = np.zeros(n_var, np.float32)
            # nth gradient
            for i in range(n_var):
                # 初始化左點和右點，同時不改變原來的展開點
                left_point = point.copy()
                right_point = point.copy()
                left_point[i] = point[i] - eps[i]
                right_point[i] = point[i] + eps[i]
                gradient[i] = (func(right_point) - func(left_point))/(2*eps[i])
            return gradient
        return gradient_func

    grad_func = grad_func_generator(func)
    hessian_matrix = np.zeros((n_var, n_var), np.float32)
    for i in range(n_var):
        for j in range(n_var):
            # 第一項
            left_point_j = point.copy()
            right_point_j = point.copy()
            right_point_j[j] = point[j] + eps[j]
            left_point_j[j] = point[j] - eps[j]
            diff_i = (grad_func(right_point_j)[i] - grad_func(left_point_j)[i])/(4*eps[j])
            # 第二項
            left_point_i = point.copy()
            right_point_i = point.copy()
            right_point_i[i] = point[i] + eps[i]
            left_point_i[i] = point[i] - eps[i]
            diff_j = (grad_func(right_point_i)[j] - grad_func(left_point_i)[j])/(4*eps[i])

            hessian_matrix[i, j] = diff_i + diff_j

    return hessian_matrix

可以通過輸入函數(shù)func和求取二階導(dǎo)數(shù)的點x，就可以輸出該點處的Hessian矩陣。

hessian(func, [0,0])

得到結(jié)果

array([[2., 0.],
       [0., 2.]], dtype=float32)

如果和numdifftools的結(jié)果對照，可以發(fā)現(xiàn)一樣。但是numdifftools非常難用，總是報錯，而且速度奇慢，如果需要循環(huán)中算，更是龜速。我們的程序只需要numpy包就能實現(xiàn)，非常方便好用，速度非?？?。文章來源地址http://www.zghlxwxcb.cn/news/detail-778492.html

到了這里，關(guān)于有限差法(Finite Difference)求梯度和Hessian Matrix(海森矩陣)的python實現(xiàn)的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！