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
FinTech

Covarianza

Recordemos cómo se define la varianza de una variable aleatoria X:
$$
Cov(X, Y) = \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
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
C#

Entran una matriz y un vector en un bar

… y claro, al rato sale un vector «transformado».

Esta entrada no es, aunque pueda parecerlo, un ripio de la anterior. Algorítmicamente, transformar un vector con una matriz se parece mucho a una sucesión de productos escalares. Pero resulta que el producto escalar, al menos hasta AVX2, tiene su truco. Vamos a comenzar por la implementación más tonta:

public static double[] Mult(double[,] a, double[] x)
{
    int m = a.GetLength(0);
    int n = a.GetLength(1);
    double[] b = new double[m];
    for (int = 0; i < m; i++)
    {
        double d = 0;
        for (int j = 0; j < n; j++)
            d += a[i, j] * x[j];
        b[i] = d;
    }
    return b;
}

Recordemos que tenemos un «handicap» autoimpuesto por representar las matrices como arrays bidimensionales de C#. Pero esta vez no voy a dar la brasa con los punteros, que ya sabemos que resuelven este problema sin pestañear. Esta es la implementación final que necesitamos, con soporte opcional de AVX para cuando esté disponible y merezca la pena:

public static unsafe double[] Mult(double[,] a, double[] x)
{
    int m = a.GetLength(0);
    int n = a.GetLength(1);
    double[] b = new double[m];
    int lastBlockIndex = n - (n % 4);
    fixed (double* pA = a)
    fixed (double* pX = x)
    fixed (double* pB = b)
    {
        double* pA1 = pA;
        double* pB1 = pB;
        if (n >= 12 && Avx2.IsSupported)
            for (int i = 0; i < m; i++)
            {
                int j = 0;
                var v = Vector256<double>.Zero;
                while (j < lastBlockIndex)
                {
                    v = Avx.Add(
                        v,
                        Avx.Multiply(
                            Avx.LoadVector256(pA1 + j),
                            Avx.LoadVector256(pX + j)));
                    j += 4;
                }
                v = Avx.HorizontalAdd(v, v);
                double d = v.ToScalar() + v.GetElement(2);
                for (; j < n; j++)
                    d += pA1[j] * pX[j];
                *pB1 = d;
                pA1 += n;
                pB1++;
            }
        else
            for (int i = 0; i < m; i++)
            {
                int j = 0;
                double d = 0;
                while (j < lastBlockIndex)
                {
                    d += (*(pA1 + j) * *(pX + j)) +
                        (*(pA1 + j + 1) * *(pX + j + 1)) +
                        (*(pA1 + j + 2) * *(pX + j + 2)) +
                        (*(pA1 + j + 3) * *(pX + j + 3));
                    j += 4;
                }
                for (; j < n; j++)
                     d += pA1[j] * pX[j];
                *pB1 = d;
                pA1 += n;
                pB1++;
            }
    }
    return b;
}

Esta vez, el código SIMD sólo se usa cuando hay doce o más elementos en el vector. La cifra la he elegido experimentando en mi i7-4770. Puede que en otros ordenadores, el umbral sea más bajo incluso.

Tengo que explicar cómo se implementa un producto escalar con SIMD, porque no es muy evidente. Uno diría que hay que acumular un escalar en una variable global al bucle… pero no hay ninguna operación SIMD que calcule directamente la suma de las cuatro multiplicaciones necesarias. La explicación oficial es que una suma de ese tipo destrozaría el paralelismo de la CPU. Y yo me lo creo, de veras. La consecuencia es que necesitamos acumular las multiplicaciones en cuatro variables; es decir, en un vector que hace de acumulador.

Las cosas se ponen de color hormiga cuando terminamos el bucle y tenemos entonces que sumar los cuatro elementos del vector acumulador. ¿Ves la zona de código resaltada en verde (lo siento si eres daltónico)? Según mis experimentos, es la forma más rápida de conseguirlo. HorizontalAdd, cuando se trata de Vector256<double>, suma el primer elemento con el segundo, y lo almacena por partida doble en el primer y segundo elemento. A la vez, suma el tercero y el cuarto y hace lo mismo para guardar el resultado. Los métodos de extensión ToScalar() y GetElement() acceden entonces directamente al primer y tercer elemento y los suma. Mantengo la llamada inicial a HorizontalAdd porque, teóricamente, puede hacer dos de las sumas en paralelo, pero puedes experimentar a ver qué pasa si accedes directamente a los cuatro elementos y los sumas como toda la vida. A mí ya se me ha acabado la partida de tiempo libre para este experimento.

