Caminos mínimos con programación dinámica

La semana pasada vimos que los subproblemas de una programación dinámica y sus dependencias forman un DAG, y llenamos tablas sobre intervalos y sobre árboles. Había un orden natural: primero los intervalos cortos, primero los hijos. Esta semana el problema es de otro tipo. Queremos caminos mínimos en un grafo que puede tener ciclos, así que los nodos no tienen ningún orden natural.

El problema ya lo conocemos. En la semana 5 demostramos Dijkstra, y vimos exactamente en qué línea se rompe la demostración si hay pesos negativos. Dejamos un pendiente: ¿qué se hace entonces? La respuesta es una idea de la semana 7. Si los nodos no dan un orden, se agrega un parámetro que sí lo dé: cuántas aristas se permite usar. Con eso, un problema con ciclos se vuelve una tabla que se llena fila por fila. Eso es Bellman-Ford.

La Sesión A desarrolla Bellman-Ford como programación dinámica, y responde qué pasa cuando hay un ciclo de peso negativo: cómo se detecta y cómo se recupera. La Sesión B calcula caminos mínimos entre todos los pares, con otro parámetro: qué nodos se permite usar como intermedios. Eso es Floyd-Warshall. Termina con un cuadro para elegir algoritmo y, como tema que se puede cortar, una programación dinámica cuyo parámetro es un subconjunto.

Ruta de lectura.

Antes de clase (15–20 minutos): la introducción, Por qué con aristas negativas hace falta otra idea, Sin ciclos negativos bastan V − 1 aristas y Todos los pares: V veces un algoritmo de un origen.

En clase: Sesión A, de Paso 1: el subproblema con a lo sumo i aristas a Empezar en “infinito”: dos errores. Sesión B, Floyd-Warshall: el subproblema con nodos intermedios permitidos, Cuadro de decisión y ✂️ Máscaras de bits y Held-Karp si hay tiempo.

Después de clase: el Resumen y los Ejercicios, para repasar. Antes del laboratorio, Del producto a la suma: los tipos de cambio y El código de la semana. Para ampliar, Dónde se usa: el vector de distancias, ✂️ Máscaras de bits y Held-Karp si no se vio en clase, y las lecturas de apoyo del recuadro siguiente.

Lectura. Sesión A: Erickson, secs. 8.5 (Directed Acyclic Graphs: Depth-First Search) y 8.7 (Relax ALL the Edges: Bellman-Ford). CLRS, capítulo Single-Source Shortest Paths, 4.ª edición. Sesión B: Erickson, cap. 9. CLRS, capítulo All-Pairs Shortest Paths, 4.ª edición. Estos apuntes se pueden estudiar sin ellos.

Código. Los programas de este capítulo están en code/week-10. Dentro del entorno del curso, entren a esa carpeta y corran make para compilar los ejemplos y make test para comprobarlos. Los comandos de este capítulo se corren desde ahí. El problema de laboratorio está en lab/arbitrage y tiene su propio make test.

Por qué con aristas negativas hace falta otra idea

Tenemos un grafo dirigido con V nodos, E aristas, un peso w(u, v) en cada arista, que puede ser negativo, y un origen s. Usamos las palabras de la semana 5. Un camino de s a v es una secuencia de aristas consecutivas, que puede repetir nodos, y su largo es la suma de sus pesos. Un camino simple no repite nodos. δ(v) es el menor largo de un camino de s a v, si ese mínimo existe.

¿Por qué alguien querría un peso negativo? Porque no todo peso es una distancia. Puede ser un costo que a veces es una ganancia: un tramo de un viaje en el que se recibe un pago, una reacción química que libera energía, una operación de cambio de moneda. Al final de este capítulo vamos a ver ese último caso con detalle.

Dijkstra no sirve. Este es el grafo de la semana 5, con el algoritmo de esa semana:

$ build/bf_demo
Aristas: 0->1 (1), 0->2 (4), 2->1 (-5), 1->3 (2)
Dijkstra:        0   -1    4    3

El valor correcto para el nodo 3 es 1, por el camino 0 → 2 → 1 → 3, de largo 4 − 5 + 2. Dijkstra fija el 1 con distancia 1, relaja sus aristas con ese valor y no las vuelve a mirar. Cuando después d[1] baja a −1, el nodo 3 ya no se entera. La idea de fijar un nodo es lo que falla: con un negativo, el futuro todavía puede mejorar lo que ya se fijó. Y en la semana 5 también vimos que sumar una constante a todos los pesos tampoco sirve, porque castiga más a los caminos con más aristas.

Hay además un problema más serio, que Dijkstra nunca tuvo. Un ciclo negativo (negative cycle) es un ciclo cuyos pesos suman menos que cero. Si desde s se puede llegar a un ciclo negativo, y desde él se puede llegar a v, entonces δ(v) no existe: cada vuelta al ciclo acorta el camino, y no hay un mínimo. Por eso todo lo que sigue tiene dos partes: qué calcular cuando no hay ciclos negativos, y cómo darse cuenta de que los hay.

Paso 1: el subproblema con a lo sumo i aristas

Seguimos el método de la semana 7. Los nodos no están ordenados, pero el número de aristas de un camino sí es un parámetro que crece.

Subproblema. d_i(v) es el menor largo de un camino de s a v con a lo sumo i aristas, o ∞ si no hay ninguno.

Con i aristas hay una cantidad finita de caminos, así que este mínimo siempre existe, incluso con ciclos negativos. Esa es la ventaja de la frase: no hace falta suponer nada para que tenga sentido. Más adelante veremos qué fila contiene la respuesta.

Paso 2: la recurrencia

Miremos un camino de s a v con a lo sumo i aristas, y hagámonos la pregunta de siempre: ¿qué pasa con lo último? Hay dos casos. O bien tiene a lo sumo i − 1 aristas, o bien tiene exactamente i, y entonces su última arista es alguna (u, v).

d_0(s) = 0
d_0(v) = ∞                                                        para v ≠ s
d_i(v) = mín( d_(i-1)(v),  mín sobre (u, v) de d_(i-1)(u) + w(u, v) )     para i ≥ 1

Con la convención ∞ + w = ∞, para cualquier w. Eso parece un detalle y es el origen de un error común, como veremos.

Teorema. Para todo i ≥ 0 y todo nodo v, la recurrencia calcula el menor largo de un camino de s a v con a lo sumo i aristas. Vale para pesos de cualquier signo, y con ciclos negativos.

Técnica. Inducción sobre i, con las dos direcciones. Sea OPT_i(v) el valor que dice la frase del subproblema.

Caso base. Con 0 aristas, el único camino es el vacío, de s a s, y mide 0.

Paso. Supongamos que d_(i−1) = OPT_(i−1) en todos los nodos.

(≥) Cada opción de la recurrencia es el largo de un camino permitido. d_(i−1)(v) lo es, porque un camino con a lo sumo i − 1 aristas también tiene a lo sumo i. Y si d_(i−1)(u) es finito, hay un camino de s a u con a lo sumo i − 1 aristas y ese largo, y al agregarle (u, v) resulta uno de s a v con a lo sumo i. Por lo tanto el mínimo de las opciones no baja de lo que un camino permitido puede lograr: OPT_i(v) ≤ d_i(v).

(≤) Tomemos un camino P de s a v con a lo sumo i aristas y largo OPT_i(v), si existe alguno. Si tiene a lo sumo i − 1 aristas, entonces su largo es al menos OPT_(i−1)(v) = d_(i−1)(v). Si tiene exactamente i ≥ 1, sea (u, v) la última arista. Lo que queda es un camino de s a u con i − 1 aristas, de largo al menos d_(i−1)(u), así que P mide al menos d_(i−1)(u) + w(u, v). En los dos casos, d_i(v) ≤ OPT_i(v). Y si no hay ningún camino, ninguna opción es finita, por la hipótesis de inducción, y las dos cantidades son ∞. ∎

Fíjense en lo que no se usó: nunca dijimos que el camino sea simple, ni que los pesos tengan un signo.

Paso 3: el orden y el costo

La fila i solo lee la fila i − 1, así que se llena de arriba hacia abajo. Este es el mismo ejemplo:

La tabla d_i(v):
  i    0    1    2    3
  0    0  inf  inf  inf
  1    0    1    4  inf
  2    0   -1    4    3
  3    0   -1    4    1
  4    0   -1    4    1

Sigamos dos casillas. d_1(1) = 1 usa la arista directa. d_2(1) mira d_1(1) = 1, d_1(0) + 1 = 1 y d_1(2) − 5 = −1, y se queda con −1. Luego d_3(3) = d_2(1) + 2 = 1: aquí aparece el camino que Dijkstra no vio. A partir de la fila 3, nada cambia. Y el 3 aparece con largo 3 en la fila 2 porque con dos aristas lo mejor es 0 → 1 → 3: esa fila no se equivoca, solo responde otra pregunta.

