Diagonalização quântica de Krylov de hamiltonianos de rede

Estimativa de uso: 20 minutos em um Heron r2 (OBSERVAÇÃO: essa é apenas uma estimativa. Seu tempo de execução pode variar)

Segundo plano

Este tutorial demonstra como implementar o Algoritmo de Diagonalização Quântica de Krylov (KQD) no contexto dos padrões Qiskit. Primeiro, você aprenderá sobre a teoria por trás do algoritmo e, em seguida, verá uma demonstração de sua execução em uma QPU.

Em todas as disciplinas, estamos interessados em aprender as propriedades do estado fundamental dos sistemas quânticos. Os exemplos incluem a compreensão da natureza fundamental das partículas e forças, a previsão e a compreensão do comportamento de materiais complexos e a compreensão das interações e reações bioquímicas. Devido ao crescimento exponencial do espaço de Hilbert e à correlação que surge em sistemas emaranhados, o algoritmo clássico tem dificuldade para resolver esse problema em sistemas quânticos de tamanho cada vez maior. Em uma extremidade do espectro está a abordagem existente que aproveita o foco do hardware quântico em métodos quânticos variacionais (por exemplo, eigensolver quântico variacional ). Essas técnicas enfrentam desafios com os dispositivos atuais devido ao alto número de chamadas de função necessárias no processo de otimização, o que adiciona uma grande sobrecarga de recursos quando são introduzidas técnicas avançadas de atenuação de erros, limitando assim sua eficácia a sistemas pequenos. No outro extremo do espectro, há métodos quânticos tolerantes a falhas com garantias de desempenho (por exemplo, estimativa de fase quântica ), que exigem circuitos profundos que só podem ser executados em um dispositivo tolerante a falhas. Por esses motivos, apresentamos aqui um algoritmo quântico baseado em métodos de subespaço (conforme descrito neste artigo de revisão ), o algoritmo de diagonalização quântica de Krylov (KQD). Esse algoritmo tem bom desempenho em larga escala [1] no hardware quântico existente, compartilha garantias de desempenho semelhantes às da estimativa de fase, é compatível com técnicas avançadas de atenuação de erros e pode fornecer resultados que são classicamente inacessíveis.

Requisitos

Antes de iniciar este tutorial, verifique se você tem os seguintes itens instalados:

• Qiskit SDK v2.0 ou posterior, com suporte para visualização
• Qiskit Runtime v0.22 ou mais tarde ( pip install qiskit-ibm-runtime )

Instalação

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

Passo 1: Mapear entradas clássicas para um problema quântico

O espaço de Krylov

O espaço de Krylov $\mathcal{K}^r$ de ordem $r$ é o espaço abrangido pelos vetores obtidos pela multiplicação das potências mais altas de uma matriz $A$, até $r-1$, com um vetor de referência $\vert v \rangle$.

Se a matriz $A$ for o Hamiltoniano $H$, nos referiremos ao espaço correspondente como o espaço de Krylov de potência $\mathcal{K}_P$. No caso em que $A$ é o operador de evolução temporal gerado pelo Hamiltoniano $U=e^{-iHt}$, vamos nos referir ao espaço como o espaço de Krylov unitário $\mathcal{K}_U$. O subespaço de Krylov de potência que usamos classicamente não pode ser gerado diretamente em um computador quântico, pois $H$ não é um operador unitário. Em vez disso, podemos usar o operador de evolução temporal $U = e^{-iHt}$, que pode ser mostrado como garantia de convergência semelhante ao método de potência. As potências de $U$ tornam-se, então, etapas de tempo diferentes $U^k = e^{-iH(kt)}$.

Consulte o Apêndice para obter uma derivação detalhada de como o espaço unitário de Krylov permite representar com precisão os estados próprios de baixa energia.

Algoritmo de diagonalização quântica de Krylov

