Produto cartesiano de pontos de matriz xey em matriz única de pontos 2D

147

Eu tenho duas matrizes numpy que definem os eixos xey de uma grade. Por exemplo:

x = numpy.array([1,2,3])
y = numpy.array([4,5])

Eu gostaria de gerar o produto cartesiano dessas matrizes para gerar:

array([[1,4],[2,4],[3,4],[1,5],[2,5],[3,5]])

De uma maneira que não é terrivelmente ineficiente, pois preciso fazer isso várias vezes seguidas. Estou assumindo que convertê-los para uma lista Python e usar itertools.producte voltar para uma matriz numpy não é a forma mais eficiente.

Rico
fonte
Percebi que o passo mais caro na abordagem de ferramentas é a conversão final de lista em matriz. Sem esse último passo, é duas vezes mais rápido que o exemplo de Ken.
Alexey Lebedev

Respostas:

88
>>> numpy.transpose([numpy.tile(x, len(y)), numpy.repeat(y, len(x))])
array([[1, 4],
       [2, 4],
       [3, 4],
       [1, 5],
       [2, 5],
       [3, 5]])

Consulte Usando numpy para criar uma matriz de todas as combinações de duas matrizes para obter uma solução geral para calcular o produto cartesiano de N matrizes.

kennytm
fonte
1
Uma vantagem dessa abordagem é que ela produz resultados consistentes para matrizes do mesmo tamanho. A abordagem meshgrid+ dstack, embora mais rápida em alguns casos, pode levar a erros se você espera que o produto cartesiano seja construído na mesma ordem para matrizes do mesmo tamanho.
tlnagy
3
@tlnagy, não notei nenhum caso em que essa abordagem produz resultados diferentes daqueles produzidos por meshgrid+ dstack. Você poderia postar um exemplo?
Senderle
148

Um canônico cartesian_product(quase)

Existem muitas abordagens para esse problema com propriedades diferentes. Alguns são mais rápidos que outros, e outros são de uso geral. Após muitos testes e ajustes, descobri que a função a seguir, que calcula uma dimensão n cartesian_product, é mais rápida que a maioria das outras para muitas entradas. Para um par de abordagens um pouco mais complexas, mas ainda mais rápidas em muitos casos, veja a resposta de Paul Panzer .

Dada essa resposta, essa não é mais a implementação mais rápida do produto cartesiano em numpyque estou ciente. No entanto, acho que sua simplicidade continuará sendo uma referência útil para melhorias futuras:

def cartesian_product(*arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[...,i] = a
    return arr.reshape(-1, la)

Vale ressaltar que essa função usa de ix_maneira incomum; enquanto o uso documentado de ix_é gerar índices em uma matriz, acontece que matrizes com a mesma forma podem ser usadas para atribuição transmitida. Muito obrigado a mgilson , que me inspirou a tentar usar ix_esse caminho, e a unutbu , que forneceu alguns comentários extremamente úteis sobre essa resposta, incluindo a sugestão de uso numpy.result_type.

Alternativas notáveis

Às vezes, é mais rápido gravar blocos de memória contíguos na ordem do Fortran. Essa é a base dessa alternativa, cartesian_product_transposeque se mostrou mais rápida em alguns hardwares do que cartesian_product(veja abaixo). No entanto, a resposta de Paul Panzer, que usa o mesmo princípio, é ainda mais rápida. Ainda assim, incluo isso aqui para leitores interessados:

def cartesian_product_transpose(*arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
    dtype = numpy.result_type(*arrays)

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

Depois de entender a abordagem de Panzer, escrevi uma nova versão quase tão rápida quanto a dele e quase tão simples quanto cartesian_product:

def cartesian_product_simple_transpose(arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([la] + [len(a) for a in arrays], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[i, ...] = a
    return arr.reshape(la, -1).T

Parece haver alguma sobrecarga de tempo constante que a torna mais lenta que a do Panzer para pequenas entradas. Mas para entradas maiores, em todos os testes que eu executei, ele executa tão bem quanto sua implementação mais rápida ( cartesian_product_transpose_pp).

Nas seções a seguir, incluo alguns testes de outras alternativas. Agora eles estão um pouco desatualizados, mas, em vez de um esforço duplicado, decidi deixá-los aqui fora de interesse histórico. Para testes atualizados, consulte a resposta de Panzer, bem como a de Nico Schlömer .

Testes contra alternativas

Aqui está uma bateria de testes que mostram o aumento de desempenho que algumas dessas funções fornecem em relação a várias alternativas. Todos os testes mostrados aqui foram realizados em uma máquina quad-core, executando o Mac OS 10.12.5, Python 3.6.1 e numpy1.12.1. Sabe-se que variações no hardware e software produzem resultados diferentes, portanto, YMMV. Execute esses testes para ter certeza!

Definições:

import numpy
import itertools
from functools import reduce

### Two-dimensional products ###

def repeat_product(x, y):
    return numpy.transpose([numpy.tile(x, len(y)), 
                            numpy.repeat(y, len(x))])

def dstack_product(x, y):
    return numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)

### Generalized N-dimensional products ###

def cartesian_product(*arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[...,i] = a
    return arr.reshape(-1, la)

def cartesian_product_transpose(*arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
    dtype = numpy.result_type(*arrays)

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

# from https://stackoverflow.com/a/1235363/577088

def cartesian_product_recursive(*arrays, out=None):
    arrays = [numpy.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = numpy.prod([x.size for x in arrays])
    if out is None:
        out = numpy.zeros([n, len(arrays)], dtype=dtype)

    m = n // arrays[0].size
    out[:,0] = numpy.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian_product_recursive(arrays[1:], out=out[0:m,1:])
        for j in range(1, arrays[0].size):
            out[j*m:(j+1)*m,1:] = out[0:m,1:]
    return out

def cartesian_product_itertools(*arrays):
    return numpy.array(list(itertools.product(*arrays)))

### Test code ###

name_func = [('repeat_product',                                                 
              repeat_product),                                                  
             ('dstack_product',                                                 
              dstack_product),                                                  
             ('cartesian_product',                                              
              cartesian_product),                                               
             ('cartesian_product_transpose',                                    
              cartesian_product_transpose),                                     
             ('cartesian_product_recursive',                           
              cartesian_product_recursive),                            
             ('cartesian_product_itertools',                                    
              cartesian_product_itertools)]

def test(in_arrays, test_funcs):
    global func
    global arrays
    arrays = in_arrays
    for name, func in test_funcs:
        print('{}:'.format(name))
        %timeit func(*arrays)

def test_all(*in_arrays):
    test(in_arrays, name_func)

# `cartesian_product_recursive` throws an 
# unexpected error when used on more than
# two input arrays, so for now I've removed
# it from these tests.

def test_cartesian(*in_arrays):
    test(in_arrays, name_func[2:4] + name_func[-1:])

x10 = [numpy.arange(10)]
x50 = [numpy.arange(50)]
x100 = [numpy.arange(100)]
x500 = [numpy.arange(500)]
x1000 = [numpy.arange(1000)]

Resultado dos testes:

In [2]: test_all(*(x100 * 2))
repeat_product:
67.5 µs ± 633 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
dstack_product:
67.7 µs ± 1.09 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product:
33.4 µs ± 558 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product_transpose:
67.7 µs ± 932 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product_recursive:
215 µs ± 6.01 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_itertools:
3.65 ms ± 38.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [3]: test_all(*(x500 * 2))
repeat_product:
1.31 ms ± 9.28 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
dstack_product:
1.27 ms ± 7.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product:
375 µs ± 4.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_transpose:
488 µs ± 8.88 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_recursive:
2.21 ms ± 38.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
105 ms ± 1.17 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [4]: test_all(*(x1000 * 2))
repeat_product:
10.2 ms ± 132 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
dstack_product:
12 ms ± 120 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product:
4.75 ms ± 57.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_transpose:
7.76 ms ± 52.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_recursive:
13 ms ± 209 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
422 ms ± 7.77 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Em todos os casos, cartesian_productconforme definido no início desta resposta, é mais rápido.

Para aquelas funções que aceitam um número arbitrário de matrizes de entrada, vale a pena verificar o desempenho len(arrays) > 2também. (Até que eu possa determinar por que causa cartesian_product_recursiveum erro nesse caso, eu o removi desses testes.)

In [5]: test_cartesian(*(x100 * 3))
cartesian_product:
8.8 ms ± 138 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_transpose:
7.87 ms ± 91.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
518 ms ± 5.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [6]: test_cartesian(*(x50 * 4))
cartesian_product:
169 ms ± 5.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_transpose:
184 ms ± 4.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_itertools:
3.69 s ± 73.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [7]: test_cartesian(*(x10 * 6))
cartesian_product:
26.5 ms ± 449 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_transpose:
16 ms ± 133 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
728 ms ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [8]: test_cartesian(*(x10 * 7))
cartesian_product:
650 ms ± 8.14 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
cartesian_product_transpose:
518 ms ± 7.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
cartesian_product_itertools:
8.13 s ± 122 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Como esses testes mostram, cartesian_productpermanece competitivo até que o número de matrizes de entrada ultrapasse (aproximadamente) quatro. Depois disso, cartesian_product_transposetem uma ligeira vantagem.

Vale a pena reiterar que usuários com outros hardwares e sistemas operacionais podem ter resultados diferentes. Por exemplo, o unutbu reporta os seguintes resultados para esses testes usando o Ubuntu 14.04, Python 3.4.3 e numpy1.14.0.dev0 + b7050a9:

>>> %timeit cartesian_product_transpose(x500, y500) 
1000 loops, best of 3: 682 µs per loop
>>> %timeit cartesian_product(x500, y500)
1000 loops, best of 3: 1.55 ms per loop

Abaixo, mostro alguns detalhes sobre os testes anteriores que realizei nesse sentido. O desempenho relativo dessas abordagens mudou ao longo do tempo, para diferentes hardwares e diferentes versões do Python e numpy. Embora não seja imediatamente útil para pessoas que usam versões atualizadas numpy, ele ilustra como as coisas mudaram desde a primeira versão desta resposta.

Uma alternativa simples: meshgrid+dstack

A resposta atualmente aceita usa tilee repeatpara transmitir duas matrizes juntas. Mas a meshgridfunção faz praticamente a mesma coisa. Aqui está a saída tilee repeatantes de ser passada para transposição:

In [1]: import numpy
In [2]: x = numpy.array([1,2,3])
   ...: y = numpy.array([4,5])
   ...: 

In [3]: [numpy.tile(x, len(y)), numpy.repeat(y, len(x))]
Out[3]: [array([1, 2, 3, 1, 2, 3]), array([4, 4, 4, 5, 5, 5])]

E aqui está a saída de meshgrid:

In [4]: numpy.meshgrid(x, y)
Out[4]: 
[array([[1, 2, 3],
        [1, 2, 3]]), array([[4, 4, 4],
        [5, 5, 5]])]

Como você pode ver, é quase idêntico. Precisamos apenas remodelar o resultado para obter exatamente o mesmo resultado.

In [5]: xt, xr = numpy.meshgrid(x, y)
   ...: [xt.ravel(), xr.ravel()]
Out[5]: [array([1, 2, 3, 1, 2, 3]), array([4, 4, 4, 5, 5, 5])]

Em vez de remodelar neste ponto, porém, poderíamos passar a saída meshgridpara dstacke remodelar depois, o que economiza algum trabalho:

In [6]: numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)
Out[6]: 
array([[1, 4],
       [2, 4],
       [3, 4],
       [1, 5],
       [2, 5],
       [3, 5]])

Ao contrário da afirmação deste comentário , não vi nenhuma evidência de que entradas diferentes produzirão saídas de formas diferentes e, como demonstrado acima, elas fazem coisas muito semelhantes, portanto seria muito estranho se o fizessem. Entre em contato se você encontrar um contra-exemplo.

Testando meshgrid+ dstackvs. repeat+transpose

O desempenho relativo dessas duas abordagens mudou ao longo do tempo. Em uma versão anterior do Python (2.7), o resultado usando meshgrid+ dstackfoi notavelmente mais rápido para pequenas entradas. (Observe que esses testes são de uma versão antiga desta resposta.) Definições:

>>> def repeat_product(x, y):
...     return numpy.transpose([numpy.tile(x, len(y)), 
                                numpy.repeat(y, len(x))])
...
>>> def dstack_product(x, y):
...     return numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)
...     

Para entrada de tamanho médio, vi uma aceleração significativa. Mas tentei novamente esses testes com versões mais recentes do Python (3.6.1) e numpy(1.12.1), em uma máquina mais nova. As duas abordagens são quase idênticas agora.

Teste antigo

>>> x, y = numpy.arange(500), numpy.arange(500)
>>> %timeit repeat_product(x, y)
10 loops, best of 3: 62 ms per loop
>>> %timeit dstack_product(x, y)
100 loops, best of 3: 12.2 ms per loop

Novo teste

In [7]: x, y = numpy.arange(500), numpy.arange(500)
In [8]: %timeit repeat_product(x, y)
1.32 ms ± 24.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [9]: %timeit dstack_product(x, y)
1.26 ms ± 8.47 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Como sempre, YMMV, mas isso sugere que nas versões recentes do Python e numpy, elas são intercambiáveis.

Funções generalizadas do produto

Em geral, podemos esperar que o uso de funções internas seja mais rápido para entradas pequenas, enquanto que para entradas grandes, uma função criada para fins específicos pode ser mais rápida. Além disso, para um produto n-dimensional generalizado, tilee repeatnão ajudará, porque eles não possuem análogos claros de alta dimensão. Portanto, vale a pena investigar o comportamento das funções criadas especificamente para esse fim.

A maioria dos testes relevantes aparece no início desta resposta, mas aqui estão alguns dos testes realizados em versões anteriores do Python e numpypara comparação.

A cartesianfunção definida em outra resposta costumava ter um desempenho muito bom para entradas maiores. (É o mesmo que a função chamada cartesian_product_recursiveacima.) A fim de comparar cartesiana dstack_prodct, usamos apenas duas dimensões.

Aqui, novamente, o teste antigo mostrou uma diferença significativa, enquanto o novo teste mostra quase nenhum.

Teste antigo

>>> x, y = numpy.arange(1000), numpy.arange(1000)
>>> %timeit cartesian([x, y])
10 loops, best of 3: 25.4 ms per loop
>>> %timeit dstack_product(x, y)
10 loops, best of 3: 66.6 ms per loop

Novo teste

In [10]: x, y = numpy.arange(1000), numpy.arange(1000)
In [11]: %timeit cartesian([x, y])
12.1 ms ± 199 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [12]: %timeit dstack_product(x, y)
12.7 ms ± 334 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Como antes, dstack_productainda bate cartesianem escalas menores.

Novo teste ( teste antigo redundante não mostrado )

In [13]: x, y = numpy.arange(100), numpy.arange(100)
In [14]: %timeit cartesian([x, y])
215 µs ± 4.75 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [15]: %timeit dstack_product(x, y)
65.7 µs ± 1.15 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Acho que essas distinções são interessantes e merecem ser registradas; mas eles são acadêmicos no final. Como mostraram os testes no início desta resposta, todas essas versões são quase sempre mais lentas que cartesian_product, definidas no início desta resposta - o que é um pouco mais lento que as implementações mais rápidas entre as respostas a essa pergunta.

remetente
fonte
1
e adicionando dtype=objectem arr = np.empty( )permitiria a utilização de tipos diferentes no produto, por exemplo arrays = [np.array([1,2,3]), ['str1', 'str2']].
user3820991
Muito obrigado por suas soluções inovadoras. Só pensei que você gostaria de saber que alguns usuários podem encontrar cartesian_product_tranposemais rápido do que cartesian_productdependendo do SO da máquina, python ou versão numpy. Por exemplo, no Ubuntu 14.04, python3.4.3, numpy 1.14.0.dev0 + b7050a9, %timeit cartesian_product_transpose(x500,y500)produz 1000 loops, best of 3: 682 µs per loopenquanto %timeit cartesian_product(x500,y500)produz 1000 loops, best of 3: 1.55 ms per loop. Também estou descobrindo que cartesian_product_transposepode ser mais rápido quando len(arrays) > 2.
Unutbu 17/07
Além disso, cartesian_productretorna uma matriz do tipo de ponto flutuante enquanto cartesian_product_transposeretorna uma matriz do mesmo tipo da primeira matriz (transmitida). A capacidade de preservar o dtype ao trabalhar com matrizes inteiras pode ser um motivo para os usuários favorecerem cartesian_product_transpose.
Unutbu 17/07
@unutbu obrigado novamente - como eu deveria saber, a clonagem do tipo não apenas adiciona conveniência; acelera o código em outros 20 a 30% em alguns casos.
Senderle
1
@ remetente: Uau, isso é bom! Além disso, ocorreu-me que algo como dtype = np.find_common_type([arr.dtype for arr in arrays], [])poderia ser usado para encontrar o dtype comum de todas as matrizes, em vez de forçar o usuário a colocar a matriz que controla o dtype primeiro.
Unutbu 18/07
44

Você pode apenas entender a lista normal em python

x = numpy.array([1,2,3])
y = numpy.array([4,5])
[[x0, y0] for x0 in x for y0 in y]

o que deve lhe dar

[[1, 4], [1, 5], [2, 4], [2, 5], [3, 4], [3, 5]]
ozooxo
fonte
28

Eu também estava interessado nisso e fiz uma pequena comparação de desempenho, talvez um pouco mais clara do que na resposta do @ senderle.

Para duas matrizes (o caso clássico):

insira a descrição da imagem aqui

Para quatro matrizes:

insira a descrição da imagem aqui

(Observe que o comprimento das matrizes é de apenas algumas dezenas de entradas aqui.)


Código para reproduzir as parcelas:

from functools import reduce
import itertools
import numpy
import perfplot


def dstack_product(arrays):
    return numpy.dstack(numpy.meshgrid(*arrays, indexing="ij")).reshape(-1, len(arrays))


# Generalized N-dimensional products
def cartesian_product(arrays):
    la = len(arrays)
    dtype = numpy.find_common_type([a.dtype for a in arrays], [])
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[..., i] = a
    return arr.reshape(-1, la)


def cartesian_product_transpose(arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = reduce(numpy.multiply, broadcasted[0].shape), len(broadcasted)
    dtype = numpy.find_common_type([a.dtype for a in arrays], [])

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T


# from https://stackoverflow.com/a/1235363/577088
def cartesian_product_recursive(arrays, out=None):
    arrays = [numpy.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = numpy.prod([x.size for x in arrays])
    if out is None:
        out = numpy.zeros([n, len(arrays)], dtype=dtype)

    m = n // arrays[0].size
    out[:, 0] = numpy.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian_product_recursive(arrays[1:], out=out[0:m, 1:])
        for j in range(1, arrays[0].size):
            out[j * m : (j + 1) * m, 1:] = out[0:m, 1:]
    return out


def cartesian_product_itertools(arrays):
    return numpy.array(list(itertools.product(*arrays)))


perfplot.show(
    setup=lambda n: 2 * (numpy.arange(n, dtype=float),),
    n_range=[2 ** k for k in range(13)],
    # setup=lambda n: 4 * (numpy.arange(n, dtype=float),),
    # n_range=[2 ** k for k in range(6)],
    kernels=[
        dstack_product,
        cartesian_product,
        cartesian_product_transpose,
        cartesian_product_recursive,
        cartesian_product_itertools,
    ],
    logx=True,
    logy=True,
    xlabel="len(a), len(b)",
    equality_check=None,
)
Nico Schlömer
fonte
17

Com base no trabalho exemplar do @ senderle, criei duas versões - uma para C e outra para layouts de Fortran - que geralmente são um pouco mais rápidas.

  • cartesian_product_transpose_ppé - ao contrário do @ senderle, cartesian_product_transposeque usa uma estratégia completamente diferente - uma versão docartesion_product que usa o layout de memória de transposição mais favorável + algumas otimizações muito menores.
  • cartesian_product_ppfica com o layout de memória original. O que o torna mais rápido é o uso de cópias contíguas. As cópias contíguas acabam sendo muito mais rápidas que copiar um bloco inteiro de memória, embora apenas parte dela contenha dados válidos, é preferível a copiar apenas os bits válidos.

Alguns perfplots. Criei separações para layouts C e Fortran, porque essas são tarefas diferentes da IMO.

Os nomes que terminam em 'pp' são as minhas abordagens.

1) muitos fatores minúsculos (2 elementos cada)

insira a descrição da imagem aquiinsira a descrição da imagem aqui

2) muitos fatores pequenos (4 elementos cada)

insira a descrição da imagem aquiinsira a descrição da imagem aqui

3) três fatores de igual comprimento

insira a descrição da imagem aquiinsira a descrição da imagem aqui

4) dois fatores de igual comprimento

insira a descrição da imagem aquiinsira a descrição da imagem aqui

Código (é necessário executar execuções separadas para cada gráfico b / c Não consegui descobrir como redefinir; também preciso editar / comentar in / out adequadamente):

import numpy
import numpy as np
from functools import reduce
import itertools
import timeit
import perfplot

def dstack_product(arrays):
    return numpy.dstack(
        numpy.meshgrid(*arrays, indexing='ij')
        ).reshape(-1, len(arrays))

def cartesian_product_transpose_pp(arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty((la, *map(len, arrays)), dtype=dtype)
    idx = slice(None), *itertools.repeat(None, la)
    for i, a in enumerate(arrays):
        arr[i, ...] = a[idx[:la-i]]
    return arr.reshape(la, -1).T

def cartesian_product(arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[...,i] = a
    return arr.reshape(-1, la)

def cartesian_product_transpose(arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
    dtype = numpy.result_type(*arrays)

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

from itertools import accumulate, repeat, chain

def cartesian_product_pp(arrays, out=None):
    la = len(arrays)
    L = *map(len, arrays), la
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty(L, dtype=dtype)
    arrs = *accumulate(chain((arr,), repeat(0, la-1)), np.ndarray.__getitem__),
    idx = slice(None), *itertools.repeat(None, la-1)
    for i in range(la-1, 0, -1):
        arrs[i][..., i] = arrays[i][idx[:la-i]]
        arrs[i-1][1:] = arrs[i]
    arr[..., 0] = arrays[0][idx]
    return arr.reshape(-1, la)

def cartesian_product_itertools(arrays):
    return numpy.array(list(itertools.product(*arrays)))


# from https://stackoverflow.com/a/1235363/577088
def cartesian_product_recursive(arrays, out=None):
    arrays = [numpy.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = numpy.prod([x.size for x in arrays])
    if out is None:
        out = numpy.zeros([n, len(arrays)], dtype=dtype)

    m = n // arrays[0].size
    out[:, 0] = numpy.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian_product_recursive(arrays[1:], out=out[0:m, 1:])
        for j in range(1, arrays[0].size):
            out[j*m:(j+1)*m, 1:] = out[0:m, 1:]
    return out

### Test code ###
if False:
  perfplot.save('cp_4el_high.png',
    setup=lambda n: n*(numpy.arange(4, dtype=float),),
                n_range=list(range(6, 11)),
    kernels=[
        dstack_product,
        cartesian_product_recursive,
        cartesian_product,
#        cartesian_product_transpose,
        cartesian_product_pp,
#        cartesian_product_transpose_pp,
        ],
    logx=False,
    logy=True,
    xlabel='#factors',
    equality_check=None
    )
else:
  perfplot.save('cp_2f_T.png',
    setup=lambda n: 2*(numpy.arange(n, dtype=float),),
    n_range=[2**k for k in range(5, 11)],
    kernels=[
#        dstack_product,
#        cartesian_product_recursive,
#        cartesian_product,
        cartesian_product_transpose,
#        cartesian_product_pp,
        cartesian_product_transpose_pp,
        ],
    logx=True,
    logy=True,
    xlabel='length of each factor',
    equality_check=None
    )
Paul Panzer
fonte
Obrigado por compartilhar esta excelente resposta. Quando o tamanho de arraysem cartesian_product_transpose_pp (matrizes) exceder um determinado tamanho, MemoryErrorocorrerá. Nesta situação, eu gostaria que essa função produzisse pequenos pedaços de resultados. Eu postei uma pergunta sobre esse assunto. Você pode responder à minha pergunta? Obrigado.
Sun Bear
13

A partir de outubro de 2017, o numpy agora tem uma np.stackfunção genérica que usa um parâmetro de eixo. Utilizando-o, podemos ter um "produto cartesiano generalizado" usando a técnica "dstack and meshgrid":

import numpy as np
def cartesian_product(*arrays):
    ndim = len(arrays)
    return np.stack(np.meshgrid(*arrays), axis=-1).reshape(-1, ndim)

Nota sobre o axis=-1parâmetro Este é o último eixo (mais interno) do resultado. É equivalente a usar axis=ndim.

Outro comentário: como os produtos cartesianos explodem muito rapidamente, a menos que precisemos realizar a matriz na memória por algum motivo, se o produto for muito grande, convém usar itertoolse usar os valores imediatamente.

MrDrFenner
fonte
8

Eu usei a resposta @kennytm por um tempo, mas ao tentar fazer o mesmo no TensorFlow, mas descobri que o TensorFlow não tem equivalente numpy.repeat(). Após um pouco de experimentação, acho que encontrei uma solução mais geral para vetores arbitrários de pontos.

Para numpy:

import numpy as np

def cartesian_product(*args: np.ndarray) -> np.ndarray:
    """
    Produce the cartesian product of arbitrary length vectors.

    Parameters
    ----------
    np.ndarray args
        vector of points of interest in each dimension

    Returns
    -------
    np.ndarray
        the cartesian product of size [m x n] wherein:
            m = prod([len(a) for a in args])
            n = len(args)
    """
    for i, a in enumerate(args):
        assert a.ndim == 1, "arg {:d} is not rank 1".format(i)
    return np.concatenate([np.reshape(xi, [-1, 1]) for xi in np.meshgrid(*args)], axis=1)

e para o TensorFlow:

import tensorflow as tf

def cartesian_product(*args: tf.Tensor) -> tf.Tensor:
    """
    Produce the cartesian product of arbitrary length vectors.

    Parameters
    ----------
    tf.Tensor args
        vector of points of interest in each dimension

    Returns
    -------
    tf.Tensor
        the cartesian product of size [m x n] wherein:
            m = prod([len(a) for a in args])
            n = len(args)
    """
    for i, a in enumerate(args):
        tf.assert_rank(a, 1, message="arg {:d} is not rank 1".format(i))
    return tf.concat([tf.reshape(xi, [-1, 1]) for xi in tf.meshgrid(*args)], axis=1)
Sean McVeigh
fonte
6

O pacote Scikit-learn possui uma rápida implementação exatamente disso:

from sklearn.utils.extmath import cartesian
product = cartesian((x,y))

Observe que a convenção desta implementação é diferente do que você deseja, se você se importa com a ordem da saída. Para o seu pedido exato, você pode fazer

product = cartesian((y,x))[:, ::-1]
jmd_dk
fonte
Isso é mais rápido que a função do @ senderle?
cs95
@ cᴏʟᴅsᴘᴇᴇᴅ Eu não testei. Eu esperava que isso fosse implementado, por exemplo, C ou Fortran e, portanto, praticamente imbatível, mas parece ter sido escrito usando o NumPy. Como tal, essa função é conveniente, mas não deve ser significativamente mais rápida do que se pode construir usando o NumPy.
jmd_dk
4

De maneira mais geral, se você possui duas matrizes 2d numpy aeb, e deseja concatenar todas as linhas de a a cada linha de b (um produto cartesiano de linhas, como uma junção em um banco de dados), você pode usar esse método :

import numpy
def join_2d(a, b):
    assert a.dtype == b.dtype
    a_part = numpy.tile(a, (len(b), 1))
    b_part = numpy.repeat(b, len(a), axis=0)
    return numpy.hstack((a_part, b_part))
Jonathan
fonte
3

O mais rápido que você pode obter é combinando uma expressão de gerador com a função map:

import numpy
import datetime
a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = (item for sublist in [list(map(lambda x: (x,i),a)) for i in b] for item in sublist)

print (list(foo))

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Saídas (na verdade, toda a lista resultante é impressa):

[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.253567 s

ou usando uma expressão de gerador duplo:

a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = ((x,y) for x in a for y in b)

print (list(foo))

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Saídas (lista completa impressa):

[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.187415 s

Leve em consideração que a maior parte do tempo de computação é direcionada ao comando de impressão. Os cálculos do gerador são decentemente eficientes. Sem imprimir, os tempos de cálculo são:

execution time: 0.079208 s

para expressão de gerador + função de mapa e:

execution time: 0.007093 s

para a expressão de gerador duplo.

Se o que você realmente deseja é calcular o produto real de cada um dos pares de coordenadas, o mais rápido é resolvê-lo como um produto matricial numpy:

a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = np.dot(np.asmatrix([[i,0] for i in a]), np.asmatrix([[i,0] for i in b]).T)

print (foo)

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Saídas:

 [[     0      0      0 ...,      0      0      0]
 [     0      1      2 ...,    197    198    199]
 [     0      2      4 ...,    394    396    398]
 ..., 
 [     0    997   1994 ..., 196409 197406 198403]
 [     0    998   1996 ..., 196606 197604 198602]
 [     0    999   1998 ..., 196803 197802 198801]]
execution time: 0.003869 s

e sem impressão (neste caso, não economiza muito, pois apenas um pequeno pedaço da matriz é realmente impresso):

execution time: 0.003083 s
mosegui
fonte
Para o cálculo do produto, a transmissão externa do produto foo = a[:,None]*bé mais rápida. Usando seu método de temporização sem print(foo), é 0,001103 s vs 0,002225 s. Usando timeit, são 304 μs vs 1,6 ms. Sabe-se que o Matrix é mais lento que o ndarray, então tentei seu código com o np.array, mas ainda é mais lento (1,57 ms) do que a transmissão.
syockit 9/03
2

Isso também pode ser feito facilmente usando o método itertools.product

from itertools import product
import numpy as np

x = np.array([1, 2, 3])
y = np.array([4, 5])
cart_prod = np.array(list(product(*[x, y])),dtype='int32')

Resultado: matriz ([
[1, 4],
[1, 5],
[2, 4],
[2, 5],
[3, 4],
[3, 5]], dtype = int32)

Tempo de execução: 0.000155 s

Muhammad Umar Amanat
fonte
1
você não precisa ligar para entorpecido. matrizes python antigas simples também funcionam com isso.
Coddy 22/01
0

No caso específico em que você precisa executar operações simples, como adição em cada par, é possível introduzir uma dimensão extra e permitir que a transmissão faça o trabalho:

>>> a, b = np.array([1,2,3]), np.array([10,20,30])
>>> a[None,:] + b[:,None]
array([[11, 12, 13],
       [21, 22, 23],
       [31, 32, 33]])

Não tenho certeza se existe alguma maneira similar de obter os pares.

Caagr98
fonte
Se dtypeé floatque você pode fazer (a[:, None, None] + 1j * b[None, :, None]).view(float)que é surpreendentemente rápido.
Paul Panzer