Categoría: FinTech

Categorías
C# FinTech

Secuencias

Esto es una perogrullada: en un lenguaje de programación funcional, las cosas se suelen hacer de forma diferente. Y probablemente escribir «perogrullada» sea una pedantería. Dejo el párrafo vivo por pereza, pero no me lo tenga en cuenta.

Factorial

Austra no pretende ser un lenguaje de programación funcional. Ni siquiera pretende ser «Turing-completo». En algún punto del camino, he pensado en diseñar un lenguaje «de verdad», e incluso en saltarme la maquinaria de .NET para compilarlo a código nativo directamente. Pero falta un poco para eso.

De todas maneras, hay que intentar que se puedan hacer todas las cosas posibles en Austra. Por ejemplo, ¿cómo calculas un factorial, si no tienes todavía funciones recursivas, y no quieres introducir bucles? Hasta ahora, la solución era parecida a ésta:

vec(10, i => i + 1).prod

No estamos definiendo una función, sino que estamos usando un parámetro a dedo, pero funciona. Creamos un vector de 10 entradas, lo llenamos con los números del 1 al 10, y multiplicamos todos los elementos del vector. Incluso podemos «optimizar» un poco la expresión:

vec(8, i => i + 2).prod * 2

Es absurdo dejar el 1 en el vector, y ya que tenemos diez elementos, en vez de dejarlos nueve, quito también el 2. Eso lo hago porque sé que internamente el producto mete los ocho elementos en dos registros AVX y no tengo que manejar el elemento nuevo aparte. El problema es que, de todas maneras, estamos pidiendo memoria para el vector.

En Austra 2.0 (que ya tiene soporte para .NET 8 y AVX512, dicho sea de paso), ya están implementadas las secuencias de valores reales, y puedes hacer esto otro, que es más natural:

seq(2, 10).prod

La nueva clase se llama seq, y es una copia descarada del diseño de enumerables y LINQ (o de los streams de Java, o de las «list comprehensions» de tantos otros lenguajes funcionales). Ahora tenemos un seq, pero luego vendrán iseq y cseq, para secuencias de enteros y complejos. Observe también que he acortado vector a vec, para usar cvec en vez de complexvector, y para que un futuro vector de enteros se pueda llamar ivec a secas.

Las secuencias pueden hacer casi todas las cosas que hacen los vectores, y en la mayoría de los casos, generarlas es más sencillo: con el vector, usamos una función lambda. También permiten ahorrar código porque aplican los mismos trucos que LINQ para objetos. Por ejemplo:

-- Esta es una secuencia que divide el intervalo
-- de 0 a 2*pi en 360 trozos.
-- Multiplicamos cada número por dos...
seq(0, τ, 360) * 2
-- .. pero internamente, lo transformamos en esto:
seq(0, 2τ, 360)

Recursividad eficiente

Y esto último todavía no está implementado, pero es bueno pensar en estas cosas antes de lanzarse a programarlas. ¿Cómo definiríamos una función para el factorial que fuese medianamente eficiente, pero usando recursividad? Lo que se suele hacer en un lenguaje funcional, pero con una notación al estilo Austra:

def fact(n: int) =
  let f = (n, acc: int) =>
    if n = 1 then acc else f(n - 1, acc * n) in
      f(n, 1)

El truco está en usar una función auxiliar, que al estilo Austra sería una lambda en una cláusula let. La función auxiliar es recursiva por la cola, por lo que el propio JIT de .NET la puede convertir en un bucle. La recursividad en lambdas va a ser problemática, porque estoy usando directamente Linq.Expressions para generar código, pero hay un truco sencillo, que es el que voy a usar en las primeras pruebas: declarar una variable de tipo lambda, asignarle un puntero vacío, y a continuación, asignarle una expresión lambda que ya podrá usar la variable recursivamente. La alternativa es usar un combinador Y. Es un truco que viene, precisamente, del cálculo lambda, sobre el que escribiré en algún momento. No es que sea cómodo o eficiente, pero es interesante.

De todas maneras, me he dado cuenta, revisando la documentación, que en F# hay que anunciar cuando una función va a ser recursiva. Una alternativa que tengo que pensar es generar un «método dinámico» en estos casos, con el esfuerzo adicional de generar directamente el código IL. Pero tengo ya experiencia de generar IL con Freya, y puedo reutilizar código.