Cada fila revisa cada arista una vez y copia V valores. Con i de 1 a V − 1, el costo es O(VE), si suponemos que todo nodo es alcanzable, para que V ≤ E + 1. Falta justificar por qué V − 1 filas bastan.

Sin ciclos negativos bastan V − 1 aristas

Lema. Si desde s no se alcanza ningún ciclo negativo, entonces para todo nodo v alcanzable, δ(v) existe y es el largo de un camino simple. Ese camino tiene a lo sumo V − 1 aristas, así que δ(v) = d_(V−1)(v).

Técnica. Quitar ciclos de un camino, sin alargarlo.

Tomemos cualquier camino P de s a v que repita un nodo x. Entre las dos apariciones de x hay un ciclo. Ese ciclo se alcanza desde s, porque P pasa por él, así que por hipótesis pesa 0 o más. Si lo quitamos, sigue siendo un camino de s a v, con menos aristas y de largo menor o igual. Repitiendo esto, que termina porque cada vez hay menos aristas, de P queda un camino simple que no mide más que P.

Entonces todo camino tiene un camino simple que lo iguala o lo mejora. Los caminos simples de s a v son una cantidad finita, así que hay uno mínimo, y ese es δ(v). Un camino simple no repite nodos, así que visita a lo sumo V, y por lo tanto tiene a lo sumo V − 1 aristas. Es uno de los que cuenta d_(V−1)(v), de modo que d_(V−1)(v) ≤ δ(v). Y d_(V−1)(v) es el largo de algún camino, así que no baja de δ(v). ∎

Esta es la pregunta de repaso 2, y la hipótesis importa. Un ciclo de peso 0 no estorba: se quita y el largo no cambia. Un ciclo negativo que no se alcanza desde s tampoco, porque ningún camino desde s pasa por él.

Idea clave. Bellman-Ford no necesita que los nodos tengan un orden. Se lo fabrica: la fila i contiene las respuestas para i aristas, y el lema dice en qué fila ya está la respuesta al problema original.

Por qué basta con una sola fila

La tabla usa memoria O(V²), y en la práctica nadie la guarda. La versión que se programa tiene un solo arreglo d, y relaja cada arista, una palabra que ya usamos en Dijkstra: si llegar por (u, v) mejora lo que se tiene, se actualiza.

algoritmo bellman_ford(grafo, s)
    d[v] ← ∞ y pred[v] ← NULO para todo nodo v
    d[s] ← 0
    para vuelta desde 1 hasta V - 1 hacer
        para cada arista (u, v) hacer
            si d[u] ≠ ∞ y d[u] + w(u, v) < d[v] entonces
                d[v] ← d[u] + w(u, v)
                pred[v] ← u
            fin si
        fin para
    fin para
    devolver d, pred
fin algoritmo

pred[v], el predecesor de v, es el nodo desde el que se relajó v por última vez. Aquí NULO significa que todavía no tiene ninguno. Lo vamos a necesitar pronto.

Ahora, algo que debería incomodarles. En la semana 7 vimos un programa que usaba una sola fila en la dirección equivocada, y sin avisar resolvía otro problema. ¿Pasa lo mismo aquí? Dentro de una vuelta, d[u] puede haber cambiado ya, así que d[v] se calcula con un valor de “esta” fila y no de la anterior. Miren este ejemplo:

Una sola fila, con las aristas en el orden de arriba:
vuelta  0    1    2    3
  1     0   -1    4    1
  2     0   -1    4    1

Una sola fila, con las aristas al revés:
vuelta  0    1    2    3
  1     0    1    4  inf
  2     0   -1    4    3
  3     0   -1    4    1
  4     0   -1    4    1

Con el primer orden, tras una sola vuelta d[3] ya vale 1, y sin embargo d_1(3) es ∞. d[v] después de la vuelta i no es d_i(v). Esa igualdad es falsa, y una demostración que la afirme está mal. Con el segundo orden, en cambio, cada vuelta avanza justo una arista, como la tabla. (La última vuelta de cada corrida es la que no cambia nada, y por eso el programa se detiene.) La afirmación correcta es más débil, y alcanza.

Teorema. Durante toda la corrida se cumple lo siguiente:

1. Si d[v] ≠ ∞, entonces d[v] es el largo de algún camino de s a v.
2. Los valores de d nunca suben.
3. Después de la vuelta i, d[v] ≤ d_i(v), para todo v.

En consecuencia, si no se alcanza ningún ciclo negativo, después de la vuelta V − 1 se cumple d[v] = δ(v) para todo v alcanzable, y d[v] = ∞ para los demás.

Técnica. Invariante, por inducción sobre las relajaciones para 1, y sobre las vueltas para 3. La conclusión aprieta a d[v] entre dos cosas iguales.

(1) y (2). Al empezar, d[s] = 0 es el largo del camino vacío. Una relajación solo cambia d[v] si lo baja, así que vale (2). Y el valor nuevo es d[u] + w(u, v). Como d[u] es el largo de un camino a u, ese camino más la arista (u, v) es un camino a v con ese largo.

(3). Para i = 0, d[s] = 0 = d_0(s), y los demás d_0(v) son ∞. Supongamos que vale para i − 1, y fijemos v. Por la recurrencia, d_i(v) es d_(i−1)(v) o d_(i−1)(u) + w(u, v) para alguna arista (u, v). En el primer caso, d[v] ≤ d_(i−1)(v) al terminar la vuelta i − 1, y (2) dice que no sube. En el segundo, si d_(i−1)(u) = ∞, entonces d_i(v) = ∞ y no hay nada que probar. Si no, al llegar a la arista (u, v) en la vuelta i, se cumple d[u] ≤ d_(i−1)(u). Esto es lo importante: que d[u] ya haya cambiado en esta vuelta solo ayuda, porque solo pudo haber bajado. Después de mirar esa arista se cumple d[v] ≤ d[u] + w(u, v) ≤ d_i(v), y (2) dice que eso se mantiene.

La conclusión. Sin ciclos negativos alcanzables, si v es alcanzable, (1) dice que d[v] es el largo de un camino, así que d[v] ≥ δ(v) una vez que es finito. Y (3) con i = V − 1, junto con el lema, dice que d[v] ≤ d_(V−1)(v) = δ(v). Las dos juntas dan la igualdad. Si v no es alcanzable, (1) dice que d[v] sigue en ∞. ∎

Compárenlo con la mochila. Allá, reusar un valor “de esta fila” era tomar el mismo objeto dos veces, y eso cambiaba el problema. Aquí, reusarlo es usar una arista más de lo que la vuelta permite, y tener más aristas nunca es ilegal: al final se quiere el mínimo sin límite de aristas. La fila i era un andamio, no una restricción del problema.

Error común. Decir que “después de la vuelta i, d[v] es el mejor camino con i aristas”. Con una sola fila eso es falso, como acabamos de ver. Es a lo sumo ese valor, y nunca menos que la distancia verdadera.

El orden de las aristas importa para el tiempo, no para el resultado

El programa bf_demo también cuenta las vueltas que hacen falta en una cadena, hasta que una vuelta no cambia nada. Agregar esa parada temprana no rompe nada: si una vuelta no cambia nada, las siguientes tampoco.

Cuántas vueltas hacen falta, hasta una que no cambia nada, en una cadena de n nodos:
    n  aristas en orden  aristas al revés
   10                 2                10
  100                 2               100
 1000                 2              1000

Con las aristas en el orden de la cadena, toda la información avanza en una sola vuelta. Al revés, avanza una arista por vuelta, y se necesitan todas. Así que Θ(VE) sigue siendo el peor caso con parada temprana, y este grafo lo alcanza.

En C++. Así queda en shortest_paths.h, con un valor INF que marca lo no alcanzado:

for (int round = 1; round <= n; round++) {
    result.rounds = round;
    int last_relaxed = -1;
    for (const Edge& e : edges) {
        if (result.dist[idx(e.from)] == INF) continue;  // no alcanzado: no hay nada que propagar
        if (result.dist[idx(e.from)] + e.weight < result.dist[idx(e.to)]) {
            result.dist[idx(e.to)] = result.dist[idx(e.from)] + e.weight;
            result.pred[idx(e.to)] = e.from;
            last_relaxed = e.to;
        }
    }
    if (after_round) after_round(round, result);
    if (last_relaxed == -1) break;
    if (round == n) result.witness = last_relaxed;
}

Las distancias van en long long: V − 1 aristas de peso grande no caben en int con facilidad. INF es el máximo de long long dividido entre 4, y se le puede sumar un peso sin overflow. Fíjense en que el ciclo llega hasta la vuelta n, una más de las que necesita el pseudocódigo. La razón viene ahora.

Detectar un ciclo negativo

Sin ciclos negativos, la vuelta V − 1 ya deja todo en su valor final. Entonces hagamos una vuelta más, la vuelta V, y miremos si algo cambia.

Teorema. La vuelta V relaja alguna arista si y solo si desde s se alcanza un ciclo negativo.

