EN
$ cat ~/articles/optimizacion-multiobjetivo-pareto.md ← volver a artículos

Diseñas un sistema de caché distribuido. Quieres baja latencia, alta consistencia y disponibilidad máxima. El teorema CAP te dice que solo puedes elegir dos. Entrenas una red neuronal para producción. Necesitas alta precisión, baja latencia de inferencia y modelo pequeño para edge devices. No puedes tener los tres simultáneamente. Optimizas una API REST. Quieres bajo tiempo de respuesta, alto throughput y uso mínimo de memoria. Cada objetivo conflictúa con al menos otro.

La optimización multiobjetivo formaliza estos problemas donde mejorar un objetivo necesariamente empeora otro. No existe "la mejor solución" única. Existe un conjunto de soluciones Pareto-óptimas donde no puedes mejorar ningún objetivo sin sacrificar al menos otro. Este artículo implementa las métricas fundamentales para trabajar con estas soluciones: dominancia de Pareto, frontera óptima, hipervolumen, spacing, coverage y crowding distance.

Por qué no existe la mejor solución

Un problema de optimización clásico tiene un objetivo único que maximizar o minimizar. Encontrar el camino más corto entre dos puntos en un grafo. Minimizar el error cuadrático medio de un modelo. Maximizar el throughput de un pipeline de procesamiento. Estos problemas tienen solución óptima única o un conjunto de soluciones con el mismo valor objetivo.

Los problemas multiobjetivo tienen múltiples funciones objetivo que queremos optimizar simultáneamente. Formalmente, dado un vector de decisión x en el espacio de búsqueda, queremos optimizar:

minimizar f₁(x), f₂(x), ..., fₖ(x)

Donde cada fᵢ es una función objetivo diferente. Escribimos "minimizar" por convención, pero maximizar f(x) equivale a minimizar -f(x).

El problema surge cuando los objetivos están en conflicto. Mejorar f₁ empeora f₂. Un ejemplo concreto: diseñar un compresor de imágenes.

  • Objetivo 1: Minimizar tamaño del archivo comprimido
  • Objetivo 2: Minimizar pérdida de calidad (medida con PSNR)

Aumentar la tasa de compresión reduce el tamaño pero degrada la calidad. Preservar calidad perfecta requiere compresión lossless con tamaños mayores. No existe ninguna configuración que simultáneamente minimice ambos objetivos. La solución no es un punto, es una curva de trade-offs.

Otro ejemplo real: configurar un cluster de Kubernetes.

  • Objetivo 1: Minimizar latencia p99 de requests
  • Objetivo 2: Minimizar coste de infraestructura (número de nodos)
  • Objetivo 3: Maximizar utilización de CPU (evitar desperdicio)

Escalar horizontalmente reduce latencia pero aumenta coste y reduce utilización por nodo. Reducir nodos aumenta utilización pero puede saturar CPU causando latencias altas. Cada configuración (número de nodos, límites de recursos, políticas de autoscaling) representa un punto diferente en el espacio de objetivos.

Dominancia de Pareto: comparar sin ponderar

Ejemplo: Elegir servidor cloud

Tienes que elegir un tipo de instancia para tu aplicación. Dos opciones:

  • Servidor A: 100ms latencia, 50€/mes
  • Servidor B: 150ms latencia, 80€/mes

El servidor A domina a B porque es mejor en ambos aspectos: más rápido Y más barato. Nadie elegiría B racionalmente. Pero si comparas:

  • Servidor A: 100ms latencia, 50€/mes
  • Servidor C: 80ms latencia, 70€/mes

Ninguno domina al otro. A es más barato pero C es más rápido. Dependes de tus prioridades: si el presupuesto es crítico eliges A, si la velocidad es vital eliges C. No hay respuesta única correcta.


Cuando tienes objetivos en conflicto, no puedes simplemente sumarlos. La expresión latencia_ms + coste_eur no tiene sentido dimensional. Sumar milisegundos con euros es inválido. Ponderar con pesos arbitrarios w₁·latencia + w₂·coste requiere decidir de antemano la importancia relativa, lo cual es exactamente lo que queremos evitar.

La dominancia de Pareto define una relación de orden parcial entre soluciones sin necesidad de ponderar objetivos. Una solución a domina a otra solución b (escrito ab) si y solo si:

  1. a es igual o mejor que b en todos los objetivos: fᵢ(a) ≤ fᵢ(b) para todo i
  2. a es estrictamente mejor que b en al menos un objetivo: existe algún j donde fⱼ(a) < fⱼ(b)

La primera condición garantiza que a no es peor en ningún aspecto. La segunda garantiza que a es genuinamente mejor, no simplemente igual. Si ambas soluciones tienen exactamente los mismos valores en todos los objetivos, ninguna domina a la otra aunque sean idénticas.

