numpy.amax () encontrará el valor máximo en una matriz, y numpy.amin () hace lo mismo con el valor mínimo. Si quiero encontrar tanto el máximo como el mínimo, tengo que llamar a ambas funciones, lo que requiere pasar dos veces la matriz (muy grande), lo que parece lento.

¿Existe una función en la API numpy que encuentre tanto el máximo como el mínimo con solo un paso a través de los datos?

¿Existe una función en la API numpy que encuentre tanto el máximo como el mínimo con solo un paso a través de los datos?

No. En el momento de escribir este artículo, no existe tal función. (Y sí, si no eran tal función, su rendimiento sería significativamente mejor que llamar numpy.amin()y numpy.amax()sucesivamente sobre una gran variedad.)

No creo que pasar por encima de la matriz dos veces sea un problema. Considere el siguiente pseudocódigo:

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

Si bien solo hay 1 bucle aquí, todavía hay 2 controles. (En lugar de tener 2 bucles con 1 cheque cada uno). Realmente, lo único que guarda es la sobrecarga de 1 bucle. Si las matrices son realmente grandes como dice, esa sobrecarga es pequeña en comparación con la carga de trabajo real del ciclo. (Tenga en cuenta que todo esto está implementado en C, por lo que los bucles son más o menos libres de todos modos).

EDITAR Lo siento por los 4 que votaron a favor y tuvieron fe en mí. Definitivamente puedes optimizar esto.

Aquí hay un código de fortran que se puede compilar en un módulo de Python a través de f2py(tal vez un Cythongurú pueda venir y comparar esto con una versión C optimizada ...):

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

Compílelo a través de:

f2py -m untitled -c fortran_code.f90

Y ahora estamos en un lugar donde podemos probarlo:

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'

Los resultados son un poco asombrosos para mí:

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

Tengo que decir que no lo entiendo del todo. Comparar solo np.minversus minmax1y minmax2sigue siendo una batalla perdida, por lo que no es solo un problema de memoria ...

notas : aumentar el tamaño en un factor de 10**ay disminuir la repetición en un factor de 10**a(mantener constante el tamaño del problema) cambia el rendimiento, pero no de una manera aparentemente consistente, lo que muestra que hay alguna interacción entre el rendimiento de la memoria y la sobrecarga de llamadas de función en pitón. Incluso comparar una minimplementación simple en fortran supera a la de numpy por un factor de aproximadamente 2 ...

Hay una función para encontrar (max-min) llamada numpy.ptp si le resulta útil:

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

pero no creo que haya una manera de encontrar tanto el mínimo como el máximo con un recorrido.

EDITAR: ptp solo llama mínimo y máximo bajo el capó

Puede usar Numba , que es un compilador dinámico de Python compatible con NumPy que usa LLVM. La implementación resultante es bastante simple y clara:

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

También debería ser más rápido que la min() & max()implementación de Numpy . Y todo ello sin tener que escribir una sola línea de código C / Fortran.

Haz tus propias pruebas de rendimiento, ya que siempre depende de tu arquitectura, tus datos, las versiones de tus paquetes ...

En general, puede reducir la cantidad de comparaciones para un algoritmo minmax procesando dos elementos a la vez y solo comparando el más pequeño con el mínimo temporal y el más grande con el máximo temporal. En promedio, solo se necesitan 3/4 de las comparaciones que un enfoque ingenuo.

Esto podría implementarse en co fortran (o cualquier otro lenguaje de bajo nivel) y debería ser casi imbatible en términos de rendimiento. Estoy usandonumba para ilustrar el principio y obtener una implementación muy rápida, independiente 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

Definitivamente es más rápido que el enfoque ingenuo que presentó 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

Como se esperaba, la nueva implementación de minmax solo toma aproximadamente 3/4 del tiempo que tomó la implementación ingenua ( 2.1 / 2.75 = 0.7636363636363637)

Solo para obtener algunas ideas sobre los números que uno podría esperar, dados los siguientes enfoques:

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]

(los extrema_loop_*()enfoques son similares a los que se proponen aquí , mientras que los extrema_while_*()enfoques se basan en el código de aquí )

Los siguientes horarios:

indican que extrema_while_*()son los más rápidos, con los extrema_while_nb()más rápidos. En cualquier caso, también las soluciones extrema_loop_nb()y extrema_loop_cy()superan el enfoque de solo NumPy (usando np.max()y por np.min()separado).

Finalmente, tenga en cuenta que ninguno de estos es tan flexible como np.min()/ np.max()(en términos de soporte n-dim, axisparámetro, etc.).

(el código completo está disponible aquí )

Nadie mencionó numpy.percentile , así que pensé que lo haría. Si pregunta por [0, 100]percentiles, le dará una matriz de dos elementos, el mínimo (percentil 0) y el máximo (percentil 100).

Sin embargo, no satisface el propósito del OP: no es más rápido que el mínimo y el máximo por separado. Eso probablemente se deba a alguna maquinaria que permitiría percentiles no extremos (un problema más difícil, que debería llevar más tiempo).

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'

Una versión futura de Numpy podría incluir un caso especial para omitir el cálculo del percentil normal si solo [0, 100]se solicita. Sin agregar nada a la interfaz, hay una manera de pedirle a Numpy el mínimo y el máximo en una llamada (al contrario de lo que se dijo en la respuesta aceptada), pero la implementación estándar de la biblioteca no aprovecha este caso para hacerlo vale la pena.

