Diagonalización cuántica de Krylov

En esta lección sobre la diagonalización cuántica de Krylov (KQD) responderemos a lo siguiente:

¿Qué es, en general, el método de Krylov?
¿Por qué funciona el método de Krylov y en qué condiciones?
¿Qué papel desempeña la informática cuántica?

La parte cuántica de los cálculos se basa en gran medida en el trabajo de la Ref [1].

El siguiente vídeo ofrece una visión general de los métodos de Krylov en computación clásica, motiva su uso y explica cómo la computación cuántica puede desempeñar un papel en esa corriente de trabajo. El texto siguiente ofrece más detalles e implementa un método de Krylov tanto de forma clásica como utilizando un ordenador cuántico.

1. Introducción a los métodos de Krylov

Un método de subespacio de Krylov puede referirse a cualquiera de varios métodos creados alrededor de lo que se denomina el subespacio de Krylov. Una revisión completa de estos está más allá del alcance de esta lección, pero Ref [2-4] todos pueden dar sustancialmente más antecedentes. Aquí nos centraremos en qué es un subespacio de Krylov, cómo y por qué es útil para resolver problemas de valores propios y, por último, cómo puede implementarse en un ordenador cuántico.

Definición: Dada una matriz $N\times N$ simétrica, semidefinida positiva $A$ , el espacio de Krylov $\mathcal{K}^r$ de orden $r$ es el espacio abarcado por los vectores obtenidos multiplicando las potencias superiores de una matriz $A$ , hasta $r-1\leq N$ , por un vector de referencia $\vert v \rangle$ .

\mathcal{K}^r = \text{span}\left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}

Aunque los vectores anteriores abarcan lo que llamamos un subespacio de Krylov, no hay razón para pensar que serán ortogonales. A menudo se utiliza un proceso iterativo de ortonormalización similar a la ortogonalización de Gram-Schmidt. En este caso, el proceso es ligeramente distinto, ya que cada nuevo vector se hace ortogonal a los demás a medida que se genera. En este contexto se denomina iteración de Arnoldi. Partiendo del vector inicial $|v\rangle$ , se genera el siguiente vector $A|v\rangle$ , y luego se asegura que este segundo vector es ortogonal al primero restando su proyección sobre $|v\rangle$ . Es decir

\begin{aligned} |v_0\rangle &=\frac{|v\rangle}{\left|\left| |v\rangle \right|\right|}\\ |v_1\rangle &=\frac{A|v\rangle-\langle v_0|A|v\rangle |v_0\rangle}{\left|\left|A|v\rangle-\langle v_0|A|v\rangle |v_0\rangle \right|\right|} \end{aligned}

Ahora es fácil ver que $|v_0\rangle \perp |v_1\rangle,$ ya que

\langle v_0 | v_1\rangle=\frac{\langle v_0 | A|v\rangle-\langle v_0 |A|v\rangle\langle v_0|v_0\rangle}{\left|\left| A|v\rangle-\langle A|v\rangle|v_0\rangle |v_0\rangle \right|\right|}=0

Hacemos lo mismo con el siguiente vector, asegurándonos de que es ortogonal a los dos anteriores:

|v_2\rangle=\frac{A |v_1\rangle-\langle v_0|A |v_1\rangle |v_0\rangle-\langle v_1|A |v_1\rangle |v_1\rangle}{\left|\left| A |v_1\rangle-\langle v_0|A |v_1\rangle |v_0\rangle-\langle v_1|A |v_1\rangle |v_1\rangle\right|\right|}

Si repetimos este proceso para todos los vectores de $r$ , tendremos una base ortonormal completa para un espacio de Krylov. Nótese que el proceso de ortogonalización aquí dará cero una vez $r>m$ , ya que $m$ vector ortogonal necesariamente abarca todo el espacio. El proceso también dará cero si algún vector es un vector propio de $A$ , ya que todos los vectores posteriores serán múltiplos de ese vector.

1.1 Un ejemplo sencillo: Krylov a mano

Veamos paso a paso la generación de un subespacio de Krylov sobre una matriz trivialmente pequeña, para que podamos ver el proceso. Partimos de una matriz inicial $A$ que nos interesa:

A=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}

Para este pequeño ejemplo, podemos determinar los vectores y valores propios fácilmente incluso a mano. Mostramos aquí la solución numérica.

# One might use linalg.eigh here, but later matrices may not be Hermitian. So we use
# linalg.eig in this lesson.

import numpy as np

A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]])
eigenvalues, eigenvectors = np.linalg.eig(A)
print("The eigenvalues are ", eigenvalues)
print("The eigenvectors are ", eigenvectors)

Output:

The eigenvalues are  [2.58578644 4.         5.41421356]
The eigenvectors are  [[ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]
 [ 7.07106781e-01  1.37464400e-16 -7.07106781e-01]
 [ 5.00000000e-01  7.07106781e-01  5.00000000e-01]]

Los registramos aquí para su posterior comparación:

\begin{aligned} a_0&=2.59,&|0\rangle&=&\begin{pmatrix}1/2\\-\sqrt{2}/2\\1/2\end{pmatrix}\\ \\ a_1&=4,&|1\rangle&=&\begin{pmatrix}\sqrt{2}/2\\0\\-\sqrt{2}/2\end{pmatrix}\\ \\ a_2&=5.41,&|2\rangle&=&\begin{pmatrix}1/2\\\sqrt{2}/2\\1/2\end{pmatrix} \end{aligned}

Nos gustaría estudiar cómo funciona (o falla) este proceso a medida que aumentamos la dimensión de nuestro subespacio de Krylov, $r$ . Para ello, aplicaremos este proceso:

Generar un subespacio del espacio vectorial completo a partir de un vector elegido al azar $|v\rangle$ (llamarlo $|v_0\rangle$ si ya está normalizado, como arriba).
Proyecte la matriz completa $A$ en ese subespacio y encuentre los valores propios de esa matriz proyectada $\tilde{A}$ .
Aumente el tamaño del subespacio generando más vectores, asegurándose de que son ortonormales, mediante un proceso similar a la ortogonalización Gram-Schmidt.
Proyecte $A$ en el subespacio mayor y encuentre los valores propios de la matriz resultante, $\tilde{A}$ .
Repita esta operación hasta que los valores propios converjan (o, en este caso de juguete, hasta que haya generado vectores que abarquen todo el espacio vectorial de la matriz original $A$ ).

Una implementación normal del método de Krylov no necesitaría resolver el problema de valores propios para la matriz proyectada en cada subespacio de Krylov a medida que se construye. Se podría construir el subespacio de la dimensión deseada, proyectar la matriz sobre ese subespacio y diagonalizar la matriz proyectada. La proyección y diagonalización en cada dimensión del subespacio sólo se realiza para comprobar la convergencia.

$r=1$ de las dimensiones:

Elegimos un vector aleatorio, por ejemplo

|v_0\rangle=\begin{pmatrix}1\\0\\0\end{pmatrix}

Si aún no está normalizado, normalícelo.

Ahora proyectamos nuestra matriz $A$ en el subespacio de este único vector:

\tilde{A}_0=\langle v_0| A|v_0\rangle=\begin{pmatrix}1&0&0\end{pmatrix}\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}=(4)

Esta es nuestra proyección de la matriz sobre nuestro subespacio de Krylov cuando contiene un único vector, $|v_0\rangle$ . El valor propio de esta matriz es trivialmente 4. Podemos pensar en esto como nuestra estimación de orden cero de los valores propios (en este caso sólo uno) de $A$ . Aunque es una estimación pobre, es el orden de magnitud correcto.

$r=2$ de las dimensiones:

Ahora generamos el siguiente vector en nuestro subespacio mediante la operación con $A$ sobre el vector anterior:

A|v_0\rangle=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}=\begin{pmatrix}4\\-1\\0\end{pmatrix}

Ahora restamos la proyección de este vector sobre nuestro vector anterior para asegurar la ortogonalidad.

|v_1\rangle=A|v_0\rangle-\langle v_0 |A|v_0\rangle|v_0\rangle

|v_1\rangle=\begin{pmatrix}4\\-1\\0\end{pmatrix}-\begin{pmatrix}1& 0& 0\end{pmatrix}\begin{pmatrix}4\\-1\\0\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}=\begin{pmatrix}0\\-1\\0\end{pmatrix}

Si aún no está normalizado, normalícelo. En este caso, el vector ya estaba normalizado, por lo que

|v_1\rangle=\begin{pmatrix}0\\-1\\0\end{pmatrix}

Ahora proyectamos nuestra matriz A en el subespacio de estos dos vectores:

\tilde{A}_1= \begin{pmatrix} 1&0&0\\0&-1&0 \end{pmatrix} \begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1&0\\0&-1\\0&0\end{pmatrix}=\begin{pmatrix}1&0&0\\0&-1&0\end{pmatrix}\begin{pmatrix}4&1\\-1&-4\\0&1\end{pmatrix}=\begin{pmatrix}4&1\\1&4\end{pmatrix}

Nos queda el problema de determinar los valores propios de esta matriz. Pero esta matriz es ligeramente más pequeña que la matriz completa. En problemas que implican matrices muy grandes, trabajar con este subespacio más pequeño puede resultar muy ventajoso.

\det(\tilde{A_1}-\lambda I)=0

\begin{vmatrix} 4-\lambda&1\\1&4-\lambda\end{vmatrix} =(4-\lambda)^2-1=0

4-\lambda=±1→\lambda=3,5

Aunque sigue sin ser una buena estimación, es mejor que la estimación de orden cero. Realizaremos una iteración más para asegurarnos de que el proceso está claro. Sin embargo, esto debilita el objetivo del método, ya que acabaremos diagonalizando una matriz 3x3 en la siguiente iteración, lo que significa que no hemos ahorrado tiempo ni potencia de cálculo.

$r=3$ de las dimensiones:

Ahora generamos el siguiente vector en nuestro subespacio mediante la operación con A sobre el vector anterior:

A|v_1\rangle=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}0\\-1\\0\end{pmatrix}=\begin{pmatrix}1\\-4\\1\end{pmatrix}

Ahora restamos la proyección de este vector sobre nuestros dos vectores anteriores para asegurar la ortogonalidad.

\begin{aligned} |v_2\rangle&=A|v_1\rangle-\langle v_0 |A|v_1\rangle|v_0\rangle-\langle v_1 |A|v_1\rangle|v_1\rangle\\ |v_2\rangle&=\begin{pmatrix}1\\-4\\1\end{pmatrix}-\begin{pmatrix}1& 0& 0 \end{pmatrix}\begin{pmatrix}1\\-4\\1\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}-\begin{pmatrix}0&-1& 0\end{pmatrix}\begin{pmatrix}1\\-4\\1\end{pmatrix}\begin{pmatrix}0\\-1\\0\end{pmatrix}=\begin{pmatrix}0\\0\\1\end{pmatrix} \end{aligned}

Si aún no está normalizado, normalícelo. En este caso, el vector ya estaba normalizado, por lo que

|v_2 \rangle=\begin{pmatrix}0\\0\\1\end{pmatrix}

Ahora proyectamos nuestra matriz $A$ en el subespacio de estos vectores:

\tilde{A}_2=\begin{pmatrix}1&0&0\\0&-1&0\\0&0&1\end{pmatrix}\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1&0&0\\0&-1&0\\0&0&1\end{pmatrix}=\begin{pmatrix}4&-1&0\\1&-4&1\\0&-1&4\end{pmatrix}\begin{pmatrix}1&0&0\\0&-1&0\\0&0&1\end{pmatrix}=\begin{pmatrix}4&1&0\\1&4&1\\0&1&4\end{pmatrix}

Ahora determinamos los valores propios:

\det(\tilde{A}_2-\lambda I)=0

\begin{vmatrix}4-\lambda&1&0\\1&4-\lambda&1\\0&1&4-\lambda\end{vmatrix} = (4-\lambda)((4-\lambda)^2-1)-(4-\lambda)=0\\

4-\lambda=0,4-\lambda=±2^{1/2}→\lambda=4-2^{1/2},4,4+2^{1/2}≈2.59,4,5.41

Estos valores propios son exactamente los valores propios de la matriz original $A$ . Este debe ser el caso, ya que hemos ampliado nuestro subespacio de Krylov para que abarque todo el espacio vectorial de la matriz original $A$ .

En este ejemplo, el método de Krylov puede no parecer particularmente más fácil que la diagonalización directa. De hecho, como veremos en secciones posteriores, el método de Krylov sólo es ventajoso a partir de una determinada dimensión de matriz; con ello se pretende ayudarnos a resolver problemas de valores propios/vectores propios de matrices extremadamente grandes.

Imagen que muestra una matriz muy grande proyectada en un subespacio de Krylov, es decir, filas de vectores de Krylov formando una matriz a la izquierda, un Hamiltoniano, y luego columnas de vectores de Krylov a la derecha.

Este es el único ejemplo que mostraremos trabajado "a mano", pero en la sección 2 se muestran ejemplos computacionales.

