Usando numpy para construir una matriz de todas las combinaciones de dos matrices

Estoy tratando de ejecutar el espacio de parámetros de una función de 6 parámetros para estudiar su comportamiento numérico antes de intentar hacer algo complejo con él, así que estoy buscando una forma eficiente de hacerlo.

Mi función toma valores flotantes a los que se les asigna una matriz de numpy de 6 dim como entrada. Lo que traté de hacer inicialmente fue esto:

Primero creé una función que toma 2 matrices y genera una matriz con todas las combinaciones de valores de las dos matrices.

from numpy import * def comb(a,b): c = [] for i in a: for j in b: c.append(r_[i,j]) return c 

Luego usé reduce() para aplicar eso a m copias de la misma matriz:

 def combs(a,m): return reduce(comb,[a]*m) 

Y luego evalúo mi función así:

 values = combs(np.arange(0,1,0.1),6) for val in values: print F(val) 

Esto funciona pero es demasiado lento. Sé que el espacio de parámetros es enorme, pero esto no debería ser tan lento. Solo he muestreado 10 6 (un millón) puntos en este ejemplo y me tomó más de 15 segundos crear los values la matriz.

¿Conoces alguna forma más eficiente de hacer esto con numpy?

Puedo modificar la forma en que la función F toma sus argumentos si es necesario.

En la versión más reciente de numpy (> 1.8.x), np.meshgrid proporciona una implementación mucho más rápida:

la solución de @ pv

 In [113]: %timeit cartesian(([1, 2, 3], [4, 5], [6, 7])) 10000 loops, best of 3: 135 µs per loop In [114]: cartesian(([1, 2, 3], [4, 5], [6, 7])) Out[114]: array([[1, 4, 6], [1, 4, 7], [1, 5, 6], [1, 5, 7], [2, 4, 6], [2, 4, 7], [2, 5, 6], [2, 5, 7], [3, 4, 6], [3, 4, 7], [3, 5, 6], [3, 5, 7]]) 

numpy.meshgrid ser solo en 2D, ahora es capaz de ND. En este caso, 3D:

 In [115]: %timeit np.array(np.meshgrid([1, 2, 3], [4, 5], [6, 7])).T.reshape(-1,3) 10000 loops, best of 3: 74.1 µs per loop In [116]: np.array(np.meshgrid([1, 2, 3], [4, 5], [6, 7])).T.reshape(-1,3) Out[116]: array([[1, 4, 6], [1, 5, 6], [2, 4, 6], [2, 5, 6], [3, 4, 6], [3, 5, 6], [1, 4, 7], [1, 5, 7], [2, 4, 7], [2, 5, 7], [3, 4, 7], [3, 5, 7]]) 

Tenga en cuenta que el orden de la resultante final es ligeramente diferente.

Aquí hay una implementación de números puros. Es ca. 5 veces más rápido que el uso de itertools.

 import numpy as np def cartesian(arrays, out=None): """ Generate a cartesian product of input arrays. Parameters ---------- arrays : list of array-like 1-D arrays to form the cartesian product of. out : ndarray Array to place the cartesian product in. Returns ------- out : ndarray 2-D array of shape (M, len(arrays)) containing cartesian products formed of input arrays. Examples -------- >>> cartesian(([1, 2, 3], [4, 5], [6, 7])) array([[1, 4, 6], [1, 4, 7], [1, 5, 6], [1, 5, 7], [2, 4, 6], [2, 4, 7], [2, 5, 6], [2, 5, 7], [3, 4, 6], [3, 4, 7], [3, 5, 6], [3, 5, 7]]) """ arrays = [np.asarray(x) for x in arrays] dtype = arrays[0].dtype n = np.prod([x.size for x in arrays]) if out is None: out = np.zeros([n, len(arrays)], dtype=dtype) m = n / arrays[0].size out[:,0] = np.repeat(arrays[0], m) if arrays[1:]: cartesian(arrays[1:], out=out[0:m,1:]) for j in xrange(1, arrays[0].size): out[j*m:(j+1)*m,1:] = out[0:m,1:] return out 

itertools.combinations es en general la forma más rápida de obtener combinaciones de un contenedor de Python (si de hecho quiere combinaciones, es decir, arreglos SIN repeticiones e independientemente del orden; eso no es lo que su código parece estar haciendo, pero no puedo diga si eso es porque su código tiene errores o porque está usando la terminología incorrecta).

Si desea algo diferente a las combinaciones, quizás otros iteradores en itertools, product o permutations , podrían servirle mejor. Por ejemplo, parece que su código es aproximadamente lo mismo que:

 for val in itertools.product(np.arange(0, 1, 0.1), repeat=6): print F(val) 