Ejemplos concretos con dos objetivos (minimizar ambos):

Solución A: f₁=10, f₂=20
Solución B: f₁=15, f₂=25

A domina a B porque 10 < 15 y 20 < 25. A es mejor en ambos objetivos.

Solución A: f₁=10, f₂=30
Solución B: f₁=15, f₂=25

Ninguna domina a la otra. A es mejor en f₁ pero peor en f₂. B es mejor en f₂ pero peor en f₁. Estas soluciones son incomparables bajo dominancia de Pareto.

Solución A: f₁=10, f₂=20
Solución B: f₁=10, f₂=25

A domina a B porque ambos son iguales en f₁ pero A es estrictamente mejor en f₂.

Implementación en Python:

import numpy as np
from typing import List

def dominates(a: np.ndarray, b: np.ndarray) -> bool:
    """
    Verifica si la solución a domina a la solución b.
    Asume minimización en todos los objetivos.

    Args:
        a: Array de objetivos de la solución a
        b: Array de objetivos de la solución b

    Returns:
        True si a domina a b, False en caso contrario
    """
    # Verificar que a es igual o mejor en todos los objetivos
    all_better_or_equal = np.all(a <= b)

    # Verificar que a es estrictamente mejor en al menos un objetivo
    any_strictly_better = np.any(a < b)

    return all_better_or_equal and any_strictly_better


# Ejemplo de uso
sol_a = np.array([10, 20])  # f1=10, f2=20
sol_b = np.array([15, 25])  # f1=15, f2=25
sol_c = np.array([10, 30])  # f1=10, f2=30

print(f"A domina B: {dominates(sol_a, sol_b)}")  # True
print(f"A domina C: {dominates(sol_a, sol_c)}")  # True
print(f"C domina A: {dominates(sol_c, sol_a)}")  # False
print(f"B domina C: {dominates(sol_b, sol_c)}")  # False

Esta función asume que todos los objetivos se minimizan. Si tienes objetivos que se maximizan, multiplica esos valores por -1 antes de llamar a la función. Por ejemplo, si f₂ se maximiza, convierte el problema usando f₂' = -f₂ y minimiza f₂'.

Frontera de Pareto: soluciones no dominadas

Ejemplo: Configurar tu base de datos

Pruebas 10 configuraciones diferentes de tu PostgreSQL: diferentes tamaños de buffer, conexiones máximas, políticas de vacuum. Mides queries por segundo y memoria usada para cada una. De esas 10 configuraciones, solo 4 están en la frontera de Pareto:

  • Config 1: 1000 qps, 2GB RAM (extremo velocidad)
  • Config 2: 800 qps, 1.5GB RAM (balance)
  • Config 3: 600 qps, 1GB RAM (balance)
  • Config 4: 400 qps, 0.5GB RAM (extremo memoria)

Las otras 6 configuraciones son dominadas: hay siempre una config que da más qps con igual o menos RAM. Esas 6 puedes descartarlas automáticamente. Solo las 4 de la frontera merecen consideración porque cada una representa un trade-off óptimo diferente.


El conjunto de Pareto es el subconjunto de todas las soluciones candidatas donde ninguna solución es dominada por otra del mismo conjunto. Formalmente:

P = { x ∈ S | no existe y ∈ S tal que y ≺ x }

Donde S es el conjunto de todas las soluciones candidatas. En otras palabras, una solución está en el conjunto de Pareto si y solo si ninguna otra solución la domina.

La frontera de Pareto es la imagen de este conjunto en el espacio de objetivos. Si graficamos todas las soluciones con ejes correspondientes a cada objetivo, la frontera de Pareto es la curva formada por los puntos no dominados.

El algoritmo más directo para calcular el conjunto de Pareto es comparar cada solución contra todas las demás:

def pareto_frontier(solutions: np.ndarray) -> np.ndarray:
    """
    Calcula la frontera de Pareto de un conjunto de soluciones.

    Args:
        solutions: Array de forma (n_solutions, n_objectives)

    Returns:
        Array booleano de tamaño n_solutions donde True indica
        que la solución está en la frontera de Pareto
    """
    n_solutions = solutions.shape[0]
    is_pareto = np.ones(n_solutions, dtype=bool)

    for i in range(n_solutions):
        if not is_pareto[i]:
            # Si ya sabemos que i está dominada, saltarla
            continue

        for j in range(n_solutions):
            if i == j:
                continue

            # Si j domina a i, marcar i como dominada
            if dominates(solutions[j], solutions[i]):
                is_pareto[i] = False
                break

    return is_pareto


# Generar soluciones de ejemplo
np.random.seed(42)
solutions = np.random.rand(100, 2) * 100

# Calcular frontera
pareto_mask = pareto_frontier(solutions)
pareto_solutions = solutions[pareto_mask]

