{"id":1301,"date":"2023-08-07T13:11:26","date_gmt":"2023-08-07T11:11:26","guid":{"rendered":"https:\/\/intsight.com\/?p=1301"},"modified":"2023-08-07T13:12:26","modified_gmt":"2023-08-07T11:12:26","slug":"el-algoritmo-de-welford","status":"publish","type":"post","link":"https:\/\/intsight.com\/index.php\/2023\/08\/07\/el-algoritmo-de-welford\/","title":{"rendered":"El algoritmo de Welford"},"content":{"rendered":"

En la entrada sobre<\/span> la varianza<\/a>, vimos que pod\u00edamos tener problemas de estabilidad num\u00e9rica si intent\u00e1bamos calcular la varianza en una sola pasada sobre los datos, usando inocentemente la definici\u00f3n matem\u00e1tica. La soluci\u00f3n de entonces fue usar un algoritmo de dos pasos: calcular la media en el primer paso, y en el segundo, calcular la varianza de la muestra menos la media. Hab\u00eda otra posibilidad: usar el primer valor de la secuencia como estimado malo de la media, y restar ese valor a las sucesivas muestras.<\/p>\n

\u00bfPodr\u00edamos hacer algo mejor si corrigi\u00e9semos el estimado de la media sobre la marcha? Resulta que se puede, y el primero en darse cuenta fue Welford<\/a>, all\u00e1 por el 1962. Donald Knuth incluy\u00f3 el algoritmo en el segundo tomo de The Art of Computer Programming<\/a>. El algoritmo original s\u00f3lo calculaba la media y la varianza sobre la marcha, pero Timothy Terriberry<\/a>, en 2007, lo ampli\u00f3 para que calculase momentos superiores. Este algoritmo est\u00e1 implementado, por ejemplo, en Math.Net Numerics<\/a>, aunque la implementaci\u00f3n es mejorable.<\/p>\n

\u00bfQu\u00e9 nos da?<\/h4>\n

En AUSTRA, la clase que implementa este algoritmo se llama Accumulator<\/code>. Hay tambi\u00e9n una clase simplificada, SimpleAccumulator<\/code>, que s\u00f3lo calcula los dos primeros momentos, con el beneficio evidente de tener que ejecutar menos trabajo.<\/p>\n

La definici\u00f3n de Accumulator<\/code>, junto con su constructor principal y los campos y propiedades de almacenamiento, es la siguiente:<\/p>\n

\n\/\/\/ Calculates statistics by adding samples.<\/summary>\npublic sealed class Accumulator\n{\n    \/\/\/Minimum value.<\/summary>\n    private double min = double.PositiveInfinity;\n    \/\/\/Maximum value.<\/summary>\n    private double max = double.NegativeInfinity;\n    \/\/\/Estimated mean.<\/summary>\n    private double m1;\n    \/\/\/Accumulated second moment.<\/summary>\n    private double m2;\n    \/\/\/Accumulated third moment.<\/summary>\n    private double m3;\n    \/\/\/Accumulated fourth moment.<\/summary>\n    private double m4;\n\n    \/\/\/Gets the total number of samples.<\/summary>\n    public long Count { get; private set; }\n\n    \/\/\/Creates an empty accumulator.<\/summary>\n    public Accumulator() { }\n\n    \/* … *\/\n}\n<\/pre>\n
La informaci\u00f3n que nos interesa se obtiene a trav\u00e9s de estos campos, por medio de propiedades calculadas:<\/p>\n
\n\/\/\/Returns the minimum value.<\/summary>\npublic double Minimum => Count > 0 ? min : double.NaN;\n\/\/\/Returns the maximum value.<\/summary>\npublic double Maximum => Count > 0 ? max : double.NaN;\n\/\/\/Gets the sample mean.<\/summary>\npublic double Mean => Count > 0 ? m1 : double.NaN;\n\/\/\/Gets the unbiased variance.<\/summary>\npublic double Variance =>\n    Count < 2 ? double.NaN : m2 \/ (Count - 1);\n\/\/\/Gets the unbiased standard deviation.<\/summary>\npublic double StandardDeviation =>\n    Count < 2 ? double.NaN : Sqrt(m2 \/ (Count - 1));\n\/\/\/Gets the unbiased population skewness.<\/summary>\npublic double Skewness =>\n    Count < 3\n    ? double.NaN\n    : Count * m3 * Sqrt(m2 \/ (Count - 1))\n        \/ (m2 * m2 * (Count - 2)) * (Count - 1);\n\/\/\/Gets the unbiased population kurtosis.<\/summary>\npublic double Kurtosis =>\n    Count < 4\n    ? double.NaN\n    : ((double)Count * Count - 1)\n        \/ ((Count - 2) * (Count - 3))\n        * (Count * m4 \/ (m2 * m2) - 3 + 6.0 \/ (Count + 1));\n<\/pre>\n
He omitido, por brevedad, el c\u00e1lculo de propiedades como PopulationVariance<\/code>, PopulationSkewness<\/code> y dem\u00e1s. De todas maneras, est\u00e1n disponibles en el c\u00f3digo de Austra, y en el propio c\u00f3digo de Math.NET Numerics.<\/p>\n
\u00bfQu\u00e9 le damos?<\/h4>\n
Para alimentar al acumulador, hay que pasarle las muestras, al menos en principio, de una en una, por medio de un m\u00e9todo que hemos nombrado Add<\/code>:<\/p>\n
\n\/\/\/ Adds a sample to this accumulator.<\/summary>\n\/\/\/ The new sample.<\/param>\npublic void Add(double sample)\n{\n    ++Count;\n    double d = sample - m1, s = d \/ Count;\n    double t = d * s * (Count - 1);\n    m1 += s;\n    m4 += (t * s * (Count * (Count - 3 + 3)\n        + 6 * s * m2 - 4 * m3) * s;\n    m3 += (t * (Count - 2) - 3 * m2) * s;\n    m2 += t;\n    if (sample < min) min = sample;\n    if (sample > max) max = sample;\n}\n<\/pre>\n
Hay algunas peque\u00f1as mejoras en el c\u00f3digo anterior, respecto al original. Hay algunas multiplicaciones menos, y est\u00e1 todo preparado por si quisi\u00e9ramos usar alguna instrucci\u00f3n de fusi\u00f3n de multiplicaci\u00f3n y suma. No las he usado porque tengo algunas dudas sobre la eficiencia en .NET Core. Es cierto que siempre tienes la ventaja de la mayor exactitud, pero ya veremos d\u00f3nde s\u00ed se usan (en pocas palabras: donde realmente importan).<\/p>\n
Combinando acumuladores<\/h4>\n
Lo mejor de todo es que podemos combinar los valores en dos acumuladores independientes y generar un acumulador de los datos conjuntos. Esto nos permitir\u00eda, por ejemplo, dividir una muestra grande en cuatro partes, calcular cuatro acumuladores en paralelo, y luego mezclarlos en el resultado final.<\/p>\n
\npublic static Accumulator operator +(\n    Accumulator a1, Accumulator a2)\n{\n    if (a1.Count == 0) return a2;\n    if (a2.Count == 0) return a1;\n\n    long n = a1.Count + a2.Count, n2 = n * n;\n    double d = a2.m1 - a1.m1, d2 = d * d;\n    double d3 = d2 * d, d4 = d2 * d2;\n    double m1 = (a1.Count * a1.m1 + a2.Count * a2.m1) \/ n;\n    double m2 = a1.m2 + a2.m2 + d2 * a1.Count * a2.Count \/ n;\n    double m3 = a1.m3 + a2.m3\n        + d3 * a1.Count * a2.Count * (a1.Count - a2.Count) \/ n2\n        + 3 * d * (a1.Count * a2.m2 - a2.Count * a1.m2) \/ n;\n    double m4 = a1.m4 + a2.m4 + d4 * a1.Count * a2.Count\n            * (a1.Count * (a1.Count - a2.Count)\n                + a2.Count * a2.Count) \/ (n2 * n)\n        + 6 * d2 * (a1.Count * a1.Count * a2.m2\n            + a2.Count * a2.Count * a1.m2) \/ n2\n        + 4 * d * (a1.Count * a2.m3 - a2.Count * a1.m3) \/ n;\n    return new() {\n        Count = n,\n        m1 = m1, m2 = m2, m3 = m3, m4 = m4,\n        min = Min(a1.min, a2.min),\n        max = Max(a1.max, a2.max),\n    };\n}\n<\/pre>\n
El c\u00f3digo es complicado, y es f\u00e1cil equivocarse copiando y pegando (ya me ha pasado). De todas maneras, es una clase que es f\u00e1cil de testear.<\/p>\n
El primer impulso es dejarlo aqu\u00ed, y confiar en el paralelismo con tareas para cuando queremos acelerar el c\u00f3digo. Mi problema con esto es que, cuando escribo c\u00f3digo para una librer\u00eda, prefiero realizar la aceleraci\u00f3n b\u00e1sica con c\u00f3digo vectorial (AVX o lo que est\u00e9 disponible). \u00bfPor qu\u00e9? Pues porque por experiencia, el programador que usa luego la biblioteca prefiere tener la opci\u00f3n del paralelismo por tareas para su propio c\u00f3digo<\/em>. Es cierto que las tareas se combinan m\u00e1s o menos bien en .NET, gracias al thread-pool<\/em> que obtienes del entorno de ejecuci\u00f3n, sin esforzarte demasiado; en Java, todo es m\u00e1s complicado con sus malditos executors<\/em>.<\/p>\n
Drum roll, please...<\/h4>\n
Prefiero, por lo tanto, ganar todo lo que pueda en paralelismo en una librer\u00eda a golpe de instrucciones SIMD. Y esto es precisamente lo que hacemos en Accumulator<\/code> con el siguiente m\u00e9todo y algunos m\u00e1s que lo llaman:<\/p>\n
\npublic unsafe void Add(double* samples, int size)\n{\n    int i = 0;\n    if (Avx.IsSupported && size >= 16)\n    {\n        var vMin = Vector256.Create(double.PositiveInfinity);\n        var vMax = Vector256.Create(double.NegativeInfinity);\n        var vM1 = Vector256<double>.Zero;\n        var vM2 = Vector256<double>.Zero;\n        var vM3 = Vector256<double>.Zero;\n        var vM4 = Vector256<double>.Zero;\n        var v3 = Vector256.Create(3.0);\n        var v4 = Vector256.Create(4.0);\n        var v6 = Vector256.Create(6.0);\n        long c = 0;\n        for (int top = size & CommonMatrix.AVX_MASK;\n             i < top; i += 4)\n        {\n            c++;\n            var vSample = Avx.LoadVector256(samples + i);\n            vMin = Avx.Min(vMin, vSample);\n            vMax = Avx.Max(vMax, vSample);\n            var vd = Avx.Subtract(vSample, vM1);\n            var vs = Avx.Divide(vd,\n                Vector256.Create((double)c));\n            var vt = Avx.Multiply(Avx.Multiply(vd, vs),\n                Vector256.Create((double)(c - 1)));\n            vM1 = Avx.Add(vM1, vs);\n            var t1 = Avx.Multiply(Avx.Multiply(vt, vs),\n                Vector256.Create((double)(c * (c - 3) + 3)));\n            var t2 = Avx.Multiply(Avx.Multiply(vs, vM2), v6);\n            var t3 = Avx.Multiply(v4, vM3);\n            vM4 = vM4.MultiplyAdd(Avx.Subtract(\n                Avx.Add(t1, t2), t3), vs);\n            t1 = Avx.Multiply(vt,\n                Vector256.Create((double)(c - 2)));\n            t2 = Avx.Multiply(vM2, v3);\n            vM3 = vM3.MultiplyAdd(Avx.Subtract(t1, t2), vs);\n            vM2 = Avx.Add(vM2, vt);\n        }\n        var acc01 = Mix(c,\n            vM1.ToScalar(), vM2.ToScalar(),\n            vM3.ToScalar(), vM4.ToScalar(),\n            vM1.GetElement(1), vM2.GetElement(1),\n            vM3.GetElement(1), vM4.GetElement(1));\n        var acc23 = Mix(c,\n            vM1.GetElement(2), vM2.GetElement(2),\n            vM3.GetElement(2), vM4.GetElement(2),\n            vM1.GetElement(3), vM2.GetElement(3),\n            vM3.GetElement(3), vM4.GetElement(3));\n        var a = Mix(c + c,\n            acc01.m1, acc01.m2, acc01.m3, acc01.m4,\n            acc23.m1, acc23.m2, acc23.m3, acc23.m4);\n        if (Count == 0)\n            (Count, m1, m2, m3, m4)\n                = (4 * c, a.m1, a.m2, a.m3, a.m4);\n        else\n        {\n            long acCnt = 4 * c, n = Count + acCnt, n2 = n * n;\n            double d = a.m1 - m1, d2 = d * d;\n            double d3 = d2 * d, d4 = d2 * d2;\n\n            double nm1 = (Count * m1 + acCnt * a.m1) \/ n;\n            double nm2 = m2 + a.m2 + d2 * Count * acCnt \/ n;\n            double nm3 = m3 + a.m3\n                + d3 * Count * acCnt * (Count - acCnt) \/ n2\n                + 3 * d * (Count * a.m2 - acCnt * m2) \/ n;\n            m4 += a.m4 + d4 * Count * acCnt\n                    * (Count * (Count - acCnt) \n                        + acCnt * acCnt) \/ (n2 * n)\n                + 6 * d2 * (Count * Count * a.m2\n                    + acCnt * acCnt * m2) \/ n2\n                + 4 * d * (Count * a.m3 - acCnt * m3) \/ n;\n            (m1, m2, m3, Count) = (nm1, nm2, nm3, n);\n        }\n        min = Min(min, vMin.Min());\n        max = Max(max, vMax.Max());\n    }\n    for (; i < size; ++i)\n        Add(samples[i]);\n\n    static (double m1, double m2, double m3, double m4) Mix(\n        long c,\n        double a1, double a2, double a3, double a4,\n        double b1, double b2, double b3, double b4)\n    {\n        long n = c + c, n2 = n * n;\n        double d = b1 - a1, d2 = d * d, d4 = d2 * d2;\n        return (\n            (a1 + b1) \/ 2,\n            a2 + b2 + d2 * c \/ 2,\n            a3 + b3 + 3 * d * (b2 - a2) \/ 2,\n            a4 + b4 + d4 * c \/ 8 + 3 * d2 * (b2 + a2) \/ 2\n               + 2 * d * (b3 - a3));\n    }\n}\n<\/pre>\n
Observe que la precondici\u00f3n para aprovechar las instrucciones vectoriales es tener todo un array<\/em> de muestras a nuestra disposici\u00f3n. Si nos diesen un IEnumerable<double><\/code>, tendr\u00edamos que hacer maniobras como materializar las muestras en grupos de cuatro, en un array<\/em>, y alimentar as\u00ed al animalito vectorial.<\/p>\n
El c\u00f3digo es relativamente sencillo, si miramos con atenci\u00f3n. La parte AVX pr\u00e1cticamente repite el c\u00f3digo del Add<\/code> escalar. Por cada campo de Accumulator<\/code> hay un vector de doble precisi\u00f3n. La excepci\u00f3n es la propiedad Count<\/code>, y la tratamos diferente porque para los cuatro acumuladores virtuales que maneja el m\u00e9todo, la cantidad de muestras es siempre la misma.<\/p>\n
Esto es una ventaja cuando tenemos que mezclar los resultados de los cuatro acumuladores. La funci\u00f3n interna est\u00e1tica Mix<\/code> aprovecha la igualdad de los contadores para simplificar algebraicamente algunas f\u00f3rmula. Observe, por ejemplo, que la f\u00f3rmula para el m3<\/code> combinado es m\u00e1s sencilla, al anularse uno de los t\u00e9rminos.<\/p>\n
Una vez que hemos mezclado los cuatro acumuladores parciales, mezclamos el resultado, a su vez, con los valores que pueda haber ya en el propio acumulador (si los hubiera). Aqu\u00ed no podemos simplificar tanto, porque los contadores nuevos y antiguos pueden ser muy diferentes, aunque en el caso en el que el acumulador inicial no tuviese muestras, es todo m\u00e1s simple.<\/p>\n
Si quiere hacerse una idea de cu\u00e1nto mejora este tipo de procesamiento vectorial, los benchmarks<\/em> que he ejecutado me dan casi cinco veces m\u00e1s velocidad. Es extra\u00f1o, porque yo esperar\u00eda una mejora de 4x, pero puede deberse a que aqu\u00ed s\u00ed<\/em> hacemos uso de las instrucciones FMA vectoriales, cuando est\u00e1n disponibles. Las instrucciones FMA est\u00e1n escondidas en los m\u00e9todos de extensi\u00f3n MultiplyAdd<\/code> que present\u00e9 en esta entrada<\/a>.<\/p>\n
Por cierto, la ni\u00f1a de la imagen de la entrada tiene poco que ver con el algoritmo, pero estoy usando im\u00e1genes generadas por AI, entre otros motivos, para evitar problemas de derechos de autor. En este caso, le ped\u00ed a la AI que generase una ni\u00f1a perdida e indefensa en un universo digital simulado. En parte, la AI me hizo caso; en parte, ignor\u00f3 la petici\u00f3n. Pero el resultado me gusta, y ah\u00ed lo tiene.<\/p>\n","protected":false},"excerpt":{"rendered":"
En la entrada sobre la varianza, vimos que pod\u00edamos tener problemas de estabilidad num\u00e9rica si intent\u00e1bamos calcular la varianza en una sola pasada sobre los datos, usando inocentemente la definici\u00f3n matem\u00e1tica. La soluci\u00f3n de entonces fue usar un algoritmo de dos pasos: calcular la media en el primer paso, y en el segundo, calcular la […]<\/p>\n","protected":false},"author":1,"featured_media":1307,"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":[65,15,73,10,23,76,30,29],"class_list":["post-1301","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-c","tag-net","tag-algorithms","tag-austra","tag-optimization","tag-simd","tag-stability","tag-statistics","tag-variance"],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/intsight.com\/wp-content\/uploads\/2023\/08\/welford.png?fit=444%2C448&ssl=1","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1301","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=1301"}],"version-history":[{"count":30,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1301\/revisions"}],"predecessor-version":[{"id":1336,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1301\/revisions\/1336"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media\/1307"}],"wp:attachment":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media?parent=1301"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/categories?post=1301"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/tags?post=1301"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}