微动弹性带方法——续鞍点

1. 鞍点
2. 微动弹性带方法
3. 例题
- 3.1 SciPy求解最小值
- 3.2 IPyVolume可视化NumPy数组
- 3.3 梯度
- 3.4 寻找最小值
- 3.5 寻找鞍点
- 3.6 Climbing Image Nudged Elastic Band

1. 鞍点

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

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

2. 微动弹性带方法

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

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

参考资料：Nudged Elastic Band Method

3. 例题

考虑一个三变量函数（见下方代码），坐标区间为[-1, 1]。

寻找这个函数的在(0.5, 0.5, 0.5)和(-0.5, -0.5, -0.5)附近的两个局部极小值，以及两个极小值之间最低能量路径上的鞍点。data是目标函数

import numpy as npdef gaussian(pos, pos0):return np.exp(-np.sum((pos-pos0)**2))def data(pos):return gaussian(pos, np.array([0.1, -0.1, -0.1])) - gaussian(pos, np.array([-0.5, -0.5, -0.5])) - gaussian(pos, np.array([0.5, 0.5, 0.5]))

3.1 SciPy求解最小值

采用BFGS求解两点附近最小值。

from scipy import optimizepos = np.array([0.5, 0.5, 0.5])
print(optimize.minimize(data, pos, method='BFGS'))pos = np.array([-0.5, -0.5, -0.5])
print(optimize.minimize(data, pos, method='BFGS'))

3.2 IPyVolume可视化NumPy数组

打印梯度

N = 50
X = Y = Z = np.linspace(-1, 1, N)
data_grid = np.zeros([N, N, N])for i in range(N):for j in range(N):for k in range(N):data_grid[k, j, i] = data(np.array([X[i], Y[j], Z[k]]))print(np.min(data_grid))

3D Volume
三维可视化

import ipyvolume as ipvipv.figure()
ipv.volshow(data_grid)
ipv.view(-40)
ipv.show()

2D Slice
二维展平

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

3.3 梯度

def gradient(pos, delta = 0.01):gradient_x = (data(pos + np.array([delta, 0, 0])) - data(pos - np.array([delta, 0, 0]))) / 2 / deltagradient_y = (data(pos + np.array([0, delta, 0])) - data(pos - np.array([0, delta, 0]))) / 2 / deltagradient_z = (data(pos + np.array([0, 0, delta])) - data(pos - np.array([0, 0, delta]))) / 2 / deltareturn np.array([gradient_x, gradient_y, gradient_z])

3.4 寻找最小值

随机生成三维点

np.random.seed(0)
pos = np.random.rand(3) * 2 - 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_minimal(pos, rate, err):pos_diff = 1data_diff = 1while np.abs(data_diff) > err:pos_new = pos - np.array(gradient(pos)) * ratedata_diff = data(pos_new) - data(pos)pos_diff = np.max(np.abs(np.array(gradient(pos)) * rate))pos = pos_newreturn pos
pos = find_minimal(pos, 0.01,  1e-10)
print(pos, data(pos))

for i in range(5):pos = np.random.rand(3) * 2 - 1pos = find_minimal(pos, 0.001,  1e-10)print(pos, data(pos))

记录两个局部最优点

pos_A = np.array([0.60486467, 0.67087103, 0.67164885])
pos_B = np.array([-0.73136024, -0.64616776, -0.64521874])

3.5 寻找鞍点

n = 11
images = np.zeros([n, 3])for i in range(n):images[i] = (pos_B - pos_A)/(n - 1) * i + pos_Aprint(images)

可视化等间距

ipv.quickscatter(images[:,0], images[:,1], images[:,2], size=5, marker="sphere")

定义弹力的合力大小和方向

dist_AB = np.sqrt(np.sum((pos_A - pos_B)**2))
dist_image = dist_AB / (n - 1)def spring_force(image_before, image, image_next, k = 2.0):dist_before = np.sqrt(np.sum((image - image_before)**2))force_before = (dist_before - dist_image) * k direction_before = (image - image_before)/dist_beforedist_next = np.sqrt(np.sum((image_next - image)**2))force_next = (dist_image - dist_next) * kdirection_next = (image_next - image)/dist_nextforce = force_before*direction_before + force_next*direction_nextdirection = (image_next - image_before) / np.sqrt(np.sum((image_next - image_before)**2))return force, direction

采用弹力重力的合力进行迭代，返回鞍点位置索引，判断条件是珠子位置导数和目标函数导数的最大值

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)rate = 0.1
err = 1e-8def NEB(rate, err):n_step = 0pos_diff = 1data_diff = 1while pos_diff > err or data_diff > err:old_pos = imagesold_saddle = np.max([data(images[i]) for i in range(n)])for i in range(1, n-1):s_force[i], direction[i] = spring_force(images[i-1], images[i], images[i+1])s_force[i] = np.dot(s_force[i], direction[i]) * direction[i]  # Vector decompositiong_force[i] = gradient(images[i])g_force[i] = g_force[i] - np.dot(g_force[i], direction[i]) * direction[i]  # Vector decompositionimages[i] -= (s_force[i]+g_force[i]) * ratenew_pos = imagesnew_saddle = np.max([data(images[i]) for i in range(n)]) idx_saddle = np.argmax([data(images[i]) for i in range(n)])pos_diff = np.max(np.abs(new_pos - old_pos))data_diff = np.abs(new_saddle - old_saddle)n_step += 1return n_step, idx_saddlen_step, idx_saddle = NEB(rate, err)
print("n_step    ", n_step)
print("idx_saddle", idx_saddle)
print("images    ", images[idx_saddle])
print("data      ", data(images[idx_saddle]))ipv.quickscatter(images[:,0], images[:,1], images[:,2], size=5, marker="sphere")
gradient(images[idx_saddle])