Técnica. Las dos direcciones, por contrarrecíproca. Una usa el teorema anterior, y la otra suma desigualdades alrededor del ciclo.

Si no se alcanza ningún ciclo negativo, la vuelta V no relaja nada. Después de la vuelta V − 1, cada nodo alcanzable tiene d[v] = δ(v), y los demás siguen en ∞. Una relajación dejaría un d[v] menor que δ(v). Por la parte (1) del teorema anterior, ese valor sería el largo de un camino de s a v, más corto que el mínimo. Es imposible.

Si la vuelta V no relaja nada, no se alcanza ningún ciclo negativo. En ese caso d no cambió en toda esa vuelta, así que al terminarla, toda arista (u, v) con d[u] ≠ ∞ cumple

d[v] ≤ d[u] + w(u, v)

Sea C un ciclo que se alcanza desde s, con nodos x₀, x₁, …, x_(L−1) y vuelta a x₀. Aquí se usa que se alcanza: cada nodo de C se alcanza desde s por un camino simple, que tiene a lo sumo V − 1 aristas. Entonces d_(V−1) es finito en ese nodo, y por la parte (3) del teorema, d también. Así que todos los d del ciclo son finitos, y la desigualdad se aplica a cada arista del ciclo. Sumándolas todas:

d[x₀] + d[x₁] + ⋯ + d[x_(L-1)]  ≤  d[x₀] + d[x₁] + ⋯ + d[x_(L-1)]  +  peso(C)

A cada lado está la misma suma, porque cada nodo del ciclo aparece una vez como inicio y una vez como fin de una arista. Como es finita, se cancela, y queda 0 ≤ peso(C). C no es negativo. ∎

El mismo algoritmo responde las dos preguntas: V − 1 vueltas para las distancias, y una más para saber si podemos creer en ellas.

Error común. Olvidar lo de “se alcanza”. Bellman-Ford desde s no ve un ciclo negativo al que s no llega. Si la pregunta es “¿hay un ciclo negativo en alguna parte del grafo?”, hay que hacer otra cosa, y esa es la que sigue.

Buscar en todo el grafo

Si la pregunta es sobre todo el grafo, se pueden poner todos los nodos en 0 al empezar, con pred en NULO. Es como agregar un origen nuevo con aristas de peso 0 hacia todos: ese origen alcanza cualquier ciclo, y las aristas nuevas no forman ningún ciclo, porque nadie entra a él.

Las demostraciones se repiten con un solo cambio: d_i(v) es ahora el menor largo de un camino que termina en v y empieza donde sea, con a lo sumo i aristas, y d_0(v) = 0 para todos. Un camino simple sigue teniendo a lo sumo V − 1 aristas. Así que la cuenta no cambia: V vueltas, y la última es la que decide. Háganlo como ejercicio, porque es fácil equivocarse en una vuelta.

Recuperar el ciclo

Saber que existe no siempre basta. Queremos ver el ciclo. Para eso guardamos pred. La tentación es esta: tomar un nodo que se relajó en la vuelta V y seguir sus predecesores hasta volver a él. No funciona, porque ese nodo puede no estar en el ciclo, sino colgando de él:

$ build/cycle_demo
Aristas: 0->1 (2), 1->2 (1), 2->3 (1), 3->1 (-3), 3->4 (1), 4->5 (1)
vuelta: d[0..5]
1: 0 1 3 4 5 6
2: 0 0 2 3 4 5
3: 0 -1 1 2 3 4
4: 0 -2 0 1 2 3
5: 0 -3 -1 0 1 2
6: 0 -4 -2 -1 0 1
nodo relajado al final de la vuelta 6: 5
pred: -1 3 1 2 3 4

Leer directamente desde ese nodo, hasta volver a él (se corta a los 12 pasos):
  5 4 3 2 1 3 2 1 3 2 1 3 ...

Retroceder 6 veces y recién entonces leer:
  5 4 3 2 1 3 -> 2
ciclo: 3 1 2 -> 3, peso -1

El ciclo 1 → 2 → 3 → 1 pesa −1, y cada vez que baja arrastra a la cola 4 y 5. El último nodo relajado es el 5, que no está en el ciclo. Leer desde él da un ciclo infinito en el programa: nunca se vuelve al 5.

La solución es retroceder V veces y recién entonces leer. Hay que demostrar que eso aterriza en un ciclo, y que ese ciclo es negativo. Lo que sigue vale en el grafo original, con s como origen.

Lema 1. En todo momento, para todo nodo v con pred[v] = u ≠ NULO, se cumple d[v] ≥ d[u] + w(u, v).

Demostración. Al relajar se cumple con igualdad. Después, d[u] solo puede bajar, y eso no rompe la desigualdad. Y si d[v] cambia, es por otra relajación que vuelve a poner igualdad con su nuevo pred. ∎

Lema 2. En todo momento, todo ciclo que formen los predecesores tiene peso negativo.

Demostración. Miremos el instante en que ese ciclo se cerró, es decir, la última relajación que cambió uno de sus predecesores: pred[v] ← u. Justo antes, los demás predecesores del ciclo ya estaban puestos y forman un camino v = x₀, x₁, …, x_k = u, con pred[x_(j+1)] = x_j. El lema 1 da d[x_(j+1)] ≥ d[x_j] + w(x_j, x_(j+1)) para cada j. Sumando, d[u] ≥ d[v]

• peso del camino. Como la relajación ocurrió, d[u] + w(u, v) < d[v] de manera estricta. Juntando las dos:

d[v] + peso del camino + w(u, v)  ≤  d[u] + w(u, v)  <  d[v]

Es decir, peso del camino + w(u, v) < 0, y esa suma es el peso del ciclo. Los pesos no cambian después, así que sigue siendo negativo mientras exista. ∎

Teorema. Si x se relajó en la vuelta V, entonces al terminar esa vuelta, retroceder V veces por pred desde x lleva a un nodo que está en un ciclo de predecesores, y ese ciclo es negativo.

Técnica. Por contradicción: si la cadena de predecesores termina en lugar de ciclarse, x no podría haberse relajado.

Seguir pred desde x da una secuencia en la que cada nodo define el siguiente. Hay dos posibilidades: llega a un nodo con pred = NULO, o sigue para siempre.

Supongamos que llega a uno. Mientras tanto no repitió nodos, porque si repite uno, ya no sale de ese ciclo. El nodo final tiene d finito, porque es predecesor de alguien, y no tiene predecesor, así que nunca se relajó: es s, y d[s] = 0. Hay entonces un camino simple P de s a x, y el lema 1, sumado a lo largo de P, da d[x] ≥ d[s] + peso(P) = peso(P).

Un camino simple tiene a lo sumo V − 1 aristas, así que peso(P) ≥ d_(V−1)(x). Pero la parte (3) del teorema de una fila dice que, al terminar la vuelta V − 1, d[x] ≤ d_(V−1)(x). Y en la vuelta V, x se relajó, así que ahora d[x] es estrictamente menor que eso. Tenemos d[x] < d_(V−1)(x) ≤ peso(P) ≤ d[x]. Es una contradicción.

Entonces la secuencia nunca termina. Como hay V nodos, entre los primeros V + 1 de la secuencia se repite uno, y a partir de ahí da vueltas en un ciclo. O sea que en a lo sumo V − 1 pasos ya está en el ciclo, y V pasos lo dejan dentro de él. Por el lema 2, ese ciclo es negativo. ∎

La demostración sirve con cualquier nodo relajado en la vuelta V, y también si se empezó con todos en 0: lo único que cambia es que un nodo sin pred es cualquiera que nunca se relajó, con d = 0, y que d_(V−1) es el de la versión con todos en 0.

algoritmo recuperar_ciclo(pred, x, V)
    para paso desde 1 hasta V hacer
        x ← pred[x]
    fin para
    ciclo ← lista con x
    y ← pred[x]
    mientras y ≠ x hacer
        agregar y a ciclo
        y ← pred[y]
    fin mientras
    invertir ciclo
    devolver ciclo
fin algoritmo

Se invierte porque se leyó de fin a inicio. En el ejemplo, esto da 3 → 1 → 2 → 3, que es el ciclo 1 → 2 → 3 → 1 escrito desde otro nodo.

Error común. Tomar directamente el nodo relajado en la vuelta V como parte del ciclo. En muchos grafos pequeños está en el ciclo, y el programa pasa las primeras pruebas. Luego aparece una cola, y queda un ciclo infinito.

Empezar en “infinito”: dos errores

El pseudocódigo pregunta d[u] ≠ ∞ antes de sumar. Esa línea es fácil de olvidar. En C++ no existe el ∞ de la recurrencia, y hay que elegir un número que haga de ∞. Las dos formas de hacerlo mal son distintas:

$ build/infinity_demo
infinity_demo.cpp:23:32: runtime error: signed integer overflow: 2147483647 + 5 cannot be represented in type 'int'
infinity_demo.cpp:24:47: runtime error: signed integer overflow: 2147483647 + 5 cannot be represented in type 'int'
Primer error: INT_MAX como infinito. Aristas 2->3 (5), 0->1 (7).
d con INT_MAX: 0 7 2147483647 -2147483644
correcto:      0 7 inf inf

Segundo error: infinito + un negativo. Aristas 0->1 (7), 2->3 (-5).
sin revisar:   0 7 1000000000000 999999999995
revisando:     0 7 inf inf

Y un falso ciclo: 0->1 (7), y el ciclo 2->3 (-5), 3->2 (1), que el origen no alcanza.
sin revisar: la vuelta 4 cambió algo: sí
revisando:   la vuelta 4 cambió algo: no

El primer error es usar INT_MAX y sumarle algo. INT_MAX + 5 es un overflow de int, y en C++ eso es comportamiento indefinido: el lenguaje no dice qué pasa. El programa está compilado con el sanitizer, que lo señala en las dos primeras líneas. Después de avisar sigue, y en esta corrida el resultado dio la vuelta y quedó en un número muy negativo: el nodo 3, que nadie alcanza, termina con distancia −2147483644. Otro compilador u otras opciones pueden hacer otra cosa. Por eso esa línea no es un resultado que se pueda esperar: es un bug.

El segundo error no tiene overflow. Con un “infinito” de 10¹² todo cabe, pero el nodo 2 no alcanzado tiene d = 10¹², y 10¹² − 5 sí es menor que 10¹². El programa “mejora” al nodo 3, con una distancia sin sentido, que no corresponde a ningún camino desde s. Se rompe la parte (1) del teorema, y con ella todo lo que sigue. La última parte de la salida lo muestra: un ciclo negativo que el origen no alcanza se detecta como si lo alcanzara. El teorema de detección decía otra cosa, y el programa dejó de cumplirlo.

Lo correcto es un valor sin overflow y la pregunta explícita antes de sumar: de un nodo no alcanzado no se propaga nada.

Error común. Partir de un valor “infinito” y sumarle a un nodo no alcanzado. Con pesos no negativos, como en Dijkstra, el segundo error no se notaba: ∞ + w nunca es menor que ∞. Con negativos, sí.

Dónde se usa: el vector de distancias

Mirar solo la recurrencia sugiere algo curioso. Para actualizar d[v], un nodo necesita los valores de sus vecinos y el peso de la arista, nada más. No hace falta conocer el grafo completo. Dijkstra, en cambio, necesita un heap global.

Eso es lo que usan los routers con un vector de distancias (distance vector). Cada router guarda su estimación de la distancia a cada destino, y se la envía a sus vecinos. Al recibir el vector de un vecino, relaja: si llegar a un destino por ese vecino es mejor, actualiza su estimación y anota al vecino como siguiente salto. RIP es un protocolo de ruteo que usa esta idea. Los detalles, incluido lo que pasa cuando un enlace se cae, los verán en Redes de Computadoras. Lo que importa aquí es que la recurrencia de este capítulo es la misma, y se puede calcular sin que nadie vea el mapa entero.

Todos los pares: V veces un algoritmo de un origen

Sesión B. Ahora queremos δ(i, j) para todo par de nodos: una tabla de distancias entre todas las ciudades, por ejemplo. La primera idea es correr V veces lo que ya tenemos, una vez por origen.

En un grafo denso, con E cercano a V², V × Bellman-Ford es O(V⁴). Si hay pesos negativos, parece que ese es el precio. La idea de esta sesión baja eso a Θ(V³), con un programa de pocas líneas.

Floyd-Warshall: el subproblema con nodos intermedios permitidos

En Bellman-Ford el parámetro era cuántas aristas usar. Probar lo mismo aquí da V capas de tablas de V². Floyd-Warshall tiene otra idea: un parámetro que limite por dónde puede pasar el camino.

Los nodos se numeran de 0 a V − 1. Los nodos intermedios de un camino son todos, menos el primero y el último.

Subproblema. Sin ciclos negativos, d_k(i, j) es el menor largo de un camino de i a j cuyos nodos intermedios están todos entre los k primeros nodos, 0, 1, …, k − 1. Es ∞ si no hay ninguno. La respuesta es d_V(i, j), donde ya se permite todo.

Ojo con los índices. d_0 no permite ningún intermedio, y d_(k+1) es el que agrega al nodo k.

d_0(i, i) = mín(0, w(i, i)), donde w(i, i) es el peso del lazo, si lo hay
d_0(i, j) = w(i, j), si hay arista de i a j, y ∞ si no                         para i ≠ j
d_(k+1)(i, j) = mín( d_k(i, j),  d_k(i, k) + d_k(k, j) )

Teorema. Si no hay ciclos negativos, la recurrencia calcula d_k(i, j) para todo k de 0 a V, y d_V(i, j) = δ(i, j).

Técnica. Inducción sobre k, con las dos direcciones. La decisión es sobre un nodo: si el camino pasa o no por k.

Caso base. Sin intermedios, un camino de i a j es la arista directa o, si i = j, el vacío. Un ciclo que vuelva a i sin intermedios sería un lazo, una arista de i a i. Si hay uno de peso no negativo no mejora al camino vacío.

Paso. Supongamos que la recurrencia es correcta para k.

(≥) Las dos opciones son largos de caminos permitidos en d_(k+1). La primera, porque tener menos intermedios permitidos sigue siendo permitido. La segunda, porque pegar un camino de i a k y uno de k a j, ambos con intermedios entre los k primeros, da un camino de i a j cuyos intermedios son esos y el propio k.

(≤) Tomemos un camino P de i a j con intermedios entre los k + 1 primeros. Como no hay ciclos negativos, quitarle ciclos no lo alarga y no agrega intermedios, igual que en el lema de V − 1 aristas, así que podemos suponer que P es simple. Si k no es intermedio de P, P ya cuenta en d_k(i, j). Si lo es, aparece una sola vez, porque P es simple. Parte a P en dos: de i a k y de k a j. Ninguno tiene a k como intermedio, así que los dos cuentan en d_k, y el largo de P es al menos d_k(i, k) + d_k(k, j). ∎

Esa frase en negrita es el centro. Es lo que hace falta para partir en dos, y es la razón por la que se pide que no haya ciclos negativos: con uno, el mejor camino no tendría por qué ser simple. Veremos abajo qué se hace entonces.

El programa, y el costo

algoritmo floyd_warshall(grafo)
    d[i][j] ← d_0(i, j) para todo par (i, j)
    para k desde 0 hasta V - 1 hacer
        si d[k][k] < 0 entonces
            devolver "hay un ciclo negativo"
        fin si
        para i desde 0 hasta V - 1 hacer
            para j desde 0 hasta V - 1 hacer
                si d[i][k] ≠ ∞ y d[k][j] ≠ ∞ y d[i][k] + d[k][j] < d[i][j] entonces
                    d[i][j] ← d[i][k] + d[k][j]
                fin si
            fin para
        fin para
    fin para
    devolver d
fin algoritmo

Tres ciclos anidados de V vueltas: Θ(V³), sin importar E. Es la misma cuenta para un grafo disperso o uno denso.

Sigamos un ejemplo con un peso negativo. Cada matriz tiene el origen en la fila y el destino en la columna:

$ build/fw_demo
Aristas: 0->1 (4), 0->3 (9), 1->2 (-2), 2->3 (1), 3->1 (5), 2->0 (3)
d_0, sin nodos intermedios:
    0    4  inf    9
  inf    0   -2  inf
    3  inf    0    1
  inf    5  inf    0
d_1, tras permitir el nodo 0:
    0    4  inf    9
  inf    0   -2  inf
    3    7    0    1
  inf    5  inf    0
d_2, tras permitir el nodo 1:
    0    4    2    9
  inf    0   -2  inf
    3    7    0    1
  inf    5    3    0
d_3, tras permitir el nodo 2:
    0    4    2    3
    1    0   -2   -1
    3    7    0    1
    6    5    3    0
d_4, tras permitir el nodo 3:
    0    4    2    3
    1    0   -2   -1
    3    6    0    1
    6    5    3    0

Sigan las celdas que cambian. En d_1 aparece el camino 2 → 0 → 1 = 3 + 4 = 7. En d_2, 0 → 1 → 2 = 4 − 2 = 2, y también 3 → 1 → 2 = 5 − 2 = 3. En d_3, permitir el nodo 2 mejora 0 → 3, de 9 a 2 + 1 = 3, y conecta a 1 con 0 y con 3. En d_4, la casilla de 2 a 1 baja de 7 a 1 + 5 = 6. Cada matriz responde una pregunta distinta, y solo la última es la de distancias.

Por qué se puede actualizar en el mismo lugar

La recurrencia lee d_k, pero el programa guarda una sola matriz: cuando lee d[i][k], esa casilla pudo haberse actualizado ya en esta vuelta. ¿Lee de la capa k o de la k + 1? Igual que con Bellman-Ford, hay que demostrar que no importa, y aquí el motivo es aún más limpio.

