Al fin después de mucho luchar pudimos hallar el resultado de esta aprieto que muchos usuarios de esta web tienen. Si quieres compartir algún detalle no dudes en aportar tu comentario.
Solución:
Esta solución es 5 veces más rápida para n=100:
coms = np.fromiter(itertools.combinations(np.arange(n), 3), 'i,i,i').view(('i', 3))
best = A[coms].min(1).max(1)
at = best.argmin()
global_best = best[at]
save_rows = coms[at]
La primera línea es un poco complicada pero convierte el resultado de itertools.combinations en un numpy array que contiene todo lo posible [i,j,k] combinaciones de índices.
A partir de ahí, es una simple cuestión de indexar en A usando todas las combinaciones de índices posibles, luego reduciendo a lo largo de los ejes apropiados.
Esta solución consume mucha más memoria a medida que construye el concreto array de todas las combinaciones posibles A[coms]. Ahorra tiempo para los más pequeños. ndigamos menos de 250, pero para grandes n el tráfico de memoria será muy alto y puede ser más lento que el código original.
El trabajo por fragmentos permite combinar la velocidad del cálculo vectorizado y evitar los errores de memoria. A continuación hay un ejemplo de cómo convertir los bucles anidados en vectorización por fragmentos.
A partir de las mismas variables que la pregunta, se define una longitud de fragmento para vectorizar los cálculos dentro del fragmento y hacer un bucle solo sobre fragmentos en lugar de sobre combinaciones.
chunk = 2000 # define chunk length, if to small, the code won't take advantage
# of vectorization, if it is too large, excessive memory usage will
# slow down execution, or Memory Error will be risen
combinations = itertools.combinations(range(n),3) # generate iterator containing
# all possible combinations of 3 columns
N = n*(n-1)*(n-2)//6 # number of combinations (length of combinations cannot be
# retrieved because it is an iterator)
# generate a list containing how many elements of combinations will be retrieved
# per iteration
n_chunks, remainder = divmod(N,chunk)
counts_list = [chunk for _ in range(n_chunks)]
if remainder:
counts_list.append(remainder)
# Iterate one chunk at a time, using vectorized code to treat the chunk
for counts in counts_list:
# retrieve combinations in current chunk
current_comb = np.fromiter(combinations,dtype='i,i,i',count=counts)
.view(('i',3))
# maximum of element-wise minimum in current chunk
chunk_best = np.minimum(np.minimum(A[current_comb[:,0],:],A[current_comb[:,1],:]),
A[current_comb[:,2],:]).max(axis=1)
ravel_save_row = chunk_best.argmin() # minimum of maximums in current chunk
# check if current chunk contains global minimum
if chunk_best[ravel_save_row] < global_best:
global_best = chunk_best[ravel_save_row]
save_rows = current_comb[ravel_save_row]
print(global_best,save_rows)
Realicé algunas comparaciones de rendimiento con los bucles anidados, obteniendo los siguientes resultados (chunk_length = 1000):
- n=100
- Bucles anidados: 1,13 s ± 16,6 ms
- Trabajo por trozos: 108 ms ± 565 µs
- n=150
- Bucles anidados: 4,16 s ± 39,3 ms
- Trabajo por trozos: 523 ms ± 4,75 ms
- n=500
- Bucles anidados: 3 min 18 s ± 3,21 s
- Trabajo por bloques: 1min 12s ± 1,6 s
Nota
Después de perfilar el código, descubrí que el np.min fue lo que más tardó en llamar np.maximum.reduce. Lo convertí directamente a np.maximum lo que mejoró un poco el rendimiento.
Puedes usar combinaciones de itertoolsque es una biblioteca estándar de Python, y lo ayudará a eliminar todos esos bucles anidados.
from itertools import combinations
import numpy as np
n = 100
np.random.seed(2)
A = np.random.rand(n,n)
global_best = 1000000000000000.0
for i, j, k in combinations(range(n), 3):
local_best = np.amax(np.array([A[i,:], A[j,:], A[k,:]]).min(0))
if local_best < global_best:
global_best = local_best
save_rows = [i, j, k]
print global_best, save_rows
Comentarios y puntuaciones
Si entiendes que te ha resultado útil nuestro artículo, sería de mucha ayuda si lo compartieras con otros juniors así nos ayudas a extender esta información.