numpy.amax () trouvera la valeur max dans un tableau et numpy.amin () fait de même pour la valeur min. Si je veux trouver à la fois max et min, je dois appeler les deux fonctions, ce qui nécessite de passer deux fois sur le (très grand) tableau, ce qui semble lent.

Y a-t-il une fonction dans l'API numpy qui trouve à la fois max et min en un seul passage dans les données?

Y a-t-il une fonction dans l'API numpy qui trouve à la fois max et min en un seul passage dans les données?

Non. Au moment d'écrire ces lignes, une telle fonction n'existe pas. (Et oui, s'il y avait une telle fonction, ses performances seraient nettement meilleures que celles d'appeler numpy.amin()et numpy.amax()successivement sur un grand tableau.)

Je ne pense pas que passer deux fois sur le tableau soit un problème. Considérez le pseudo-code suivant:

minval = array[0]
maxval = array[0]
for i in array:
    if i < minval:
       minval = i
    if i > maxval:
       maxval = i

Bien qu'il n'y ait qu'une seule boucle ici, il y a encore 2 vérifications. (Au lieu d'avoir 2 boucles avec 1 chèque chacune). En réalité, la seule chose que vous économisez est la surcharge d'une boucle. Si les tableaux sont vraiment gros comme vous le dites, cette surcharge est faible par rapport à la charge de travail réelle de la boucle. (Notez que tout cela est implémenté en C, donc les boucles sont de toute façon plus ou moins libres).

EDIT Désolé pour les 4 d'entre vous qui ont voté pour et ont eu confiance en moi. Vous pouvez certainement optimiser cela.

Voici un code fortran qui peut être compilé dans un module python via f2py(peut-être qu'un Cythongourou peut venir le comparer avec une version C optimisée ...):

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

Compilez-le via:

f2py -m untitled -c fortran_code.f90

Et maintenant, nous sommes dans un endroit où nous pouvons le tester:

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'

Les résultats sont un peu stupéfiants pour moi:

8.61869883537 # numpy min/max
1.60417699814 # minmax1
2.30169081688 # minmax2

Je dois dire que je ne comprends pas complètement. Comparer juste np.mincontre minmax1et minmax2est toujours une bataille perdue, donc ce n'est pas seulement un problème de mémoire ...

notes - L'augmentation de la taille d'un facteur de 10**aet la diminution de la répétition d'un facteur de 10**a(en gardant la taille du problème constante) modifie les performances, mais pas d'une manière apparemment cohérente, ce qui montre qu'il existe une interaction entre les performances de la mémoire et la surcharge des appels de fonction dans python. Même en comparant une minimplémentation simple dans fortran bat numpy par un facteur d'environ 2 ...

Il existe une fonction de recherche (max-min) appelée numpy.ptp si cela vous est utile:

>>> import numpy
>>> x = numpy.array([1,2,3,4,5,6])
>>> x.ptp()
5

mais je ne pense pas qu'il y ait un moyen de trouver à la fois min et max avec un seul parcours.

EDIT: ptp appelle juste min et max sous le capot

Vous pouvez utiliser Numba , qui est un compilateur Python dynamique prenant en charge NumPy utilisant LLVM. L'implémentation qui en résulte est assez simple et claire:

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

Il devrait également être plus rapide que l' min() & max()implémentation d'un Numpy . Et tout cela sans avoir à écrire une seule ligne de code C / Fortran.

Faites vos propres tests de performances, car cela dépend toujours de votre architecture, de vos données, de vos versions de packages ...

En général, vous pouvez réduire le nombre de comparaisons pour un algorithme minmax en traitant deux éléments à la fois et en comparant uniquement le plus petit au minimum temporaire et le plus grand au maximum temporaire. En moyenne, il suffit de 3/4 des comparaisons qu'une approche naïve.

Cela pourrait être implémenté en c ou fortran (ou tout autre langage de bas niveau) et devrait être presque imbattable en termes de performances. j'utilisenumba pour illustrer le principe et obtenir une implémentation très rapide, indépendante de 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

C'est définitivement plus rapide que l'approche naïve présentée par 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

Comme prévu, la nouvelle implémentation minmax ne prend qu'environ 3/4 du temps de mise en œuvre naïve ( 2.1 / 2.75 = 0.7636363636363637)

Juste pour avoir quelques idées sur les chiffres auxquels on peut s'attendre, étant donné les approches suivantes:

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]

