python中的双重反导计算

Question

我有以下问题。我使用 numpy 函数在 python 中定义了一个函数f。该函数在正实数上是平滑且可积的。我想构造函数的双反导数（假设反导数在 0 处的值和斜率均为 0），以便我可以在任何小于 100 的正实数上对其进行评估。

f 的反导数在x 的定义：

integrate f(s) with s from 0 to x

fx 的双反导数的定义：

integrate (integrate f(t) with t from 0 to s) with s from 0 to x

f 的实际形式并不重要，所以为了方便，我将使用一个简单的形式。但请注意，即使我的示例具有已知的封闭形式，但我的实际函数却没有。

import numpy as np

f = lambda x: np.exp(-x)*x

我的解决方案是使用简单的数值积分将反导数构造为数组：

N = 10000
delta = 100/N

xs = np.linspace(0,100,N+1)

vs = f(xs)
avs = np.cumsum(vs)*delta
aavs = np.cumsum(avs)*delta

这当然有效，但它给了我数组而不是函数。但这不是一个大问题，因为我可以使用样条线插入 aavs 以获得函数并摆脱数组。

from scipy.interpolate import UnivariateSpline

aaf = UnivariateSpline(xs, aavs)

函数aaf 近似为f 的双重反导数。

问题是即使它有效，在我获得我的功能之前还有相当多的开销，而且精度很昂贵。

我的另一个想法是通过样条插值f 并取其反导数，但这会引入对于我想要使用该函数来说太大的数值错误。

有没有更好的方法来做到这一点？更好是指在不牺牲准确性的情况下更快。

编辑：我希望可以使用某种傅里叶变换来避免积分两次。我希望vs 有一些方便的转换，允许将值与xs 逐个相乘并转换回来以获得双重反导数。我玩了一下这个，但我迷路了。

编辑：我发现通过使用梯形规则而不是简单求和，可以大大提高准确性。使用 Simpson 规则应该会进一步提高准确度，但是对于 numpy 数组来说有点繁琐。

编辑：正如@user202729 理所当然地抱怨的那样，这似乎不对。它似乎关闭的原因是因为我跳过了一些细节。我在这里解释了为什么我说的有道理，但这并不影响我的问题。

我的实际目标不是找到f 的双重反导数，而是找到它的变换。我已经跳过了，因为我认为它只会混淆问题。

函数f 在 x 接近 0 或无穷大时呈指数衰减。我通过从 0 开始求和到大约 f 的峰值来最小化积分中的数值误差。这确保了相对误差大致恒定。然后我从一些非常大的 x 的相反方向开始，然后回到峰值。然后我对反微分值做同样的事情。

然后我通过另一个对数值错误敏感的函数来转换aavs。然后我找到错误很大的区域（值剧烈振荡）并删除这些值。最后，我用样条线近似我认为是好的值。

现在，如果我使用样条逼近f，它会引入一个绝对误差，这是一个相当大的区间内的主导项。这被“集成”了两次，最终在aavs 中成为一个相当大的相对错误。然后一旦我改造aavs，我发现'好区域'已经大大缩小了。

编辑：f 的实际形式是我仍在研究的东西。但是，它将是对数正态分布的推广。现在我正在和下面的家人一起玩。

我首先定义正态分布的概括：

def pdf_n(params, center=0.0, slope=8):
    scale, min, diff = params
    if diff > 0:
        r = min
        l = min + diff
    else:
        r = min - diff
        l = min
    def retfun(m):
        x = (m - center)/scale
        E = special.expit(slope*x)*(r - l) + l
        return np.exp( -np.power(1 + x*x, E)/2 )
    return np.vectorize(retfun)

这里发生的事情可能并不明显，但结果很简单。该函数在左侧衰减为 exp(-x^(2l))，在右侧衰减为 exp(-x^(2r))。对于 min=1 和 diff=0，这是正态分布。请注意，这不是标准化的。然后我定义

g = pdf(params)
f = np.vectorize(lambda x:g(np.log(x))/x/area)

其中area 是归一化常数。

请注意，这不是我使用的实际代码。我将其剥离到最低限度。

Answer 1

您可以使用 Numba 更有效地同时计算两个 np.cumsum（和除法）。这明显更快，因为不需要分配、填充、再次读取和释放多个 临时数组。这是一个幼稚的实现：

import numba as nb

@nb.njit('float64[::1](float64[::1], float64)')  # Assume vs is contiguous
def doubleAntiderivative_naive(vs, delta):
    res = np.empty(vs.size, dtype=np.float64)
    sum1, sum2 = 0.0, 0.0
    for i in range(vs.size):
        sum1 += vs[i] * delta
        sum2 += sum1 * delta
        res[i] = sum2
    return res

