有限差法(Finite Difference)求梯度和Hessian Matrix(海森矩阵)的python实现
数学参考
有限差方法求导,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 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矩阵。
# -*- coding: utf-8 -*-
# @author: Dasheng Fan
# @email: fandasheng1999@163.comdef 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 gradientfor 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 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] + 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_jreturn hessian_matrix
可以通过输入函数func
和求取二阶导数的点x
,就可以输出该点处的Hessian矩阵。
hessian(func, [0,0])
得到结果
array([[2., 0.],[0., 2.]], dtype=float32)
如果和numdifftools
的结果对照,可以发现一样。但是numdifftools
非常难用,总是报错,而且速度奇慢,如果需要循环中算,更是龟速。我们的程序只需要numpy
包就能实现,非常方便好用,速度非常快。
有限差法(Finite Difference)求梯度和Hessian Matrix(海森矩阵)的python实现相关推荐
- 局部最优、梯度消失、鞍点、海森矩阵(Hessian Matric)、批梯度下降算法(btach批梯度下降法BGD、小批量梯度下降法Mini-Batch GD、随机梯度下降法SGD)
日萌社 人工智能AI:Keras PyTorch MXNet TensorFlow PaddlePaddle 深度学习实战(不定时更新) BATCH_SIZE大小设置对训练耗时的影响:1.如果当设置B ...
- 梯度、雅克比矩阵、海森矩阵、多元泰勒公式
梯度向量的表达式为: [∂f∂x1∂f∂x2...∂f∂xn]=[∂f∂x1∂f∂x2..∂f∂xn]T\left[ \begin{array} { c c } {\frac {\partial{ ...
- 【转载】矩阵求导、几种重要的矩阵及常用的矩阵求导公式
一.矩阵求导 一般来讲,我们约定x=(x1,x2,-xN)Tx=(x1,x2,-xN)T,这是分母布局.常见的矩阵求导方式有:向量对向量求导,标量对向量求导,向量对标量求导. 1.向量对向量求导 Nu ...
- 【机器学习中的矩阵求导】(六)Jacobian矩阵和Hessian矩阵
学习总结 (0)回顾矩阵向量化,和 克罗内克积的主要运算法则. (1)梯度向量是雅克比矩阵的特例. (2)Hessian矩阵是梯度向量g(x)对自变量x的Jacobian矩阵,描述了函数的局部曲率. ...
- python求向量函数的雅可比矩阵_使用python,pytorch求海森Hessian矩阵
考虑一个函数$y=f(\textbf{x}) (R^n\rightarrow R)$,y的Hessian矩阵定义如下: 考虑一个函数:$$f(x)=b^Tx+\frac{1}{2}x^{T}Ax\\其 ...
- 梯度、Hessian矩阵、Jacobian矩阵的计算
x表示为如下列向量: 一.f(x)为一维 此时其一阶导数构成的向量为梯度向量g(x),其二阶导数构成的矩阵为Hessian(海森/黑塞)矩阵G(x): 导数可以表示为: 梯度向量 Hessian矩阵 ...
- 13-反向传播法求梯度
文章目录 反向传播法求梯度 梯度下降求最小值 反向传播法求梯度 利用计算图求梯度是一种比较方便又快速的方法,如何利用计算图求梯度?先回忆一下计算图: 以 z=x2+y2z=x^2+y^2z=x2+y2 ...
- 机器学习-算法-有监督学习:EM(最大期望值算法)<=> MLE(最大似然估计法)【关系类似“梯度下降法”<=>“直接求导法”】【EM“梯度下降”:先初始化一个随机值,然后通过迭代不断靠近真实值】
机器学习-算法-有监督学习:EM(最大期望值算法)<=> MLE(最大似然估计法)[关系类似"梯度下降法"<=>"直接求导法"][EM& ...
- 二阶偏微分方程组 龙格库塔法_有限单元法(Finite Element Method)实现声波方程模拟(Part 2)...
2.1 前言 承接上一篇文章,前面我们已经介绍了一维声波方程有限元求解: 蓝不是蓝:有限单元法(Finite Element Method)实现声波方程模拟(Part 1)zhuanlan.zhih ...
最新文章
- DI容器是代码污染者
- 手机 服务器 推送消息推送消息,推送信息到手机的pushover使用方法及sample code
- 前端学习(3268):js中this在类中的表现
- 2021快手奢侈品行业数据价值报告
- 详解Python生成器函数和生成器对象的原理和用法
- 牛逼! IDEA 2020 要本土化,真的是全中文了!中国开发者话语权越来越大了
- linux系统查看java安装路径
- 《高质量程序设计指南》读书笔记
- vs运行c语言代码快捷键,VS2010快捷键
- C语言基础题练习10道
- SP namespace (sp.js)
- 『开发技巧』MacbookM1芯片深度学习环境配置最全教程:简明安装开发TensorFlow与PyTorch
- Application Virtualization 4.5 部署【3】
- Arduino案例实操 -- OLED中文显示(IIC)
- 手把手教你React Native接入聊天IM即时通讯功能-源码分享
- 【总结】二手书网站开发总结(业余时间开发)
- Fiddler 抓包工具使用详解
- 疯子网页采集器之提取内容教程
- python输入,Python中的基本输入和输出
- 做影视后期需要学习哪些行业软件?