Algoritmos quânticos: algoritmos quânticos variacionais

nota

Takashi Imamichi (24 de maio de 2024)

Baixe o PDF da aula original. Note que alguns trechos de código podem estar desatualizados, pois são imagens estáticas.

O tempo aproximado de QPU para executar este experimento é de 9 minutos (testado em um processador Eagle).

(este notebook pode não ser avaliado no tempo permitido pelo Plano Aberto. Use os recursos de computação quântica com sabedoria.)

1. Introdução

Este tutorial apresenta uma visão geral de um algoritmo híbrido quântico-clássico, com foco específico no resolvedor variacional de autovalores (VQE) e no algoritmo de otimização aproximada quântica (QAOA). O objetivo principal desses algoritmos é resolver problemas de otimização utilizando circuitos quânticos com portas quânticas parametrizadas.

Apesar dos avanços na computação quântica, a presença de ruído nos dispositivos quânticos atuais torna difícil extrair resultados significativos de circuitos quânticos profundos. Para superar esse desafio, o VQE e o QAOA adotam uma abordagem híbrida quântico-clássica, que consiste em executar iterativamente circuitos quânticos relativamente curtos usando computação quântica e otimizar os parâmetros dos circuitos quânticos parametrizados alvo usando computação clássica.

O QAOA tem potencial para fornecer soluções ótimas para os problemas alvo em escala de utilidade, graças à aplicação de diversas técnicas de mitigação e supressão de erros. O VQE tem muitas aplicações (como química quântica) nas quais é menos escalável. Mas surgiram várias abordagens relacionadas a autovalores para complementar e ampliar o VQE, incluindo a diagonalização por subespaço de Krylov e a diagonalização quântica baseada em amostragem (SQD). Compreender o VQE é um primeiro passo importante para entender a ampla gama de algoritmos híbridos clássico-quânticos que surgiram.

Este módulo descreve os conceitos fundamentais e a implementação do VQE e do QAOA. Tutoriais adicionais explorarão tópicos avançados e técnicas para escalar esses algoritmos. Você precisa da seguinte biblioteca no seu ambiente para executar este notebook. Se ainda não instalou, pode fazê-lo descomentando e executando a célula a seguir.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-ibm-runtime rustworkx scipy

# % pip install 'qiskit[visualization]' qiskit-ibm-runtime

2. Calculando o autovalor mínimo de um Hamiltoniano simples

Vamos começar aplicando o VQE a um caso muito simples, para ver como ele funciona. Vamos calcular o autovalor mínimo da matriz de Pauli $Z$ com o VQE. Começaremos importando alguns pacotes gerais.

import numpy as np
from qiskit.circuit import ParameterVector, QuantumCircuit
from qiskit.primitives import StatevectorEstimator, StatevectorSampler
from qiskit.quantum_info import SparsePauliOp
from scipy.optimize import minimize

Agora definimos o operador de interesse e o visualizamos na forma matricial.

op = SparsePauliOp("Z")
op.to_matrix()

array([[ 1.+0.j,  0.+0.j],
       [ 0.+0.j, -1.+0.j]])

É fácil obter os autovalores de forma clássica, então podemos verificar nosso trabalho. Isso pode se tornar difícil à medida que escalamos em direção à utilidade. Aqui usamos o numpy.

# compute eigenvalues with numpy
result = np.linalg.eigh(op.to_matrix())
print("Eigenvalues:", result.eigenvalues)

Eigenvalues: [-1.  1.]

Para obter autovalores usando um algoritmo quântico variacional, construímos um circuito com portas que recebem parâmetros variacionais:

# define a variational form
param = ParameterVector("a", 3)
qc = QuantumCircuit(1, 1)
qc.u(param[0], param[1], param[2], 0)
qc_estimator = qc.copy()
qc.measure(0, 0)
qc.draw("mpl")

Saída da célula de código anterior

Se quisermos estimar o valor esperado de um operador (como $Z$ ), devemos usar o Estimator. Se quisermos observar os estados do sistema, usamos o Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Podemos calcular as contagens de cadeias de bits 0 e 1 com valores de parâmetros aleatórios [1, 2, 3] usando o Sampler.

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3])]).result()
counts = result[0].data.c.get_counts()
counts

{'0': 783, '1': 241}

