快速对数计算
Fast logarithm calculation
所有代码都在linux上的同一台机器上运行。
在python中:
import numpy as np
drr = abs(np.random.randn(100000,50))
%timeit np.log2(drr)
10个环路,3个最佳:每个环路 77.9 ms
在C++中(用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但这只是浮动的,不是双倍的,我不知道如何将其转换为双倍。
我也尝试过:http://hackage.haskell.org/package/approximate-0.2.2.1/src/cbits/fast.c它具有快速日志功能,当编译为numpy ufunc时,运行时间为20ms,这很好,但准确性损失很大。
关于如何实现MATLAB神奇的log2速度,有什么想法吗?
更新
感谢大家的评论,这非常快,非常有帮助!事实上,答案是并行化,即将负载分散在几个线程上。遵循@morningsun的建议,
%timeit numexpr.evaluate('log(drr)')
给出5.6毫秒,与MATLAB不相上下,谢谢!numexpr是启用MKL的
请注意,下面的ALL都是float32,而不是双精度。
更新:我已经完全抛弃了gcc,转而支持英特尔的icc。当性能至关重要时,以及当你没有时间微调你的"编译器提示"以强制执行gcc矢量化时,这一切都会有所不同(参见,例如这里)
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/comporter/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笔记本中,因此有%%magic):
%%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
- 为什么"do while"循环不断退出,即使条件计算结果为 false?
- 递归函数计算序列中的平方和(并输出过程)
- (C++)分析树以计算返回错误值的简单算术表达式
- 我的字符计数代码计算错误.为什么
- 在计算中使用二的幂有多有利可图
- 如何计算文件中的"columns"数?
- 计算排序向量的向量中唯一值的计数
- 如何使用 std::累积在 C++ 中计算总和立方体
- 使用Qt C++计算类似Git的SHA1哈希
- OpenCV C++.快速计算混淆矩阵
- cpp二进制搜索问题,计算给定数组中输入元素的出现次数
- C++如何计算用户输入的数字中的偶数位数
- 如何计算数据类型的范围,例如int
- 类似枚举的计算常量
- 计算每个节点的树高,帮助我解释这个代码解决方案
- 多个If语句与使用逻辑运算符计算条件的单个语句的比较
- 计算缩放多边形的比例,得到给定的多边形面积
- 在C++中如何在没有pow的情况下进行基础计算
- 计算平均值,不包括上次得分
- 如何计算多映射中重复对的数量