1. 简介

2. 作业题目

2.1 要求

2.2 提示

3. 解答

3.1 建立数据来源

3.2 设计当前点梯度

3.3 寻找最小值

4. 寻找鞍点

4.1 珠子等间距连接

4.2 设置弹簧弹力

4.3 微动弹性带

5. 总结

1. 简介

鞍点是一种特殊的驻点。对于多变量函数，在鞍点位置，函数沿任意方向的导数都为0，但函数并不是最大值或者最小值。我们关注一类特殊的鞍点，在这个位置，函数在某一方向上是最大值，但是在剩余所有方向上是极小值。

寻找鞍点在科学和工程研究中有很多应用。一个常用的例子是地形图，地势高度取决于水平坐标，因此这是一个双变量函数。假设在起伏的地势中有两个盆地（对应于函数的局部极小值）A和B。一个人想要从A出发到达B，在连接A和B的所有可能的路径中，哪一条路径走过的地势最高点最低？这个问题的实质就是寻找这个双变量函数的鞍点（或者一个更常见的名称，过渡态）。

2. 作业题目

考虑一个三变量函数（见下方代码），寻找这个函数的在(0.5, 0.5, 0.5)和(-0.5, -0.5, -0.5)附近的两个局部极小值，以及两个极小值之间路径上(0.7, 0.3, 0.5)附近的鞍点。

2.1 要求

数值结果误差小于1e-3。
不允许使用numpy和scipy之外的数学库，需手动实现算法。
不允许对提供的data的生成函数做符号分析（例如解析求导），data只能作为一个黑盒子使用。
为了模拟真实场景（data中数据点需要借助实验或者模拟仿真高成本获得），不允许对data进行遍历或暴力搜索。
代码优越性的评价取决于data中数据的读取次数，越低越好（读取次数越低，对应于真实场景中解决问题速度越快）。

为了方便大家对函数值的分布有直观认识，这里使用ipyvolume进行了3D可视化。

import numpy as np
import ipyvolume as ipvx = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
z = np.linspace(-1, 1, 100)
Y, Z, X = np.meshgrid(x, y, z)def gaussian(x0, y0, z0):return np.exp(-(X - x0)**2 - (Y - y0)**2 - (Z - z0)**2)data = gaussian(0.1, -0.1, -0.1) - gaussian(-0.5, -0.5, -0.5) - gaussian(0.5, 0.5, 0.5)ipv.figure()
ipv.volshow(data)
ipv.view(-40)
ipv.show()

借助IPython提供的交互控件，可以使用下方的拖动条查看data的在xy方向上的二维数据切片。

from ipywidgets import interactive
import matplotlib.pyplot as pltdef slice_z(i):fig, ax = plt.subplots(figsize=(8, 8))im = ax.imshow(data[i,:,:], vmin=-1, vmax=1, cmap=plt.get_cmap('gist_rainbow'))ct = ax.contour(data[i,:,:])bar = plt.colorbar(im)plt.show()interact_plot = interactive(slice_z, i=(0, 99))
interact_plot

2.2 提示

寻找过渡态的一个常用算法是微动弹性带（Nudged Elastic Band）。它的核心思想是，将初始坐标和终态坐标用若干个中间态（例如11个）连接起来，然后让这些中间态沿着函数梯度的反方向移动（类似于小球在地形图中沿着山坡向下移动）；为了避免这些中间态都收敛到附近的局部极小（类似于小球都落入了盆地），相邻中间态之间用一根虚拟的弹簧连接，从而迫使相邻中间态有一定的间距。当这个小球弹簧链（微动弹性带）移动停止时，其所在位置就是所谓的最低能量路径（minimal energy path），其中间函数值最大的位置就是鞍点或者过渡态。

在迭代计算过程中，中间态的移动同时受函数梯度和弹簧弹力的调控。为了保持中间态的间距尽量不改变，以及虚拟弹簧不影响路径正确性，可以对函数梯度和弹簧弹力进行矢量分解。其中，函数梯度只保留垂直于路径的分量；弹簧弹力只保留沿着路径的分量。

