{"id":159,"date":"2020-03-20T17:06:04","date_gmt":"2020-03-20T16:06:04","guid":{"rendered":"https:\/\/intsight.com\/?p=159"},"modified":"2021-01-21T17:50:03","modified_gmt":"2021-01-21T16:50:03","slug":"multiplicacion-de-matrices","status":"publish","type":"post","link":"https:\/\/intsight.com\/index.php\/2020\/03\/20\/multiplicacion-de-matrices\/","title":{"rendered":"Multiplicaci\u00f3n de matrices"},"content":{"rendered":"

Supongamos que queremos<\/span> multiplicar un par de matrices, $A$ y $B$. Digamos que la primera tiene dimensiones $m\\times n$ y que la segunda es $n\\times p$. La coincidencia entre columnas de la primera y filas de la segunda es condici\u00f3n necesaria para que podamos multiplicarlas.<\/p>\n

Si me piden que escriba de carrerilla un m\u00e9todo para esta multiplicaci\u00f3n, esto es lo que se me ocurre:<\/p>\n

public static double[,] Mult(double[,] a, double[,] b)\n{\n    int m = a.GetLength(0);\n    int n = a.GetLength(1);\n    int p = b.GetLength(1);\n    double[,] result = new double[m, p];\n    for (int i = 0; i < m; i++)\n        for (int j = 0; j < p; j++)\n        {\n            double d = 0;\n            for (int k = 0; k < n; k++)\n                d += a[i, k] * b[k, j];\n            result[i, j] = d;\n        }\n    return result;\n}<\/pre>\n
He utilizado matrices bidimensionales de C# porque acceder a sus elementos individuales es sencillo. Internamente, C# las almacena en una sola memoria contigua de memoria, fila por fila.<\/p>\n
El c\u00f3digo que he mostrado no es una maravilla. Para empezar, cada vez que decimos algo como a[i, k]<\/code>, el compilador tiene que multiplicar la variable i<\/code> por el n\u00famero de columnas y por los ocho bytes que tiene un flotante de doble precisi\u00f3n. Hacerlo una vez no es problema… pero tenemos tres bucles anidados. Eso tiene que doler. Si en vez de C# escribi\u00e9semos esto en C++, el compilador podr\u00eda sustituir un mont\u00f3n de multiplicaciones por sumas. RyuJIT ha mejorado much\u00edsimo, pero no tanto.<\/p>\n
C#, adem\u00e1s, es un lenguaje mucho m\u00e1s seguro que C++, pero esta seguridad nos cuesta un mont\u00f3n de verificaciones de rango para poder indexar. Recordemos, adem\u00e1s, que cada acceso necesita dos \u00edndices.<\/p>\n
Y hay un tercer problema, mucho m\u00e1s sutil: cuando las matrices son grandes, el c\u00f3digo anterior machaca la cach\u00e9 de la CPU sin piedad. Toma un folio de papel y haz el experimento: dibuja dos matrices, y ve numerando las celdas siguiendo el orden en que las usa el algoritmo.<\/p>\n
La clase Unsafe<\/h4>\n
Llegados a este punto, tenemos dos alternativas: o marcamos el m\u00e9todo como unsafe<\/strong> y usamos directamente punteros de C#, o intentamos evitarlo haciendo uso de la clase Unsafe<\/code>, de System.Runtime.CompilerServices<\/code>. Vamos a comenzar por esta \u00faltima. De paso, voy a invertir el orden de los dos bucles m\u00e1s internos, para ver qu\u00e9 conseguimos con ello. Este es el c\u00f3digo modificado, y suele funcionar el doble de r\u00e1pido, o un poco m\u00e1s:<\/p>\n
public static double[,] Mult(double[,] a, double[,] b)\n{\n    int m = a.GetLength(0);\n    int n = a.GetLength(1);\n    int p = b.GetLength(1);\n    double[,] c= new double[m, p];\n    ref double rA = ref a[0, 0];\n    ref double rB = ref b[0, 0];\n    ref double rC = ref c[0, 0];\n    for (int i = 0; i < m; i++)\n    {\n        ref double rAi = ref Unsafe.Add(ref rA, i * n);\n        ref double rCi = ref Unsafe.Add(ref rC, i * n);\n        for (int k = 0; k < n; k++)\n        {\n            double d = Unsafe.Add(ref rAi, k);\n            int kp = k * p;\n            for (int j = 0; j < p; j++)\n                Unsafe.Add(ref rCi, j) +=\n                    d * Unsafe.Add(ref rB, kp + j);\n        }\n    }\n    return c;\n}<\/pre>\n
La regla principal del uso de Unsafe.Add<\/code> es que si inicializamos as\u00ed:<\/p>\n
ref double rA = ref a[0, 0];<\/pre>\n
entonces el acceso a a[i, j]<\/code> debe parecerse a esto:<\/p>\n
Unsafe.Add(ref rA, i * n + j) = 42;<\/pre>\n
Esa multiplicaci\u00f3n es un problema del que ya advertimos. En nuestro c\u00f3digo lo paliamos moviendo la multiplicaci\u00f3n al inicio del bucle donde se le da valor al \u00edndice de la fila. Mi apa\u00f1o no es la palabra definitiva: le dejo como ejercicio la eliminaci\u00f3n total de esas multiplicaciones.<\/p>\n
Ahora hay que prestar atenci\u00f3n, sobre todo, al patr\u00f3n de acceso a memoria que se produce en el bucle m\u00e1s interno. En el algoritmo inicial, acumul\u00e1bamos todos los t\u00e9rminos de un elemento de la matriz final en el bucle interno, y asign\u00e1bamos su suma de golpe a la celda del resultado. Esta variante, sin embargo, no parece tan buena. Tenemos que asumir que, al reservar memoria para la matriz, todas sus entradas valen cero (y es as\u00ed). Luego, cada celda del resultado se va rellenando por pasos, no de una vez. Puede que esto sea bueno para la cach\u00e9 de la CPU, pero no me queda tan claro que sea bueno para el compilador de C#.<\/p>\n
Pero lo que nos interesa realmente es que ahora ejecutamos el siguiente patr\u00f3n de c\u00e1lculo:<\/p>\n
    \n
  1. Tenemos dos zonas de memoria consecutiva.<\/li>\n
  2. Leemos algo de la primera zona.<\/li>\n
  3. Lo transformamos como sea.<\/li>\n
  4. Lo asignamos a la celda equivalente en la segunda zona de memoria.<\/li>\n<\/ol>\n
    Instrucciones SIMD<\/h4>\n
    Ese patr\u00f3n de actividad es el t\u00edpico algoritmo \u00abvectorial\u00bb que podemos acelerar utilizando operaciones SIMD. Tenemos dos opciones:<\/p>\n
      \n
    • Utilizar System.Numerics.Vector<\/code>, que se adapta autom\u00e1ticamente a cualquier m\u00e1quina que soporte SIMD, e incluso ofrece una alternativa cuando no existe ese soporte. Este tipo funciona tambi\u00e9n para .NET Framework, a trav\u00e9s de un paquete.<\/li>\n
    • Si podemos usar .NET Core 3.1, podemos ir directamente a las clases declaradas en System.Runtime.Intrinsics<\/code> y System.Runtime.Intrinsics.X86<\/code>. Es un poco m\u00e1s complicado y no est\u00e1 bien documentado, pero da resultados ligeramente mejores.<\/li>\n<\/ul>\n
      Vamos a ir directamente por la segunda v\u00eda. Vamos a optimizar las CPUs que soporten el conjunto de instrucciones AVX, haremos algo m\u00e1s en el caso en que soporte el conjunto FMA (que mezcla multiplicaciones y sumas en una misma operaci\u00f3n) y, de todas maneras, habilitaremos c\u00f3digo de respaldo para cuando el procesador no soporte SIMD.<\/p>\n
      Cuando hay soporte para instrucciones AVX, podemos procesar hasta cuatro variables de tipo\u00a0double<\/code> de una tacada. Para ello tenemos que utilizar el tipo de estructura Vector256<\/code>, que tiene capacidad para cuatro elementos. La forma m\u00e1s sencilla de inicializar estos vectores es utilizando punteros, por lo que vamos a tener que declarar nuestro m\u00e9todo\u00a0unsafe<\/strong> y pasarnos directamente a los punteros.<\/p>\n
      public static unsafe double[,] Mult(double[,] a, double[,] b)\n{\n    int m = a.GetLength(0);\n    int n = a.GetLength(1);\n    int p = b.GetLength(1);\n    double[,] c = new double[m, p];\n    int lastBlockIndex = p - (p % 4);\n    fixed (double* pA = a)\n    fixed (double* pB = b)\n    fixed (double* pC = c)\n    {\n        double* pAi = pA;\n        double* pCi = pC;\n        for (int i = 0; i < m; i++)\n        {\n            double* pBk = pB;\n            for (int k = 0; k < n; k++)\n            {\n                double d = *(pAi + k);\n                if (Avx.IsSupported)\n                {\n                    int j = 0;\n                    var vd = Vector256.Create(d);\n                    while (j < lastBlockIndex)\n                    {\n                        if (Fma.IsSupported)\n                            Avx.Store(pCi + j,\n                                Fma.MultiplyAdd(\n                                Avx.LoadVector256(pBk + j),\n                                vd,\n                                Avx.LoadVector256(pCi + j)));\n                        else\n                            Avx.Store(pCi + j,\n                                Avx.Add(\n                                Avx.LoadVector256(pCi + j),\n                                Avx.Multiply(\n                                Avx.LoadVector256(pBk + j),\n                                vd)));\n                        j += 4;\n                    }\n                    while (j < p)\n                    {\n                        pCi[j] += d * pBk[j];\n                        j++;\n                    }\n                }\n                else\n                {\n                    for (int j = 0; j < p; j++)\n                        pCi[j] += d * pBk[j];\n                }\n                pBk += p;\n            }\n            pAi += n;\n            pCi += p;\n        }\n    }\n    return c;\n}<\/pre>\n
      Observaciones:<\/p>\n
        \n
      1. Lo peor de trabajar con SIMD es tener que lidiar con vectores que no son m\u00faltiplos exactos del tama\u00f1o del vector b\u00e1sico. Nuestros vectores b\u00e1sicos tienen cuatro elementos. Si tenemos un vector de 75 elementos, necesitaremos un bucle de 18 repeticiones que procese cuatro elementos por vez, para una mierdecilla de bucle final que maneje los 3 elementos que nos sobran.<\/li>\n
      2. Aunque la llamada a Avx.IsSupported<\/code> est\u00e1 metida dentro de dos bucles anidados, no se preocupe: el compilador JIT la trata como una constante en tiempo de generaci\u00f3n de c\u00f3digo nativo, y no cuesta nada. Si no se soporta AVX, el compilador JIT solamente genera el c\u00f3digo de la cl\u00e1usula else<\/strong>, que funciona sobre cualquier arquitectura.<\/li>\n
      3. Ojo: ese c\u00f3digo \u00abpara cualquier m\u00e1quina\u00bb podr\u00eda optimizarse echando mano de la t\u00e9cnica de\u00a0loop unrolling<\/em>. Pero mi pol\u00edtica en estos casos es: si no tienes una m\u00e1quina decente, j\u00f3dete.<\/li>\n
      4. En el ejemplo anterior, cuando intercambiamos el orden de los bucles m\u00e1s internos, ten\u00edamos un valor escalar que sac\u00e1bamos fuera del tercer bucle. Pero SIMD no ofrece instrucciones para multiplicar un vector por un escalar: tenemos que convertir ese escalar en todo un vector y utilizar la instrucci\u00f3n de multiplicaci\u00f3n m\u00e1s general. No es grave, de todos modos.<\/li>\n
      5. Si, adem\u00e1s de AVX, la m\u00e1quina soporta el conjunto FMA de instrucciones, podemos utilizar el m\u00e9todo MultiplyAdd<\/code> para acelerar un poco el algoritmo. Pero con esto hay que tener cuidado: a * b + c<\/code> puede dar resultados diferentes si se hacen las dos operaciones por separado o a la vez. Si se hacen a la vez, aumenta la exactitud de la operaci\u00f3n al existir menos redondeos. Pero el efecto secundario es que los c\u00e1lculos con y sin esa opci\u00f3n dan resultados ligeramente diferentes. Tenemos que decidir cu\u00e1ndo es aceptable que exista esa diferencia y cu\u00e1ndo no. En cualquier caso, tengamos presente que el resultado de MultiplyAdd<\/code> es m\u00e1s preciso.<\/li>\n<\/ol>\n
        Benchmark.NET<\/h4>\n
        Para estar seguro de las ganancias en velocidad, he utilizado el package Benchmark.NET para generar las pruebas. Estos son los resultados:<\/p>\n\n\n\n\n\n\n\n
        Method<\/th>\nMean<\/th>\nError<\/th>\nStdDev<\/th>\n<\/tr>\n<\/thead>\n
        MultMatrix<\/td>\n4,482.3 \u03bcs<\/td>\n88.75 \u03bcs<\/td>\n138.17 \u03bcs<\/td>\n<\/tr>\n
        UMultMatrix<\/td>\n1,895.2 \u03bcs<\/td>\n37.87 \u03bcs<\/td>\n63.26 \u03bcs<\/td>\n<\/tr>\n
        FMultMatrix<\/td>\n506.3 \u03bcs<\/td>\n3.44 \u03bcs<\/td>\n2.87 \u03bcs<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n
        La mejora por el uso de SIMD es cerca de cuatro veces, porque es el n\u00famero de operaciones simult\u00e1neas que permite esta arquitectura en particular. Con AVX512 tendr\u00edamos vectores de ocho valores, pero necesitar\u00edamos procesadores mucho m\u00e1s modernos, y de momento .NET Core no lo soporta.<\/p>\n
        Para esta prueba, he utilizado matrices de 128×128. He probado tambi\u00e9n con matrices de 8×8 e incluso de 4×4. La ganancia no es tan espectacular, pero en total se consigue una cuarta parte del tiempo de ejecuci\u00f3n respecto al algoritmo m\u00e1s sencillo.<\/p>\n","protected":false},"excerpt":{"rendered":"
        Supongamos que queremos multiplicar un par de matrices, $A$ y $B$. Digamos que la primera tiene dimensiones $m\\times n$ y que la segunda es $n\\times p$. La coincidencia entre columnas de la primera y filas de la segunda es condici\u00f3n necesaria para que podamos multiplicarlas. Si me piden que escriba de carrerilla un m\u00e9todo para […]<\/p>\n","protected":false},"author":1,"featured_media":161,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[4],"tags":[20,10,25,23,24],"class_list":["post-159","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-c","tag-matrices","tag-optimization","tag-pointers","tag-simd","tag-unsafe"],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/intsight.com\/wp-content\/uploads\/2020\/03\/img_0555.jpg?fit=540%2C275&ssl=1","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/159","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/comments?post=159"}],"version-history":[{"count":40,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/159\/revisions"}],"predecessor-version":[{"id":577,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/159\/revisions\/577"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media\/161"}],"wp:attachment":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media?parent=159"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/categories?post=159"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/tags?post=159"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}