Sabemos que podemos calcular o valor esperado de Z por $\langle Z \rangle = p_0 - p_1$ com probabilidades $\{0: p_0, 1: p_1\}$ .

# compute the expectation value of Z based on the counts
(counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())

0.529296875

Este circuito funcionou, mas os valores de parâmetros escolhidos não corresponderam a um estado de baixa energia (ou baixo autovalor). O autovalor obtido é consideravelmente maior que o mínimo. O resultado é similar quando se usa o estimator.

Note que o Estimator recebe circuitos quânticos sem medições.

result = estimator.run([(qc_estimator, op, [1, 2, 3])]).result()
result[0].data.evs

array(0.54030231)

Precisamos buscar pelos parâmetros que produzem o autovalor mais baixo. Criamos uma função para receber os valores de parâmetros da forma variacional e retornar o valor esperado $\langle Z \rangle$ .

# define a cost function to look for the minimum eigenvalue of Z
def cost(x):
    result = sampler.run([(qc, x)]).result()
    counts = result[0].data.c.get_counts()
    expval = (counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())
    # the following line shows the trajectory of the optimization
    print(expval, counts)
    return expval

Vamos aplicar a função minimize do SciPy para encontrar o autovalor mínimo de Z.

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0], method="COBYLA", tol=1e-8)
min_result

1.0 {'0': 1024}
0.494140625 {'0': 765, '1': 259}
0.466796875 {'0': 751, '1': 273}
0.564453125 {'0': 801, '1': 223}
-0.4296875 {'1': 732, '0': 292}
-0.984375 {'1': 1016, '0': 8}
-0.8984375 {'1': 972, '0': 52}
-0.990234375 {'1': 1019, '0': 5}
-0.892578125 {'1': 969, '0': 55}
-0.986328125 {'1': 1017, '0': 7}
-0.861328125 {'1': 953, '0': 71}
-1.0 {'1': 1024}
-0.982421875 {'1': 1015, '0': 9}
-0.99609375 {'1': 1022, '0': 2}
-0.986328125 {'1': 1017, '0': 7}
-1.0 {'1': 1024}
-0.990234375 {'1': 1019, '0': 5}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.994140625 {'1': 1021, '0': 3}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.0
       x: [ 3.182e+00  1.338e+00  1.664e-01]
    nfev: 63
   maxcv: 0.0

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()

{'0': 1, '1': 1023}

2.1 Exercício

Calcule o autovalor mínimo de $Z \otimes Z$ com o VQE.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

SparsePauliOp(['ZZ'],
              coeffs=[1.+0.j])
[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]

# compute eigenvalues with numpy

# define a variational form
# qc = ...

# compute counts of bitstrings with a random parameter values by Sampler
# result = sampler.run(...)
# result

# compute the expectation value of ZZ based on the counts

# verify the expectation value of ZZ with Estimator

# define a cost function to look for the minimum eigenvalue of ZZ
# def cost(x):
#    expval = ...
#    return expval

# minimize the cost function with scipy's minimize
# min_result = minimize(cost, [...], method="COBYLA", tol=1e-8)
# min_result

# check counts of bitstrings with the optimal parameter values
# result = sampler.run(qc, min_result.x).result()
# result

Soluções do exercício

Definimos o operador de interesse e o visualizamos na forma matricial.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

SparsePauliOp(['ZZ'],
              coeffs=[1.+0.j])
[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]

Para obter autovalores usando um algoritmo quântico variacional, construímos um circuito com portas que recebem parâmetros variacionais:

# define a variational form
param = ParameterVector("a", 6)
qc = QuantumCircuit(2, 2)
qc.u(param[0], param[1], param[2], 0)
qc.u(param[3], param[4], param[5], 1)
qc_estimator = qc.copy()
qc.measure([0, 1], [0, 1])
qc.draw("mpl")

Saída da célula de código anterior

Se quisermos estimar o valor esperado de um operador (como $Z \otimes Z$ ), usaríamos o Estimator. Se quisermos observar os estados do sistema, usamos o Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3, 4, 5, 6])]).result()
counts = result[0].data.c.get_counts()
counts

{'10': 661, '11': 203, '01': 47, '00': 113}

