numpy.amax () encontrará o valor máximo em uma matriz, e numpy.amin () faz o mesmo para o valor mínimo. Se eu quiser encontrar max e min, tenho que chamar as duas funções, o que requer a passagem da matriz (muito grande) duas vezes, o que parece lento.

Existe uma função na API numpy que encontra o máximo e o mínimo com apenas uma única passagem pelos dados?

Existe uma função na API numpy que encontra o máximo e o mínimo com apenas uma única passagem pelos dados?

Não. No momento da redação deste artigo, essa função não existia. (E sim, se não fosse essa função, seu desempenho seria significativamente melhor do que chamar numpy.amin()e numpy.amax()sucessivamente em uma grande variedade.)

Não acho que passar sobre o array duas vezes seja um problema. Considere o seguinte pseudocódigo:

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

Embora haja apenas 1 loop aqui, ainda há 2 verificações. (Em vez de ter 2 loops com 1 cheque cada). Na verdade, a única coisa que você salva é a sobrecarga de 1 loop. Se os arrays realmente forem grandes, como você diz, a sobrecarga é pequena em comparação com a carga de trabalho do loop real. (Observe que tudo isso é implementado em C, então os loops são mais ou menos livres de qualquer maneira).

EDITAR Sinto muito para os 4 de vocês que votaram positivamente e tiveram fé em mim. Você definitivamente pode otimizar isso.

Aqui está um código fortran que pode ser compilado em um módulo Python via f2py(talvez um Cythonguru possa vir e comparar isso com uma versão C otimizada ...):

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

Compile-o via:

f2py -m untitled -c fortran_code.f90

E agora estamos em um lugar onde podemos testá-lo:

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'

Os resultados são um pouco surpreendentes para mim:

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

Devo dizer que não entendo completamente. Comparar apenas np.minversus minmax1e minmax2ainda é uma batalha perdida, então não é apenas um problema de memória ...

notas - Aumentar o tamanho por um fator de 10**ae diminuir a repetição por um fator de 10**a(mantendo o tamanho do problema constante) muda o desempenho, mas não de uma forma aparentemente consistente, o que mostra que há alguma interação entre o desempenho da memória e a sobrecarga da chamada de função em Pitão. Mesmo comparando uma minimplementação simples em fortran bate numpy por um fator de aproximadamente 2 ...

Existe uma função para encontrar (max-min) chamada numpy.ptp se for útil para você:

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

mas não acho que haja uma maneira de encontrar o mínimo e o máximo com uma travessia.

EDIT: ptp apenas chama min e max sob o capô

Você poderia usar o Numba , que é um compilador Python dinâmico compatível com o NumPy usando LLVM. A implementação resultante é muito simples e 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()))

Também deve ser mais rápido do que a min() & max()implementação de um Numpy . E tudo isso sem precisar escrever uma única linha de código C / Fortran.

Faça seus próprios testes de desempenho, pois é sempre dependente de sua arquitetura, seus dados, suas versões de pacote ...

Em geral, você pode reduzir a quantidade de comparações para um algoritmo minmax processando dois elementos por vez e comparando apenas o menor com o mínimo temporário e o maior com o máximo temporário. Em média, são necessários apenas 3/4 das comparações do que uma abordagem ingênua.

Isso pode ser implementado em c ou fortran (ou qualquer outra linguagem de baixo nível) e deve ser quase imbatível em termos de desempenho. estou a usarnumba para ilustrar o princípio e obter uma implementação muito rápida e independente do tipo:

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 mais rápido do que a abordagem ingênua que Peque apresentou:

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 esperado, a nova implementação de minmax leva apenas cerca de 3/4 do tempo que a implementação ingênua levou ( 2.1 / 2.75 = 0.7636363636363637)

Apenas para ter algumas idéias sobre os números que você pode esperar, dadas as seguintes abordagens:

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]

(as extrema_loop_*()abordagens são semelhantes ao que é proposto aqui , enquanto as extrema_while_*()abordagens são baseadas no código a partir daqui )

Os seguintes horários:

indicam que extrema_while_*()são os mais rápidos, extrema_while_nb()sendo os mais rápidos. Em qualquer caso, também as soluções extrema_loop_nb()e extrema_loop_cy()superam a abordagem apenas NumPy (usando np.max()e np.min()separadamente).

Finalmente, observe que nenhum desses é tão flexível quanto np.min()/ np.max()(em termos de suporte n-dim, axisparâmetro, etc.).

(o código completo está disponível aqui )

Ninguém mencionou numpy.percentile , então pensei que faria. Se você solicitar os [0, 100]percentis, receberá uma matriz de dois elementos, o mínimo (0º percentil) e o máximo (100º percentil).

No entanto, não satisfaz o propósito do OP: não é mais rápido do que mín e máx separadamente. Isso provavelmente se deve a algum mecanismo que permitiria percentis não extremos (um problema mais difícil, que deve demorar mais).

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'

Uma versão futura do Numpy poderia colocar em um caso especial para pular o cálculo de percentil normal se apenas [0, 100]for solicitado. Sem adicionar nada à interface, há uma maneira de pedir a Numpy o mínimo e o máximo em uma chamada (ao contrário do que foi dito na resposta aceita), mas a implementação padrão da biblioteca não aproveita este caso para fazê-lo que vale a pena.

