Categorías
FinTech

La distribución normal multivariante

La distribución normal multivariante es la generalización más inmediata de la distribución normal a un espacio multidimensional. Esto es: cada vez que tiremos los dados, queremos obtener, en vez de un número flotante, un vector de $N$ dimensiones.

La manera más sencilla de definir, y a la vez explicar, esta distribución es constructivamente. Primero tenemos que definir a qué llamaremos un «vector aleatorio normal estándar». Esto es simplemente un vector cuyos elementos son variables aleatorias normales independientes, cada una con media cero y varianza uno… como las que genera nuestro iterador BoxMuller de la entrada anterior.

Ahora supongamos que $Z$ es uno de estos vectores aleatorios normales y estándares, que $A$ es una matriz de dimensiones compatibles con $Z$, y que $\mu$ es un vector que, para simplificar, asumiremos que tiene las mismas dimensiones que $Z$. Entonces, los vectores aleatorios $X$ definidos mediante la siguiente ecuación pertenecen a una distribución normal multivariante:
$$
X = A \times Z + \mu
$$Para nosotros, los programadores, esto simplemente quiere decir que podemos generar vectores aleatorios normales multivariantes generando primero vectores gaussianos independientes y luego transformándolos con una multiplicación matricial seguida de una suma vectorial.

Intuitivamente, es más o menos claro que la suma vectorial nos sirve para mover la esperanza de la distribución, pero no es tan sencillo ver para qué multiplicamos por una matriz. La respuesta es que así conseguimos que las distintas dimensiones de la distribución no sean independientes. La matriz $\Sigma = A \times A^T$ sería entonces la matriz de covarianza entre las dimensiones de la distribución.

Una distribución muy general

La definición constructiva anterior es muy general, con toda intención. De hecho, en la definición más general, los vectores $X$ y $Z$ no tienen necesariamente que tener la misma dimensión, y la matriz $A$ puede ser, en consecuencia, una matriz rectangular.

De hecho, nuestra definición no garantiza que $\Sigma$ sea una matriz de covarianza razonable. Para ello, todos sus elementos tendrían que ser no negativos, y los elementos de la diagonal, en particular, tendrían que ser positivos. Eso no se cumple para cualquier $A$, y cuando no se cumple, no se puede definir una función de densidad para la distribución. Pero cuando la matriz de covarianza está bien definida, ocurre algo interesante, porque la función de densidad asociada se puede escribir de esta manera:
$$
{1 \over \sqrt{(2\pi)^k\vert\Sigma\vert}}e^{-{1\over 2}(x – \mu)^T \Sigma ^{-1}(x – \mu)}
$$Esta definición es casi idéntica a la de una gaussiana escalar. Las diferencias son que utilizamos vectores para el argumento y la media, y que en vez de tener la varianza en el denominador de la exponencial, utilizamos la inversa de la matriz de covarianza (la variable misteriosa k del factor de escala es simplemente el número de dimensiones de la distribución).

Monsieur Cholesky

¿Y si partimos del extremo contrario? En vez de plantearnos la distribución más general posible, teóricamente, podemos partir de una función de densidad ya asumida. Esto es: tenemos una distribución multivariante, y ya conocemos (o podemos calcular) su media y su matriz de covarianza. Tenemos la matriz $\Sigma$, y lo que queremos es encontrar qué matriz $A$ multiplicada por su traspuesta genera la matriz de covarianza…

Permettez-moi de vous présenter M. Cholesky. André-Louis Cholesky fue un militar y matemático francés, muerto en combate pocos meses antes de que terminase la Primera Guerra Mundial. Durante el conflicto, se dedicó a la geodesia y, para facilitar la confección de mapas, inventó eso que ahora llamamos «descomposición matricial de Cholesky», y que podemos entender intuitivamente como una forma de calcular la raíz cuadrada de una matriz.

La descomposición puede aplicarse a matrices hermitianas definidas positivas; si sabemos que la matriz sólo contiene valores reales, esto es equivalente a pedir que la matriz sea simétrica y que la expresión $x^T M x$ sea estrictamente positiva para cualquier vector no nulo. Y, vaya, esto lo cumple cualquier matriz de covarianza decente. Con esta premisa, se cumple entonces que existe una matriz triangular inferior $L$ tal que $M=L \times L^T$. Como ejemplo sencillo:
$$
\pmatrix{1&0.5\cr 0.5&1} = \pmatrix{1&0\cr 0.5&0.866} \times \pmatrix{1&0.5\cr 0&0.866}
$$No voy a describir en esta entrada el algoritmo para calcular la factorización (quizás más adelante), pero es un algoritmo sencillo, que ya implementan la casi totalidad de las librerías numéricas.