Este es un hilo antiguo, pero de todos modos, si alguien vuelve a mirar esto ...

Al buscar el mínimo y el máximo simultáneamente, es posible reducir el número de comparaciones. Si está comparando flotadores (lo que supongo), esto podría ahorrarle algo de tiempo, aunque no complejidad computacional.

En lugar de (código 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]

primero puede comparar dos valores adyacentes en la matriz, y luego solo comparar el más pequeño con el mínimo actual y el más grande con el máximo actual:

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

El código aquí está escrito en Python, claramente para la velocidad usarías C o Fortran o Cython, pero de esta manera haces 3 comparaciones por iteración, con len (ar) / 2 iteraciones, dando 3/2 * len (ar) comparaciones. A diferencia de eso, al hacer la comparación "de la manera obvia" se hacen dos comparaciones por iteración, lo que lleva a comparaciones 2 * len (ar). Le ahorra un 25% de tiempo de comparación.

Tal vez alguien algún día lo encuentre útil.

A primera vista, parece hacer el truco:numpy.histogram

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

... pero si nos fijamos en la fuente de esa función, simplemente se llama a.min()y a.max()de forma independiente, y por lo tanto no puede evitar los problemas de rendimiento abordan en esta pregunta. :-(

Del mismo modo, scipy.ndimage.measurements.extremaparece una posibilidad, pero también, simplemente llama a.min()y de forma a.max()independiente.

De todos modos, valió la pena el esfuerzo para mí, así que propondré aquí la solución más difícil y menos elegante para quien pueda estar interesado. Mi solución es implementar un mínimo-máximo de subprocesos múltiples en un algoritmo de un paso en C ++, y usarlo para crear un módulo de extensión de Python. Este esfuerzo requiere un poco de sobrecarga para aprender a usar las API de Python y NumPy C / C ++, y aquí mostraré el código y daré algunas pequeñas explicaciones y referencias para quien desee seguir este camino.

Min / Max multiproceso

No hay nada demasiado interesante aquí. La matriz se divide en trozos de tamaño length / workers. El mínimo / máximo se calcula para cada fragmento en a future, que luego se escanea para el mínimo / máximo 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

El módulo de extensión de Python

Aquí es donde las cosas empiezan a ponerse feas ... Una forma de usar el código C ++ en Python es implementar un módulo de extensión. Este módulo se puede construir e instalar utilizando el distutils.coremódulo estándar. Una descripción completa de lo que esto implica se cubre en la documentación de Python: https://docs.python.org/3/extending/extending.html . NOTA: ciertamente hay otras formas de obtener resultados similares, por citar https://docs.python.org/3/extending/index.html#extending-index :

Esta guía solo cubre las herramientas básicas para crear extensiones proporcionadas como parte de esta versión de CPython. Las herramientas de terceros como Cython, cffi, SWIG y Numba ofrecen enfoques más simples y sofisticados para crear extensiones C y C ++ para Python.

Esencialmente, esta ruta es probablemente más académica que práctica. Habiendo dicho eso, lo que hice a continuación fue, manteniéndome bastante cerca del tutorial, crear un archivo de módulo. Esto es esencialmente un texto estándar para que los distutils sepan qué hacer con su código y creen un módulo de Python a partir de él. Antes de hacer algo de esto, probablemente sea aconsejable crear un entorno virtual de Python para no contaminar los paquetes de su sistema (consulte https://docs.python.org/3/library/venv.html#module-venv ).

Aquí está el archivo del módulo:

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

En este archivo hay un uso significativo tanto de Python como de la API de NumPy, para más información consultar: https://docs.python.org/3/c-api/arg.html#c.PyArg_ParseTuple , y para NumPy : https://docs.scipy.org/doc/numpy/reference/c-api.array.html .

Instalación del módulo

Lo siguiente que debe hacer es utilizar distutils para instalar el módulo. Esto requiere un archivo de instalación:

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

Para finalmente instalar el módulo, ejecute python3 setup.py install desde su entorno virtual.

Prueba del módulo

Finalmente, podemos probar para ver si la implementación de C ++ realmente supera al uso ingenuo de NumPy. Para hacerlo, aquí hay un script de prueba 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))

Estos son los resultados que obtuve al hacer todo esto:

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  

Estos son mucho menos alentadores de lo que los resultados indican anteriormente en el hilo, que indicaron en algún lugar una aceleración de alrededor de 3.5x, y no incorporaron subprocesos múltiples. Los resultados que obtuve son algo razonables, esperaría que la sobrecarga de subprocesos y dominaría el tiempo hasta que las matrices se volvieran muy grandes, momento en el que el aumento de rendimiento comenzaría a acercarse a std::thread::hardware_concurrencyx aumento.

Conclusión

Ciertamente, parece que hay espacio para optimizaciones específicas de aplicaciones para algunos códigos de NumPy, en particular con respecto al subproceso múltiple. No tengo claro si vale la pena el esfuerzo o no, pero ciertamente parece un buen ejercicio (o algo así). Creo que quizás aprender algunas de esas "herramientas de terceros" como Cython puede ser un mejor uso del tiempo, pero quién sabe.

La forma más corta que he encontrado es esta:

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

Pero dado que ordena la matriz, no es la más eficiente.

Otra forma corta sería:

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

Esto debería ser más eficiente, pero el resultado se calcula y se devuelve un flotante.