Lema. Sin ciclos negativos, durante la vuelta de k no cambia ninguna casilla de la fila k ni de la columna k.

Demostración. d[k][k] = 0: el camino vacío mide 0, y ningún camino cerrado mide menos porque no hay ciclos negativos. Entonces, para la fila k, la recurrencia mira d[k][k] + d[k][j] = d[k][j], que no es menor que d[k][j]. Para la columna k, mira d[i][k] + d[k][k] = d[i][k], y tampoco. ∎

Entonces los valores d[i][k] y d[k][j] que el programa lee en la vuelta de k son exactamente d_k(i, k) y d_k(k, j), sin importar en qué orden se recorren i y j. Por eso, al terminar la vuelta, cada casilla vale d_(k+1). El lema depende de que no haya ciclos negativos, y por eso el programa revisa d[k][k] antes de cada vuelta.

Por qué el ciclo de k va por fuera

Esta es la pregunta de repaso 4, y el error más característico de esta semana. Si se escriben los tres ciclos con k adentro,

para i ...
    para j ...
        para k ...

el programa compila, termina y da números creíbles. Pero está mal, y no se puede ver mirando la salida:

El ciclo de k por dentro. Aristas: 1->3 (5), 3->2 (5), 2->0 (9), 1->0 (30)
k por fuera:
    0  inf  inf  inf
   19    0   10    5
    9  inf    0  inf
   14  inf    5    0
k por dentro:
    0  inf  inf  inf
   30    0   10    5
    9  inf    0  inf
   14  inf    5    0

Grafos al azar en los que k por dentro se equivoca (20000 por fila, 35 % de aristas, pesos 1 a 9):
  n = 3: 0
  n = 4: 1534
  n = 5: 5272
  n = 8: 17060

Miremos la casilla de 1 a 0. El camino bueno es 1 → 3 → 2 → 0, de largo 19, y el programa malo da 30, la arista directa. ¿Por qué? Con k adentro, cuando se calcula d[1][0] se prueba una sola vez cada intermedio:

Después, d[1][2] y d[3][0] se vuelven finitos, pero nadie vuelve a d[1][0]. Con k por fuera, las casillas que hacen falta se llenan en la vuelta anterior. La recurrencia dice cuál es el orden: toda la capa k antes de la capa k + 1, y la capa está marcada por k. Es el error de la semana 9, leer una casilla antes de llenarla, escondido en el orden de tres ciclos.

La última parte de la salida es la razón para desconfiar. Con 3 nodos, nuestra prueba no encontró ningún error en 20000 grafos, y con 4 lo encontró en menos de 1 de cada 10. Con pruebas chiquitas, un programa así pasa.

Error común. Poner el ciclo de k por dentro. Compila, da respuestas plausibles y está mal. Una prueba con grafos muy pequeños puede no notarlo. Lo que hay que recordar es la frase: k es el parámetro del subproblema, así que va por fuera.

Detectar un ciclo negativo con la diagonal

La matriz también responde esta pregunta, sin nada más.

Teorema. El grafo tiene un ciclo negativo si y solo si en algún momento el programa ve d[k][k] < 0 al empezar la vuelta de k.

Técnica. Las dos direcciones. Para una, cada casilla es el largo de algún camino. Para la otra, una cota que vale siempre, aunque haya ciclos negativos, aplicada al nodo mayor del ciclo.

Si ve d[k][k] < 0, hay un ciclo negativo. Por inducción, cada casilla finita es el largo de algún camino: al empezar, lo son, y un valor nuevo pega dos caminos. Entonces d[k][k] < 0 es un camino cerrado de largo negativo. Un camino cerrado se descompone en ciclos simples, porque cada vez que repite un nodo se le puede separar un ciclo. Sus pesos suman un negativo, así que alguno es negativo.

Si hay un ciclo negativo, el programa lo ve. Primero, una cota que no necesita suponer nada sobre ciclos:

Tras procesar los nodos 0, …, k − 1, para todo par (i, j), incluso con i = j, d[i][j] es menor o igual que el largo de cualquier camino de i a j cuyos nodos intermedios sean distintos entre sí y estén entre los k primeros.

Para k = 0 es la definición de d_0, que en la diagonal toma el mínimo entre 0 y el peso del lazo, si lo hay. Para pasar de k a k + 1, si el camino no usa a k como intermedio, vale la hipótesis, y los valores solo bajan. Si lo usa, lo usa una vez, y se parte en dos caminos a los que se aplica la hipótesis. Al mirar el par (i, j) en la vuelta de k, las casillas d[i][k] y d[k][j] no valen más que al empezar la vuelta, y por eso d[i][j] queda por debajo de la suma.

Ahora, sea C un ciclo negativo simple, que existe porque un ciclo negativo se descompone en simples, y sea h su nodo de mayor número. Si el programa se detuvo antes de la vuelta de h, ya lo vio. Si no, llega a la vuelta de h. Allí, C es un camino de h a h, con intermedios distintos, todos menores que h. La cota dice d[h][h] ≤ peso(C) < 0. Lo ve. ∎

Aquí está la pregunta de repaso 5, y el programa en acción:

La diagonal. Aristas: 0->1 (1), 1->2 (-3), 2->1 (2), 2->3 (1)
  antes de k = 0, diagonal: 0 0 0 0
  antes de k = 1, diagonal: 0 0 0 0
  antes de k = 2, diagonal: 0 0 -1 0
  d[2][2] < 0: hay un ciclo negativo, y se detiene

El ciclo es 1 → 2 → 1, de peso −1. Su nodo mayor es el 2, y tras permitir al nodo 1, la casilla d[2][2] ya lo muestra.

Por qué se detiene: los valores pueden crecer muchísimo

¿Y si no se detiene, y se revisa la diagonal al final? En principio funciona, pero con un ciclo negativo los valores pueden crecer de forma exponencial, porque cada vuelta suma dos valores ya negativos:

Sin detenerse: 40 nodos, todas las aristas de peso -1.
  tras la vuelta k =  0: menor valor -2.000e+00
  tras la vuelta k =  1: menor valor -1.200e+01
  tras la vuelta k =  2: menor valor -7.200e+01
  tras la vuelta k =  3: menor valor -4.320e+02
  tras la vuelta k =  4: menor valor -2.592e+03
  tras la vuelta k =  5: menor valor -1.555e+04
  tras la vuelta k = 24: menor valor -9.477e+18  (ya no cabe en long long)

En este grafo se multiplica por 6 en cada vuelta. Con 40 nodos, sin ningún peso grande, en la vuelta 24 ya no cabe en long long. El programa de arriba lo mide con double, justamente para poder mostrarlo sin overflow. En un programa con enteros, eso sería un overflow, y la diagonal dejaría de significar algo.

La revisión de d[k][k] al empezar cada vuelta lo evita. Mientras el programa siga, no hay ningún ciclo negativo entre los k primeros nodos, porque si lo hubiera, ya se habría detenido en su nodo mayor. Entonces d_k tiene su significado de siempre, y sus casillas son largos de caminos con intermedios distintos: a lo sumo V aristas por el mayor peso, en valor absoluto. Y como d[k][k] ≥ 0, de hecho es 0, y el lema de la fila y la columna k sigue valiendo en esa vuelta.

En C++. En all_pairs.h, la vuelta y la guardia:

inline void floyd_round(Matrix& d, int n, int k) {
    for (int i = 0; i < n; i++) {
        if (d[idx(i)][idx(k)] == INF) continue;
        for (int j = 0; j < n; j++) {
            if (d[idx(k)][idx(j)] == INF) continue;
            const long long through = d[idx(i)][idx(k)] + d[idx(k)][idx(j)];
            if (through < d[idx(i)][idx(j)]) d[idx(i)][idx(j)] = through;
        }
    }
}

Los dos continue son el mismo cuidado de Bellman-Ford: nada de ∞ + un negativo.

Qué conviene, según las mediciones

make bench compara las tres formas, con -O2, midiendo una sola corrida con std::chrono::steady_clock. Estos tiempos dependen de la máquina y de lo que estuviera corriendo. Los pesos del grafo disperso y del completo son enteros de 1 a 1000. Los tres métodos dan matrices iguales en todas las filas.

$ build/all_pairs_bench
grafo                V        E V Dijkstra      V B-F        F-W
disperso, E=4V     100      399      0.002      0.001      0.001  iguales
disperso, E=4V     200      797      0.009      0.003      0.006  iguales
disperso, E=4V     400     1595      0.030      0.010      0.042  iguales
completo           100     9900      0.010      0.004      0.001  iguales
completo           200    39800      0.066      0.040      0.009  iguales
completo           400   159600      0.447      0.336      0.061  iguales
cadena al revés   200      199      0.001      0.004      0.003  iguales
cadena al revés   400      399      0.005      0.025      0.025  iguales
cadena al revés   800      799      0.018      0.193      0.194  iguales

