Calcolo logaritmo veloce

Question

Tutto il codice è stato eseguito sulla stessa macchina su linux.

In pitone:

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

10 cicli, più di 3: 77,9 ms per loop

In C++ (compilato con g ++ -o log ./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(); 

Viene eseguito in 60 ms

in MATLAB:

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

Il tempo trascorso è 7,8 ms.

Riesco a capire la differenza di velocità tra C++ e numpy, ma MATLAB batte tutto. Mi sono imbattuto in ma questo è solo mobile, non doppio, e non sono sicuro di come convertirlo in doppio.

Ho anche provato questo: http://hackage.haskell.org/package/approximate-0.2.2.1/src/cbits/fast.c che ha funzioni di registro veloci, e quando compilato come ufunc NumPy, viene eseguito in 20 ms, che è grande, ma la perdita di precisione è significativo.

Qualche idea su come ottenere la magica velocità log2 che MATLAB ottiene?

UPDATE

Grazie a tutti per i commenti, che era molto veloce e molto utile! In effetti, la risposta è la parallelizzazione, cioè la diffusione del carico su più thread. Seguendo il suggerimento @morningsun,

% timeit numexpr.evaluate ('log (RRC)')

dà 5.6 ms, che è alla pari con MATLAB, grazie! numexpr è abilitato per MKL

Answer 1

Si noti che TUTTI di seguito sono float32, non doppia precisione.

UPDATE: Ho abbandonato gcc completamente a favore di icc di Intel. Si fa la differenza quando le prestazioni sono critiche e quando non si ha il tempo per mettere a punto i vostri suggerimenti "compilatore" per far rispettare vettorializzazione gcc (vedi, ad esempio, 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.c -lm -O3 -fPIC -shared -std = c99 -vec-relazione1 -xAVX -I/opt/intel/compositore/mkl/include

#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); 
} 

in 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 

codice Cython (in notebook ipython, da qui il %% magia):

%%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 

Tempi:

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 

Quindi numexpr sembra fare un lavoro migliore rispetto al codice gcc (scarsamente) compilato, ma vince icc.

Alcune risorse che ho trovato codice utile e vergognosamente usato da:

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

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