# compute the expectation value of ZZ based on the counts
(
    counts.get("00", 0)
    - counts.get("01", 0)
    - counts.get("10", 0)
    + counts.get("11", 0)
) / sum(counts.values())

-0.3828125

# verify the expectation value of ZZ with Estimator
result = estimator.run([(qc_estimator, z2, [1, 2, 3, 4, 5, 6])]).result()
result[0].data.evs

array(-0.35316516)

Precisamos buscar pelos parâmetros que produzem o autovalor mais baixo.

# define a cost function to look for the minimum eigenvalue of ZZ
def cost(x):
    result = sampler.run([(qc, x)]).result()
    counts = result[0].data.c.get_counts()
    expval = (
        counts.get("00", 0)
        - counts.get("01", 0)
        - counts.get("10", 0)
        + counts.get("11", 0)
    ) / sum(counts.values())
    print(expval, counts)
    return expval

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0, 0, 0, 0], method="COBYLA", tol=1e-8)
min_result

1.0 {'00': 1024}
0.578125 {'00': 808, '01': 216}
0.5234375 {'00': 780, '01': 244}
0.548828125 {'00': 793, '01': 231}
0.3515625 {'00': 637, '10': 164, '11': 55, '01': 168}
0.3359375 {'00': 638, '11': 46, '10': 174, '01': 166}
0.283203125 {'00': 602, '10': 181, '01': 186, '11': 55}
-0.087890625 {'01': 414, '00': 184, '10': 143, '11': 283}
0.236328125 {'10': 27, '11': 623, '01': 364, '00': 10}
-0.0625 {'11': 261, '01': 403, '00': 219, '10': 141}
0.248046875 {'01': 366, '11': 628, '00': 11, '10': 19}
-0.0625 {'10': 145, '11': 254, '01': 399, '00': 226}
0.228515625 {'01': 373, '11': 609, '00': 20, '10': 22}
0.0546875 {'11': 376, '10': 273, '01': 211, '00': 164}
-0.447265625 {'01': 731, '10': 10, '11': 267, '00': 16}
-0.71484375 {'01': 871, '11': 99, '00': 47, '10': 7}
-0.46484375 {'01': 741, '00': 253, '10': 9, '11': 21}
-0.87890625 {'01': 962, '00': 39, '11': 23}
-0.640625 {'00': 176, '01': 837, '11': 8, '10': 3}
-0.88671875 {'01': 966, '00': 41, '11': 17}
-0.994140625 {'01': 1021, '11': 3}
-0.91796875 {'01': 982, '11': 35, '00': 7}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.939453125 {'01': 993, '00': 31}
-0.990234375 {'01': 1019, '11': 5}
-0.90234375 {'01': 974, '00': 21, '11': 29}
-0.98046875 {'01': 1014, '11': 10}
-0.994140625 {'01': 1021, '00': 3}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.98828125 {'01': 1018, '11': 6}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '00': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '00': 1, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '11': 1, '00': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.994140625 {'01': 1021, '00': 3}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.998046875
       x: [ 3.167e+00  6.940e-01  1.033e+00 -2.894e-02  8.933e-01
            1.885e+00]
    nfev: 128
   maxcv: 0.0

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.99609375
       x: [ 3.098e+00 -5.402e-01  1.091e+00 -1.004e-02  3.615e-01
            6.913e-01]
    nfev: 115
   maxcv: 0.0

Obtivemos um autovalor extremamente próximo do mínimo fornecido pelo numpy.

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()

{'01': 1024}

3. Otimização Quântica com padrões do Qiskit

Nesta seção você vai aprender sobre padrões do Qiskit e otimização aproximada quântica. Um padrão do Qiskit é um conjunto intuitivo e repetível de etapas para implementar um fluxo de trabalho de computação quântica: "Qiskit function" Vamos aplicar esses padrões ao contexto de otimização combinatória e mostrar como resolver o problema Max-Cut usando o Quantum Approximate Optimization Algorithm (QAOA), um método iterativo híbrido (quântico-clássico).

Observe que esta parte sobre QAOA é baseada em "Part 1: Small-scale QAOA" do tutorial Quantum approximate optimization algorithm. Consulte o tutorial para aprender como escalar essa abordagem.

3.1 Padrão Qiskit (pequena escala) para otimização

