फास्ट लॉगरिथम गणना

Question

सभी कोड लिनक्स पर एक ही मशीन पर चलाया गया था।

अजगर में:

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

10 छोरों, 3 का सबसे अच्छा: पाश प्रति 77.9 एमएस

C++ (छ साथ ++ ./log.cpp लोग इन -ओ -std = C++ 11 संकलित -O3):

01,235,164: 60 एमएस

MATLAB में

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

रन

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

बीता हुआ समय 7.8 एमएस है।

मैं सेल्सियस के बीच गति अंतर ++ और 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 एमएस, जो बहुत अच्छा है में चलाता है, लेकिन सटीकता में नुकसान महत्वपूर्ण है।

MATLAB को जादुई लॉग 2 गति को प्राप्त करने के तरीके पर कोई विचार है?

अद्यतन

आप सभी टिप्पणियों के लिए धन्यवाद, यह बहुत तेजी से और काफी मददगार था! दरअसल, जवाब समांतरता है, यानी कई धागे पर भार फैलाना। @morningsun सुझाव,

% timeit numexpr.evaluate ('लॉग (DRR)')

5.6 एमएस देता है, जो MATLAB के समान है के बाद, धन्यवाद! numexpr एमकेएल सक्षम है

Answer 1

ध्यान दें कि नीचे सभी फ्लोट 32 हैं, डबल परिशुद्धता नहीं हैं।

अद्यतन: मैंने इंटेल के आईसीसी के पक्ष में पूरी तरह से जीसीसी को हटा दिया है। यह है, जब प्रदर्शन के लिए महत्वपूर्ण है और आप को फ़ाइन-ट्यून करने के लिए अपने "संकलक संकेत" जीसीसी vectorization लागू करने के लिए समय नहीं है जब (देखें, जैसे here)

log_omp.c सभी फर्क नहीं पड़ता

जीसीसी: जीसीसी -ओ log_omp.so -fopenmp log_omp.c -lm -O3 -fPIC कमरा साझा -std = c99

आईसीसी: आईसीसी -ओ log_omp.so -openmp loge_omp।ग -lm -O3 -fPIC कमरा साझा -std = c99 -vec-Report1 -xAVX मैं/opt/इंटेल/संगीतकार/MKL/शामिल

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

अजगर में :

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 

तो numxpr संकलित जीसीसी कोड (खराब) संकलित जीसीसी कोड से बेहतर काम कर रहा है, लेकिन आईसीसी जीतता है।

कुछ संसाधन मैं से उपयोगी और शर्मनाक तरीके से इस्तेमाल किया कोड मिला:

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

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