Dado um Hamiltoniano $H$ que desejamos diagonalizar, primeiro consideramos o espaço de Krylov unitário correspondente $\mathcal{K}_U$. O objetivo é encontrar uma representação compacta do Hamiltoniano em $\mathcal{K}_U$, que chamaremos de $\tilde{H}$. Os elementos da matriz de $\tilde{H}$, a projeção do Hamiltoniano no espaço de Krylov, podem ser calculados pelos seguintes valores de expectativa

Onde $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ são os vetores do espaço unitário de Krylov e $t_n = n dt$ são os múltiplos da etapa de tempo $dt$ escolhida. Em um computador quântico, o cálculo de cada elemento da matriz pode ser feito com qualquer algoritmo que permita obter a sobreposição entre os estados quânticos. Este tutorial se concentra no teste de Hadamard. Como o $\mathcal{K}_U$ tem dimensão $r$, o Hamiltoniano projetado no subespaço terá dimensões $r \times r$. Com $r$ suficientemente pequeno (geralmente $r<<100$ é suficiente para obter a convergência das estimativas das energias próprias), podemos então diagonalizar facilmente o Hamiltoniano projetado $\tilde{H}$. Entretanto, não podemos diagonalizar diretamente $\tilde{H}$ devido à não ortogonalidade dos vetores do espaço de Krylov. Teremos que medir suas sobreposições e construir uma matriz $\tilde{S}$

Isso nos permite resolver o problema do valor próprio em um espaço não ortogonal (também chamado de problema do valor próprio generalizado)

Pode-se, então, obter estimativas dos autovalores e dos estados próprios de $H$ observando os de $\tilde{H}$. Por exemplo, a estimativa da energia do estado fundamental é obtida tomando-se o menor autovalor de $c$ e o estado fundamental do vetor próprio correspondente $\vec{c}$. Os coeficientes em $\vec{c}$ determinam a contribuição dos diferentes vetores que abrangem $\mathcal{K}_U$.

A figura mostra uma representação do circuito do teste de Hadamard modificado, um método usado para calcular a sobreposição entre diferentes estados quânticos. Para cada elemento da matriz $\tilde{H}_{i,j}$, é realizado um teste de Hadamard entre os estados $\vert \psi_i \rangle$, $\vert \psi_j \rangle$. Isso é destacado na figura pelo esquema de cores dos elementos da matriz e das operações $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ correspondentes. Assim, é necessário um conjunto de testes de Hadamard para todas as combinações possíveis de vetores do espaço de Krylov para calcular todos os elementos da matriz do Hamiltoniano projetado $\tilde{H}$. O fio superior no circuito de teste de Hadamard é um qubit ancilla que é medido na base X ou Y, e seu valor de expectativa determina o valor da sobreposição entre os estados. O fio inferior representa todos os qubits do sistema Hamiltoniano. A operação $\text{Prep} \; \psi_i$ prepara o qubit do sistema no estado $\vert \psi_i \rangle$ controlado pelo estado do qubit ancilla (da mesma forma para $\text{Prep} \; \psi_j$ ) e a operação $P$ representa a decomposição de Pauli do sistema Hamiltoniano $H = \sum_i P_i$. Uma derivação mais detalhada das operações calculadas pelo teste de Hadamard é apresentada abaixo.

Defina Hamiltoniano

Vamos considerar o Hamiltoniano de Heisenberg para $N$ qubits em uma cadeia linear: $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)]

Definir parâmetros para o algoritmo

Escolhemos heuristicamente um valor para a etapa de tempo dt (com base nos limites superiores da norma Hamiltoniana). A Ref [2] mostrou que um intervalo de tempo suficientemente pequeno é $\pi/\vert \vert H \vert \vert$, e que é preferível, até certo ponto, subestimar esse valor em vez de superestimá-lo, pois a superestimação pode permitir que as contribuições de estados de alta energia corrompam até mesmo o estado ideal no espaço de Krylov. Por outro lado, escolher $dt$ para ser muito pequeno leva a um condicionamento pior do subespaço de Krylov, já que os vetores da base de Krylov diferem menos de um intervalo de tempo para outro.

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