3. 解答

3.1 建立数据来源

import numpy as np
import ipyvolume as ipvx = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
z = np.linspace(-1, 1, 100)
Y, Z, X = np.meshgrid(x, y, z)def gaussian(x,y,z,x0, y0, z0):return np.exp(-(x- x0)**2 - (y - y0)**2 - (z - z0)**2)
def data(x,y,z):return gaussian(x,y,z,0.1, -0.1, -0.1) - gaussian(x,y,z,-0.5, -0.5, -0.5) - gaussian(x,y,z,0.5, 0.5, 0.5)
data(0,0,0)

返回查找的data数据

3.2 设计当前点梯度

三个方向分别采用相邻两点的差求解对应的偏导

##梯度
def gradient(x,y,z,delta=0.01):gradient_x = (data(x+delta,y,z) - data(x-delta,y,z))/2/deltagradient_y = (data(x,y+delta,z) - data(x,y-delta,z))/2/deltagradient_z = (data(x,y,z+delta) - data(x,y,z-delta))/2/deltareturn gradient_x,gradient_y,gradient_z
gradient(0,0,0)

返回对应的点的梯度

3.3 寻找最小值

np.random.seed(0)#定义随机种子数
pos = np.random.rand(3)*2-1#生成点的坐标（-1,1）之间三个随机数
print(pos,data(*pos),gradient(*pos))
rate=0.1#下降速率
pos_new=pos-np.array(gradient(*pos))*rate#采用梯度下降法
print(pos_new,data(*pos_new))
def find_minimum(pos,rate,err):#寻找局部最小值的点data_diff = 1.0while np.abs(data_diff) > err:pos_new = pos - np.array(gradient(*pos))*ratedata_diff = data(*pos_new) - data(*pos)pos = pos_newreturn pos
pos = find_minimum(pos,rate,err = 1e-6)
print(pos,data(*pos))

可以随机寻找多次得到局部最优点，原则采用相似的取最小值的位置

for i in range(10):pos = np.random.rand(3)*2-1pos = find_minimum(pos,rate,err = 1e-6)print(pos,data(*pos))

选择第四个和第五个最小值点，然后保存A和B两点

pos_A = np.array([0.60501534,0.67242841,0.67002837])#三个坐标用逗号
pos_B = np.array([-0.73177095,-0.64662581,-0.64406625])

4. 寻找鞍点

4.1 珠子等间距连接

首先采用9个珠子等间距连接两点

import ipyvolume as ipy#导入画图库
n = 9#采用9粒珠子将两个局部最小值点圈连起来
images = np.zeros([n,3])#设置9个位置为零清空内存for i in range(n):#连接A\B两点images[i] = (pos_B - pos_A)/(n - 1) * i + pos_Aipy.quickscatter(images[:,0], images[:,1], images[:,2], size = 5, marker = "sphere")#可视化9个珠子

显示

4.2 设置弹簧弹力

两端珠子已经固定，中间的1 一个珠子分别连接两个，先算单个弹力（与弹簧长度相关）和方向（采用单位方向），direction = np.zeros_like(force)其维度与矩阵force一致，并为其初始化为全0，后期叠加。

dist_AB = np.sqrt(np.sum((pos_A - pos_B)**2))#AB两点之间距离
dist_avg = dist_AB / (n - 1)#两珠子之间平均距离
n = 9
images = np.zeros([n,3])
for i in range(n):#用9粒珠子均匀连接A\B两点images[i] = (pos_B - pos_A)/(n - 1) * i + pos_Adef spring_force(image_before, image, image_next, k = 2.0):#前三个三参数是相邻的位置点，k为弹力增益系数dist_before = np.sqrt(np.sum((image - image_before)**2))force_before = (dist_before - dist_avg)*kdirection_before = (image - image_before)/dist_beforedist_next = np.sqrt(np.sum((image - image_next)**2))force_next = (dist_next - dist_avg)*kdirection_next = (image - image_next)/dist_nextforce = force_before*direction_before + force_next*direction_nextif np.sum(force**2)< 1e-8:direction = np.zeros_like(force)else:direction = force / np.sqrt(np.sum(force**2))return force,directionspring_force(images[4], images[5], images[6], k = 2.0)#测试4/5/6点力和方向