Categorías
C# FinTech

Ostara

Todavía es work-in-progress, pero ya tenemos una versión open-source de una aplicación de escritorio, en WPF, para Austra. El código está ya en GitHub.

El editor de código es AvalonEdit. Reconoce la sintaxis del lenguaje de fórmulas e implementa una versión bastante decente de completamiento de código. Esto es lo que aporta la versión WPF respecto a la aplicación de consola: es mejor para aprender a usar el lenguaje. Los gráficos están hechos con OxyPlot, de momento.

Faltan cosas, tanto en la librería como en la aplicación, para darme por satisfecho, pero todo es cuestión de tiempo. Me gustaría, sobre todo, terminar de definir un mecanismo genérico de carga de series desde fuentes de datos externas.

Quiero también actualizar el código de la función de autocorrelación, e implementarla usando la transformada de Fourier que ya viene incluida. E incluir, de una puñetera vez, el modelo ARIMA completo y algo de GARCH y familia. De momento, sólo está incluido el modelo de series autoregresivas. También está pendiente una implementación con AVX de generación de números aleatorios.

Nota: autocorrelación ya actualizada. Como hay que añadir ceros al final de la serie para evitar que se cuelen las correlaciones cíclicas, la FFT se aplica a un número de muestras que es potencia de dos, y se utiliza el sub-algoritmo más eficiente.

Categorías
FinTech

Austra

Ostara, the Spring Goddess

Otro pequeño interludio: si os interesan estas cosas, echad un vistazo a la ayuda online de Austra. Es un lenguaje de expresiones orientado al manejo de series financieras y modelos econométricos. Está todavía en desarrollo, pero ya va tomando forma. Entre otras cosas interesantes, permite usar funciones lambda, y la librería, Austra.Library, utiliza AVX2 para optimizar las operaciones sobre matrices, vectores y series temporales. Hasta 200/300 dimensiones, es generalmente más eficiente que las alternativas, incluso las que usan Intel MKL. A partir de ese punto, naturalmente, Intel MKL es más eficiente. En cualquier caso, tengo la posibilidad de mejorar la implementación para ese volumen de información. Y como última alternativa, puedo incorporar Intel MKL como una opción más; opcional, por supuesto, porque esta librería tiene licencias de pago.

Mi idea es, en el momento en que esté suficientemente maduro, convertir parte del código en open source. Aparte de añadir más algoritmos a la librería básica, quiero dar alternativas al analizador del lenguaje. Ahora mismo, el lenguaje genera árboles de expresiones, que a su vez generan código nativo sobre la marcha. Para una línea de comandos, sin embargo, puede que esto sea excesivo. No voy a destruir el compilador actual, pero quiero subir, antes de abrir el código, un intérprete más modesto que, al mismo tiempo, seguramente será más eficiente para el propósito de la herramienta.

Categorías
FinTech

Covarianza

Recordemos cómo se define la varianza de una variable aleatoria X:
$$
Var(X) = \mathbb{E}[(X – \mathbb{E}[X])^2]
$$Vamos a suponer que tenemos, en vez de una, dos variables aleatorias, $X$ e $Y$. Si manipulamos un poco la definición de varianza, obtendremos la definición de la covarianza entre dos variables:
$$
Cov(X, Y) = \mathbb{E}[(X – \mathbb{E}[X])*(Y – \mathbb{E}[Y])]
$$Según esta definición, la varianza de una variable es la covarianza de la variable consigo misma, por lo que la definición parece tener sentido:
$$
Var(X) = Cov(X, X)
$$Además, la covarianza es simétrica respecto a los argumentos:
$$Cov(X,Y)=Cov(Y,X)
$$La interpretación de la covarianza no es del todo inmediata. En este caso sencillo, si la covarianza es positiva, cuando la $X$ aumenta, también aumenta la $Y$: las variaciones van coordinadas en la misma dirección. Si la covarianza es negativa, cuando $X$ aumenta, la $Y$ disminuye. Si la covarianza es cero, son variables independientes. Sin embargo, para medir el grado de relación entre dos variables aleatorias, existe una medida mejor, llamada precisamente correlación. Dedicaré, en otro momento, una entrada a la correlación. De momento, puedo adelantar que la correlación es una forma de covarianza «normalizada» para que sus valores vayan siempre entre menos uno y uno.