Aclaración de términos

Un error común es creer que existe un único subespacio de Krylov para un problema determinado. Pero claro, como hay muchos vectores iniciales a los que se podría aplicar nuestra matriz, hay muchos subespacios de Krylov posibles. Sólo utilizaremos la expresión "el subespacio de Krylov" para referirnos a un subespacio de Krylov específico ya definido para un ejemplo concreto. Para los planteamientos generales de resolución de problemas nos referiremos a "un subespacio de Krylov". Una última aclaración es que es válido referirse a un " espacio de Krylov". A menudo se le denomina " subespacio de Krylov" por su uso en el contexto de la proyección de matrices de un espacio inicial a un subespacio. En consonancia con ese contexto, aquí nos referiremos principalmente a él como subespacio.

Comprueba tu comprensión

Explique por qué no es (a) útil, y (b) posible extender la dimensión del subespacio de Krylov $r$ más allá de la dimensión $N$ de la matriz de interés.

(a) Dado que estamos ortonormalizando los vectores a medida que los generamos, un conjunto de $N$ tales vectores formará una base completa, lo que significa que una combinación lineal de ellos puede utilizarse para crear cualquier vector del espacio.
(b) El proceso de ortogonalización consiste en restar la proyección de un nuevo vector sobre todos los vectores anteriores. Si todos los vectores anteriores generan el espacio vectorial completo, al restar las proyecciones en el subespacio completo siempre obtendremos un vector nulo.

Imaginemos que un compañero investigador está mostrando cómo se aplica el método de Krylov a una pequeña matriz de ejemplo. ¿Hay algún problema con la elección de la matriz $A$ y el vector inicial $|\psi\rangle$ ?

A=\begin{pmatrix}2&1&3\\1&2&3\\3&3&5\end{pmatrix}

|\psi\rangle=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\-1\\0\end{pmatrix}.

Su colega ha elegido accidentalmente un vector propio para su vector inicial. Actuar con la matriz sobre el vector inicial simplemente devolverá el mismo vector, escalado por el valor propio. Esto no generará un subespacio de dimensión creciente. Aconseje a su colega que seleccione un vector inicial diferente, asegurándose de que no sea un vector propio.

Aplica el método de Krylov a la matriz dada, seleccionando un nuevo vector inicial adecuado. Anota las estimaciones del valor propio mínimo en los órdenes 0 y 1 de tu subespacio de Krylov.

A=\begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix}

