Producto cartesiano de puntos de matriz x e y en una sola matriz de puntos 2D

Tengo dos matrices numpy que definen los ejes x e y de una cuadrícula. Por ejemplo:

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

Me gustaría generar el producto cartesiano de estas matrices para generar:

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

De una manera que no es terriblemente ineficiente, ya que necesito hacer esto muchas veces en un bucle. Supongo que convertirlos a una lista de Python y usar itertools.producty volver a una matriz numpy no es la forma más 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 Uso de numpy para crear una matriz de todas las combinaciones de dos matrices para obtener una solución general para calcular el producto cartesiano de N matrices.

Un canónico cartesian_product(casi)

Hay muchos enfoques para este problema con diferentes propiedades. Algunos son más rápidos que otros, y algunos son más generales. Después de muchas pruebas y ajustes, descubrí que la siguiente función, que calcula un n-dimensional cartesian_product, es más rápida que la mayoría de las otras entradas. Para un par de enfoques que son un poco más complejos, pero que son incluso un poco más rápidos en muchos casos, vea la respuesta de Paul Panzer .

Dada esa respuesta, esta ya no es la implementación más rápida del producto cartesiano en lo numpyque sé. Sin embargo, creo que su simplicidad continuará haciéndolo un punto de referencia útil para futuras mejoras:

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 la pena mencionar que esta función se usa ix_de una manera inusual; Mientras que el uso documentado de ix_es generar índices en una matriz, resulta que las matrices con la misma forma pueden usarse para la asignación emitida. Muchas gracias a mgilson , que me inspiró a intentar usarlo de ix_esta manera, y a unutbu , que me proporcionó algunos comentarios extremadamente útiles sobre esta respuesta, incluida la sugerencia de uso numpy.result_type.

Alternativas notables

A veces es más rápido escribir bloques de memoria contiguos en el orden Fortran. Esa es la base de esta alternativa, cartesian_product_transposeque ha demostrado ser más rápida en algunos equipos que cartesian_product(ver más abajo). Sin embargo, la respuesta de Paul Panzer, que utiliza el mismo principio, es aún más rápida. Aún así, incluyo esto aquí para los lectores interesados:

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

Después de comprender el enfoque de Panzer, escribí una nueva versión que es casi tan rápida como la suya, y es casi tan simple como 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

Esto parece tener una sobrecarga de tiempo constante que lo hace correr más lento que el de Panzer para entradas pequeñas. Pero para entradas más grandes, en todas las pruebas que ejecuté, funciona tan bien como su implementación más rápida ( cartesian_product_transpose_pp).

En las siguientes secciones, incluyo algunas pruebas de otras alternativas. Ahora están algo desactualizados, pero en lugar de duplicar el esfuerzo, he decidido dejarlos aquí por interés histórico. Para ver las pruebas actualizadas, vea la respuesta de Panzer, así como la de Nico Schlömer .

Pruebas contra alternativas

Aquí hay una batería de pruebas que muestran el aumento de rendimiento que algunas de estas funciones proporcionan en relación con una serie de alternativas. Todas las pruebas que se muestran aquí se realizaron en una máquina de cuatro núcleos, con Mac OS 10.12.5, Python 3.6.1 y numpy1.12.1. Se sabe que las variaciones en hardware y software producen resultados diferentes, por lo que YMMV. ¡Ejecute estas pruebas por sí mismo para estar seguro!

Definiciones:

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

Resultados de la prueba:

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)

En todos los casos, cartesian_productcomo se define al comienzo de esta respuesta, es más rápido.

Para aquellas funciones que aceptan un número arbitrario de matrices de entrada, también vale la pena verificar el rendimiento len(arrays) > 2. (Hasta que pueda determinar por qué cartesian_product_recursivearroja un error en este caso, lo he eliminado de estas pruebas).

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 muestran estas pruebas, cartesian_productsigue siendo competitivo hasta que el número de matrices de entrada supera (aproximadamente) cuatro. Después de eso, cartesian_product_transposetiene una ligera ventaja.

