numpy.amax() will find the max value in an array, and numpy.amin() does the same for the min value. If I want to find both max and min, I have to call both functions, which requires passing over the (very big) array twice, which seems slow.
Is there a function in the numpy API that finds both max and min with only a single pass through the data?
numpy APIに、データを1回通過するだけで最大値と最小値の両方を検出する関数はありますか？
いいえ。この記事の執筆時点では、そのような機能はありません。（はい、そのような関数があれば、そのパフォーマンスは、大きな配列で連続して呼び出すよりも大幅に優れています。）numpy.amin()numpy.amax()
配列を2回渡すことは問題ではないと思います。 次の擬似コードを検討してください。
minval = array[0]
maxval = array[0]
for i in array:
if i < minval:
minval = i
if i > maxval:
maxval = i
ここにはループが1つしかありませんが、それでも2つのチェックがあります。（それぞれ1つのチェックで2つのループを持つ代わりに）。実際に節約できるのは、1ループのオーバーヘッドだけです。あなたが言うように配列が本当に大きい場合、そのオーバーヘッドは実際のループの作業負荷と比較して小さいです。（これはすべてCで実装されているため、ループは多かれ少なかれ自由であることに注意してください）。
編集賛成して私を信じてくれた4人に申し訳ありません。あなたは間違いなくこれを最適化することができます。
これは、PythonモジュールにコンパイルできるいくつかのFortranコードです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。fortranでの単純な実装を比較しても、minnumpyを約2倍上回っています...
LLVMを使用するNumPy対応の動的PythonコンパイラであるNumbaを使用できます。結果として得られる実装は非常に単純で明確です。
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()))
min() & max()また、Numpyの実装よりも高速である必要があります。そして、すべてC/Fortranの1行のコードを記述する必要はありません。
アーキテクチャ、データ、パッケージのバージョンに常に依存するため、独自のパフォーマンステストを実行してください...
それがあなたにとって有用であるならば、numpy.ptpと呼ばれる（max-min）を見つけるための関数があります：
>>> import numpy
>>> x = numpy.array([1,2,3,4,5,6])
>>> x.ptp()
5
しかし、1回のトラバーサルで最小値と最大値の両方を見つける方法はないと思います。
次のアプローチを前提として、予想される数についていくつかのアイデアを得るだけです。
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()最後に、これらのどれも/ほど柔軟ではないことに注意してくださいnp.max()（n-dimサポート、axisパラメーターなどの点で）。
（完全なコードはここにあります）
一般に、一度に2つの要素を処理し、小さい方を一時的な最小値と比較し、大きい方を一時的な最大値と比較するだけで、minmaxアルゴリズムの比較の量を減らすことができます。平均すると、単純なアプローチよりも比較の3/4しか必要ありません。
これは、cまたはfortran（またはその他の低水準言語）で実装でき、パフォーマンスの点でほとんど無敵であるはずです。私はnumbaを使用して原理を説明し、非常に高速でdtypeに依存しない実装を取得しています。
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）
numpy.percentileについては誰も言及していなかったので、そうだと思いました。パーセンタイルを要求する[0, 100]と、最小（0パーセンタイル）と最大（100パーセンタイル）の2つの要素の配列が表示されます。
ただし、OPの目的を満たしていません。最小値と最大値を別々に高速化することはできません。これはおそらく、極端でないパーセンタイルを可能にするいくつかの機械によるものです（より難しい問題であり、より時間がかかるはずです）。
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'
[0, 100]Numpyの将来のバージョンでは、要求された場合にのみ通常のパーセンタイル計算をスキップする特別なケースを設定する可能性があります。インターフェイスに何も追加せずに、1回の呼び出しでNumpyに最小値と最大値を要求する方法があります（受け入れられた回答で述べられたこととは異なります）が、ライブラリの標準実装はこのケースを利用してそれを作成しません価値があります。
これは古いスレッドですが、とにかく、誰かがこれをもう一度見ると...
最小値と最大値を同時に検索する場合、比較の数を減らすことができます。比較しているフロートの場合（おそらくそうだと思います）、計算の複雑さはありませんが、時間を節約できる可能性があります。
（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]
最初に配列内の2つの隣接する値を比較してから、小さい方の値を現在の最小値と比較し、大きい方の値を現在の最大値と比較することができます。
## 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）の比較を行います。それとは対照的に、「明白な方法」で比較を行うと、反復ごとに2つの比較が行われ、2 * len（ar）の比較になります。比較時間を25％節約できます。
多分誰かがいつかこれが役に立つと思うでしょう。
一見すると、トリックを実行しているように見えます。numpy.histogram
count, (amin, amax) = numpy.histogram(a, bins=1)
...しかし、その関数のソースを見ると、それは単に独立して呼び出すためa.min()、a.max()この質問で取り上げたパフォーマンスの懸念を回避できません。:-(
同様に、scipy.ndimage.measurements.extrema可能性のように見えますが、それも単に呼び出しa.min()てa.max()独立しています。
とにかく努力する価値があったので、興味のある人のために、ここで最も難しく、最もエレガントでない解決策を提案します。私の解決策は、C ++のワンパスアルゴリズムでマルチスレッドのmin-maxを実装し、これを使用してPython拡張モジュールを作成することです。この作業には、PythonおよびNumPy C / C ++ APIの使用方法を学習するためのオーバーヘッドが少し必要です。ここでは、コードを示し、このパスをたどりたい人のために簡単な説明とリファレンスを示します。
ここにはあまり興味深いものはありません。配列はサイズのチャンクに分割されますlength / workers。最小/最大は、のチャンクごとに計算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でC++コードを使用する1つの方法は、拡張モジュールを実装することです。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_concurrencyxの増加に近づき始めると予想されます。
一部のNumPyコードには、アプリケーション固有の最適化の余地が確かにあります。特にマルチスレッドに関してはそう思われます。努力する価値があるかどうかは私にはわかりませんが、それは確かに良い運動（または何か）のようです。Cythonのような「サードパーティツール」のいくつかを学ぶことは、時間のより良い使い方かもしれないと思いますが、誰が知っていますか。
前の回答に触発されて、2次元配列からaxis=0のminmaxを返すnumba実装を作成しました。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)
私が思いついた最短の方法はこれです：
mn, mx = np.sort(ar)[[0, -1]]
ただし、配列を並べ替えるため、最も効率的ではありません。
もう1つの短い方法は次のとおりです。
mn, mx = np.percentile(ar, [0, 100])
これはより効率的であるはずですが、結果が計算され、floatが返されます。