Hay muchas respuestas posibles en función de la elección del vector inicial. Vamos a elegir:
$|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}.$
Para obtener $|v_1\rangle$ aplicamos $A$ una vez a $|v_0\rangle$ , y luego hacemos $|v_1\rangle$ ortogonal a $|v_0\rangle.$
$A|v_0\rangle=\begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}$ $A|v_0\rangle - \langle v_0|A|v_0\rangle |v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix} - \frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}-\frac{7}{3}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}=\sqrt{\frac{3}{2}}\begin{pmatrix}-1/3\\2/3\\-1/3\end{pmatrix}$
En orden 0, la proyección sobre nuestro subespacio de Krylov es
$\langle v_0|A|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{7}{3}$
En 1er orden, la proyección sobre este subespacio de Krylov es
$\langle V^1|A|V^1\rangle=\begin{pmatrix}\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\-\sqrt{\frac{1}{6}}&\sqrt{\frac{2}{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \begin{pmatrix}\frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\\\frac{1}{\sqrt{3}}& \sqrt{\frac{2}{3}} \\ \frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix}$
Esto se puede hacer a mano, pero es más fácil con numpy:
```
import numpy as np
vstar = np.array([[1/np.sqrt(3),1/np.sqrt(3),1/np.sqrt(3)],[-1/np.sqrt(6),np.sqrt(2/3),-1/np.sqrt(6)]]
)
A = np.array([[1, 1, 0],
              [1, 1, 1],
              [0, 1, 1]])
v = np.array([[1/np.sqrt(3),-1/np.sqrt(6)],[1/np.sqrt(3),np.sqrt(2/3)],[1/np.sqrt(3),-1/np.sqrt(6)]])
proj = vstar@A@v
print(proj)
eigenvalues, eigenvectors = np.linalg.eig(proj)
print("The eigenvalues are ", eigenvalues)
print("The eigenvectors are ", eigenvectors)
```
outputs:
```
[[ 2.33333333  0.47140452]
 [ 0.47140452 -0.33333333]]
The eigenvalues are  [ 2.41421356 -0.41421356]
The eigenvectors are  [[ 0.98559856 -0.16910198]
 [ 0.16910198  0.98559856]]
```
La estimación del valor propio mínimo es -0.414.

1.2 Tipos de métodos de Krylov

los "métodos de subespacios de Krylov" pueden referirse a cualquiera de las diversas técnicas iterativas utilizadas para resolver grandes sistemas lineales y problemas de valores propios. Todas ellas tienen en común que construyen una solución aproximada a partir de un subespacio de Krylov

$\mathcal{K}^r(A,|v\rangle ) = \text{span}\{|v\rangle, A|v\rangle, A^2|v\rangle, ..., A^{r-1}|v\rangle\},$

donde $|v\rangle$ es la conjetura inicial (ver Ref [5] ). Se diferencian en cómo eligen la mejor aproximación de este subespacio, equilibrando factores como la velocidad de convergencia, el uso de memoria y el coste computacional global. El objetivo de esta lección es aprovechar la computación cuántica en el contexto de los métodos de subespacios de Krylov; una discusión exhaustiva de estos métodos queda fuera de su alcance. Las breves definiciones que figuran a continuación sólo sirven de contexto e incluyen algunas referencias para investigar más a fondo estos métodos.

El método del gradiente conjugado (CG ): Este método se utiliza para resolver sistemas lineales simétricos y definidos positivamente [6]. Minimiza la norma A del error en cada iteración, por lo que es particularmente eficaz para sistemas derivados de EDP elípticas discretizadas [7]. Utilizaremos este enfoque en la siguiente sección para explicar por qué un subespacio de Krylov sería un subespacio eficaz para buscar soluciones mejoradas a los sistemas lineales.

El método del residuo mínimo generalizado (GMRES ): Está diseñado para resolver sistemas lineales generales no simétricos. Minimiza la norma residual sobre un espacio de Krylov en cada iteración, lo que lo hace robusto pero potencialmente intensivo en memoria para sistemas grandes [7].

El método del residuo mínimo (MINRES ): Este método se utiliza para resolver sistemas lineales indefinidos simétricos. Es similar al GMRES pero aprovecha la simetría matricial para reducir el coste computacional [8].

Otros enfoques destacables son el método de ortogonalización completa (FOM), estrechamente relacionado con el método de Arnoldi para problemas de valores propios, el método de gradiente bi-conjugado ( BiCG ) y el método de reducción de dimensión inducida (IDR).

1.3 ¿Por qué funciona el método del subespacio de Krylov?

Aquí motivaremos que el método del subespacio de Krylov debería ser una forma eficiente de aproximar los valores propios de la matriz mediante el refinamiento iterativo de las aproximaciones del vector propio, a través de la lente del descenso más pronunciado. Argumentaremos que dada una conjetura inicial de un estado base, el espacio de correcciones sucesivas a esa conjetura inicial que produce la convergencia más rápida es un subespacio de Krylov. No llegamos a una prueba rigurosa del comportamiento de convergencia.

Supongamos que nuestra matriz de interés $A$ es simétrica y definida positiva. Esto hace que nuestro argumento sea más relevante para el método CG anterior. Aquí no hacemos suposiciones sobre la dispersión; tampoco afirmamos que $A$ deba ser hermitiano (lo que tiene que ser si es un hamiltoniano).

Típicamente deseamos resolver un problema de la forma

A|x\rangle=|b\rangle.

Uno podría imaginar que $|b\rangle=c|x\rangle$ donde $c$ es alguna constante, como en un problema de valor propio. Pero el enunciado de nuestro problema sigue siendo más general por ahora.

Partimos de un vector $|x_0\rangle$ que es una solución aproximada. Aunque existen paralelismos entre esta conjetura $|x_0\rangle$ y $|v_0\rangle$ en la sección 1.1, aquí no los aprovechamos. Nuestra conjetura $|x_0\rangle$ tiene error, que llamamos $|e_0\rangle:$

|e_0\rangle:=|x\rangle−|x_0\rangle.

También definimos el residuo $R_0:$

|R_0\rangle=|b\rangle−A|x_0\rangle.

Aquí utilizamos la mayúscula $R$ para distinguir el residuo de la dimensión de nuestro subespacio de Krylov $r$ .

Un vector propio verdadero etiquetado como x, una estimación etiquetada como x₀ y una representación gráfica del error entre ambos.

Ahora queremos hacer un paso de corrección de la forma

|x_1\rangle=|x_0\rangle+|p_0\rangle,

lo que esperamos que mejore nuestra aproximación. Aquí $|p_0\rangle$ es algún vector aún por determinar. Sea $|e_1\rangle$ el error después de la corrección. Entonces

|e_1\rangle=|x\rangle−|x_1\rangle=|x\rangle−(|x_0\rangle+|p_0\rangle)=|e_0\rangle−|p_0\rangle.

Un vector propio verdadero y una actualización de la estimación inicial. La estimación actualizada se acerca más al vector propio real.

Nos interesa saber cómo se comporta nuestro error cuando es transformado por nuestra matriz. Así pues, calculemos la $A$ -norma del error. Es decir

\begin{aligned} ∥|e_0\rangle−|p_0\rangle∥_A^2&=\left(\langle e_0|A−\langle p_0|A\right)\left(|e_0\rangle−|p_0\rangle\right)\\ & = \langle e_0|A|e_0 \rangle − \langle e_0|A|p_0\rangle − \langle p_0|A|e_0\rangle+\langle p_0|A|p_0\rangle\\ & = \langle e_0|A|e_0\rangle−2\langle e_0|A|p_0\rangle+\langle p_0|A|p_0\rangle\\ & = d−2\langle R_0|p_0\rangle +\langle p_0|A|p_0\rangle, \end{aligned}

donde hemos utilizado la simetría de $A$ y también que $A |e_0\rangle = |R_0\rangle.$ Aquí $d$ es alguna constante independiente de $|p_0\rangle$ . Como se mencionó en la Sección 1.2, la $A$ -norma del error no es la única cantidad que podríamos elegir para minimizar, pero es una buena. Queremos ver cómo varía esta cantidad con nuestra elección de vectores de corrección $|p_0\rangle.$ Así que definimos la función $f$ estableciendo

f(|p_0\rangle)=\langle p_0|A|p_0\rangle−2\langle R_0|p_0\rangle+d.

$f$ no es más que el error $|e_1\rangle$ en función de la corrección $|p_0\rangle$ medida en la $A$ -norma. Por lo tanto, queremos elegir $|p_0\rangle$ de forma que $f(|p_0\rangle)$ sea lo más pequeño posible. Para ello, calculamos el gradiente de $f$ . Utilizando la simetría de $A$ tenemos

\nabla f(|p_0\rangle) = 2(A|p_0\rangle−|R_0\rangle).

El gradiente apunta en la dirección del ascenso más pronunciado, lo que significa que su opuesto nos da la dirección en la que la función disminuye más: la dirección del descenso más pronunciado. En nuestra conjetura inicial $|x_0\rangle$ , donde $|p_0\rangle=0$ , tenemos que $\nabla f(0) = -2|R_0\rangle.$ Por lo tanto, la función $f$ es la que más decrece en la dirección del residuo $|R_0\rangle.$ Así que nuestra elección inicial se beneficiaría más de la adición del vector $|p_0\rangle=\alpha_0 |R_0\rangle$ por algún escalar $\alpha_0$ .

En el siguiente paso, elegimos, de nuevo, un vector $|p_1\rangle$ y añadimos su valor a la aproximación actual. Usando el mismo argumento que antes elegimos $|p_1\rangle = \alpha_1 |R_1\rangle$ para algún escalar $\alpha_1$ . Continuamos de esta manera, de forma que la iteración $k^\text{th}$ de nuestro vector es

|x_{k+1}\rangle=|x_0\rangle+\alpha_0 |R_0\rangle+\alpha_1 |R_1\rangle+⋯+\alpha_k |R_k\rangle.

Equivalentemente, queremos construir el espacio del que elegimos nuestras estimaciones mejoradas añadiendo $|R_0\rangle$ , $|R_1\rangle$ , y así sucesivamente, en orden. El vector estimado $k^\text{th}$ se encuentra en

|x_{k+1}\rangle\in |x_0\rangle+\text{span}\{|R_0\rangle,|R_1\rangle,…,|R_k\rangle \}.

Ahora, utilizando la relación que

|R_{k+1}\rangle=|b\rangle−A |x_{k+1}\rangle=|b\rangle−A(|x_k\rangle+\alpha_k |R_k\rangle)=|R_k\rangle−\alpha_k A |R_k\rangle,

vemos que

\text{span} \{|R_0\rangle,|R_1\rangle,…,|R_k\rangle \}=\text{span} \{|R_0\rangle,A|R_0\rangle,…,A^{k}|R_0\rangle \}.

Es decir, el espacio que construimos que se aproxima más eficientemente a la solución correcta $|x\rangle$ es exactamente el espacio construido por la operación sucesiva de la matriz $A$ en $|R_0\rangle.$ . Un subespacio de Krylov es el espacio abarcado por los vectores de las direcciones sucesivas del descenso más pronunciado.

Por último, reiteramos que no hemos hecho afirmaciones numéricas sobre el escalado de este enfoque, ni hemos discutido el beneficio comparativo para matrices dispersas. Esto sólo pretende motivar el uso de los métodos de subespacios de Krylov, y añadir algo de sentido intuitivo para ellos. A continuación exploraremos numéricamente el comportamiento de estos métodos.

Comprueba tu comprensión

En el flujo de trabajo anterior, propusimos minimizar la $A$ -norma del error. ¿Qué otras cantidades se podrían minimizar al buscar el estado fundamental y su valor propio?

Se podría imaginar utilizar el vector residual en lugar de la $A$ -norma del error. Puede haber casos en los que sea útil considerar el propio vector de error.

2. Métodos de Krylov en la computación clásica

En esta sección implementamos computacionalmente las iteraciones de Arnoldi para poder aprovechar un subespacio de Krylov en la resolución de problemas de valores propios. Primero lo aplicaremos a un ejemplo a pequeña escala y, a continuación, examinaremos cómo se escala el tiempo de cálculo a medida que aumenta el tamaño de la matriz de interés. Una idea clave aquí será que la generación de los vectores que abarcan el espacio de Krylov contribuirá en gran medida al tiempo total de cálculo necesario. La memoria necesaria varía según el método de Krylov. Pero las restricciones de memoria pueden limitar el uso de los métodos tradicionales de Krylov.

2.1 Ejemplo sencillo a pequeña escala

En el proceso de creación de un subespacio de Krylov, necesitaremos ortonormalizar los vectores de nuestro subespacio. Definamos una función que tome un vector establecido de nuestro subespacio vknown (que no se supone normalizado) y un vector candidato para añadirlo a nuestro subespacio vnext y hacer que vnext sea ortogonal a vknown y normalizado. Definamos además una función que recorra este proceso para todos los vectores establecidos en nuestro subespacio de Krylov para garantizar un conjunto totalmente ortonormal.

# vknown is some established vector in our subspace. vnext is one we wish to add,
# which must be orthogonal to vknown.


def orthog_pair(vknown, vnext):
    vknown = vknown / np.sqrt(vknown.T @ vknown)
    diffvec = vknown.T @ vnext * vknown
    vnext = vnext - diffvec
    return vnext


# v is the candidate vector to be added to our subspace. s is the existing subspace.


def orthoset(v, s):
    v = v / np.sqrt(v.T @ v)
    temp = v
    for i in range(len(s)):
        temp = orthog_pair(s[i], temp)
    v = temp / np.sqrt(temp.T @ temp)
    return v

Definamos ahora una función que construya iterativamente un subespacio de Krylov cada vez mayor, hasta que el espacio de vectores de Krylov abarque el espacio completo de la matriz original. Esto nos permitirá ver hasta qué punto los valores propios obtenidos mediante nuestro método del subespacio de Krylov coinciden con los valores exactos, en función de la dimensión del subespacio de Krylov. Es importante destacar que nuestra función krylov_full_build devuelve los vectores de Krylov, los hamiltonianos proyectados, los valores propios y el tiempo necesario.

# Necessary imports and definitions to track time in microseconds
import time


def time_mus():
    return int(time.time() * 1000000)


# This function constructs a Krylov subspace that spans the whole space of the original matrix.
#     Input:
#       v0          : initial vector
#       matrix      : original matrix to be diagonalized
#     Output:
#       ks          : Krylov vectors
#       Hs          : projected Hamiltonians
#       eigs        : eigenvalues
#       k_tot_times : time required for the operation


def krylov_full_build(v0, matrix):
    t0 = time_mus()
    b = v0 / np.sqrt(v0 @ v0.T)
    A = matrix
    ks = []
    ks.append(b)
    Hs = []
    eigs = []
    Hs.append(b.T @ A @ b)
    eigs.append(np.array([b.T @ A @ b]))
    k_tot_times = []

    for j in range(len(A) - 1):
        vec = A @ ks[j].T
        ortho = orthoset(vec, ks)
        ks.append(ortho)
        ksarray = np.array(ks)
        Hs.append(ksarray @ A @ ksarray.T)
        eigs.append(np.linalg.eig(Hs[j + 1]).eigenvalues)
        k_tot_times.append(time_mus() - t0)

    # Return the Krylov vectors, the projected Hamiltonians, the eigenvalues,
    # and the total time required.
    return (ks, Hs, eigs, k_tot_times)

Probaremos esto en una matriz que sigue siendo bastante pequeña, pero más grande de lo que querríamos hacer a mano.

# Define our small test matrix
test_matrix = np.array(
    [
        [4, -1, 0, 1, 0],
        [-1, 4, -1, 2, 1],
        [0, -1, 4, 3, 3],
        [1, 2, 3, 4, 0],
        [0, 1, 3, 0, 4],
    ]
)

# Give the test matrix and an initial guess as arguments in the function defined above.
# Calculate outputs.
test_ks, test_Hs, test_eigs, text_k_tot_times = krylov_full_build(
    np.array([0.5, 0.5, 0, 0.5, 0.5]), test_matrix
)

Podemos comprobar nuestras funciones asegurándonos de que en el último paso (cuando el espacio de Krylov es el espacio vectorial completo de la matriz original) los valores propios del método de Krylov coinciden exactamente con los de la diagonalización numérica exacta:

print(np.linalg.eig(test_matrix).eigenvalues)
print(test_eigs[len(test_matrix) - 1])

Output:

[-1.36956923  8.43756009  2.9040308   5.34436028  4.68361806]
[-1.36956923  8.43756009  2.9040308   4.68361806  5.34436028]

Eso fue un éxito. Por supuesto, lo que realmente importa es lo buena que es nuestra aproximación en función de la dimensión de nuestro subespacio de Krylov. Dado que a menudo nos preocupamos por encontrar estados básicos y otros valores propios mínimos (y por otras razones más algebraicas que se explican más adelante), veamos nuestra estimación del valor propio más bajo en función de la dimensión del subespacio de Krylov. Es decir

def errors(matrix, krylov_eigs):
    targ_min = min(np.linalg.eig(matrix).eigenvalues)
    err = []
    for i in range(len(matrix)):
        err.append(min(krylov_eigs[i]) - targ_min)
    return err

import matplotlib.pyplot as plt

krylov_error = errors(test_matrix, test_eigs)

plt.plot(krylov_error)
plt.axhline(y=0, color="red", linestyle="--")  # Add dashed red line at y=0
plt.xlabel("Order of Krylov subspace")  # Add x-axis label
plt.ylabel("Error in minimum eigenvalue")  # Add y-axis label
plt.show()

Output:

Vemos que el valor propio mínimo se alcanza con bastante precisión una vez que el subespacio de Krylov ha crecido hasta $\mathcal{K}^2,$ y es perfecto por $\mathcal{K}^3.$

2.2 Escalado temporal con dimensión matricial

Convenzámonos de que el método de Krylov puede resultar ventajoso frente a los eigensolvers numéricos exactos de la siguiente manera:

Construir matrices aleatorias (no dispersas, no es la aplicación ideal para KQD)
Determine los valores propios utilizando dos métodos: directamente utilizando NumPy y utilizando un subespacio de Krylov.
Elegimos un límite para la precisión de nuestros valores propios antes de aceptar las estimaciones de Krylov.
Compara el tiempo de pared necesario para resolver de estas dos maneras.

Advertencias: Como discutiremos en detalle más adelante, la diagonalización cuántica de Krylov se aplica mejor a operadores cuyas representaciones matriciales son dispersas y/o pueden escribirse utilizando un pequeño número de grupos de operadores de Pauli conmutantes. Las matrices aleatorias que utilizamos aquí no se ajustan a esa descripción. Sólo son útiles para sondear la escala a la que los métodos clásicos de Krylov podrían ser útiles. En segundo lugar, al utilizar el método de Krylov calcularemos los valores propios utilizando muchos subespacios de Krylov de distintos tamaños. Informaremos del tiempo requerido para el subespacio de Krylov de dimensión mínima que alcanza nuestra precisión requerida para el valor propio del estado fundamental. De nuevo, esto es un poco diferente de resolver un problema que es intratable para eigensolvers exactos, ya que estamos utilizando la solución exacta para evaluar la dimensión necesaria.

Comenzamos generando nuestro conjunto de matrices aleatorias.

import numpy as np

# Set the random seed
np.random.seed(42)

# how many random matrices will we make
num_matrix = 200

matrices = []
for m in range(1, num_matrix):
    matrices.append(np.random.rand(m, m))

Ahora diagonalizamos cada matriz directamente, usando numpy. Calculamos el tiempo necesario para la diagonalización para su posterior comparación.

matrix_numpy_times = []
matrix_numpy_eigs = []
for mm in range(num_matrix - 1):
    t0 = time_mus()
    matrix_numpy_eigs.append(min(np.linalg.eig(matrices[mm]).eigenvalues))
    matrix_numpy_times.append(time_mus() - t0)

plt.plot(matrix_numpy_times)
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Obsérvese que, en la imagen anterior, el tiempo anómalamente elevado en torno a una dimensión de 125 puede deberse a la naturaleza aleatoria de las matrices o a la implementación en el procesador clásico utilizado, pero no es reproducible. Si se vuelve a ejecutar el código, se obtendrá un perfil diferente con distintos picos anómalos.

Ahora, para cada matriz, construiremos un subespacio de Krylov y calcularemos los valores propios por pasos. En cada paso, comprobaremos si el valor propio más bajo se ha obtenido dentro de nuestro error absoluto especificado. El subespacio que primero nos da valores propios dentro de nuestro error especificado es el subespacio para el que registraremos los tiempos de cálculo. La ejecución de esta célula puede tardar varios minutos, dependiendo de la velocidad del procesador. Puede omitir la evaluación o reducir la dimensión máxima de las matrices diagonalizadas. Basta con mirar los resultados precalculados.

# Choose the absolute error you can tolerate, and make a list for tracking the Krylov subspace size
# at which that error is achieved.
abserr = 0.05
accept_subspace_size = []

# Lists to store total time spent on the Krylov method, and the subset of that time spent on
# diagonalizing the projected matrix.
matrix_krylov_tot_times = []
matrix_krylov_dim = []

# Step through all our random matrices
for mm in range(0, num_matrix - 1):
    test_ks, test_Hs, test_eigs, test_k_tot_times = krylov_full_build(
        np.ones(len(matrices[mm])), matrices[mm]
    )
    # We have not yet found a Krylov subspace that produces our minimum eigenvalue to
    # within the required error.
    found = 0
    for j in range(0, len(matrices[mm]) - 1):
        # If we still haven't found the desired subspace...
        if found == 0:
            # ...but if this one satisfies the requirement, then record everything
            if (
                abs((min(test_eigs[j]) - matrix_numpy_eigs[mm]) / matrix_numpy_eigs[mm])
                < abserr
            ):
                accept_subspace_size.append(j)
                matrix_krylov_tot_times.append(test_k_tot_times[j])
                matrix_krylov_dim.append(mm)
                found = 1

Comparemos los tiempos obtenidos con estos dos métodos:

plt.plot(matrix_numpy_times, color="blue")
plt.plot(matrix_krylov_dim, matrix_krylov_tot_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Estos son los tiempos reales requeridos, pero a efectos de discusión, vamos a suavizar estas curvas promediando sobre unos pocos puntos adyacentes / dimensiones de la matriz. Esto se hace a continuación:

smooth_numpy_times = []
smooth_krylov_times = []

# Choose the number of adjacent points over which to average forward;
# the same will be used backward.
smooth_steps = 10

# We will do this smoothing for all points/matrix dimensions
for i in range(len(matrix_krylov_tot_times)):
    # Ensure we don't exceed the boundaries of our lists
    start = max(0, i - smooth_steps)
    end = min(len(matrix_krylov_tot_times) - 1, i + smooth_steps)

    # Dummy variables for accumulating an average over adjacent points. This is done for both Krylov
    # and the NumPy calculations.
    smooth_count = 0
    smooth_numpy_sum = 0
    smooth_krylov_sum = 0

    for j in range(start, end):
        smooth_numpy_sum = smooth_numpy_sum + matrix_numpy_times[j]
        smooth_krylov_sum = smooth_krylov_sum + matrix_krylov_tot_times[j]
        smooth_count = smooth_count + 1

    # Appending the averaged adjacent values to our new smooth lists
    smooth_numpy_times.append(smooth_numpy_sum / smooth_count)
    smooth_krylov_times.append(smooth_krylov_sum / smooth_count)

plt.plot(smooth_numpy_times, color="blue")
plt.plot(smooth_krylov_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (smoothed, microsec)")  # Add y-axis label
plt.show()

Output:

Nótese que el tiempo requerido para la construcción de un subespacio de Krylov inicialmente excede el tiempo requerido para la diagonalización completa de numpy. Pero a medida que aumenta el tamaño de la matriz, el método de Krylov se vuelve ventajoso. Esto es cierto incluso si reducimos nuestro error aceptable, pero la ventaja se establece en un tamaño de matriz mayor. Merece la pena analizarlo.

La complejidad temporal de la diagonalización numérica es $O(n^3)$ (con alguna variación entre algoritmos). La complejidad temporal de generar una base ortonormal de vectores $n$ también es $O(n^3)$ . Así que la ventaja del método de Krylov no está relacionada con el uso de $\it{some}$ base ortonormal, sino con el uso de una base ortonormal particular que efectivamente escoge los valores propios de interés. Ya hemos visto esto en el esbozo de una prueba en la primera sección de esta lección, y esto es crítico para las garantías de convergencia en los métodos de Krylov.

Repasemos nuestros progresos hasta ahora:

Para matrices muy grandes, el método del subespacio de Krylov puede producir valores propios aproximados dentro de las tolerancias requeridas más rápidamente que los algoritmos tradicionales de diagonalización.
Para matrices tan grandes, la generación de un subespacio de Krylov es la parte que más tiempo consume del método de subespacios de Krylov.
Por lo tanto, sería muy valiosa una forma eficaz de generar un subespacio de Krylov. Aquí es donde entra en escena el ordenador cuántico.

Comprueba tu comprensión

Consulte el gráfico suavizado de los tiempos de diagonalización en función de la dimensión de la matriz que se muestra más arriba.

(a) Según este gráfico, ¿a partir de qué dimensión aproximada de la matriz el método de Krylov empezó a ser más rápido?

b) ¿Qué aspectos del cálculo podrían modificar el orden de magnitud a partir del cual el método de Krylov se vuelve más rápido?

(a) Los resultados pueden variar si se vuelve a realizar el cálculo, pero el método de Krylov se vuelve más rápido a partir de una dimensión de aproximadamente 80-85.
(b) Hay muchas respuestas posibles. Algunos factores importantes son la precisión que necesitamos y la dispersión de las matrices que se van a diagonalizar.

3. Krylov mediante evolución temporal

Todo lo que hemos descrito hasta ahora puede hacerse de forma clásica. Entonces, ¿cómo y cuándo utilizaríamos un ordenador cuántico? Para matrices muy grandes, el método de Krylov puede requerir largos tiempos de cálculo y grandes cantidades de memoria. El tiempo necesario para la operación matricial de $H$ en $|v\rangle$ escala como $O(N^2)$ en el peor de los casos. Incluso la multiplicación de matrices dispersas en un vector (el caso típico de los solucionadores clásicos de tipo Krylov) tiene una complejidad temporal que escala como $O(N)$ . Esto se hace para cada vector que queramos en nuestro subespacio. La dimensión del subespacio $r$ no suele ser una fracción significativa de $N$ , y a menudo escala como $\log(N)$ . Así que generar todos los vectores se escala como $O(N^2 \log(N))$ en el peor de los casos. Aunque hay otros pasos, como la ortogonalización, éste es el escalado dominante que hay que tener en cuenta.

La computación cuántica nos permite cambiar qué atributos del problema determinan el escalado del tiempo y los recursos necesarios. En lugar de depender del tamaño de la matriz $N$ de forma generalizada, veremos cosas como el número de disparos y el número de términos de Pauli no conmutativos que componen el Hamiltoniano. Veamos cómo funciona.

3.1 Evolución temporal

Recordemos que el operador que evoluciona en el tiempo un estado cuántico es $e^{-iHt/\hbar}$ (y es muy común, especialmente en computación cuántica, eliminar el $\hbar$ de la notación). Una forma de entender e incluso realizar dicha función exponencial de un operador es observar su expansión en serie de Taylor. Nótese que esta operación actuando sobre algún vector inicial $|v\rangle$ produce una suma de términos con potencias crecientes de $H$ aplicados al estado inicial. Parece que podemos crear nuestro subespacio de Krylov evolucionando en el tiempo nuestro estado inicial

\begin{aligned} e^{-iHt/\hbar}→e^{-iHt}&≈1-iHt-\frac{(H^2 t^2)}{2}+⋯\\ e^{-iHt} |v\rangle &≈ |v\rangle-iHt|v\rangle-\frac{(H^2 t^2)}{2}|v\rangle+⋯ \end{aligned}

La salvedad está en realizar la evolución temporal en un ordenador cuántico real. Muchos de los términos del Hamiltoniano no conmutarán entre sí. Así, mientras que algunos operadores exponenciales simples como $e^{-iZ}$ corresponden a circuitos simples, los hamiltonianos generales no. Y como contienen términos no conmutativos, no podemos descomponer simplemente la exponencial en un producto de simples, como podemos hacer con los números.

e^{-iHt}=e^{-i(H_1+H_2+⋯+H_n)t}\neq e^{-iH_1 t} e^{-iH_2 t}... e^{-iH_n t}

Así que esto no es trivial, pero es un proceso bien estudiado en la computación cuántica. Llevamos a cabo la evolución temporal en ordenadores cuánticos mediante un proceso llamado trotterización, que en sí mismo es un tema rico [10]. Pero a un nivel muy alto, rompiendo la evolución temporal en pasos muy pequeños, digamos $m$ pasos de tamaño $dt$ , limitamos los efectos de la no conmutatividad de los términos.

e^{-iHt}=e^{-i(H_1+H_2+⋯+H_n )t} = (e^{-i(H_1+H_2+⋯+H_n )t/m} )^m ≈(e^{-iH_1 dt} e^{-iH_2 dt} …e^{-iH_n dt} )^m

donde $dt = t/m$ .

Llamemos "subespacio de Krylov de potencia" a un subespacio de Krylov de orden r que hayamos generado en el contexto clásico utilizando potencias de H directamente.

\mathcal{K}_P^r (H,|v\rangle)=\text{span}\{|v\rangle,H|v\rangle,H^2 |v\rangle… H^{r-1} |v\rangle\}

Ahora generamos un espacio similar utilizando el operador unitario de evolución temporal $U \equiv e^{-iHt}$ ; nos referiremos a éste como el "espacio unitario de Krylov" $\mathcal{K}_U^r$ . El subespacio de Krylov de potencia $\mathcal{K}_P^r$ que utilizamos clásicamente no puede generarse directamente en un ordenador cuántico, ya que $H$ no es un operador unitario. Se puede demostrar que el uso del subespacio unitario de Krylov ofrece garantías de convergencia similares a las del subespacio de Krylov de potencia, es decir, el error del estado fundamental converge eficientemente siempre y cuando el estado inicial $|v\rangle$ tenga un solapamiento con el estado fundamental verdadero que no sea exponencialmente evanescente, y siempre y cuando exista una brecha suficiente entre los valores propios. Véase la Ref [1] para un análisis más preciso de la convergencia.

Aquí, las potencias de $U$ se convierten en diferentes pasos de tiempo (la potencia $k^\text{th}$ de $U$ se adelanta un paso de tiempo $k \times dt$ ). Podemos etiquetar el elemento del subespacio que evoluciona en el tiempo para el tiempo total $k dt$ como $|\psi_k\rangle$ .

\begin{aligned} U&=e^{-iHdt}\\ U^k&=e^{-iH(kdt)}\\ \mathcal{K}_U^r&=\text{span}\{|\psi\rangle,U|\psi\rangle,U^2 |\psi\rangle… U^{r-1} |\psi\rangle\} \end{aligned}

Podemos proyectar nuestro Hamiltoniano H en el subespacio unitario de Krylov, $\mathcal{K}_U^r$ . En otras palabras, calculamos cada elemento matricial de $H$ en la base $\mathcal{K}_U^r$ . Denominaremos a esta matriz proyectada $\tilde{H}$ .

3.2 Cómo implementar en un ordenador cuántico

Los elementos de la matriz de $\tilde{H}$ vienen dados por los valores de expectativa $\langle \psi_m |H| \psi_n\rangle$ , que pueden estimarse utilizando el ordenador cuántico. Hay que tener en cuenta que $H$ puede escribirse como una suma de operadores de Pauli en diferentes qubits, y que no todos los operadores de Pauli pueden medirse simultáneamente. Podemos clasificar los términos de Pauli en grupos de términos conmutativos y medirlos todos a la vez. Pero puede que necesitemos muchos grupos de este tipo para abarcar todos los términos. Por lo tanto, el número de grupos de conmutación distintos en los que se pueden dividir los términos, $N_\text{GCP}$ , es importante.

H=\sum_{\alpha=1}^{N_\text{GCP}} c_\alpha P_\alpha

Aquí, $P_\alpha$ es una cadena de Pauli de la forma $P_\alpha \sim IZIXII...YZXIX$ o un conjunto de tales cadenas de Pauli que conmutan entre sí. Dado que podemos escribir « $H$ » como una suma de operadores medibles, las siguientes expresiones para los elementos de matriz de « $\tilde{H}$ » pueden obtenerse utilizando el estimador primitivo « IBM Quantum ».

\begin{aligned} \tilde{H}_{mn}&=\langle \psi_m |H| \psi_n\rangle\\ &=\langle \psi e^{iHt_m} |H| \psi e^{-iHt_n}\rangle\\ &=\langle \psi e^{iHmdt} |H|\psi e^{-iHndt}\rangle \end{aligned}

Donde $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ son los vectores del espacio unitario de Krylov y $t_n = n dt$ son los múltiplos del paso de tiempo $dt$ elegidos. En un ordenador cuántico, el cálculo de los elementos de cada matriz puede realizarse con cualquier algoritmo que permita obtener solapamientos entre estados cuánticos. En esta lección nos centraremos en la prueba de Hadamard. Dado que $\mathcal{K}_U$ tiene dimensión $r$ , el Hamiltoniano proyectado en el subespacio tendrá dimensiones $r \times r$ . Con $r$ suficientemente pequeño (generalmente $r<<100$ es suficiente para obtener la convergencia de las estimaciones de los valores propios) podemos entonces diagonalizar fácilmente el Hamiltoniano proyectado $\tilde{H},$ clásicamente. Sin embargo, no podemos diagonalizar directamente $\tilde{H}$ debido a la no ortogonalidad de los vectores del espacio de Krylov. Tendremos que medir sus solapamientos y construir una matriz $\tilde{S}$

\tilde{S}_{mn} = \langle \psi_m \vert \psi_n \rangle

Esto nos permite resolver el problema de valores propios en un espacio no ortogonal (también llamado problema de valores propios generalizado)

\tilde{H} \ \vec{c} = E \ \tilde{S} \ \vec{c}

A continuación, se pueden obtener estimaciones de los valores propios y los estados propios de $H$ observando las soluciones de este problema de valores propios generalizado. Por ejemplo, la estimación de la energía del estado fundamental se obtiene tomando el valor propio más pequeño $E$ y el estado fundamental del vector propio correspondiente $\vec{c}$ . Los coeficientes en $\vec{c}$ determinan la contribución de los distintos vectores que abarcan $\mathcal{K}_U$ .

Problema generalizado de valores propios

¿Por qué no podemos simplemente diagonalizar $\tilde{H}$ ? Dado que $\tilde{S}$ contiene la información sobre la geometría de la base de Krylov (que es no ortogonal en todos los casos excepto en casos muy especiales), $\tilde{H}$ por sí solo no describe una proyección del Hamiltoniano completo, por lo que sus valores propios no tienen ninguna relación particular con los del Hamiltoniano completo -- podrían ser cualquier valor aleatorio. La resolución del problema generalizado de valores propios es necesaria para obtener los valores propios y vectores propios aproximados correspondientes a la proyección del Hamiltoniano completo en el espacio de Krylov....

Un diagrama de circuito con muchas capas que indica que el circuito debe utilizarse muchas veces con diferentes estados para realizar la prueba de Hadamard modificada.

La figura muestra una representación en circuito de la prueba de Hadamard modificada, un método que se utiliza para calcular el solapamiento entre diferentes estados cuánticos. Para cada elemento de la matriz $\tilde{H}_{i,j}$ , se realiza una prueba Hadamard entre el estado $\vert \psi_i \rangle$ , $\vert \psi_j \rangle$ . Esto se destaca en la figura por el esquema de colores para los elementos de la matriz y las operaciones $\text{Prep} \; \psi_i$ , $\text{Prep} \; \psi_j$ correspondientes. Así, se requiere un conjunto de pruebas de Hadamard para todas las combinaciones posibles de vectores del espacio de Krylov para calcular todos los elementos de la matriz del Hamiltoniano proyectado $\tilde{H}$ . El hilo superior del circuito de pruebas de Hadamard es un qubit ancilla que se mide en la base X o Y, su valor de expectativa determina el valor del solapamiento entre los estados. El hilo inferior representa todos los qubits del Hamiltoniano del sistema. La operación $\text{Prep} \; \psi_i$ prepara el qubit del sistema en el estado $\vert \psi_i \rangle$ controlado por el estado del qubit ancilla (de forma similar para $\text{Prep} \; \psi_j$ ) y la operación $P$ representa la descomposición de Pauli del Hamiltoniano del sistema $H = \sum_i P_i$ . La implementación de esto en un ordenador cuántico se muestra con más detalle a continuación.

4. Diagonalización cuántica de Krylov en un ordenador cuántico

Ahora implementaremos la diagonalización cuántica de Krylov en un ordenador cuántico real. Empecemos por importar algunos paquetes útiles.

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import warnings

from qiskit.quantum_info import SparsePauliOp, Pauli
from qiskit.circuit import Parameter
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

# from qiskit.providers.fake_provider import Fake20QV1
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2 as Estimator, Batch

import itertools as it

warnings.filterwarnings("ignore")

Definimos la siguiente función para resolver el problema generalizado de valores propios que acabamos de explicar.

def solve_regularized_gen_eig(
    h: np.ndarray,
    s: np.ndarray,
    threshold: float,
    k: int = 1,
    return_dimn: bool = False,
) -> Union[float, List[float]]:
    """
    Method for solving the generalized eigenvalue problem with regularization

    Args:
        h (numpy.ndarray):
            The effective representation of the matrix in our Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of our Krylov subspace
        threshold (float):
            Cut-off value for the eigenvalue of s
        k (int):
            Number of eigenvalues to return
        return_dimn (bool):
            Whether to return the size of the regularized subspace

    Returns:
        lowest k-eigenvalue(s) that are the solution of the regularized generalized eigenvalue problem


    """
    s_vals, s_vecs = sp.linalg.eigh(s)
    s_vecs = s_vecs.T
    good_vecs = np.array([vec for val, vec in zip(s_vals, s_vecs) if val > threshold])
    h_reg = good_vecs.conj() @ h @ good_vecs.T
    s_reg = good_vecs.conj() @ s @ good_vecs.T
    if k == 1:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][0]
    else:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k]

