数学参考
有限差方法求导,Finite Difference Approximations of Derivatives
,是数值计算中常用的求导方法。数学上也比较简单易用。本文主要针对的是向量值函数,也就是
f
(
x
)
:
R
n
→
R
f(x):\mathbb{R^n}\rightarrow \mathbb{R}
f(x):Rn→R当然,普通的标量值函数是向量值函数的一种特例。
本文采用的数学参考是:有限差方法
参考的主要是Central Difference Approximations
小节中的Second-order derivatives based on gradient calls
的那个公式。
代码
用法
将下面代码中的Hessian矩阵一节中的Hessian函数直接复制到你的代码中,然后就可以按照用法示例使用。
特别要注意,eps的选择比较关键,直接决定了有限差方法的精度。建议大家根据函数参数的数量级动态的设置,例如某参数变化范围1-10,就可以设置为0.001;而某参数变化范围为0-0.0001,则可设置为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])
得到结果:
array([[2., 0.],
[0., 2.]], dtype=float32)
函数主体
准备
本文的方法只需要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
该函数是一个
R
2
→
R
\mathbb{R^2}\rightarrow \mathbb{R}
R2→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 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) # 生成梯度函数
grad_f([1,1])
可以得到结果:
array([2., 2.], dtype=float32)
Hessian矩阵
利用已经实现的梯度函数,可以实现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
可以通过输入函数func
和求取二阶导数的点x
,就可以输出该点处的Hessian矩阵。
hessian(func, [0,0])
得到结果文章来源:https://www.toymoban.com/news/detail-778492.html
array([[2., 0.],
[0., 2.]], dtype=float32)
如果和numdifftools
的结果对照,可以发现一样。但是numdifftools
非常难用,总是报错,而且速度奇慢,如果需要循环中算,更是龟速。我们的程序只需要numpy
包就能实现,非常方便好用,速度非常快。文章来源地址https://www.toymoban.com/news/detail-778492.html
到了这里,关于有限差法(Finite Difference)求梯度和Hessian Matrix(海森矩阵)的python实现的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!