NumPy:同时实现 max() 和 min() 的函数
- 2025-04-16 08:58:00
- admin 原创
- 20
问题描述:
numpy.amax()会查找数组中的最大值,而numpy.amin()则会查找最小值。如果我想同时查找最大值和最小值,则必须调用这两个函数,这需要两次传递(非常大的)数组,这似乎很慢。
numpy API 中是否有一个函数,只需传递一次数据即可找到最大值和最小值?
解决方案 1:
numpy API 中是否有一个函数,只需传递一次数据即可找到最大值和最小值?
没有。截至撰写本文时,还没有这样的函数。(当然,如果有这样的函数,其性能会比在一个大型数组上依次调用和好得多。)numpy.amin()
`numpy.amax()`
解决方案 2:
你可以使用Numba,它是一个使用 LLVM 且支持 NumPy 的动态 Python 编译器。最终的实现非常简单明了:
import numpy
import numba
@numba.jit
def minmax(x):
maximum = x[0]
minimum = x[0]
for i in x[1:]:
if i > maximum:
maximum = i
elif i < minimum:
minimum = i
return (minimum, maximum)
numpy.random.seed(1)
x = numpy.random.rand(1000000)
print(minmax(x) == (x.min(), x.max()))
它也应该比 Numpy 的min() & max()
实现更快。而且无需编写任何 C/Fortran 代码。
做你自己的性能测试,因为它总是依赖于你的架构、你的数据、你的包版本……
解决方案 3:
我认为两次传递数组不会有问题。 考虑以下伪代码:
minval = array[0]
maxval = array[0]
for i in array:
if i < minval:
minval = i
if i > maxval:
maxval = i
虽然这里只有 1 个循环,但仍然有 2 个检查。(而不是 2 个循环,每个循环检查 1 次)。实际上,您唯一节省的就是 1 个循环的开销。如果数组真的像您说的那样大,那么与实际循环的工作负载相比,这部分开销就很小了。(请注意,这都是用 C 语言实现的,所以循环本身或多或少是空闲的。)
编辑:对于点赞并信任我的四位朋友,我深感抱歉。你们肯定可以优化这一点。
这里有一些 Fortran 代码,可以通过以下方式编译成 Python 模块f2py
(也许有位Cython
大师可以来将它与优化的 C 版本进行比较……):
subroutine minmax1(a,n,amin,amax)
implicit none
!f2py intent(hidden) :: n
!f2py intent(out) :: amin,amax
!f2py intent(in) :: a
integer n
real a(n),amin,amax
integer i
amin = a(1)
amax = a(1)
do i=2, n
if(a(i) > amax)then
amax = a(i)
elseif(a(i) < amin) then
amin = a(i)
endif
enddo
end subroutine minmax1
subroutine minmax2(a,n,amin,amax)
implicit none
!f2py intent(hidden) :: n
!f2py intent(out) :: amin,amax
!f2py intent(in) :: a
integer n
real a(n),amin,amax
amin = minval(a)
amax = maxval(a)
end subroutine minmax2
通过以下方式编译:
f2py -m untitled -c fortran_code.f90
现在我们可以测试它了:
import timeit
size = 100000
repeat = 10000
print timeit.timeit(
'np.min(a); np.max(a)',
setup='import numpy as np; a = np.arange(%d, dtype=np.float32)' % size,
number=repeat), " # numpy min/max"
print timeit.timeit(
'untitled.minmax1(a)',
setup='import numpy as np; import untitled; a = np.arange(%d, dtype=np.float32)' % size,
number=repeat), '# minmax1'
print timeit.timeit(
'untitled.minmax2(a)',
setup='import numpy as np; import untitled; a = np.arange(%d, dtype=np.float32)' % size,
number=repeat), '# minmax2'
结果让我有点震惊:
8.61869883537 # numpy min/max
1.60417699814 # minmax1
2.30169081688 # minmax2
我不得不说,我不太明白。仅仅np.min
比较minmax1
和minmax2
仍然是徒劳的,所以这不仅仅是一个记忆问题……
注释——将规模增加 1 倍10**a
,将重复次数减少 1 倍10**a
(保持问题规模不变)确实会改变性能,但这种改变似乎并不一致,这表明 Python 中的内存性能和函数调用开销之间存在一些相互作用。即使min
与 Fortran 中的简单实现进行比较,其性能也比 NumPy 的要好大约 2 倍……
解决方案 4:
如果这对您有用的话,有一个名为numpy.ptp的函数可以用于查找(最大值-最小值):
>>> import numpy
>>> x = numpy.array([1,2,3,4,5,6])
>>> x.ptp()
5
但我不认为有办法通过一次遍历找到最小值和最大值。
编辑: ptp 只是在后台调用 min 和 max
解决方案 5:
仅就以下方法得出一些关于预期数字的想法:
import numpy as np
def extrema_np(arr):
return np.max(arr), np.min(arr)
import numba as nb
@nb.jit(nopython=True)
def extrema_loop_nb(arr):
n = arr.size
max_val = min_val = arr[0]
for i in range(1, n):
item = arr[i]
if item > max_val:
max_val = item
elif item < min_val:
min_val = item
return max_val, min_val
import numba as nb
@nb.jit(nopython=True)
def extrema_while_nb(arr):
n = arr.size
odd = n % 2
if not odd:
n -= 1
max_val = min_val = arr[0]
i = 1
while i < n:
x = arr[i]
y = arr[i + 1]
if x > y:
x, y = y, x
min_val = min(x, min_val)
max_val = max(y, max_val)
i += 2
if not odd:
x = arr[n]
min_val = min(x, min_val)
max_val = max(x, max_val)
return max_val, min_val
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True
import numpy as np
cdef void _extrema_loop_cy(
long[:] arr,
size_t n,
long[:] result):
cdef size_t i
cdef long item, max_val, min_val
max_val = arr[0]
min_val = arr[0]
for i in range(1, n):
item = arr[i]
if item > max_val:
max_val = item
elif item < min_val:
min_val = item
result[0] = max_val
result[1] = min_val
def extrema_loop_cy(arr):
result = np.zeros(2, dtype=arr.dtype)
_extrema_loop_cy(arr, arr.size, result)
return result[0], result[1]
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True
import numpy as np
cdef void _extrema_while_cy(
long[:] arr,
size_t n,
long[:] result):
cdef size_t i, odd
cdef long x, y, max_val, min_val
max_val = arr[0]
min_val = arr[0]
odd = n % 2
if not odd:
n -= 1
max_val = min_val = arr[0]
i = 1
while i < n:
x = arr[i]
y = arr[i + 1]
if x > y:
x, y = y, x
min_val = min(x, min_val)
max_val = max(y, max_val)
i += 2
if not odd:
x = arr[n]
min_val = min(x, min_val)
max_val = max(x, max_val)
result[0] = max_val
result[1] = min_val
def extrema_while_cy(arr):
result = np.zeros(2, dtype=arr.dtype)
_extrema_while_cy(arr, arr.size, result)
return result[0], result[1]
(这些方法与这里extrema_loop_*()
提出的方法类似,而方法基于这里的代码)extrema_while_*()
时间安排如下:
表明 是extrema_while_*()
最快的, 是extrema_while_nb()
最快的。无论如何,extrema_loop_nb()
和 的extrema_loop_cy()
解决方案也确实优于仅使用 NumPy 的方法(分别使用np.max()
和np.min()
)。
最后,请注意,这些都不像np.min()
/那样灵活(在 n-dim 支持、参数等np.max()
方面)。axis
(完整代码可在此处获取)
解决方案 6:
没人提到numpy.percentile,所以我想试试。如果你求[0, 100]
百分位数,它会返回一个包含两个元素的数组,分别是最小值(0 百分位数)和最大值(100 百分位数)。
然而,它并没有满足原帖作者的目的:它的速度并不比分别使用 min 和 max 更快。这可能是因为某些机制允许处理非极端百分位数(这是一个更难的问题,需要更长的时间)。
In [1]: import numpy
In [2]: a = numpy.random.normal(0, 1, 1000000)
In [3]: %%timeit
...: lo, hi = numpy.amin(a), numpy.amax(a)
...:
100 loops, best of 3: 4.08 ms per loop
In [4]: %%timeit
...: lo, hi = numpy.percentile(a, [0, 100])
...:
100 loops, best of 3: 17.2 ms per loop
In [5]: numpy.__version__
Out[5]: '1.14.4'
Numpy 的未来版本可能会添加一个特殊情况,以便在仅[0, 100]
需要时跳过正常的百分位数计算。无需在接口中添加任何内容,就可以在一次调用中向 Numpy 询问最小值和最大值(与已接受答案中所述相反),但该库的标准实现并未利用这种情况来使其变得有价值。
解决方案 7:
一般来说,可以通过一次处理两个元素,并仅将较小的元素与临时最小值进行比较,将较大的元素与临时最大值进行比较,来减少 minmax 算法的比较次数。平均而言,与简单方法相比,这种方法只需要 3/4 的比较次数。
这可以用 C 或 Fortran(或任何其他低级语言)实现,并且性能几乎无与伦比。我正在使用numba为了说明原理并获得非常快速、与数据类型无关的实现:
import numba as nb
import numpy as np
@nb.njit
def minmax(array):
# Ravel the array and return early if it's empty
array = array.ravel()
length = array.size
if not length:
return
# We want to process two elements at once so we need
# an even sized array, but we preprocess the first and
# start with the second element, so we want it "odd"
odd = length % 2
if not odd:
length -= 1
# Initialize min and max with the first item
minimum = maximum = array[0]
i = 1
while i < length:
# Get the next two items and swap them if necessary
x = array[i]
y = array[i+1]
if x > y:
x, y = y, x
# Compare the min with the smaller one and the max
# with the bigger one
minimum = min(x, minimum)
maximum = max(y, maximum)
i += 2
# If we had an even sized array we need to compare the
# one remaining item too.
if not odd:
x = array[length]
minimum = min(x, minimum)
maximum = max(x, maximum)
return minimum, maximum
它绝对比Peque提出的简单方法要快:
arr = np.random.random(3000000)
assert minmax(arr) == minmax_peque(arr) # warmup and making sure they are identical
%timeit minmax(arr) # 100 loops, best of 3: 2.1 ms per loop
%timeit minmax_peque(arr) # 100 loops, best of 3: 2.75 ms per loop
正如预期的那样,新的 minmax 实现仅花费原始实现所需时间的大约 3/4(2.1 / 2.75 = 0.7636363636363637
)
解决方案 8:
这是一个老话题了,但无论如何,如果有人再次看到这个……
当同时查找最小值和最大值时,可以减少比较次数。如果你比较的是浮点数(我猜是的),这可能会节省一些时间,尽管计算复杂度不会降低。
而不是(Python代码):
_max = ar[0]
_min= ar[0]
for ii in xrange(len(ar)):
if _max > ar[ii]: _max = ar[ii]
if _min < ar[ii]: _min = ar[ii]
您可以先比较数组中的两个相邻值,然后仅将较小的值与当前最小值进行比较,将较大的值与当前最大值进行比较:
## for an even-sized array
_max = ar[0]
_min = ar[0]
for ii in xrange(0, len(ar), 2)): ## iterate over every other value in the array
f1 = ar[ii]
f2 = ar[ii+1]
if (f1 < f2):
if f1 < _min: _min = f1
if f2 > _max: _max = f2
else:
if f2 < _min: _min = f2
if f1 > _max: _max = f1
这里的代码是用 Python 编写的,显然为了提高速度,你会使用 C、Fortran 或 Cython,但这样每次迭代你都会进行 3 次比较,迭代次数为 len(ar)/2,因此比较次数为 3/2 len(ar)。与此相反,用“显而易见的方式”进行比较,每次迭代你都会进行两次比较,因此比较次数为 2len(ar)。这样可以节省 25% 的比较时间。
也许有一天有人会发现这很有用。
解决方案 9:
乍一看,似乎可以达到这样的效果:numpy.histogram
count, (amin, amax) = numpy.histogram(a, bins=1)
...但是如果你查看该函数的源代码a.min()
,它只是独立调用a.max()
,因此无法避免该问题中解决的性能问题。:-(
同样,scipy.ndimage.measurements.extrema
看起来像是一种可能性,但它也只是独立地调用a.min()
和a.max()
。
解决方案 10:
无论如何,对我来说,这一切都是值得的,所以我将在这里为感兴趣的人提出一个最困难、最不优雅的解决方案。我的解决方案是用 C++ 实现一个多线程的单遍最小最大算法,并用它来创建一个 Python 扩展模块。这项工作需要一些学习 Python 和 NumPy C/C++ API 的开销,在这里我将展示代码,并为想要尝试这条路线的人提供一些简单的解释和参考。
多线程最小/最大
这里没什么特别的。数组被分成大小为 的块length / workers
。计算 a 中每个块的最小值/最大值future
,然后扫描这些块以查找全局最小值/最大值。
// mt_np.cc
//
// multi-threaded min/max algorithm
#include <algorithm>
#include <future>
#include <vector>
namespace mt_np {
/*
* Get {min,max} in interval [begin,end)
*/
template <typename T> std::pair<T, T> min_max(T *begin, T *end) {
T min{*begin};
T max{*begin};
while (++begin < end) {
if (*begin < min) {
min = *begin;
continue;
} else if (*begin > max) {
max = *begin;
}
}
return {min, max};
}
/*
* get {min,max} in interval [begin,end) using #workers for concurrency
*/
template <typename T>
std::pair<T, T> min_max_mt(T *begin, T *end, int workers) {
const long int chunk_size = std::max((end - begin) / workers, 1l);
std::vector<std::future<std::pair<T, T>>> min_maxes;
// fire up the workers
while (begin < end) {
T *next = std::min(end, begin + chunk_size);
min_maxes.push_back(std::async(min_max<T>, begin, next));
begin = next;
}
// retrieve the results
auto min_max_it = min_maxes.begin();
auto v{min_max_it->get()};
T min{v.first};
T max{v.second};
while (++min_max_it != min_maxes.end()) {
v = min_max_it->get();
min = std::min(min, v.first);
max = std::max(max, v.second);
}
return {min, max};
}
}; // namespace mt_np
Python扩展模块
事情开始变得糟糕起来……在 Python 中使用 C++ 代码的一种方法是实现扩展模块。该模块可以使用distutils.core
标准模块构建和安装。Python 文档中详细介绍了扩展模块的具体功能:https://docs.python.org/3/extending/extending.html。 注意:当然还有其他方法可以获得类似的结果,例如:https://docs.python.org/3/extending/index.html#extending-index:
本指南仅涵盖此版本 CPython 提供的用于创建扩展的基本工具。Cython、cffi、SWIG 和 Numba 等第三方工具提供了更简单和更复杂的方法来为 Python 创建 C 和 C++ 扩展。
本质上,这条路线可能更偏向学术性而非实用性。话虽如此,我接下来做的是,严格按照教程操作,创建一个模块文件。这本质上是 distutils 的样板,它知道如何处理你的代码,并用它来创建一个 Python 模块。在执行任何操作之前,最好先创建一个 Python虚拟环境,这样就不会污染你的系统包(参见https://docs.python.org/3/library/venv.html#module-venv)。
这是模块文件:
// mt_np_forpy.cc
//
// C++ module implementation for multi-threaded min/max for np
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <python3.6/numpy/arrayobject.h>
#include "mt_np.h"
#include <cstdint>
#include <iostream>
using namespace std;
/*
* check:
* shape
* stride
* data_type
* byteorder
* alignment
*/
static bool check_array(PyArrayObject *arr) {
if (PyArray_NDIM(arr) != 1) {
PyErr_SetString(PyExc_RuntimeError, "Wrong shape, require (1,n)");
return false;
}
if (PyArray_STRIDES(arr)[0] != 8) {
PyErr_SetString(PyExc_RuntimeError, "Expected stride of 8");
return false;
}
PyArray_Descr *descr = PyArray_DESCR(arr);
if (descr->type != NPY_LONGLTR && descr->type != NPY_DOUBLELTR) {
PyErr_SetString(PyExc_RuntimeError, "Wrong type, require l or d");
return false;
}
if (descr->byteorder != '=') {
PyErr_SetString(PyExc_RuntimeError, "Expected native byteorder");
return false;
}
if (descr->alignment != 8) {
cerr << "alignment: " << descr->alignment << endl;
PyErr_SetString(PyExc_RuntimeError, "Require proper alignement");
return false;
}
return true;
}
template <typename T>
static PyObject *mt_np_minmax_dispatch(PyArrayObject *arr) {
npy_intp size = PyArray_SHAPE(arr)[0];
T *begin = (T *)PyArray_DATA(arr);
auto minmax =
mt_np::min_max_mt(begin, begin + size, thread::hardware_concurrency());
return Py_BuildValue("(L,L)", minmax.first, minmax.second);
}
static PyObject *mt_np_minmax(PyObject *self, PyObject *args) {
PyArrayObject *arr;
if (!PyArg_ParseTuple(args, "O", &arr))
return NULL;
if (!check_array(arr))
return NULL;
switch (PyArray_DESCR(arr)->type) {
case NPY_LONGLTR: {
return mt_np_minmax_dispatch<int64_t>(arr);
} break;
case NPY_DOUBLELTR: {
return mt_np_minmax_dispatch<double>(arr);
} break;
default: {
PyErr_SetString(PyExc_RuntimeError, "Unknown error");
return NULL;
}
}
}
static PyObject *get_concurrency(PyObject *self, PyObject *args) {
return Py_BuildValue("I", thread::hardware_concurrency());
}
static PyMethodDef mt_np_Methods[] = {
{"mt_np_minmax", mt_np_minmax, METH_VARARGS, "multi-threaded np min/max"},
{"get_concurrency", get_concurrency, METH_VARARGS,
"retrieve thread::hardware_concurrency()"},
{NULL, NULL, 0, NULL} /* sentinel */
};
static struct PyModuleDef mt_np_module = {PyModuleDef_HEAD_INIT, "mt_np", NULL,
-1, mt_np_Methods};
PyMODINIT_FUNC PyInit_mt_np() { return PyModule_Create(&mt_np_module); }
在此文件中,Python 和 NumPy API 得到了大量使用,有关更多信息,请参阅:https://docs.python.org/3/c-api/arg.html#c.PyArg_ParseTuple,有关 NumPy 的信息,请参阅:https://docs.scipy.org/doc/numpy/reference/c-api.array.html。
安装模块
接下来要做的是利用 distutils 安装模块。这需要一个安装文件:
# setup.py
from distutils.core import setup,Extension
module = Extension('mt_np', sources = ['mt_np_module.cc'])
setup (name = 'mt_np',
version = '1.0',
description = 'multi-threaded min/max for np arrays',
ext_modules = [module])
要最终安装模块,python3 setup.py install
请从虚拟环境中执行。
测试模块
最后,我们可以测试一下 C++ 实现是否真的比 NumPy 的简单使用效果更好。为此,这里有一个简单的测试脚本:
# timing.py
# compare numpy min/max vs multi-threaded min/max
import numpy as np
import mt_np
import timeit
def normal_min_max(X):
return (np.min(X),np.max(X))
print(mt_np.get_concurrency())
for ssize in np.logspace(3,8,6):
size = int(ssize)
print('********************')
print('sample size:', size)
print('********************')
samples = np.random.normal(0,50,(2,size))
for sample in samples:
print('np:', timeit.timeit('normal_min_max(sample)',
globals=globals(),number=10))
print('mt:', timeit.timeit('mt_np.mt_np_minmax(sample)',
globals=globals(),number=10))
以下是我完成这些工作后得到的结果:
8
********************
sample size: 1000
********************
np: 0.00012079699808964506
mt: 0.002468645994667895
np: 0.00011947099847020581
mt: 0.0020772050047526136
********************
sample size: 10000
********************
np: 0.00024697799381101504
mt: 0.002037393998762127
np: 0.0002713389985729009
mt: 0.0020942929986631498
********************
sample size: 100000
********************
np: 0.0007130410012905486
mt: 0.0019842900001094677
np: 0.0007540129954577424
mt: 0.0029724110063398257
********************
sample size: 1000000
********************
np: 0.0094779249993735
mt: 0.007134920000680722
np: 0.009129883001151029
mt: 0.012836456997320056
********************
sample size: 10000000
********************
np: 0.09471094200125663
mt: 0.0453535050037317
np: 0.09436299200024223
mt: 0.04188535599678289
********************
sample size: 100000000
********************
np: 0.9537652180006262
mt: 0.3957935369980987
np: 0.9624398809974082
mt: 0.4019058070043684
这些结果远不如之前线程中显示的结果那么鼓舞人心,之前线程中显示的速度大约是 3.5 倍,而且没有使用多线程。我得到的结果还算合理,我预计线程的开销会占据主要时间,直到数组变得非常大,此时性能提升将开始接近std::thread::hardware_concurrency
10 倍。
结论
看来,某些 NumPy 代码确实有针对特定应用进行优化的空间,尤其是在多线程方面。我不清楚这是否值得付出努力,但这确实看起来是一个很好的练习(或者说某种方式)。我觉得学习一些像 Cython 这样的“第三方工具”或许能更好地利用时间,但谁知道呢。
解决方案 11:
受上一个答案的启发,我编写了一个numba实现,从二维数组中返回axis=0时的minmax。它比调用numpy min/max快约5倍。也许有人会觉得它有用。
from numba import jit
@jit
def minmax(x):
"""Return minimum and maximum from 2D array for axis=0."""
m, n = len(x), len(x[0])
mi, ma = np.empty(n), np.empty(n)
mi[:] = ma[:] = x[0]
for i in range(1, m):
for j in range(n):
if x[i, j]>ma[j]: ma[j] = x[i, j]
elif x[i, j]<mi[j]: mi[j] = x[i, j]
return mi, ma
x = np.random.normal(size=(256, 11))
mi, ma = minmax(x)
np.all(mi == x.min(axis=0)), np.all(ma == x.max(axis=0))
# (True, True)
%timeit x.min(axis=0), x.max(axis=0)
# 15.9 µs ± 9.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit minmax(x)
# 2.62 µs ± 31.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
解决方案 12:
我想到的最短方法是这样的:
mn, mx = np.sort(ar)[[0, -1]]
但由于它对数组进行排序,所以它不是最有效的。
另一种简便方法是:
mn, mx = np.percentile(ar, [0, 100])
这应该更有效率,但结果被计算出来,并返回一个浮点数。
解决方案 13:
我编写了一个小型 Python 包,用于对 float32 数组执行此操作。与 np.amax 和 np.amin 相比,速度提高了约 2.3 倍。
用法:
pip install numpy-minmax
import numpy_minmax
min_val, max_val = numpy_minmax.minmax(arr)
该算法用 C 语言编写,并使用 SIMD 指令进行了优化。代码仓库位于:https://github.com/nomonosound/numpy-minmax
解决方案 14:
也许可以使用numpy.unique
? 像这样:
min_, max_ = numpy.unique(arr)[[0, -1]]
只是为了多样性而将它添加到这里:)它和排序一样慢。
扫码咨询,免费领取项目管理大礼包!