하나의 변수로 많은 함수의 근원을 찾는다

Question

나는 작은 광선 추적기를 작성하기 위해 Python/numpy/scipy로 작업하고있다. 곡면은 2 차원 함수로 모델링되어 일반 평면 위에 높이를 제공합니다. 나는 광선과 표면의 교차점을 찾는 문제를 해결하여 하나의 변수로 함수의 근원을 찾아 냈습니다. 기능은 지속적이고 차별화됩니다.

단순히 scipy root finder를 사용하여 모든 기능을 루핑하는 것보다 효율적으로 수행하는 방법이 있습니까?

편집 : 함수는 광선을 나타내는 선형 함수와 교차 평면에 구속 된 표면 함수의 차이입니다.

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()의 조합으로 대체됩니다. 모든 뿌리가 필요한 정밀도로 발견 될 때까지 또는 반복적으로 충분한 뿌리가 발견 될 때까지 문제를 해결하고 그 뿌리를 제외한보다 작은 문제를 계속 수행 할 수 있습니다.

해결하기 위해 선택한 기능은 임의적이지만 함수가 정상적으로 작동하면 도움이됩니다. 이 경우 패밀리의 모든 함수는 단조롭고 정확히 하나의 양수근을가집니다. 또한, 2 분법의 경우 함수의 다른 부호를주는 변수에 대한 추측이 필요합니다. 여기서는 여기서도 쉽게 얻을 수 있습니다 (x0 및 x1의 초기 값).

위의 코드는 아마도 가장 단순한 루트 파인더 (bisection)를 사용하지만 동일한 기법을 Newton-Raphson, Ridder 's 등에 쉽게 적용 할 수 있습니다. 루트 찾기 방법에있는 조건 수가 적을수록 적합합니다 이에. 그러나 원하는 알고리즘을 다시 구현해야하므로 기존의 라이브러리 루트 파인더 함수를 직접 사용할 방법이 없습니다.

위의 코드 스 니펫은 속도가 아닌 명료하게 작성되었습니다. 비교를 위해

... 
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()를 사용하여 다음과 같이 한 번만 대신 3 번 반복 당 기능을 평가하는 특히, 일부 계산의 반복을 방지, 9 초이 속도 한번에 루트가 ~ 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)