Con todos estos elementos en la mano, ya tenemos una receta para generar vectores aleatorios normales con dimensiones correlacionadas:

  1. Necesitamos conocer o calcular tanto la media como la matriz de covarianza de la distribución deseada.
  2. Calculamos la descomposición de Cholesky de la matriz de covarianza.
  3. Podemos entonces usar la fórmula $X=L\times Z + \mu$, donde $Z$ es un vector aleatorio normal estándar que podemos generar con un algoritmo sencillo como el de Box-Muller o el del zigurat.
Categorías
FinTech

La transformación de Box-Muller

En casi todos los fenómenos aleatorios, ya pertenezcan a la física, la genética o las finanzas, la distribución normal, o de Gauss-Laplace (la de la famosa curva de la campana) juega un papel importante. Sin embargo, .NET no ofrece de serie una clase, o un método, que genere valores aleatorios pertenecientes a esta distribución. Podemos utilizar una librería de terceros, por supuesto. Pero no está de más conocer alternativas, sobre todo para aplicaciones pequeñas o pruebas de concepto, en los que no merezca la pena usar algo más completo.

El problema a resolver es: teniendo como punto de partida un generador de números aleatorios que utilice una distribución uniforme, como la clase Random, ¿cómo podemos transformarlos para obtener la distribución normal? Lo primero es ponernos de acuerdo sobre los parámetros de la distribución normal que generaremos. Hay dos parámetros: la media y la varianza. Pero podemos ceñirnos a una distribución con media igual a cero y varianza igual a uno. Es fácil cambiar de parámetros desplazando y estirando los números que vamos a generar.

¿Cuál es el algoritmo adecuado para transformar una distribución uniforme en una normal? La respuesta es el llamado algoritmo del zigurat, que realiza un muestreo por regiones. El enlace anterior incluye código en C#. Pero existe un algoritmo mucho más sencillo, que se conoce como la transformación de Box-Muller. Esta transformación convierte dos valores aleatorios $u$ y $v$, pertenecientes a una distribución uniforme sobre el intervalo [0, 1], en otros dos valores aleatorios, a los que llamaremos $x$ e $y$, pertenecientes a una normal con media cero y varianza uno. Las fórmulas necesarias son estas:
$$
\eqalign{x&=\sqrt{-2 \ln u} \cos 2\pi v\cr
y&=\sqrt{-2 \ln u} \sin 2\pi v}
$$Existen métodos alternativos, como el de Marsaglia, que evitan las funciones trigonométricas, pero al precio de descartar algunas muestras. Antes de recomendar el método original de Box-Muller, he hecho la prueba en un Core i7-4770, y no he encontrado diferencias significativas entre ambos métodos:

  1. Probablemente, los procesadores más o menos modernos (el mío es un Intel Core de cuarta generación, que ya tiene su edad) penalicen más los saltos que las funciones trigonométricas.
  2. Además, la función Random de .NET utiliza internamente un algoritmo relativamente bueno, pero que tiene su propio coste.

La manera más sencilla de implementar un generador de números aleatorios con las características anteriores sea probablemente utilizar un iterador basado en un bucle infinito.

public static IEnumerable<double> BoxMuller()
{
    Random rnd = new Random();
    while (true)
    {
        double u = Math.Log(1 - rnd.NextDouble());
        double r = Math.Sqrt(-u - u);
        double v = 2 * Math.PI * rnd.NextDouble();
        yield return Math.Cos(v) * r;
        yield return Math.Sin(v) * r;
    }
}

Por supuesto, esta es la implementación más tonta posible: la instancia que contiene las variables de estado de la iteración pertenece a una clase y ocupa memoria dinámica. Además, es bastante probable que el compilador llame a la propiedad Current y al método MoveNext a través del tipo de interfaz IEnumerator, con lo que se trataría de llamadas virtuales. Pero existen técnicas sencillas para resolver estos dos problemas, aunque las explicaré en otro momento. Si tiene prisa, puede mirar como la clase List implementa internamente su iterador (se utiliza una estructura). He hecho la prueba y, al menos en .NET Core, la ganancia en velocidad no es significativa.

