Etiqueta: matrices

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

Transformaciones unitarias

Un algoritmo cuántico funciona partiendo de un estado inicial en el que todos los qubits están en el estado $\vert 0\rangle$. A partir de ese punto, el estado se transforma mediante matrices, que se implementan físicamente como circuitos lógicos. El objetivo del algoritmo es llevar el estado cuántico a determinado valor. Entonces se realiza una medición.

La «explicación» anterior no «explica» mucho en este momento, pero necesitaba incluirla para saber en qué parte nos encontramos, y hacia dónde vamos. En esta entrada comenzaremos a ver las transformaciones más sencillas que afectan al estado cuántico: aquellas que solamente afectan un qubit cada vez.

Rotaciones complejas

Sabemos que el estado cuántico es un vector complejo. En el caso de un solo qubit, el vector tiene dos componentes. Veremos más adelante que, para $N$ qubits, tendremos $2^N$ componentes. Es fácil ver que, si queremos transformar estos vectores, necesitamos matrices complejas. Lo importante es qué tipo de matrices nos valen: necesitamos matrices unitarias, que son el equivalente de las rotaciones (matrices ortogonales) en un espacio vectorial complejo.

Matemáticamente, si tenemos una matriz compleja $U$, y llamamos $U^*$ a su transpuesta conjugada (cambiamos filas y columnas e invertimos el signo de las partes imaginarias), entonces la matriz es unitaria si se cumple:
$$UU^* = I$$donde $I$ es la matriz identidad. Es decir: toda matriz unitaria tiene una inversa (que es su transpuesta conjugada, en concreto). Este es un detalle importante: un estado cuántico regido por la ecuación de Schrödinger es simétrico respecto al paso del tiempo. Toda transformación cuántica que no implique una medición tiene que ser reversible.

Una propiedad interesante que tienen las matrices unitarias es que preservan el producto escalar entre dos vectores arbitrarios:

$$\langle U\phi\vert U\psi\rangle = \langle\phi\vert\psi\rangle$$Como el producto escalar es proporcional al ángulo entre dos vectores (ángulos complejos, en este caso), tenemos entonces que las matrices unitarias respetan el ángulo entre vectores. La experiencia en tres dimensiones reales nos dice que las rotaciones, precisamente, son las transformaciones que respetan los ángulos. En particular, si se preserva el producto escalar, también se preservará la longitud de los vectores transformados:

$$\vert\vert U\psi\vert\vert = \vert\vert\psi\vert\vert$$

Matrices de Pauli

Ahora tenemos que aprendernos un pequeño repertorio de matrices que actúan sobre un qubit. Las tres primeras matrices que veremos se conocen como matrices de Pauli, y la forma más fácil de recordarlas es utilizar la esfera de Bloch como recurso mnemotécnico:

Hay tres ejes espaciales en esta representación, y a cada uno de ellos le corresponde una matriz: $X$, $Y$ y $Z$. Cada una de ellas representa una rotación de 180 grados alrededor del eje correspondiente. Ojo: estoy hablando de rotaciones en el espacio de la esfera de Bloch, no en el espacio vectorial complejo.

Por ejemplo, la matriz $X$ gira 180 grados el estado cuántico alrededor del eje $X$. Como da la casualidad de que nuestra base computacional está alineada según el eje $Z$, esta transformación transforma el polo norte, el estado $\vert 0\rangle$, en el polo sur $\vert 1\rangle$ y viceversa. En el espacio vectorial complejo, esto equivale a intercambiar los dos componentes complejos, por lo que no nos debe extrañar que la matriz $X$ se escriba de esta manera en la base computacional:

$$X=\pmatrix{0&1\cr 1&0}$$Machaquemos esta información para estar seguro de entenderla. Imaginemos que el estado cuántico en un momento dado es, en la notación de Dirac:

$$\alpha\vert 0\rangle + \beta\vert 1\rangle$$Si le aplicamos la matriz $X$ a este vector, obtenemos un estado con los componentes intercambiados:

$$\beta\vert 0\rangle + \alpha\vert 1\rangle$$Ahora repitamos el cálculo, pero usando directamente la notación matricial:

$$\pmatrix{0&1\cr 1&0} \times \pmatrix{\alpha\cr\beta} = \pmatrix{\beta\cr\alpha}$$Son dos maneras equivalentes de expresar la misma operación.