返回，没有变换之前力和方向都几乎为零

4.3 微动弹性带

珠子在重力和弹力限制下运动，判断所有合内力是否在允许范围，允许结束，反而言之。

n = 9
images = np.zeros([n,3])
for i in range(n):images[i] = (pos_B - pos_A)/(n - 1) * i + pos_As_force = np.zeros_like(images)
direction = np.zeros_like(images)
g_force = np.zeros_like(images)
Saddle_Point = np.zeros([n])rate = 0.1#滚动速率
err =1e-5print(data(*images[5]))
def Saddle(rate, err):n_step = 0data_diff = 1.0while np.abs(data_diff) > err:#判断所有合内力是否在允许范围data_diff = 0for i in range(1,n-1):old_data = data(*images[i])s_force[i],direction[i] = spring_force(images[i-1], images[i], images[i+1], k = 2.0)s_force[i] = np.dot(s_force[i], direction[i]) * direction[i]#力矢量分解g_force[i] = gradient(*images[i])#重力与梯度相关g_force[i] = g_force[i] - np.dot(g_force[i], direction[i]) * direction[i]images[i] -= (s_force[i] + g_force[i]) * rate#珠子在重力和弹力限制下运动new_data = data(*images[i])data_diff=data_diff + (new_data - old_data)#计算所有合力n_step+= 1return n_step
n_step = Saddle(rate, err)
for i in range(n):Saddle_Point[i]=data(*images[i])print(Saddle_Point[i])
print(n_step)
print('Saddle Point value:'+str(max(Saddle_Point)))#打印鞍点值
ipy.quickscatter(images[:,0], images[:,1], images[:,2], size = 5, marker = "sphere")

经过2158次迭代，得到“Saddle Point value:-0.09851946983280313”。返回如下

完整代码，需要修改地方是pos_A,pos_B的坐标