La imagen de la entrada, por cierto, es una representación ficticia de la famosa torre de Babel. Quizás habría sido más apropiado usar una imagen de un zigurat, pero pensándolo mejor, la forma de la torre se parece un poco a la campana de Gauss.

Categorías
C#

Un ejemplo de estructura

Para ver algo más de las nuevas posibilidades de las estructuras de C#, voy a mostrar una clase sencilla que utilizo en Iridium, un motor de valoración de swaps desarrollado para Comunytek, en su versión nativa para .NET Core. La estructura (¡qué sorpresa!) sirve para representar fechas pero, a diferencia de DateTime, sin la parte de la hora.

Internamente, una fecha se representa como el número de días transcurridos desde el uno de enero del año 1. Como son fechas para software financiero, me la trae al viento los problemas de cambio de calendario. Con estas premisas, puedo representar el número de días con un valor entero de 32 bits, con lo que mi tipo Date ocupa la mitad del espacio que un DateTime. Esta es la declaración de la estructura y de su único campo de estado:

/// <summary>A date with efficient operations.</summary>
public readonly struct Date : IEquatable<Date>, IComparable<Date>
{
    /// <summary>Number of days since Jan 1st, 1.</summary>
    private readonly int date;
    public Date(int year, int month, int day) { ... }
}

Para ganar velocidad, además, el constructor no verifica que los componentes de una fecha sean correctos: en Iridium, eso es responsabilidad del generador de cupones y otras partes del código que generan fechas.

Ahora viene la parte más interesante del tipo de datos:

public void Deconstruct(
     out int year, out int month, out int day) { ... }

Un deconstructor es un método que permite, precisamente, extraer en una sola operación las partes integrantes de una instancia de un tipo. Se introdujeron pensando en las tuplas, pero en realidad se pueden usar con cualquier clase o estructura, y es una pena que DateTime no cuente con uno de estos métodos. Gracias al deconstructor, podemos ejecutar instrucciones como la siguiente:

var (y1, m1, day1) = fromDate;

Si tuviese que usar DateTime, tendría que llamar por separado a las tres propiedades Year, Month y Day de la fecha. ¿El problema? Pues que para recuperar el año a partir de la representación interna hay que ejecutar un pequeño algoritmo que lleva su tiempo. Si luego quiero el mes, no importa: tengo que volver a ejecutar la parte que extrae el año. Y lo mismo pasa al pedir el día del mes. Con el deconstructor, en cambio, sólo tengo que descomponer la fecha en partes una sola vez, y obtengo los tres componentes. De hecho, la implementación de mi propiedad Day se permite el lujo de usar internamente el deconstructor:

public int Day
{
    get
    {
        Deconstruct(out _, out _, out int d);
        return d;
    }
}

Los subrayados son comodines para descartar el año y el mes.

¿Más cosas que pueden interesar? Por ejemplo, se permiten las conversiones de tipo entre Date y DateTime, en ambos sentidos, y permito convertir un valor Date en un número (la conversión inversa se consigue más elegantemente con un constructor adicional):

public static explicit operator DateTime(Date d) =>
    new DateTime(d.date * TicksPerDay);
public static explicit operator Date(DateTime d) =>
    new Date((int)(d.Ticks / TicksPerDay));
public static explicit operator int(Date d) => d.date;

El código completo de la clase puede descargar desde este enlace.

Categorías
C#

Estructuras en C#

Hay un antes y un después en C# en lo que atañe a los tipos struct, y tiene que ver con la llegada de RyuJIT, más que con .NET Core. El compilador de código nativo ha aprendido a manejar eficientemente estos tipos de datos.

Todos sabemos que las estructuras se diferencian de las clases en que son tipos con semántica de asignación por valor, en vez de la semántica habitual de asignación por referencia. Una variable de tipo struct contiene directamente los datos, mientras que una variable de clase es realmente un puntero a la zona de memoria donde residen los datos.

Lo quiero explicar en esta entrada, sin embargo, es cuándo es eficiente y aconsejable utilizar un tipo de estructura en vez de un tipo de clase. Hay un par de casos, y el primero de ellos es el más evidente: cuando el espacio que ocupan los campos de la estructura es suficientemente pequeño como caber en un registro de hardware.

