Categorías
C#

Reuniendo piezas de un vector

Cuando hay que utilizar instrucciones vectoriales, es siempre aconsejable que se cargue en memoria un trozo consecutivo de memoria para crear vectores. Para un programador ya entrenado en AVX/AVX2/AVX512, por ejemplo, es fácil darse cuenta de que el siguiente bucle (un producto escalar de dos vectores) es un candidato ideal para usar instrucciones SIMD:

double sum = 0.0;
for (int i = 0; i < size; i++)
    sum += p[i] * q[i];

Hay dos zonas de memoria, p y q, y el índice se mueve ascendentemente de uno en uno, en ambos casos. Por lo tanto, para convertir este código en instrucciones vectoriales de AVX, hacemos que cada paso del código cargue cuatro entradas de cada vector, en cada paso del bucle:

double sum = 0;
int i = 0;
if (Fma.IsSupported)
{
    var acc = Vector256<double>.Zero;
    for (; i < size - 4; i += 4)
        acc = Fma.MultiplyAdd(
            Avx.LoadVector256(p + i),
            Avx.LoadVector256(q + i),
            acc);
    acc = Avx.HorizontalAdd(acc, acc);
    sum = acc.ToScalar() + acc.GetElement(2);
}
for (; i < size; i++)
    sum += p[i] * q[i];

Gracias a este sencillo truco, se aprovecha mejor la caché de memoria del procesador. Es cierto que p y q pueden estar en direcciones distantes, pero al menos la carga de cada trozo de cuatro valores no provoca invalidaciones de la caché. En cambio, este otro fragmento de código, extraído del código de Austra que realiza la descomposición de una matriz según la factorización LU, no cumple con esta regla:

double m = c[k];
for (int i = k + 1; i < size; i++)
    c[i] -= m * a[i * size + k];

El almacenamiento en c es consecutivo, pero si reunimos cuatro pasos del bucle, las cargas desde el vector a van a venir de zonas dispersas de la memoria: esto es porque la variable de bucle en el indexador está multiplicada por size. Sería una pena no poder utilizar AVX aquí, porque podríamos reducir por cuatro las restas y multiplicaciones del bucle.

Por fortuna, Intel agregó al conjunto de instrucciones de AVX2 una instrucción que, precisamente, sirve para reunir en una misma carga fragmentos de memoria no consecutivos. Se trata de la instrucción vgatherdpd, que .NET implementa mediante el método Avx2.GatherVector:

public static Vector256<double> GatherVector256(
    double* baseAddress,
    Vector128<int> index,
    byte scale);

A primera vista, parece arameo antiguo. A segunda y tercera vistas, se asemeja más al sumerio de chabolas. Pero no es tan complicado explicar cómo funciona, sobre todo si lo acompañamos con un ejemplo. A este método le pasamos un puntero a un array de números de doble precisión, esa es la parte obvia. Además, le pasamos cuatro índices de tipo entero, en el parámetro de tipo Vector128<int>. Como este tipo contiene 128 bits, es evidente que ahí podemos pasar cuatro enteros de 32 bits; ya veremos cómo. ¿Y el factor de escala? Aquí es un poco redundante, y Microsoft podría habernos ahorrado el trabajo: cuando se trabaja con vectores de double, tenemos que pasar un 8, o lo que es igual, sizeof(double).

Juntando las piezas, he aquí como quedaría el código transformado del fragmento original:

double m = c[k];
int i = k + 1;
if (Avx2.IsSupported)
{
    var vm = Vector256.Create(m);
    var vx = Vector128.Create(0, size, 2 * size, 3 * size);
    for (double* p = a + i * size + k;
        i < size - 4;
        i += 4, p += size * 4)
        Avx.Store(c + i,
            Avx.Subtract(Avx.LoadVector256(c + i),
            Avx.Multiply(Avx2.GatherVector256(p, vx, 8), vm)));
}
for (; i < size; i++)
    c[i] -= m * a[i * size + k];

Expliquémonos: antes de entrar en el bucle, creamos un vector de números reales a partir del valor de la variable m. Lo que sigue es lo que nos interesa: creamos por adelantado, un Vector128<int> con cuatro offsets o desplazamientos diferentes, múltiplos del tamaño de la fila de la matriz, que viene en size. Y ahora viene el detalle artesanal delicado: al iniciar el bucle, creo una variable adicional de puntero, p, que apunte i * size + k entradas a partir de la dirección inicial, que está en a. Recuerde que leíamos el valor de a[i * size + k] en cada paso del bucle original, y que i es la variable de bucle. Si hubiese querido ser más claro, debería haber escrito:

double* p = a + (k + 1) * size + k

Pero me puede la pereza, y al fin y al cabo, en ese preciso instante, i es equivalente a k + 1, ¿no?

El detalle es que, en cada paso del bucle AVX, la variable p se incrementa en 4 * size entradas (es decir, en cuatro filas enteras de size entradas). Digo que es un detallazo elegante porque la alternativa obvia sería incrementar directamente los cuatro índices enteros que guardamos en el vector vx. ¿Por qué no lo he hecho así? Pues porque, al trabajar con vectores de 128 bits, la instrucción de suma pertenecería al conjunto de instrucciones SSE2… y resulta que el cambio entre AVX y SSE tiene un coste pequeño, pero no desdeñable.

… y este es el momento más apropiado para exclamar ¡pero coño! (con perdón), ¿esto no tendría que hacerlo el compilador, o el JIT, o quien narices sea?. Pues sí. De hecho, teóricamente, Hotspot, el JIT de Java, tendría que ocuparse de estas cosas. En la práctica, la vectorización automática en Hotspot es un truño como un puño. Y el famoso escape analysis tampoco hace lo que se supone que tendría que hacer. O si lo hace, no se nota. Al final, Java ha añadido una librería de vectores basados en instrucciones SIMD, al estilo de .NET. Los compiladores de C/C++ hacen mejor papel en estos casos, pero no son omnipotentes, y a veces tienes que echarles una mano.

NOTA: AVX2, creo recordar, se introdujo en la cuarta generación de Intel, Haswell, en 2013. La primera implementación en hardware de vgatherdpd no era para tirar cohetes, y todavía hay páginas y comentarios en Internet que lo recuerdan. La implementación mejoró notablemente en Skylake, la sexta generación (2015). Yo tuve un Haswell durante muchos años, y empecé a probar AVX2 con esa máquina. Mi penúltimo portátil fue un Skylake. Ahora mismo tengo un PC con CPU de oncena generación, Tiger Lake, que es un pepino. Lo menciono porque, por casualidad, Tiger Lake es la última generación de procesadores de «consumo» que incluye las instrucciones AVX512: con la llegada de las CPU con núcleos heterogéneos, estas instrucciones se deshabilitan porque los núcleos «eficientes» no las implementan. Ahora mismo, si quieres AVX512, tienes que irte a Xeon, o a AMD.

Categorías
C#

Multiplicar por la traspuesta

Una operación muy frecuente con matrices consiste en multiplicar una matriz por la traspuesta de otra. La implementación más sencilla de esto sería trasponer la segunda matriz y luego multiplicar la primera matriz por la traspuesta. Lamentablemente, trasponer una matriz es poco eficiente. No hay un patrón sencillo de acceso a la caché que sea bueno al mismo tiempo: hay que descomponer la matriz en trozos pequeños… en pocas palabras, un lío. Y resulta que hay una forma más sencilla de conseguirlo: invertir la posición de los índices al acceder a la segunda matriz. El código «naïve» para la multiplicación directa de dos matrices se resumía en un triple bucle:

for (int i = 0; i < m; i++)
    for (int j = 0; j < p; j++)
    {
        double d = 0;
        for (int k = 0; k < n; k++)
            d += a[i, k] * b[k, j];
        result[i, j] = d;
    }

En la entrada sobre multiplicación de matrices, invertíamos el orden de los dos bucles interiores para mejorar el rendimiento de la caché de memoria. Además, añadíamos punteros e instrucciones AVX. Para la nueva rutina, lo que invertimos son los índices de la segunda matrix y, en consecuencia, renunciamos al truco de intercambiar los bucles internos. Perdemos algo de velocidad pero, en cualquier caso, el resultado es mejor que trasponer explícitamente y luego multiplicar.

La rutina que finalmente está incluida en Austra es la siguiente:

public unsafe Matrix MultiplyTranspose(Matrix m)
{
    Contract.Requires(IsInitialized);
    Contract.Requires(m.IsInitialized);
    if (Cols != m.Cols)
        throw new MatrixSizeException();
    Contract.Ensures(Contract.Result<Matrix>().Rows == Rows);
    Contract.Ensures(Contract.Result<Matrix>().Cols == m.Rows);

    int r = Rows;
    int n = Cols;
    int c = m.Rows;
    double[,] result = new double[r, c];
    fixed (double* pA = values, pB = m.values, pC = result)
    {
        int lastBlockIndex = n & CommonMatrix.AVX_MASK;
        double* pAi = pA;
        double* pCi = pC;
        for (int i = 0; i < r; i++)
        {
            double* pBj = pB;
            for (int j = 0; j < c; j++)
            {
                int k = 0;
                double acc = 0;
                if (Avx.IsSupported)
                {
                    Vector256<double> sum = Vector256<double>.Zero;
                    for (; k < lastBlockIndex; k += 4)
                        sum = Avx.Add(
                            left: Avx.Multiply(
                                left: Avx.LoadVector256(pAi + k),
                                right: Avx.LoadVector256(pBj + k)),
                            right: sum);
                    sum = Avx.HorizontalAdd(sum, sum);
                    acc = sum.ToScalar() + sum.GetElement(2);
                }
                for (; k < n; k++)
                    acc += pAi[k] * pBj[k];
                pCi[j] = acc;
                pBj += n;
            }
            pAi += n;
            pCi += c;
        }
    }
    return result;
}

Una parte importante del truco es hacer que el lenguaje que implementa Austra reconozca la multiplicación por la traspuesta y sepa usar la operación más eficiente. Pero eso es fácil de conseguir, y lo explicaré en otro momento.

Categorías
Quantum

Mediciones

Con esta entrada puedo ganarme la enemistad de algún físico fanático. En cuanto alguno lea la frase «colapso de la función de onda» sufrirá un ataque de epilepsia. Esto tiene que ver con las interpretaciones filosóficas de la mecánica cuántica, de las que hay muchas. Por fortuna, no hay ambigüedad la forma de realizar cálculos. Y por suerte para nosotros, en Computación Cuántica tendremos que lidiar casi siempre con un único observable. Siento mucho inundarle con más contenido que el estrictamente necesario, pero mi teoría es que en conocimientos, más es siempre más.

Dos principios

En la entrada anterior presenté tres axiomas. Ahora nos toca ver dos principios adicionales.

  • La aplicación de operador observable a un sistema cuántico provoca un colapso del estado a uno de los vectores propios (eigenvectors) de dicho operador.
  • La probabilidad de que un vector propio concreto sea el elegido depende de la magnitud al cuadrado de su amplitud asociada en el estado cuántico anterior a la observación. Este principio se conoce como «regla de Born».

Estos dos principios se tratan por separado por estar asociados a un proceso físico diferente del de los axiomas. La visión del mundo (o Weltanschauung, porque me gusta el palabro alemán) que se infiere de los tres axiomas es muy diferente a la que surge de estos dos nuevos principios. En el mundo de los tres axiomas, el estado cuántico evoluciona de forma determinista y reversible. No hay nada probabilístico en la ecuación de Schrödinger. A este modelo se le conoce como «evolución unitaria», por razones que veremos más adelante.

En cambio, los dos nuevos principios corresponden al proceso conocido como medición. Se trata de un proceso irreversible y probabilístico, en el que se pierde buena parte de la información existente en la función de onda anterior a la medición.

La mayoría de los físicos creen (en el sentido de creer en los dioses del Olimpo o en Nyarlathotep) que el proceso de medición y la regla de Born podrían derivarse de los tres axiomas. Nadie ha podido demostrarlo, de momento. En el otro bando, hay quienes creen que hace falta ampliar la ecuación de Schrödinger, o que el colapso tiene que ver con discordancias en el espaciotiempo, y que hace falta una teoría cuántica de la gravitación para explicar esta parte. Personalmente, no me parece descabellado. Al fin y al cabo, la gravedad destruye información en los agujeros negros (nadie ha observado experimentalmente la radiación de Hawking). En cualquier caso, mi opinión sobre estos temas es la de un absoluto profano.

El estado como superposición

Vamos a ocuparnos del primer principio. Recordemos que, de acuerdo a los axiomas, un «observable» es un operador auto-adjunto que actúa sobre los vectores de un espacio vectorial complejo. Exigimos que sea auto-adjunto, además, para que sus valores propios sean números reales.