La razón para la controversia es que, en realidad, Internet está lleno de recomendaciones para hacer esta suma final de esta otra manera:

v = Avx2.Permute4x64(
    Avx.HorizontalAdd(v, v),
    0b00_10_01_11);
double d = Avx.HorizontalAdd(v, v).ToScalar();
// v = Avx.HorizontalAdd(v, v);
// double d = v.ToScalar() + v.GetElement(2);

Es decir: se llama dos veces a HorizontalAdd, pasando entre medias por una permutación entre escalares. En la arquitectura Haswell, al menos, esto funciona más lento que mi solución.

Si multiplico una matriz aleatoria de 64×64 por un vector de 64 elementos, obtengo estas cifras:

Method Mean Error StdDev Median
MultVector 5.762 μs 0.1142 μs 0.2227 μs 5.646 μs
FMultVector 1.814 μs 0.0320 μs 0.0416 μs 1.818 μs

No está mal, aunque no conseguimos tanta ventaja como con la multiplicación entre matrices. La versión con punteros y sin SIMD tampoco va mal, pero queda muy claro que el SIMD acelera este código. De paso, ya tenemos un patrón de código para productos escalares (y para cosas más raras como multiplicar un vector de sensibilidad delta-gamma por un escenario histórico: cosas de la valoración de productos financieros).

Por cierto, el mejor chiste que conozco sobre gente que entra en un bar tiene que ver con la Mecánica Cuántica. Dice así: entra el Gato de Schrödinger en un bar… y no entra.

Categorías
C#

Multiplicación de matrices

Supongamos que queremos multiplicar un par de matrices, $A$ y $B$. Digamos que la primera tiene dimensiones $m\times n$ y que la segunda es $n\times p$. La coincidencia entre columnas de la primera y filas de la segunda es condición necesaria para que podamos multiplicarlas.

Si me piden que escriba de carrerilla un método para esta multiplicación, esto es lo que se me ocurre:

public static double[,] Mult(double[,] a, double[,] b)
{
    int m = a.GetLength(0);
    int n = a.GetLength(1);
    int p = b.GetLength(1);
    double[,] result = new double[m, p];
    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;
        }
    return result;
}

He utilizado matrices bidimensionales de C# porque acceder a sus elementos individuales es sencillo. Internamente, C# las almacena en una sola memoria contigua de memoria, fila por fila.

El código que he mostrado no es una maravilla. Para empezar, cada vez que decimos algo como a[i, k], el compilador tiene que multiplicar la variable i por el número de columnas y por los ocho bytes que tiene un flotante de doble precisión. Hacerlo una vez no es problema… pero tenemos tres bucles anidados. Eso tiene que doler. Si en vez de C# escribiésemos esto en C++, el compilador podría sustituir un montón de multiplicaciones por sumas. RyuJIT ha mejorado muchísimo, pero no tanto.

C#, además, es un lenguaje mucho más seguro que C++, pero esta seguridad nos cuesta un montón de verificaciones de rango para poder indexar. Recordemos, además, que cada acceso necesita dos índices.

Y hay un tercer problema, mucho más sutil: cuando las matrices son grandes, el código anterior machaca la caché de la CPU sin piedad. Toma un folio de papel y haz el experimento: dibuja dos matrices, y ve numerando las celdas siguiendo el orden en que las usa el algoritmo.

La clase Unsafe

Llegados a este punto, tenemos dos alternativas: o marcamos el método como unsafe y usamos directamente punteros de C#, o intentamos evitarlo haciendo uso de la clase Unsafe, de System.Runtime.CompilerServices. Vamos a comenzar por esta última. De paso, voy a invertir el orden de los dos bucles más internos, para ver qué conseguimos con ello. Este es el código modificado, y suele funcionar el doble de rápido, o un poco más:

public static double[,] Mult(double[,] a, double[,] b)
{
    int m = a.GetLength(0);
    int n = a.GetLength(1);
    int p = b.GetLength(1);
    double[,] c= new double[m, p];
    ref double rA = ref a[0, 0];
    ref double rB = ref b[0, 0];
    ref double rC = ref c[0, 0];
    for (int i = 0; i < m; i++)
    {
        ref double rAi = ref Unsafe.Add(ref rA, i * n);
        ref double rCi = ref Unsafe.Add(ref rC, i * n);
        for (int k = 0; k < n; k++)
        {
            double d = Unsafe.Add(ref rAi, k);
            int kp = k * p;
            for (int j = 0; j < p; j++)
                Unsafe.Add(ref rCi, j) +=
                    d * Unsafe.Add(ref rB, kp + j);
        }
    }
    return c;
}

La regla principal del uso de Unsafe.Add es que si inicializamos así:

ref double rA = ref a[0, 0];

entonces el acceso a a[i, j] debe parecerse a esto:

Unsafe.Add(ref rA, i * n + j) = 42;

Esa multiplicación es un problema del que ya advertimos. En nuestro código lo paliamos moviendo la multiplicación al inicio del bucle donde se le da valor al índice de la fila. Mi apaño no es la palabra definitiva: le dejo como ejercicio la eliminación total de esas multiplicaciones.

Ahora hay que prestar atención, sobre todo, al patrón de acceso a memoria que se produce en el bucle más interno. En el algoritmo inicial, acumulábamos todos los términos de un elemento de la matriz final en el bucle interno, y asignábamos su suma de golpe a la celda del resultado. Esta variante, sin embargo, no parece tan buena. Tenemos que asumir que, al reservar memoria para la matriz, todas sus entradas valen cero (y es así). Luego, cada celda del resultado se va rellenando por pasos, no de una vez. Puede que esto sea bueno para la caché de la CPU, pero no me queda tan claro que sea bueno para el compilador de C#.

Pero lo que nos interesa realmente es que ahora ejecutamos el siguiente patrón de cálculo:

  1. Tenemos dos zonas de memoria consecutiva.
  2. Leemos algo de la primera zona.
  3. Lo transformamos como sea.
  4. Lo asignamos a la celda equivalente en la segunda zona de memoria.

Instrucciones SIMD

Ese patrón de actividad es el típico algoritmo «vectorial» que podemos acelerar utilizando operaciones SIMD. Tenemos dos opciones:

  1. Utilizar System.Numerics.Vector<T>, que se adapta automáticamente a cualquier máquina que soporte SIMD, e incluso ofrece una alternativa cuando no existe ese soporte. Este tipo funciona también para .NET Framework, a través de un paquete.
  2. Si podemos usar .NET Core 3.1, podemos ir directamente a las clases declaradas en System.Runtime.Intrinsics y System.Runtime.Intrinsics.X86. Es un poco más complicado y no está bien documentado, pero da resultados ligeramente mejores.

Vamos a ir directamente por la segunda vía. Vamos a optimizar las CPUs que soporten el conjunto de instrucciones AVX, haremos algo más en el caso en que soporte el conjunto FMA (que mezcla multiplicaciones y sumas en una misma operación) y, de todas maneras, habilitaremos código de respaldo para cuando el procesador no soporte SIMD.

Cuando hay soporte para instrucciones AVX, podemos procesar hasta cuatro variables de tipo double de una tacada. Para ello tenemos que utilizar el tipo de estructura Vector256<double>, que tiene capacidad para cuatro elementos. La forma más sencilla de inicializar estos vectores es utilizando punteros, por lo que vamos a tener que declarar nuestro método unsafe y pasarnos directamente a los punteros.