(les extrema_loop_*()approches sont similaires à ce qui est proposé ici , tandis que les extrema_while_*()approches sont basées sur le code d' ici )

Les horaires suivants:

indiquent que les extrema_while_*()sont les plus rapides, les plus extrema_while_nb()rapides Dans tous les cas, les solutions extrema_loop_nb()et extrema_loop_cy()surpassent également l'approche NumPy uniquement (en utilisant np.max()et np.min()séparément).

Enfin, notez qu'aucun de ceux-ci n'est aussi flexible que np.min()/ np.max()(en termes de support n-dim, axisparamètre, etc.).

(le code complet est disponible ici )

Personne n'a mentionné numpy.percentile , alors j'ai pensé que je le ferais. Si vous demandez des [0, 100]centiles, cela vous donnera un tableau de deux éléments, le min (0e centile) et le max (100e centile).

Cependant, cela ne répond pas à l'objectif de l'OP: ce n'est pas plus rapide que min et max séparément. Cela est probablement dû à certaines machines qui permettraient des percentiles non extrêmes (un problème plus difficile, qui devrait prendre plus de temps).

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'

Une future version de Numpy pourrait mettre dans un cas spécial pour ignorer le calcul normal du percentile si seulement [0, 100]sont demandés. Sans rien ajouter à l'interface, il existe un moyen de demander à Numpy le minimum et le maximum en un seul appel (contrairement à ce qui a été dit dans la réponse acceptée), mais l'implémentation standard de la bibliothèque ne profite pas de ce cas pour le faire digne d'intérêt.

C'est un vieux fil de discussion, mais de toute façon, si quelqu'un regarde à nouveau cela ...

Lors de la recherche simultanée du min et du max, il est possible de réduire le nombre de comparaisons. Si vous comparez des flotteurs (ce que je suppose), cela pourrait vous faire gagner du temps, mais pas en complexité de calcul.

Au lieu de (code 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]

vous pouvez d'abord comparer deux valeurs adjacentes dans le tableau, puis comparer uniquement la plus petite au minimum actuel et la plus grande au maximum actuel:

## 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

Le code ici est écrit en Python, clairement pour la vitesse, vous utiliseriez C ou Fortran ou Cython, mais de cette façon, vous faites 3 comparaisons par itération, avec len (ar) / 2 itérations, ce qui donne des comparaisons 3/2 * len (ar). Par opposition à cela, faire la comparaison "de manière évidente" vous faites deux comparaisons par itération, conduisant à des comparaisons 2 * len (ar). Économise 25% du temps de comparaison.

Peut-être que quelqu'un trouvera un jour cela utile.

À première vue, semble faire l'affaire:numpy.histogram

count, (amin, amax) = numpy.histogram(a, bins=1)

... mais si vous regardez la source de cette fonction, elle appelle simplement a.min()et a.max()indépendamment, et ne parvient donc pas à éviter les problèmes de performances abordés dans cette question. :-(

De même, cela scipy.ndimage.measurements.extremaressemble à une possibilité, mais il appelle aussi simplement a.min()et a.max()indépendamment.

Cela valait la peine pour moi de toute façon, donc je proposerai ici la solution la plus difficile et la moins élégante pour quiconque pourrait être intéressé. Ma solution consiste à implémenter un algorithme min-max multithread en un seul passage en C ++ et à l'utiliser pour créer un module d'extension Python. Cet effort nécessite un peu de temps système pour apprendre à utiliser les API Python et NumPy C / C ++, et ici je vais montrer le code et donner quelques petites explications et références pour quiconque souhaite emprunter cette voie.

Multi-threadé Min / Max

Il n'y a rien de trop intéressant ici. Le tableau est divisé en morceaux de taille length / workers. Le min / max est calculé pour chaque morceau dans a future, qui sont ensuite analysés pour le min / max global.

    // 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

Le module d'extension Python

C'est là que les choses commencent à devenir moche ... Une façon d'utiliser le code C ++ en Python est d'implémenter un module d'extension. Ce module peut être construit et installé à l'aide du distutils.coremodule standard. Une description complète de ce que cela implique est couverte dans la documentation Python: https://docs.python.org/3/extending/extending.html . REMARQUE: il existe certainement d'autres moyens d'obtenir des résultats similaires, pour citer https://docs.python.org/3/extending/index.html#extending-index :

Ce guide ne couvre que les outils de base pour la création d'extensions fournis dans le cadre de cette version de CPython. Des outils tiers tels que Cython, cffi, SWIG et Numba offrent des approches à la fois plus simples et plus sophistiquées pour créer des extensions C et C ++ pour Python.

Essentiellement, cette voie est probablement plus académique que pratique. Cela étant dit, ce que j'ai fait ensuite, c'était, en m'en tenant assez près du didacticiel, de créer un fichier de module. C'est essentiellement un passe-partout pour que les distutils sachent quoi faire avec votre code et en créer un module Python. Avant de faire quoi que ce soit, il est probablement sage de créer un environnement virtuel Python afin de ne pas polluer vos packages système (voir https://docs.python.org/3/library/venv.html#module-venv ).

Voici le fichier du module:

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

Dans ce fichier, il y a une utilisation significative de Python ainsi que de l'API NumPy, pour plus d'informations, consultez: https://docs.python.org/3/c-api/arg.html#c.PyArg_ParseTuple , et pour NumPy : https://docs.scipy.org/doc/numpy/reference/c-api.array.html .

Installation du module

La prochaine chose à faire est d'utiliser distutils pour installer le module. Cela nécessite un fichier d'installation:

# 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])

Pour enfin installer le module, exécutez à python3 setup.py installpartir de votre environnement virtuel.

Test du module

Enfin, nous pouvons tester pour voir si l'implémentation C ++ surpasse réellement l'utilisation naïve de NumPy. Pour ce faire, voici un script de test simple:

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

Voici les résultats que j'ai obtenus en faisant tout cela:

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  

Celles-ci sont beaucoup moins encourageantes que les résultats indiquent plus tôt dans le fil, qui indiquaient une accélération d'environ 3,5x, et n'intégraient pas le multi-threading. Les résultats que j'ai obtenus sont quelque peu raisonnables, je m'attendrais à ce que la surcharge de threading et domine le temps jusqu'à ce que les tableaux deviennent très grands, auquel point l'augmentation des performances commencerait à se rapprocher de std::thread::hardware_concurrencyx augmenter.

Conclusion

Il y a certainement de la place pour des optimisations spécifiques aux applications pour certains codes NumPy, semble-t-il, en particulier en ce qui concerne le multi-threading. Je ne sais pas si cela en vaut la peine ou non, mais cela semble certainement être un bon exercice (ou quelque chose du genre). Je pense que peut-être apprendre certains de ces "outils tiers" comme Cython peut être une meilleure utilisation du temps, mais qui sait.

Le moyen le plus court que j'ai trouvé est le suivant:

mn, mx = np.sort(ar)[[0, -1]]

Mais comme il trie le tableau, ce n'est pas le plus efficace.

Un autre chemin court serait:

mn, mx = np.percentile(ar, [0, 100])

Cela devrait être plus efficace, mais le résultat est calculé et un flottant est renvoyé.