Al menos en la evaluación comparativa inicial, es útil conocer una solución clásica exacta para comprobar el comportamiento de la convergencia. La siguiente función calcula la energía del estado fundamental de un Hamiltoniano, utilizando el Hamiltoniano y el número de qubits como argumentos.

def single_particle_gs(H_op, n_qubits):
    """
    Find the ground state of the single particle(excitation) sector
    """
    H_x = []
    for p, coeff in H_op.to_list():
        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))

    H_z = []
    for p, coeff in H_op.to_list():
        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))

    H_c = H_op.coeffs

    print("n_sys_qubits", n_qubits)

    n_exc = 1
    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))
    print("n_exc", n_exc, ", subspace dimension", sub_dimn)

    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)

    sparse_vecs = [
        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)
    ]  # list all of the possible sets of n_exc indices of 1s in n_exc-particle states

    m = 0
    for i, i_set in enumerate(sparse_vecs):
        for j, j_set in enumerate(sparse_vecs):
            m += 1

            if len(i_set.symmetric_difference(j_set)) <= 2:
                for p_x, p_z, coeff in zip(H_x, H_z, H_c):
                    if i_set.symmetric_difference(j_set) == p_x:
                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (
                            (-1) ** len(i_set.intersection(p_z))
                        )
                    else:
                        sgn = 0

                    few_particle_H[i, j] += sgn * coeff

    gs_en = min(np.linalg.eigvalsh(few_particle_H))
    print("single particle ground state energy: ", gs_en)
    return gs_en

