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*xx

He generalizado y metido un multiplicador λ, pero para una rotación podemos dejar que este λ 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 λ? 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*xx, 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 - λ*I) * 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 – λ*I valga cero.

Tomemos el caso más sencillo: una matriz de rotación en 2D, y veamos cómo queda el determinante.

characteristic polynomial

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:

Q = [v1, v2, v3 ... vn]

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 = [Av1, Av2, Av3 ... Avn]

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 vn son vectores propios, podemos simplificar la ecuación anterior de esta manera:

AQ = [λ1*v1, λ2*v2, λ3*v3 ... λn*vn]

Usando una de esas simplificaciones no muy evidentes, pero que se pueden comprobar fácilmente, la ecuación se puede reducir a esto:

AQ = 

La nueva matriz Λ 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Λ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:

AA = QΛQ-1QΛQ-1 = 2Q-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 μ 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 * Z + μ

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 Σ=A*AT 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 Σ 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:

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 Σ, 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 xT*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=A*AT. Como ejemplo sencillo:

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 * Z + μ, 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:
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<double>, 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<T> 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.