用一个变量找出大量函数的根 [英] Finding the roots of a large number of functions with one variable

问题描述

我正在使用Python/numpy/scipy编写一个小型射线跟踪器.将曲面建模为二维函数，从而给出法线以上的高度.我将寻找射线与曲面之间的交点的问题简化为寻找具有一个变量的函数的根的问题.功能是连续的和连续可微的.

I'm working with Python/numpy/scipy to write a small ray tracer. Surfaces are modelled as two-dimensional functions giving a height above a normal plane. I reduced the problem of finding the point of intersection between ray and surface to finding the root of a function with one variable. The functions are continuous and continuously differentiable.

有没有一种方法比使用scipy根查找器(甚至可能使用多个进程)简单地遍历所有功能更有效?

Is there a way to do this more efficiently than simply looping over all the functions, using scipy root finders (and maybe using multiple processes)?

这些函数是表示射线的线性函数与表面函数之间的差，并限制在相交平面上.

The functions are the difference between a linear function representing the ray and the surface function, constrained to a plane of intersection.

推荐答案

以下示例显示了在其中计算函数x **(a + 1)-b(均带有不同的a和b)的一百万个副本的根的方法使用二等分方法并行.这里大约需要12秒.

The following example shows calculating the roots for 1 million copies of the function x**(a+1) - b (all with different a and b) in parallel using the bisection method. Takes about ~12 seconds here.

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

基本思想是使用可以在变量向量和定义单个变量的参数的等效向量上求值的函数，简单地对变量向量并行运行寻根器的所有常规步骤组件功能.有条件的替换为masks和numpy.where()的组合.这可以继续进行，直到找到所有根源都达到所需的精度为止，或者直到找到足够的根源值得将其从问题中删除，然后再继续进行一个较小的问题(将那些根排除在外).

The basic idea is to simply run all the usual steps of a root finder in parallel on a vector of variables, using a function that can be evaluated on a vector of variables and equivalent vector(s) of parameters that define the individual component functions. Conditionals are replaced with a combination of masks and numpy.where(). This can continue until all roots have been found to the required precision, or alternately until enough roots have been found that it is worth to remove them from the problem and continue with a smaller problem that excludes those roots.

我选择要解决的功能是任意的，但是如果功能良好，它会有所帮助；在这种情况下，该族中的所有函数都是单调的，并且具有一个正根.此外，对于二等分方法，我们需要猜测给出该函数不同符号的变量，并且在这里也很容易想到这些变量(x0和x1的初始值).

The functions I chose to solve are arbitrary, but it helps if the functions are well-behaved; in this case all functions in the family are monotonic and have exactly one positive root. Additionally, for the bisection method we need guesses for the variable that give different signs of the function, and those happen to be quite easy to come up with here as well (the initial values of x0 and x1).

上面的代码可能使用了最简单的根查找器(二分法)，但是相同的技术可以轻松地应用于Newton-Raphson，Ridder's等.根查找方法中的条件数越少，它越适合于这.但是，您将不得不重新实现所需的任何算法，无法直接使用现有的库根查找器功能.

The above code uses perhaps the simplest root finder (bisection), but the same technique could be easily applied to Newton-Raphson, Ridder's, etc. The fewer conditionals there are in a root finding method, the better suited it is to this. However, you will have to reimplement any algorithm you want, there is no way to use an existing library root finder function directly.

上面的代码段在编写时要牢记清晰，而不是速度.避免重复某些计算，特别是每次迭代仅评估函数一次，而不是3次，因此将计算速度提高了9秒，如下所示:

The above code snippet is written with clarity in mind, not speed. Avoiding the repetition of some calculations, in particular evaluating the function only once per iteration instead of 3 times, speeds this up to 9 seconds, as follows:

...
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 comparison, using scipy.bisect() to find one root at a time takes ~94 seconds:

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)

这篇关于用一个变量找出大量函数的根的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！