Esta parte usará um problema Max-Cut de pequena escala para ilustrar as etapas necessárias para resolver um problema de otimização em um computador quântico. O problema Max-Cut é um problema de otimização difícil de resolver (mais especificamente, é um problema NP-difícil) com diversas aplicações em clustering, ciência de redes e física estatística. Este tutorial considera um grafo de nós conectados por arestas e tem como objetivo particionar os nós em dois conjuntos "cortando" arestas, de forma que o número de arestas cortadas seja maximizado. "Maxcut" Para contextualizar antes de mapear esse problema para um algoritmo quântico, você pode entender melhor como o problema Max-Cut se torna um problema de otimização combinatória clássico considerando primeiro a minimização de uma função $f(x)$

\min_{x\in \{0, 1\}^n}f(x),

onde a entrada $x$ é um vetor cujos componentes correspondem a cada nó de um grafo. Em seguida, restrinja cada um desses componentes a ser $0$ ou $1$ (que representam estar incluído ou não no corte). Este exemplo de pequena escala usa um grafo com $n=5$ nós.

Você pode escrever uma função para um par de nós $i,j$ que indica se a aresta correspondente $(i,j)$ está no corte. Por exemplo, a função $x_i + x_j - 2 x_i x_j$ é 1 somente se um dos valores $x_i$ ou $x_j$ for 1 (o que significa que a aresta está no corte) e zero caso contrário. O problema de maximizar as arestas no corte pode ser formulado como

\max_{x\in \{0, 1\}^n} \sum_{(i,j)} x_i + x_j - 2 x_i x_j,

que pode ser reescrito como uma minimização da forma

\min_{x\in \{0, 1\}^n} \sum_{(i,j)} 2 x_i x_j - x_i - x_j.

O mínimo de $f(x)$ neste caso ocorre quando o número de arestas cruzadas pelo corte é máximo. Como você pode ver, ainda não há nada relacionado à computação quântica. Você precisa reformular esse problema em algo que um computador quântico possa entender. Inicialize o problema criando um grafo com $n=5$ nós.

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 5

graph = rx.PyGraph()
graph.add_nodes_from(range(1, n + 1))
edge_list = [
    (0, 1, 1.0),
    (0, 2, 1.0),
    (1, 2, 1.0),
    (1, 3, 1.0),
    (2, 4, 1.0),
    (3, 4, 1.0),
]
graph.add_edges_from(edge_list)
pos = rx.spring_layout(graph, seed=2)
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str)

Saída da célula de código anterior

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

O primeiro passo do padrão é mapear o problema clássico (grafo) em circuitos e operadores quânticos. Para isso, há três etapas principais a seguir:

Utilizar uma série de reformulações matemáticas para representar esse problema usando a notação de problemas de Otimização Binária Irrestrita Quadrática (QUBO).
Reescrever o problema de otimização como um Hamiltoniano cujo estado fundamental corresponde à solução que minimiza a função de custo.
Criar um circuito quântico que preparará o estado fundamental desse Hamiltoniano por meio de um processo semelhante ao quantum annealing.

Observação: Na metodologia QAOA, o objetivo final é ter um operador (Hamiltoniano) que represente a função de custo do nosso algoritmo híbrido, bem como um circuito parametrizado (Ansatz) que represente estados quânticos com soluções candidatas para o problema. Você pode amostrar a partir desses estados candidatos e então avaliá-los usando a função de custo.

Grafo → problema de otimização

O primeiro passo do mapeamento é uma mudança de notação. O seguinte expressa o problema na notação QUBO:

\min_{x\in \{0, 1\}^n}x^T Q x,

onde $Q$ é uma matriz $n\times n$ de números reais, $n$ corresponde ao número de nós no seu grafo, $x$ é o vetor de variáveis binárias introduzido acima e $x^T$ indica a transposta do vetor $x$ .

Problem name: maxcut

Minimize
  2*x_1*x_2 + 2*x_1*x_3 + 2*x_2*x_3 + 2*x_2*x_4 + 2*x_3*x_5 + 2*x_4*x_5 - 2*x_1
  - 3*x_2 - 3*x_3 - 2*x_4 - 2*x_5

Subject to
  No constraints

  Binary variables (5)
    x_1 x_2 x_3 x_4 x_5

Problema de otimização → Hamiltoniano