Matrices de covarianza

Los valores de las varianzas y la covarianza de dos variables se pueden organizar en una tabla o matriz:
$$\pmatrix{Cov(X,X)&Cov(X,Y)\cr Cov(Y,X)&Cov(Y,Y)}
$$
Esto nos interesa porque el siguiente paso es extender la definición de covarianza a cualquier número de variables aleatorias. Por ejemplo, con tres variables aleatorias:
$$
\pmatrix{c_{1,1}&c_{1,2}&c_{1,3} \cr c_{2,1}&c_{2,2}&c_{2,3} \cr c_{3,1}&c_{3,2}&c_{2,3}}
$$Como es fácil de entender, una matriz de covarianza es una matriz simétrica, por definición.

Matrices semidefinidas positivas

Una matriz de covarianza tiene siempre la interesante propiedad de ser semidefinida positiva. Para todo vector x se debe cumplir lo siguiente:
$$
\forall x : x^T \cdot \Sigma \cdot x \ge 0
$$Vamos a interpretar geométricamente la propiedad en cuestión. Lo que estamos haciendo es transformar el vector x con la matriz de covarianza. Luego calculamos el producto escalar del vector respecto al vector transformado. Si decimos que ese producto escalar es mayor o igual a cero, estamos diciendo que, sin importar el número de dimensiones del vector y la matriz, el ángulo entre el vector y el vector transformado está siempre entre -π/2 y +π/2. En otras palabras, la matriz nunca «retuerce» demasiado los vectores que transforma.

¿Nos sirve de algo esta imagen visual? Pues no lo sé. Pero me gusta tener presente este tipo de interpretaciones gráficas. Por experiencia, terminas encontrándole un uso más tarde o más temprano.

Esto sí es importante: una matriz semidefinida positiva tiene, obligatoriamente, un determinante mayor o igual que cero. Este criterio puede servir para descartar rápidamente matrices de covarianza mal construidas. De hecho, si el determinante es cero, es porque existen variables aleatorias redundantes.

Categorías
C# FinTech

Cholesky

Esta va a ser una entrada «utilitaria»: poca explicación, pero aportando código fuente, por si alguien lo necesita. ¿Recuerda la entrada en la que presenté la descomposición de Cholesky? En aquel momento, no incluí el algoritmo de la descomposición, porque quería experimentar un poco con la implementación. Ya lo he hecho, y estoy más o menos satisfecho con el resultado.

Antes de ver el código, le explico el contexto. Este es un método de una estructura Matrix, que encapsula solamente una matriz bidimensional de valores flotantes (campo values). Dentro de la estructura se definen todos los operadores y métodos que podéis imaginar. En paralelo, defino otra estructura, LowerMatrix, para representar matrices triangulares inferiores. No hay relaciones de herencia, al tratarse de estructura, pero la clase de matrices triangulares permite definir métodos y operadores más eficientes. Es importante tener en cuenta que LowerMatrix gasta exactamente la misma memoria que Matrix. La ventaja está en esos métodos que evitan procesar las dos mitades de la matriz.

También hace falta saber que he definido operadores de conversión implícitos que transforman un array bidimensional en uno u otro tipo de matrices. Este operador es el que me permite evitar un constructor explícito para devolver un valor en el método:

public unsafe LowerMatrix Cholesky()
{
    int n = Rows;
    double[,] dest = new double[n, n];
    double[,] src = values;
    double* tmp = stackalloc double[n + n];

    // First column is special.
    double ajj = src[0, 0];
    if (ajj <= 0)
    {
        dest[0, 0] = double.NaN;
        return dest;
    }
    dest[0, 0] = ajj = Math.Sqrt(ajj);
    double r = 1 / ajj;
    n--;
    for (int i = 1; i <= n; i++)
        dest[i, 0] = src[i, 0] * r;
    for (int j = 1; j <= n; j++)
    {
        // Compute the diagonal cell.
        double v = 0.0;
        for (int i = 0; i < j; i++)
        {
            double a = dest[j, i];
            v += a * a;
        }
        ajj = src[j, j] - v;
        if (ajj <= 0)
        {
            dest[j, j] = double.NaN;
            return dest;
        }
        dest[j, j] = ajj = Math.Sqrt(ajj);

        // Compute the other cells of column J.
        if (j < n)
        {
            r = 1 / ajj;
            for (int i = 0; i < j; i++)
                tmp[i] = dest[j, i];
            for (int i = j; i < n; i++)
            {
                v = 0.0;
                for (int k = 0; k < j; k++)
                    v += dest[i + 1, k] * tmp[k];
                tmp[i] = v;
            }
            for (int i = j; i < n; i++)
                dest[i + 1, j] = (src[i + 1, j] - tmp[i]) * r;
        }
    }
    return dest;
}