中间构型仍然存在沿着最低能量路径方向的梯度，说明其并不在鞍点位置。我们需要进一步精修其位置。“生活就是这样给你一个不怎么精确的解，还需要优化迭代出最优解”

3.6 Climbing Image Nudged Elastic Band

采用单纯的重力更新位置，判断条件是珠子受重力方向的最大值

def cNEB(image, direction):g_force = 1while np.max(np.abs(g_force)) > err:g_force = gradient(image)g_force = np.dot(g_force, direction) * directionimage += g_force * ratereturn imagesaddle_direction = (images[idx_saddle+1] - images[idx_saddle-1]) / np.sqrt(np.sum((images[idx_saddle+1] - images[idx_saddle-1])**2))
saddle_point = cNEB(images[idx_saddle], saddle_direction)print(saddle_point)
print(gradient(saddle_point))
print(data(saddle_point))

精度还是不错的，梯度非常接近零啦！干饭喽

参考文献来自桑鸿乾老师的课件：科学计算和人工智能

【微动弹性带方法——续鞍点】相关推荐

的微波感知_上海交大彭志科教授团队研发：微波微动监测与智能感知技术
上海9月18日电(葛俊俊) 准确监测方舱医院大量感染患者的生命状况,精确"诊断"大桥工程结构是否存在安全隐患,随时随地获取独居老人在家的健康体征--上海交通大学彭志科教授团队研发的 ...
使用html语言检测鼠标微动是否发生双击
寒假以来打了几天游戏,却发现有时候枪会走火,便得知是鼠标微动出现了双击现象.可以采用以下方法检测微动是否产生双击.将其中的html代码保存至记事本,然后修改拓展名为html再打开即可,最终比对实际点击 ...
为什么通过微服务的方法构建应用程序？
作为软件开发人员,我们已知道思考如何将应用程序因数分解成组件部分. 这是对象导向.软件抽象和组件化的中心模式. 现在,这种因数分解往往以共享库和技术层之间的类与接口呈现. 通常采用一种分层方法,有后端 ...
电脑ppt录制微课软件哪个好电脑ppt录制微课的方法
如今线上课程已经逐渐成为线下课程的补充,拓宽知识面,让学生能够学到更多知识.微课是线上课程里比较方便观看的一类,制作起来也很便捷,很多人会直接用ppt来制作微课,简单快速又能传播知识.今天就来分享一下 ...
性价比哪家强？富勒G93S光磁微动鼠标深度评测
光磁微动打破传统壁垒,更好的解决了由于弹簧片老化带来的双击问题,快速响应,5ms触发,寿命长.富勒黑科技G93S将光磁微动再次演绎. 要玩就玩高清无码版每次用逆天的4K显示器看视频,我总是担心他们会 ...
4am永远鼠标按键设置_可换微动的三模游戏鼠标：华硕ROG烈刃2体验
双十一期间购买的商品基本已经开始陆陆续续收到货了,也是时候开启疯狂晒单模式了.今天跟大家分享的是华硕ROG烈刃2三模电竞鼠标,初上市的时候就放在购物车了,也算是种草很久的一款外设产品,双十一期间价格优 ...
工业机器人、工艺夹具、送料机械手、电火花镗磨机床、半自动钻床、机械手、套筒、铣床升降台、精密播种机、卧式组合钻床、六自由度微动机器人、花生收获机、山茶采摘平台、车载起重机、锤式破碎机、螺旋输送机……
机械手-四自由度的工业机器人[毕业设计(论文)+CAD图纸] 工艺夹具-减速器箱体零件工艺及加工Φ120外圆的夹具设计车床-数控车床系统XY工作台与控制系统设计机械手-送料机械手设计及Solidw ...
钉钉企业微应用调试方法
解决钉钉企业微应用需要反复部署调试的方法钉钉微应用调试方法启动你的本地项目(前提要后端允许本地的id地址访问[关于后端如何允许前端id地址访问的操作请看最下方~]) 首先下载钉钉RC版进到RC版 ...
快速微课制作方法和技巧
很多老师常常说:我也想做微课,可是我对电脑操作不太熟练,做不成.其实微课最重要的是设计思路,有了巧妙的构思,总能找到简单的计算完成它,不会使用Camtasia Studio复杂的录屏工具,可以使用Co ...
智能探测静止微动，云望爱希活体存在感应器，雷达触发技术应用
活体存在感应器具有极高的灵敏度和准确性,活体存在感应器能够实现对人在正常工作生活中行走.移动.抬头.转身等微小动作的识别,以及呼吸信号.非睡眠状态下的人体存在的探测,具有极高的灵敏度和准确性. 同时用 ...

【微动弹性带方法——续鞍点】