#导入基本库
import numpy as np
import ipyvolume as ipv#定义三坐标线性矩阵系列
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
z = np.linspace(-1, 1, 100)np.random.seed(0)#定义随机种子数
pos = np.random.rand(3)*2-1#生成点的坐标（-1,1）之间三个随机数
rate=0.1#下降速率#定义高斯函数
def gaussian(x,y,z,x0, y0, z0):return np.exp(-(x- x0)**2 - (y - y0)**2 - (z - z0)**2)
#定义数据库
def data(x,y,z):return gaussian(x,y,z,0.1, -0.1, -0.1) - gaussian(x,y,z,-0.5, -0.5, -0.5) - gaussian(x,y,z,0.5, 0.5, 0.5)
data(0,0,0)
##定义梯度函数
def gradient(x,y,z,delta=0.01):gradient_x = (data(x+delta,y,z) - data(x-delta,y,z))/2/deltagradient_y = (data(x,y+delta,z) - data(x,y-delta,z))/2/deltagradient_z = (data(x,y,z+delta) - data(x,y,z-delta))/2/deltareturn gradient_x,gradient_y,gradient_z
gradient(0,0,0)#寻找局部最小值函数
def find_minimum(pos,rate,err):data_diff = 1.0while np.abs(data_diff) > err:pos_new = pos - np.array(gradient(*pos))*ratedata_diff = data(*pos_new) - data(*pos)pos = pos_newreturn pos#随机搜索10次寻找局部最小值
for i in range(10):pos = np.random.rand(3)*2-1pos = find_minimum(pos,rate,err = 1e-6)print(pos,data(*pos))#注意此次需要根据自己计算的差异性修改为局部最小值，保存两个最小值
pos_A = np.array([0.60501534,0.67242841,0.67002837])#三个坐标用逗号
pos_B = np.array([-0.73177095,-0.64662581,-0.64406625])#准备计算弹力
dist_AB = np.sqrt(np.sum((pos_A - pos_B)**2))#AB两点之间距离
dist_avg = dist_AB / (n - 1)#两珠子之间平均距离
n = 9
images = np.zeros([n,3])
for i in range(n):#用9粒珠子均匀连接A\B两点images[i] = (pos_B - pos_A)/(n - 1) * i + pos_A#定义弹力大小和方向
def spring_force(image_before, image, image_next, k = 2.0):#前三个三参数是相邻的位置点，k为弹力增益系数dist_before = np.sqrt(np.sum((image - image_before)**2))force_before = (dist_before - dist_avg)*kdirection_before = (image - image_before)/dist_beforedist_next = np.sqrt(np.sum((image - image_next)**2))force_next = (dist_next - dist_avg)*kdirection_next = (image - image_next)/dist_nextforce = force_before*direction_before + force_next*direction_nextif np.sum(force**2)< 1e-8:direction = np.zeros_like(force)else:direction = force / np.sqrt(np.sum(force**2))return force,direction#定义珠子之间的弹力、方向、重力、珠子目标函数值
s_force = np.zeros_like(images)
direction = np.zeros_like(images)
g_force = np.zeros_like(images)
Saddle_Point = np.zeros([n])rate = 0.1#滚动速率
err =1e-5print(data(*images[5]))
def Saddle(rate, err):n_step = 0data_diff = 1.0while np.abs(data_diff) > err:#判断所有合内力是否在允许范围data_diff = 0for i in range(1,n-1):old_data = data(*images[i])s_force[i],direction[i] = spring_force(images[i-1], images[i], images[i+1], k = 2.0)s_force[i] = np.dot(s_force[i], direction[i]) * direction[i]#力矢量分解g_force[i] = gradient(*images[i])#重力与梯度相关g_force[i] = g_force[i] - np.dot(g_force[i], direction[i]) * direction[i]images[i] -= (s_force[i] + g_force[i]) * rate#珠子在重力和弹力限制下运动new_data = data(*images[i])data_diff=data_diff + (new_data - old_data)#计算所有合力n_step+= 1return n_step
n_step = Saddle(rate, err)
for i in range(n):Saddle_Point[i]=data(*images[i])print(Saddle_Point[i])
print(n_step)
print('Saddle Point value:'+str(max(Saddle_Point)))#打印鞍点值
ipy.quickscatter(images[:,0], images[:,1], images[:,2], size = 5, marker = "sphere")

升级版，不需要手动更新，直接一键运行