4.1 Paso 1: Asignar el problema a circuitos y operadores cuánticos

Ahora definiremos un Hamiltoniano. Se diferencia de la función anterior en que la función anterior toma un Hamiltoniano como argumento y devuelve sólo el estado fundamental, y lo hace de forma clásica. Este Hamiltoniano que definimos aquí determina los niveles de energía de todos los eigenestados energéticos, y este Hamiltoniano puede construirse utilizando operadores de Pauli e implementarse en un ordenador cuántico.

Elegimos un Hamiltoniano correspondiente a una cadena de espines que puede tener cualquier orientación en el espacio, llamada "cadena de Heisenberg". Asumimos que el espín $i^\text{th}$ puede ser influenciado por sus vecinos más cercanos (los espines $(i-1)^\text{th}$ y $(i+1)^\text{th}$ ) pero no por vecinos más distantes. También tenemos en cuenta la posibilidad de que la interacción entre los espines sea diferente cuando los espines apuntan a lo largo de ejes diferentes. A veces, esta asimetría se debe, por ejemplo, a la estructura de la red cristalina en la que están incrustados los espines.

# Define problem Hamiltonian.
n_qubits = 10
# coupling strength for XX, YY, and ZZ interactions
JX = 1
JY = 3
JZ = 2

# Define the Hamiltonian:
H_int = [["I"] * n_qubits for _ in range(3 * (n_qubits - 1))]
for i in range(n_qubits - 1):
    H_int[i][i] = "Z"
    H_int[i][i + 1] = "Z"
for i in range(n_qubits - 1):
    H_int[n_qubits - 1 + i][i] = "X"
    H_int[n_qubits - 1 + i][i + 1] = "X"
for i in range(n_qubits - 1):
    H_int[2 * (n_qubits - 1) + i][i] = "Y"
    H_int[2 * (n_qubits - 1) + i][i + 1] = "Y"
H_int = ["".join(term) for term in H_int]
H_tot = [
    (term, JZ)
    if term.count("Z") == 2
    else (term, JY)
    if term.count("Y") == 2
    else (term, JX)
    for term in H_int
]

# Get operator
H_op = SparsePauliOp.from_list(H_tot)
print(H_tot)

Output:

[('ZZIIIIIIII', 2), ('IZZIIIIIII', 2), ('IIZZIIIIII', 2), ('IIIZZIIIII', 2), ('IIIIZZIIII', 2), ('IIIIIZZIII', 2), ('IIIIIIZZII', 2), ('IIIIIIIZZI', 2), ('IIIIIIIIZZ', 2), ('XXIIIIIIII', 1), ('IXXIIIIIII', 1), ('IIXXIIIIII', 1), ('IIIXXIIIII', 1), ('IIIIXXIIII', 1), ('IIIIIXXIII', 1), ('IIIIIIXXII', 1), ('IIIIIIIXXI', 1), ('IIIIIIIIXX', 1), ('YYIIIIIIII', 3), ('IYYIIIIIII', 3), ('IIYYIIIIII', 3), ('IIIYYIIIII', 3), ('IIIIYYIIII', 3), ('IIIIIYYIII', 3), ('IIIIIIYYII', 3), ('IIIIIIIYYI', 3), ('IIIIIIIIYY', 3)]

El código siguiente restringe el Hamiltoniano a estados de una sola partícula, y utiliza la norma espectral para establecer un buen tamaño para nuestro paso de tiempo $dt$ . Elegimos heurísticamente un valor para el paso temporal dt (basándonos en los límites superiores de la norma hamiltoniana). Ref [9] demostró que un paso de tiempo suficientemente pequeño es $\pi/\vert \vert H \vert \vert$ , y que es preferible hasta cierto punto subestimar este valor en lugar de sobreestimarlo, ya que la sobreestimación puede permitir que las contribuciones de los estados de alta energía corrompan incluso el estado óptimo en el espacio de Krylov. Por otro lado, elegir $dt$ demasiado pequeño conduce a un peor acondicionamiento del subespacio de Krylov, ya que los vectores base de Krylov difieren menos de un paso temporal a otro.

# Get Hamiltonian restricted to single-particle states
single_particle_H = np.zeros((n_qubits, n_qubits))
for i in range(n_qubits):
    for j in range(i + 1):
        for p, coeff in H_op.to_list():
            p_x = Pauli(p).x
            p_z = Pauli(p).z
            if all(p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)):
                sgn = ((-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))) * (
                    (-1) ** p_z[i]
                )
            else:
                sgn = 0
            single_particle_H[i, j] += sgn * coeff
for i in range(n_qubits):
    for j in range(i + 1, n_qubits):
        single_particle_H[i, j] = np.conj(single_particle_H[j, i])

# Set dt according to spectral norm
dt = np.pi / np.linalg.norm(single_particle_H, ord=2)
dt

Output:

np.float64(0.17453292519943295)

Especificamos el número de pasos de Trotter que se utilizarán en la evolución temporal. También especificamos una dimensión máxima de Krylov de 4. Esta dimensión de Krylov no es lo suficientemente grande para aplicaciones realistas. Pero es suficiente para este ejemplo. Además, comprobaremos la convergencia en dimensiones aún más pequeñas. En lecciones posteriores exploraremos métodos que nos permiten escalar y proyectar nuestros hamiltonianos en subespacios más grandes.

# Set parameters for quantum Krylov algorithm
krylov_dim = 4  # size of krylov subspace
num_trotter_steps = 4
dt_circ = dt / num_trotter_steps

Preparación del estado

Elija un estado de referencia $\vert \psi \rangle$ que tenga cierto solapamiento con el estado de tierra. Para este Hamiltoniano, Usamos el estado a con una excitación en el qubit del medio $\vert 00..010...00 \rangle$ como nuestro estado de referencia.

qc_state_prep = QuantumCircuit(n_qubits)
qc_state_prep.x(int(n_qubits / 2) + 1)
qc_state_prep.draw("mpl", scale=0.5)

Output:

Evolución temporal

Podemos realizar el operador de evolución temporal generado por un Hamiltoniano dado: $U=e^{-iHt}$ mediante la aproximación de Lie-Trotter. Para simplificar, utilizamos la dirección PauliEvolutionGate integrada en el circuito de evolución temporal. La sintaxis general es la siguiente

t = Parameter("t")

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)
)

qr = QuantumRegister(n_qubits)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

Output:

<qiskit.circuit.instructionset.InstructionSet at 0x7ccaa4664250>

Utilizaremos una versión de esto a continuación en la prueba de Hadamard, pero dando un paso adelante para los tiempos $dt$ .

Prueba de Hadamard

Recordemos que deseamos calcular los elementos matriciales tanto de $\tilde{H}$ como de la matriz de Gram $\tilde{S}$ utilizando la prueba de Hadamard. Repasemos cómo funciona en este contexto, centrándonos primero en la construcción de $\tilde{H}.$ . El proceso general se representa gráficamente a continuación. Las capas de bloques de colores de preparación de estados $\text{Prep}|\psi_i\rangle$ sirven para recordar que este proceso se lleva a cabo para todas las combinaciones de $|\psi_i\rangle$ y $|\psi_j\rangle$ en nuestro subespacio.

Imagen de un diagrama de circuito cuántico con muchas capas que indica que el circuito debe evaluarse para muchos estados diferentes con el fin de realizar la prueba de Hadamard.

Los estados del sistema en los pasos indicados son:

\begin{aligned} \text{Step 0:}\qquad|\Psi\rangle & = |0\rangle|0\rangle^N \\ \text{Step 1:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle + |1\rangle \Big)|0\rangle^N \\ \text{Step 2:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle|0\rangle^N+|1\rangle |\psi_i\rangle\Big)\\ \text{Step 3:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle |0\rangle^N+|1\rangle P |\psi_i\rangle\Big) \\ \text{Step 4:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle |\psi_j\rangle+|1\rangle P|\psi_i\rangle\Big) \end{aligned}

Aquí $P$ es un término de Pauli en la descomposición del Hamiltoniano (nótese que no puede ser una combinación lineal de múltiples términos de Pauli conmutativos ya que eso no sería unitario -- la agrupación es posible usando una construcción diferente que mostraremos más adelante) $\text{Prep} \; \psi_i$ , $\text{Prep} \; \psi_j$ son operaciones controladas que preparan $|\psi_i\rangle$ , $|\psi_j\rangle$ vectores del espacio unitario de Krylov, con $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$ . Aplicando medidas de $X$ y $Y$ a este circuito se calculan las partes real e imaginaria, respectivamente, de los elementos matriciales que necesitamos.

Empezando por el paso 4 anterior, aplique la puerta Hadamard $H$ al qubit zeroth.

\begin{equation*} |\Psi\rangle \longrightarrow\quad\frac{1}{2}|0\rangle\Big( |\psi_j\rangle + P|\psi_i\rangle\Big) + \frac{1}{2}|1\rangle\Big(|\psi_j\rangle - P|\psi_i\rangle\Big) \end{equation*}

A continuación, mida $X$ o $Y$ .

\begin{equation*} \begin{split} \Rightarrow\quad\langle X\rangle &= \frac{1}{4}\Bigg(\Big\|| \psi_j\rangle + P|\psi_i\rangle \Big\|^2-\Big\||\psi_j\rangle - P|\psi_i\rangle\Big\|^2\Bigg) \\ &= \text{Re}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{split} \end{equation*}