Você pode então reformular o problema QUBO como um Hamiltoniano (aqui, uma matriz que representa a energia de um sistema):

H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.

Etapas de reformulação do problema QAOA para o Hamiltoniano

Para demonstrar como o problema QAOA pode ser reescrito dessa forma, primeiro substitua as variáveis binárias $x_i$ por um novo conjunto de variáveis $z_i\in\{-1, 1\}$ via

x_i = \frac{1-z_i}{2}.

Aqui você pode ver que se $x_i$ é $0$ , então $z_i$ deve ser $1$ . Quando os $x_i$ 's são substituídos pelos $z_i$ 's no problema de otimização ( $x^TQx$ ), uma formulação equivalente pode ser obtida.

x^TQx=\sum_{ij}Q_{ij}x_ix_j \\ =\frac{1}{4}\sum_{ij}Q_{ij}(1-z_i)(1-z_j) \\=\frac{1}{4}\sum_{ij}Q_{ij}z_iz_j-\frac{1}{4}\sum_{ij}(Q_{ij}+Q_{ji})z_i + \frac{n^2}{4}.

Agora, se definirmos $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$ , removermos o fator de prefixo e o termo constante $n^2$ , chegamos às duas formulações equivalentes do mesmo problema de otimização.

min_{x\in\{0,1\}^n} x^TQx\Longleftrightarrow \min_{z\in\{-1,1\}^n}z^TQz + b^Tz

Aqui, $b$ depende de $Q$ . Observe que para obter $z^TQz + b^Tz$ descartamos o fator 1/4 e uma constante de deslocamento $n^2$ , que não desempenham papel na otimização.

Agora, para obter uma formulação quântica do problema, promova as variáveis $z_i$ a uma matriz de Pauli $Z$ , como uma matriz $2\times 2$ da forma

Z_i = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}.

Quando você substitui essas matrizes no problema de otimização acima, obtém o seguinte Hamiltoniano

H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.

Lembre-se também de que as matrizes $Z$ são incorporadas no espaço computacional do computador quântico, ou seja, um espaço de Hilbert de tamanho $2^n\times 2^n$ . Portanto, você deve entender termos como $Z_iZ_j$ como o produto tensorial $Z_i\otimes Z_j$ incorporado no espaço de Hilbert $2^n\times 2^n$ . Por exemplo, em um problema com cinco variáveis de decisão, o termo $Z_1Z_3$ é entendido como $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ onde $I$ é a matriz identidade $2\times 2$ .

Este Hamiltoniano é chamado de Hamiltoniano da função de custo. Ele tem a propriedade de que seu estado fundamental corresponde à solução que minimiza a função de custo $f(x)$ . Portanto, para resolver seu problema de otimização, você agora precisa preparar o estado fundamental de $H_C$ (ou um estado com grande sobreposição com ele) no computador quântico. Em seguida, amostrar a partir desse estado irá, com alta probabilidade, fornecer a solução para $\min~f(x)$ .

def build_max_cut_operator(graph: rx.PyGraph) -> tuple[SparsePauliOp, float]:
    sp_list = []
    constant = 0
    for s, t in graph.edge_list():
        w = graph.get_edge_data(s, t)
        sp_list.append(("ZZ", [s, t], w / 2))
        constant -= 1 / 2
    return SparsePauliOp.from_sparse_list(
        sp_list, num_qubits=graph.num_nodes()
    ), constant

cost_hamiltonian, constant = build_max_cut_operator(graph)
print("Cost Function Hamiltonian:", cost_hamiltonian)
print("Constant:", constant)

Cost Function Hamiltonian: SparsePauliOp(['IIIZZ', 'IIZIZ', 'IIZZI', 'IZIZI', 'ZIZII', 'ZZIII'],
              coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])
Constant: -3.0

Hamiltoniano → circuito quântico

O Hamiltoniano $H_C$ contém a definição quântica do seu problema. Agora você pode criar um circuito quântico que ajudará a amostrar boas soluções do computador quântico. O QAOA é inspirado no quantum annealing e aplica camadas alternadas de operadores no circuito quântico.