#导入基本库
import numpy as np
import ipyvolume as ipy#定义三坐标线性矩阵系列
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
z = np.linspace(-1, 1, 100)n = 9#珠子个数
#np.random.seed(0)#定义随机种子数
#pos = np.random.rand(3)*2-1#生成点的坐标（-1,1）之间三个随机数
poss = np.zeros([n+1, 4])#
gradient_rate=0.1#梯度下降速率
find_minimum_err = 1e-6#局部最小值误差#定义高斯函数
def gaussian(x,y,z,x0, y0, z0):return np.exp(-(x- x0)**2 - (y - y0)**2 - (z - z0)**2)#定义数据库
def data(x,y,z):return gaussian(x,y,z,0.1, -0.1, -0.1) - gaussian(x,y,z,-0.5, -0.5, -0.5) - gaussian(x,y,z,0.5, 0.5, 0.5)
data(0,0,0)##定义梯度函数
def gradient(x,y,z,delta=0.01):gradient_x = (data(x+delta,y,z) - data(x-delta,y,z))/2/deltagradient_y = (data(x,y+delta,z) - data(x,y-delta,z))/2/deltagradient_z = (data(x,y,z+delta) - data(x,y,z-delta))/2/deltareturn gradient_x,gradient_y,gradient_z#寻找局部最小值函数
def find_minimum(pos,gradient_rate,find_minimum_err):data_diff = 1.0while np.abs(data_diff) > find_minimum_err:pos_new = pos - np.array(gradient(*pos))*gradient_ratedata_diff = data(*pos_new) - data(*pos)pos = pos_newreturn pos#定义弹力大小和方向
def spring_force(image_before, image, image_next, k = 2.0):#前三个三参数是相邻的位置点，k为弹力增益系数dist_before = np.sqrt(np.sum((image - image_before)**2))force_before = (dist_before - dist_avg)*kdirection_before = (image - image_before)/dist_beforedist_next = np.sqrt(np.sum((image - image_next)**2))force_next = (dist_next - dist_avg)*kdirection_next = (image - image_next)/dist_nextforce = force_before*direction_before + force_next*direction_nextif np.sum(force**2)< 1e-8:direction = np.zeros_like(force)else:direction = force / np.sqrt(np.sum(force**2))return force,direction#采用微动弹性带寻找过渡态
def Saddle(roll_rate, roll_err):n_step = 0data_diff = 1.0while np.abs(data_diff) > roll_err:#判断所有合内力是否在允许范围data_diff = 0for i in range(1,n-1):old_data = data(*images[i])s_force[i],direction[i] = spring_force(images[i-1], images[i], images[i+1], k = 2.0)s_force[i] = np.dot(s_force[i], direction[i]) * direction[i]#力矢量分解g_force[i] = gradient(*images[i])#重力与梯度相关g_force[i] = g_force[i] - np.dot(g_force[i], direction[i]) * direction[i]images[i] -= (s_force[i] + g_force[i]) * roll_rate#珠子在重力和弹力限制下运动new_data = data(*images[i])data_diff=data_diff + (new_data - old_data)#计算所有合力n_step+= 1return n_stepnp.random.seed(0)#定义随机种子数
#pos = np.random.rand(3)*2-1
#随机搜索10次寻找局部最小值
for i in range(10):#np.random.seed(0)#定义随机种子数pos = np.random.rand(3)*2-1pos = find_minimum(pos,gradient_rate,find_minimum_err)print(pos,data(*pos))poss[i, 0:3] = posposs[i, 3] = data(*pos)#注意此次需要根据自己计算的差异性修改为局部最小值，保存两个最小值
index_a = np.squeeze(np.where(poss[:, 3] == min(poss[:, 3])))
pos_A = poss[index_a, 0:3]
pos_A = np.array([round(i,8) for i in pos_A])
for i in range(10):if abs(data(*pos_A) - poss[i, 3]) < 1e-3:poss[i, 3] = 0
index_b = np.squeeze(np.where(poss[:, 3] == min(poss[:, 3])))
pos_B = poss[index_b, 0:3]
pos_B = np.array([round(i,8) for i in pos_B])
print(pos_A,pos_B)
#pos_A = np.array([0.60501534,0.67242841,0.67002837])#三个坐标用逗号
#pos_B = np.array([-0.73177095,-0.64662581,-0.64406625])
#print(pos_A,pos_B)dist_AB = np.sqrt(np.sum((pos_A - pos_B)**2))#AB两点之间距离
dist_avg = dist_AB / (n - 1)#两珠子之间平均距离
images = np.zeros([n,3])#清空所以珠子的位置信息s_force = np.zeros_like(images)#定义珠子之间的弹力
direction = np.zeros_like(images)#定义珠子之间的方向
g_force = np.zeros_like(images)#定义珠子之间的重力
Saddle_Point = np.zeros([n])#定义珠子目标函数值roll_rate = 0.1#滚动速率
roll_err =1e-5#准备计算弹力
for i in range(n):#用9粒珠子均匀连接A\B两点images[i] = (pos_B - pos_A)/(n - 1) * i + pos_A
n_step = Saddle(roll_rate, roll_err)
for i in range(n):Saddle_Point[i]=data(*images[i])print(Saddle_Point[i])
print(n_step)
print('Saddle Point value:'+str(max(Saddle_Point)))#打印鞍点值
ipy.quickscatter(images[:,0], images[:,1], images[:,2], size = 5, marker = "sphere")