Sobre el algoritmo, en sí: el algoritmo de Cholesky puede fallar cuando la matriz de origen no es positiva semidefinida. Mi implementación detecta ese caso al calcular las raíces cuadradas… y se limita a parar, poniendo un NaN en la celda diagonal donde se ha detectado el problema. Esto significa que el método, en la práctica, asume que la matriz es positiva semidefinida. En mi código, tengo un segundo método, TryCholesky, que devuelve un valor lógico para ver si la conversión fue posible, y retorna la matriz transformada como parámetro de salida.

Desde el punto de vista de un programador de C#, el único detalle interesante es el uso de stackalloc para reservar un array de memoria en la pila, en vez de usar memoria dinámica. Esto es lo que obliga a declarar el método con unsafe.

En rendimiento, el método es más rápido que la «versión base» de la librería que he visto que es más rápida usando sólo C# (usando Intel MKL, no hay competencia posible). Me refiero a la versión base porque, para matrices grandes, las librerías serias suelen dividir la matriz en bloques que se pueden procesar en paralelo. Este código, como ve, no usa threads, instrucciones SIMD y sólo utiliza punteros para la caché. En menos palabras: todo es mejorable.

Categorías
FinTech

Varianza

Todo el mundo sabe lo que es una media (me refiero a la media estadística, por supuesto, no a las otras), por lo que no tiene mucho sentido escribir una entrada sobre medias. Conceptualmente, sin embargo, la media es el primer elemento de una serie de valores que caracterizan a las series aleatorias, y que se conocen como momentos. Si nos saltamos el primer momento, aterrizamos en el segundo momento estadístico, que es la varianza o, cuando tenemos una distribución multivariante, la matriz de covarianza.

Definición

Casi todo el mundo tiene también claro cómo se define la varianza, pero refresquémoslo, por si acaso. Supongamos que la media de una serie, o de una variable aleatoria, se representa como $\mathbb{E}[X]$. Entonces, la varianza se puede definir así:
$$
Var(X)=\mathbb{E}[(X-\mathbb{E}[X])^2]
$$Sí, estamos usando la media dentro de otra media. Según esta definición, literalmente, tenemos primero que conocer la media de la variable aleatoria. Entonces, tenemos que ver cuánto se desvía la variable aleatoria de esa media. Esto es, se mide el grado de «desparrame» de la variable. No podemos medir directamente la media de la resta que hemos mencionado, porque los valores negativos se compensarían con los valores positivos, y obtendríamos siempre cero. La manera de evitar esta compensación es elevar la resta al cuadrado.

Si masajeamos un poco la definición, podemos transformarla de esta manera:
$$
Var(X)=\mathbb{E}[X^2] – \mathbb{E}[X]^2
$$La transformación es fácil de realizar, y no la incluyo por brevedad. La nueva fórmula viene a decir que la varianza es la diferencia entre la media del cuadrado y el cuadrado de la media. Enseguida veremos la gran utilidad de esta definición alternativa.

Por cierto, a partir de la varianza se define la desviación estándar, que es simplemente la raíz cuadrada de la varianza. ¿Para qué necesitamos la raíz cuadrada? Principalmente, por las unidades de medida. Si la variable aleatoria X representa euros, la media estará expresada como euros, pero la varianza estará en euros al cuadrado. La desviación estándar, en cambio, vuelve a estar en euros. Probablemente tenga otros usos y beneficios, pero ahora mismo no los recuerdo.

Más definiciones

Si un matemático leyese lo que acabo de escribir, seguramente me pegaría una somanta de palos. Y con razón. He manejado con demasiada alegría los términos «muestra aleatoria» y «variable aleatoria», pero si hubiese sido más riguroso, la introducción habría sido infinita. Ahora aclaro lo que hace falta aclarar desde el punto de vista práctico.