A ideia geral é começar no estado fundamental de um sistema conhecido, $H^{\otimes n}|0\rangle$ acima, e então conduzir o sistema para o estado fundamental do operador de custo de interesse. Isso é feito aplicando os operadores $\exp\{-i\gamma_k H_C\}$ e $\exp\{-i\beta_k H_m\}$ com ângulos $\gamma_1,...,\gamma_p$ e $\beta_1,...,\beta_p~$ .

O circuito quântico gerado é parametrizado por $\gamma_i$ e $\beta_i$ , portanto você pode testar diferentes valores de $\gamma_i$ e $\beta_i$ e amostrar a partir do estado resultante. "QAOA circuit diagram" Neste caso, vamos usar um exemplo com 1 camada QAOA contendo dois parâmetros: $\gamma_1$ e $\beta_1$ .

from qiskit.circuit.library import QAOAAnsatz

circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")

Saída da célula de código anterior

circuit.decompose(reps=3).draw("mpl", fold=-1)

Saída da célula de código anterior

circuit.parameters

ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])

3.3 Passo 2. Otimizar circuitos para execução em hardware quântico

O circuito acima contém uma série de abstrações úteis para pensar sobre algoritmos quânticos, mas que não podem ser executadas diretamente no hardware. Para executar em uma QPU, o circuito precisa passar por uma série de operações que compõem a etapa de transpilação ou otimização de circuito do padrão.

A biblioteca Qiskit oferece uma série de passes de transpilação que atendem a uma ampla variedade de transformações de circuito. Você precisa garantir que seu circuito esteja otimizado para o seu objetivo.

A transpilação pode envolver várias etapas, como:

Mapeamento inicial dos qubits no circuito (como variáveis de decisão) para qubits físicos no dispositivo.
Expansão das instruções no circuito quântico para as instruções nativas do hardware que o backend entende.
Roteamento de quaisquer qubits no circuito que interagem para qubits físicos adjacentes entre si.
Supressão de erros pela adição de gates de qubit único para suprimir ruído com desacoplamento dinâmico.

Mais informações sobre transpilação estão disponíveis em nossa documentação.

O código a seguir transforma e otimiza o circuito abstrato em um formato pronto para execução em um dos dispositivos acessíveis pela nuvem usando o serviço Qiskit IBM® Runtime.

Observe que você pode testar seus programas localmente usando o "modo de teste local" antes de enviá-los para computadores quânticos reais. Mais informações sobre o modo de teste local estão disponíveis na documentação.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Use a quantum device
service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
# backend = service.backend("ibm_kingston")

# You can test your programs locally with a fake backend (local testing mode)
# backend = FakeBrisbane()

print(backend)

# Create pass manager for transpilation
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)

candidate_circuit = pm.run(circuit)
candidate_circuit.draw("mpl", fold=False, idle_wires=False)

service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

Saída da célula de código anterior

3.4 Passo 3. Executar usando primitivas do Qiskit

No fluxo de trabalho do QAOA, os parâmetros ótimos do QAOA são encontrados em um laço de otimização iterativo, que executa uma série de avaliações de circuito e usa um otimizador clássico para encontrar os parâmetros ótimos $\beta_k$ e $\gamma_k$ . Esse laço de execução é realizado pelas seguintes etapas:

Definir os parâmetros iniciais
Instanciar uma nova Session contendo o laço de otimização e a primitiva usada para amostrar o circuito
Uma vez encontrado o conjunto ótimo de parâmetros, executar o circuito uma última vez para obter uma distribuição final que será usada na etapa de pós-processamento.

Definir o circuito com parâmetros iniciais

Começamos com parâmetros escolhidos arbitrariamente.

initial_gamma = np.pi
initial_beta = np.pi / 2
init_params = [initial_gamma, initial_beta]

Definir o backend e a primitiva de execução

Use as primitivas do Qiskit Runtime para interagir com backends da IBM®. As duas primitivas são Sampler e Estimator, e a escolha da primitiva depende do tipo de medição que você deseja executar no computador quântico. Para a minimização de $H_C$ , use o Estimator, pois a medição da função de custo é simplesmente o valor esperado de $\langle H_C \rangle$ .

Executar

As primitivas oferecem uma variedade de modos de execução para agendar cargas de trabalho em dispositivos quânticos, e um fluxo de trabalho QAOA é executado de forma iterativa em uma sessão. "execution mode" Você pode conectar a função de custo baseada no Sampler à rotina de minimização do SciPy para encontrar os parâmetros ótimos.