¿Qué operación clásica convierte los ceros en unos y viceversa? La negación lógica, por supuesto. Podemos entonces considerar que $X$ es la negación dentro de nuestro repertorio de operaciones cuánticas.

Hemos elegido como negación la rotación alrededor del eje $X$, pero podríamos haber elegido también la rotación alrededor de $Y$. Pero no lo hicimos, porque la representación de $Y$ en la base computacional es un poco más complicada:

$$Y=\pmatrix{0&-i\cr i&0}$$Observe que yo no le estoy explicando que la rotación alrededor del eje $Y$ en el espacio de Bloch se corresponde realmente a la matriz que estoy presentando. No es difícil de demostrar, pero para abreviar, se lo dejo de momento a su cargo. Lo que ahora nos interesa es si de verdad esta forma tan rara con números imaginarios es realmente una negación. En realidad, es equivalente solamente cuando se trata de transformar los polos. Pero cualquier otro vector dará un vector con una fase rotada. Simplemente, vamos a preferir $X$ como forma estándar de negación.

La tercera matriz de Pauli es, naturalmente, la matriz $Z$:

$$Z=\pmatrix{1&0\cr 0&-1}$$Obviamente, si aplicamos una rotación de 180 grados alrededor del eje $Z$ a uno de los polos, nos quedaremos como al principio.

La transformación de Hadamard

Las tres matrices de Pauli sólo nos permiten, de momento, movernos de un polo al otro. ¿Qué tal si quisiéramos mover el estado cuántico al ecuador de la esfera de Bloch? Para eso necesitaremos una de las transformaciones más populares en computación cuántica: la transformación de Hadamard, o de Walsh-Hadamard.

$$H=\frac{1}{\sqrt{2}}\pmatrix{1&1\cr 1&-1}$$Veamos qué efecto tiene esta transformación sobre los vectores de la base computacional:

$$\eqalign{
\frac{1}{\sqrt{2}}\pmatrix{1&1\cr 1&-1}\times\pmatrix{1\cr 0} =& \frac{1}{\sqrt{2}}\pmatrix{1\cr 1} \cr
\frac{1}{\sqrt{2}}\pmatrix{1&1\cr 1&-1}\times\pmatrix{0\cr -1} =& \frac{1}{\sqrt{2}}\pmatrix{1\cr -1}}
$$Utilizando la notación de Dirac:

$$\eqalign{
H\vert 0\rangle =& \frac{1}{\sqrt{2}}(\vert0 \rangle + \vert1 \rangle)\cr
H\vert 1\rangle =& \frac{1}{\sqrt{2}}(\vert0 \rangle – \vert1 \rangle)
}$$O, si recordamos la definición de la base de Hadamard:

$$\eqalign{
H\vert 0\rangle =& \vert+ \rangle\cr
H\vert 1\rangle =& \vert- \rangle
}$$Es decir: la transformación de Hadamard mueve los polos a puntos situados en el ecuador de la esfera de Bloch y alineados con el eje $X$.

Como ejercicio, puede comprobar qué ocurre cuando se le aplica la matriz $X$ a los vectores de la base de Hadamard. No olvide la existencia del fenómeno de la equivalencia de fase: si multiplica los dos componentes por un vector unitario, el estado cuántico que se obtiene es indistinguible del original.

¿Unitarias o hermitianas?

Las cuatro matrices que hemos visto comparten una característica que, aunque no tiene mayor importancia, puede provocar confusión. Las hemos presentado aquí porque son matrices unitarias. Por ejemplo, se cumple que:

$$X^* X = I$$Pero al mismo tiempo, ¡las cuatro matrices son hermitianas o auto-adjuntas! Esto es:

$$X^* = X$$… y, por lo tanto:

$$XX = I$$Recuerde que los operadores hermitianos, en Mecánica Cuántica, se utilizan para definir observables, al ser operadores con valores propios siempre reales. Esto tiene alguna importancia en Física. Pero en Computación Cuántica casi siempre utilizaremos la base computacional, y la coincidencia no tiene mayor consecuencia. Sólo la menciono para no confundir el concepto de operador unitario con el de operador hermitiano.

Categorías
Quantum

Historia de un qubit

