Breve introducción a Algoritmos Genéticos - Parte 2

En la primera parte vimos qué eran los algoritmos genéticos y cómo aplicarlos para encontrar el máximo (o el mínimo) de una función. Sin embargo, el ejemplo que dimos era una función continua, diferenciable y con un máximo muy claro:

y = 8x -2x^2

Veremos en este post cómo aplicar el mismo algoritmo pero para buscar el resultado óptimo de una función no continua más compleja, como ser el conjunto de ecuaciones que presentó el Ministro de Educación hace unas semanas:

Definiendo el sistema de ecuaciones

Los algoritmos genéticos intentan aproximarse a la solución de un problema dado. Hacen esto tratando de minimizar o maximizar una función. ¿Cómo podemos representar este problema como algo a minimizar?

Una forma posible es la siguiente:

def f(xyz):

  # calculator => x
  # book => y
  # computer => z

  x,y,z = xyz

  # 1
  s_1 = x + x + x
  s_2 = x + y + y
  s_3 = y - z * 2

  # 2
  err_1 = np.abs(s_1 - 30)
  err_2 = np.abs(s_2 - 18)
  err_3 = np.abs(s_3 - 2)

  # 3
  return np.sum([err_1, err_2, err_3])

1. escribimos las 3 ecuaciones en base a los argumentos x, y, z, que serán los que nuestro algoritmo genético produzca

2. calculamos el "error" de cada ecuación, es decir, qué tan lejos está la suma con los parámetros que generamos del resultado esperado

3. el "fitness" de nuestro x, y, z es la suma de todos los errores. Ésta es la función que intentaremos minimizar

Población inicial

Al igual que en el post anterior creamos una población inicial, pero esta vez:

1. los individuos de la población son tres valores en vez de uno (x, y, z)

2. asumimos que las soluciones están entre 1000 y -1000
   

np.random.seed(1234) # make results reproducible.


def _r():
  return np.random.random() * 2000 - 1000


def random_individual(): 
  return np.array([_r(), _r(), _r()])


initial_generation_size = 1000
first_generation = [random_individual() for _ in range(initial_generation_size)]
pop = Population(chromosomes=first_generation, eval_function=f, maximize=False)

Esta vez seteamos maximize=False porque queremos minimizar, o sea hacer el error lo más pequeño posible.

Definiendo la Evolución

Previamente hicimos todos los pasos de a uno para que sea más fácil entender el proceso pero la librería evol tiene una función para definir de forma más concisa el pipeline de evolución, el objeto Evolution:

def pick_parents(population):
  mom, dad = np.random.choice(population, size=2)
  return mom, dad


def make_child(mom, dad):
  mom_p = np.random.random() # percent of mom contribution
  dad_p = 1 - mom_p # percent of dad contribution
  mom_x, mom_y, mom_z = mom
  dad_x, dad_y, dad_z = dad
  x = mom_x * mom_p + dad_x * dad_p
  y = mom_y * mom_p + dad_y * dad_p
  z = mom_z * mom_p + dad_z * dad_p
  return np.array([x, y, z])


# add gaussian noise with mean zero.
def random_noise(individual, sigma):
  noise = np.random.normal(loc=0, scale=sigma, size=3)
  return individual + noise


best_of = []

def add_best_of(population):
  b = population.current_best
  best_of.append((b.chromosome, f(b.chromosome)))

evo = (Evolution()
      .survive(fraction=0.1)
      .breed(parent_picker=pick_parents, combiner=make_child)
      .mutate(mutate_function=random_noise, sigma=2)
      .callback(add_best_of))

Algunos detalles que valen la pena mencionar:

• la función make_child en vez de tomar el promedio elige un valor al azar entre 0 y 1, toma ese porcentaje de la madre y el restante del padre.

• usamos un array (best_of) y el hook callback del objeto Evolution para ir acumulando los mejores cromosomas de cada generación.

Con el pipeline listo, sólo es cuestión de pasar la población inicial y decirle cuántas generaciones queremos:

new_generation = pop.evolve(evolution=evo, n=25)

Listo! ya tenemos en best_of a los mejores exponentes de cada generación y su valor de "fitness", es decir, qué tanto minimizan el error de nuestra función y por lo tanto encuadran mejor en el sistema de ecuaciones.

Como la última vez, observemos cómo se reduce el error a medida que se van creando nuevas generaciones.

Tomamos entonces como solución los valores de la última generación:

x,y,z = best_of[-1][0]

10.155879542288524, 4.02024292752548, 0.92945755333642

Pero la solución no es la óptima 🤨

No, no lo es, el algoritmo puede no encontrar el mínimo (o máximo) global. En este caso, como indica el Ministro la respuesta correcta era:

x = 10, y = 4, z = 1.

Sin embargo, a diferencia de este sistema sencillo de ecuaciones, hay otros donde la respuesta 100% correcta no es obvia.

Nuestro algoritmo puede generalizar bien a estos otros problemas como veremos a continuación.

Modificando un poco la ecuación

Como muchos conocerán, este enunciado fue blanco de críticas por quienes pensaban que la última ecuación no es y - z * 2 sino y - z^z:

Para resolver esta variante sólo tenemos que cambiar la función a minimizar y correr el programa otra vez:

def f2(xyz):

  # calculator => x
  # book => y
  # computer => z

  x,y,z = xyz

  # 1
  s_1 = x + x + x
  s_2 = x + y + y
  s_3 = y - z ** z # no longer y - 2 * z

  # 2
  err_1 = np.abs(s_1 - 30)
  err_2 = np.abs(s_2 - 18)
  err_3 = np.abs(s_3 - 2)

  # 3
  return np.sum([err_1, err_2, err_3])

initial_generation_size = 1000
first_generation = [random_individual() for _ in range(initial_generation_size)]
pop2 = Population(chromosomes=first_generation, eval_function=f2, maximize=False)

new_generation = pop2.evolve(evolution=evo, n=100)

Analizamos los errores:

Como vemos, los errores no convergen esta vez de forma tan directa. Sin embargo muchas observaciones están con un error bajo. Nos quedamos con los valores que tienen el mejor fit:

fits = [b[1] for b in best_of]
best_fit = np.argmin(fits)
x,y,z = best_of[best_fit][0]

x = 10.60, y = 9.13, z = 2.27
error: 15.27

Un error de "solamente" 15 puntos.

Conclusiones

• Pudimos adaptar nuestro algoritmo de optimización a una nueva función de manera trivial

• Estos algoritmos no llegan siempre al óptimo global sino que intentan aproximarse a una solución "good enough"

• Es posible iterar la solución. En ambos casos al ver los resultados nos damos cuenta de que nuestra intuición general del rango de valores óptimos (-1.000, 1.000) era exagerada. Podríamos volver a intentar con un rango más razonable, de hecho es posible llegar a un error muy bajo, de casi 1 punto, tomando un rango más acotado (esto queda como ejercicio para el lector)