def cost_func_estimator(params, ansatz, hamiltonian, estimator):
    # transform the observable defined on virtual qubits to
    # an observable defined on all physical qubits
    isa_hamiltonian = hamiltonian.apply_layout(ansatz.layout)

    pub = (ansatz, isa_hamiltonian, params)
    job = estimator.run([pub])

    results = job.result()[0]
    cost = results.data.evs

    objective_func_vals.append(cost)

    return cost

from qiskit_ibm_runtime import Session, EstimatorV2
from scipy.optimize import minimize

objective_func_vals = []  # Global variable
with Session(backend=backend) as session:
    # If using qiskit-ibm-runtime<0.24.0, change `mode=` to `session=`
    estimator = EstimatorV2(mode=session)
    estimator.options.default_shots = 1000

    # Set simple error suppression/mitigation options
    estimator.options.dynamical_decoupling.enable = True
    estimator.options.dynamical_decoupling.sequence_type = "XY4"
    estimator.options.twirling.enable_gates = True
    estimator.options.twirling.num_randomizations = "auto"

    result = minimize(
        cost_func_estimator,
        init_params,
        args=(candidate_circuit, cost_hamiltonian, estimator),
        method="COBYLA",
        tol=1e-2,
    )
    print(result)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.6557925874481715
       x: [ 2.873e+00  9.414e-01]
    nfev: 21
   maxcv: 0.0

O otimizador conseguiu reduzir o custo e encontrar melhores parâmetros para o circuito.

plt.figure(figsize=(12, 6))
plt.plot(objective_func_vals)
plt.xlabel("Iteration")
plt.ylabel("Cost")
plt.show()

Saída da célula de código anterior

Depois de encontrar os parâmetros ótimos para o circuito, você pode atribuir esses parâmetros e amostrar a distribuição final obtida com os parâmetros otimizados. É aqui que a primitiva Sampler deve ser usada, pois é a distribuição de probabilidade das medições de bitstrings que corresponde ao corte ótimo do grafo.

Observação: Isso significa preparar um estado quântico $\psi$ no computador e então medi-lo. Uma medição irá colapsar o estado para um único estado da base computacional — por exemplo, 010101110000... — que corresponde a uma solução candidata $x$ para o nosso problema de otimização inicial ( $\max f(x)$ ou $\min f(x)$ dependendo da tarefa).

optimized_circuit = candidate_circuit.assign_parameters(result.x)
optimized_circuit.draw("mpl", fold=False, idle_wires=False)

Saída da célula de código anterior

from qiskit_ibm_runtime import SamplerV2

# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `backend=`
sampler = SamplerV2(mode=backend)

# Set simple error suppression/mitigation options
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XY4"
sampler.options.twirling.enable_gates = True
sampler.options.twirling.num_randomizations = "auto"

pub = (optimized_circuit,)
job = sampler.run([pub], shots=int(1e4))
counts_int = job.result()[0].data.meas.get_int_counts()
counts_bin = job.result()[0].data.meas.get_counts()
shots = sum(counts_int.values())
final_distribution_int = {key: val / shots for key, val in counts_int.items()}
final_distribution_bin = {key: val / shots for key, val in counts_bin.items()}
print(final_distribution_int)

{12: 0.0652, 31: 0.0089, 4: 0.0085, 13: 0.0731, 26: 0.0256, 28: 0.0246, 17: 0.0405, 25: 0.0591, 20: 0.031, 15: 0.0221, 8: 0.017, 21: 0.0371, 14: 0.0461, 16: 0.0229, 19: 0.0723, 23: 0.0199, 22: 0.0478, 18: 0.0708, 24: 0.0165, 6: 0.0525, 7: 0.0155, 5: 0.0245, 3: 0.0231, 29: 0.0121, 30: 0.0062, 10: 0.0363, 1: 0.0097, 9: 0.042, 27: 0.0094, 11: 0.0349, 0: 0.0129, 2: 0.0119}

3.5 Passo 4. Pós-processar e retornar o resultado em formato clássico