Vale la pena reiterar que los usuarios con otro hardware y sistemas operativos pueden ver resultados diferentes. Por ejemplo, unutbu informa haber visto los siguientes resultados para estas pruebas con Ubuntu 14.04, Python 3.4.3 y 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

A continuación, entro en algunos detalles sobre las pruebas anteriores que he realizado en este sentido. El rendimiento relativo de estos enfoques ha cambiado con el tiempo, para diferentes hardware y diferentes versiones de Python y numpy. Si bien no es inmediatamente útil para las personas que usan versiones actualizadas numpy, ilustra cómo han cambiado las cosas desde la primera versión de esta respuesta.

Una alternativa simple: meshgrid+dstack

La respuesta actualmente aceptada utiliza tiley repeatpara transmitir dos matrices juntas. Pero la meshgridfunción hace prácticamente lo mismo. Aquí está la salida de tiley repeatantes de pasar a transponer:

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

Y aquí está el resultado 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 puede ver, es casi idéntico. Solo necesitamos remodelar el resultado para obtener exactamente el mismo 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])]

Sin embargo, en lugar de remodelar en este punto, podríamos pasar la salida de meshgrida dstacky remodelar después, lo que ahorra algo de trabajo:

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

Contrariamente a lo que se afirma en este comentario , no he visto evidencia de que diferentes entradas produzcan salidas de formas diferentes, y como lo demuestra lo anterior, hacen cosas muy similares, por lo que sería bastante extraño si lo hicieran. Avíseme si encuentra un contraejemplo.

Prueba meshgrid+ dstackvs. repeat+transpose

El rendimiento relativo de estos dos enfoques ha cambiado con el tiempo. En una versión anterior de Python (2.7), el resultado usando meshgrid+ dstackfue notablemente más rápido para entradas pequeñas. (Tenga en cuenta que estas pruebas son de una versión anterior de esta respuesta). Definiciones:

>>> 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 una entrada de tamaño moderado, vi una aceleración significativa. Pero volví a probar estas pruebas con versiones más recientes de Python (3.6.1) y numpy(1.12.1), en una máquina más nueva. Los dos enfoques son casi idénticos ahora.

Prueba antigua

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

Nueva prueba

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 siempre, YMMV, pero esto sugiere que en versiones recientes de Python y numpy, estas son intercambiables.

Funciones de producto generalizadas

En general, podríamos esperar que el uso de funciones integradas sea más rápido para entradas pequeñas, mientras que para entradas grandes, una función especialmente diseñada podría ser más rápida. Además, para un producto n-dimensional generalizado, tiley repeatno ayudará, porque no tienen análogos claros de dimensiones superiores. Por lo tanto, vale la pena investigar el comportamiento de las funciones especialmente diseñadas también.

La mayoría de las pruebas relevantes aparecen al comienzo de esta respuesta, pero aquí hay algunas de las pruebas realizadas en versiones anteriores de Python y numpypara comparación.

La cartesianfunción definida en otra respuesta solía funcionar bastante bien para entradas más grandes. (Es la misma que la función llamada cartesian_product_recursiveanteriormente.) Con el fin de comparar cartesiana dstack_prodct, utilizamos sólo dos dimensiones.

Aquí nuevamente, la prueba anterior mostró una diferencia significativa, mientras que la nueva prueba casi no muestra ninguna.

Prueba antigua

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

Nueva prueba

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_producttodavía late cartesiana escalas más pequeñas.

Nueva prueba ( prueba antigua redundante no mostrada )

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)

Estas distinciones son, creo, interesantes y vale la pena registrarlas; pero son académicos al final. Como lo demostraron las pruebas al comienzo de esta respuesta, todas estas versiones son casi siempre más lentas que las cartesian_productdefinidas al comienzo de esta respuesta, que en sí es un poco más lenta que las implementaciones más rápidas entre las respuestas a esta pregunta.

Puedes hacer una comprensión normal de la lista en Python

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

que debería darte

[[1, 4], [1, 5], [2, 4], [2, 5], [3, 4], [3, 5]]

