Descubriendo la estabilidad numérica

Recientemente tuve necesidad de calcular la desviación estándar de una serie de valores y tropecé con una serie de problemas interesantes, a pesar de lo sencillo que pareció a primera vista.

Comenzando por el principio:

La desviación estándar es una medida que se utiliza para cuantificar la variación o la dispersión de un conjunto de datos numéricos. Mide el grado de dispersión de los datos con respecto a la media. Mientras más pequeña es, son más homogéneos, el incremento de los valores de la desviación estándar indica una mayor variabilidad. Utiliza la misma unidad de medida que los datos.

Estrechamente relacionada con la desviación estándar se encuentra la varianza, que es otra medida de dispersión de los datos. La desviación estándar es la raíz cuadrada de la varianza. Comparten las mismas fórmulas incluyendo la raíz cuadrada para la desviación estándar.

La desviación estándar y la varianza se calculan de forma diferente cuando se quiere evaluar una población completa o solo una muestra, en este caso solo utilizaremos las fórmulas para una muestra de población. Por otra parte, diversas fuentes bibliográficas utilizan fórmulas aparentemente diferentes, pero no son más que transformaciones algebraicamente equivalentes, aunque pueden causar confusión al verlas por primera vez.

La fórmula para calcular la varianza es:

Donde X es la media, que se calcula como:

Estas fórmulas son sencillas, no suponen ninguna dificultad de implementación, unas sumatorias, unas multiplicaciones, unas divisiones, una raíz cuadrada, un estudiante de primer año podría hacerlo, ¿qué podría salir mal? Alerta de revelación (spoiler): todo.

El primer problema es que el algoritmo basado en las fórmulas anteriores, conocido como algoritmo de dos pasadas requiere realizar dos recorridos a los datos, lo que implica mantenerlos en memoria. Esto puede ser indeseable en muchas aplicaciones, por ejemplo, cuando el juego de datos es muy grande para ser almacenado en memoria o cuando es necesario calcular la varianza dinamicamente a medida que los datos son recolectados.

Afortunadamente, manipulando las fórmulas podemos llegar a la siguiente:

Esta forma es generalmente recomendada en los libros de texto como algoritmo de una pasada. Solo requiere mantener varios valores con los cuales se puede calcular la varianza o desviación estándar en cualquier momento. Podemos comenzar a implementar.

Al realizar pruebas para comprobar el funcionamiento del algoritmo, los resultados no se acercaban a los esperados. Los casos de prueba fueron generados usando las funciones de LibreOffice Calc, por lo que sería un poco fantasioso creer que el equipo de Calc se equivocaba y mi grandiosa implementación era correcta (probablemente alguna deidad algorítmica se arrastró de la risa ante ese pensamiento).

Luego de días de lucha infructuosa, descubro el artículo “Algoritmos para calcular la varianza de una muestra: análisis y recomendaciones”, escrito por Tony F. Chan, Gene H. Golub y Randall J. LeVeque para la revista “American Statistician”, en el año 1983. Realmente es difícil decidir qué fue más triste, descubrir que mi implementación sufría de una inestabilidad numérica de risa o que me lo dijera un artículo de 1983, escrito en TeX, dos años antes de que Tex fuera finalizado, ni siquiera era LaTeX.

En el subcampo matemático del análisis numérico, la estabilidad numérica es una propiedad de los algoritmos numéricos. Describe cómo los errores en los datos de entrada se propagan a través del algoritmo. En un método estable, los errores debidos a las aproximaciones se atenúan a medida que la computación procede. En un método inestable, cualquier error en el procesamiento se magnifica conforme el cálculo procede. Métodos inestables generan rápidamente anomalías y son inútiles para el procesamiento numérico.

Algunas veces un solo cálculo puede ser logrado de varias maneras, que pueden ser algebraicamente idénticas en términos de números reales o complejos, pero que en la práctica producen resultados diferentes según varían los niveles de estabilidad numérica.

Resulta que mis sumas y multiplicaciones tenían la mala costumbre de hacerse muy grandes muy rápido y mis restas y divisiones muy pequeñas, causando todo tipo de problemas, desde desbordamiento aritmético (que ocurre cuando el valor es mayor de lo que puede contener el tipo) hasta errores de redondeo y cancelaciones catastróficas, resultando en pérdida de precisión y un margen de error inaceptable.

Llegado a este punto podía dedicar (un largo) tiempo y esfuerzo a implementar un algoritmo numéricamente estable, con un correcto tratamiento de tipos para evitar pérdida de precisión y desbordamientos aritméticos o podía utilizar el arma secreta de los programadores. Legada por una larga tradición de desarrolladores e investigadores de todo el mundo, esta herramienta casi mágica ha sido vital en el desarrollo de esta ciencia: reutilización de código. Escogí la segunda.

Reutilizar código es un fino arte, como sabrán todos aquellos seguidores de Stackoverflow, saber qué copiar o utilizar y donde ponerlo, no es una tarea trivial. En este caso no fue Stackoverflow la fuente de conocimiento, sino otra forma de reutilización de código más antigua, las bibliotecas de código. Uno de los proyectos de la Apache Software Foundation es una biblioteca matemática y estadística que implementa todo lo necesario para resolver el problema de forma numéricamente estable: Commons Math. A partir de aquí, todo marchó sobre ruedas.

Las moralejas del día:

• La reutilización de código es una herramienta poderosa que debemos usar sabiamente para no reinventar la rueda.
• Se debe ser muy cuidadoso al realizar operaciones aritméticas con la computadora, esta es un dispositivo digital, capaz de realizar cálculos muy complejos pero en última instancia operan con números binarios y operaciones matemáticas simples.

Nota: Esta historia está escrita de forma lineal y sencilla, la realidad fue mucho más siniestra.