A etapa de pós-processamento interpreta a saída de amostragem para retornar uma solução para o seu problema original. Neste caso, você está interessado no bitstring com a maior probabilidade, pois ele determina o corte ótimo. As simetrias no problema permitem quatro possíveis soluções, e o processo de amostragem retornará uma delas com uma probabilidade ligeiramente maior, mas você pode ver na distribuição plotada abaixo que quatro dos bitstrings são distintamente mais prováveis do que os demais.

# auxiliary functions to sample most likely bitstring
def to_bitstring(integer, num_bits):
    result = np.binary_repr(integer, width=num_bits)
    return [int(digit) for digit in result]

keys = list(final_distribution_int.keys())
values = list(final_distribution_int.values())
most_likely = keys[np.argmax(np.abs(values))]
most_likely_bitstring = to_bitstring(most_likely, len(graph))
most_likely_bitstring.reverse()

print("Result bitstring:", most_likely_bitstring)

Result bitstring: [1, 0, 1, 1, 0]

import matplotlib.pyplot as plt

matplotlib.rcParams.update({"font.size": 10})
final_bits = final_distribution_bin
values = np.abs(list(final_bits.values()))
top_4_values = sorted(values, reverse=True)[:4]
positions = []
for value in top_4_values:
    positions.append(np.where(values == value)[0])
fig = plt.figure(figsize=(11, 6))
ax = fig.add_subplot(1, 1, 1)
plt.xticks(rotation=45)
plt.title("Result Distribution")
plt.xlabel("Bitstrings (reversed)")
plt.ylabel("Probability")
ax.bar(list(final_bits.keys()), list(final_bits.values()), color="tab:grey")
for p in positions:
    ax.get_children()[p[0].item()].set_color("tab:purple")
plt.show()

Saída da célula de código anterior

Visualizar o melhor corte

A partir do bitstring ótimo, você pode visualizar esse corte no grafo original.

colors = ["tab:grey" if i == 0 else "tab:purple" for i in most_likely_bitstring]
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str, node_color=colors)

Saída da célula de código anterior

E calcular o valor do corte. A solução não é ótima devido ao ruído (o valor de corte da solução ótima é 5).

from typing import Sequence

def evaluate_sample(x: Sequence[int], graph: rx.PyGraph) -> float:
    assert len(x) == len(
        list(graph.nodes())
    ), "The length of x must coincide with the number of nodes in the graph."
    return sum(
        x[u] * (1 - x[v]) + x[v] * (1 - x[u]) for u, v in list(graph.edge_list())
    )

cut_value = evaluate_sample(most_likely_bitstring, graph)
print("The value of the cut is:", cut_value)

The value of the cut is: 5

Este tutorial pequena escala de QAOA está concluído. Você aprenderá como adaptar o QAOA em escala utilitária em "Part 2: scale it up!" do tutorial Quantum approximate optimization algorithm.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'

1. Introdução​

2. Calculando o autovalor mínimo de um Hamiltoniano simples​

2.1 Exercício​

Soluções do exercício​

3. Otimização Quântica com padrões do Qiskit​

3.1 Padrão Qiskit (pequena escala) para otimização​

3.2 Passo 1. Mapear entradas clássicas para um problema quântico​

Grafo → problema de otimização​

Problema de otimização → Hamiltoniano​

Hamiltoniano → circuito quântico​

3.3 Passo 2. Otimizar circuitos para execução em hardware quântico​

3.4 Passo 3. Executar usando primitivas do Qiskit​

Definir o circuito com parâmetros iniciais​

Definir o backend e a primitiva de execução​

Executar​

3.5 Passo 4. Pós-processar e retornar o resultado em formato clássico​

Visualizar o melhor corte​

1. Introdução

2. Calculando o autovalor mínimo de um Hamiltoniano simples

2.1 Exercício

Soluções do exercício

3. Otimização Quântica com padrões do Qiskit

3.1 Padrão Qiskit (pequena escala) para otimização

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

Grafo → problema de otimização

Problema de otimização → Hamiltoniano

Hamiltoniano → circuito quântico

3.3 Passo 2. Otimizar circuitos para execução em hardware quântico

3.4 Passo 3. Executar usando primitivas do Qiskit

Definir o circuito com parâmetros iniciais

Definir o backend e a primitiva de execução

Executar

3.5 Passo 4. Pós-processar e retornar o resultado em formato clássico

Visualizar o melhor corte