Todos nosotros resolvemos un problema de “ camino más corto ” muchas veces al día. Sin querer, por supuesto. Lo resolvemos de camino al trabajo, cuando navegamos por Internet o cuando ordenamos nuestras cosas en el escritorio.
¿Suena un poco… demasiado? Averigüémoslo.
Imagina que has decidido reunirte con amigos, bueno... digamos en una cafetería. En primer lugar, tienes que encontrar una ruta (o camino) para llegar a la cafetería, así que empiezas a buscar transporte público disponible (si vas a pie) o rutas y calles (si vas en coche). Estás buscando una ruta desde tu ubicación actual hasta la cafetería, pero no “cualquier” ruta, sino la más corta o la más rápida.
He aquí otro ejemplo, que no es tan obvio como el anterior. Durante el trabajo, decides tomar un “atajo” por una calle lateral porque, bueno, es un “atajo” y es “más rápido” ir por ese camino. Pero, ¿cómo sabes que este “atajo” es más rápido? Basándote en tu experiencia personal, resolviste un problema de “camino más corto” y seleccionaste una ruta que pasa por la calle lateral.
En ambos ejemplos, la ruta “más corta” se determina en función de la distancia o el tiempo necesarios para llegar de un lugar a otro. Los ejemplos de viajes son aplicaciones muy naturales de la teoría de grafos y, en particular, del problema de la “ruta más corta”. Sin embargo, ¿qué sucede con el ejemplo de organizar cosas en un escritorio? En este caso, la “ruta más corta” puede representar una cantidad o complejidad de acciones que se deben realizar para obtener, por ejemplo, una hoja de papel: abrir el escritorio, abrir una carpeta, tomar una hoja de papel frente a tomar una hoja de papel directamente del escritorio.
Todos los ejemplos anteriores representan un problema de búsqueda de un camino más corto entre dos vértices en un grafo (una ruta o camino entre dos lugares, una cantidad de acciones o la complejidad de obtener una hoja de papel de un lugar u otro). Esta clase de problemas de camino más corto se denomina SSSP (Single Source Shortest Path) y el algoritmo base para resolver estos problemas es el algoritmo de Dijkstra , que tiene una complejidad computacional O(n^2)
.
Pero, a veces, necesitamos encontrar todos los caminos más cortos entre todos los vértices. Considere el siguiente ejemplo: está creando un mapa para sus movimientos regulares entre casa , trabajo y teatro . En este escenario, terminará con 6 rutas: work ⭢ home
, home ⭢ work
, work ⭢ theatre
, theatre ⭢ work
, home ⭢ theatre
y theatre ⭢ home
(las rutas inversas pueden ser diferentes debido a las carreteras de un solo sentido, por ejemplo).
Añadir más lugares al mapa dará como resultado un crecimiento significativo de combinaciones: según las permutaciones de la fórmula n de combinatoria, se puede calcular como:
A(k, n) = n! / (n - m)! // where n - is a number of elements, // k - is a number of elements in permutation (in our case k = 2)
Lo que nos da 12 rutas para 4 lugares y 90 rutas para 10 lugares (lo cual es impresionante). Nota… esto es sin considerar los puntos intermedios entre los lugares, es decir, para ir de casa al trabajo tienes que cruzar 4 calles, caminar a lo largo del río y cruzar el puente. Si imaginamos, algunas rutas pueden tener puntos intermedios comunes… bueno… como resultado tendremos un grafo muy grande, con muchos vértices, donde cada vértice representará un lugar o un punto intermedio significativo. La clase de problemas, donde necesitamos encontrar todos los caminos más cortos entre todos los pares de vértices en el grafo, se llama APSP (All Pairs Shortest Paths) y el algoritmo base para resolver estos problemas es el algoritmo de Floyd-Warshall , que tiene una complejidad computacional de O(n^3)
.
Y este es el algoritmo que implementaremos hoy 🙂
El algoritmo Floyd-Warshall encuentra todos los caminos más cortos entre cada par de vértices de un grafo. El algoritmo fue publicado por Robert Floyd en [1] (ver la sección “Referencias” para más detalles). Ese mismo año, Peter Ingerman en [2] describió una implementación moderna del algoritmo en forma de tres bucles for
anidados:
algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm // where W - is a weight matrix of N x N size, // min() - is a function which returns lesser of it's arguments
Si nunca ha tenido la oportunidad de trabajar con un gráfico representado en forma de matriz, puede resultarle difícil comprender qué hace el algoritmo anterior. Por lo tanto, para asegurarnos de que estamos en la misma página, veamos cómo se puede representar un gráfico en forma de matriz y por qué dicha representación es beneficiosa para resolver el problema de la ruta más corta.
La siguiente imagen ilustra un gráfico dirigido y ponderado de 5 vértices. A la izquierda, el gráfico se presenta en una forma visual, que está hecha de círculos y flechas, donde cada círculo representa un vértice y la flecha representa una arista con dirección. El número dentro de un círculo corresponde a un número de vértice y el número sobre una arista corresponde al peso de la arista. A la derecha, el mismo gráfico se presenta en forma de matriz de pesos. La matriz de pesos es una forma de matriz de adyacencia donde cada celda de la matriz contiene un "peso": una distancia entre el vértice i
(fila) y el vértice j
(columna). La matriz de pesos no incluye información sobre un "camino" entre vértices (una lista de vértices a través de los cuales se llega de i
a j
), solo un peso (distancia) entre estos vértices.
En una matriz de pesos podemos ver que los valores de las celdas son iguales a los pesos entre vértices adyacentes . Por eso, si inspeccionamos las rutas desde el vértice 0
(fila 0
), veremos que... solo hay una ruta: de 0
a 1
Pero, en una representación visual podemos ver claramente las rutas desde el vértice 0
a los vértices 2
y 3
(a través del vértice 1
). La razón de esto es simple: en un estado inicial, la matriz de pesos contiene solo la distancia entre vértices adyacentes. Sin embargo, esta información por sí sola es suficiente para encontrar el resto.
Veamos cómo funciona. Preste atención a la celda W[0, 1]
. Su valor indica que hay una ruta desde el vértice 0
al vértice 1
con un peso igual a 1
(en resumen, podemos escribirlo como: 0 ⭢ 1: 1
). Con este conocimiento, ahora podemos escanear todas las celdas de la fila 1
(que contiene todos los pesos de todas las rutas desde el vértice 1
) y volver a transferir esta información a la fila 0
, aumentando el peso en 1
(valor de W[0, 1]
).
Siguiendo los mismos pasos podemos encontrar caminos desde el vértice 0
a través de otros vértices. Durante la búsqueda, puede suceder que haya más de un camino que conduzca al mismo vértice y, lo que es más importante, los pesos de estos caminos pueden ser diferentes. Un ejemplo de una situación de este tipo se ilustra en la siguiente imagen, donde la búsqueda desde el vértice 0
hasta el vértice 2
reveló un camino más hacia el vértice 3
de un peso menor.
Tenemos dos caminos: un camino original – 0 ⭢ 3: 4
y un nuevo camino que acabamos de descubrir – 0 ⭢ 2 ⭢ 3: 3
(tenga en cuenta que la matriz de peso no contiene caminos, por lo que no sabemos qué vértices están incluidos en el camino original). Este es un momento en el que tomamos la decisión de mantener el más corto y escribir 3
en la celda W[0, 3]
.
¡Parece que acabamos de encontrar nuestro primer camino más corto!
Los pasos que acabamos de ver son la esencia del algoritmo Floyd-Warshall. Observemos una vez más el pseudocódigo del algoritmo:
algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm
El ciclo más externo for
en k
itera sobre todos los vértices en un grafo y en cada iteración, la variable k
representa un vértice a través del cual estamos buscando caminos . El ciclo interno for
en i
también itera sobre todos los vértices en el grafo y en cada iteración, i
representa un vértice desde el cual estamos buscando caminos . Y por último, un ciclo más interno for
en j
itera sobre todos los vértices en el grafo y en cada iteración, j
representa un vértice hacia el cual estamos buscando caminos. En combinación, nos da lo siguiente: en cada iteración k
estamos buscando caminos desde todos los vértices i
hacia todos los vértices j
a través del vértice k
. Dentro de un ciclo comparamos el camino i ⭢ j
(representado por W[i, j]
) con un camino i ⭢ k ⭢ j
(representado por la suma de W[I, k]
y W[k, j]
) y escribimos el más corto nuevamente en W[i, j]
.
Ahora, cuando entendemos la mecánica, es hora de implementar el algoritmo.
El código fuente y los gráficos experimentales están disponibles en el repositorio de GitHub. Los gráficos experimentales se pueden encontrar en el directorio Data/Sparse-Graphs.zip
. Todos los puntos de referencia de esta publicación se implementan en el archivo APSP01.cs .
Antes de sumergirnos en la implementación necesitamos aclarar algunos aspectos técnicos:
NO_EDGE = (int.MaxValue / 2) - 1
.
Ahora veamos por qué ocurre esto.
Respecto al punto 1, cuando hablamos de matrices definimos las celdas en términos de “filas” y “columnas”. Por ello, parece natural imaginar una matriz en forma de “cuadrado” o “rectángulo” (Imagen 4a).
Sin embargo, esto no significa necesariamente que tengamos que representar la matriz en forma de un arreglo de arreglos (Imagen 4b) para ceñirnos a nuestra imaginación. En lugar de esto, podemos representar una matriz en forma de un arreglo lineal (Imagen 4c) donde el índice de cada celda se calcula utilizando la siguiente fórmula:
i = row * row_size + col; // where row - cell row index, // col - cell column index, // row_size - number of cells in a row.
La matriz lineal de pesos es una condición previa para la ejecución eficaz del algoritmo Floyd-Warshall. La razón no es sencilla y una explicación detallada merecería una publicación aparte... o varias publicaciones. Sin embargo, actualmente es importante mencionar que dicha representación mejora significativamente la localidad de los datos , lo que de hecho tiene un gran impacto en el rendimiento del algoritmo.
Aquí les pido que crean en mí y que tengan esta información en mente como condición previa, sin embargo, al mismo tiempo les recomiendo que dediquen algún tiempo a estudiar la pregunta y, por cierto, no crean a la gente en Internet.
- Nota del autor
Respecto al punto n.° 2, si observas más de cerca el pseudocódigo del algoritmo, no encontrarás ninguna comprobación relacionada con la existencia de una ruta entre dos vértices. En su lugar, el pseudocódigo simplemente utiliza la función min()
. La razón es simple: originalmente, si no hay una ruta entre dos vértices, el valor de una celda se establece en infinity
y en todos los lenguajes, excepto quizás JavaScript, todos los valores son menores que infinity
. En el caso de los números enteros, puede resultar tentador utilizar int.MaxValue
como un valor de “sin ruta”. Sin embargo, esto provocará un desbordamiento de números enteros en los casos en que los valores de las rutas i ⭢ k
y k ⭢ j
sean iguales a int.MaxValue
. Por eso utilizamos un valor que es uno menos que la mitad de int.MaxValue
.
¡Oye! Pero, ¿por qué no podemos comprobar si la ruta existe antes de realizar ningún cálculo? Por ejemplo, comparando ambas rutas con 0 (si tomamos un cero como valor de “sin ruta”).
- Lector reflexivo
De hecho, es posible, pero lamentablemente esto conllevará una pérdida de rendimiento significativa. En resumen, la CPU mantiene estadísticas de los resultados de la evaluación de las ramas, por ejemplo, cuando algunas de las declaraciones if
se evalúan como true
o false
. Utiliza estas estadísticas para ejecutar el código de la “rama predicha estadísticamente” por adelantado, antes de que se evalúe la declaración if
real (esto se denomina ejecución especulativa ) y, por lo tanto, ejecutar el código de manera más eficiente. Sin embargo, cuando la predicción de la CPU es inexacta, causa una pérdida de rendimiento significativa en comparación con la predicción correcta y la ejecución incondicional porque la CPU tiene que detenerse y calcular la rama correcta.
Debido a que en cada iteración k
actualizamos una parte significativa de la matriz de ponderación, las estadísticas de las ramas de la CPU se vuelven inútiles porque no hay un patrón de código y todas las ramas se basan exclusivamente en datos. Por lo tanto, dicha verificación dará como resultado una cantidad significativa de predicciones erróneas de las ramas .
Aquí también les pido que me crean (por ahora) y luego dediquen algún tiempo a estudiar el tema.
Eh, parece que hemos terminado con la parte teórica: implementemos el algoritmo (denotamos esta implementación como Baseline
):
public void Baseline(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }
El código anterior es casi una copia idéntica del pseudocódigo mencionado anteriormente con una única excepción: en lugar de Math.Min()
usamos if
para actualizar una celda solo cuando sea necesario.
¡Oye! Espera un momento, ¿no fuiste tú quien escribió un montón de palabras sobre por qué los "si" no son buenos aquí y unas líneas más adelante nosotros mismos presentamos un "si"?
- Lector reflexivo
La razón es simple. En el momento de escribir, JIT emite un código casi equivalente para las implementaciones if
y Math.Min
. Puedes comprobarlo en detalle en sharplab.io , pero aquí tienes fragmentos de los cuerpos de los bucles principales:
// // == Assembly code of implementation of innermost loop for of j using if. // 53: movsxd r14, r14d // // compare matrix[i * sz + j] and distance (Condition) // 56: cmp [rdx+r14*4+0x10], ebx 5b: jle short 64 // // if matrix[i * sz + j] greater than distance write distance to matrix // 5d: movsxd rbp, ebp 60: mov [rdx+rbp*4+0x10], ebx 64: // end of loop // // == Assembly code of implementation of innermost loop for of j using Math.Min. // 4f: movsxd rbp, ebp 52: mov r14d, [rdx+rbp*4+0x10] // // compare matrix[i * sz + j] and distance (Condition) // 57: cmp r14d, ebx. // // if matrix[i * sz + j] less than distance write value to matrix // 5a: jle short 5e // // otherwise write distance to matrix // 5c: jmp short 61 5e: mov ebx, r14d 61: mov [rdx+rbp*4+0x10], ebx 65: // end of loop
Como puede ver, independientemente de si usamos if
o Math.Min
todavía hay una verificación condicional. Sin embargo, en el caso de if
no hay escritura innecesaria.
Ahora que hemos terminado con la implementación, ¡es hora de ejecutar el código y ver qué tan rápido es!
Puede verificar usted mismo la corrección del código ejecutando pruebas en un repositorio .
Utilizo Benchmark.NET (versión 0.12.1) para realizar pruebas comparativas de código. Todos los gráficos utilizados en las pruebas comparativas son gráficos acíclicos dirigidos de 300, 600, 1200, 2400 y 4800 vértices. La cantidad de aristas en los gráficos es de alrededor del 80 % del máximo posible, que para los gráficos acíclicos dirigidos se puede calcular como:
var max = v * (v - 1)) / 2; // where v - is a number of vertexes in a graph.
¡Vamos a rockear!
Estos son los resultados de los puntos de referencia ejecutados en mi PC (Windows 10.0.19042, CPU Intel Core i7-7700 3.60Ghz (Kaby Lake) / 8 procesadores lógicos / 4 núcleos):
Método | Tamaño | Significar | Error | Desviación estándar |
---|---|---|---|---|
Base | 300 | 27,525 ms | 0,1937 ms | 0,1617 ms |
Base | 600 | 217,897 ms | 1,6415 ms | 1,5355 ms |
Base | 1200 | 1.763,335 ms | 7,4561 ms | 6,2262 ms |
Base | 2400 | 14.533,335 ms | 63,3518 ms | 52,9016 ms |
Base | 4800 | 119.768,219 ms | 181,5514 ms | 160,9406 ms |
De los resultados podemos ver que el tiempo de cálculo crece dramáticamente comparado con el tamaño de un gráfico: para un gráfico de 300 vértices tomó 27 milisegundos, para un gráfico de 2400 vértices, tomó 14.5 segundos y para un gráfico de 4800, 119 segundos, lo que es casi 2 minutos .
Mirando el código del algoritmo puede ser difícil imaginar que hay algo que podemos hacer para acelerar los cálculos… porque bueno… hay TRES bucles, sólo TRES bucles.
Sin embargo, como suele ocurrir, las posibilidades se esconden en los detalles 🙂
El algoritmo de Floyd-Warshall es un algoritmo base para resolver problemas de ruta más corta de todos los pares, especialmente cuando se trata de gráficos densos o completos (porque el algoritmo busca rutas entre todos los pares de vértices).
Sin embargo, en nuestros experimentos utilizamos grafos acíclicos dirigidos , que tienen una propiedad maravillosa: si hay un camino desde el vértice 1
al vértice 2
, entonces no hay camino desde el vértice 2
al vértice 1
Para nosotros, significa que hay muchos vértices no adyacentes que podemos omitir si no hay un camino desde i
hasta k
(denotamos esta implementación como SpartialOptimisation
).
public void SpartialOptimisation(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }
Aquí se muestran los resultados de las implementaciones anteriores ( Baseline
) y actuales ( SpartialOptimisation
) en el mismo conjunto de gráficos (los resultados más rápidos están resaltados en negrita ):
Método | Tamaño | Significar | Error | Desviación estándar | Relación |
---|---|---|---|---|---|
Base | 300 | 27,525 ms | 0,1937 ms | 0,1617 ms | 1.00 |
Optimización parcial | 300 | 12,399 ms | 0,0943 ms | 0,0882 ms | 0,45 |
Base | 600 | 217,897 ms | 1,6415 ms | 1,5355 ms | 1.00 |
Optimización parcial | 600 | 99,122 ms | 0,8230 ms | 0,7698 ms | 0,45 |
Base | 1200 | 1.763,335 ms | 7,4561 ms | 6,2262 ms | 1.00 |
Optimización parcial | 1200 | 766,675 ms | 6,1147 ms | 5,7197 ms | 0,43 |
Base | 2400 | 14.533,335 ms | 63,3518 ms | 52,9016 ms | 1.00 |
Optimización parcial | 2400 | 6.507,878 ms | 28,2317 ms | 26,4079 ms | 0,45 |
Base | 4800 | 119.768,219 ms | 181,5514 ms | 160,9406 ms | 1.00 |
Optimización parcial | 4800 | 55.590,374 ms | 414,6051 ms | 387,8218 ms | 0,46 |
¡Bastante impresionante!
¡Hemos reducido el tiempo de cálculo a la mitad! Por supuesto, cuanto más denso sea el gráfico, menor será la velocidad. Sin embargo, esta es una de las optimizaciones que resulta útil si sabes de antemano con qué clase de datos vas a trabajar.
¿Podemos hacer más? 🙂
Mi PC está equipada con un procesador Inter Core i7-7700 CPU 3.60GHz (Kaby Lake)
con 8 procesadores lógicos ( HW ) o 4 núcleos con tecnología Hyper-Threading . Tener más de un núcleo es como tener más “manos libres” que podemos poner a trabajar. Así que veamos qué parte del trabajo se puede dividir eficientemente entre varios trabajadores y luego, paralelizarlo.
Los bucles son siempre los candidatos más obvios para la paralelización porque en la mayoría de los casos las iteraciones son independientes y pueden ejecutarse simultáneamente. En el algoritmo, tenemos tres bucles que debemos analizar y ver si hay dependencias que nos impidan paralelizarlos.
Empecemos con el bucle for of k
. En cada iteración, el bucle calcula las rutas desde cada vértice a cada vértice a través del vértice k
. Las rutas nuevas y actualizadas se almacenan en la matriz de pesos. Al observar las iteraciones, es posible que notes que se pueden ejecutar en cualquier orden: 0, 1, 2, 3 o 3, 1, 2, 0 sin afectar el resultado. Sin embargo, aún deben ejecutarse en secuencia; de lo contrario, algunas de las iteraciones no podrán usar rutas nuevas o actualizadas porque aún no se escribirán en la matriz de pesos. Tal inconsistencia definitivamente arruinará los resultados. Por lo tanto, debemos seguir buscando.
El siguiente candidato es for of i
. En cada iteración, el bucle lee una ruta desde el vértice i
al vértice k
(celda: W[i, k]
), una ruta desde el vértice k
al vértice j
(celda: W[k, j
]) y luego verifica si la ruta conocida desde i
a j
(celda: W[i, j]
) es más corta que la ruta i ⭢ k ⭢ j
(suma de: W[i, k]
+ W[k, j]
). Si observa más de cerca el patrón de acceso, puede notar que, en cada iteración, el bucle i
lee datos de la fila k
y actualiza la fila i
lo que básicamente significa que las iteraciones son independientes y se pueden ejecutar no solo en cualquier orden, sino también en paralelo.
Esto parece prometedor, así que vamos a implementarlo (denotamos esta implementación como SpartialParallelOptimisations
).
El bucle for of j
también se puede paralelizar. Sin embargo, la paralelización del ciclo más interno en este caso es muy ineficiente. Puedes comprobarlo tú mismo haciendo unos pocos cambios simples en el código fuente .
- Nota del autor
public void SpartialParallelOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } }
Aquí están los resultados de las implementaciones de Baseline
, SpartialOptimisation
y SpartialParallelOptimisations
en el mismo conjunto de gráficos (la paralelización se realiza utilizando la clase Parallel ):
Método | Tamaño | Significar | Error | Desviación estándar | Relación |
---|---|---|---|---|---|
Base | 300 | 27,525 ms | 0,1937 ms | 0,1617 ms | 1.00 |
Optimización parcial | 300 | 12,399 ms | 0,0943 ms | 0,0882 ms | 0,45 |
Optimizaciones paralelas parciales | 300 | 6,252 ms | 0,0211 ms | 0,0187 ms | 0,23 |
Base | 600 | 217,897 ms | 1,6415 ms | 1,5355 ms | 1.00 |
Optimización parcial | 600 | 99,122 ms | 0,8230 ms | 0,7698 ms | 0,45 |
Optimizaciones paralelas parciales | 600 | 35,825 ms | 0,1003 ms | 0,0837 ms | 0,16 |
Base | 1200 | 1.763,335 ms | 7,4561 ms | 6,2262 ms | 1.00 |
Optimización parcial | 1200 | 766,675 ms | 6,1147 ms | 5,7197 ms | 0,43 |
Optimizaciones paralelas parciales | 1200 | 248,801 ms | 0,6040 ms | 0,5043 ms | 0,14 |
Base | 2400 | 14.533,335 ms | 63,3518 ms | 52,9016 ms | 1.00 |
Optimización parcial | 2400 | 6.507,878 ms | 28,2317 ms | 26,4079 ms | 0,45 |
Optimizaciones paralelas parciales | 2400 | 2.076,403 ms | 20,8320 ms | 19,4863 ms | 0,14 |
Base | 4800 | 119.768,219 ms | 181,5514 ms | 160,9406 ms | 1.00 |
Optimización parcial | 4800 | 55.590,374 ms | 414,6051 ms | 387,8218 ms | 0,46 |
Optimizaciones paralelas parciales | 4800 | 15.614,506 ms | 115,6996 ms | 102,5647 ms | 0,13 |
De los resultados podemos ver que la paralelización del bucle for of i
ha reducido el tiempo de cálculo entre 2 y 4 veces en comparación con la implementación anterior ( SpartialOptimisation
). Esto es muy impresionante, sin embargo, es importante recordar que la paralelización de cálculos puros consume todos los recursos computacionales disponibles, lo que puede provocar la escasez de recursos de otras aplicaciones en el sistema.
Como puedes imaginar, este no es el final 🙂
La vectorización es una transformación de un código que opera sobre un solo elemento en un código que opera sobre múltiples elementos simultáneamente.
Puede parecer un asunto complejo, así que veamos cómo funciona con un ejemplo sencillo:
var a = new [] { 5, 7, 8, 1 }; var b = new [] { 4, 2, 2, 6 }; var c = new [] { 0, 0, 0, 0 }; for (var i = 0; i < 4; ++i) c[i] = a[i] + b[i];
De manera muy simplificada , la ejecución de la iteración 0
del bucle for
anterior a nivel de CPU se puede ilustrar de la siguiente manera:
Esto es lo que sucede. La CPU carga los valores de las matrices a
y b
desde la memoria a los registros de la CPU (cada registro puede contener exactamente un valor). Luego, la CPU ejecuta la operación de suma escalar en estos registros y escribe el resultado nuevamente en la memoria principal, directamente en una matriz c
. Este proceso se repite cuatro veces antes de que finalice el bucle.
Vectorización significa el uso de registros especiales de CPU: registros vectoriales o SIMD (instrucción única, múltiples datos) e instrucciones de CPU correspondientes para ejecutar operaciones en múltiples valores de matriz a la vez:
Esto es lo que sucede. La CPU carga los valores de las matrices a
y b
desde la memoria a los registros de la CPU (sin embargo, esta vez, un registro vectorial puede contener dos valores de matriz). Luego, la CPU ejecuta la operación de suma vectorial en estos registros y escribe el resultado nuevamente en la memoria principal, directamente en una matriz c
. Debido a que operamos con dos valores a la vez, este proceso se repite dos veces en lugar de cuatro.
Para implementar la vectorización en .NET podemos usar los tipos Vector y Vector<T> (también podemos usar tipos del espacio de nombres System.Runtime.Intrinsics , sin embargo son un poco avanzados para la implementación actual, por lo que no los usaré, pero siéntete libre de probarlos tú mismo):
public void SpartialVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }
El código vectorizado puede parecer un poco extraño, así que veámoslo paso a paso.
Creamos un vector de i ⭢ k
path: var ik_vec = new Vector<int>(matrix[i * sz + k])
. Como resultado, si vector puede contener cuatro valores de tipo int
y el peso de i ⭢ k
path es igual a 5, crearemos un vector de cuatro 5: [5, 5, 5, 5]. A continuación, en cada iteración, calculamos simultáneamente las rutas desde el vértice i
hasta los vértices j, j + 1, ..., j + Vector<int>.Count
.
La propiedad Vector<int>.Count
devuelve la cantidad de elementos de tipo int
que caben en los registros vectoriales. El tamaño de los registros vectoriales depende del modelo de CPU, por lo que esta propiedad puede devolver valores diferentes en distintas CPU.
- Nota del autor
Todo el proceso de cálculo se puede dividir en tres pasos:
ij_vec
e ikj_vec
.ij_vec
e ikj_vec
y seleccione los valores más pequeños en un nuevo vector r_vec
.r_vec
nuevamente en la matriz de peso.
Si bien los puntos 1 y 3 son bastante sencillos, el punto 2 requiere una explicación. Como se mencionó anteriormente, con los vectores manipulamos múltiples valores al mismo tiempo. Por lo tanto, el resultado de algunas operaciones no puede ser singular, es decir, la operación de comparación Vector.LessThan(ij_vec, ikj_vec)
no puede devolver true
o false
porque compara múltiples valores. Por lo tanto, devuelve un nuevo vector que contiene los resultados de la comparación entre los valores correspondientes de los vectores ij_vec
e ikj_vec
( -1
, si el valor de ij_vec
es menor que el valor de ikj_vec
y 0
en caso contrario). El vector devuelto (por sí solo) no es realmente útil, sin embargo, podemos usar la operación de vector Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec)
para extraer los valores requeridos de los vectores ij_vec
e ikj_vec
en un nuevo vector: r_vec
. Esta operación devuelve un nuevo vector donde los valores son iguales al menor de dos valores correspondientes de los vectores de entrada, es decir, si el valor del vector lt_vec
es igual a -1
, entonces la operación selecciona un valor de ij_vec
; de lo contrario, selecciona un valor de ikj_vec
.
Además del punto 2 , hay una parte más que requiere explicación: el segundo bucle no vectorizado. Es necesario cuando el tamaño de la matriz de peso no es un producto del valor de Vector<int>.Count
. En ese caso, el bucle principal no puede procesar todos los valores (porque no se puede cargar parte del vector) y el segundo bucle no vectorizado calculará las iteraciones restantes.
Aquí están los resultados de esta y todas las implementaciones anteriores:
Método | tal | Significar | Error | Desviación estándar | Relación |
---|---|---|---|---|---|
Base | 300 | 27,525 ms | 0,1937 ms | 0,1617 ms | 1.00 |
Optimización parcial | 300 | 12,399 ms | 0,0943 ms | 0,0882 ms | 0,45 |
Optimizaciones paralelas parciales | 300 | 6,252 ms | 0,0211 ms | 0,0187 ms | 0,23 |
Optimizaciones de vectores parciales | 300 | 3,056 ms | 0,0301 ms | 0,0281 ms | 0,11 |
Base | 600 | 217,897 ms | 1,6415 ms | 1,5355 ms | 1.00 |
Optimización parcial | 600 | 99,122 ms | 0,8230 ms | 0,7698 ms | 0,45 |
Optimizaciones paralelas parciales | 600 | 35,825 ms | 0,1003 ms | 0,0837 ms | 0,16 |
Optimizaciones de vectores parciales | 600 | 24,378 ms | 0,4320 ms | 0,4041 ms | 0,11 |
Base | 1200 | 1.763,335 ms | 7,4561 ms | 6,2262 ms | 1.00 |
Optimización parcial | 1200 | 766,675 ms | 6,1147 ms | 5,7197 ms | 0,43 |
Optimizaciones paralelas parciales | 1200 | 248,801 ms | 0,6040 ms | 0,5043 ms | 0,14 |
Optimizaciones de vectores parciales | 1200 | 185,628 ms | 2,1240 ms | 1,9868 ms | 0,11 |
Base | 2400 | 14.533,335 ms | 63,3518 ms | 52,9016 ms | 1.00 |
Optimización parcial | 2400 | 6.507,878 ms | 28,2317 ms | 26,4079 ms | 0,45 |
Optimizaciones paralelas parciales | 2400 | 2.076,403 ms | 20,8320 ms | 19,4863 ms | 0,14 |
Optimizaciones de vectores parciales | 2400 | 2.568,676 ms | 31,7359 ms | 29,6858 ms | 0,18 |
Base | 4800 | 119.768,219 ms | 181,5514 ms | 160,9406 ms | 1.00 |
Optimización parcial | 4800 | 55.590,374 ms | 414,6051 ms | 387,8218 ms | 0,46 |
Optimizaciones paralelas parciales | 4800 | 15.614,506 ms | 115,6996 ms | 102,5647 ms | 0,13 |
Optimizaciones de vectores parciales | 4800 | 18.257,991 ms | 84,5978 ms | 79,1329 ms | 0,15 |
Del resultado se desprende claramente que la vectorización ha reducido significativamente el tiempo de cálculo: de 3 a 4 veces en comparación con la versión no paralelizada ( SpartialOptimisation
). Lo interesante aquí es que la versión vectorizada también supera a la versión paralela ( SpartialParallelOptimisations
) en gráficos más pequeños (menos de 2400 vértices).
Y por último, pero no menos importante: ¡combinemos vectorización y paralelismo!
Si está interesado en la aplicación práctica de la vectorización, le recomiendo leer la serie de publicaciones de Dan Shechter donde vectorizó Array.Sort
(estos resultados luego aparecieron en la actualización de Garbage Collector en .NET 5 ).
La última implementación combina esfuerzos de paralelización y vectorización y… para ser honestos, carece de individualidad 🙂 Porque… bueno, acabamos de reemplazar un cuerpo de Parallel.For
de SpartialParallelOptimisations
con un bucle vectorizado de SpartialVectorOptimisations
:
public void SpartialParallelVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } }
A continuación se presentan todos los resultados de este artículo. Como era de esperar, el uso simultáneo de paralelismo y vectorización demostró los mejores resultados, reduciendo el tiempo de cálculo hasta 25 veces (para gráficos de 1200 vértices) en comparación con la implementación inicial.
Método | tal | Significar | Error | Desviación estándar | Relación |
---|---|---|---|---|---|
Base | 300 | 27,525 ms | 0,1937 ms | 0,1617 ms | 1.00 |
Optimización parcial | 300 | 12,399 ms | 0,0943 ms | 0,0882 ms | 0,45 |
Optimizaciones paralelas parciales | 300 | 6,252 ms | 0,0211 ms | 0,0187 ms | 0,23 |
Optimizaciones de vectores parciales | 300 | 3,056 ms | 0,0301 ms | 0,0281 ms | 0,11 |
Optimizaciones de vectores paralelos parciales | 300 | 3,008 ms | 0,0075 ms | 0,0066 ms | 0,11 |
Base | 600 | 217,897 ms | 1,6415 ms | 1,5355 ms | 1.00 |
Optimización parcial | 600 | 99,122 ms | 0,8230 ms | 0,7698 ms | 0,45 |
Optimizaciones paralelas parciales | 600 | 35,825 ms | 0,1003 ms | 0,0837 ms | 0,16 |
Optimizaciones de vectores parciales | 600 | 24,378 ms | 0,4320 ms | 0,4041 ms | 0,11 |
Optimizaciones de vectores paralelos parciales | 600 | 13,425 ms | 0,0319 ms | 0,0283 ms | 0,06 |
Base | 1200 | 1.763,335 ms | 7,4561 ms | 6,2262 ms | 1.00 |
Optimización parcial | 1200 | 766,675 ms | 6,1147 ms | 5,7197 ms | 0,43 |
Optimizaciones paralelas parciales | 1200 | 248,801 ms | 0,6040 ms | 0,5043 ms | 0,14 |
Optimizaciones de vectores parciales | 1200 | 185,628 ms | 2,1240 ms | 1,9868 ms | 0,11 |
Optimizaciones de vectores paralelos parciales | 1200 | 76,770 ms | 0,3021 ms | 0,2522 ms | 0,04 |
Base | 2400 | 14.533,335 ms | 63,3518 ms | 52,9016 ms | 1.00 |
Optimización parcial | 2400 | 6.507,878 ms | 28,2317 ms | 26,4079 ms | 0,45 |
Optimizaciones paralelas parciales | 2400 | 2.076,403 ms | 20,8320 ms | 19,4863 ms | 0,14 |
Optimizaciones de vectores parciales | 2400 | 2.568,676 ms | 31,7359 ms | 29,6858 ms | 0,18 |
Optimizaciones de vectores paralelos parciales | 2400 | 1.281,877 ms | 25,1127 ms | 64,8239 ms | 0,09 |
Base | 4800 | 119.768,219 ms | 181,5514 ms | 160,9406 ms | 1.00 |
Optimización parcial | 4800 | 55.590,374 ms | 414,6051 ms | 387,8218 ms | 0,46 |
Optimizaciones paralelas parciales | 4800 | 15.614,506 ms | 115,6996 ms | 102,5647 ms | 0,13 |
Optimizaciones de vectores parciales | 4800 | 18.257,991 ms | 84,5978 ms | 79,1329 ms | 0,15 |
Optimizaciones de vectores paralelos parciales | 4800 | 12.785,824 ms | 98,6947 ms | 87,4903 ms | 0,11 |
En esta publicación, analizamos en profundidad un problema de ruta más corta entre pares e implementamos un algoritmo Floyd-Warshall en C# para resolverlo. También actualizamos nuestra implementación para respetar los datos y utilizar características de bajo nivel de .NET y hardware.
En esta publicación, jugamos a lo grande. Sin embargo, en aplicaciones de la vida real es importante recordar que no estamos solos. La paralelización extrema puede degradar significativamente el rendimiento del sistema y causar más daño que beneficio. La vectorización, por otro lado, es un poco más segura si se realiza de manera independiente de la plataforma. Recuerde que algunas de las instrucciones vectoriales podrían estar disponibles solo en ciertas CPU.
Espero que hayas disfrutado la lectura y te hayas divertido “exprimiendo” un poco el rendimiento de un algoritmo con solo TRES ciclos 🙂