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 una matriz NumPy que contiene todos los posibles [i,j,k] combinaciones de índices.
A partir de ahí, es una simple cuestión de indexar en A utilizando todas las posibles combinaciones de índices, luego reduciendo a lo largo de los ejes apropiados.
Esta solución consume mucha más memoria ya que construye la matriz concreta de todas las combinaciones posibles. A[coms]. Ahorra tiempo para los más pequeños. n, digamos 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.
Trabajar por fragmentos permite combinar la velocidad del cálculo vectorizado y evitar que se produzcan errores de memoria. A continuación, se muestra un ejemplo de conversión de bucles anidados a 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 realizar 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)
Ejecuté 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
- Trabajar por trozos: 1 min 12 s ± 1,6 s
Nota
Después de perfilar el código, encontré que el np.min fue lo que más tardó en llamar np.maximum.reduce. Lo convertí directamente a np.maximum que mejoró un poco el rendimiento.
Puede utilizar combinaciones de itertools, que es una biblioteca estándar de Python y le 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