A partir de la identidad $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$ . Del mismo modo, midiendo $Y$ se obtiene

\begin{equation*} \langle Y\rangle = \text{Im}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{equation*}

Añadiendo estos pasos a la evolución temporal que establecimos anteriormente escribimos lo siguiente.

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)

## Create the time-evo op dagger circuit
evol_gate_d = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)
evol_gate_d = evol_gate_d.inverse()

# Put pieces together
qc_reg = QuantumRegister(n_qubits)
qc_temp = QuantumCircuit(qc_reg)
qc_temp.compose(qc_state_prep, inplace=True)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate, qargs=qc_reg)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate_d, qargs=qc_reg)
qc_temp.compose(qc_state_prep.inverse(), inplace=True)

# Create controlled version of the circuit
controlled_U = qc_temp.to_gate().control(1)

# Create hadamard test circuit for real part
qr = QuantumRegister(n_qubits + 1)
qc_real = QuantumCircuit(qr)
qc_real.h(0)
qc_real.append(controlled_U, list(range(n_qubits + 1)))
qc_real.h(0)

print("Circuit for calculating the real part of the overlap in S via Hadamard test")
qc_real.draw("mpl", fold=-1, scale=0.5)

Output:

Circuit for calculating the real part of the overlap in S via Hadamard test

Ya advertimos de la profundidad que suponen los circuitos Trotter. Realizar la prueba de Hadamard en estas condiciones puede dar lugar a un circuito aún más profundo, especialmente una vez que descomponemos a puertas nativas. Esto aumentará aún más si tenemos en cuenta la topología del dispositivo. Así que antes de utilizar cualquier tiempo en el ordenador cuántico, es una buena idea comprobar la profundidad de 2 qubits de nuestro circuito.

print(
    "Number of layers of 2Q operations",
    qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),
)

Output:

Number of layers of 2Q operations 14401

Un circuito de esta profundidad no puede devolver resultados utilizables en los modernos ordenadores cuánticos. Si queremos construir $\tilde{H}$ y $\tilde{S},$ necesitamos una forma mejor. Esta es la razón de la prueba de Hadamard eficiente que se presenta a continuación.

4. Paso 2. Optimizar circuitos y operadores para el hardware de destino

Prueba de Hadamard eficiente

Podemos optimizar los circuitos profundos para la prueba de Hadamard que hemos obtenido introduciendo algunas aproximaciones y basándonos en alguna suposición sobre el Hamiltoniano del modelo. Por ejemplo, considere el siguiente circuito para la prueba de Hadamard:

Imagen de un diagrama de circuito cuántico con muchas capas que indica que el circuito debe evaluarse para muchos operadores unitarios diferentes con el fin de realizar la prueba de Hadamard modificada y eficiente.

Supongamos que podemos calcular clásicamente $E_0$ , el valor propio de $|0\rangle^N$ bajo el Hamiltoniano $H$ . Esto se cumple cuando el Hamiltoniano preserva la simetría U(1). Aunque esto puede parecer una suposición fuerte, hay muchos casos en los que es seguro asumir que existe un estado de vacío (en este caso se mapea al estado $|0\rangle^N$ ) que no se ve afectado por la acción del Hamiltoniano. Este es el caso, por ejemplo, de los hamiltonianos químicos que describen moléculas estables (en las que se conserva el número de electrones). Dado que la puerta $\text{Prep} \; \psi_0$ , prepara el estado de referencia deseado $\ket{\psi_0} = \text{Prep} \; \psi_0 \ket{0} = e^{-i H 0 dt} U_{\psi_0} \ket{0}$ , por ejemplo, preparar el estado HF para la química $\text{Prep} \; \psi_0$ sería un producto de NOTs de un solo qubit, por lo que controlada- $\text{Prep} \; \psi_0$ es sólo un producto de CNOTs. Entonces el circuito anterior implementa el siguiente estado antes de la medición:

\begin{aligned} \text{Step 0:}\qquad|\Psi\rangle & = \ket{0} \ket{0}^{N}\\ \text{Step 1:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(\ket{0}\ket{0}^N+ \ket{1} \ket{0}^N\right)\\ \text{Step 2:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(|0\rangle|0\rangle^N+|1\rangle|\psi_0\rangle\right)\\ \text{Step 3:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0}\ket{0}^N+\ket{1} U\ket{\psi_0}\right)\\ \text{Step 4:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0} \ket{\psi_0}+\ket{1} U\ket{\psi_0}\right)\\ & = \frac{1}{2}\left(\ket{+}\left(e^{i\phi}\ket{\psi_0}+U\ket{\psi_0}\right)+\ket{-}\left(e^{i\phi}\ket{\psi_0}-U\ket{\psi_0}\right)\right)\\ & = \frac{1}{2}\left(\ket{+i}\left(e^{i\phi}\ket{\psi_0}-iU\ket{\psi_0}\right)+\ket{-i}\left(e^{i\phi}\ket{\psi_0}+iU\ket{\psi_0}\right)\right) \end{aligned}

donde hemos utilizado el desfase simulable clásico $U\ket{0}^N = e^{i\phi}\ket{0}^N$ del paso 2 al 3. Por lo tanto, los valores esperados son

\begin{aligned} \langle X\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi_0}+\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi}\ket{\psi_0}+U\ket{\psi_0}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi_0}-\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi}\ket{\psi_0}-U\ket{\psi_0}\right) \Big)\\ &=\text{Re}\left[e^{-i\phi}\bra{\psi_0}PU\ket{\psi_0}\right], \end{aligned}

\begin{aligned} \langle Y\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi_0}+i\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi_0}\ket{\psi_0}-iU\ket{\psi_0}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi_0}-i\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi}\ket{\psi_0}+iU\ket{\psi_0}\right) \Big)\\ &=\text{Im}\left[e^{-i\phi}\bra{\psi_0}PU\ket{\psi_0}\right]. \end{aligned}

Utilizando estos supuestos pudimos escribir los valores de las expectativas de los operadores de interés con menos operaciones controladas. De hecho, sólo necesitamos implementar la preparación controlada de estados $\text{Prep} \; \psi_0$ y no evoluciones controladas en el tiempo. Si reformulamos nuestro cálculo de la forma anterior, podremos reducir en gran medida la profundidad de los circuitos resultantes.

Obsérvese que, como ventaja adicional, dado que el operador de Pauli aparece ahora como una medida al final del circuito en lugar de como una puerta controlada en el medio, puede medirse junto con otros operadores de Pauli conmutantes como en la descomposición $H=\sum_{\alpha = 1}^{N_\text{GCP}}c_\alpha P_\alpha$ dada anteriormente.

Descomponer el operador de evolución temporal con la descomposición de Trotter

En lugar de implementar el operador de evolución temporal de forma exacta, podemos utilizar la descomposición de Trotter para implementar una aproximación del mismo. Repitiendo varias veces una determinada orden de descomposición de Trotter, obtenemos una mayor reducción del error introducido por la aproximación. A continuación, construimos directamente la implementación de Trotter de la manera más eficiente para el gráfico de interacción del hamiltoniano que estamos considerando (solo interacciones entre vecinos más cercanos). En la práctica, insertamos rotaciones de Pauli $R_{xx}$ , $R_{yy}$ , $R_{zz}$ con fuerzas de acoplamiento $J_x,$ $J_y,$ y $J_z$ y un ángulo parametrizado $t$ , que corresponden a la implementación aproximada de $e^{-i (J_x XX + J_y YY + J_z ZZ) t}$ . Dada la diferencia en la definición de las rotaciones de Pauli y la evolución temporal que estamos tratando de implementar, tendremos que utilizar el parámetro $2*dt$ para lograr una evolución temporal de $dt$ . Además, invertimos el orden de las operaciones para un número impar de repeticiones de los pasos de Trotter, lo que es funcionalmente equivalente pero permite sintetizar operaciones adyacentes en una única unidad $SU(2)$ . Esto da como resultado un circuito mucho menos profundo que el que se obtiene utilizando la funcionalidad PauliEvolutionGate() genérica.

t = Parameter("t")

# Create instruction for rotation about XX+YY-ZZ:
Rxyz_circ = QuantumCircuit(2)
Rxyz_circ.rxx(2 * JX * t, 0, 1)
Rxyz_circ.ryy(2 * JY * t, 0, 1)
Rxyz_circ.rzz(2 * JZ * t, 0, 1)
Rxyz_instr = Rxyz_circ.to_instruction(label="R J_x XX + J_y YY + J_z ZZ")

interaction_list = [
    [[i, i + 1] for i in range(0, n_qubits - 1, 2)],
    [[i, i + 1] for i in range(1, n_qubits - 1, 2)],
]  # linear chain

qr = QuantumRegister(n_qubits)
trotter_step_circ = QuantumCircuit(qr)
for i, color in enumerate(interaction_list):
    for interaction in color:
        trotter_step_circ.append(Rxyz_instr, interaction)
    if i < len(interaction_list) - 1:
        trotter_step_circ.barrier()
reverse_trotter_step_circ = trotter_step_circ.reverse_ops()

qc_evol = QuantumCircuit(qr)
for step in range(num_trotter_steps):
    if step % 2 == 0:
        qc_evol = qc_evol.compose(trotter_step_circ)
    else:
        qc_evol = qc_evol.compose(reverse_trotter_step_circ)

qc_evol.decompose().draw("mpl", fold=-1, scale=0.5)

Output:

Preparamos de nuevo un estado inicial para esta prueba Hadamard eficiente.

control = 0
excitation = int(n_qubits / 2) + 1
controlled_state_prep = QuantumCircuit(n_qubits + 1)
controlled_state_prep.cx(control, excitation)
controlled_state_prep.draw("mpl", fold=-1, scale=0.5)

Output:

Circuitos modelo para calcular elementos matriciales de $\tilde{S}$ y $\tilde{H}$ mediante la prueba de Hadamard

La única diferencia entre los circuitos utilizados en la prueba de Hadamard será la fase en el operador de evolución temporal y los observables medidos. Por lo tanto, podemos preparar un circuito plantilla que represente el circuito genérico para la prueba Hadamard, con marcadores de posición para las puertas que dependen del operador de evolución temporal.

# Parameters for the template circuits
parameters = []
for idx in range(1, krylov_dim):
    parameters.append(dt_circ * (idx))

# Create modified hadamard test circuit
qr = QuantumRegister(n_qubits + 1)
qc = QuantumCircuit(qr)
qc.h(0)
qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)
qc.barrier()
qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)
qc.barrier()
qc.x(0)
qc.compose(controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True)
qc.x(0)

qc.decompose().draw("mpl", fold=-1)

Output:

print(
    "The optimized circuit has 2Q gates depth: ",
    qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),
)

Output:

The optimized circuit has 2Q gates depth:  50

Esta profundidad se reduce sustancialmente en comparación con la prueba de Hadamard original. Esta profundidad es manejable en los ordenadores cuánticos modernos, aunque sigue siendo bastante elevada. Tendremos que recurrir a la mitigación de errores más avanzada para obtener resultados útiles.

Selecciona un backend en el que ejecutar nuestro cálculo cuántico de Krylov, de modo que podamos transpilar nuestro circuito para ejecutarlo en ese ordenador cuántico.

# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)

# Or you may choose a specify backend and channel if necessary for your workflow.
# service = QiskitRuntimeService(channel="ibm_quantum_platform")
# backend = service.backend("ibm_fez")

Ahora transpilemos nuestros circuitos y operadores.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

target = backend.target
basis_gates = list(target.operation_names)
pm = generate_preset_pass_manager(
    optimization_level=3, backend=backend, basis_gates=basis_gates
)

qc_trans = pm.run(qc)

print(qc_trans.depth(lambda x: x[0].num_qubits == 2))
print(qc_trans.count_ops())
qc_trans.draw("mpl", fold=-1, idle_wires=False, scale=0.5)

Output:

36
OrderedDict([('rz', 410), ('sx', 361), ('cz', 156), ('x', 18), ('barrier', 6)])

Tras la optimización, nuestra profundidad transpilada de dos qubits se reduce aún más.

4.3 Paso 3. Ejecutar utilizando una primitiva « IBM Quantum »

Ahora creamos PUBs para su ejecución con Estimator.

# Define observables to measure for S
observable_S_real = "I" * (n_qubits) + "X"
observable_S_imag = "I" * (n_qubits) + "Y"

observable_op_real = SparsePauliOp(
    observable_S_real
)  # define a sparse pauli operator for the observable
observable_op_imag = SparsePauliOp(observable_S_imag)

layout = qc_trans.layout  # get layout of transpiled circuit
observable_op_real = observable_op_real.apply_layout(
    layout
)  # apply physical layout to the observable
observable_op_imag = observable_op_imag.apply_layout(layout)
observable_S_real = (
    observable_op_real.paulis.to_labels()
)  # get the label of the physical observable
observable_S_imag = observable_op_imag.paulis.to_labels()