La diferencia que tenemos que conocer es la de «varianza» contra «varianza de una muestra». La definición anterior es válida para una variable aleatoria definida analíticamente, o para fenómenos en los que disponemos de todos los datos. En la vida real, lo que solemos tener es una muestra de una serie, no todos los elementos. En ese caso, lo que podemos hacer pragmáticamente es calcular un «estimado» de la varianza a partir de la muestra.

El ajuste, de todas maneras, es sencillo:
$$
Var_M[X] = Var(X) \cdot N / (N – 1)
$$donde N es el tamaño de la muestra. Este es el motivo por el que Excel tiene dos funciones, VAR y VAR.P, para la varianza. La primera es la varianza de una muestra, y la segunda es la varianza «completa» de una población.

Implementación

Escribo esta entrada, además de para que sirva de referencia a entradas posteriores, porque quiero explicar un pequeño truco de implementación que he visto pasar por alto muchas veces.

Vamos a suponer que tenemos una muestra aleatoria en un array. Para calcular la media, hacemos un bucle y vamos sumando las entradas. Si utilizo la primera definición que hemos dado de la varianza, necesitamos hacer dos pasadas sobre el array. La primera, para calcular la media, y la segunda, para calcular la media de los cuadrados de las diferencias respecto a la media. ¿No sería mejor hacer una sola pasada? Si tenemos los datos en un array, no es tan acuciante, pero si los datos son grandes, o vienen de un enumerador, nos interesa hacerlo todo en una pasada.

Es en estos casos en los que la segunda fórmula equivalente es importante: podemos calcular simultáneamente, en una sola pasada, la media y la media de los cuadrados, y luego restar al segundo valor el cuadrado del primero. Algo así:

double sumX = 0, sumX2 = 0;
int total = 0;
foreach (double v in source)
{
    sumX += v;
    sumX2 += v * v;
    total++;
}
double mean = sumX / total;
double variance = (sumX2 / total - mean * mean);

Si quisiéramos la varianza de una muestra (es decir, si no tuviésemos todos los valores de la población) tendríamos que ajustar la varianza calculada multiplicándola por total y dividiéndola por total - 1.

De todos modos, el problema principal con el código anterior es otro, mucho más sutil. Las variables sumX y sumX2 pueden llegar a tener valores muy altos. Y al restar los dos valores se puede producir una cancelación catastrófica que nos haga perder la precisión del resultado. Ese problema no lo tiene el algoritmo en dos pasadas.

La solución, sin embargo, es muy sencilla. Resulta que:
$$
Var[X] = Var[X – a]
$$Es decir: si le restamos una constante arbitraria a la variable aleatoria, la «dispersión» o varianza de la misma no varía. Lo ideal sería que la constante en cuestión fuese la media de la muestra, pero para eso necesitaríamos dos pasadas. Lo que haremos es llegar a un compromiso: como no podemos tener la media, nos conformaremos con un valor representativo de la muestra. ¿Qué tal si elegimos el primero?

double sumX = 0, sumX2 = 0, mean = 0, first = 0;
bool hasFirst = false;
int total = 0;
foreach (double v in source)
{
    if (!hasFirst)
    {
        hasFirst = true;
        first = v;
    }
    else
    {
        double v1 = v - first;
        sumX += v1;
        sumX2 += v1 * v1;
    }
    mean += v; 
    total++;
}
mean /= total;
double variance = (sumX2 - sumX * sumX / total) / total;

He puesto un else dentro del bucle porque la primera suma no aporta nada a los acumuladores. Por esto, también he tenido que modificar el cálculo de la media al final. Idealmente, el primer elemento de la muestra debería estar lo más cercano posible a la media. Si sabemos que la media va a estar alrededor de cero, además, podemos olvidar esta precaución.

Quizás alguien crea que la instrucción condicional dentro del bucle lo puede ralentizar un poco. No lo he medido para ver si ocurre, pero si eso fuese importante, lo que podríamos hacer es desarrollar el bucle en dos trozos:

var enumerator = source.GetEnumerator();
if (enumerator.MoveNext()
{
    double first = enumerator.Current;
    double sumX = 0, sumX2 = 0, mean = first;
    int total = 1;
    while (enumerator.MoveNext())
    {
        double v = enumerator.Current - first;
        sumX += v;
        sumX2 += v * v;
        mean += enumerator.Current;
        total++;
    }
    mean /= total;
    double variance = (sumX2 - sumX * sumX / total) / total;
}

Me falta un else en el fragmento anterior, para cuando la secuencia está vacía, pero me da pereza ponerlo. También habría que llamar a Dispose().

Concurrencia

Los otros trucos interesantes, que no voy a tratar en esta entrada, tienen que ver con la posibilidad de repartir el cálculo entre varios hilos, cuando el tamaño de la serie lo amerite. No hay grandes complicaciones a la vista, pero hay que tener cuidado al mezclar los valores obtenidos en cada hilo. Si la serie está en un array, es quizás más sencillo, porque se pueden repartir a priori los rangos entre tareas. Lo recomendable, en cualquier caso, es utilizar una de las sobrecargas de Parallel.ForEach que permita el uso de variables de estados. Permítame que me haga el sueco y pase de página en este punto.

Categorías
FinTech

Valores y vectores propios

Though this be madness, yet there is method in’t.
Polonius

Esta va a ser, probablemente, la entrada más esotérica de esta serie. Yo mismo no tengo claro si esto me lo enseñaron en el primer año de la carrera. Supongo que sí, pero no tuve que usar estas cosas hasta mucho tiempo después.

Todo el mundo tiene una idea más o menos intuitiva sobre qué es un vector. Las intuiciones sobre las matrices no son tan populares, pero una que nos valdrá es considerar que una matriz representa una transformación sobre un vector. Esa transformación puede ser una rotación, un cambio de escala, una traslación (con ciertas modificaciones) o una combinación de todas estas cosas. Supongamos que, en un caso concreto, la transformación es una rotación. Entonces tiene que haber un eje de rotación, ¿no? Y ese eje de rotación va a estar determinado por un vector que debe cumplir la siguiente ecuación:
$$
A \times x = \lambda x
$$He generalizado y metido un multiplicador $\lambda$, pero para una rotación podemos dejar que este $\lambda$ valga 1. ¿Qué quiere decir entonces la ecuación anterior? Pues que existe un vector que se transforma en sí mismo. En realidad, cualquier múltiplo de ese vector se va a transformar en sí mismo. ¿Y para qué quiero entonces el multiplicador $\lambda$? Muy sencillo: imaginemos que la transformación es un cambio de escala uniforme en todas las direcciones. Cualquier vector cumple entonces la igualdad anterior de forma trivial.

Compliquemos un poco la transformación, entonces: vamos a estirar todos los vectores en una dirección. En este caso, si $A \times x = \lambda x$, entonces el vector x apunta en la dirección del estiramiento:

En general, esos vectores que se transforman en sí mismos, se conocen como «vectores propios» de la matriz o transformación, y los multiplicadores se conocen como «valores propios» de la transformación. En inglés: eigenvector y eigenvalue, respectivamente.

Cómo se calculan

¿Cómo se pueden calcular valores y vectores propios? Los algoritmos prácticos son relativamente complicados. Pero hay una forma relativamente sencilla cuando las matrices son pequeñas. Si manipulamos los términos de la definición, podemos agruparlos así:
$$
(A – \lambda I) \times x = 0
$$En este caso, $I$ es la matriz identidad (toda la diagonal con unos y el resto de las celdas con ceros). Por lo tanto, los valores propios, es decir, los valores que puede tomar λ son aquellos para los que el determinante de la matriz $A – \lambda I$ valga cero.

Tomemos el caso más sencillo: una matriz de rotación en 2D, y veamos cómo queda el determinante.
$$
\displaylines{\pmatrix{\cos \theta \,- \lambda & – \sin \theta \cr \sin \theta & \cos \theta \,- \lambda }
\cr
\lambda ^2 – 2 \lambda \cos ^2 \theta + 1 = 0}
$$El polinomio sobre λ que se genera se conoce como «polinomio característico» y, oops, en este caso tiene un discriminante negativo, lo que quiere decir que sus dos raíces son complejas (excepto si el seno del ángulo es igual a cero, en cuyo caso hay dos raíces reales idénticas). ¿No habíamos quedado en que el vector propio de una matriz de rotación era el eje de rotación? Sí, pero en dos dimensiones no existe un eje de rotación, porque quedaría siempre fuera del plano. Por este motivo, los valores propios son complejos y lo mismo ocurre con los vectores propios. Otras matrices 2D sí tienen valores propios reales, pero no las de rotación. Si, por el contrario, se tratase de una matriz de rotación en 3D, el polinomio característico sería de tercer grado. Y da la casualidad que todo polinomio de tercer grado (o de grado impar, en general) tiene al menos una raíz real. Si quieres, haz la prueba.

Power iteration

Para matrices pequeñas, los polinomios característicos son manejables, pero en estadísticas se suele trabajar con matrices enormes. Hay varios métodos «serios» para cubrir estos casos, pero todos son complicadillos de implementar. Es mejor tirar de librerías probadas que intentar reinventar la rueda. No obstante, existe un método muy sencillo que nos puede valer cuando sólo necesitamos un vector propio, el asociado al valor propio de más magnitud:

  1. Selecciona un vector aleatorio, y asegúrate de que su longitud sea uno.
  2. Multiplica el vector con la matriz.
  3. Normaliza el vector. Esto es, divídelo por su longitud para que el vector tenga nuevamente longitud uno.
  4. Repetir desde el paso dos, hasta que el vector converja a algo, o te aburras de esperar a la convergencia.

Este es el algoritmo conocido como «power iteration», y no siempre converge. Cuando lo hace, la velocidad de la convergencia depende de la magnitud de la diferencia entre el mayor de los valores propios y el siguiente. Tiene el defecto adicional de que sólo calcula ese valor propio. Para calcular los siguientes, hay que transformar la matriz.

¿Aplicaciones?

Las hay a montones… pero estoy siguiendo a rajatabla la táctica de hacer entradas pequeñas para evitar la tentación de abandonar el blog cuando tarde mucho en escribir cada entrada. Lo que sí puedo es adelantar algunos de los usos de estas cosas.

Por ejemplo, los «observables» en Mecánica Cuántica son valores propios de operadores hermitianos. Ojo: estoy hablando ahora de operadores en vez de matrices, pero la mayoría de estos operadores pueden representarse mediante matrices.

En Estadística, los vectores propios son la base de un algoritmo conocido como Principal Component Analysis, o PCA. Para ir haciendo boca, le adelanto una de las propiedades que personalmente me molan más. Imaginemos que formamos una matriz a partir de los vectores propios:
$$
A = [v_1, v_2, v_3 \cdots v_n]
$$Si todos los vectores propios son diferentes o, más bien, independientes, resulta que esta es una matriz ortogonal, que representa un giro en algún número de direcciones. Ahora transformaremos esta matriz con la matriz original:
$$
AQ = [Av_1, Av_2, Av_3 \cdots Av_n]
$$Si no ves inmediatamente lo que ocurre en el lado derecho, tranquilo, que es la falta de práctica: yo estas cosas las aprendí hace muchos años, y cuesta resucitarlas. Si es tu caso, aplica la fórmula de multiplicación de matrices y desarróllala. El caso es que, sabiendo que los $v_n$ son vectores propios, podemos simplificar la ecuación anterior de esta manera:
$$
AQ = [\lambda _1 v_1, \lambda _2 v_2, \lambda _3 v_3 \cdots \lambda _n v_n]
$$Usando una de esas simplificaciones no muy evidentes, pero que se pueden comprobar fácilmente, la ecuación se puede reducir a esto:
$$
AQ = Q \Lambda
$$La nueva matriz $\Lambda$ es simplemente una matriz diagonal con un valor propio en cada uno de los elementos de la diagonal. El último paso es multiplicar ambos lados de la igualdad por la inversa de $Q$, la matriz ortogonal:
$$
AQQ^{-1} = A = Q \Lambda Q^{-1}
$$En otras palabras, podemos descomponer la matriz original en una matriz ortogonal y una matriz diagonal, debidamente combinadas.

¿Qué tiene esto de interesante para que yo diga que me mola? Vamos a multiplicar la matriz $A$ por sí misma, es decir, vamos a elevarla al cuadrado:
$$
A \cdot A = Q \Lambda Q^{-1} \cdot Q \Lambda Q^{-1} = Q \Lambda ^2 Q^{-1}
$$Si ya tenemos la descomposición de la matriz, las potencias de la raíz se obtienen fácilmente elevando la matriz diagonal a la potencia deseada… que como se puede comprobar, es una operación muy sencilla.

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.