Para simplificar, usaré como ejemplo el sistema más simple: un ordenador cuántico con sólo 1 qubit. El estado cuántico, en este caso, es un vector complejo en $\mathbb{C}^2$. Es decir, el estado se describe mediante dos números complejos:
$$
\alpha\hat e_1 + \beta\hat e_2
$$La clave, en este punto, es averiguar qué son esos $\hat e_1$ y $\hat e_2$ que han aparecido de la nada. Son dos vectores unitarios (longitud igual a uno) que definen una «base» en el espacio vectorial. ¿Cuál base, en concreto? Necesitamos algo más de información para poder dar respuesta a esta pregunta…

La respuesta consiste en que este ordenador de un qubit debe ofrecer un «operador observable», cuya implementación exacta es cosa del hardware. Lo que nos importa es que ese operador va a definir implícitamente una base formada por sus vectores propios. ¿Cuántos vectores propios tiene un operador auto-adjunto en $\mathbb{C}^2$? Tiene dos, por supuesto: es un espacio de dos dimensiones, ¿no?. Podemos seguir llamándolos $\hat e_1$ y $\hat e_2$… o podemos ser más prácticos y llamarlos $\vert 0\rangle$ y $\vert 1\rangle$. Con esto, ya podemos afinar un poco más la definición del estado cuántico del qubit. Siempre podrá describirse como una combinación lineal de estos dos vectores o estados especiales:
$$
\alpha\vert 0\rangle + \beta\vert 1\rangle
$$La base formada por los vectores $\vert 0\rangle$ y $\vert 1\rangle$ se conoce como base computacional o, en inglés, computational basis. Al observable cuyo operador tiene esta base como vectores propios, voy a llamarlo, por motivos que desvelaremos a su debido tiempo, $M_z$.

La Regla de Born

Vamos a darle valores concretos a $\alpha$ y $\beta$. Supongamos que el estado del qubit es el siguiente:
$$
\vert 0\rangle + \vert 1\rangle
$$La regla de Born nos dice que, si aplicamos la medición $M_z$ a este estado, obtendremos como resultado el vector $\vert 0\rangle$ en la mitad de los experimentos, y $\vert 1\rangle$ en la otra mitad de las veces. Supongamos, por el contrario, que el estado es:
$$
3\vert 0\rangle + 4\vert 1\rangle
$$En este caso, obtendremos $\vert 0\rangle$ con una probabilidad de $\frac{9}{25}$, esto es, un 36% de los casos. Y obtendremos $\vert 1\rangle$ con una probabilidad de $\frac{16}{25}$, que es el 64%.

¿Y si el estado inicial es directamente $\vert 0\rangle$? Pues en este caso, siempre saldrá $\vert 0\rangle$ como resultado. En general, si el estado cuántico es éste:
$$
\alpha\vert 0\rangle + \beta\vert 1\rangle
$$la probabilidad de que obtengamos $\vert 0\rangle$ es
$$
\frac{\alpha^2}{\alpha^2+\beta^2}
$$y la probabilidad de que salga $\vert 1\rangle$ es la complementaria.

Normalización

Un médico no debe hacer daño, y un escritor técnico nunca debe confundir al lector. A quienes ya conocen algo de Mecánica Cuántica les extrañará que haya puesto el siguiente estado como ejemplo:
$$
\vert 0\rangle + \vert 1\rangle
$$¿Por qué? Pues porque los estados suelen representarse de manera que su longitud sea la unidad: $\vert\vert\psi\vert\vert =1$. Se trata de un convenio, simplemente. El estado anterior normalmente se escribe así:
$$
\frac{1}{\sqrt 2}\vert 0\rangle + \frac{1}{\sqrt 2}\vert 1\rangle
$$En general, como convenio se pide lo siguiente:
$$
\alpha\vert 0\rangle + \beta\vert 1\rangle\quad\alpha\alpha^* + \beta\beta^* = 1
$$Como $\alpha$ y $\beta$ son números complejos, hemos tenido que multiplicarlos por sus respectivos conjugados. La diferencia entre normalizar los estados y no normalizarlos consiste en que, si no los normalizamos, tendremos que hacer malabares con la norma del estado en algunas fórmulas. Por lo tanto, de ahora en adelante, el estado que representábamos como $3\vert 0\rangle + 4\vert 1\rangle$ lo escribiremos como $0.6\vert 0\rangle + 0.8\vert 1\rangle$ para evitar problemas.