observables_S = [[observable_S_real], [observable_S_imag]]


# Define observables to measure for H
# Hamiltonian terms to measure
observable_list = []
for pauli, coeff in zip(H_op.paulis, H_op.coeffs):
    # print(pauli)
    observable_H_real = pauli[::-1].to_label() + "X"
    observable_H_imag = pauli[::-1].to_label() + "Y"
    observable_list.append([observable_H_real])
    observable_list.append([observable_H_imag])

layout = qc_trans.layout

observable_trans_list = []
for observable in observable_list:
    observable_op = SparsePauliOp(observable)
    observable_op = observable_op.apply_layout(layout)
    observable_trans_list.append([observable_op.paulis.to_labels()])

observables_H = observable_trans_list


# Define a sweep over parameter values
params = np.vstack(parameters).T


# Estimate the expectation value for all combinations of
# observables and parameter values, where the pub result will have
# shape (# observables, # parameter values).
pub = (qc_trans, observables_S + observables_H, params)

Los circuitos de $t=0$ se pueden calcular de forma clásica. Realizamos esta operación antes de pasar al caso $t\neq 0$ utilizando un ordenador cuántico.

from qiskit.quantum_info import StabilizerState, Pauli


qc_cliff = qc.assign_parameters({t: 0})


# Get expectation values from experiment
S_expval_real = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "X")
)
S_expval_imag = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "Y")
)

# Get expectation values
S_expval = S_expval_real + 1j * S_expval_imag

H_expval = 0
for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Get expectation values from experiment
    expval_real = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "X")
    )
    expval_imag = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "Y")
    )
    expval = expval_real + 1j * expval_imag

    # Fill-in matrix elements
    H_expval += coeff * expval


print(H_expval)

Output:

(10+0j)

Aunque pudimos reducir nuestra profundidad de compuerta en órdenes de magnitud utilizando la prueba Hadamard eficiente, la profundidad sigue siendo suficiente para requerir una mitigación de errores de última generación. A continuación, especificamos los atributos de la mitigación utilizada. Todos los métodos utilizados son importantes, pero merece la pena destacar específicamente la amplificación probabilística de errores (PEA). Esta potente técnica conlleva una gran sobrecarga cuántica. El cálculo realizado aquí puede tardar 20 minutos o más en ejecutarse en un ordenador cuántico real. Si lo desea, puede jugar con los parámetros siguientes para aumentar o disminuir la precisión y, en consecuencia, la sobrecarga. Los ajustes por defecto que se indican a continuación ofrecen resultados de alta fidelidad.

# Experiment options
num_randomizations = 300
num_randomizations_learning = 20
max_batch_circuits = 20
shots_per_randomization = 100
learning_pair_depths = [0, 4, 24]
noise_factors = [1, 1.3, 1.6]

# Base option formatting
options = {
    # Builtin resilience settings for ZNE
    "resilience": {
        "measure_mitigation": True,
        "zne_mitigation": True,
        "zne": {"noise_factors": noise_factors},
        # TREX noise learning configuration
        "measure_noise_learning": {
            "num_randomizations": num_randomizations_learning,
            "shots_per_randomization": shots_per_randomization,
        },
        # PEA noise model configuration
        "layer_noise_learning": {
            "max_layers_to_learn": 10,
            "layer_pair_depths": learning_pair_depths,
            "shots_per_randomization": shots_per_randomization,
            "num_randomizations": num_randomizations_learning,
        },
    },
    # Randomization configuration
    "twirling": {
        "num_randomizations": num_randomizations,
        "shots_per_randomization": shots_per_randomization,
        "strategy": "all",
    },
    # Experimental settings for PEA method
    "experimental": {
        # # Just in case, disable any further qiskit transpilation not related to twirling / DD
        # "skip_transpilation": True,
        # Execution configuration
        "execution": {
            "max_pubs_per_batch_job": max_batch_circuits,
            "fast_parametric_update": True,
        },
        # Error Mitigation configuration
        "resilience": {
            # ZNE Configuration
            "zne": {
                "amplifier": "pea",
                "return_all_extrapolated": True,
                "return_unextrapolated": True,
                "extrapolated_noise_factors": [0] + noise_factors,
            }
        },
    },
}

Por último, ejecutamos los circuitos para $\tilde{S}$ y $\tilde{H}$ con Estimator.

# This job required 17 minutes of QPU time to run on a Heron r2 processor. This is only an estimate.
# Your execution time may vary.

with Batch(backend=backend) as batch:
    # Estimator
    estimator = Estimator(mode=batch, options=options)

    job = estimator.run([pub], precision=1)

4.4 Paso 4. Procesar y analizar los resultados

Lo que hemos obtenido del ordenador cuántico son los elementos matriciales individuales de $\tilde{S}$ y los grupos conmutativos de Pauli que componen los elementos matriciales de $\tilde{H}$ . Estos términos deben combinarse para recuperar nuestras matrices, de modo que podamos resolver el problema generalizado de valores propios.

# Store the outputs as 'results'.
results = job.result()[0]

Calcular el hamiltoniano efectivo y las matrices de superposición

En primer lugar, calcule la fase acumulada por el estado $\vert 0 \rangle$ durante la evolución temporal no controlada

prefactors = [
    np.exp(-1j * sum([c for p, c in H_op.to_list() if "Z" in p]) * i * dt)
    for i in range(1, krylov_dim)
]

Una vez que tenemos los resultados de las ejecuciones de los circuitos podemos post-procesar los datos para calcular los elementos de la matriz de $S$

# Assemble S, the overlap matrix of dimension D:
S_first_row = np.zeros(krylov_dim, dtype=complex)
S_first_row[0] = 1 + 0j

# Add in ancilla-only measurements:
for i in range(krylov_dim - 1):
    # Get expectation values from experiment
    expval_real = results.data.evs[0][0][i]  # automatic extrapolated evs if ZNE is used
    expval_imag = results.data.evs[1][0][i]  # automatic extrapolated evs if ZNE is used

    # Get expectation values
    expval = expval_real + 1j * expval_imag
    S_first_row[i + 1] += prefactors[i] * expval

S_first_row_list = S_first_row.tolist()  # for saving purposes


S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        S_circ[j, i] = S_first_row[i - j]
    else:
        S_circ[j, i] = np.conj(S_first_row[j - i])

from sympy import Matrix

Matrix(S_circ)

Output:

\displaystyle \left[\begin{matrix}1.0 & 0.149322296177984 - 0.283023058106896 i & 0.185815978760175 - 0.0910521940394691 i & 0.0940509850777074 - 0.094154537369141 i\\0.149322296177984 + 0.283023058106896 i & 1.0 & 0.149322296177984 - 0.283023058106896 i & 0.185815978760175 - 0.0910521940394691 i\\0.185815978760175 + 0.0910521940394691 i & 0.149322296177984 + 0.283023058106896 i & 1.0 & 0.149322296177984 - 0.283023058106896 i\\0.0940509850777074 + 0.094154537369141 i & 0.185815978760175 + 0.0910521940394691 i & 0.149322296177984 + 0.283023058106896 i & 1.0\end{matrix}\right]

Y los elementos de la matriz de $\tilde{H}$

import itertools

# Assemble S, the overlap matrix of dimension D:
H_first_row = np.zeros(krylov_dim, dtype=complex)
H_first_row[0] = H_expval

for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Add in ancilla-only measurements:
    for i in range(krylov_dim - 1):
        # Get expectation values from experiment
        expval_real = results.data.evs[2 + 2 * obs_idx][0][
            i
        ]  # automatic extrapolated evs if ZNE is used
        expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][
            i
        ]  # automatic extrapolated evs if ZNE is used

        # Get expectation values
        expval = expval_real + 1j * expval_imag
        H_first_row[i + 1] += prefactors[i] * coeff * expval

H_first_row_list = H_first_row.tolist()

H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in itertools.product(range(krylov_dim), repeat=2):
    if i >= j:
        H_eff_circ[j, i] = H_first_row[i - j]
    else:
        H_eff_circ[j, i] = np.conj(H_first_row[j - i])

from sympy import Matrix

Matrix(H_eff_circ)

Output:

\displaystyle \left[\begin{matrix}10.0 & -3.02044405310714 - 2.80721615865252 i & 0.496487054782717 + 0.188101957039621 i & 1.0770511571923 + 0.104340737159455 i\\-3.02044405310714 + 2.80721615865252 i & 10.0 & -3.02044405310714 - 2.80721615865252 i & 0.496487054782717 + 0.188101957039621 i\\0.496487054782717 - 0.188101957039621 i & -3.02044405310714 + 2.80721615865252 i & 10.0 & -3.02044405310714 - 2.80721615865252 i\\1.0770511571923 - 0.104340737159455 i & 0.496487054782717 - 0.188101957039621 i & -3.02044405310714 + 2.80721615865252 i & 10.0\end{matrix}\right]

Por último, podemos resolver el problema de valores propios generalizado para $\tilde{H}$ :

$\tilde{H} \vec{c} = c S \vec{c}$

y obtener una estimación de la energía del estado fundamental $c_{min}$

gnd_en_circ_est_list = []
for d in range(1, krylov_dim + 1):
    # Solve generalized eigenvalue problem
    gnd_en_circ_est = solve_regularized_gen_eig(
        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=1e-1
    )
    gnd_en_circ_est_list.append(gnd_en_circ_est)
    print("The estimated ground state energy is: ", gnd_en_circ_est)

Output:

The estimated ground state energy is:  10.0
The estimated ground state energy is:  5.933953916292923
The estimated ground state energy is:  4.4101773995740645
The estimated ground state energy is:  3.921288588521255

Para un sector de una sola partícula, podemos calcular eficientemente el estado fundamental de este sector del Hamiltoniano clásicamente

gs_en = single_particle_gs(H_op, n_qubits)

Output:

n_sys_qubits 10
n_exc 1 , subspace dimension 11
single particle ground state energy:  2.391547869638771

len(H_op)

Output:

plt.plot(
    range(1, krylov_dim + 1),
    gnd_en_circ_est_list,
    color="blue",
    linestyle="-.",
    label="KQD estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.title("Estimating Ground state energy with Krylov Quantum Diagonalization")
plt.show()

Output:

5. Debate y ampliación

Para recapitular, empezamos con un estado de referencia y luego lo evolucionamos durante distintos periodos de tiempo para generar el subespacio unitario de Krylov. Proyectamos nuestro Hamiltoniano sobre ese subespacio. También estimamos los solapamientos de los vectores del subespacio. Por último, resolvemos clásicamente el problema de valores propios generalizado de dimensión inferior.

Un diagrama de flujo general de QKD: empezar con un estado de referencia, evolucionar el estado para aproximar los vectores de Krylov, proyectar en el subespacio de Krylov, diagonalizar el subespacio proyectado clásicamente y determinar las propiedades del estado base.

Comparemos lo que determina los costes computacionales de utilizar la técnica de Krylov de forma clásica y mecánica cuántica. No existen analogías perfectas entre los enfoques clásico y cuántico para todos los pasos. Esta tabla recoge algunas escalas de diferentes pasos para su consideración.

Una tabla que describe el escalado de diferentes procesos de forma clásica y en el enfoque cuántico de los métodos de Krylov. Algunos pasos cuánticos no tienen análogo. Las escalas son las mismas que las indicadas en el texto.

Recordemos que los hamiltonianos suelen tener términos que no pueden medirse simultáneamente (porque no conmutan entre sí). Clasificamos los términos del hamiltoniano en grupos de operadores de Pauli conmutables que pueden medirse simultáneamente, y es posible que necesitemos muchos de estos grupos para tener en cuenta todos los términos que no conmutan entre sí. Para construir $\tilde{H}$ en un ordenador cuántico se requieren mediciones separadas para cada grupo de cadenas de Pauli conmutativas en el Hamiltoniano, y cada una de ellas requiere muchos disparos. Debemos hacerlo para $r^2$ elementos de matriz diferentes, correspondientes a $r^2$ combinaciones de diferentes factores de evolución temporal. A veces hay maneras de reducir esto, pero en este tratamiento aproximado, el tiempo requerido para esto escala como $N_\text{shots}\times N_\text{GCP} \times r^2.$ Los elementos de $S$ deben ser estimados, que escala como $O(N_\text{shots}\times r^2)$ . Por último, la resolución del problema generalizado de valores propios en el espacio proyectado, clásicamente, toma $O(r^3).$

Vemos que la diagonalización cuántica de Krylov puede ser útil en casos en los que el número de grupos de Pauli conmutantes en el Hamiltoniano es relativamente pequeño. Estas dependencias de escala sugieren algunas aplicaciones en las que el método de Krylov puede ser útil, y otras en las que probablemente no lo sea. Algunos Hamiltonianos tienen una alta complejidad cuando se mapean a qubits, involucrando muchas cadenas de Pauli no conmutativas que no pueden ser fácilmente particionadas en unos pocos grupos conmutativos. Esto suele ocurrir, por ejemplo, con los problemas de química cuántica. Esta complejidad presenta dos retos principales para los ordenadores cuánticos a corto plazo:

La estimación de cada elemento de $\tilde{H}$ resulta costosa desde el punto de vista computacional debido al gran número de términos.
Los circuitos Trotter requeridos se vuelven prohibitivamente profundos.

Los dos puntos anteriores serán menos problemáticos cuando los ordenadores cuánticos alcancen la tolerancia a fallos, pero hay que tenerlos en cuenta a corto plazo. Incluso los sistemas con mapeos "más sencillos" que los de la química cuántica pueden experimentar los mismos impedimentos, si los hamiltonianos tienen demasiados términos no conmutativos. El método de Krylov es más útil cuando el Hamiltoniano puede dividirse en relativamente pocos grupos Pauli conmutativos, y cuando $H$ es fácil de implementar en circuitos trotadores. Ambas condiciones se cumplen, por ejemplo, para muchos modelos reticulares de interés en física. KQD es especialmente útil si se sabe muy poco sobre el estado básico. Esto se debe a sus garantías de convergencia inherentes y a su aplicabilidad en escenarios en los que los métodos alternativos son insostenibles debido a un conocimiento insuficiente del estado del suelo.

Aunque el KQD es una herramienta potente, los aspectos del protocolo que consumen mucho tiempo, en particular la estimación de cada elemento del Hamiltoniano proyectado y la superposición de los estados de Krylov, representan oportunidades de mejora. Un enfoque alternativo consiste en aprovechar los métodos de Krylov junto con los métodos basados en el muestreo, que son el tema de la siguiente lección.

6. Apéndices

Apéndice I: Subespacio de Krylov a partir de evoluciones en tiempo real

El espacio unitario de Krylov se define como

\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ |\psi\rangle, e^{-iH\,dt} |\psi\rangle, \dots, e^{-irH\,dt} |\psi\rangle \right\}

