有限差法(Finite Difference)求梯度和Hessian Matrix(海森矩阵)的python实现

数学参考

有限差方法求导，Finite Difference Approximations of Derivatives，是数值计算中常用的求导方法。数学上也比较简单易用。本文主要针对的是向量值函数，也就是$$f(x):{\mathbb{R}}^n\to \mathbb{R}$$当然，普通的标量值函数是向量值函数的一种特例。

本文采用的数学参考是：有限差方法
参考的主要是Central Difference Approximations小节中的Second-order derivatives based on gradient calls的那个公式。

代码

准备

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

梯度函数

为了求Hessian矩阵，本文采用的方法需要首先求取梯度。首先需要有一个函数func，示例的func如下：

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

该函数是一个${\mathbb{R}}^2\to \mathbb{R}$的函数。将该函数输入进下面的函数grad_func_generator中之后，就可以返回梯度函数，支持在任何一点求取梯度。这里输入x应该是一个列表，是各个维度的输入。例如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 gradientfor i in range(n_var):# 初始化左点和右点，同时不改变原来的展开点left_point = point.copy()right_point = point.copy()left_point[i] = point[i] - epsright_point[i] = point[i] + epsgradient[i] = (func(right_point) - func(left_point))/(2*eps)return gradientreturn gradient_func

求取梯度：

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

可以得到结果：

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

Hessian矩阵

利用已经实现的梯度函数，可以实现Hessian矩阵。

def hessian(func, point = [0, 0], eps = 0.00001):"""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 gradientfor i in range(n_var):# 初始化左点和右点，同时不改变原来的展开点left_point = point.copy()right_point = point.copy()left_point[i] = point[i] - epsright_point[i] = point[i] + epsgradient[i] = (func(right_point) - func(left_point))/(2*eps)return gradientreturn gradient_funcgrad_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] + epsleft_point_j[j] = point[j] - epsdiff_i = (grad_func(right_point_j)[i] - grad_func(left_point_j)[i])/(4*eps)# 第二项left_point_i = point.copy()right_point_i = point.copy()right_point_i[i] = point[i] + epsleft_point_i[i] = point[i] - epsdiff_j = (grad_func(right_point_i)[j] - grad_func(left_point_i)[j])/(4*eps)hessian_matrix[i, j] = diff_i + diff_jreturn hessian_matrix

可以通过输入函数func和求取二阶导数的点x，就可以输出该点处的Hessian矩阵。

hessian(func, [0,0])

得到结果

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

如果和numdifftools的结果对照，可以发现一样。但是numdifftools非常难用，总是报错。我们的程序只需要numpy包就能实现，非常方便好用。