{ "cells": [ { "cell_type": "markdown", "id": "509f7bd9-b597-4d23-b3af-0a76a7b4d33d", "metadata": {}, "source": [ "---\n", "title: \"Diagonalización cuántica de Krylov de hamiltonianos de red\"\n", "description: \"Implementar el algoritmo de diagonalización cuántica de Krylov (KQD) en el contexto de los patrones de Qiskit.\"\n", "---\n", "\n", "{/* cspell:ignore prefactors */}\n", "\n", "\n", "\n", "# Diagonalización cuántica de Krylov de hamiltonianos de red\n", "\n", "*Estimación de uso: 20 minutos en un Heron r2 (NOTA: Esto es sólo una estimación. Su tiempo de ejecución puede variar)*\n", "\n" ] }, { "cell_type": "markdown", "id": "921c7b04-5b5d-4cfc-aba0-15a61334e619", "metadata": {}, "source": [ "\n", "\n", "## En segundo plano\n", "\n", "Este tutorial muestra cómo implementar el Algoritmo de Diagonalización Cuántica de Krylov (KQD) en el contexto de los patrones Qiskit. Primero aprenderás la teoría que hay detrás del algoritmo y luego verás una demostración de su ejecución en una QPU.\n", "\n", "En todas las disciplinas nos interesa conocer las propiedades del estado básico de los sistemas cuánticos. Algunos ejemplos son la comprensión de la naturaleza fundamental de las partículas y las fuerzas, la predicción y comprensión del comportamiento de materiales complejos y la comprensión de las interacciones y reacciones bioquímicas. Debido al crecimiento exponencial del espacio de Hilbert y a la correlación que surge en los sistemas enredados, los algoritmos clásicos tienen dificultades para resolver este problema en sistemas cuánticos de tamaño creciente. En un extremo del espectro está el enfoque existente que aprovecha el hardware cuántico centrado en métodos cuánticos variacionales (por ejemplo, [el eigensolver cuántico variacional](/docs/tutorials/spin-chain-vqe) ). Estas técnicas se enfrentan a retos con los dispositivos actuales debido al elevado número de llamadas a funciones necesarias en el proceso de optimización, lo que añade una gran sobrecarga de recursos una vez introducidas las técnicas avanzadas de mitigación de errores, limitando así su eficacia a sistemas pequeños. En el otro extremo del espectro, están los métodos cuánticos tolerantes a fallos con garantías de rendimiento (por ejemplo, la [estimación cuántica de fase](https://arxiv.org/abs/quant-ph/0604193) ), que requieren circuitos profundos que sólo pueden ejecutarse en un dispositivo tolerante a fallos. Por estas razones, introducimos aquí un algoritmo cuántico basado en métodos de subespacios (como se describe en este [artículo de revisión](https://arxiv.org/abs/2312.00178) ), el algoritmo de diagonalización cuántica de Krylov (KQD). Este algoritmo funciona bien a gran escala \\[[1](#references) ] en el hardware cuántico existente, comparte [garantías de rendimiento](https://arxiv.org/abs/2110.07492) similares a las de la estimación de fase, es compatible con técnicas avanzadas de mitigación de errores y podría proporcionar resultados clásicamente inaccesibles.\n", "\n" ] }, { "cell_type": "markdown", "id": "5f698c82-95ca-4dc4-b9fa-d6e741e2c02c", "metadata": {}, "source": [ "\n", "\n", "## Requisitos\n", "\n", "Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:\n", "\n", "* Qiskit SDK v2.0 o posterior, con soporte [de visualización](/docs/api/qiskit/visualization)\n", "* Qiskit Runtime v0.22 o posterior ( `pip install qiskit-ibm-runtime` )\n", "\n" ] }, { "cell_type": "markdown", "id": "c44956e4-47ab-4b0f-9d6d-553080110062", "metadata": {}, "source": [ "\n", "\n", "## Configuración\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ce7d5adc-3ef2-4654-b865-14d5141ce41a", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy as sp\n", "import matplotlib.pylab as plt\n", "from typing import Union, List\n", "import itertools as it\n", "import copy\n", "from sympy import Matrix\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "from qiskit.quantum_info import SparsePauliOp, Pauli, StabilizerState\n", "from qiskit.circuit import Parameter, IfElseOp\n", "from qiskit import QuantumCircuit, QuantumRegister\n", "from qiskit.circuit.library import PauliEvolutionGate\n", "from qiskit.synthesis import LieTrotter\n", "from qiskit.transpiler import Target, CouplingMap\n", "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", "\n", "\n", "from qiskit_ibm_runtime import (\n", " QiskitRuntimeService,\n", " EstimatorV2 as Estimator,\n", ")\n", "\n", "\n", "def solve_regularized_gen_eig(\n", " h: np.ndarray,\n", " s: np.ndarray,\n", " threshold: float,\n", " k: int = 1,\n", " return_dimn: bool = False,\n", ") -> Union[float, List[float]]:\n", " \"\"\"\n", " Method for solving the generalized eigenvalue problem with regularization\n", "\n", " Args:\n", " h (numpy.ndarray):\n", " The effective representation of the matrix in the Krylov subspace\n", " s (numpy.ndarray):\n", " The matrix of overlaps between vectors of the Krylov subspace\n", " threshold (float):\n", " Cut-off value for the eigenvalue of s\n", " k (int):\n", " Number of eigenvalues to return\n", " return_dimn (bool):\n", " Whether to return the size of the regularized subspace\n", "\n", " Returns:\n", " lowest k-eigenvalue(s) that are the solution of the\n", " regularized generalized eigenvalue problem\n", "\n", "\n", " \"\"\"\n", " s_vals, s_vecs = sp.linalg.eigh(s)\n", " s_vecs = s_vecs.T\n", " good_vecs = np.array(\n", " [vec for val, vec in zip(s_vals, s_vecs) if val > threshold]\n", " )\n", " h_reg = good_vecs.conj() @ h @ good_vecs.T\n", " s_reg = good_vecs.conj() @ s @ good_vecs.T\n", " if k == 1:\n", " if return_dimn:\n", " return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)\n", " else:\n", " return sp.linalg.eigh(h_reg, s_reg)[0][0]\n", " else:\n", " if return_dimn:\n", " return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)\n", " else:\n", " return sp.linalg.eigh(h_reg, s_reg)[0][:k]\n", "\n", "\n", "def single_particle_gs(H_op, n_qubits):\n", " \"\"\"\n", " Find the ground state of the single particle(excitation) sector\n", " \"\"\"\n", " H_x = []\n", " for p, coeff in H_op.to_list():\n", " H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))\n", "\n", " H_z = []\n", " for p, coeff in H_op.to_list():\n", " H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))\n", "\n", " H_c = H_op.coeffs\n", "\n", " print(\"n_sys_qubits\", n_qubits)\n", "\n", " n_exc = 1\n", " sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))\n", " print(\"n_exc\", n_exc, \", subspace dimension\", sub_dimn)\n", "\n", " few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)\n", "\n", " # list all of the possible sets of n_exc indices of 1s in\n", " # n_exc-particle states\n", " sparse_vecs = [\n", " set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)\n", " ]\n", "\n", " m = 0\n", " for i, i_set in enumerate(sparse_vecs):\n", " for j, j_set in enumerate(sparse_vecs):\n", " m += 1\n", "\n", " if len(i_set.symmetric_difference(j_set)) <= 2:\n", " for p_x, p_z, coeff in zip(H_x, H_z, H_c):\n", " if i_set.symmetric_difference(j_set) == p_x:\n", " sgn = ((-1j) ** len(p_x.intersection(p_z))) * (\n", " (-1) ** len(i_set.intersection(p_z))\n", " )\n", " else:\n", " sgn = 0\n", "\n", " few_particle_H[i, j] += sgn * coeff\n", "\n", " gs_en = min(np.linalg.eigvalsh(few_particle_H))\n", " print(\"single particle ground state energy: \", gs_en)\n", " return gs_en" ] }, { "cell_type": "markdown", "id": "76a1c616-a749-48a2-b5fd-8beff760760b", "metadata": {}, "source": [ "\n", "\n", "## Paso 1: Asignar entradas clásicas a un problema cuántico\n", "\n" ] }, { "cell_type": "markdown", "id": "a166e2a1-3003-4799-9c15-59f5e78b6ea8", "metadata": {}, "source": [ "\n", "\n", "### El espacio de Krylov\n", "\n", "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$, por un vector de referencia $\\vert v \\rangle$.\n", "\n", "$$\n", "\\mathcal{K}^r = \\left\\{ \\vert v \\rangle, A \\vert v \\rangle, A^2 \\vert v \\rangle, ..., A^{r-1} \\vert v \\rangle \\right\\}\n", "$$\n", "\n", "Si la matriz $A$ es el Hamiltoniano $H$, nos referiremos al espacio correspondiente como el espacio de Krylov de potencia $\\mathcal{K}_P$. En el caso de que $A$ sea el operador de evolución temporal generado por el hamiltoniano $U=e^{-iHt}$, nos referiremos al espacio como el espacio unitario de Krylov $\\mathcal{K}_U$. El subespacio de Krylov de potencia que utilizamos clásicamente no puede generarse directamente en un ordenador cuántico, ya que $H$ no es un operador unitario. En su lugar, podemos utilizar el operador de evolución en el tiempo $U = e^{-iHt}$ que puede demostrarse que ofrece [garantías de convergencia](https://arxiv.org/abs/2110.07492) similares a las del método de la potencia. Las potencias de $U$ se convierten entonces en diferentes pasos temporales $U^k = e^{-iH(kt)}$.\n", "\n", "$$\n", "\\mathcal{K}_U^r = \\left\\{ \\vert \\psi \\rangle, U \\vert \\psi \\rangle, U^2 \\vert \\psi \\rangle, ..., U^{r-1} \\vert \\psi \\rangle \\right\\}\n", "$$\n", "\n", "Véase el Apéndice para una derivación detallada de cómo el espacio unitario de Krylov permite representar con precisión los estados propios de baja energía.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "5573ca7e-ab16-4488-88b3-d8d1eba9e20c", "metadata": {}, "source": [ "\n", "\n", "### Algoritmo de diagonalización cuántica de Krylov\n", "\n", "Dado un Hamiltoniano $H$ que deseamos diagonalizar, primero consideramos el correspondiente espacio unitario de Krylov $\\mathcal{K}_U$. El objetivo es encontrar una representación compacta del Hamiltoniano en $\\mathcal{K}_U$, al que nos referiremos como $\\tilde{H}$. Los elementos matriciales de $\\tilde{H}$, la proyección del hamiltoniano en el espacio de Krylov, pueden calcularse calculando los siguientes valores de expectativa\n", "\n", "$$\n", "\\tilde{H}_{mn} = \\langle \\psi_m \\vert H \\vert \\psi_n \\rangle =\n", "$$\n", "\n", "$$\n", "= \\langle \\psi \\vert e^{i H t_m} H e^{-i H t_n} \\vert \\psi \\rangle\n", "$$\n", "\n", "$$\n", "= \\langle \\psi \\vert e^{i H m dt} H e^{-i H n dt} \\vert \\psi \\rangle\n", "$$\n", "\n", "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. Este tutorial se centra 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 las eigenenergías) podemos entonces diagonalizar fácilmente el Hamiltoniano proyectado $\\tilde{H}$. 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}$\n", "\n", "$$\n", "\\tilde{S}_{mn} = \\langle \\psi_m \\vert \\psi_n \\rangle\n", "$$\n", "\n", "Esto nos permite resolver el problema de valores propios en un espacio no ortogonal (también llamado problema de valores propios generalizado)\n", "\n", "$$\n", "\\tilde{H} \\ \\vec{c} = E \\ \\tilde{S} \\ \\vec{c}\n", "$$\n", "\n", "A continuación, se pueden obtener estimaciones de los valores propios y los estados propios de $H$ observando los de $\\tilde{H}$. Por ejemplo, la estimación de la energía del estado fundamental se obtiene tomando el valor propio más pequeño $c$ y el estado fundamental del vector propio correspondiente $\\vec{c}$. Los coeficientes de $\\vec{c}$ determinan la contribución de los distintos vectores que abarcan $\\mathcal{K}_U$.\n", "\n", "![fig1.png](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-subspace-diagonalization/fc662b76-8ad7-4a6c-8c49-5f08c125aee8.avif)\n", "\n", "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$. A continuación se ofrece una derivación más detallada de las operaciones calculadas por la prueba de Hadamard.\n", "\n" ] }, { "cell_type": "markdown", "id": "1a6d7a4f-6c1f-4069-93d1-f5b670645d7f", "metadata": {}, "source": [ "\n", "\n", "#### Definir hamiltoniano\n", "\n", "Consideremos el Hamiltoniano de Heisenberg para $N$ qubits en una cadena lineal: $H= \\sum_{i,j}^N X_i X_j + Y_i Y_j - J Z_i Z_j$\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "82163249-bafd-4bc6-9741-20a4077971a7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[('ZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 1), ('XXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 1), ('YYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 1)]\n" ] } ], "source": [ "# Define problem Hamiltonian.\n", "n_qubits = 30\n", "J = 1 # coupling strength for ZZ interaction\n", "\n", "# Define the Hamiltonian:\n", "H_int = [[\"I\"] * n_qubits for _ in range(3 * (n_qubits - 1))]\n", "for i in range(n_qubits - 1):\n", " H_int[i][i] = \"Z\"\n", " H_int[i][i + 1] = \"Z\"\n", "for i in range(n_qubits - 1):\n", " H_int[n_qubits - 1 + i][i] = \"X\"\n", " H_int[n_qubits - 1 + i][i + 1] = \"X\"\n", "for i in range(n_qubits - 1):\n", " H_int[2 * (n_qubits - 1) + i][i] = \"Y\"\n", " H_int[2 * (n_qubits - 1) + i][i + 1] = \"Y\"\n", "H_int = [\"\".join(term) for term in H_int]\n", "H_tot = [(term, J) if term.count(\"Z\") == 2 else (term, 1) for term in H_int]\n", "\n", "# Get operator\n", "H_op = SparsePauliOp.from_list(H_tot)\n", "print(H_tot)" ] }, { "cell_type": "markdown", "id": "f26578ba-24ad-42f3-baad-c348d3c05699", "metadata": {}, "source": [ "\n", "\n", "#### Establecer los parámetros del algoritmo\n", "\n" ] }, { "cell_type": "markdown", "id": "8d376c0e-fb7a-4a41-838a-ae041d5c9afa", "metadata": {}, "source": [ "Elegimos heurísticamente un valor para el paso temporal `dt` (basado en límites superiores de la norma hamiltoniana). Ref [\\[2\\]](#references) 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.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "963dd2e9-f4ea-456d-b4ed-e3d05a114082", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(0.10833078115826875)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get Hamiltonian restricted to single-particle states\n", "single_particle_H = np.zeros((n_qubits, n_qubits))\n", "for i in range(n_qubits):\n", " for j in range(i + 1):\n", " for p, coeff in H_op.to_list():\n", " p_x = Pauli(p).x\n", " p_z = Pauli(p).z\n", " if all(\n", " p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)\n", " ):\n", " sgn = (\n", " (-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))\n", " ) * ((-1) ** p_z[i])\n", " else:\n", " sgn = 0\n", " single_particle_H[i, j] += sgn * coeff\n", "for i in range(n_qubits):\n", " for j in range(i + 1, n_qubits):\n", " single_particle_H[i, j] = np.conj(single_particle_H[j, i])\n", "\n", "# Set dt according to spectral norm\n", "dt = np.pi / np.linalg.norm(single_particle_H, ord=2)\n", "dt" ] }, { "cell_type": "markdown", "id": "d2bcc2b9-ca7c-4147-8e6b-0f15208297ce", "metadata": {}, "source": [ "Y establecer otros parámetros del algoritmo. Por el bien de este tutorial, nos limitaremos a utilizar un espacio de Krylov con sólo cinco dimensiones, lo cual es bastante limitante.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "ddde76b8-446f-4cd5-bc4b-8f00e2ec726c", "metadata": {}, "outputs": [], "source": [ "# Set parameters for quantum Krylov algorithm\n", "krylov_dim = 5 # size of Krylov subspace\n", "num_trotter_steps = 6\n", "dt_circ = dt / num_trotter_steps" ] }, { "cell_type": "markdown", "id": "3c0dcc08-6963-4e14-a50c-21c5dd7f8ade", "metadata": {}, "source": [ "\n", "\n", "#### Preparación del estado\n", "\n", "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.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "70161afe-8ace-4642-894a-cd21ed77a3b9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "qc_state_prep = QuantumCircuit(n_qubits)\n", "qc_state_prep.x(int(n_qubits / 2) + 1)\n", "qc_state_prep.draw(\"mpl\", scale=0.5)" ] }, { "cell_type": "markdown", "id": "1786de6b-476f-47ea-9c20-3464d3cfbbdb", "metadata": {}, "source": [ "\n", "\n", "#### Evolución temporal\n", "\n", "Podemos realizar el operador de evolución temporal generado por un Hamiltoniano dado: $U=e^{-iHt}$ mediante la [aproximación de Lie-Trotter](/docs/api/qiskit/qiskit.synthesis.LieTrotter).\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "a23b9e9c-5dc8-447f-8c73-b4fa01630c8f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = Parameter(\"t\")\n", "\n", "## Create the time-evo op circuit\n", "evol_gate = PauliEvolutionGate(\n", " H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)\n", ")\n", "\n", "qr = QuantumRegister(n_qubits)\n", "qc_evol = QuantumCircuit(qr)\n", "qc_evol.append(evol_gate, qargs=qr)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "41c2cc43-51a2-42c3-88ee-a3b7b30eabb4", "metadata": {}, "source": [ "\n", "\n", "#### Prueba de Hadamard\n", "\n", "![fig2.png](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-subspace-diagonalization/c5263851-6067-4ca2-8e0c-a835631cdc7f.avif)\n", "\n", "$$\n", "\\begin{equation*}\n", " |0\\rangle|0\\rangle^N \\quad\\longrightarrow\\quad \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle + |1\\rangle \\Big)|0\\rangle^N \\quad\\longrightarrow\\quad \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle|0\\rangle^N+|1\\rangle |\\psi_i\\rangle\\Big) \\quad\\longrightarrow\\quad \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle |0\\rangle^N+|1\\rangle P |\\psi_i\\rangle\\Big) \\quad\\longrightarrow\\quad\\frac{1}{\\sqrt{2}}\\Big(|0\\rangle |\\psi_j\\rangle+|1\\rangle P|\\psi_i\\rangle\\Big)\n", "\\end{equation*}\n", "$$\n", "\n", "Donde $P$ es uno de los términos de la descomposición del Hamiltoniano $H=\\sum P$ y $\\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$. Para medir $X$, primero se aplica $H$...\n", "\n", "$$\n", "\\begin{equation*}\n", " \\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)\n", "\\end{equation*}\n", "$$\n", "\n", "... y luego medir:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\begin{split}\n", " \\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) \\\\\n", " &= \\text{Re}\\Big[\\langle\\psi_j| P|\\psi_i\\rangle\\Big].\n", "\\end{split}\n", "\\end{equation*}\n", "$$\n", "\n", "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\n", "\n", "$$\n", "\\begin{equation*}\n", " \\langle Y\\rangle = \\text{Im}\\Big[\\langle\\psi_j| P|\\psi_i\\rangle\\Big].\n", "\\end{equation*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "7c1efca7-7db9-43a9-bcb7-053ba274d6f6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Circuit for calculating the real part of the overlap in S via Hadamard test\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Create the time-evo op circuit\n", "evol_gate = PauliEvolutionGate(\n", " H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)\n", ")\n", "\n", "## Create the time-evo op dagger circuit\n", "evol_gate_d = PauliEvolutionGate(\n", " H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)\n", ")\n", "evol_gate_d = evol_gate_d.inverse()\n", "\n", "# Put pieces together\n", "qc_reg = QuantumRegister(n_qubits)\n", "qc_temp = QuantumCircuit(qc_reg)\n", "qc_temp.compose(qc_state_prep, inplace=True)\n", "for _ in range(num_trotter_steps):\n", " qc_temp.append(evol_gate, qargs=qc_reg)\n", "for _ in range(num_trotter_steps):\n", " qc_temp.append(evol_gate_d, qargs=qc_reg)\n", "qc_temp.compose(qc_state_prep.inverse(), inplace=True)\n", "\n", "# Create controlled version of the circuit\n", "controlled_U = qc_temp.to_gate().control(1)\n", "\n", "# Create hadamard test circuit for real part\n", "qr = QuantumRegister(n_qubits + 1)\n", "qc_real = QuantumCircuit(qr)\n", "qc_real.h(0)\n", "qc_real.append(controlled_U, list(range(n_qubits + 1)))\n", "qc_real.h(0)\n", "\n", "print(\n", " \"Circuit for calculating the real part of the overlap in S via Hadamard test\"\n", ")\n", "qc_real.draw(\"mpl\", fold=-1, scale=0.5)" ] }, { "cell_type": "markdown", "id": "c72fd5e5-a586-498e-bbc1-0f48ecea49d5", "metadata": {}, "source": [ "El circuito de prueba Hadamard puede ser un circuito profundo una vez que descomponemos a puertas nativas (que aumentará aún más si tenemos en cuenta la topología del dispositivo)\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "06f22c05-6e1c-4540-9dac-884b3400a50e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of layers of 2Q operations 112753\n" ] } ], "source": [ "print(\n", " \"Number of layers of 2Q operations\",\n", " qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "id": "8e15903d-d4a8-4fbe-8f9e-e960e686629e", "metadata": {}, "source": [ "\n", "\n", "## Paso 2: Optimizar el problema para la ejecución en hardware cuántico\n", "\n", "\n", "\n", "### Prueba de Hadamard eficiente\n", "\n", "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:\n", "\n", "![fig3.png](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-subspace-diagonalization/35b13797-5a46-486c-b50e-97c205cc9747.avif)\n", "\n", "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).\n", "Dado que la puerta $\\text{Prep} \\; \\psi$, prepara el estado de referencia deseado $\\ket{psi} = \\text{Prep} \\; \\psi \\ket{0} = e^{-i H 0 dt} U_{\\psi} \\ket{0}$, por ejemplo, preparar el estado HF para la química $\\text{Prep} \\; \\psi$ sería un producto de NOTs de un solo qubit, por lo que controlada- $\\text{Prep} \\; \\psi$ es sólo un producto de CNOTs.\n", "Entonces el circuito anterior implementa el siguiente estado antes de la medición:\n", "\n", "$$\n", "\\begin{equation}\n", "\\begin{split}\n", " \\ket{0} \\ket{0}^N\\xrightarrow{H}&\\frac{1}{\\sqrt{2}}\n", " \\left(\n", " \\ket{0}\\ket{0}^N+ \\ket{1} \\ket{0}^N\n", " \\right)\\\\\n", " \\xrightarrow{\\text{1-ctrl-init}}&\\frac{1}{\\sqrt{2}}\\left(|0\\rangle|0\\rangle^N+|1\\rangle|\\psi\\rangle\\right)\\\\\n", " \\xrightarrow{U}&\\frac{1}{\\sqrt{2}}\\left(e^{i\\phi}\\ket{0}\\ket{0}^N+\\ket{1} U\\ket{\\psi}\\right)\\\\\n", " \\xrightarrow{\\text{0-ctrl-init}}&\\frac{1}{\\sqrt{2}}\n", " \\left(\n", " e^{i\\phi}\\ket{0} \\ket{\\psi}\n", " +\\ket{1} U\\ket{\\psi}\n", " \\right)\\\\\n", " =&\\frac{1}{2}\n", " \\left(\n", " \\ket{+}\\left(e^{i\\phi}\\ket{\\psi}+U\\ket{\\psi}\\right)\n", " +\\ket{-}\\left(e^{i\\phi}\\ket{\\psi}-U\\ket{\\psi}\\right)\n", " \\right)\\\\\n", " =&\\frac{1}{2}\n", " \\left(\n", " \\ket{+i}\\left(e^{i\\phi}\\ket{\\psi}-iU\\ket{\\psi}\\right)\n", " +\\ket{-i}\\left(e^{i\\phi}\\ket{\\psi}+iU\\ket{\\psi}\\right)\n", " \\right)\n", "\\end{split}\n", "\\end{equation}\n", "$$\n", "\n", "donde hemos utilizado el desfase simulable clásico $ U\\ket{0}^N = e^{i\\phi}\\ket{0}^N$ en la tercera línea. Por lo tanto, los valores de las expectativas se obtienen como\n", "\n", "$$\n", "\\begin{equation}\n", "\\begin{split}\n", " \\langle X\\otimes P\\rangle&=\\frac{1}{4}\n", " \\Big(\n", " \\left(e^{-i\\phi}\\bra{\\psi}+\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}+U\\ket{\\psi}\\right)\n", " \\\\\n", " &\\qquad-\\left(e^{-i\\phi}\\bra{\\psi}-\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}-U\\ket{\\psi}\\right)\n", " \\Big)\\\\\n", " &=\\text{Re}\\left[e^{-i\\phi}\\bra{\\psi}PU\\ket{\\psi}\\right],\n", "\\end{split}\n", "\\end{equation}\n", "$$\n", "\n", "$$\n", "\\begin{equation}\n", "\\begin{split}\n", " \\langle Y\\otimes P\\rangle&=\\frac{1}{4}\n", " \\Big(\n", " \\left(e^{-i\\phi}\\bra{\\psi}+i\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}-iU\\ket{\\psi}\\right)\n", " \\\\\n", " &\\qquad-\\left(e^{-i\\phi}\\bra{\\psi}-i\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}+iU\\ket{\\psi}\\right)\n", " \\Big)\\\\\n", " &=\\text{Im}\\left[e^{-i\\phi}\\bra{\\psi}PU\\ket{\\psi}\\right].\n", "\\end{split}\n", "\\end{equation}\n", "$$\n", "\n", "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$ 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.\n", "\n" ] }, { "cell_type": "markdown", "id": "2e29da5f-b5f5-4b9e-96fc-3fa7ab698398", "metadata": {}, "source": [ "\n", "\n", "### Descomponer el operador de evolución temporal con la descomposición de Trotter\n", "\n", "En lugar de implementar exactamente el operador de evolución temporal, podemos utilizar la descomposición de Trotter para implementar una aproximación del mismo. La repetición varias veces de una descomposición de Trotter de cierto orden nos proporciona una mayor reducción del error introducido por la aproximación. A continuación, construimos directamente la implementación de Trotter de la forma más eficiente para el grafo de interacciones del Hamiltoniano que estamos considerando (sólo interacciones de vecino más cercano). En la práctica insertamos rotaciones Pauli $R_{xx}$, $R_{yy}$, $R_{zz}$ con un ángulo parametrizado $t$ que corresponden a la implementación aproximada de $e^{-i (XX + YY + ZZ) t}$. Dada la diferencia de definición de las rotaciones Pauli y la evolución temporal que intentamos 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 el número impar de repeticiones de los pasos de Trotter, lo que es funcionalmente equivalente pero permite sintetizar operaciones adyacentes en un único $SU(2)$ unitario. Esto da un circuito mucho menos profundo que el que se obtiene utilizando la funcionalidad genérica `PauliEvolutionGate()` .\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "267716dc-fa23-41bd-abe4-6d4e0499a0f4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = Parameter(\"t\")\n", "\n", "# Create instruction for rotation about XX+YY-ZZ:\n", "Rxyz_circ = QuantumCircuit(2)\n", "Rxyz_circ.rxx(t, 0, 1)\n", "Rxyz_circ.ryy(t, 0, 1)\n", "Rxyz_circ.rzz(t, 0, 1)\n", "Rxyz_instr = Rxyz_circ.to_instruction(label=\"RXX+YY+ZZ\")\n", "\n", "interaction_list = [\n", " [[i, i + 1] for i in range(0, n_qubits - 1, 2)],\n", " [[i, i + 1] for i in range(1, n_qubits - 1, 2)],\n", "] # linear chain\n", "\n", "qr = QuantumRegister(n_qubits)\n", "trotter_step_circ = QuantumCircuit(qr)\n", "for i, color in enumerate(interaction_list):\n", " for interaction in color:\n", " trotter_step_circ.append(Rxyz_instr, interaction)\n", " if i < len(interaction_list) - 1:\n", " trotter_step_circ.barrier()\n", "reverse_trotter_step_circ = trotter_step_circ.reverse_ops()\n", "\n", "qc_evol = QuantumCircuit(qr)\n", "for step in range(num_trotter_steps):\n", " if step % 2 == 0:\n", " qc_evol = qc_evol.compose(trotter_step_circ)\n", " else:\n", " qc_evol = qc_evol.compose(reverse_trotter_step_circ)\n", "\n", "qc_evol.decompose().draw(\"mpl\", fold=-1, scale=0.5)" ] }, { "cell_type": "markdown", "id": "0a9c3a4d-a678-41e4-b9dd-995fe34341fb", "metadata": {}, "source": [ "\n", "\n", "### Utilice un circuito optimizado para la preparación del estado\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "70411715-eed3-4cf5-961d-06a6f1e04efc", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "control = 0\n", "excitation = int(n_qubits / 2) + 1\n", "controlled_state_prep = QuantumCircuit(n_qubits + 1)\n", "controlled_state_prep.cx(control, excitation)\n", "controlled_state_prep.draw(\"mpl\", fold=-1, scale=0.5)" ] }, { "cell_type": "markdown", "id": "2421ff11-d97a-488a-a382-a652df8c94d6", "metadata": {}, "source": [ "\n", "\n", "### Circuitos modelo para calcular elementos matriciales de $\\tilde{S}$ y $\\tilde{H}$ mediante la prueba de Hadamard\n", "\n", "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.\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "2f102112-4ddc-41ea-999c-db5863bc77ac", "metadata": {}, "outputs": [], "source": [ "# Parameters for the template circuits\n", "parameters = []\n", "for idx in range(1, krylov_dim):\n", " parameters.append(2 * dt_circ * (idx))" ] }, { "cell_type": "code", "execution_count": 13, "id": "33ec7c29-904e-4445-a654-405214349a4d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create modified hadamard test circuit\n", "qr = QuantumRegister(n_qubits + 1)\n", "qc = QuantumCircuit(qr)\n", "qc.h(0)\n", "qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)\n", "qc.barrier()\n", "qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)\n", "qc.barrier()\n", "qc.x(0)\n", "qc.compose(\n", " controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True\n", ")\n", "qc.x(0)\n", "\n", "qc.decompose().draw(\"mpl\", fold=-1)" ] }, { "cell_type": "code", "execution_count": 14, "id": "157356c9-06bb-411c-87ef-cd1d2f73be6f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The optimized circuit has 2Q gates depth: 74\n" ] } ], "source": [ "print(\n", " \"The optimized circuit has 2Q gates depth: \",\n", " qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),\n", ")" ] }, { "cell_type": "markdown", "id": "5a59bbc4-72c7-4a55-9a9e-57bb7b469021", "metadata": {}, "source": [ "Hemos reducido considerablemente la profundidad de la prueba de Hadamard con una combinación de aproximación de Trotter y unitarios no controlados\n", "\n" ] }, { "cell_type": "markdown", "id": "ca78bfe3-684f-4d3a-b8c2-3d4e55c9ec30", "metadata": {}, "source": [ "\n", "\n", "## Paso 3: Ejecutar utilizando Qiskit primitives\n", "\n" ] }, { "cell_type": "markdown", "id": "a082f023-d752-40e4-b693-c4d0da5ba102", "metadata": {}, "source": [ "Instanciar el backend y establecer los parámetros de ejecución\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0d90e4df-e262-4852-a811-ff6a1d3232ae", "metadata": {}, "outputs": [], "source": [ "service = QiskitRuntimeService()\n", "backend = service.least_busy(operational=True, simulator=False)\n", "if (\n", " \"if_else\" not in backend.target.operation_names\n", "): # Needed as \"op_name\" could be \"if_else\"\n", " backend.target.add_instruction(IfElseOp, name=\"if_else\")\n", "print(backend.name)" ] }, { "cell_type": "markdown", "id": "2b3f7227-e6ef-4a11-93d6-d3696054b9cc", "metadata": {}, "source": [ "\n", "\n", "### Transpiling a una QPU\n", "\n", "En primer lugar, seleccionemos subconjuntos del mapa de acoplamiento con qubits de «buen» rendimiento (donde «buen» es bastante arbitrario en este caso, ya que lo que queremos principalmente es evitar los qubits de muy mal rendimiento) y creemos un nuevo objetivo para la transpilación\n", "\n" ] }, { "cell_type": "code", "execution_count": 32, "id": "e123fda1-6454-4893-9612-2b8591e8cfb9", "metadata": {}, "outputs": [], "source": [ "target = backend.target\n", "cmap = target.build_coupling_map(filter_idle_qubits=True)\n", "cmap_list = list(cmap.get_edges())\n", "\n", "cust_cmap_list = copy.deepcopy(cmap_list)\n", "for q in range(target.num_qubits):\n", " meas_err = target[\"measure\"][(q,)].error\n", " t2 = target.qubit_properties[q].t2 * 1e6\n", " if meas_err > 0.02 or t2 < 100:\n", " for q_pair in cmap_list:\n", " if q in q_pair:\n", " try:\n", " cust_cmap_list.remove(q_pair)\n", " except:\n", " continue\n", "\n", "for q in cmap_list:\n", " op_name = list(target.operation_names_for_qargs(q))[0]\n", " twoq_gate_err = target[f\"{op_name}\"][q].error\n", " if twoq_gate_err > 0.005:\n", " for q_pair in cmap_list:\n", " if q == q_pair:\n", " try:\n", " cust_cmap_list.remove(q)\n", " except:\n", " continue\n", "\n", "\n", "cust_cmap = CouplingMap(cust_cmap_list)\n", "cust_target = Target.from_configuration(\n", " basis_gates=backend.configuration().basis_gates,\n", " coupling_map=cust_cmap,\n", ")" ] }, { "cell_type": "markdown", "id": "96f79e35-a265-4de4-b98c-f9c84f58f0d5", "metadata": {}, "source": [ "A continuación, transpile el circuito virtual a la mejor disposición física en este nuevo objetivo\n", "\n" ] }, { "cell_type": "code", "execution_count": 36, "id": "62109c95-79fe-4076-b591-9bd920fd51f4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "depth 52\n", "num 2q ops OrderedDict([('rz', 2058), ('sx', 1703), ('cz', 728), ('x', 84), ('barrier', 8)])\n", "physical qubits [91, 92, 93, 94, 95, 98, 99, 108, 109, 110, 111, 113, 114, 115, 119, 127, 132, 133, 134, 135, 137, 139, 147, 148, 149, 150, 151, 152, 153, 154, 155]\n" ] } ], "source": [ "basis_gates = list(target.operation_names)\n", "pm = generate_preset_pass_manager(\n", " optimization_level=3,\n", " target=cust_target,\n", " basis_gates=basis_gates,\n", ")\n", "\n", "qc_trans = pm.run(qc)\n", "\n", "print(\"depth\", qc_trans.depth(lambda x: x[0].num_qubits == 2))\n", "print(\"num 2q ops\", qc_trans.count_ops())\n", "print(\n", " \"physical qubits\",\n", " sorted(\n", " [\n", " idx\n", " for idx, qb in qc_trans.layout.initial_layout.get_physical_bits().items()\n", " if qb._register.name != \"ancilla\"\n", " ]\n", " ),\n", ")" ] }, { "cell_type": "markdown", "id": "9b9b0170-b96a-47b9-b5ab-5112d02338a3", "metadata": {}, "source": [ "\n", "\n", "### Crear PUB para ejecución con Estimator\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "bcfd6693-d1fb-44e4-9b06-c70e6766e877", "metadata": {}, "outputs": [], "source": [ "# Define observables to measure for S\n", "observable_S_real = \"I\" * (n_qubits) + \"X\"\n", "observable_S_imag = \"I\" * (n_qubits) + \"Y\"\n", "\n", "observable_op_real = SparsePauliOp(\n", " observable_S_real\n", ") # define a sparse pauli operator for the observable\n", "observable_op_imag = SparsePauliOp(observable_S_imag)\n", "\n", "layout = qc_trans.layout # get layout of transpiled circuit\n", "observable_op_real = observable_op_real.apply_layout(\n", " layout\n", ") # apply physical layout to the observable\n", "observable_op_imag = observable_op_imag.apply_layout(layout)\n", "observable_S_real = (\n", " observable_op_real.paulis.to_labels()\n", ") # get the label of the physical observable\n", "observable_S_imag = observable_op_imag.paulis.to_labels()\n", "\n", "observables_S = [[observable_S_real], [observable_S_imag]]\n", "\n", "\n", "# Define observables to measure for H\n", "# Hamiltonian terms to measure\n", "observable_list = []\n", "for pauli, coeff in zip(H_op.paulis, H_op.coeffs):\n", " # print(pauli)\n", " observable_H_real = pauli[::-1].to_label() + \"X\"\n", " observable_H_imag = pauli[::-1].to_label() + \"Y\"\n", " observable_list.append([observable_H_real])\n", " observable_list.append([observable_H_imag])\n", "\n", "layout = qc_trans.layout\n", "\n", "observable_trans_list = []\n", "for observable in observable_list:\n", " observable_op = SparsePauliOp(observable)\n", " observable_op = observable_op.apply_layout(layout)\n", " observable_trans_list.append([observable_op.paulis.to_labels()])\n", "\n", "observables_H = observable_trans_list\n", "\n", "\n", "# Define a sweep over parameter values\n", "params = np.vstack(parameters).T\n", "\n", "\n", "# Estimate the expectation value for all combinations of\n", "# observables and parameter values, where the pub result will have\n", "# shape (# observables, # parameter values).\n", "pub = (qc_trans, observables_S + observables_H, params)" ] }, { "cell_type": "markdown", "id": "8bd27e84-d673-4326-9004-f3be5a65aed3", "metadata": {}, "source": [ "\n", "\n", "### Circuitos de carrera\n", "\n", "Los circuitos de $t=0$ se pueden calcular de forma clásica\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "768cd443-2ed8-4ba2-ba9d-98491b36fa54", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(25+0j)\n" ] } ], "source": [ "qc_cliff = qc.assign_parameters({t: 0})\n", "\n", "\n", "# Get expectation values from experiment\n", "S_expval_real = StabilizerState(qc_cliff).expectation_value(\n", " Pauli(\"I\" * (n_qubits) + \"X\")\n", ")\n", "S_expval_imag = StabilizerState(qc_cliff).expectation_value(\n", " Pauli(\"I\" * (n_qubits) + \"Y\")\n", ")\n", "\n", "# Get expectation values\n", "S_expval = S_expval_real + 1j * S_expval_imag\n", "\n", "H_expval = 0\n", "for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):\n", " # Get expectation values from experiment\n", " expval_real = StabilizerState(qc_cliff).expectation_value(\n", " Pauli(pauli[::-1].to_label() + \"X\")\n", " )\n", " expval_imag = StabilizerState(qc_cliff).expectation_value(\n", " Pauli(pauli[::-1].to_label() + \"Y\")\n", " )\n", " expval = expval_real + 1j * expval_imag\n", "\n", " # Fill-in matrix elements\n", " H_expval += coeff * expval\n", "\n", "\n", "print(H_expval)" ] }, { "cell_type": "markdown", "id": "794edf4c-d539-4864-865d-52a88a450b5c", "metadata": {}, "source": [ "Ejecuta circuitos para $S$ y $\\tilde{H}$ con Estimator\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9059dc9e-9203-4572-92c4-76e9c85dfed9", "metadata": {}, "outputs": [], "source": [ "# Experiment options\n", "num_randomizations = 300\n", "num_randomizations_learning = 30\n", "shots_per_randomization = 100\n", "noise_factors = [1, 1.2, 1.4]\n", "learning_pair_depths = [0, 4, 24, 48]\n", "\n", "\n", "experimental_opts = {}\n", "experimental_opts[\"resilience\"] = {\n", " \"measure_mitigation\": True,\n", " \"measure_noise_learning\": {\n", " \"num_randomizations\": num_randomizations_learning,\n", " \"shots_per_randomization\": shots_per_randomization,\n", " },\n", " \"zne_mitigation\": True,\n", " \"zne\": {\"noise_factors\": noise_factors},\n", " \"layer_noise_learning\": {\n", " \"max_layers_to_learn\": 10,\n", " \"layer_pair_depths\": learning_pair_depths,\n", " \"shots_per_randomization\": shots_per_randomization,\n", " \"num_randomizations\": num_randomizations_learning,\n", " },\n", " \"zne\": {\n", " \"amplifier\": \"pea\",\n", " \"extrapolated_noise_factors\": [0] + noise_factors,\n", " },\n", "}\n", "experimental_opts[\"twirling\"] = {\n", " \"num_randomizations\": num_randomizations,\n", " \"shots_per_randomization\": shots_per_randomization,\n", " \"strategy\": \"all\",\n", "}\n", "\n", "estimator = Estimator(mode=backend, options=experimental_opts)\n", "\n", "\n", "job = estimator.run([pub])" ] }, { "cell_type": "markdown", "id": "9524e609-b457-42e9-9de7-8dbd4464ac74", "metadata": {}, "source": [ "\n", "\n", "## Paso 4: Procesamiento posterior y devolución del resultado en el formato clásico deseado\n", "\n" ] }, { "cell_type": "code", "execution_count": 46, "id": "30d72032-02e4-40e7-af28-cadda2f7f3bd", "metadata": {}, "outputs": [], "source": [ "results = job.result()[0]" ] }, { "cell_type": "markdown", "id": "8149d431-f069-4567-a0be-0cc4ed8d523b", "metadata": {}, "source": [ "\n", "\n", "### Calcular el hamiltoniano efectivo y las matrices de superposición\n", "\n", "En primer lugar, calcule la fase acumulada por el estado $\\vert 0 \\rangle$ durante la evolución temporal no controlada\n", "\n" ] }, { "cell_type": "code", "execution_count": 47, "id": "34ba8b05-d52c-4b0f-8d28-acaf7a5ddffc", "metadata": {}, "outputs": [], "source": [ "prefactors = [\n", " np.exp(-1j * sum([c for p, c in H_op.to_list() if \"Z\" in p]) * i * dt)\n", " for i in range(1, krylov_dim)\n", "]" ] }, { "cell_type": "markdown", "id": "36875c8b-36c8-464b-9f0e-7a79e17d5fef", "metadata": {}, "source": [ "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$\n", "\n" ] }, { "cell_type": "code", "execution_count": 48, "id": "16dd2534-8d5f-40f7-a938-5d519475fd8d", "metadata": {}, "outputs": [], "source": [ "# Assemble S, the overlap matrix of dimension D:\n", "S_first_row = np.zeros(krylov_dim, dtype=complex)\n", "S_first_row[0] = 1 + 0j\n", "\n", "# Add in ancilla-only measurements:\n", "for i in range(krylov_dim - 1):\n", " # Get expectation values from experiment\n", " expval_real = results.data.evs[0][0][\n", " i\n", " ] # automatic extrapolated evs if ZNE is used\n", " expval_imag = results.data.evs[1][0][\n", " i\n", " ] # automatic extrapolated evs if ZNE is used\n", "\n", " # Get expectation values\n", " expval = expval_real + 1j * expval_imag\n", " S_first_row[i + 1] += prefactors[i] * expval\n", "\n", "S_first_row_list = S_first_row.tolist() # for saving purposes\n", "\n", "\n", "S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)\n", "\n", "# Distribute entries from first row across matrix:\n", "for i, j in it.product(range(krylov_dim), repeat=2):\n", " if i >= j:\n", " S_circ[j, i] = S_first_row[i - j]\n", " else:\n", " S_circ[j, i] = np.conj(S_first_row[j - i])" ] }, { "cell_type": "code", "execution_count": 49, "id": "01c6563d-87ca-4bd6-b487-dcdaece2d8c2", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\n", "\\displaystyle \\left[\\begin{matrix}1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i & 0.0012070853532697 + 0.312052218182462 i\\\\-0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i\\\\0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i\\\\-0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i\\\\0.0012070853532697 - 0.312052218182462 i & -0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0\\end{matrix}\\right]\n", "$$" ], "text/plain": [ "Matrix([\n", "[ 1.0, -0.723052998582984 - 0.345085413575966*I, 0.467051960502366 + 0.516197865254034*I, -0.180546747798251 - 0.492624093654174*I, 0.0012070853532697 + 0.312052218182462*I],\n", "[-0.723052998582984 + 0.345085413575966*I, 1.0, -0.723052998582984 - 0.345085413575966*I, 0.467051960502366 + 0.516197865254034*I, -0.180546747798251 - 0.492624093654174*I],\n", "[ 0.467051960502366 - 0.516197865254034*I, -0.723052998582984 + 0.345085413575966*I, 1.0, -0.723052998582984 - 0.345085413575966*I, 0.467051960502366 + 0.516197865254034*I],\n", "[-0.180546747798251 + 0.492624093654174*I, 0.467051960502366 - 0.516197865254034*I, -0.723052998582984 + 0.345085413575966*I, 1.0, -0.723052998582984 - 0.345085413575966*I],\n", "[0.0012070853532697 - 0.312052218182462*I, -0.180546747798251 + 0.492624093654174*I, 0.467051960502366 - 0.516197865254034*I, -0.723052998582984 + 0.345085413575966*I, 1.0]])" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Matrix(S_circ)" ] }, { "cell_type": "markdown", "id": "a34c72c0-2984-47fb-8579-089ea795ef39", "metadata": {}, "source": [ "Y los elementos de la matriz de $\\tilde{H}$\n", "\n" ] }, { "cell_type": "code", "execution_count": 50, "id": "9cde2419-8a29-4d7a-b5b1-9ef2d4c126d2", "metadata": {}, "outputs": [], "source": [ "# Assemble S, the overlap matrix of dimension D:\n", "H_first_row = np.zeros(krylov_dim, dtype=complex)\n", "H_first_row[0] = H_expval\n", "\n", "for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):\n", " # Add in ancilla-only measurements:\n", " for i in range(krylov_dim - 1):\n", " # Get expectation values from experiment\n", " expval_real = results.data.evs[2 + 2 * obs_idx][0][\n", " i\n", " ] # automatic extrapolated evs if ZNE is used\n", " expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][\n", " i\n", " ] # automatic extrapolated evs if ZNE is used\n", "\n", " # Get expectation values\n", " expval = expval_real + 1j * expval_imag\n", " H_first_row[i + 1] += prefactors[i] * coeff * expval\n", "\n", "H_first_row_list = H_first_row.tolist()\n", "\n", "H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)\n", "\n", "# Distribute entries from first row across matrix:\n", "for i, j in it.product(range(krylov_dim), repeat=2):\n", " if i >= j:\n", " H_eff_circ[j, i] = H_first_row[i - j]\n", " else:\n", " H_eff_circ[j, i] = np.conj(H_first_row[j - i])" ] }, { "cell_type": "code", "execution_count": 51, "id": "5800d892-5554-4a92-bfe0-3750ef0a3ed7", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\n", "\\displaystyle \\left[\\begin{matrix}25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i & 1.98818301405581 + 5.8897614762563 i\\\\-14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i\\\\10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i\\\\-5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i\\\\1.98818301405581 - 5.8897614762563 i & -5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0\\end{matrix}\\right]\n", "$$" ], "text/plain": [ "Matrix([\n", "[ 25.0, -14.2437089383409 - 6.50486277982165*I, 10.2857217968584 + 9.0431912203186*I, -5.15587257589417 - 8.88280836036843*I, 1.98818301405581 + 5.8897614762563*I],\n", "[-14.2437089383409 + 6.50486277982165*I, 25.0, -14.2437089383409 - 6.50486277982165*I, 10.2857217968584 + 9.0431912203186*I, -5.15587257589417 - 8.88280836036843*I],\n", "[ 10.2857217968584 - 9.0431912203186*I, -14.2437089383409 + 6.50486277982165*I, 25.0, -14.2437089383409 - 6.50486277982165*I, 10.2857217968584 + 9.0431912203186*I],\n", "[-5.15587257589417 + 8.88280836036843*I, 10.2857217968584 - 9.0431912203186*I, -14.2437089383409 + 6.50486277982165*I, 25.0, -14.2437089383409 - 6.50486277982165*I],\n", "[ 1.98818301405581 - 5.8897614762563*I, -5.15587257589417 + 8.88280836036843*I, 10.2857217968584 - 9.0431912203186*I, -14.2437089383409 + 6.50486277982165*I, 25.0]])" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Matrix(H_eff_circ)" ] }, { "cell_type": "markdown", "id": "896553fa-3c28-44cc-89ad-5031d2f8d35f", "metadata": {}, "source": [ "Por último, podemos resolver el problema de valores propios generalizado para $\\tilde{H}$ :\n", "\n", "$\\tilde{H} \\vec{c} = c S \\vec{c}$\n", "\n", "y obtener una estimación de la energía del estado fundamental $c_{min}$\n", "\n" ] }, { "cell_type": "code", "execution_count": 58, "id": "8b997d15-ee25-40eb-80e8-1f054a298ee9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The estimated ground state energy is: 25.0\n", "The estimated ground state energy is: 22.572154819954875\n", "The estimated ground state energy is: 21.691509219286587\n", "The estimated ground state energy is: 21.23882298756386\n", "The estimated ground state energy is: 20.965499325470294\n" ] } ], "source": [ "gnd_en_circ_est_list = []\n", "for d in range(1, krylov_dim + 1):\n", " # Solve generalized eigenvalue problem for different size of the Krylov space\n", " gnd_en_circ_est = solve_regularized_gen_eig(\n", " H_eff_circ[:d, :d], S_circ[:d, :d], threshold=9e-1\n", " )\n", " gnd_en_circ_est_list.append(gnd_en_circ_est)\n", " print(\"The estimated ground state energy is: \", gnd_en_circ_est)" ] }, { "cell_type": "markdown", "id": "b39943f8-b75f-4bd4-b246-833ea4f9b710", "metadata": {}, "source": [ "Para un sector de una sola partícula, podemos calcular eficientemente el estado fundamental de este sector del Hamiltoniano clásicamente\n", "\n" ] }, { "cell_type": "code", "execution_count": 59, "id": "fcfe07e5-99e4-4276-a4d6-6f1c7e28c5f2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n_sys_qubits 30\n", "n_exc 1 , subspace dimension 31\n", "single particle ground state energy: 21.021912418526906\n" ] } ], "source": [ "gs_en = single_particle_gs(H_op, n_qubits)" ] }, { "cell_type": "code", "execution_count": 60, "id": "4bc52594-0376-497f-8a61-0949415a1fe0", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(\n", " range(1, krylov_dim + 1),\n", " gnd_en_circ_est_list,\n", " color=\"blue\",\n", " linestyle=\"-.\",\n", " label=\"KQD estimate\",\n", ")\n", "plt.plot(\n", " range(1, krylov_dim + 1),\n", " [gs_en] * krylov_dim,\n", " color=\"red\",\n", " linestyle=\"-\",\n", " label=\"exact\",\n", ")\n", "plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))\n", "plt.legend()\n", "plt.xlabel(\"Krylov space dimension\")\n", "plt.ylabel(\"Energy\")\n", "plt.title(\n", " \"Estimating Ground state energy with Krylov Quantum Diagonalization\"\n", ")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "2ef2fac9-b4b9-4457-a551-c98999a88710", "metadata": {}, "source": [ "\n", "\n", "## Apéndice: Subespacio de Krylov a partir de evoluciones en tiempo real\n", "\n", "El espacio unitario de Krylov se define como\n", "\n", "$$\n", "\\mathcal{K}_U(H, |\\psi\\rangle) = \\text{span}\\left\\{ |\\psi\\rangle, e^{-iH\\,dt} |\\psi\\rangle, \\dots, e^{-irH\\,dt} |\\psi\\rangle \\right\\}\n", "$$\n", "\n", "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\n", "\n", "$$\n", "\\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\\},\n", "$$\n", "\n", "es decir, donde todas las evoluciones temporales se desplazan hacia atrás $d$ timesteps.\n", "La razón por la que es indistinguible es porque los elementos de la matriz\n", "\n", "$$\n", "\\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\n", "$$\n", "\n", "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$.\n", "\n", "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\\]](#references) :\n", "\n", "**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)...\n", "\n", "1. $f(E_0)=1$\n", "2. $|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\n", "3. $f(E)$ es una combinación lineal de $e^{ijE\\,dt}$ para $j=-d,-d+1,...,d-1,d$\n", "\n", "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:\n", "\n", "$$\n", "|\\psi\\rangle = \\sum_{k=0}^{N}\\gamma_k|E_k\\rangle,\n", "$$\n", "\n", "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\n", "\n", "$$\n", "f(H)|\\psi\\rangle = \\sum_{k=0}^{N}\\gamma_kf(E_k)|E_k\\rangle,\n", "$$\n", "\n", "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\n", "\n", "$$\n", "\\text{energy error} = \\frac{\\langle\\psi|f(H)(H-E_0)f(H)|\\psi\\rangle}{\\langle\\psi|f(H)^2|\\psi\\rangle}\n", "$$\n", "\n", "$$\n", "= \\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}.\n", "$$\n", "\n", "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$ :\n", "\n", "$$\n", "\\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}.\n", "$$\n", "\n", "Podemos acotar el primer término en $\\delta$,\n", "\n", "$$\n", "\\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,\n", "$$\n", "\n", "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\n", "\n", "$$\n", "\\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).\n", "$$\n", "\n", "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, el límite superior de $E_k-E_0<2\\|H\\|$ y el límite superior de $\\sum_{E_k>E_0+\\delta}|\\gamma_k|^2<1$ dan como resultado\n", "\n", "$$\n", "\\text{energy error} \\le \\delta + \\frac{8}{|\\gamma_0|^2}\\|H\\|\\left(1 + \\delta\\right)^{-2d}.\n", "$$\n", "\n", "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\n", "\n", "## Apéndice: prueba de la reclamación 1\n", "\n", "Lo siguiente se deriva en su mayor parte de [\\[3\\]](#references), 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\n", "\n", "$$\n", "\\beta(a, b, d) = \\min_{p \\in \\Pi^*_d} \\max_{x \\in [a, b]} |p(x)| \\quad\n", "$$\n", "\n", "es\n", "\n", "$$\n", "p^*(x) = \\frac{T_d\\left(\\frac{b + a - 2x}{b - a}\\right)}{T_d\\left(\\frac{b + a}{b - a}\\right)}, \\quad\n", "$$\n", "\n", "y el valor mínimo correspondiente es\n", "\n", "$$\n", "\\beta(a, b, d) = T_d^{-1}\\left(\\frac{b + a}{b - a}\\right).\n", "$$\n", "\n", "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.\n", "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\n", "\n", "$$\n", "g(E) = \\frac{1-\\cos\\big((E-E_0)dt\\big)}{2},\n", "$$\n", "\n", "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$.\n", "\n", "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:\n", "\n", "$$\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)}\n", "$$\n", "\n", "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\n", "\n", "$$\n", "|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)\n", "$$\n", "\n", "$$\n", "\\leq 2\\left(1 + \\delta\\right)^{-d} = 2\\left(1 + \\delta\\right)^{-\\lfloor k/2\\rfloor}.\n", "$$\n", "\n" ] }, { "cell_type": "markdown", "id": "ccd6c21d-ed1e-4894-9cbb-8d3a9d91afe2", "metadata": {}, "source": [ "\n", "\n", "## Referencias\n", "\n", "\\[1] N. Yoshioka, M. Amico, W. Kirby et al. \"Diagonalization of large many-body Hamiltonians on a quantum processor\". [arXiv:2407.14431](https://arxiv.org/abs/2407.14431)\n", "\n", "\\[2] 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).\n", "\n", "\\[3] Å. Björck. \"Métodos numéricos en el cálculo de matrices\". Textos de Matemática Aplicada. Springer International Publishing. (2014).\n", "\n", "\\[4] William Kirby. \"Análisis de algoritmos cuánticos de Krylov con errores\". Quantum 8, 1457 (2024).\n", "\n" ] }, { "cell_type": "markdown", "id": "74948cc7-041f-412c-ab16-2554bf164061", "metadata": {}, "source": [ "\n", "\n", "## Encuesta tutorial\n", "\n", "Responda a esta breve encuesta para darnos su opinión sobre este tutorial. Su opinión nos ayudará a mejorar nuestra oferta de contenidos y la experiencia de los usuarios.\n", "\n", "[Enlace a la encuesta](https://your.feedback.ibm.com/jfe/form/SV_82nennpKIjjD8rQ)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" }, "hours": 1.5, "qpuSeconds": 1200 }, "nbformat": 4, "nbformat_minor": 4 }