使用一个变量查找大量函数的根

Question

我正在使用Python/numpy/scipy编写一个小型光线跟踪器.曲面被建模为二维函数,其高度高于法线平面.我减少了找到光线和曲面之间的交点以找到具有一个变量的函数的根的问题.功能是连续的,可连续区分的.

有没有办法比使用scipy根查找器(并且可能使用多个进程)简单地循环遍历所有函数更有效地执行此操作？

编辑:函数是表示光线的线性函数和表面函数之间的差异,约束到交叉平面.

Answer 1

以下示例显示使用二分法并行计算函数 x**(a+1) - b（均具有不同的 a 和 b）的 100 万个副本的根。这里大约需要 12 秒。

import numpy

def F(x, a, b):
    return numpy.power(x, a+1.0) - b

N = 1000000

a = numpy.random.rand(N)
b = numpy.random.rand(N)

x0 = numpy.zeros(N)
x1 = numpy.ones(N) * 1000.0

max_step = 100
for step in range(max_step):
    x_mid = (x0 + x1)/2.0
    F0 = F(x0, a, b)
    F1 = F(x1, a, b)
    F_mid = F(x_mid, a, b)
    x0 = numpy.where( numpy.sign(F_mid) == numpy.sign(F0), x_mid, x0 )
    x1 = numpy.where( numpy.sign(F_mid) == numpy.sign(F1), x_mid, x1 )
    error_max = numpy.amax(numpy.abs(x1 - x0))
    print "step=%d error max=%f" % (step, error_max)
    if error_max < 1e-6: break

基本思想是简单地在变量向量上并行运行求根器的所有常用步骤，使用可以在变量向量和定义各个分量函数的等效参数向量上进行评估的函数。条件语句被替换为掩码和 numpy.where() 的组合。这可以继续下去，直到找到所需精度的所有根，或者直到找到足够的根，值得将它们从问题中删除，并继续排​​除这些根的较小问题。

我选择求解的函数是任意的，但如果函数表现良好就会有所帮助；在这种情况下，族中的所有函数都是单调的并且只有一个正根。此外，对于二分法，我们需要猜测给出函数不同符号的变量，而这些猜测在这里也很容易得出（x0 和 x1 的初始值）。

上面的代码可能使用了最简单的求根器（二分法），但是相同的技术可以很容易地应用于 Newton-Raphson、Ridder 等。求根方法中的条件越少，就越适合这种情况。但是，您将必须重新实现您想要的任何算法，无法直接使用现有的库根查找函数。

上面的代码片段是为了清晰而写的，而不是为了速度。避免重复某些计算，特别是每次迭代仅评估函数一次而不是 3 次，可将计算速度加快至 9 秒，如下所示：

...
F0 = F(x0, a, b)
F1 = F(x1, a, b)

max_step = 100
for step in range(max_step):
    x_mid = (x0 + x1)/2.0
    F_mid = F(x_mid, a, b)
    mask0 = numpy.sign(F_mid) == numpy.sign(F0)
    mask1 = numpy.sign(F_mid) == numpy.sign(F1)
    x0 = numpy.where( mask0, x_mid, x0 )
    x1 = numpy.where( mask1, x_mid, x1 )
    F0 = numpy.where( mask0, F_mid, F0 )
    F1 = numpy.where( mask1, F_mid, F1 )
...

作为比较，使用 scipy.bisect() 一次查找一个根大约需要 94 秒：

for i in range(N):
    x_root = scipy.optimize.bisect(lambda x: F(x, a[i], b[i]), x0[i], x1[i], xtol=1e-6)