是否可以矢量化 scipy.optimize.fminbound?

2024-04-26

我有一些按时间参数化的轨迹数据点,我想知道每个点到拟合它们的曲线的最短距离。似乎有几种方法可以解决这个问题(例如here https://kitchingroup.cheme.cmu.edu/blog/2013/02/14/Find-the-minimum-distance-from-a-point-to-a-curve/ or here https://stackoverflow.com/a/19869154/5472354),但我选择了scipy.optimize.fminbound https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fminbound.html?highlight=scipy%20optimize%20fminbound#scipy.optimize.fminbound,如提议的here https://stackoverflow.com/a/19102112/5472354(因为这似乎是最巧妙的方法,而且因为我实际上让它发挥了作用):

import numpy as np
import pandas as pd
from scipy.optimize import fminbound


data = np.array(
    [
        (212.82275865, 1650.40828168, 0., 0),
        (214.22056952, 1649.99898924, 10.38, 0),
        (212.86786868, 1644.25228805, 116.288, 0),
        (212.78680031, 1643.87461108, 122.884, 0),
        (212.57489485, 1642.01124032, 156.313, 0),
        (212.53483954, 1641.61858242, 162.618, 0),
        (212.43922274, 1639.58782771, 196.314, 0),
        (212.53726315, 1639.13842423, 202.619, 0),
        (212.2888428, 1637.24641296, 236.306, 0),
        (212.2722447, 1636.92307229, 242.606, 0),
        (212.15559302, 1635.0529813, 276.309, 0),
        (212.17535631, 1634.60618711, 282.651, 0),
        (212.02545613, 1632.72139574, 316.335, 0),
        (211.99988779, 1632.32053329, 322.634, 0),
        (211.33419846, 1631.07592039, 356.334, 0),
        (211.58972239, 1630.21971902, 362.633, 0),
        (211.70332122, 1628.2088542, 396.316, 0),
        (211.74610735, 1627.67591368, 402.617, 0),
        (211.88819518, 1625.67310022, 436.367, 0),
        (211.90709414, 1625.39410321, 442.673, 0),
        (212.00090348, 1623.42655008, 476.332, 0),
        (211.9249017, 1622.94540583, 482.63, 0),
        (212.34321938, 1616.32949197, 597.329, 0),
        (213.09638942, 1615.2869643, 610.4, 0),
        (219.41313491, 1580.22234313, 1197.332, 0),
        (220.38660128, 1579.20043302, 1210.37, 0),
        (236.35472669, 1542.30863041, 1798.267, 0),
        (237.41755384, 1541.41679119, 1810.383, 0),
        (264.08373622, 1502.66620597, 2398.244, 0),
        (265.65655239, 1501.64308908, 2410.443, 0),
        (304.66999824, 1460.94068336, 2997.263, 0),
        (306.22550945, 1459.75817211, 3010.38, 0),
        (358.88879764, 1416.472238, 3598.213, 0),
        (361.14046402, 1415.40942931, 3610.525, 0),
        (429.96379858, 1369.7972467, 4198.282, 0),
        (432.06565776, 1368.22265539, 4210.505, 0),
        (519.30493383, 1319.01141844, 4798.277, 0),
        (522.12134083, 1317.68234967, 4810.4, 0),
        (630.00294242, 1265.05368942, 5398.236, 0),
        (633.67624272, 1263.63633508, 5410.431, 0),
        (766.29767476, 1206.91262814, 5997.266, 0),
        (770.78300649, 1205.48393374, 6010.489, 0),
        (932.92308019, 1145.87780431, 6598.279, 0),
        (937.54373403, 1141.55438694, 6609.525, 0),
    ], dtype=[
        ('x', 'f8'), ('y', 'f8'), ('t', 'f8'), ('dmin', 'f8'),
    ]
)

# fyi my data comes as a structured array; unfortunately, simply passing
# data[['x', 'y']] to np.polyfit does not work. using
# pd.DataFrame(data[['x', y']]).values seems to be the easiest solution:
# https://stackoverflow.com/a/26175750/5472354
coeffs = np.polyfit(
    data['t'], pd.DataFrame(data[['x', 'y']]).values, 3
)


def curve(t):
    # this can probably also be done in one statement, but I don't know how
    x = np.polyval(coeffs[:, 0], t)
    y = np.polyval(coeffs[:, 1], t)

    return x, y


def f(t, p):
    x, y = curve(t)
    return np.hypot(x - p['x'], y - p['y'])


# instead of this:
for point in data:
    tmin = fminbound(f, -50, 6659.525, args=(point, ))
    point['dmin'] = f(tmin, point)


# do something like this:
# tmin = fminbound(f, -50, 6659.525, args=(data, ))
# data['dmin'] = f(tmin, data)

但正如你所看到的,我使用for-循环来计算每个数据点的最短距离,这会显着减慢我的程序速度,因为这会执行数千次。因此我想向量化操作/提高性能,但还没有找到方法。有与此相关的帖子(例如here https://stackoverflow.com/q/59371930/5472354 or here https://stackoverflow.com/q/10828477/5472354),但我不知道建议的解决方案如何适用于我的情况。


不,不可能矢量化fminbound因为它期望单变量的标量函数。但是,您仍然可以通过重新表述底层数学优化问题来向量化循环:

The N标量优化问题

min f_1(t) s.t. t_l <= t <= t_u
min f_2(t) s.t. t_l <= t <= t_u
.
.
.
min f_N(t) s.t. t_l <= t <= t_u

对于标量函数f_i等价于一个优化问题 在N变量:

min f_1(t_1)**2 + ... + f_N(t_N)**2  s.t. t_l <= t_i <= t_u for all i = 1, .., N

可以通过以下方式解决scipy.optimize.minimize https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize。根据您的整个算法,您可以使用这种方法来进一步消除更多循环,即您只解决一个大规模优化问题,而不是数千个标量优化问题。

清理代码后,可以按如下方式完成:

import numpy as np
from scipy.optimize import minimize

data = np.array([
    (212.82275865, 1650.40828168, 0., 0),
    (214.22056952, 1649.99898924, 10.38, 0),
    (212.86786868, 1644.25228805, 116.288, 0),
    (212.78680031, 1643.87461108, 122.884, 0),
    (212.57489485, 1642.01124032, 156.313, 0),
    (212.53483954, 1641.61858242, 162.618, 0),
    (212.43922274, 1639.58782771, 196.314, 0),
    (212.53726315, 1639.13842423, 202.619, 0),
    (212.2888428, 1637.24641296, 236.306, 0),
    (212.2722447, 1636.92307229, 242.606, 0),
    (212.15559302, 1635.0529813, 276.309, 0),
    (212.17535631, 1634.60618711, 282.651, 0),
    (212.02545613, 1632.72139574, 316.335, 0),
    (211.99988779, 1632.32053329, 322.634, 0),
    (211.33419846, 1631.07592039, 356.334, 0),
    (211.58972239, 1630.21971902, 362.633, 0),
    (211.70332122, 1628.2088542, 396.316, 0),
    (211.74610735, 1627.67591368, 402.617, 0),
    (211.88819518, 1625.67310022, 436.367, 0),
    (211.90709414, 1625.39410321, 442.673, 0),
    (212.00090348, 1623.42655008, 476.332, 0),
    (211.9249017, 1622.94540583, 482.63, 0),
    (212.34321938, 1616.32949197, 597.329, 0),
    (213.09638942, 1615.2869643, 610.4, 0),
    (219.41313491, 1580.22234313, 1197.332, 0),
    (220.38660128, 1579.20043302, 1210.37, 0),
    (236.35472669, 1542.30863041, 1798.267, 0),
    (237.41755384, 1541.41679119, 1810.383, 0),
    (264.08373622, 1502.66620597, 2398.244, 0),
    (265.65655239, 1501.64308908, 2410.443, 0),
    (304.66999824, 1460.94068336, 2997.263, 0),
    (306.22550945, 1459.75817211, 3010.38, 0),
    (358.88879764, 1416.472238, 3598.213, 0),
    (361.14046402, 1415.40942931, 3610.525, 0),
    (429.96379858, 1369.7972467, 4198.282, 0),
    (432.06565776, 1368.22265539, 4210.505, 0),
    (519.30493383, 1319.01141844, 4798.277, 0),
    (522.12134083, 1317.68234967, 4810.4, 0),
    (630.00294242, 1265.05368942, 5398.236, 0),
    (633.67624272, 1263.63633508, 5410.431, 0),
    (766.29767476, 1206.91262814, 5997.266, 0),
    (770.78300649, 1205.48393374, 6010.489, 0),
    (932.92308019, 1145.87780431, 6598.279, 0),
    (937.54373403, 1141.55438694, 6609.525, 0)])

# the coefficients
coeffs = np.polyfit(data[:, 2], data[:, 0:2], 3)

# the points
points = data[:, :2]

# vectorized version of your objective function
# i.e. evaluates f_1, ..., f_N
def f(t, points):
    poly_x = np.polyval(coeffs[:, 0], t)
    poly_y = np.polyval(coeffs[:, 1], t)
    return np.hypot(poly_x - points[:, 0], poly_y - points[:, 1])

# the scalar objective function we want to minimize
def obj_vec(t, points):
    vals = f(t, points)
    return np.sum(vals**2)

# variable bounds
bnds = [(-50, 6659.525)]*len(points)

# solve the optimization problem
res = minimize(lambda t: obj_vec(t, points), x0=np.zeros(len(points)), bounds=bnds)
dmins = f(res.x, points)

为了进一步加速优化,强烈建议将目标函数的精确梯度传递给minimize。目前,梯度是通过有限差分来近似的,速度相当慢:

In [7]: %timeit res = minimize(lambda t: obj_vec(t, points), x0=np.zeros(len(points)), bounds=bnds)
91.5 ms ± 868 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Since the gradient can easily be computed by the chain rule, I'll leave it as an exercise for the reader :). By the chain rule, the gradient reads as

def grad(t, points):
    poly_x = np.polyval(coeffs[:, 0], t)
    poly_y = np.polyval(coeffs[:, 1], t)
    poly_x_deriv = np.polyval(np.polyder(coeffs[:, 0], m=1), t)
    poly_y_deriv = np.polyval(np.polyder(coeffs[:, 1], m=1), t)
    return 2*poly_x_deriv*(poly_x - points[:, 0]) + 2*(poly_y_deriv)*(poly_y - points[:, 1])

并将其传递给minimize显着减少运行时间:

In [9]: %timeit res = minimize(lambda t: obj_vec(t, points), jac=lambda t: grad(t, points), x0=np.zeros(len(points)),
     ...: bounds=bnds)
6.13 ms ± 63.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

最后,正如您在评论中已经指出的那样,设置另一个起点也可能会导致迭代次数减少。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

是否可以矢量化 scipy.optimize.fminbound? 的相关文章

  • 尽管极其懒惰,但如何在 Python 中模拟 IMAP 服务器?

    我很好奇是否有一种简单的方法来模拟 IMAP 服务器 例如imaplib模块 在Python中 without做很多工作 是否有预先存在的解决方案 理想情况下 我可以连接到现有的 IMAP 服务器 进行转储 并让模拟服务器在真实的邮箱 电子
  • Python BigQuery 存储。并行读取多个流

    我有以下玩具代码 import pandas as pd from google cloud import bigquery storage v1beta1 import os import google auth os environ G
  • Django REST序列化器:创建对象而不保存

    我已经开始使用 Django REST 框架 我想做的是使用一些 JSON 发布请求 从中创建一个 Django 模型对象 然后使用该对象而不保存它 我的 Django 模型称为 SearchRequest 我所拥有的是 api view
  • DreamPie 不适用于 Python 3.2

    我最喜欢的 Python shell 是DreamPie http dreampie sourceforge net 我想将它与 Python 3 2 一起使用 我使用了 添加解释器 DreamPie 应用程序并添加了 Python 3 2
  • 如何使用 Scrapy 从网站获取所有纯文本?

    我希望在 HTML 呈现后 可以从网站上看到所有文本 我正在使用 Scrapy 框架使用 Python 工作 和xpath body text 我能够获取它 但是带有 HTML 标签 而且我只想要文本 有什么解决办法吗 最简单的选择是ext
  • keras加载模型错误尝试将包含17层的权重文件加载到0层的模型中

    我目前正在使用 keras 开发 vgg16 模型 我用我的一些图层微调 vgg 模型 拟合我的模型 训练 后 我保存我的模型model save name h5 可以毫无问题地保存 但是 当我尝试使用以下命令重新加载模型时load mod
  • 在循环中每次迭代开始时将变量重新分配给原始值(在循环之前定义)

    在Python中 你使用 在每次迭代开始时将变量重新分配给原始值 在循环之前定义 时 也就是说 original 1D o o o for i in range 0 3 new original 1D revert back to orig
  • 在 NumPy 中获取 ndarray 的索引和值

    我有一个 ndarrayA任意维数N 我想创建一个数组B元组 数组或列表 其中第一个N每个元组中的元素是索引 最后一个元素是该索引的值A 例如 A array 1 2 3 4 5 6 Then B 0 0 1 0 1 2 0 2 3 1 0
  • python pandas 中的双端队列

    我正在使用Python的deque 实现一个简单的循环缓冲区 from collections import deque import numpy as np test sequence np array range 100 2 resha
  • 循环中断打破tqdm

    下面的简单代码使用tqdm https github com tqdm tqdm在循环迭代时显示进度条 import tqdm for f in tqdm tqdm range 100000000 if f gt 100000000 4 b
  • Python - 按月对日期进行分组

    这是一个简单的问题 起初我认为很简单而忽略了它 一个小时过去了 我不太确定 所以 我有一个Python列表datetime对象 我想用图表来表示它们 x 值是年份和月份 y 值是此列表中本月发生的日期对象的数量 也许一个例子可以更好地证明这
  • Python - 在窗口最小化或隐藏时使用 pywinauto 控制窗口

    我正在尝试做的事情 我正在尝试使用 pywinauto 在 python 中创建一个脚本 以在后台自动安装 notepad 隐藏或最小化 notepad 只是一个示例 因为我将编辑它以与其他软件一起使用 Problem 问题是我想在安装程序
  • 如何改变Python中特定打印字母的颜色?

    我正在尝试做一个简短的测验 并且想将错误答案显示为红色 欢迎来到我的测验 您想开始吗 是的 祝你好运 法国的首都是哪里 法国 随机答案不正确的答案 我正在尝试将其显示为红色 我的代码是 print Welcome to my Quiz be
  • 如何将 PIL 图像转换为 NumPy 数组?

    如何转换 PILImage来回转换为 NumPy 数组 这样我就可以比 PIL 进行更快的像素级转换PixelAccess允许 我可以通过以下方式将其转换为 NumPy 数组 pic Image open foo jpg pix numpy
  • 为美国东部以外地区的 Cloudwatch 警报发送短信?

    AWS 似乎没有为美国东部以外的 SNS 主题订阅者提供 SMS 作为协议 我想连接我的 CloudWatch 警报并在发生故障时接收短信 但无法将其发送到 SMS YES 经过一番挖掘后 我能够让它发挥作用 它比仅仅选择一个主题或输入闹钟
  • 在Python中重置生成器对象

    我有一个由多个yield 返回的生成器对象 准备调用该生成器是相当耗时的操作 这就是为什么我想多次重复使用生成器 y FunctionWithYield for x in y print x here must be something t
  • 用于运行可执行文件的python多线程进程

    我正在尝试将一个在 Windows 上运行可执行文件并管理文本输出文件的 python 脚本升级到使用多线程进程的版本 以便我可以利用多个核心 我有四个独立版本的可执行文件 每个线程都知道要访问它们 这部分工作正常 我遇到问题的地方是当它们
  • 对输入求 Keras 模型的导数返回全零

    所以我有一个 Keras 模型 我想将模型的梯度应用于其输入 这就是我所做的 import tensorflow as tf from keras models import Sequential from keras layers imp
  • Python 分析:“‘select.poll’对象的‘poll’方法”是什么?

    我已经使用 python 分析了我的 python 代码cProfile模块并得到以下结果 ncalls tottime percall cumtime percall filename lineno function 13937860 9
  • PyAudio ErrNo 输入溢出 -9981

    我遇到了与用户相同的错误 Python 使用 Pyaudio 以 16000Hz 录制音频时出错 https stackoverflow com questions 12994981 python error audio recording

随机推荐