Veamos ahora cómo representar un qubit utilizando la menor redundancia posible. Sabemos que el estado de un qubit es, por definición, un vector en $\mathbb{C}^2$. Eso significa que necesitamos dos componentes complejos para representar su estado. Por motivos histórico, se suele usar la letra $\psi$ para denotar el estado cuántico. Por lo tanto y resumiendo:
$$\displaylines{
\vert\psi\rangle = \alpha\vert 0\rangle + \beta\vert 1\rangle \cr
\alpha, \beta \in \mathbb{C}}
$$A los componentes complejos del estado cuántico se les suele llamar amplitudes, como recuerdo de que en la metodología original de Schrödinger, $\psi$ se conoce como «función de onda».

Observad que, sin decir nada, he hecho que los componentes vayan automáticamente asociados a la base formada por los vectores $\vert 0\rangle$ y $\vert 1\rangle$. Los componentes de un vector siempre tienen que referirse a una base concreta, y la base universalmente utilizada en Quantum Computing es la base anterior, conocida como «base computacional». Estos dos vectores son especiales porque son los vectores propios, o eigenvectors, del operador de medición que transforma un qubit en un bit clásico.

Ahora bien, ya he mencionado que, en realidad, el estado cuántico es un rayo dentro del espacio vectorial complejo. Es decir, todos los estados que apuntan en la misma dirección son físicamente equivalentes. Por convenio, y para simplificar la forma de calcular, nos interesa elegir el vector dentro del rayo cuya longitud es la unidad. Por lo tanto, vamos a exigir que nuestras amplitudes cumplan con esta condición adicional:
$$\alpha\alpha^* +\beta\beta^*=1$$No olvidemos, ni por un momento, que tanto $\alpha$ como $\beta$ son números complejos. Por este motivo, no elevamos simplemente al cuadrado cada amplitud, sino que tenemos que multiplicarlas por sus valores conjugados, si queremos que la suma de los cuadrados de las amplitudes sea un valor real.

Notación matricial

Hay programadores ninjas que se buscan directamente la vida con la notación de Dirac y pueden hacer todos los cálculos necesarios directamente con ella. La mayoría de los mortales, en cambio, lo vemos todo más claro si utilizamos la notación matricial de toda la vida. ¡Ojo!, de momento estamos tratando con vectores, y las matrices cuadradas propiamente dichas aparecerán en breve, pero en esta notación, un vector puede representarse como una matriz vertical, de dos filas y una columna:
$$\vert\psi\rangle = \pmatrix{\alpha \cr \beta}$$La ventaja de esta notación es que nos permite omitir la base: implícitamente, se refiere a la base computacional. Por lo tanto, es una notación mucho más compacta. Otra ventaja es que nos permite entender por otra vía de qué va todo ese rollo de los bras, los kets y eso del «espacio vectorial dual». Hemos visto, por ejemplo, que un ket es una matriz vertical. Resulta, entonces, que un bra es simplemente una matriz horizontal:
$$\langle\psi\vert = \pmatrix{\alpha^* & \beta^* }$$Por supuesto, para pasar del ket al bra hemos tenido que conjugar las amplitudes. Con esta notación, además, el producto escalar se calcula muy fácilmente:
$$\eqalign{
\langle\psi\vert\psi\rangle =& \pmatrix{\alpha \cr \beta} \times \pmatrix{\alpha^* & \beta^* } \cr
=& \pmatrix{ \alpha\alpha^* +\beta\beta^* }
}$$Aquí tenemos que recordar las reglas del álgebra lineal: si multiplicamos una matriz $m\times n$ por otra $n\times p$, el resultado es una matriz $m\times p$. Nosotros hemos multiplicado una matriz $1 \times 2$ por una $2\times 1$, por lo que el resultado es una matriz de $1\times 1$. Una sola celda, es decir, un escalar o valor numérico a secas.

El cero no es cero-cero

A esta parte le pongo un subtítulo separado porque es importante. ¿Cómo se representan, usando componentes, los dos vectores de la base computacional?
$$\vert 0\rangle \equiv \pmatrix{1 \cr 0}\quad\vert 1\rangle \equiv \pmatrix{0 \cr 1}$$Esto es extremadamente simple, pero la experiencia personal me dice que es fácil confundirse y creer que el vector del estado $\vert 0\rangle$ tiene sus componentes a cero.

De hecho, el vector $\pmatrix{0 \cr 0}$ no es un estado cuántico correcto.

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

Matrices de covarianza

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

Matrices semidefinidas positivas

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

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

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

Categorías
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

Lilith

… 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. Analicemos las líneas 27 y 28 del listado anterior. 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:

  • Utilizar System.Numerics.Vector, 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.
  • 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, 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.