E definir outros parâmetros do algoritmo. Para fins deste tutorial, vamos nos limitar a usar um espaço de Krylov com apenas cinco dimensões, o que é bastante limitado.

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

Preparação do estado

Escolha um estado de referência $\vert \psi \rangle$ que tenha alguma sobreposição com o estado fundamental. Para esse Hamiltoniano, usamos o estado a com uma excitação no qubit do meio $\vert 00..010...00 \rangle$ como nosso estado de referência.

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

Output:

Evolução temporal

Podemos realizar o operador de evolução temporal gerado por um determinado Hamiltoniano: $U=e^{-iHt}$ por meio da aproximação 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>

teste de Hadamard

Onde $P$ é um dos termos na decomposição do Hamiltoniano $H=\sum P$ e $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ são operações controladas que preparam $|\psi_i\rangle$, $|\psi_j\rangle$ vetores do espaço unitário de Krylov, com $|\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$, primeiro aplique $H$...

... e depois medir:

A partir da identidade $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Da mesma forma, a medição de $Y$ resulta em

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

O circuito de teste Hadamard pode ser um circuito profundo quando decomposto em portas nativas (o que aumentará ainda mais se levarmos em conta a topologia do 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

Etapa 2: Otimizar o problema para execução em hardware quântico

Teste de Hadamard eficiente

Podemos otimizar os circuitos profundos para o teste de Hadamard que obtivemos introduzindo algumas aproximações e confiando em algumas suposições sobre o modelo Hamiltoniano. Por exemplo, considere o seguinte circuito para o teste de Hadamard:

Suponha que possamos calcular classicamente $E_0$, o valor próprio de $|0\rangle^N$ sob o Hamiltoniano $H$. Isso é satisfeito quando o Hamiltoniano preserva a simetria U(1). Embora isso possa parecer uma suposição forte, há muitos casos em que é seguro supor que há um estado de vácuo (nesse caso, ele é mapeado para o estado $|0\rangle^N$ ) que não é afetado pela ação do Hamiltoniano. Isso é verdadeiro, por exemplo, para Hamiltonianos químicos que descrevem moléculas estáveis (em que o número de elétrons é conservado). Considerando que a porta $\text{Prep} \; \psi$ prepara o estado de referência desejado $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, por exemplo, preparar o estado HF para a química $\text{Prep} \; \psi$ seria um produto de NOTs de um único qubit, portanto, controlado- $\text{Prep} \; \psi$ é apenas um produto de CNOTs. Então, o circuito acima implementa o seguinte estado antes da medição:

em que usamos o deslocamento de fase simulável clássico $U\ket{0}^N = e^{i\phi}\ket{0}^N$ na terceira linha. Portanto, os valores de expectativa são obtidos como

Usando essas premissas, conseguimos escrever os valores de expectativa dos operadores de interesse com menos operações controladas. Na verdade, só precisamos implementar a preparação do estado controlado $\text{Prep} \; \psi$ e não as evoluções de tempo controladas. Reenquadrar nosso cálculo como acima nos permitirá reduzir bastante a profundidade dos circuitos resultantes.

Decomponha o operador de evolução temporal com a decomposição de Trotter

Em vez de implementar exatamente o operador de evolução temporal, podemos usar a decomposição de Trotter para implementar uma aproximação dele. Repetir várias vezes uma determinada ordem de decomposição de Trotter nos proporciona uma redução adicional do erro introduzido pela aproximação. A seguir, criamos diretamente a implementação do Trotter da maneira mais eficiente para o gráfico de interação do Hamiltoniano que estamos considerando (somente interações de vizinhos mais próximos). Na prática, inserimos as rotações de Pauli $R_{xx}$, $R_{yy}$, $R_{zz}$ com um ângulo parametrizado $t$ que corresponde à implementação aproximada de $e^{-i (XX + YY + ZZ) t}$. Dada a diferença na definição das rotações de Pauli e a evolução temporal que estamos tentando implementar, teremos de usar o parâmetro $2*dt$ para obter uma evolução temporal de $dt$. Além disso, invertemos a ordem das operações para um número ímpar de repetições das etapas de Trotter, o que é funcionalmente equivalente, mas permite sintetizar operações adjacentes em uma única $SU(2)$ unitária. Isso proporciona um circuito muito mais raso do que o obtido com o uso da funcionalidade genérica do 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:

Use um circuito otimizado para a preparação do 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 o cálculo de elementos matriciais de $\tilde{S}$ e $\tilde{H}$ através do teste de Hadamard

A única diferença entre os circuitos usados no teste de Hadamard será a fase no operador de evolução temporal e os observáveis medidos. Portanto, podemos preparar um circuito modelo que represente o circuito genérico para o teste Hadamard, com espaços reservados para as portas que dependem do operador de evolução 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

Reduzimos consideravelmente a profundidade do teste de Hadamard com uma combinação de aproximação de Trotter e unidades não controladas

Passo 3: Execute usando Qiskit primitives

Instanciar o backend e definir os parâmetros de tempo de execução

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 para um QPU

Primeiro, vamos selecionar subconjuntos do mapa de acoplamento com qubits de “bom” desempenho (onde “bom” é bastante arbitrário aqui, queremos principalmente evitar qubits de desempenho realmente ruim) e criar um novo alvo para transpilagem

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

Em seguida, transpile o circuito virtual para o melhor layout físico nesse novo destino

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]