public static unsafe double[,] Mult(double[,] a, double[,] b)
{
    int m = a.GetLength(0);
    int n = a.GetLength(1);
    int p = b.GetLength(1);
    double[,] c = new double[m, p];
    int lastBlockIndex = p - (p % 4);
    fixed (double* pA = a)
    fixed (double* pB = b)
    fixed (double* pC = c)
    {
        double* pAi = pA;
        double* pCi = pC;
        for (int i = 0; i < m; i++)
        {
            double* pBk = pB;
            for (int k = 0; k < n; k++)
            {
                double d = *(pAi + k);
                if (Avx.IsSupported)
                {
                    int j = 0;
                    var vd = Vector256.Create(d);
                    while (j < lastBlockIndex)
                    {
                        if (Fma.IsSupported)
                            Avx.Store(pCi + j,
                                Fma.MultiplyAdd(
                                Avx.LoadVector256(pBk + j),
                                vd,
                                Avx.LoadVector256(pCi + j)));
                        else
                            Avx.Store(pCi + j,
                                Avx.Add(
                                Avx.LoadVector256(pCi + j),
                                Avx.Multiply(
                                Avx.LoadVector256(pBk + j),
                                vd)));
                        j += 4;
                    }
                    while (j < p)
                    {
                        pCi[j] += d * pBk[j];
                        j++;
                    }
                }
                else
                {
                    for (int j = 0; j < p; j++)
                        pCi[j] += d * pBk[j];
                }
                pBk += p;
            }
            pAi += n;
            pCi += p;
        }
    }
    return c;
}

Observaciones:

  1. Lo peor de trabajar con SIMD es tener que lidiar con vectores que no son múltiplos exactos del tamaño del vector básico. Nuestros vectores básicos tienen cuatro elementos. Si tenemos un vector de 75 elementos, necesitaremos un bucle de 18 repeticiones que procese cuatro elementos por vez, para una mierdecilla de bucle final que maneje los 3 elementos que nos sobran.
  2. Aunque la llamada a Avx.IsSupported está metida dentro de dos bucles anidados, no se preocupe: el compilador JIT la trata como una constante en tiempo de generación de código nativo, y no cuesta nada. Si no se soporta AVX, el compilador JIT solamente genera el código de la cláusula else, que funciona sobre cualquier arquitectura.
  3. Ojo: ese código «para cualquier máquina» podría optimizarse echando mano de la técnica de loop unrolling. Pero mi política en estos casos es: si no tienes una máquina decente, jódete.
  4. En el ejemplo anterior, cuando intercambiamos el orden de los bucles más internos, teníamos un valor escalar que sacábamos fuera del tercer bucle. Pero SIMD no ofrece instrucciones para multiplicar un vector por un escalar: tenemos que convertir ese escalar en todo un vector y utilizar la instrucción de multiplicación más general. No es grave, de todos modos.
  5. Si, además de AVX, la máquina soporta el conjunto FMA de instrucciones, podemos utilizar el método MultiplyAdd para acelerar un poco el algoritmo. Pero con esto hay que tener cuidado: a * b + c puede dar resultados diferentes si se hacen las dos operaciones por separado o a la vez. Si se hacen a la vez, aumenta la exactitud de la operación al existir menos redondeos. Pero el efecto secundario es que los cálculos con y sin esa opción dan resultados ligeramente diferentes. Tenemos que decidir cuándo es aceptable que exista esa diferencia y cuándo no. En cualquier caso, tengamos presente que el resultado de MultiplyAdd es más preciso.

Benchmark.NET

Para estar seguro de las ganancias en velocidad, he utilizado el package Benchmark.NET para generar las pruebas. Estos son los resultados:

Method Mean Error StdDev
MultMatrix 4,482.3 μs 88.75 μs 138.17 μs
UMultMatrix 1,895.2 μs 37.87 μs 63.26 μs
FMultMatrix 506.3 μs 3.44 μs 2.87 μs

La mejora por el uso de SIMD es cerca de cuatro veces, porque es el número de operaciones simultáneas que permite esta arquitectura en particular. Con AVX512 tendríamos vectores de ocho valores, pero necesitaríamos procesadores mucho más modernos, y de momento .NET Core no lo soporta.

Para esta prueba, he utilizado matrices de 128×128. He probado también con matrices de 8×8 e incluso de 4×4. La ganancia no es tan espectacular, pero en total se consigue una cuarta parte del tiempo de ejecución respecto al algoritmo más sencillo.