Con todas las aristas, Floyd-Warshall gana con claridad. En el grafo disperso de 400 nodos es el más lento. V × Bellman-Ford le gana a V × Dijkstra en estos grafos al azar, cosa que la cota no sugiere. bellman_ford tiene parada temprana, y contando las vueltas aparte, en grafos generados igual, con 400 nodos fueron unas 8 en promedio, muy lejos de las 399. Por eso esas filas no muestran su peor caso. La cadena con las aristas al revés sí lo acerca: duplicar V multiplicó el tiempo de V × Bellman-Ford por 6 a 8. No saquen de esto más conclusiones de las que permite una sola corrida en tan pocos tamaños: la cota O(V²E) es un peor caso, y estos grafos no son una red real.

Cuadro de decisión

Ahora tenemos cinco herramientas. Cada una tiene una condición, y fuera de esa condición, o se equivoca, o se paga de más.

La tabla responde la pregunta de repaso 3. Sin pesos: BFS. DAG con pesos negativos: orden topológico. Red de carreteras: Dijkstra, con distancias no negativas. Todos los pares en un grafo pequeño y denso: Floyd-Warshall.

El programa decision_demo lo comprueba en un ejemplo de cada fila, frente a Bellman-Ford. Esta es la parte más interesante, con el mismo grafo que rompió a Dijkstra, que resulta ser un DAG:

$ build/decision_demo
DAG con un peso negativo: 0->1 (1), 0->2 (4), 2->1 (-5), 1->3 (2)
  orden topológico: 0 -1 4 1
  Bellman-Ford:     0 -1 4 1
  Dijkstra:         0 -1 4 3

Con pesos negativos, un DAG no necesita Bellman-Ford. Procesar los nodos en orden topológico es la programación dinámica de la semana 9: cuando se llega a un nodo, todos sus predecesores ya están listos. Es una de las razones para mirar la forma del grafo antes de elegir.

En la práctica. Con ciclos negativos, ninguno de estos algoritmos da “el camino simple más corto”. Ese problema es otro: es el camino más largo de la semana pasada, con los pesos cambiados de signo, que retomaremos en la semana 14. Lo que sí hacen Bellman-Ford y Floyd-Warshall es avisar.

✂️ Máscaras de bits y Held-Karp

Esta sección se puede cortar, o dejar como lectura si falta tiempo. Quien la salte tendrá que leerla antes del laboratorio de la semana 15, que usa Held-Karp. Muestra un parámetro de otra clase: un subconjunto.

Un subconjunto como un entero

Un subconjunto de {0, …, n − 1} se puede guardar en un entero: el bit v vale 1 si v está en el subconjunto. Por ejemplo, {0, 2, 3} con n = 4 es 1101 en binario, o sea 13. A ese entero se le llama máscara de bits (bitmask). Las operaciones cuestan O(1):

Hay una propiedad muy útil para una tabla: quitar un bit da un entero menor. Así que recorrer las máscaras de 0 a 2ⁿ − 1 visita cada subconjunto después de todos sus subconjuntos propios. El orden de cálculo sale gratis.

El problema del viajero

Hay n ciudades numeradas desde 0 y un costo c(u, v) para ir de u a v. El problema del viajero (traveling salesperson problem) pide el recorrido más barato que sale de 0, pasa una sola vez por cada ciudad y vuelve a 0. Probar todos los órdenes cuesta (n − 1)! recorridos. Con 20 ciudades eso es 19! ≈ 1,2 × 10¹⁷. Más adelante en el curso veremos por qué no se espera un algoritmo polinomial para él. Aun así, se puede hacer mucho mejor que n!.

¿Qué necesita saber el futuro de un recorrido a medias? No el orden en que se visitaron las ciudades. Solo cuáles se visitaron y dónde estamos. Es el grafo de estados de la semana 5, con un subconjunto en vez de un contador.

Subproblema. Para un subconjunto S de nodos distintos de 0, no vacío, y un v en S, C(S, v) es el costo mínimo de un camino que sale de 0, visita exactamente los nodos de S, cada uno una vez, y termina en v. La respuesta es el mínimo sobre v de C(todos, v) + c(v, 0).

C({v}, v) = c(0, v)
C(S, v) = mín sobre u en S, u ≠ v, de C(S - {v}, u) + c(u, v)         si S tiene dos nodos o más

Teorema. La recurrencia calcula C(S, v), y da el recorrido mínimo.

Técnica. Inducción sobre el tamaño de S, con las dos direcciones, mirando el penúltimo nodo.

Caso base. Con S = {v}, el único camino es la arista de 0 a v.

Paso. Sea |S| ≥ 2. (≤) Un camino válido, quitándole la última arista (u, v), visita exactamente S − {v} y termina en un u ≠ v de ese conjunto. Ese tramo cuesta al menos C(S − {v}, u), así que el camino cuesta al menos la opción de u. (≥) Cada opción es un camino válido: como v no está en S − {v}, agregar (u, v) a un camino que visitó S − {v} no repite ningún nodo. La respuesta sale igual: un recorrido es un camino que visita a todos y termina en algún v, más la vuelta c(v, 0). ∎

Esa última frase de (≥) es el motivo de tener S en el subproblema: sin él, no se podría saber si v ya se visitó.

Con n nodos hay 2ⁿ⁻¹ subconjuntos de los otros n − 1, y por cada uno a lo sumo n − 1 finales: memoria O(n·2ⁿ). Cada casilla mira a lo sumo n − 1 opciones: tiempo O(n²·2ⁿ). Es el algoritmo de Held-Karp. Sigue siendo exponencial, pero bastante menos.

algoritmo held_karp(c, n)
    para v desde 1 hasta n - 1 hacer
        C[{v}][v] ← c(0, v)
    fin para
    para cada subconjunto S de {1, ..., n - 1} con dos nodos o más, en orden creciente de máscara hacer
        para cada v en S hacer
            C[S][v] ← mín sobre u en S - {v} de C[S - {v}][u] + c(u, v)
        fin para
    fin para
    devolver mín sobre v de C[{1, ..., n - 1}][v] + c(v, 0)
fin algoritmo

En held_karp.h los bits representan a los nodos 1, …, n − 1: el bit v − 1 es el nodo v. Aquí está la tabla completa con 4 ciudades:

$ build/held_karp_demo
Costos:
   0  10  15  20
   5   0   9  10
   6  13   0  12
   8   8   9   0

máscara  subconjunto  termina en 1  en 2  en 3
      1          {1}            10     -     -
      2          {2}             -    15     -
      3        {1,2}            28    19     -
      4          {3}             -     -    20
      5        {1,3}            28     -    20
      6        {2,3}             -    29    27
      7      {1,2,3}            35    29    31

mejor recorrido, Held-Karp: 35
mejor recorrido, permutaciones: 35