Todos estos iteradores producen tuplas, no listas o matrices numpy, por lo que si su F es exigente para obtener específicamente una matriz numpy, tendrá que aceptar la sobrecarga adicional de construir o limpiar y rellenar una en cada paso.

La siguiente implementación numpy debe ser de aprox. 2 veces la velocidad de la respuesta dada:

 def cartesian2(arrays): arrays = [np.asarray(a) for a in arrays] shape = (len(x) for x in arrays) ix = np.indices(shape, dtype=int) ix = ix.reshape(len(arrays), -1).T for n, arr in enumerate(arrays): ix[:, n] = arrays[n][ix[:, n]] return ix 

Parece que quieres una cuadrícula para evaluar tu función, en cuyo caso puedes usar numpy.ogrid (abrir) o numpy.mgrid (completo):

 import numpy my_grid = numpy.mgrid[[slice(0,1,0.1)]*6] 

Puedes hacer algo como esto

 import numpy as np def cartesian_coord(*arrays): grid = np.meshgrid(*arrays) coord_list = [entry.ravel() for entry in grid] points = np.vstack(coord_list).T return points a = np.arange(4) # fake data print(cartesian_coord(*6*[a]) 

lo que da

 array([[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 2], ..., [3, 3, 3, 3, 3, 1], [3, 3, 3, 3, 3, 2], [3, 3, 3, 3, 3, 3]]) 

puede utilizar np.array(itertools.product(a, b))

Aquí hay otra manera, usando NumPy puro, sin recursión, sin comprensión de listas, y no explícito para bucles. Es aproximadamente un 20% más lento que la respuesta original, y se basa en np.meshgrid.

 def cartesian(*arrays): mesh = np.meshgrid(*arrays) # standard numpy meshgrid dim = len(mesh) # number of dimensions elements = mesh[0].size # number of elements, any index will do flat = np.concatenate(mesh).ravel() # flatten the whole meshgrid reshape = np.reshape(flat, (dim, elements)).T # reshape and transpose return reshape 

Por ejemplo,

 x = np.arange(3) a = cartesian(x, x, x, x, x) print(a) 

da

 [[0 0 0 0 0] [0 0 0 0 1] [0 0 0 0 2] ..., [2 2 2 2 0] [2 2 2 2 1] [2 2 2 2 2]] 

Para una implementación numérica pura del producto cartesiano de arrays 1D (o listas de python planas), simplemente use meshgrid() , meshgrid() rodar los ejes con transpose() , y remodelar a la salida deseada:

  def cartprod(*arrays): N = len(arrays) return transpose(meshgrid(*arrays, indexing='ij'), roll(arange(N + 1), -1)).reshape(-1, N) 

Tenga en cuenta que esto tiene la convención de que el último eje cambia más rápido (“estilo C” o “fila principal”).

 In [88]: cartprod([1,2,3], [4,8], [100, 200, 300, 400], [-5, -4]) Out[88]: array([[ 1, 4, 100, -5], [ 1, 4, 100, -4], [ 1, 4, 200, -5], [ 1, 4, 200, -4], [ 1, 4, 300, -5], [ 1, 4, 300, -4], [ 1, 4, 400, -5], [ 1, 4, 400, -4], [ 1, 8, 100, -5], [ 1, 8, 100, -4], [ 1, 8, 200, -5], [ 1, 8, 200, -4], [ 1, 8, 300, -5], [ 1, 8, 300, -4], [ 1, 8, 400, -5], [ 1, 8, 400, -4], [ 2, 4, 100, -5], [ 2, 4, 100, -4], [ 2, 4, 200, -5], [ 2, 4, 200, -4], [ 2, 4, 300, -5], [ 2, 4, 300, -4], [ 2, 4, 400, -5], [ 2, 4, 400, -4], [ 2, 8, 100, -5], [ 2, 8, 100, -4], [ 2, 8, 200, -5], [ 2, 8, 200, -4], [ 2, 8, 300, -5], [ 2, 8, 300, -4], [ 2, 8, 400, -5], [ 2, 8, 400, -4], [ 3, 4, 100, -5], [ 3, 4, 100, -4], [ 3, 4, 200, -5], [ 3, 4, 200, -4], [ 3, 4, 300, -5], [ 3, 4, 300, -4], [ 3, 4, 400, -5], [ 3, 4, 400, -4], [ 3, 8, 100, -5], [ 3, 8, 100, -4], [ 3, 8, 200, -5], [ 3, 8, 200, -4], [ 3, 8, 300, -5], [ 3, 8, 300, -4], [ 3, 8, 400, -5], [ 3, 8, 400, -4]]) 

Si desea cambiar el primer eje más rápido (“estilo FORTRAN” o “columna principal”), simplemente cambie el parámetro de order de reshape() esta manera: reshape((-1, N), order='F')