Categorías
Quantum

Estado cuántico

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

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

Tres axiomas

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

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

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

Espacios de Hilbert (I)

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

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

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

Producto interior hermitiano (I)

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

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

La notación de Dirac (I)

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

Observables (II)

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

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

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

La ecuación de Schrödinger (III)

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

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

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

Categorías
C#

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. 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.