快速对数计算

Question

所有代码都在同一台Linux机器上运行。

使用Python：

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

每轮运行10次，3轮取最佳结果：平均每轮耗时77.9毫秒。

在C++中（使用g++编译器，命令为"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();

运行时间为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 但这个只支持float类型，不支持double类型，而且我不确定如何将其转换为double类型。

我还尝试了这个： http://hackage.haskell.org/package/approximate-0.2.2.1/src/cbits/fast.c 它具有快速的log函数，并且编译为numpy ufunc时运行时间为20毫秒，这很好，但精度损失很大。

有什么想法可以达到MATLAB得到的神奇的log2速度吗？

更新

谢谢大家的评论，非常快速和有帮助！实际上，答案是并行化，即将负载分布在多个线程上。根据@morningsun的建议：

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

得到5.6毫秒，与MATLAB相当，谢谢！numexpr是MKL启用的。

- Fermion Portal

---

向量化和并行化。 - IKavanagh

请查看此问题：https://dev59.com/aG025IYBdhLWcg3wW0oc 我猜MATLAB的MAGMA和BLAST包是魔法词。这个问题涉及C++和Java中的循环，而MATLAB在幕后将所有内容向量化。 - Adriaan

Python 模块 numexpr 也值得进行基准测试。如果您启用了 MKL（VML），它可以实时执行 SIMD 和线程。它不能直接执行 log2，所以请使用 log(a)/log(2.0)。 - user2379410

有一件事 - 如果你必须处理巨大的矩阵，而且无法以密集格式放入RAM中，那么numexpr可能不是解决问题的方法。 numexp与scipy稀疏矩阵格式不兼容（在这方面遇到了问题 - http://stackoverflow.com/questions/33824617/numexpr-doesnt-recognize-float-type-sparse-matrix） - El Brutale

@ElBrutale，由于取对数是逐元素进行的操作，因此您只需要稀疏数组中存储的值，这些值可以通过其.data属性轻松访问。 - ali_m

Answer 1

请注意，以下所有内容均为float32，不是双精度。

更新： 我已经完全放弃gcc，转而使用英特尔的icc。当性能至关重要且没有时间微调您的“编译器提示”来强制执行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.c -lm -O3 -fPIC -shared -std=c99 -vec-report1 -xAVX -I/opt/intel/composer/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);
}

在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

- Fermion Portal