En una exposición más rigurosa de la Mecánica Cuántica tendríamos que haber empezado diciendo que el estado cuántico se representa mediante un «rayo» en $\mathbb{C}^n$. O más oscuramente, que es un elemento de un espacio proyectivo de $\mathbb{C}^n$. Chino mandarín, vamos, pero ya sabemos qué es lo que quieren decir.

Categorías
Quantum

Estado cuántico

Este blog se llama Quantum Insights porque mi intención inicial era dedicarlo a la Computación Cuántica. Me he distraído un poco con los preliminares, pero es hora ya de saltar a la materia que nos interesa. No nos hará falta saber mucha Física para aprender, pero un poco no nos vendrá mal tampoco.

Mi plan para las próximas entradas es el siguiente: primero, voy a explicar qué es el «estado cuántico», a grandes rasgos. Segunda entrada: el proceso de medición. En la tercera entrada ya veremos, entonces, cómo se define un ordenador cuántico. Entiendo que, por mucho que simplifique, siempre hay temas que necesitarán aclaraciones. Utilice los comentarios ad libitum, y si lo cree necesario, envíeme un correo electrónico. Mi cuenta es mi nombre, Ian, más el dominio de este blog.

Tres axiomas

La primera mitad de la Mecánica Cuántica se explica con estos tres axiomas:

  • El estado de un sistema cuántico se describe mediante un vector en un espacio vectorial complejo $\cal H$, dotado de un producto interior hermitiano.
  • Los observables del sistema se corresponden con operadores lineales auto-adjuntos en $\cal H$.
  • La evolución en el tiempo del estado cuántico está determinada por la ecuación de Schrödinger.

Hay un montón de términos técnicos, y por ello vamos a dedicar varias secciones a cada axioma.

Espacios de Hilbert (I)

Lo primero es ver qué es un espacio vectorial complejo. Llevamos unas cuantas entradas hablando de vectores a secas, por lo que imagino que el concepto es más o menos intuitivo. Los vectores más populares son los vectores euclidianos: tres valores reales, como en $(1.2,\,2,\,-33)$. No obstante, los vectores que nos interesan para la Mecánica Cuántica tienen dos diferencias importantes:

  1. Cada componente de estos vectores va a ser un número complejo, en vez de un número real.
  2. En el espacio euclidiano hay tres dimensiones. El estado cuántico puede tener, dependiendo del sistema que se estudie, un número diferente de dimensiones. Puede ser un número finito o infinito de dimensiones. Y cuando hay infinitas dimensiones, puede tratarse de infinito numerable o infinito no numerable. Sí: manda…

No se asuste: los estados cuánticos que se manejan en Computación Cuántica son espacios finitos, y el número de dimensiones es una potencia de dos, como en $2^N$, donde $N$ es el número de qubits del sistema. Por ejemplo, dedicaremos algún tiempo a estudiar el sistema de 1 qubit, por ser el más sencillo posible. El estado de un qubit, por lo antes dicho, se puede representar como un vector de dos dimensiones complejas, como estos ejemplos:
$$
\eqalign{(0,&\,1)\cr
(0.7071 + 0.7071i,&\, 0.7071 – 0.7071i)}
$$Estos espacios vectoriales complejos, una vez que definimos un producto interior hermitiano (no lo hemos hecho todavía) se conocen como espacios de Hilbert, si tienen una propiedad adicional: que sean espacios métricos completos. Los espacios que se utilizan en computación cuántica tienen dimensiones finitas, y cumplen automáticamente con esta regla. De ahí, la $\cal H$ caligráfica que se menciona en los axiomas.

Producto interior hermitiano (I)

