目录
1. 简介
2. 作业题目
2.1 要求
2.2 提示
3. 解答
3.1 建立数据来源
3.2 设计当前点梯度
3.3 寻找最小值
4. 寻找鞍点
4.1 珠子等间距连接
4.2 设置d簧d力
4.3 微动d性带
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 ipv
x = 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 plt
def 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_plot2.2 提示
寻找过渡态的一个常用算法是微动d性带（Nudged Elastic Band）。它的核心思想是，将初始坐标和终态坐标用若干个中间态（例如11个）连接起来，然后让这些中间态沿着函数梯度的反方向移动（类似于小球在地形图中沿着山坡向下移动）；为了避免这些中间态都收敛到附近的局部极小（类似于小球都落入了盆地），相邻中间态之间用一根虚拟的d簧连接，从而迫使相邻中间态有一定的间距。当这个小球d簧链（微动d性带）移动停止时，其所在位置就是所谓的最低能量路径（minimal energy path），其中间函数值最大的位置就是鞍点或者过渡态。
在迭代计算过程中，中间态的移动同时受函数梯度和d簧d力的调控。为了保持中间态的间距尽量不改变，以及虚拟d簧不影响路径正确性，可以对函数梯度和d簧d力进行矢量分解。其中，函数梯度只保留垂直于路径的分量；d簧d力只保留沿着路径的分量。
3. 解答 3.1 建立数据来源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)
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/delta
gradient_y = (data(x,y+delta,z) - data(x,y-delta,z))/2/delta
gradient_z = (data(x,y,z+delta) - data(x,y,z-delta))/2/delta
return 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.0
while np.abs(data_diff) > err:
pos_new = pos - np.array(gradient(*pos))*rate
data_diff = data(*pos_new) - data(*pos)
pos = pos_new
return pos
pos = find_minimum(pos,rate,err = 1e-6)
print(pos,data(*pos))返回
可以随机寻找多次得到局部最优点，原则采用相似的取最小值的位置
for i in range(10):
pos = np.random.rand(3)*2-1
pos = 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])首先采用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_A
ipy.quickscatter(images[:,0], images[:,1], images[:,2], size = 5, marker = "sphere")#可视化9个珠子显示
4.2 设置d簧d力两端珠子已经固定，中间的1 一个珠子分别连接两个，先算单个d力（与d簧长度相关）和方向（采用单位方向），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_A
def spring_force(image_before, image, image_next, k = 2.0):#前三个三参数是相邻的位置点，k为d力增益系数
dist_before = np.sqrt(np.sum((image - image_before)**2))
force_before = (dist_before - dist_avg)*k
direction_before = (image - image_before)/dist_before
dist_next = np.sqrt(np.sum((image - image_next)**2))
force_next = (dist_next - dist_avg)*k
direction_next = (image - image_next)/dist_next
force = force_before*direction_before + force_next*direction_next
if np.sum(force**2)< 1e-8:
direction = np.zeros_like(force)
else:
direction = force / np.sqrt(np.sum(force**2))
return force,direction
spring_force(images[4], images[5], images[6], k = 2.0)#测试4/5/6点力和方向
返回，没有变换之前力和方向都几乎为零
4.3 微动d性带珠子在重力和d力限制下运动，判断所有合内力是否在允许范围，允许结束，反而言之。
n = 9
images = np.zeros([n,3])
for i in range(n):
images[i] = (pos_B - pos_A)/(n - 1) * i + pos_A
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-5
print(data(*images[5]))
def Saddle(rate, err):
n_step = 0
data_diff = 1.0
while np.abs(data_diff) > err:#判断所有合内力是否在允许范围
data_diff = 0
for 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#珠子在重力和d力限制下运动
new_data = data(*images[i])
data_diff=data_diff + (new_data - old_data)#计算所有合力
n_step+= 1
return 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/delta
gradient_y = (data(x,y+delta,z) - data(x,y-delta,z))/2/delta
gradient_z = (data(x,y,z+delta) - data(x,y,z-delta))/2/delta
return gradient_x,gradient_y,gradient_z
gradient(0,0,0)
#寻找局部最小值函数
def find_minimum(pos,rate,err):
data_diff = 1.0
while np.abs(data_diff) > err:
pos_new = pos - np.array(gradient(*pos))*rate
data_diff = data(*pos_new) - data(*pos)
pos = pos_new
return pos
#随机搜索10次寻找局部最小值
for i in range(10):
pos = np.random.rand(3)*2-1
pos = 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])
#准备计算d力
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
#定义d力大小和方向
def spring_force(image_before, image, image_next, k = 2.0):#前三个三参数是相邻的位置点，k为d力增益系数
dist_before = np.sqrt(np.sum((image - image_before)**2))
force_before = (dist_before - dist_avg)*k
direction_before = (image - image_before)/dist_before
dist_next = np.sqrt(np.sum((image - image_next)**2))
force_next = (dist_next - dist_avg)*k
direction_next = (image - image_next)/dist_next
force = force_before*direction_before + force_next*direction_next
if 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-5
print(data(*images[5]))
def Saddle(rate, err):
n_step = 0
data_diff = 1.0
while np.abs(data_diff) > err:#判断所有合内力是否在允许范围
data_diff = 0
for 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#珠子在重力和d力限制下运动
new_data = data(*images[i])
data_diff=data_diff + (new_data - old_data)#计算所有合力
n_step+= 1
return 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/delta
gradient_y = (data(x,y+delta,z) - data(x,y-delta,z))/2/delta
gradient_z = (data(x,y,z+delta) - data(x,y,z-delta))/2/delta
return gradient_x,gradient_y,gradient_z
#寻找局部最小值函数
def find_minimum(pos,gradient_rate,find_minimum_err):
data_diff = 1.0
while np.abs(data_diff) > find_minimum_err:
pos_new = pos - np.array(gradient(*pos))*gradient_rate
data_diff = data(*pos_new) - data(*pos)
pos = pos_new
return pos
#定义d力大小和方向
def spring_force(image_before, image, image_next, k = 2.0):#前三个三参数是相邻的位置点，k为d力增益系数
dist_before = np.sqrt(np.sum((image - image_before)**2))
force_before = (dist_before - dist_avg)*k
direction_before = (image - image_before)/dist_before
dist_next = np.sqrt(np.sum((image - image_next)**2))
force_next = (dist_next - dist_avg)*k
direction_next = (image - image_next)/dist_next
force = force_before*direction_before + force_next*direction_next
if 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 = 0
data_diff = 1.0
while np.abs(data_diff) > roll_err:#判断所有合内力是否在允许范围
data_diff = 0
for 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#珠子在重力和d力限制下运动
new_data = data(*images[i])
data_diff=data_diff + (new_data - old_data)#计算所有合力
n_step+= 1
return n_step
np.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-1
pos = find_minimum(pos,gradient_rate,find_minimum_err)
print(pos,data(*pos))
poss[i, 0:3] = pos
poss[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")本文完成了鞍点Saddle Point Locator求解，后期会分享更多有趣的 *** 作从而实现对外部世界进行感知，充分认识这个有机与无机的环境，科学地合理地进行创作和发挥效益，然后为人类社会发展贡献一点微薄之力。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）