para algún paso temporal $dt$ que determinaremos más adelante. Supongamos temporalmente que $r$ es par: definamos entonces $d=r/2$ . Obsérvese que cuando proyectamos el Hamiltoniano en el espacio de Krylov anterior, es indistinguible del espacio de Krylov

\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ e^{i\,d\,H\,dt}|\psi\rangle, e^{i(d-1)H\,dt} |\psi\rangle, \dots, e^{-i(d-1)H\,dt} |\psi\rangle, e^{-i\,d\,H\,dt} |\psi\rangle \right\},

es decir, donde todas las evoluciones temporales se desplazan hacia atrás $d$ timesteps. La razón por la que es indistinguible es porque los elementos de la matriz

\tilde{H}_{j,k} = \langle\psi|e^{i\,j\,H\,dt}He^{-i\,k\,H\,dt}|\psi\rangle=\langle\psi|He^{i(j-k)H\,dt}|\psi\rangle

son invariantes bajo desplazamientos globales del tiempo de evolución, ya que las evoluciones temporales conmutan con el Hamiltoniano. Para los impares $r$ , podemos utilizar el análisis para $r-1$ .

Queremos demostrar que en algún lugar de este espacio de Krylov está garantizada la existencia de un estado de baja energía. Lo hacemos mediante el siguiente resultado, que se deriva del Teorema 3.1 de [3] :

Afirmación 1: existe una función $f$ tal que para energías $E$ en el rango espectral del Hamiltoniano (es decir, entre la energía del estado fundamental y la energía máxima)...

$f(E_0)=1$
$|f(E)|\le2\left(1 + \delta\right)^{-d}$ para todos los valores de $E$ que se encuentran a $\ge\delta$ de $E_0$ , es decir, se suprime exponencialmente
$f(E)$ es una combinación lineal de $e^{ijE\,dt}$ para $j=-d,-d+1,...,d-1,d$

A continuación ofrecemos una demostración, pero puede omitirse sin problemas a menos que se desee comprender el argumento completo y riguroso. Por ahora nos centraremos en las implicaciones de la afirmación anterior. Por la propiedad 3 anterior, podemos ver que el espacio de Krylov desplazado anterior contiene el estado $f(H)|\psi\rangle$ . Este es nuestro estado de baja energía. Para ver por qué, escriba $|\psi\rangle$ en la base propia de energía:

|\psi\rangle = \sum_{k=0}^{N}\gamma_k|E_k\rangle,

donde $|E_k\rangle$ es el k-ésimo eigenestado energético y $\gamma_k$ es su amplitud en el estado inicial $|\psi\rangle$ . Expresado en estos términos, $f(H)|\psi\rangle$ viene dado por

f(H)|\psi\rangle = \sum_{k=0}^{N}\gamma_kf(E_k)|E_k\rangle,

utilizando el hecho de que podemos sustituir $H$ por $E_k$ cuando actúa sobre el estado propio $|E_k\rangle$ . Por tanto, el error energético de este estado es

\text{energy error} = \frac{\langle\psi|f(H)(H-E_0)f(H)|\psi\rangle}{\langle\psi|f(H)^2|\psi\rangle}

= \frac{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.

Para convertir esto en un límite superior más fácil de entender, primero separamos la suma en el numerador en términos con $E_k-E_0\le\delta$ y términos con $E_k-E_0>\delta$ :

\text{energy error} = \frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} + \frac{\sum_{E_k> E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.

Podemos acotar el primer término en $\delta$ ,

\frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} < \frac{\delta\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} \le \delta,

donde el primer paso se sigue porque $E_k-E_0\le\delta$ para cada $E_k$ en la suma, y el segundo paso se sigue porque la suma en el numerador es un subconjunto de la suma en el denominador. Para el segundo término, primero acotamos a la baja el denominador en $|\gamma_0|^2$ , ya que $f(E_0)^2=1$ : sumando todo, se obtiene

\text{energy error} \le \delta + \frac{1}{|\gamma_0|^2}\sum_{E_k>E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0).

Para simplificar lo que queda, observe que para todos estos $E_k$ , por la definición de $f$ sabemos que $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$ . Además, si acotamos $E_k-E_0<2\|H\|$ y acotamos $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ , obtenemos

\text{energy error} \le \delta + \frac{8}{|\gamma_0|^2}\|H\|\left(1 + \delta\right)^{-2d}.

Esto es válido para cualquier $\delta>0$ , así que si fijamos $\delta$ igual a nuestro error objetivo, entonces el límite de error anterior converge hacia él exponencialmente con la dimensión de Krylov $2d=r$ . Obsérvese también que si $\delta<E_1-E_0$ , el término $\delta$ desaparece por completo en el límite anterior.

Para completar el argumento, primero observamos que lo anterior es sólo el error de energía del estado particular $f(H)|\psi\rangle$ , en lugar del error de energía del estado de menor energía en el espacio de Krylov. Sin embargo, por el principio variacional (Rayleigh-Ritz), el error de energía del estado de menor energía en el espacio de Krylov está acotado superiormente por el error de energía de cualquier estado en el espacio de Krylov, por lo que lo anterior es también un límite superior en el error de energía del estado de menor energía, es decir, la salida del algoritmo de diagonalización cuántica de Krylov.

Se puede realizar un análisis similar al anterior que, además, tenga en cuenta el ruido y el procedimiento de umbralización comentado en el cuaderno. Véase [2] y [4] para este análisis.

Apéndice II: prueba de la reclamación 1

Lo siguiente se deriva en su mayor parte de [3], Teorema 3.1: sea $0 < a < b$ y sea $\Pi^*_d$ el espacio de polinomios residuales (polinomios cuyo valor en 0 es 1) de grado como máximo $d$ . La solución de

\beta(a, b, d) = \min_{p \in \Pi^*_d} \max_{x \in [a, b]} |p(x)| \quad

p^*(x) = \frac{T_d\left(\frac{b + a - 2x}{b - a}\right)}{T_d\left(\frac{b + a}{b - a}\right)}, \quad

y el valor mínimo correspondiente es

\beta(a, b, d) = T_d^{-1}\left(\frac{b + a}{b - a}\right).

Queremos convertir esto en una función que pueda expresarse naturalmente en términos de exponenciales complejos, porque esas son las evoluciones temporales reales que generan el espacio cuántico de Krylov. Para ello, es conveniente introducir la siguiente transformación de energías dentro del rango espectral del Hamiltoniano a números en el rango $[0,1]$ : definir

g(E) = \frac{1-\cos\big((E-E_0)dt\big)}{2},

donde $dt$ es un paso de tiempo tal que $-\pi < E_0dt < E_\text{max}dt < \pi$ . Obsérvese que $g(E_0)=0$ y $g(E)$ crecen a medida que $E$ se aleja de $E_0$ .

Ahora, utilizando el polinomio $p^*(x)$ con los parámetros a, b, d fijados en $a = g(E_0 + \delta)$ , $b = 1$ , y d = int( r/2 ), definimos la función:

f(E) = p^* \left( g(E) \right) = \frac{T_d\left(1 + 2\frac{\cos\big((E-E_0)dt\big) - \cos\big(\delta\,dt\big)}{1 +\cos\big(\delta\,dt\big)}\right)}{T_d\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)}

donde $E_0$ es la energía del estado básico. Podemos ver insertando $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ que $f(E)$ es un polinomio trigonométrico de grado $d$ , es decir, una combinación lineal de $e^{ijE\,dt}$ para $j=-d,-d+1,...,d-1,d$ . Además, a partir de la definición de $p^*(x)$ anterior tenemos que $f(E_0)=p(0)=1$ y para cualquier $E$ en el rango espectral tal que $\vert E-E_0 \vert > \delta$ tenemos

|f(E)| \le \beta(a, b, d) = T_d^{-1}\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)

\leq 2\left(1 + \delta\right)^{-d} = 2\left(1 + \delta\right)^{-\lfloor k/2\rfloor}.

Referencias:

[1] https://arxiv.org/abs/2407.14431

[2] https://arxiv.org/abs/1811.09025

[3] https://people.math.ethz.ch/~mhg/pub/biksm.pdf

[4] https://academic.oup.com/book/36426

[5] https://en.wikipedia.org/wiki/Krylov_subespacio

[6] Métodos de subespacios de Krylov: Principios y análisis, Jörg Liesen, Zdenek Strakos https://academic.oup.com/book/36426

[7] Métodos iterativos para sistemas lineales dispersos" por Yousef Saad

[8] "MINRES-QLP: A Krylov Subspace Method for Indefinite or Singular Symmetric Systems" por Sou-Cheng Choi, Christopher Paige, y Michael Saunders ( https://epubs.siam.org/doi/10.1137/100787921 )

[9] Ethan N. Epperly, Lin Lin y Yuji Nakatsukasa. "Una teoría de diagonalización de subespacios cuánticos". SIAM Journal on Matrix Analysis and Applications 43, 1263-1290 (2022).

[10] https://link.aps.org/doi/10.1103/PRXQuantum.4.030319

¿Le ha resultado útil esta página?

Informe de un error, de una errata o solicite contenido en GitHub.

Diagonalización cuántica de Krylov

1. Introducción a los métodos de Krylov

1.1 Un ejemplo sencillo: Krylov a mano

r=1r=1r=1 de las dimensiones:

r=2r=2r=2 de las dimensiones:

r=3r=3r=3 de las dimensiones:

Aclaración de términos

Comprueba tu comprensión

1.2 Tipos de métodos de Krylov

1.3 ¿Por qué funciona el método del subespacio de Krylov?

Comprueba tu comprensión

2. Métodos de Krylov en la computación clásica

2.1 Ejemplo sencillo a pequeña escala

2.2 Escalado temporal con dimensión matricial

Comprueba tu comprensión

3. Krylov mediante evolución temporal

3.1 Evolución temporal

3.2 Cómo implementar en un ordenador cuántico

Problema generalizado de valores propios

4. Diagonalización cuántica de Krylov en un ordenador cuántico

4.1 Paso 1: Asignar el problema a circuitos y operadores cuánticos

Preparación del estado

Evolución temporal

Prueba de Hadamard

4. Paso 2. Optimizar circuitos y operadores para el hardware de destino

Prueba de Hadamard eficiente

Descomponer el operador de evolución temporal con la descomposición de Trotter

Circuitos modelo para calcular elementos matriciales de S~\tilde{S}S~ y H~\tilde{H}H~ mediante la prueba de Hadamard

4.3 Paso 3. Ejecutar utilizando una primitiva « IBM Quantum »

4.4 Paso 4. Procesar y analizar los resultados

Calcular el hamiltoniano efectivo y las matrices de superposición

5. Debate y ampliación

6. Apéndices

Apéndice I: Subespacio de Krylov a partir de evoluciones en tiempo real

Apéndice II: prueba de la reclamación 1

Referencias:

$r=1$ de las dimensiones:

$r=2$ de las dimensiones:

$r=3$ de las dimensiones:

Circuitos modelo para calcular elementos matriciales de $\tilde{S}$ y $\tilde{H}$ mediante la prueba de Hadamard