Dos vectores se pueden sumar y restar entre sí, y es fácil imaginar cómo se hace: componente a componente. También se puede multiplicar un vector por un escalar, casi como con los vectores euclidianos. La diferencia está en que el escalar ahora es un número complejo. Esta multiplicación, naturalmente, también se calcula componente a componente. Por ejemplo:
$$\eqalign{
(0,\, 1)+(i,\, 0)=&(i,\, 1)\cr
i\cdot(i,\, 1)=&(-1,\, i)}
$$Los vectores euclidianos tienen una operación de multiplicación escalar entre vectores, que recibe dos vectores y devuelve un número real. Cuando el espacio vectorial es complejo, sin embargo, esta operación se complica. El problema está en cómo se define la longitud de un vector. En un espacio euclidiano, acostumbramos a calcular el producto interior del vector consigo mismo, y aplicarle entonces la raíz cuadrada:
$$
\vert\vert v \vert\vert = \sqrt{v \cdot v}
$$Si queremos que las longitudes de los vectores complejos sean reales y positivas, entonces tenemos que ajustar la definición del producto interior para que, al menos, el producto de un vector consigo mismo sea no solamente real, sino además positivo. La generalización necesaria es sencilla. Supongamos que tenemos un par de vectores, $x$ e $y$, con componentes reales. En este caso, el producto interior se define clásicamente así:
$$
\eqalign{
x =& [x_1,x_2,x_3\cdots x_n]\cr
y =& [y_1,y_2,y_3\cdots y_n]\cr
x \cdot y =& \sum_{i=1}^n{x_i y_i}
}$$Estos son los productos escalares que hemos estado calculando a gogó en entradas anteriores.

Ahora supongamos que los componentes de estos vectores son complejos. La conjugada de un número complejo se define como un segundo número complejo con la misma parte real y la negación de la parte imaginaria del número original:
$$
(a + b i)^* = a – b i
$$El producto interior de dos vectores con componentes complejos es, entonces:
$$
x \cdot y = \sum_{i=1}^n{x_i^* y_i}
$$En el caso de los vectores reales, el producto interior es simétrico. Esto es, $x \cdot y = y \cdot x$. Pero para vectores complejos, la propiedad que se cumple es la siguiente:
$$
x \cdot y = (y \cdot x)^*
$$La igualdad anterior es fácil de demostrar si nos vamos a la definición de producto interior en espacios complejos. Lo que nos importa ahora es lo que ocurre cuando se toma el producto interior de un vector complejo consigo mismo:
$$
x \cdot x = (x \cdot x)^*
$$Esto significa que el valor del producto interior de un vector consigo mismo es, a la vez, igual a su valor conjugado. Pero esto sólo puede ocurrir cuando el valor es real, o sea, cuando la parte imaginaria es cero. Podemos incluso ir más lejos, y demostrar que el producto interior complejo que acabamos de definir es siempre no negativo cuando se multiplica cualquier vector consigo mismo. Por lo tanto, podemos seguir definiendo la longitud de un vector en un espacio complejo como antes:
$$
\vert\vert x \vert\vert = \sqrt{x \cdot x}
$$

La notación de Dirac (I)

Hagamos una pequeña pausa: resulta que para calcular un producto escalar, necesitamos una versión modificada de uno de los vectores… y no del otro. Pero el operador que hemos utilizado para el producto interior (el punto) «sugiere» que se trata de un operador simétrico (y no lo es). A Paul Adrien Maurice Dirac se le ocurrió una idea: digamos que a todo espacio vectorial $\mathbb{C}^n$ le corresponde automáticamente un espacio vectorial dual, y que entre ambos espacios hay una transformación biunívoca. El dual de un vector es un vector con los mismos componentes, pero conjugados. A los elementos del espacio vectorial original los llamaremos «kets» y los escribiremos de esta manera:
$$
\vert \psi \rangle
$$A los vectores del espacio dual los llamaremos «bras», y el dual del «ket» anterior se representa así:
$$
\langle \psi \vert
$$Lo que haremos a continuación es definir el producto interior como una operación que siempre toma un primer operando de tipo «bra» y un segundo operando de tipo «ket». Si yo quiero calcular el producto interior de $\vert \phi \rangle$ con $\vert\psi \rangle$, tengo que convertir antes el primer operando en un «bra», y sólo entonces puedo obtener el producto interior:
$$
\langle \phi \vert \psi \rangle
$$La expresión anterior es un número complejo, un escalar, y se cumple la antisimetría conjugada:
$$
\langle \phi \vert \psi \rangle = \langle \psi \vert \phi \rangle ^*
$$Naturalmente, la longitud de un vector, no importa si es bra o ket, puede definirse así:
$$
\vert\vert \psi \vert\vert = \sqrt{\langle\psi\vert\psi\rangle}
$$Siendo rigurosos, esto es un poco de abuso de notación. Pero no nos supondrá problema alguno.

