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

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.

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

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.

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

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

Para quatro matrizes:

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

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)

2) muitos fatores pequenos (4 elementos cada)

3) três fatores de igual comprimento

4) dois fatores de igual comprimento

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
    )

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.

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)

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]

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

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

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

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.