Использование серии Тейлора для ускорения вычислений

Question

В математике Taylor series важны для получения приближений функций с многочленами малой степени.

Я хотел видеть, как такое приближение может быть полезно, например, для ускорения вычислений. Давайте использовать знаменитый ряд Тейлора:

log(1+x) = x + 0.5 * x^2 + (error term) 

Морально вычисления значения многочлена степени 2 должен быть намного быстрее, чем вычисления log.

Таким образом, код, чтобы проверить это:

import numpy, time 

def f(t): 
    return t + 0.5 * t ** 2 
f = numpy.vectorize(f) 

s = time.time() 
for i in range(100): 
    x = numpy.random.rand(100000) 
    numpy.log(1 + x) 
print time.time() - s   # 0.556999921799 seconds 

s = time.time() 
for i in range(100): 
    x = numpy.random.rand(100000) 
    f(x) 
print time.time() - s   # arghh! 4.81500005722 seconds 

Почему полиномиальный метод в 10 раз медленнее, чем фактический журнал? Я ожидал обратного.

PS: Этот вопрос, вероятно, находится посреди SO и math.SE.

Answer 1

С Python + Numpy он, вероятно, оптимизирован здесь и там, и поэтому невозможно реально оценить log(1+x) vs x + 0.5 * x^2. Итак, я переехал на C++.

Результат:

Время на операции с журналом: 19.57 нс
Время на операцию с заказом-2 Тейлора расширения журнала: 3.73 ns

Так грубо x5 фактор!

---

#include <iostream> 
#include <math.h> 
#include <time.h> 
#define N (1000*1000*100) 
#define NANO (1000*1000*1000) 

int main() 
{ 
    float *x = (float*) malloc(N * sizeof(float)); 
    float y; 
    float elapsed1, elapsed2; 
    clock_t begin, end; 
    int i; 

    for (i = 0; i < N; i++) 
    x[i] = (float) (rand() + 1)/(float)(RAND_MAX); 

    begin = clock(); 
    for (i = 0; i < N; i++) 
    y = logf(x[i]); 
    end = clock(); 
    elapsed1 = float(end - begin)/CLOCKS_PER_SEC/N * NANO; 

    begin = clock(); 
    for (i = 0; i < N; i++) 
    y = x[i] + 0.5 * x[i] * x[i]; 
    end = clock(); 
    elapsed2 = float(end - begin)/CLOCKS_PER_SEC/N * NANO; 

    std::cout << "Time per operation with log: " << elapsed1 << " ns\n"; 
    std::cout << "Time per operation with order-2 Taylor epansion: " << elapsed2 << " ns"; 

    free(x); 

} 

Answer 2

Использование векторизованных операций numpy почти всегда будет быстрее любых попыток оптимизации в вашем собственном коде. Как отметил @Divakar, vectorize - это действительно удобный способ написания цикла for, поэтому ваш код будет медленнее, чем собственный код numpy.

Замена оптимизированных подпрограмм numpy стандартным кодом python показывает, что ваш метод имеет одинаковую скорость.

import math, numpy, time 


def f(t): 
    return t + 0.5 * t ** 2 

x = numpy.random.rand(1000000) 

s = time.time() 
for num in x: 
    math.log(1 + num) 
print (time.time() - s ) 

s = time.time() 
for num in x: 
    f(num) 
print (time.time() - s)  

результат:

1.1951053142547607 
1.3485901355743408 

приближение лишь немного медленнее, но экспоненцирование очень дорого. Замена t ** 2 с t*t дает хорошее улучшение, и позволяет приближение немного опережать Питон log

1.1818947792053223 
0.8402454853057861 

Edit: В качестве альтернативы, так как большой урока здесь Оптимизированные научные библиотеки будут производительнее handcoded решения практически в любой день недели, вот приближение серии Тейлора с векторизованными операциями numpy, что на сегодняшний день является самым быстрым. Обратите внимание, что единственным большим изменением является то, что vectorize не вызывается для функции аппроксимации, поэтому по умолчанию используются векторизованные операции numpy.

import numpy, time 

def f(t): 
    return t + 0.5 * t ** 2 

x = numpy.random.rand(1000000) 
s = time.time() 
numpy.log(1 + x) 
print (time.time() - s) 

s = time.time() 
x = numpy.random.rand(100000) 
f(x) 
print (time.time() - s ) 

результат:

0.07202601432800293 
0.0019881725311279297 

Там у вас есть, векторизованное приближение более чем на порядок величины быстрее, чем NumPy это векторизация log.