但是，总和在数值稳定性方面并不是很好。需要Kahan summation 来提高准确性（或者如果您对准确性和性能不那么重要，则可能是替代的 Kahan–Babuška-Klein 算法）。请注意，Numpy 使用成对算法，该算法非常好，但在准确性方面远非完美（这是性能和准确性的良好折衷）。

另外，delta在求和过程中可以因式分解（即结果只需预乘delta**2）。

这是一个使用更准确的 Kahan 求和的实现：

@nb.njit('float64[::1](float64[::1], float64)')
def doubleAntiderivative_accurate(vs, delta):
    res = np.empty(vs.size, dtype=np.float64)
    delta2 = delta * delta
    sum1, sum2 = 0.0, 0.0
    c1, c2 = 0.0, 0.0

    for i in range(vs.size):
        # Kahan summation of the antiderivative of vs
        y1 = vs[i] - c1
        t1 = sum1 + y1
        c1 = (t1 - sum1) - y1
        sum1 = t1

        # Kahan summation of the double antiderivative of vs
        y2 = sum1 - c2
        t2 = sum2 + y2
        c2 = (t2 - sum2) - y2
        sum2 = t2

        res[i] = sum2 * delta2

    return res

这是我的机器（使用 i5-9600KF 处理器）上的方法的性能：

Numpy cumsum:    51.3 us
Naive Numba:     11.6 us
Accutate Numba:  37.2 us

这是方法的相对误差（基于提供的输入函数）：

Numpy cumsum:       1e-13
Naive Numba:        5e-14
Accutate Numba:     2e-16
Perfect precision:  1e-16  (assuming 64-bit numbers are used)

如果 f 可以使用 Numba 轻松计算（这里就是这种情况），那么 vs[i] 可以替换为对 f 的调用（由 Numba 内联）。这有助于减少计算的内存消耗（N 可以很大，但不会使 RAM 饱和）。

至于插值，样条线通常会给出很好的数值结果，但它们的计算成本非常高，而且 AFAIK 需要计算整个数组（数组中的每个项目都会影响所有样条线，尽管有些项目的影响可以忽略不计独自的）。关于您的需求，您可以考虑使用Lagrange polynomials。在边缘上使用拉格朗日多项式时应该小心。在您的情况下，您可以通过使用边界值扩展数组大小来轻松解决边缘上的数值发散问题（因为您知道vs 的每个边缘上的导数为 0）。您可以使用这种方法即时应用插值，这对性能（通常是并行计算）和内存使用都有好处。

Answer 2

首先，我创建了一个我觉得更直观的代码版本。在这里，我将累积总和值乘以 bin 宽度。我相信原始版本的代码中存在与 bin 宽度问题相关的小错误。

import numpy as np
f = lambda x: np.exp(-x)*x
N = 1000
xs = np.linspace(0,100,N+1)

domainwidth = ( np.max(xs) - np.min(xs) )
binwidth = domainwidth / N
vs = f(xs)
avs = np.cumsum(vs)*binwidth
aavs = np.cumsum(avs)*binwidth

接下来，为了可视化，这里有一些非常简单的绘图代码：

import matplotlib
import matplotlib.pyplot as plt
plt.figure()
plt.scatter( xs, vs )
plt.figure()
plt.scatter( xs, avs )
plt.figure()
plt.scatter( xs, aavs )
plt.show()

第一个积分匹配示例表达式的已知结果，可以在wolfram上看到

下面是一个从二阶导数中提取元素的简单函数。请注意，int 是一个糟糕的舍入函数。我假设这是您已经实现的。

def extract_double_antideriv_value(x):
    return aavs[int(x/binwidth)]
singleresult = extract_double_antideriv_value(50.24)
print('singleresult', singleresult)

无论需要什么完整的计算步骤，我们都需要在开始优化之前了解它们。您是否有一百万种不同的功能需要集成？如果您只需要多次查询单个双反导数，您的原始解决方案应该是相当理想的。

符号近似： 您是否考虑过原始函数f 的近似值，它可以具有封闭形式的集成解决方案？您有一个函数所在的有限域。也许近似f 与泰勒级数（可以用known maximum error 构造）然后精确整合？ （考虑 Pade、Taylor、Fourier、Cheby、Lagrange（根据另一个答案的建议）等......）

日志技巧： 处理尖峰错误的另一种方法是记录原始函数的日志。 f 总是积极的吗？是否因为最大值附近的邻域非常小而导致积分错误？如果是这样，您可以改为学习ln(f) 甚至ln(ln(f))。了解f 的样子真的很有帮助。

近似积分技巧 一般来说，存在无数的积分技巧，可以对不可撤销积分做出近似封闭形式的解决方案。当涉及指数函数时（我认为你的函数是指数函数？），一个非常常见的方法是使用Laplace's Method。但是要从袋子里拿出什么技巧在很大程度上取决于f 满足的条件。