También estaba interesado en esto e hice una pequeña comparación de rendimiento, tal vez algo más claro que en la respuesta de @ senderle.

Para dos matrices (el caso clásico):

Para cuatro matrices:

(Tenga en cuenta que la longitud de las matrices es de unas pocas docenas de entradas aquí)

Código para reproducir las tramas:

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

Sobre la base del trabajo de campo ejemplar de @ senderle, he creado dos versiones, una para C y otra para diseños de Fortran, que a menudo son un poco más rápidas.

• cartesian_product_transpose_ppes, a diferencia de @ senderle, cartesian_product_transposeque usa una estrategia completamente diferente, una versión cartesion_productque usa el diseño de memoria de transposición más favorable + algunas optimizaciones menores.
• cartesian_product_ppse queda con el diseño de memoria original. Lo que lo hace rápido es usar copias contiguas. Las copias contiguas resultan ser mucho más rápidas que copiar un bloque completo de memoria a pesar de que solo una parte contiene datos válidos es preferible a copiar solo los bits válidos.

Algunas parcelas. Hice otros por separado para los diseños C y Fortran, porque estas son tareas diferentes de la OMI.

Los nombres que terminan en 'pp' son mis enfoques.

1) muchos factores pequeños (2 elementos cada uno)

2) muchos factores pequeños (4 elementos cada uno)

3) tres factores de igual longitud

4) dos factores de igual longitud

Código (necesito hacer corridas separadas para cada parcela b / c No pude encontrar la manera de restablecer; también necesito editar / comentar dentro / fuera adecuadamente):

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 octubre de 2017, numpy ahora tiene una np.stackfunción genérica que toma un parámetro de eje. Al usarlo, podemos tener un "producto cartesiano generalizado" usando la 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 el axis=-1parámetro. Este es el último eje (más interno) en el resultado. Es equivalente a usar axis=ndim.

Otro comentario, dado que los productos cartesianos explotan muy rápidamente, a menos que necesitemos realizar la matriz en la memoria por alguna razón, si el producto es muy grande, es posible que deseemos usar itertoolsy usar los valores sobre la marcha.

Utilicé @kennytm answer por un tiempo, pero cuando intenté hacer lo mismo en TensorFlow, pero descubrí que TensorFlow no tiene equivalente numpy.repeat(). Después de un poco de experimentación, creo que encontré una solución más general para vectores arbitrarios de puntos.

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)

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

El paquete Scikit-learn tiene una rápida implementación de exactamente esto:

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

Tenga en cuenta que la convención de esta implementación es diferente de lo que desea, si le importa el orden de la salida. Para su pedido exacto, puede hacer

product = cartesian((y,x))[:, ::-1]

En términos más generales, si tiene dos matrices nudosas 2d a y b, y desea concatenar cada fila de a cada fila de b (un producto cartesiano de filas, algo así como una unión en una base de datos), puede usar este 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))

Lo más rápido que puede obtener es combinando una expresión generadora con la función de mapa:

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

Salidas (en realidad, se imprime toda la lista resultante):

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

o usando una expresión de doble generador:

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

Salidas (lista completa impresa):

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

Tenga en cuenta que la mayor parte del tiempo de cálculo se destina al comando de impresión. Los cálculos del generador son, por lo demás, decentemente eficientes. Sin imprimir los tiempos de cálculo son:

execution time: 0.079208 s

para la expresión del generador + función de mapa y:

execution time: 0.007093 s

para la expresión de doble generador.

Si lo que realmente desea es calcular el producto real de cada uno de los pares de coordenadas, lo más rápido es resolverlo como un producto de matriz 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()))

Salidas:

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

y sin imprimir (en este caso no ahorra mucho ya que solo se imprime una pequeña parte de la matriz):

execution time: 0.003083 s

Esto también se puede hacer fácilmente usando el 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)

Tiempo de ejecución: 0.000155 s

En el caso específico de que necesite realizar operaciones simples, como la suma de cada par, puede introducir una dimensión adicional y dejar que la transmisión haga el trabajo:

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

No estoy seguro de si hay alguna forma similar de obtener los pares ellos mismos.