print(f"Total de soluciones: {len(solutions)}")
print(f"Soluciones en frontera de Pareto: {len(pareto_solutions)}")
print(f"Porcentaje: {100 * len(pareto_solutions) / len(solutions):.1f}%")

Este algoritmo tiene complejidad O(n²) donde n es el número de soluciones. Para cada solución i, potencialmente compara contra todas las demás n-1 soluciones. En problemas con decenas de miles de soluciones, esto puede volverse lento.

Una optimización común es ordenar las soluciones por el primer objetivo antes de comparar. Si las soluciones están ordenadas ascendentemente por f₁, cuando procesamos la solución i, solo necesitamos comparar contra soluciones anteriores que tienen f₁ ≤ f₁(i). Esta optimización reduce el número de comparaciones en la práctica, aunque la complejidad asintótica sigue siendo O(n²).

def pareto_frontier_optimized(solutions: np.ndarray) -> np.ndarray:
    """
    Versión optimizada que ordena por el primer objetivo.
    """
    # Ordenar por primer objetivo
    sorted_indices = np.argsort(solutions[:, 0])
    sorted_solutions = solutions[sorted_indices]

    n_solutions = solutions.shape[0]
    is_pareto = np.ones(n_solutions, dtype=bool)

    for i in range(n_solutions):
        if not is_pareto[i]:
            continue

        # Solo comparar contra soluciones anteriores (f1 menor o igual)
        for j in range(i):
            if is_pareto[j] and dominates(sorted_solutions[j], sorted_solutions[i]):
                is_pareto[i] = False
                break

    # Revertir el ordenamiento para retornar en orden original
    original_order = np.empty(n_solutions, dtype=bool)
    original_order[sorted_indices] = is_pareto

    return original_order

Hipervolumen: medir calidad de una frontera

Ejemplo: Comparar dos algoritmos de tuning

Ejecutas dos frameworks de hyperparameter optimization para tu modelo: Optuna y Ray Tune. Ambos te dan fronteras de Pareto optimizando accuracy vs tiempo de entrenamiento. ¿Cuál funcionó mejor?

El hipervolumen responde esto. Defines un punto de referencia "peor aceptable": accuracy 70%, tiempo 1 hora. El hipervolumen mide el área del espacio que cada frontera domina respecto a ese punto de referencia. Si Optuna tiene hipervolumen 0.45 y Ray Tune tiene 0.38, Optuna ganó: cubre más espacio de soluciones buenas.

Es como medir "cuánto terreno conquistó" cada algoritmo en el espacio de soluciones. Más área dominada = mejores soluciones encontradas.


Una vez que tienes una frontera de Pareto, necesitas medir qué tan buena es. Esto es esencial para comparar diferentes algoritmos de optimización o para monitorear convergencia durante la búsqueda.

El hipervolumen (también llamado indicador S o medida de Lebesgue) es la métrica más robusta para evaluar fronteras de Pareto. Mide el volumen del espacio de objetivos que es dominado por la frontera con respecto a un punto de referencia.

Formalmente, dado un conjunto de soluciones P y un punto de referencia r, el hipervolumen es:

HV(P, r) = λ( ⋃_{x∈P} [f(x), r] )

Donde λ es la medida de Lebesgue (volumen) y [f(x), r] es el hiperrectángulo desde el punto f(x) hasta el punto de referencia r.

En dos dimensiones, esto se visualiza fácilmente. Cada solución en la frontera define un rectángulo desde ese punto hasta el punto de referencia. El hipervolumen es el área total cubierta por la unión de estos rectángulos.

El punto de referencia debe ser peor que cualquier solución de interés en todos los objetivos. Típicamente se elige como un punto ligeramente peor que el nadir point (el punto con los peores valores de cada objetivo en la frontera).

Para dos objetivos, el cálculo es eficiente ordenando los puntos:

def hypervolume_2d(pareto_front: np.ndarray, reference_point: np.ndarray) -> float:
    """
    Calcula el hipervolumen de una frontera de Pareto en 2D.

    Args:
        pareto_front: Array de forma (n_points, 2) con puntos en la frontera
        reference_point: Array de forma (2,) con el punto de referencia

    Returns:
        Hipervolumen (área) dominada por la frontera
    """
    # Ordenar por primer objetivo (ascendente)
    sorted_front = pareto_front[np.argsort(pareto_front[:, 0])]

    # Calcular área acumulada
    hypervolume = 0.0
    prev_x = reference_point[0]

    for point in sorted_front:
        # Ancho del rectángulo: desde el punto anterior hasta este punto
        width = prev_x - point[0]

        # Altura del rectángulo: desde este punto hasta el punto de referencia
        height = reference_point[1] - point[1]

        # Acumular área
        hypervolume += width * height

        # Actualizar posición X para el siguiente rectángulo
        prev_x = point[0]

    return hypervolume