Sigamos dos casillas. C({1, 2}, 2) = C({1}, 1) + c(1, 2) = 10 + 9 = 19. C({1, 2, 3}, 2) mira dos opciones, C({1, 3},

1. • c(1, 2) = 28 + 9 = 37 y C({1, 3}, 3) + c(3, 2) = 20 + 9 = 29, y se queda con 29. La respuesta es el mínimo entre 35 + 5, 29 + 6 y 31 + 8, que es 35: el recorrido 0 → 1 → 3 → 2 → 0.

El resto de la salida cuenta el trabajo. Una transición es probar un u para una casilla, y son exactamente (n − 1)(n − 2)·2ⁿ⁻³, una fórmula que la prueba comprueba:

 n  recorridos (n-1)!  transiciones  n^2 2^n     mismo costo
 4                  6            12      256             sí
 5                 24            48      800             sí
 6                120           160     2304             sí
 7                720           480     6272             sí
 8               5040          1344    16384             sí
 9              40320          3584    41472             sí
10             362880          9216   102400             sí
11            3628800         23040   247808             sí

Desde 7 ciudades hay menos transiciones que recorridos, y un recorrido además cuesta n sumas. En tiempo, con make bench: una corrida con -O2, un solo grafo al azar por fila, y tiempos que dependen de la máquina. La memoria no está medida: es la cuenta n·2ⁿ⁻¹ × 8 bytes.

$ build/held_karp_bench
  n  permutaciones      Held-Karp memoria de tabla
  8        0.000 s        0.000 s           0.0 MB  iguales
  9        0.000 s        0.000 s           0.0 MB  iguales
 10        0.004 s        0.000 s           0.0 MB  iguales
 11        0.036 s        0.000 s           0.1 MB  iguales
 12        0.358 s        0.001 s           0.2 MB  iguales
 14              -        0.003 s           0.9 MB
 16              -        0.013 s           4.2 MB
 18              -        0.067 s          18.9 MB
 20              -        0.294 s          83.9 MB

Las permutaciones van unas 10 veces más lento por cada ciudad, y ahí se dejó de correr. Held-Karp resuelve 20 ciudades en menos de un segundo. Una extrapolación, no una medición: si cada recorrido cuesta lo que costó con n = 12, 19! recorridos tomarían más de 30 años. El límite de Held-Karp suele ser otro, la memoria: cada dos ciudades más, la tabla se multiplica por algo más de 4.

Del producto a la suma: los tipos de cambio

Esto ya es para el laboratorio. Una tabla de tipos de cambio dice cuántas unidades de una moneda se reciben por una de otra. Un arbitraje es una cadena de cambios que vuelve a la moneda inicial con más dinero: un ciclo cuyas tasas multiplicadas dan más de 1. ¿Cómo se busca eso con lo que tenemos, si todo lo nuestro suma?

Afirmación. Sea un ciclo con tasas r₁, …, r_L, todas positivas. Con el peso w = −ln(r) en cada arista, el producto de las tasas es mayor que 1 si y solo si la suma de los pesos es negativa.

Demostración. El logaritmo es creciente, y ln(r₁ ⋯ r_L) = ln r₁ + ⋯ + ln r_L. Entonces r₁ ⋯ r_L > 1 ⟺ ln r₁ + ⋯ + ln r_L > 0 ⟺ −ln r₁ − ⋯ − ln r_L < 0. ∎

El logaritmo natural ln es el que trae C++ como std::log. La base no cambia nada, porque solo cambia la escala. Un arbitraje es un ciclo negativo. Y si se quiere en cualquier parte de la tabla, no solo desde una moneda, se usa la versión con todos en 0.

Ahora, el problema. double no guarda números reales, los redondea. Miren lo que pasa con productos que valen 1 exacto:

$ build/rates_demo
Producto real y suma de -ln(tasa) en double:
tasas                    producto         suma
2 y 0,5                  1 exacto    0.000e+00
1,25 y 0,8               1 exacto   -5.551e-17
1,6 y 0,625              1 exacto   -5.551e-17
2,5 y 0,4                1 exacto   -1.110e-16
1,6, 1,25 y 0,5          1 exacto   -1.110e-16
1,25, 1,25 y 0,64        1 exacto   -5.551e-17
1,1, 1,1 y 1/1,21 (*)   menor a 1    4.996e-16
(*) 0,826446280991735 es 1/1,21 truncado a 15 decimales.

Con 1,25 y 0,8 no hay arbitraje, pero la suma sale negativa. 2 y 0,5 tienen suerte: sus logaritmos son exactamente opuestos. Y con una tasa de vuelta truncada, la suma sale positiva, porque el producto real es un poco menor que 1. No se puede saber el signo del error sin hacer la cuenta, que es justo lo que no queremos hacer para cada ciclo. Estos valores salen de nuestro entorno: otra biblioteca podría redondear distinto en la última cifra.

Una tabla que parece sana da entonces un arbitraje falso, si se compara con 0 tal cual. Hace falta una tolerancia. Y una tolerancia necesita algo que la justifique: si se aceptara cualquier ganancia, por pequeña que sea, nunca se podría distinguir de un redondeo. Por eso el problema promete un margen.

La promesa. El producto de todo ciclo simple es a lo sumo 1, o es por lo menos 1 + g, con g > 0 conocido. Basta con los simples: un camino cerrado se descompone en ciclos simples.

La tolerancia. Sumar un mismo ε > 0 al peso de cada arista: w = −ln(r) + ε. Es un grafo con otros pesos, así que todas las demostraciones de este capítulo siguen valiendo tal cual. Lo que cambia es qué ciclos son negativos: uno de L aristas ahora pesa −ln(P) + L·ε, donde P es su producto.

Teorema. En aritmética exacta, si L·ε < ln(1 + g) para todo largo de ciclo posible, entonces un ciclo es negativo con los pesos nuevos si y solo si tiene arbitraje.

Demostración. Si P ≥ 1 + g, el peso nuevo es −ln(P) + L·ε ≤ −ln(1 + g) + L·ε < 0. Si P ≤ 1, entonces −ln(P) ≥ 0, y el peso nuevo es al menos L·ε > 0. Por la promesa, no hay otros casos. ∎

Además, un ciclo recuperado es de verdad: si su peso nuevo es negativo, entonces P > 1, y por la promesa, P ≥ 1 + g.

Eso dice que ε no puede ser grande. ¿Cuán pequeño puede ser? Lo bastante para que aguante el redondeo. Hagamos la cuenta con los límites del laboratorio: g = 10⁻⁶, a lo sumo 400 monedas, tasas entre 10⁻⁶ y 10⁶.

• Cota superior. ln(1 + 10⁻⁶) es casi 10⁻⁶, y L ≤ 400, así que ε < 2,5 × 10⁻⁹.
• Cota inferior. Sin arbitraje, cada d es el largo de un camino, que no baja del de un camino simple. Cada peso mide a lo sumo ln(10⁶) ≈ 13,8, y hay a lo sumo 399 aristas, así que |d| < 5526. Una suma en double se equivoca a lo sumo en 2⁻⁵³ de su resultado, unos 6 × 10⁻¹³. El error de std::log y el de leer la tasa son mucho menores, porque el peso es mucho menor que d. Si ε supera con holgura esos 6 × 10⁻¹³, cada arista aporta más de lo que puede quitar el redondeo.

Cualquier ε entre esas cotas sirve, y ε = 10⁻⁹ queda a 2,5 veces de la de arriba y a más de mil de la de abajo. Con él, un ciclo de 400 aristas con la ganancia mínima sigue negativo: 400 × 10⁻⁹ = 4 × 10⁻⁷, menos que 10⁻⁶. Y sobra para el redondeo: harían falta valores de |d| de más de diez millones para consumir la diferencia que queda. Este programa lo mide, con la tabla exacta de cinco monedas de rates.h y un ciclo de 400 con producto exactamente 1,000001:

Con pesos -ln(tasa) + epsilon, ¿la vuelta n relaja algo?
   epsilon tabla exacta (sin arbitraje)       ciclo de 400, 1,000001
     0e+00        sí (falso arbitraje)                          sí
     1e-15                           no                          sí
     1e-12                           no                          sí
     1e-09                           no                          sí
     1e-08                           no            no (se le escapa)
     1e-06                           no            no (se le escapa)

Con ε = 0, la tabla exacta declara un arbitraje que no existe. Con ε = 10⁻⁸ se escapa uno verdadero: 400 × 10⁻⁸ = 4 × 10⁻⁶, más que la ganancia. Que 10⁻¹⁵ pase en esta tabla chica no es una garantía: la cuenta de arriba solo la da por encima de 6 × 10⁻¹³.

Error común. Elegir ε “porque se ve pequeño”. Una tolerancia por arista se acumula por todo el ciclo. Lo que tiene que quedar por debajo de la ganancia prometida es L·ε, no ε.

Error común. Comparar con ε, pero guardar el valor sin ε: si d[u] + w < d[v] - ε entonces d[v] ← d[u] + w. Se parece mucho, pero ya no es Bellman-Ford sobre otros pesos, y la demostración de la vuelta hacia atrás no se traslada: d puede quedarse hasta ε por encima de lo que dice la tabla en cada arista. Sumar ε al peso evita tener que pensarlo.

Dos advertencias más. El teorema para recuperar el ciclo está demostrado en aritmética exacta. Con double, una mejora ínfima, de puro redondeo, podría llegar a un nodo sin predecesor al retroceder, y por eso conviene revisar que no aparezca NULO. Y con arbitraje, nada garantiza esa cota de 5526. En una corrida con todas las tablas del laboratorio, el menor d con que se llegó a la última vuelta fue −704,591, dentro de ella. Pero en un ciclo de 400 monedas con todas las tasas en 10⁶ y las aristas en el orden del ciclo, cada vuelta baja d en unos 5526, y con una sola fila d llega a unos −2,2 × 10⁶.

Floyd-Warshall también resuelve este problema, con la misma tolerancia y los mismos 400 nodos, pero recuperar el ciclo pide más trabajo.

El código de la semana

$ make test
mkdir -p build
g++ -Wall -Wextra -std=c++20 -g -fsanitize=address,undefined test_examples.cpp -o build/test_examples
build/test_examples
Todas las pruebas de la semana 10 pasaron.

Estas pruebas comprueban, entre otras cosas, lo siguiente:

• La tabla coincide con un árbitro que prueba todos los caminos de hasta i aristas, con cualquier signo y con ciclos negativos, en 1500 grafos.
• Tras cada vuelta, d[v] ≤ d_i(v), los valores solo bajan y, sin ciclos negativos, nunca quedan bajo la verdad. La fila V − 1 no cambia en la fila V, y una sola fila llega a lo mismo.
• La vuelta V relaja algo si y solo si un árbitro, que recorre todos los ciclos simples, encuentra uno negativo alcanzable, en 4000 grafos. Cada ciclo recuperado es simple, usa aristas del grafo y pesa menos que cero. En algunos, el nodo relajado no estaba en el ciclo.
• Un ciclo que el origen no alcanza no se detecta, y sus nodos quedan en ∞.
• BFS, el orden topológico y Dijkstra coinciden con Bellman-Ford dentro de su condición.
• Floyd-Warshall coincide con V × Bellman-Ford en 3000 grafos, y detecta los ciclos negativos que encuentra el árbitro, y solo esos. La fila y la columna k no cambian en la vuelta k, y dos capas dan lo mismo que una. Con k por dentro, el ejemplo da 30 en lugar de 19.
• Las sumas de rates_demo, y la ventana de ε: ε = 0 falla en la tabla exacta, ε = 10⁻⁹ acierta en los dos casos, y ε = 10⁻⁸ falla en el ciclo.
• Held-Karp coincide con las permutaciones en 400 instancias, con aristas faltantes y con n = 1 y n = 2. La cantidad de transiciones es exactamente la de la fórmula.

El laboratorio

Entren a lab/arbitrage. El enunciado define con precisión qué es una cadena y qué se promete de la tabla, y qué imprimir en cada caso. El esqueleto siempre responde NO, así que falla en todos los casos con arbitraje. Puede haber varias cadenas válidas, y por eso hay un verificador: recalcula el producto con las tasas de la entrada, en long double, y no confía en la cuenta de ustedes.

Piensen antes en estas preguntas. ¿Hace falta correr desde una sola moneda? ¿Cuántas vueltas son? ¿Qué hay que imprimir si el nodo relajado no está en el ciclo? ¿Qué número van a usar como ε, y por qué? Y prueben con su propia solución lenta en tablas chicas: es la forma de saber que la rápida acierta.

Resumen

• Por qué otra idea: con aristas negativas, un nodo fijado todavía puede mejorar, y Dijkstra pierde su demostración. Con un ciclo negativo alcanzable, δ(v) puede no existir.
• Subproblema de Bellman-Ford: d_i(v) es el menor largo con a lo sumo i aristas. d_i(v) = mín(d_(i−1)(v), mín de d_(i−1)(u) + w(u, v)). Vale con cualquier signo. Costo O(VE).
• Lema de V − 1: sin ciclos negativos alcanzables, quitar ciclos no alarga, y un camino simple tiene a lo sumo V − 1 aristas.
• Una fila: d[v] tras la vuelta i no es d_i(v). Pero es el largo de algún camino, solo baja, y queda bajo d_i(v). Por eso nunca pasa de la tabla, y nunca baja de la distancia verdadera.
• Detección: la vuelta V relaja algo si y solo si se alcanza un ciclo negativo. Una dirección usa el teorema anterior, y la otra suma las desigualdades alrededor del ciclo, con d finitos porque se alcanza. Para todo el grafo: todos los d en 0.
• Recuperación: el nodo relajado en la vuelta V puede no estar en el ciclo. Se retrocede V veces y se lee. Los predecesores llevan a un ciclo, porque una cadena simple hasta s contradice que x se haya relajado, y sus ciclos son negativos, porque la relajación que los cierra es estricta.
• Infinito: INT_MAX + algo es overflow, comportamiento indefinido. Infinito + un negativo “mejora” a un nodo no alcanzado. De un nodo no alcanzado no se propaga nada.
• Vector de distancias: los routers usan esta relajación sin ver el mapa completo. RIP es un protocolo que usa esa idea.
• Todos los pares: V × Dijkstra si no hay negativos. V × Bellman-Ford, O(V²E). Floyd-Warshall, Θ(V³).
• Floyd-Warshall: d_k(i, j) con intermedios entre los k primeros nodos. d_(k+1)(i, j) = mín(d_k(i, j), d_k(i, k) + d_k(k, j)). Un camino simple usa a k una sola vez. En el mismo lugar es seguro porque sin ciclos negativos la fila y la columna k no cambian en la vuelta k.
• k por fuera: k es el parámetro del subproblema. Con k por dentro, el programa compila, da números creíbles, y está mal.
• Diagonal: hay un ciclo negativo si y solo si aparece d[k][k] < 0 al empezar la vuelta de k. Hay que detenerse ahí, porque los valores pueden crecer de forma exponencial.
• Cuadro de decisión: BFS sin pesos, orden topológico en un DAG, Dijkstra sin negativos, Bellman-Ford en cualquier grafo, Floyd-Warshall para todos los pares en un grafo pequeño o denso.
• ✂️ Held-Karp: C(S, v) guarda qué se visitó y dónde se está. O(n²·2ⁿ) en tiempo y O(n·2ⁿ) en memoria, frente a (n − 1)! recorridos. Con máscaras de bits, quitar un bit da un entero menor.
• Arbitraje: producto > 1 ⟺ suma de −ln(tasa) negativa. Con double, un producto exacto de 1 puede dar una suma negativa. Sumar ε a cada peso: L·ε bajo la ganancia prometida, y ε muy sobre el redondeo.

Ejercicios

1. A mano. Corran Bellman-Ford sobre el grafo de 4 nodos 0 → 1 (6), 0 → 2 (7), 1 → 3 (5), 2 → 1 (−3), 2 → 3 (9), con la tabla completa y con una sola fila, primero en este orden de aristas y después al revés. Comparen cuántas vueltas hacen falta hasta una que no cambia nada.
2. El lema. Den un grafo con un ciclo negativo que no se alcance desde s, en el que Bellman-Ford igual calcule bien δ(v) para todos los nodos alcanzables. Luego den uno con un ciclo negativo alcanzable y un nodo v que no se alcance desde el ciclo: ¿δ(v) existe? ¿Lo calcula bien la vuelta V − 1? Demuéstrenlo.
3. Todos en 0. Escriban la demostración completa de la versión con todos los nodos en 0: el significado de d_i(v), el teorema de una fila, y los dos sentidos del teorema de detección. Digan cuántas vueltas son, y por qué.
4. Ciclo de peso 0. Muestren que con un ciclo de peso 0 la vuelta V no relaja nada. ¿Dónde se usa que la relajación sea estricta? ¿Qué pasaría con <= en vez de <? Pruébenlo con bellman_ford.
5. El árbol de predecesores. Muestren un grafo sin ciclos negativos en el que los predecesores, en medio de una corrida, no sean todavía el árbol de caminos mínimos. Demuestren que al terminar, sin ciclos negativos, los predecesores forman un árbol con raíz en s cuyos caminos miden δ(v).
6. Floyd-Warshall, a mano. Calculen las matrices d_0, …, d_3 para las aristas 0 → 1 (3), 1 → 2 (−1), 2 → 0 (2), 2 → 1 (4). Comprueben con floyd_round.
7. Reconstruir. Agreguen a Floyd-Warshall una matriz next, con next[i][j] el primer nodo después de i en un camino mínimo a j. Demuestren que se puede mantener con una sola línea, y escriban una función que devuelva el camino. Pruébenla comparando su largo con la distancia, en grafos al azar.
8. Clausura transitiva. Con pesos booleanos, la misma recurrencia dice si hay un camino de i a j. Escriban el subproblema, demuéstrenlo, y digan por qué aquí no importa si hay ciclos. Compárenlo con correr BFS desde cada nodo: ¿cuál conviene en un grafo disperso?
9. k por dentro, más veces. Tomen floyd_warshall_k_inside. Si se repite completa dos veces sobre la misma matriz, ¿se arregla el ejemplo de este capítulo? Mídanlo con grafos al azar de 8 nodos, y reporten en cuántos falla con 1, 2 y 3 repeticiones. Un resultado así no reemplaza a una demostración: digan qué se puede concluir.
10. Potenciales. Sea h(v) el resultado de la versión con todos en 0, en un grafo sin ciclos negativos. Demuestren que w(u, v) + h(u) − h(v) ≥ 0 para toda arista, y que un camino de i a j cambia de largo exactamente en h(i) − h(j), sea cual sea. Expliquen por qué esto no es el error de “sumar una constante” de la semana 5, y por qué permite correr V × Dijkstra con pesos negativos. Den el costo total.
11. ✂️ Reconstruir el recorrido. Modifiquen held_karp_table para recuperar el recorrido, y no solo el costo, y comprueben que su costo sea el mínimo. ¿Cómo cambia la memoria si se guarda solo el costo? ¿Se puede guardar una sola “fila” como en la mochila, y seguir reconstruyendo?
12. Un viaje con ganancias. Tienen un mapa dirigido con un costo en cada tramo, que puede ser negativo si el tramo paga. Se quiere saber si desde la ciudad s se puede llegar a t ganando una cantidad tan grande como se quiera. Atención: el ciclo negativo tiene que alcanzarse desde s y tiene que llegar a t. Escriban el algoritmo, demuéstrenlo y pruébenlo con una fuerza bruta. ¿Qué hace la versión de “la vuelta V relaja algo” en un grafo con un ciclo que no llega a t?