Este é um tópico antigo, mas de qualquer forma, se alguém olhar para isso novamente ...

Ao procurar o mínimo e o máximo simultaneamente, é possível reduzir o número de comparações. Se for floats que você está comparando (o que eu acho que é), isso pode economizar algum tempo, embora não seja uma complexidade computacional.

Em vez 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]

você pode primeiro comparar dois valores adjacentes na matriz e, em seguida, comparar apenas o menor com o mínimo atual e o maior com o máximo atual:

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

O código aqui é escrito em Python, claramente para velocidade você usaria C ou Fortran ou Cython, mas desta forma você faz 3 comparações por iteração, com len (ar) / 2 iterações, dando 3/2 * len (ar) comparações. Em oposição a isso, fazendo a comparação "da maneira óbvia", você faz duas comparações por iteração, resultando em comparações 2 * len (ar). Economiza 25% do tempo de comparação.

Talvez alguém um dia ache isso útil.

À primeira vista, parece funcionar:numpy.histogram

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

... mas se você olhar a origem dessa função, ela simplesmente chama a.min()e de forma a.max()independente e, portanto, falha em evitar as preocupações de desempenho abordadas nesta questão. :-(

Da mesma forma, scipy.ndimage.measurements.extremaparece uma possibilidade, mas também simplesmente chama a.min()e de forma a.max()independente.

Valeu a pena o esforço para mim de qualquer forma, então vou propor a solução mais difícil e menos elegante aqui para quem estiver interessado. Minha solução é implementar um algoritmo mín-máx multissegmentado em uma passagem em C ++ e usar isso para criar um módulo de extensão Python. Esse esforço requer um pouco de sobrecarga para aprender a usar as APIs Python e NumPy C / C ++, e aqui vou mostrar o código e dar algumas pequenas explicações e referências para quem deseja seguir esse caminho.

Multi-threaded Min / Max

Não há nada muito interessante aqui. A matriz é dividida em pedaços de tamanho length / workers. O mín. / Máx. É calculado para cada pedaço em a future, que é então verificado para o mín. / Máx. 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

O Módulo de Extensão Python

É aqui que as coisas começam a ficar feias ... Uma maneira de usar o código C ++ em Python é implementar um módulo de extensão. Este módulo pode ser construído e instalado usando o distutils.coremódulo padrão. Uma descrição completa do que isso implica é coberta na documentação do Python: https://docs.python.org/3/extending/extending.html . NOTA: certamente existem outras maneiras de obter resultados semelhantes, para citar https://docs.python.org/3/extending/index.html#extending-index :

Este guia cobre apenas as ferramentas básicas para a criação de extensões fornecidas como parte desta versão do CPython. Ferramentas de terceiros como Cython, cffi, SWIG e Numba oferecem abordagens mais simples e sofisticadas para a criação de extensões C e C ++ para Python.

Essencialmente, esse caminho é provavelmente mais acadêmico do que prático. Com isso dito, o que fiz a seguir foi, mantendo-me bem próximo ao tutorial, criar um arquivo de módulo. Este é essencialmente um padrão para distutils saber o que fazer com seu código e criar um módulo Python a partir dele. Antes de fazer qualquer coisa, provavelmente é aconselhável criar um ambiente virtual Python para não poluir os pacotes do sistema (consulte https://docs.python.org/3/library/venv.html#module-venv ).

Aqui está o arquivo do 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); }

Neste arquivo há um uso significativo do Python, bem como da API NumPy, para mais informações consulte: https://docs.python.org/3/c-api/arg.html#c.PyArg_ParseTuple e para NumPy : https://docs.scipy.org/doc/numpy/reference/c-api.array.html .

Instalando o Módulo

A próxima coisa a fazer é utilizar distutils para instalar o módulo. Isso requer um arquivo de configuração:

# 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 o módulo, execute python3 setup.py install partir de seu ambiente virtual.

Testando o Módulo

Por fim, podemos testar para ver se a implementação C ++ realmente supera o uso ingênuo de NumPy. Para fazer isso, aqui está um script de teste simples:

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

Aqui estão os resultados que obtive com tudo isso:

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  

Estes são muito menos encorajadores do que os resultados indicam anteriormente no thread, que indicou algo em torno de 3,5x a aceleração e não incorporou multi-threading. Os resultados que obtive são um tanto razoáveis, eu esperaria que a sobrecarga de threading e dominasse o tempo até que os arrays ficassem muito grandes, ponto em que o aumento de desempenho começaria a se aproximar de std::thread::hardware_concurrencyaumento x.

Conclusão

Certamente há espaço para otimizações específicas de aplicativos para algum código NumPy, ao que parece, em particular com relação a multi-threading. Se vale a pena ou não o esforço não está claro para mim, mas certamente parece um bom exercício (ou algo assim). Acho que talvez aprender algumas dessas "ferramentas de terceiros" como o Cython possa ser um uso melhor do tempo, mas quem sabe.

O caminho mais curto que descobri é este:

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

Mas, uma vez que classifica o array, não é o mais eficiente.

Outro caminho curto seria:

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

Isso deve ser mais eficiente, mas o resultado é calculado e um float é retornado.