# Ejemplo: dos fronteras para comparar
front_a = np.array([[1, 5], [2, 3], [3, 2], [5, 1]])
front_b = np.array([[1, 6], [3, 4], [4, 3], [6, 1]])
ref_point = np.array([7, 7])

hv_a = hypervolume_2d(front_a, ref_point)
hv_b = hypervolume_2d(front_b, ref_point)

print(f"Hipervolumen frontera A: {hv_a:.2f}")
print(f"Hipervolumen frontera B: {hv_b:.2f}")
print(f"Frontera {'A' if hv_a > hv_b else 'B'} domina más espacio")

La propiedad más importante del hipervolumen es la monotonía con respecto a dominancia: si una frontera P' domina completamente a otra frontera P (es decir, para cada punto en P existe al menos un punto en P' que lo domina), entonces HV(P', r) ≥ HV(P, r). Esto hace que el hipervolumen sea una métrica teóricamente sólida para comparar fronteras.

Para problemas con tres o más objetivos, calcular el hipervolumen exacto es computacionalmente costoso (NP-hard). La aproximación más práctica es Monte Carlo:

def hypervolume_monte_carlo(pareto_front: np.ndarray,
                           reference_point: np.ndarray,
                           n_samples: int = 10000) -> float:
    """
    Aproxima el hipervolumen usando Monte Carlo para cualquier número de objetivos.

    Args:
        pareto_front: Array de forma (n_points, n_objectives)
        reference_point: Array de forma (n_objectives,)
        n_samples: Número de puntos aleatorios a samplear

    Returns:
        Estimación del hipervolumen
    """
    n_objectives = pareto_front.shape[1]

    # Encontrar límites del hiperrectángulo de sampleo
    min_bounds = np.min(pareto_front, axis=0)
    max_bounds = reference_point

    # Generar puntos aleatorios uniformemente en el hiperrectángulo
    random_points = np.random.uniform(
        low=min_bounds,
        high=max_bounds,
        size=(n_samples, n_objectives)
    )

    # Contar cuántos puntos son dominados por al menos un punto de la frontera
    dominated_count = 0

    for random_point in random_points:
        # Verificar si algún punto de la frontera domina este punto aleatorio
        is_dominated = False

        for pareto_point in pareto_front:
            # Un punto en la frontera domina al punto aleatorio si es
            # mejor o igual en todos los objetivos
            if np.all(pareto_point <= random_point):
                is_dominated = True
                break

        if is_dominated:
            dominated_count += 1

    # El hipervolumen es la proporción de puntos dominados
    # multiplicado por el volumen total del hiperrectángulo
    total_volume = np.prod(max_bounds - min_bounds)
    hypervolume = (dominated_count / n_samples) * total_volume

    return hypervolume


# Ejemplo con 3 objetivos
front_3d = np.array([
    [1, 5, 3],
    [2, 3, 4],
    [3, 2, 5],
    [4, 1, 6]
])
ref_3d = np.array([7, 7, 7])

hv_3d = hypervolume_monte_carlo(front_3d, ref_3d, n_samples=50000)
print(f"Hipervolumen 3D (Monte Carlo): {hv_3d:.2f}")

El método Monte Carlo converge al valor real con O(1/√n) donde n es el número de muestras. Con 10.000 muestras obtienes aproximadamente 1% de error. Para aplicaciones prácticas donde solo necesitas comparar fronteras relativamente (A es mejor que B), esto es suficientemente preciso.

Spacing: medir distribución uniforme

Ejemplo: Ofrecer opciones al usuario

Tu aplicación de e-commerce optimiza velocidad de entrega vs coste de envío. Encuentras 10 soluciones Pareto-óptimas, pero 8 están en la zona "muy rápido, muy caro" y solo 2 en "lento, económico". Técnicamente son óptimas, pero ofreces opciones desequilibradas: el usuario que busca equilibrio no tiene opciones intermedias.

El spacing mide esto. Un spacing bajo significa distribución uniforme: hay opciones distribuidas equitativamente a lo largo de toda la frontera. Spacing alto significa clustering: las soluciones están agrupadas dejando huecos. Para UX, quieres spacing bajo: dar al usuario opciones realmente diversas, no 8 variaciones de "express premium".


El hipervolumen mide qué tan cerca está la frontera del óptimo teórico, pero no mide qué tan uniformemente distribuidas están las soluciones a lo largo de la frontera. Una frontera puede tener hipervolumen alto pero todas las soluciones concentradas en una región, dejando otras regiones sin cobertura.

El spacing (espaciado) cuantifica la uniformidad de la distribución. Mide la desviación estándar de las distancias al vecino más cercano de cada punto:

S = √( (1/n) Σᵢ(dᵢ - d̄)² )

Donde:

  • n es el número de puntos en la frontera
  • dᵢ es la distancia mínima desde el punto i a cualquier otro punto
  • d̄ es la media de todas las distancias mínimas

Un spacing bajo indica distribución uniforme (todas las soluciones están aproximadamente a la misma distancia de sus vecinos). Un spacing alto indica clustering (algunas soluciones están muy juntas y otras muy separadas).

from scipy.spatial import distance_matrix

def spacing(pareto_front: np.ndarray) -> float:
    """
    Calcula la métrica de spacing (uniformidad de distribución).

    Args:
        pareto_front: Array de forma (n_points, n_objectives)

    Returns:
        Valor de spacing (menor es mejor = más uniforme)
    """
    n_points = len(pareto_front)

    if n_points < 2:
        return 0.0

    # Calcular matriz de distancias entre todos los puntos
    dist_matrix = distance_matrix(pareto_front, pareto_front)

    # Poner infinito en la diagonal para ignorar distancia a sí mismo
    np.fill_diagonal(dist_matrix, np.inf)

    # Encontrar distancia mínima para cada punto
    min_distances = np.min(dist_matrix, axis=1)

    # Calcular desviación estándar de estas distancias
    mean_dist = np.mean(min_distances)
    variance = np.mean((min_distances - mean_dist) ** 2)
    spacing_value = np.sqrt(variance)

    return spacing_value


# Comparar dos fronteras: una uniforme y una con clustering
uniform_front = np.array([
    [1, 5], [2, 4], [3, 3], [4, 2], [5, 1]
])

clustered_front = np.array([
    [1, 5], [1.1, 4.9], [1.2, 4.8],  # Cluster en un extremo
    [4, 2], [5, 1]                     # Puntos aislados en el otro
])

print(f"Spacing frontera uniforme: {spacing(uniform_front):.4f}")
print(f"Spacing frontera con clustering: {spacing(clustered_front):.4f}")

En aplicaciones prácticas, quieres spacing bajo porque significa que estás dando al usuario opciones diversas. Si todas tus soluciones están en una región pequeña de la frontera, el usuario solo puede elegir entre trade-offs muy similares. Una frontera bien distribuida ofrece opciones desde "optimiza casi exclusivamente objetivo 1" hasta "optimiza casi exclusivamente objetivo 2" con gradaciones intermedias.

Coverage: comparar dos fronteras

Ejemplo: Benchmark de compiladores

Optimizas flags de compilación de GCC vs Clang para tu codebase: binario más pequeño vs ejecución más rápida. GCC produce 12 configuraciones Pareto-óptimas, Clang produce 15. ¿Cuál es mejor?

Coverage(GCC, Clang) = 0.80 significa que el 80% de las soluciones de Clang son dominadas por al menos una de GCC. Coverage(Clang, GCC) = 0.25 significa que solo el 25% de las de GCC son dominadas por Clang. Conclusión: GCC encontró mejores soluciones, domina la mayoría de lo que Clang ofrece.

Es como un combate directo: "¿cuántos de tus puntos son peores que al menos uno de los míos?". Alto coverage en una dirección = superioridad clara.


Cuando comparas dos algoritmos de optimización multiobjetivo, necesitas determinar cuál produce fronteras mejores. El hipervolumen te dice cuál domina más espacio, pero no te dice cuántas soluciones de una frontera son dominadas por la otra.

La métrica de coverage (cobertura) C(A, B) mide qué porcentaje de soluciones en B son dominadas por al menos una solución en A:

C(A, B) = |{b ∈ B | ∃a ∈ A : a ≺ b}| / |B|

Esta métrica es asimétrica: C(A, B) no es igual a C(B, A). Si C(A, B) = 1, significa que todas las soluciones de B son dominadas por al menos una solución de A. Si C(A, B) = 0, ninguna solución de B es dominada por A.

def coverage(front_a: np.ndarray, front_b: np.ndarray) -> float:
    """
    Calcula C(A, B): qué porcentaje de B es dominado por A.

    Args:
        front_a: Array de forma (n_a, n_objectives)
        front_b: Array de forma (n_b, n_objectives)

    Returns:
        Valor entre 0 y 1 (mayor significa A domina más de B)
    """
    n_b = len(front_b)
    dominated_count = 0

    for sol_b in front_b:
        # Verificar si alguna solución de A domina esta solución de B
        for sol_a in front_a:
            if dominates(sol_a, sol_b):
                dominated_count += 1
                break  # No necesitamos seguir chequeando otras de A

    return dominated_count / n_b if n_b > 0 else 0.0


# Ejemplo: comparar dos fronteras
front_a = np.array([[1, 5], [2, 3], [3, 2], [5, 1]])
front_b = np.array([[2, 6], [3, 4], [4, 3], [6, 2]])

c_ab = coverage(front_a, front_b)
c_ba = coverage(front_b, front_a)

print(f"C(A, B) = {c_ab:.2f}")  # Qué porcentaje de B domina A
print(f"C(B, A) = {c_ba:.2f}")  # Qué porcentaje de A domina B

if c_ab > c_ba:
    print("La frontera A domina más soluciones de B que viceversa")
elif c_ba > c_ab:
    print("La frontera B domina más soluciones de A que viceversa")
else:
    print("Ambas fronteras tienen cobertura similar")

En benchmarking de algoritmos, típicamente ejecutas cada algoritmo múltiples veces con semillas aleatorias diferentes, obtienes una frontera por ejecución, y calculas C(A, B) entre todas las parejas de fronteras de ambos algoritmos. Luego reportas media y desviación estándar de C(A, B) y C(B, A).

Crowding distance: preservar diversidad en NSGA-II

Ejemplo: Evolucionar arquitecturas de redes neuronales

Usas un algoritmo genético para encontrar arquitecturas de redes neuronales óptimas: accuracy vs número de parámetros. Después de 100 generaciones, todas las soluciones convergen a la misma región: redes de ~50M parámetros con 92-93% accuracy. Perdiste diversidad: no hay opciones de redes pequeñas para mobile ni redes grandes para máxima precisión.

El crowding distance soluciona esto. Asigna un "score de soledad" a cada solución: soluciones en regiones densas reciben score bajo, soluciones aisladas reciben score alto. Durante selección, el algoritmo prefiere soluciones solitarias para reproducción, incentivando explorar regiones menos pobladas de la frontera.

Es como repartir recursos uniformemente en un territorio: si hay mucha gente en una zona, incentivas migración a zonas vacías para mantener cobertura equilibrada del espacio.


El algoritmo NSGA-II (Non-dominated Sorting Genetic Algorithm II) es uno de los algoritmos evolutivos más populares para optimización multiobjetivo. Uno de sus componentes clave es el crowding distance, que estima la densidad de soluciones alrededor de cada punto de la frontera.

El problema que resuelve es el siguiente: durante la evolución, el algoritmo tiende a converger hacia regiones de la frontera donde las soluciones son localmente óptimas según los operadores genéticos. Esto causa clustering, dejando otras regiones sin explorar. El crowding distance penaliza soluciones en regiones densas, incentivando al algoritmo a mantener diversidad.

El crowding distance de una solución se calcula como la suma de las distancias a sus vecinos en cada objetivo, normalizada por el rango del objetivo:

CD(i) = Σⱼ [ (fⱼ(i+1) - fⱼ(i-1)) / (fⱼᵐᵃˣ - fⱼᵐⁱⁿ) ]

Donde:

  • i+1 es el vecino inmediato superior en el objetivo j
  • i-1 es el vecino inmediato inferior en el objetivo j
  • fⱼᵐᵃˣ y fⱼᵐⁱⁿ son los valores máximo y mínimo del objetivo j en toda la frontera

Las soluciones en los extremos de cada objetivo reciben crowding distance infinito porque son únicas (no tienen vecino en esa dirección).

def crowding_distance(pareto_front: np.ndarray) -> np.ndarray:
    """
    Calcula el crowding distance para cada punto en la frontera de Pareto.

    Args:
        pareto_front: Array de forma (n_points, n_objectives)

    Returns:
        Array de crowding distances de tamaño n_points
    """
    n_points, n_objectives = pareto_front.shape

    # Inicializar todas las distancias a 0
    distances = np.zeros(n_points)

    # Calcular para cada objetivo
    for obj_idx in range(n_objectives):
        # Ordenar por este objetivo
        sorted_indices = np.argsort(pareto_front[:, obj_idx])

        # Extraer valores de este objetivo ordenados
        obj_values = pareto_front[sorted_indices, obj_idx]

        # Rango del objetivo (para normalización)
        obj_range = obj_values[-1] - obj_values[0]

        # Casos especiales: si todos los valores son iguales, el rango es 0
        if obj_range == 0:
            continue

        # Asignar distancia infinita a los extremos
        distances[sorted_indices[0]] = np.inf
        distances[sorted_indices[-1]] = np.inf

        # Calcular distancia para puntos intermedios
        for i in range(1, n_points - 1):
            # Distancia al vecino de arriba menos distancia al vecino de abajo
            distance_contribution = (obj_values[i + 1] - obj_values[i - 1]) / obj_range

            # Acumular distancia (sumar contribuciones de todos los objetivos)
            distances[sorted_indices[i]] += distance_contribution

    return distances


# Ejemplo: frontera con distribución no uniforme
front = np.array([
    [1, 10],   # Extremo objetivo 1
    [2, 8],
    [2.1, 7.9],  # Muy cerca del anterior (debería tener CD bajo)
    [5, 5],
    [10, 1]    # Extremo objetivo 2
])

cd = crowding_distance(front)

print("Crowding distances:")
for i, (point, distance) in enumerate(zip(front, cd)):
    print(f"Punto {i}: {point} -> CD = {distance:.3f if distance != np.inf else 'inf'}")

print("\nPuntos ordenados por crowding distance (de más aislado a más denso):")
sorted_indices = np.argsort(cd)[::-1]
for idx in sorted_indices:
    cd_str = "inf" if cd[idx] == np.inf else f"{cd[idx]:.3f}"
    print(f"Punto {idx}: {front[idx]} -> CD = {cd_str}")

En NSGA-II, durante la selección de padres para reproducción, cuando dos soluciones están en el mismo frente de dominancia (mismo rank), se prefiere la que tiene mayor crowding distance. Esto incentiva explorar regiones menos pobladas de la frontera.

El algoritmo completo de NSGA-II combina:

  1. Non-dominated sorting: clasificar soluciones en múltiples frentes (frente 1 = frontera de Pareto, frente 2 = frontera de Pareto del resto, etc.)
  2. Crowding distance: calcular densidad dentro de cada frente
  3. Selección por torneo: comparar soluciones por rank primero, por crowding distance si tienen mismo rank
  4. Operadores genéticos estándar: crossover y mutación

Implementación completa: clase reutilizable

Juntando todas las métricas en una clase cohesiva para usar en proyectos reales:

from dataclasses import dataclass
from typing import Tuple
import numpy as np
from scipy.spatial import distance_matrix


@dataclass
class ParetoMetrics:
    """Resultados de análisis de frontera de Pareto."""
    hypervolume: float
    spacing: float
    n_points: int
    crowding_distances: np.ndarray


class MultiObjectiveProblem:
    """
    Clase para trabajar con problemas de optimización multiobjetivo.
    """

    def __init__(self, solutions: np.ndarray, reference_point: np.ndarray = None):
        """
        Args:
            solutions: Array (n_solutions, n_objectives)
            reference_point: Punto de referencia para hipervolumen.
                           Si es None, se calcula automáticamente.
        """
        self.solutions = solutions
        self.n_solutions, self.n_objectives = solutions.shape

        if reference_point is None:
            # Punto de referencia: ligeramente peor que el peor en cada objetivo
            self.reference_point = np.max(solutions, axis=0) * 1.1
        else:
            self.reference_point = reference_point

    def get_pareto_front(self) -> Tuple[np.ndarray, np.ndarray]:
        """
        Calcula la frontera de Pareto.

        Returns:
            Tupla de (soluciones en frontera, índices originales)
        """
        pareto_mask = self._pareto_frontier()
        pareto_solutions = self.solutions[pareto_mask]
        pareto_indices = np.where(pareto_mask)[0]

        return pareto_solutions, pareto_indices

    def _pareto_frontier(self) -> np.ndarray:
        """Implementación del algoritmo de frontera de Pareto."""
        is_pareto = np.ones(self.n_solutions, dtype=bool)

        for i in range(self.n_solutions):
            if not is_pareto[i]:
                continue

            for j in range(self.n_solutions):
                if i == j:
                    continue

                if self._dominates(self.solutions[j], self.solutions[i]):
                    is_pareto[i] = False
                    break

        return is_pareto

    @staticmethod
    def _dominates(a: np.ndarray, b: np.ndarray) -> bool:
        """Verifica dominancia de Pareto."""
        return np.all(a <= b) and np.any(a < b)

    def calculate_metrics(self) -> ParetoMetrics:
        """
        Calcula todas las métricas para la frontera de Pareto.

        Returns:
            Objeto ParetoMetrics con todas las métricas
        """
        pareto_front, _ = self.get_pareto_front()

        if self.n_objectives == 2:
            hv = self._hypervolume_2d(pareto_front)
        else:
            hv = self._hypervolume_monte_carlo(pareto_front)

        sp = self._spacing(pareto_front)
        cd = self._crowding_distance(pareto_front)

        return ParetoMetrics(
            hypervolume=hv,
            spacing=sp,
            n_points=len(pareto_front),
            crowding_distances=cd
        )

    def _hypervolume_2d(self, front: np.ndarray) -> float:
        """Hipervolumen exacto para 2 objetivos."""
        sorted_front = front[np.argsort(front[:, 0])]

        hypervolume = 0.0
        prev_x = self.reference_point[0]

        for point in sorted_front:
            width = prev_x - point[0]
            height = self.reference_point[1] - point[1]
            hypervolume += width * height
            prev_x = point[0]

        return hypervolume

    def _hypervolume_monte_carlo(self, front: np.ndarray, n_samples: int = 10000) -> float:
        """Hipervolumen Monte Carlo para 3+ objetivos."""
        min_bounds = np.min(front, axis=0)
        max_bounds = self.reference_point

        random_points = np.random.uniform(
            low=min_bounds,
            high=max_bounds,
            size=(n_samples, self.n_objectives)
        )

        dominated_count = 0
        for random_point in random_points:
            for pareto_point in front:
                if np.all(pareto_point <= random_point):
                    dominated_count += 1
                    break

        total_volume = np.prod(max_bounds - min_bounds)
        return (dominated_count / n_samples) * total_volume

    def _spacing(self, front: np.ndarray) -> float:
        """Calcula spacing."""
        if len(front) < 2:
            return 0.0

        dist_matrix = distance_matrix(front, front)
        np.fill_diagonal(dist_matrix, np.inf)

        min_distances = np.min(dist_matrix, axis=1)
        mean_dist = np.mean(min_distances)
        variance = np.mean((min_distances - mean_dist) ** 2)

        return np.sqrt(variance)

    def _crowding_distance(self, front: np.ndarray) -> np.ndarray:
        """Calcula crowding distance."""
        n_points = len(front)
        distances = np.zeros(n_points)

        for obj_idx in range(self.n_objectives):
            sorted_indices = np.argsort(front[:, obj_idx])
            obj_values = front[sorted_indices, obj_idx]

            obj_range = obj_values[-1] - obj_values[0]
            if obj_range == 0:
                continue

            distances[sorted_indices[0]] = np.inf
            distances[sorted_indices[-1]] = np.inf

            for i in range(1, n_points - 1):
                distance_contribution = (obj_values[i + 1] - obj_values[i - 1]) / obj_range
                distances[sorted_indices[i]] += distance_contribution

        return distances

    @staticmethod
    def compare_fronts(front_a: np.ndarray, front_b: np.ndarray) -> Tuple[float, float]:
        """
        Compara dos fronteras usando coverage.

        Returns:
            Tupla (C(A,B), C(B,A))
        """
        def coverage(fa: np.ndarray, fb: np.ndarray) -> float:
            dominated = 0
            for sb in fb:
                for sa in fa:
                    if MultiObjectiveProblem._dominates(sa, sb):
                        dominated += 1
                        break
            return dominated / len(fb) if len(fb) > 0 else 0.0

        return coverage(front_a, front_b), coverage(front_b, front_a)


# Ejemplo de uso completo
if __name__ == "__main__":
    # Generar problema de prueba
    np.random.seed(42)
    solutions = np.random.rand(200, 2) * 100

    # Crear problema
    problem = MultiObjectiveProblem(solutions)

    # Obtener frontera
    pareto_front, indices = problem.get_pareto_front()
    print(f"Encontradas {len(pareto_front)} soluciones en frontera de Pareto")

    # Calcular métricas
    metrics = problem.calculate_metrics()
    print(f"\nMétricas de la frontera:")
    print(f"  Hipervolumen: {metrics.hypervolume:.2f}")
    print(f"  Spacing: {metrics.spacing:.4f}")
    print(f"  Puntos en frontera: {metrics.n_points}")

    # Mostrar crowding distances
    print(f"\nTop 5 puntos más aislados (mayor crowding distance):")
    sorted_cd = np.argsort(metrics.crowding_distances)[::-1]
    for i in sorted_cd[:5]:
        cd_str = "inf" if metrics.crowding_distances[i] == np.inf else f"{metrics.crowding_distances[i]:.3f}"
        print(f"  {pareto_front[i]} -> CD = {cd_str}")

Esta implementación está lista para usar en proyectos de producción. Incluye manejo de casos especiales (fronteras con un solo punto, objetivos con rango cero, soluciones duplicadas) y usa NumPy eficientemente para operaciones vectorizadas donde es posible.

Aplicaciones prácticas

Las técnicas de optimización multiobjetivo se usan en múltiples dominios:

Machine Learning: Hyperparameter tuning optimizando precisión vs tiempo de entrenamiento vs tamaño del modelo. Neural Architecture Search balanceando accuracy vs latencia vs FLOPs.

Ingeniería de software: Optimización de compiladores balanceando tiempo de compilación vs tamaño del binario vs velocidad de ejecución. Asignación de recursos en clusters balanceando latencia vs coste vs utilización.

Redes neuronales: Pruning de redes optimizando accuracy vs sparsity vs inference speed. Quantization balanceando precisión vs tamaño vs latencia en edge devices.

Sistemas distribuidos: Configuración de caches balanceando hit rate vs memoria usada vs overhead de sincronización. Diseño de pipelines de procesamiento balanceando throughput vs latencia vs coste.

Finanzas: Portfolio optimization maximizando retorno esperado mientras se minimiza riesgo y se respetan restricciones de liquidez.

En todos estos casos, la frontera de Pareto proporciona un conjunto de soluciones óptimas entre las cuales el usuario o sistema puede elegir según preferencias específicas del contexto. Las métricas implementadas aquí permiten cuantificar y comparar la calidad de estas fronteras de manera rigurosa.