Por ejemplo, .NET representa las fechas mediante el tipo de estructura DateTime. Este tipo guarda internamente un campo de tipo long, que representa el número de tics transcurridos desde cierta fecha base. Esto cabe en un registro de 64 bits. Para que todo vaya sobre ruedas, además, la estructura se declara de sólo lectura:

public readonly struct DateTime { }

En realidad, .NET Core es quien declara el tipo de sólo lectura. .NET Framework no lo hace, pero trata el tipo como si lo fuese. El cualificador readonly es relativamente nuevo en C#.

Para entender la utilidad de este caso de uso (disfrazar un tipo de datos numérico) hay que comparar con lo que ocurriría en un lenguaje como Java, que de momento no permite definir tipos de valor. O metes el entero en una clase, y cada vez que necesitas una fecha necesitas pedir memoria, o defines un montón de métodos que trabajen directamente con un entero… y te buscas la vida luego para averiguar cuándo un entero es un número de verdad o está representando una fecha.

Por supuesto, hay otra posibilidad: permitir modificaciones sobre las instancias de tu clase Java. Pero así es menos elegante y te arriesgas a todo tipo de errores. Y, de todas maneras, si necesitas crear 10.000 fechas, tienes otros tantos objetos ocupando memoria dinámica.

El segundo caso

Es más complicado explicar el segundo caso, por lo que voy a recurrir a un ejemplo en uno de mis proyectos: necesitaba un tipo para representar matrices. Las operaciones sobre ellas no formaban parte de la ruta crítica para la velocidad del código, pero la transformación de un vector por una matriz sí era crítica. Además, iba a necesitar muchas instancias de matrices.

Mi solución fue declarar la matriz como una estructura de sólo lectura. No se trata de un tipo pequeño, pues cada una necesita nueve números flotantes de doble precisión. Pero de esta manera ahorré memoria y, muy importante, una indirección innecesaria cada vez que tenía que transformar un vector.

El hacer que el tipo sea de sólo lectura es importante en estos casos, sobre todo si quiere declarar operadores. Mi multiplicación de matrices, por ejemplo, se implementa mediante un operador:

public static Matrix operator *(in Matrix m1, in Matrix m2){}

Si no utilizo los modificadores in en los parámetros, las matrices se pasan por copia. Pero si utilizo in y el tipo no es de sólo lectura, me sigo arriesgando a que el compilador haga una copia preventiva de la matriz, porque al fin y al cabo, in sólo significa que el operador no va a modificar el parámetro en su implementación. Es un poco raro, pero es lo que hay.

Hay que tener en cuenta, además, que en un método normal puedo pasar parámetros mediante ref, pero esto no es posible para un operador.

Un caso más

En mi código, por motivos históricos, hay un tercer caso intermedio: el de los vectores. Mis vectores no son de sólo lectura: quizás sería recomendable, pero me da pereza rescribir todo el código y hacer pruebas para comprobar que todo va mejor. Sin embargo, hago uso profuso del modificador in con vectores pasados como parámetros. ¿Cómo evito el desastre de la copia preventiva? He aquí la respuesta:

public readonly Vector Scale(in Vector factor) =>
    new Vector(X * factor.X, Y * factor.Y, Z * factor.Z);

C# permite declarar métodos con el modificador readonly, para garantizar que el método no modifica campos de la estructura. De este modo, si un método que recibe un vector como parámetro in ejecuta el método Scale u otro método de sólo lectura, el compilador evita crear una copia defensiva del vector.

Más novedades

Hay muchas más novedades en las últimas versiones de C# que tienen que ver con las estructuras (ref struct y ref readonly, por ejemplo), pero en esta entrada me he limitado a tratar algunas de las características que ayudan a la mejor generación de código nativo.

Categorías
Insights

Prelude

Este va a ser un blog de Informática, pero no sólo de ella. Voy a utilizarlo también para anotar ideas relacionadas: matemáticas, finanzas, música, física, y lo que se me ocurra.

¿Por qué empiezo un nuevo blog, en vez de retomar alguno de los antiguos? Pues porque mi trabajo ha cambiado mucho en los últimos años. Sigo programando, principalmente en C# y Java, pero cada vez tengo que ocuparme más de la parte matemática de mis proyectos. La gota que, en el buen sentido, ha colmado el vaso, es cierta proof of concept de algoritmos financieros para ordenadores cuánticos, en la que estoy participando. Sorprendente, pero cierto.

Así que… welcome back, my friends, to the show that never ends.