详细代码可见：考虑一个三变量函数（见下方代码），寻找这个函数的在(0.5,0.5,0.5)和(-0.5,-0.5,-0.5)附近的两个-Python文档类资源-CSDN下载

5. 总结

本文完成了鞍点Saddle Point Locator求解，后期会分享更多有趣的操作从而实现对外部世界进行感知，充分认识这个有机与无机的环境，科学地合理地进行创作和发挥效益，然后为人类社会发展贡献一点微薄之力。

鞍点Saddle Point Locator相关推荐

(一)神经网络训练不起来怎么办：局部最小值(local minia)与鞍点(saddle point)
Optimization的时候,怎么把gradient descent做的更好? 1.局部最小值(Local minima)与鞍点(saddle point) 所谓的saddle point其实就是g ...
局部最优点+鞍点+学习率的调节
深度学习系列第一篇局部最优点+鞍点+学习率的调节文章目录深度学习系列局部最优点和鞍点一样吗? 自动调整学习率我们为什么要调整学习率? 总结局部最优点和鞍点一样吗? 这是不一样的,局部最优 ...
c语言二维数组找鞍点,C语言,二维数组找鞍点
还是以前写过的东西.. 鞍点是什么?百度出来的东西鞍点(Saddle point)在微分方程中,沿着某一方向是稳定的,另一条方向是不稳定的奇点,叫做鞍点.在泛函中,既不是极大值点也不是极小值点的临界 ...
梯度下降优化算法综述与PyTorch实现源码剖析
现代的机器学习系统均利用大量的数据,利用梯度下降算法或者相关的变体进行训练.传统上,最早出现的优化算法是SGD,之后又陆续出现了AdaGrad.RMSprop.ADAM等变体,那么这些算法之间又有哪些 ...
第二章：1、函数求导
机器学习中问题的求解往往转化为优化问题. 1.函数极限函数的导数是通过极限来定义和表达的.极限理解: 向量的2范数:表示点到原点的欧氏距离,就是向量的长度.点的邻域表示以点a为心,以为半径的一个d维 ...
系列笔记 | 深度学习连载（4）：优化技巧（上）
点击上方"AI有道",选择"星标"公众号重磅干货,第一时间送达深度学习中我们总结出 5 大技巧: 1. Adaptive Learning Rate 我们先 ...
非线性动力学_非线性动力学特辑低维到高维的联通者
序言: 本文将以维度为主线, 带量大家进入非线性动力学的世界. 文章数学部分不需要全部理解, 理解思维方法为主非线性动力学,是物理学的思维进入传统方法所不能解决的问题的一座丰碑.它可以帮助我们理解不 ...
卖萌屋算法工程师思维导图part3—深度学习篇
卖萌屋的妹子们(划掉)作者团整理的算法工程师思维导图,求职/自我提升/查漏补缺神器.该手册一共分为数据结构与算法.数学基础.统计机器学习和深度学习四个部分. 下面是第三部分深度学习的内容~ 公众号后台 ...
非线性最小二乘通俗易懂解释
转https://www.cnblogs.com/leexiaoming/p/7257198.html备份用 1. 非线性最小二乘介绍 1.1. 最小二乘问题回顾: 在上一篇博客中我们知道最小二乘问题 ...

鞍点Saddle Point Locator

1. 简介

2. 作业题目

2.1 要求

2.2 提示

3. 解答

3.1 建立数据来源

3.2 设计当前点梯度

3.3 寻找最小值

4. 寻找鞍点

4.1 珠子等间距连接

4.2 设置弹簧弹力

4.3 微动弹性带

5. 总结

鞍点Saddle Point Locator相关推荐

最新文章

热门文章