Criar PUBs para execução com o Estimador

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

Os circuitos para $t=0$ são calculáveis de forma clássica

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)

Executar circuitos para $S$ e $\tilde{H}$ com o 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])

Etapa 4: Pós-processamento e retorno do resultado no formato clássico desejado

results = job.result()[0]

Calcule as matrizes Hamiltoniana Eficaz e de Sobreposição

Primeiro, calcule a fase acumulada pelo estado $\vert 0 \rangle$ durante a evolução do tempo sem controle

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

Quando tivermos os resultados das execuções do circuito, poderemos pós-processar os dados para calcular os elementos da 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:

E os elementos da 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 fim, podemos resolver o problema do valor próprio generalizado para $\tilde{H}$ :

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

e obter uma estimativa da energia do 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 um setor de partícula única, podemos calcular com eficiência o estado fundamental desse setor do Hamiltoniano de forma clássica

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: Subespaço de Krylov a partir de evoluções em tempo real

O espaço unitário de Krylov é definido como

para algum intervalo de tempo $dt$ que determinaremos posteriormente. Suponha temporariamente que $r$ seja par: em seguida, defina $d=r/2$. Observe que, quando projetamos o Hamiltoniano no espaço de Krylov acima, ele é indistinguível do espaço de Krylov

ou seja, onde todas as evoluções temporais são deslocadas para trás em $d$ passos de tempo. O motivo pelo qual é indistinguível é que os elementos da matriz

são invariantes sob mudanças gerais do tempo de evolução, uma vez que as evoluções temporais são compatíveis com o Hamiltoniano. Para $r$ ímpares, podemos usar a análise para $r-1$.

Queremos mostrar que, em algum lugar desse espaço de Krylov, é garantido que haja um estado de baixa energia. Fazemos isso por meio do seguinte resultado, que é derivado do Teorema 3.1 em [3] :

Afirmação 1: existe uma função $f$ de modo que, para as energias $E$ na faixa espectral do Hamiltoniano (ou seja, entre a energia do estado fundamental e a energia máxima)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ para todos os valores de $E$ que se encontram a uma distância de $\ge\delta$ de $E_0$, ou seja, é exponencialmente suprimido
3. $f(E)$ é uma combinação linear de $e^{ijE\,dt}$ para $j=-d,-d+1,...,d-1,d$

Apresentamos uma prova abaixo, mas ela pode ser ignorada com segurança, a menos que se queira entender o argumento completo e rigoroso. Por enquanto, vamos nos concentrar nas implicações da afirmação acima. Pela propriedade 3 acima, podemos ver que o espaço de Krylov deslocado acima contém o estado $f(H)|\psi\rangle$. Esse é o nosso estado de baixa energia. Para entender o motivo, escreva $|\psi\rangle$ na base de energia própria:

em que $|E_k\rangle$ é o k-ésimo estado próprio de energia e $\gamma_k$ é sua amplitude no estado inicial $|\psi\rangle$. Expresso em termos disso, $f(H)|\psi\rangle$ é dado por

usando o fato de que podemos substituir $H$ por $E_k$ quando ele atua no estado próprio $|E_k\rangle$. O erro de energia desse estado é, portanto

Para transformar isso em um limite superior que seja mais fácil de entender, primeiro separamos a soma no numerador em termos com $E_k-E_0\le\delta$ e termos com $E_k-E_0>\delta$ :

Podemos limitar o primeiro termo por $\delta$,

em que a primeira etapa segue porque $E_k-E_0\le\delta$ para cada $E_k$ na soma, e a segunda etapa segue porque a soma no numerador é um subconjunto da soma no denominador. Para o segundo termo, primeiro reduzimos o denominador em $|\gamma_0|^2$, já que $f(E_0)^2=1$ : somando tudo, obtemos

Para simplificar o que resta, observe que, para todos esses $E_k$, pela definição de $f$, sabemos que $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Além disso, o limite superior de $E_k-E_0<2\|H\|$ e o limite superior de $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ fornecem

Isso é válido para qualquer $\delta>0$, portanto, se definirmos $\delta$ igual ao nosso erro de meta, então o limite de erro acima converge exponencialmente para esse valor com a dimensão de Krylov $2d=r$. Observe também que, se $\delta<E_1-E_0$, o termo $\delta$ desaparece totalmente no limite acima.

Para completar o argumento, primeiro observamos que o resultado acima é apenas o erro de energia do estado específico $f(H)|\psi\rangle$, em vez do erro de energia do estado de menor energia no espaço de Krylov. No entanto, pelo princípio variacional (Rayleigh-Ritz), o erro de energia do estado de menor energia no espaço de Krylov é limitado superiormente pelo erro de energia de qualquer estado no espaço de Krylov, de modo que o acima exposto também é um limite superior do erro de energia do estado de menor energia, ou seja, o resultado do algoritmo de diagonalização quântica de Krylov.

Uma análise semelhante à anterior pode ser realizada, levando em conta, adicionalmente, o ruído e o procedimento de limiar discutido no bloco de notas. Consulte [2] e [4] para ver essa análise.

Apêndice: prova da alegação 1

O seguinte é derivado principalmente de [3], Teorema 3.1:. Seja $0 < a < b$ e $\Pi^*_d$ o espaço de polinômios residuais (polinômios cujo valor em 0 é 1) de grau no máximo $d$. A solução para

é

e o valor mínimo correspondente é

Queremos converter isso em uma função que possa ser expressa naturalmente em termos de exponenciais complexos, porque essas são as evoluções de tempo reais que geram o espaço de Krylov quântico. Para isso, é conveniente introduzir a seguinte transformação de energias dentro do intervalo espectral do Hamiltoniano para números no intervalo $[0,1]$ : definir

onde $dt$ é um intervalo de tempo tal que $-\pi < E_0dt < E_\text{max}dt < \pi$. Observe que $g(E_0)=0$ e $g(E)$ crescem à medida que $E$ se afasta de $E_0$.

Agora, usando o polinômio $p^*(x)$ com os parâmetros a, b, d definidos como $a = g(E_0 + \delta)$, $b = 1$ e d = int( r/2 ), definimos a função:

em que $E_0$ é a energia do estado fundamental. Podemos ver, inserindo $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$, que $f(E)$ é um polinômio trigonométrico de grau $d$, ou seja, uma combinação linear de $e^{ijE\,dt}$ para $j=-d,-d+1,...,d-1,d$. Além disso, a partir da definição de $p^*(x)$ acima, temos que $f(E_0)=p(0)=1$ e para qualquer $E$ na faixa espectral tal que $\vert E-E_0 \vert > \delta$ temos