Diagonalización cuántica de Krylov de hamiltonianos de red

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)

En segundo plano

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.

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 ). 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 ), 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 ), el algoritmo de diagonalización cuántica de Krylov (KQD). Este algoritmo funciona bien a gran escala [1 ] en el hardware cuántico existente, comparte garantías de rendimiento 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.

Requisitos

Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:

• Qiskit SDK v2.0 o posterior, con soporte de visualización
• Qiskit Runtime v0.22 o posterior ( pip install qiskit-ibm-runtime )

Configuración

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import itertools as it
import copy
from sympy import Matrix
import warnings

warnings.filterwarnings("ignore")

from qiskit.quantum_info import SparsePauliOp, Pauli, StabilizerState
from qiskit.circuit import Parameter, IfElseOp
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter
from qiskit.transpiler import Target, CouplingMap
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager


from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    EstimatorV2 as Estimator,
)


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 the Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of the 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]


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)

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

    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

Paso 1: Asignar entradas clásicas a un problema cuántico

El espacio de Krylov

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

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 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)}$.

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.

Algoritmo de diagonalización cuántica de Krylov

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

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}$

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

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

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.

Definir hamiltoniano

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$

# Define problem Hamiltonian.
n_qubits = 30
J = 1  # coupling strength for ZZ interaction

# 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, J) if term.count("Z") == 2 else (term, 1) for term in H_int]

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

Output:

[('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)]

Establecer los parámetros del algoritmo

Elegimos heurísticamente un valor para el paso temporal dt (basado en límites superiores de la norma hamiltoniana). Ref [2] 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.10833078115826875)

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.

# Set parameters for quantum Krylov algorithm
krylov_dim = 5  # size of Krylov subspace
num_trotter_steps = 6
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.

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 0x11eef9be0>

Prueba de Hadamard

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

... y luego medir:

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

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

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)

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 112753

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

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:

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$, 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. Entonces el circuito anterior implementa el siguiente estado antes de la medició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

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.

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

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

t = Parameter("t")

# Create instruction for rotation about XX+YY-ZZ:
Rxyz_circ = QuantumCircuit(2)
Rxyz_circ.rxx(t, 0, 1)
Rxyz_circ.ryy(t, 0, 1)
Rxyz_circ.rzz(t, 0, 1)
Rxyz_instr = Rxyz_circ.to_instruction(label="RXX+YY+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:

Utilice un circuito optimizado para la preparación del estado

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(2 * 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:  74

Hemos reducido considerablemente la profundidad de la prueba de Hadamard con una combinación de aproximación de Trotter y unitarios no controlados

Paso 3: Ejecutar utilizando Qiskit primitives

Instanciar el backend y establecer los parámetros de ejecución

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
if (
    "if_else" not in backend.target.operation_names
):  # Needed as "op_name" could be "if_else"
    backend.target.add_instruction(IfElseOp, name="if_else")
print(backend.name)

Transpiling a una QPU

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

target = backend.target
cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    t2 = target.qubit_properties[q].t2 * 1e6
    if meas_err > 0.02 or t2 < 100:
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except:
                    continue

for q in cmap_list:
    op_name = list(target.operation_names_for_qargs(q))[0]
    twoq_gate_err = target[f"{op_name}"][q].error
    if twoq_gate_err > 0.005:
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q)
                except:
                    continue


cust_cmap = CouplingMap(cust_cmap_list)
cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates,
    coupling_map=cust_cmap,
)

A continuación, transpile el circuito virtual a la mejor disposición física en este nuevo objetivo

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

qc_trans = pm.run(qc)

print("depth", qc_trans.depth(lambda x: x[0].num_qubits == 2))
print("num 2q ops", qc_trans.count_ops())
print(
    "physical qubits",
    sorted(
        [
            idx
            for idx, qb in qc_trans.layout.initial_layout.get_physical_bits().items()
            if qb._register.name != "ancilla"
        ]
    ),
)

Output:

depth 52
num 2q ops OrderedDict([('rz', 2058), ('sx', 1703), ('cz', 728), ('x', 84), ('barrier', 8)])
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]

Crear PUB para 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)

Circuitos de carrera

Los circuitos de $t=0$ se pueden calcular de forma clásica

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:

(25+0j)

Ejecuta circuitos para $S$ y $\tilde{H}$ con Estimator

# Experiment options
num_randomizations = 300
num_randomizations_learning = 30
shots_per_randomization = 100
noise_factors = [1, 1.2, 1.4]
learning_pair_depths = [0, 4, 24, 48]


experimental_opts = {}
experimental_opts["resilience"] = {
    "measure_mitigation": True,
    "measure_noise_learning": {
        "num_randomizations": num_randomizations_learning,
        "shots_per_randomization": shots_per_randomization,
    },
    "zne_mitigation": True,
    "zne": {"noise_factors": noise_factors},
    "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,
    },
    "zne": {
        "amplifier": "pea",
        "extrapolated_noise_factors": [0] + noise_factors,
    },
}
experimental_opts["twirling"] = {
    "num_randomizations": num_randomizations,
    "shots_per_randomization": shots_per_randomization,
    "strategy": "all",
}

estimator = Estimator(mode=backend, options=experimental_opts)


job = estimator.run([pub])

Paso 4: Procesamiento posterior y devolución del resultado en el formato clásico deseado

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])

Matrix(S_circ)

Output:

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

# 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 it.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])

Matrix(H_eff_circ)

Output:

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 for different size of the Krylov space
    gnd_en_circ_est = solve_regularized_gen_eig(
        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=9e-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:  25.0
The estimated ground state energy is:  22.572154819954875
The estimated ground state energy is:  21.691509219286587
The estimated ground state energy is:  21.23882298756386
The estimated ground state energy is:  20.965499325470294

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 30
n_exc 1 , subspace dimension 31
single particle ground state energy:  21.021912418526906

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:

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

El espacio unitario de Krylov se define como

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

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

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

1. $f(E_0)=1$
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
3. $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:

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

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

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$ :

Podemos acotar el primer término en $\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

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

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 energía más bajo en el espacio de Krylov está limitado 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 energía más bajo, 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: 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

es

y el valor mínimo correspondiente es

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

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:

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