Observables (II)

Toca explicar qué narices es el operador lineal auto-adjunto del segundo axioma. No es complicado: si el espacio vectorial es un espacio de dimensiones finitas, como en Computación Cuántica, un operador es simplemente una matriz con componentes complejos. Si tuviésemos que tratar con espacios de infinitas dimensiones, tendríamos que hilar un poco más fino, pero no es necesario en nuestro caso.

¿Recuerda que nuestro producto interior no es simétrico? Este detalle provoca que, en general, aplicar un operador al bra y al ket tengan resultados diferentes. En términos generales:
$$
\langle Ax\vert y\rangle\neq\langle x\vert Ay\rangle
$$Si queremos mover el operador $A$ al otro lado de la barra central, tenemos que transformar el operador $A$ en su adjunto $A^{*}$:
$$
\langle Ax\vert y\rangle = \langle x\vert A^{*}y\rangle
$$A nivel de celdas, la adjunta de una matriz es una matriz transpuesta creada a partir de los valores conjugados de sus componentes. Por ejemplo:
$$
\pmatrix{1&-2-i\cr 1+i&i}^* = \pmatrix{1&1-i\cr -2+i&-i}
$$Un operador auto-adjunto, simplemente, es un operador que no cambia al calcular su adjunto: $A = A^*$. Por ejemplo:
$$
\pmatrix{0&-i\cr i&0}^* = \pmatrix{0&-i\cr i&0}
$$Por lo tanto, si volvemos al producto interior, si el operador $A$ es auto-adjunto, se cumple que:
$$
\langle Ax\vert y\rangle = \langle x\vert Ay\rangle
$$¿Qué importancia tienen los operadores auto-adjuntos? Pues que los valores propios, o eigenvalues, de un operador auto-adjunto en un espacio vectorial complejo, son siempre valores reales. Este es el dato técnico. Pasemos a la interpretación física:

  1. Un «observable» de un sistema cuántico es simplemente una propiedad física del sistema que podemos medir. Ejemplos: la posición de una partícula, la velocidad de una partícula, la orientación del espín respecto a una dirección predeterminada, etc, etc.
  2. Nosotros no vamos a medir la velocidad de una partícula en Computación Cuántica. Lo advierto para la salud mental de todos nosotros. El «observable» que vamos a manejar prácticamente siempre en un ordenador cuántico es el estado binario de sus qubits. Esto lo veremos con más detalles en el momento adecuado.
  3. Como bien dice el axioma, cada «observable» se asocia a un operador auto-adjunto.
  4. El valor del observable se obtiene mediante un proceso de «medición», que veremos en la próxima entrada.
  5. Matemáticamente, la medición consiste en escoger probabilísticamente uno de los valores propios del operador asociado al observable.
  6. Como las cantidades físicas suelen ser magnitudes reales (las anoréxicas tienen una masa compleja, con una parte real y otra imaginaria), tenemos que exigir que los operadores observables tengan esta conveniente propiedad de ser auto-adjuntos.

La ecuación de Schrödinger (III)

Ya llegamos al tercer axioma, que es donde se menciona por primera vez la ecuación de Schrödinger:$$
i\hbar{d \over dt}\vert \psi(t) \rangle=H\vert \psi(t) \rangle
$$Tengo una buena noticia: ¡no necesitaremos resolver la ecuación de Schrödinger! Al menos, mientras no tengamos que enredarnos con el hardware a muy bajo nivel, claro. Sin embargo, menciono la dichosa ecuación porque conocerla nos va a ayudar a comprender mejor las reglas de la Computación Cuántica.

A la derecha de la igualdad tenemos la derivada temporal de la función de onda. Y a la izquierda, el operador de Hamilton del sistema. En pocas palabras: se trata de una ecuación lineal. Si $\vert \psi_0\rangle$ es una solución de la ecuación, y también lo es $\vert \psi_1\rangle$, entonces cualquier combinación
$$\alpha\vert \psi_0\rangle + \beta\vert \psi_1\rangle
$$donde $\alpha$ y $\beta$ son números reales, vale también como solución.

En la próxima entrada de este blog, trataremos el modelo de medición sobre un sistema cuántico.

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.