1

在数学中，泰勒级数对于获得函数的近似值很重要，具有小次数的多项式。

我想看看这样的近似是如何有帮助的，例如为了加快计算速度。让我们使用著名的泰勒级数：

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

从道德上讲，计算 2 次多项式的值应该比计算 a 快得多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中间。

4

2 回答 2

1

使用 numpy 的矢量化操作几乎总是比在您自己的代码中尝试的任何优化更快。正如@Divakar 所提到的，vectorize这实际上只是编写 for 循环的一种便捷方式，因此您的代码将比 numpy 的本机代码慢。

用标准 python 代码替换 numpy 的优化例程表明您的方法速度大致相同。

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提供了很好的改进，并允许近似值略微优于 pythonlog

1.1818947792053223
0.8402454853057861

编辑：或者，由于这里的重要教训是优化的科学库几乎在一周中的任何一天都将胜过手动编码的解决方案，这是使用 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 的 vectorized 快一个数量级log。

于 2016-11-07T19:47:58.760 回答
1

使用 Python+Numpy，它可能在这里和那里进行了优化，因此不可能真正log(1+x)对x + 0.5 * x^2. 所以我转向了 C++。

结果：

对数的每次操作时间：19.57 ns
对数的 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);

}
于 2016-11-09T15:19:52.530 回答