Быстрый расчет логарифма

Question

Весь код был запущен на том же компьютере в Linux.

В Python:

import numpy as np 
drr = abs(np.random.randn(100000,50)) 
%timeit np.log2(drr) 

10 петель, лучший из 3: 77,9 мс на петле

В C++ (скомпилирован с G ++ -o войти ./log.cpp -std = C++ 11 -O3):

#include <iostream> 
#include <iomanip> 
#include <string> 
#include <map> 
#include <random> 
#include <ctime> 
int main() 
{ 
std::mt19937 e2(0); 
std::normal_distribution<> dist(0, 1); 
const int n_seq = 100000; 
const int l_seq = 50; 
static double x[n_seq][l_seq]; 
for (int n = 0;n < n_seq; ++n) { 
    for (int k = 0; k < l_seq; ++k) { 
    x[n][k] = abs(dist(e2)); 
    if(x[n][k] <= 0) 
     x[n][k] = 0.1; 
    } 
    } 
clock_t begin = clock(); 

for (int n = 0; n < n_seq; ++n) { 
    for (int k = 0; k < l_seq; ++k) { 
    x[n][k] = std::log2(x[n][k]); 
     } 
    } 
    clock_t end = clock(); 

Запускается в 60 мс

В MATLAB:

abr = abs(randn(100000,50)); 
tic;abr=log2(abr);toc 

Истекшее время составляет 7,8 мс.

Я могу понять разницу в скорости между C++ и numpy, но MATLAB превосходит все. Я столкнулся с http://fastapprox.googlecode.com/svn/trunk/fastapprox/src/fastonebigheader.h , но это делает только плавающий, а не двойной, и я не уверен, как его преобразовать в двойное.

Я также попытался это: http://hackage.haskell.org/package/approximate-0.2.2.1/src/cbits/fast.c который имеет быстрые функции журнала, и при компиляции как Numpy ufunc, работает в 20 мс, что является большим, но потеря точности является значительным.

Любые идеи о том, как достичь магической скорости log2, которую получает MATLAB?

UPDATE

Спасибо всем за комментарии, которые очень быстро и очень полезно! В самом деле, ответ - это параллелизация, т. Е. Распространение нагрузки на несколько потоков. По @morningsun предложению

% timeit numexpr.evaluate ('Журнал (УОБ)')

дает 5,6 мс, которая находится на одном уровне с MATLAB, спасибо! numexpr является MKL включен

Answer 1

Обратите внимание, что ВСЕ ниже - float32, а не двойная точность.

ОБНОВЛЕНИЕ: Я отбросил gcc полностью в пользу ICC Intel. Это делает все различие, когда производительность имеет решающее значение, и когда у вас нет времени, чтобы произвести тонкую настройку «компилятор подсказку» для обеспечения GCC векторизации (см, например, here)

log_omp.c,

GCC: gcc -o log_omp.so -fopenmp log_omp.c -lm -O3 -fPIC -shared -std = c99

ICC: icc -o log_omp.so -openmp loge_omp.с -lm -O3 -fPIC -shared -std = C99 -vec-доклад1 -xAVX -I/Opt/Intel/композитор/мкл/включить

#include <math.h> 
#include "omp.h" 
#include "mkl_vml.h" 

#define restrict __restrict 

inline void log_omp(int m, float * restrict a, float * restrict c); 

void log_omp(int m, float * restrict a, float * restrict c) 
{ 
    int i; 
#pragma omp parallel for default(none) shared(m,a,c) private(i) 
    for (i=0; i<m; i++) { 
     a[i] = log(c[i]); 
    } 
} 

// VML/icc only: 
void log_VML(int m, float * restrict a, float * restrict c) 
{ 
    int i; 
    int split_to = 14; 
    int iter = m/split_to; 
    int additional = m % split_to; 

// vsLn(m, c, a); 
#pragma omp parallel for default(none) shared(m,a,c, additional, iter) private(i) num_threads(split_to) 
    for (i=0;i < (m-additional); i+=iter) 
    vsLog10(iter,c+i,a+i); 
    //vmsLn(iter,c+i,a+i, VML_HA); 

    if (additional > 0) 
    vsLog10(additional, c+m-additional, a+m-additional); 
    //vmsLn(additional, c+m-additional, a+m-additional, VML_HA); 
} 

в Python:

from ctypes import CDLL, c_int, c_void_p 
def log_omp(xs, out): 
    lib = CDLL('./log_omp.so') 
    lib.log_omp.argtypes = [c_int, np.ctypeslib.ndpointer(dtype=np.float32), np.ctypeslib.ndpointer(dtype=np.float32)] 
    lib.log_omp.restype = c_void_p 
    n = xs.shape[0] 
    out = np.empty(n, np.float32) 
    lib.log_omp(n, out, xs) 
    return out 

Cython код (в IPython ноутбука, поэтому %% магия):

%%cython --compile-args=-fopenmp --link-args=-fopenmp 
import numpy as np 
cimport numpy as np 
from libc.math cimport log 

from cython.parallel cimport prange 
import cython 

@cython.boundscheck(False) 
def cylog(np.ndarray[np.float32_t, ndim=1] a not None, 
     np.ndarray[np.float32_t, ndim=1] out=None): 
    if out is None: 
     out = np.empty((a.shape[0]), dtype=a.dtype) 
    cdef Py_ssize_t i 
    with nogil: 
     for i in prange(a.shape[0]): 
      out[i] = log(a[i]) 
    return out 

Тайминги:

numexpr.detect_number_of_cores() // 2 
28 

%env OMP_NUM_THREADS=28 
x = np.abs(np.random.randn(50000000).astype('float32')) 
y = x.copy() 

# GCC 
%timeit log_omp(x, y) 
10 loops, best of 3: 21.6 ms per loop 
# ICC 
%timeit log_omp(x, y) 
100 loops, best of 3: 9.6 ms per loop 
%timeit log_VML(x, y) 
100 loops, best of 3: 10 ms per loop 

%timeit cylog(x, out=y) 
10 loops, best of 3: 21.7 ms per loop 

numexpr.set_num_threads(28) 
%timeit out = numexpr.evaluate('log(x)') 
100 loops, best of 3: 13 ms per loop 

Так что numexpr, кажется, делает лучшую работу, чем (плохо) скомпилированный код gcc, но icc выигрывает.

Некоторые ресурсы, которые я нашел полезным и позорно используется код из:

http://people.duke.edu/~ccc14/sta-663/Optimization_Bakeoff.html

https://gist.github.com/zed/2051661