Algoritmos quânticos: Algoritmos quânticos variacionais

1. Introdução

Este tutorial fornece uma visão geral de um algoritmo híbrido quântico-clássico, concentrando-se especificamente no eigensolver quântico variacional (VQE) e no algoritmo de otimização aproximada quântica (QAOA). O principal objetivo desses algoritmos é lidar com problemas de otimização empregando 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 dificulta a extração de resultados significativos dos circuitos quânticos profundos. Para superar esse desafio, o VQE e o QAOA adotam uma abordagem híbrida quântica-clássica, que envolve a execução iterativa de circuitos quânticos relativamente curtos usando a computação quântica e a otimização dos parâmetros dos circuitos quânticos parametrizados de destino usando a computação clássica.

O QAOA tem o potencial de fornecer as soluções ideais para os problemas-alvo em escala de serviços públicos, graças à aplicação de várias técnicas de atenuação e supressão de erros. A VQE tem muitas aplicações (como a química quântica) nas quais é menos dimensionável. Mas várias abordagens relacionadas a valores próprios surgiram para complementar e aumentar a VQE, incluindo a diagonalização do subespaço de Krylov e a diagonalização quântica baseada em amostragem (SQD). Entender o VQE é uma primeira etapa importante para compreender 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. Outros tutoriais explorarão tópicos avançados e técnicas para aumentar o tamanho desses algoritmos.

Você precisa da seguinte biblioteca em seu ambiente para executar este notebook. Se ainda não o tiver instalado, você poderá instalá-lo retirando o comentário e executando a seguinte célula.

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

2. Cálculo do valor próprio mínimo de um hamiltoniano simples

Começaremos aplicando o VQE a um caso muito simples, para ver como ele funciona. Calcularemos o valor próprio mínimo da matriz Pauli $Z$ com 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 em forma de matriz.

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

Output:

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

É fácil obter os valores próprios de forma clássica, portanto, podemos verificar nosso trabalho. Isso pode se tornar difícil à medida que nos aproximamos da utilidade. Aqui usamos o numpy.

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

Output:

Eigenvalues: [-1.  1.]

Para obter valores próprios 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")

Output:

Se quisermos estimar o valor da expectativa de um operador (como $Z$ ), devemos usar o Estimator. Se quisermos examinar os estados do sistema, usaremos o Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Podemos calcular 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

Output:

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

Sabemos que podemos calcular o valor da expectativa de Z pelo site $\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())

Output:

0.529296875

Esse circuito funcionou, mas os valores dos parâmetros escolhidos não corresponderam a um estado de energia muito baixa (ou de baixo valor próprio). O valor próprio obtido é um pouco maior que o mínimo. O resultado é semelhante quando se usa o estimador.

Observe que o Estimator usa circuitos quânticos sem medições.

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

Output:

array(0.54030231)

Será necessário pesquisar os parâmetros e encontrar aqueles que produzem o menor valor próprio. Criamos uma função para receber os valores dos 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 SciPy's minimize para encontrar o valor próprio 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

Output:

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

Output:

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

2.1 Exercício

Calcule o valor próprio mínimo de $Z \otimes Z$ com VQE.

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

Output:

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 exibimos em forma de matriz.

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

Output:

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 valores próprios 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")

Output:

Se quisermos estimar o valor da expectativa de um operador (como $Z \otimes Z$ ), usaremos o Estimator. Se quisermos examinar os estados do sistema, usaremos 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

Output:

{'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())

Output:

-0.3828125

Esse circuito funcionou, mas os valores dos parâmetros escolhidos não corresponderam a um estado de energia muito baixa (ou de baixo valor próprio). O valor próprio obtido é um pouco maior que o mínimo. O resultado é semelhante quando se usa o estimador.

# 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

Output:

array(-0.35316516)

Será necessário pesquisar os parâmetros e encontrar aqueles que produzem o menor valor próprio.

# 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

Output:

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 valor próprio extremamente próximo do mínimo que nos foi fornecido pela numpy.

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

Output:

{'01': 1024}

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

Neste tutorial, aprenderemos sobre os padrões Qiskit e a otimização quântica aproximada. Um padrão Qiskit é um conjunto intuitivo e repetível de etapas para implementar um fluxo de trabalho de computação quântica:

Aplicaremos esses padrões ao contexto da otimização combinatória e mostraremos como resolver o problema do corte máximo utilizando o Algoritmo de Otimização Aproximada Quântica (QAOA), um método iterativo híbrido (quântico-clássico).

Observe que essa parte do QAOA é baseada na "Parte 1: QAOA em pequena escala" do tutorial do algoritmo de otimização aproximada quântica. Consulte o tutorial para saber como aumentar a escala.

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

Esta seção utilizará um problema de corte máximo em pequena escala para ilustrar as etapas necessárias para resolver um problema de otimização usando um computador quântico.

O problema do corte máximo é um problema de otimização difícil de resolver (mais especificamente, é um problema NP-difícil) com diversas aplicações em agrupamento, ciência das redes e física estatística. Este tutorial considera um grafo composto por nós conectados por arestas e tem como objetivo dividir os nós em dois conjuntos por meio do "corte" de arestas, de modo a maximizar o número de arestas cortadas.

Para contextualizar antes de traduzir esse problema em um algoritmo quântico, é possível compreender melhor como o problema do corte máximo se torna um problema clássico de otimização combinatória ao considerar, em primeiro lugar, a minimização de uma função $f(x)$

em que a entrada $x$ é um vetor cujos componentes correspondem a cada nó de um gráfico. Em seguida, restrinja cada um desses componentes para que sejam $0$ ou $1$ (que representam a inclusão ou não inclusão no corte). Este caso de exemplo em pequena escala usa um gráfico com $n=5$ nós.

Você poderia escrever uma função de um par de nós $i,j$ que indica se a borda 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 dois $x_i$ ou $x_j$ for 1 (o que significa que a borda está no corte) e zero caso contrário. O problema de maximizar as bordas no corte pode ser formulado como

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

O mínimo de $f(x)$ nesse caso é quando o número de bordas percorridas 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 seu problema criando um gráfico 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)

Output:

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

A primeira etapa do padrão é mapear o problema clássico (gráfico) em circuitos e operadores quânticos. Para isso, há três etapas principais a serem seguidas:

1. Utilize uma série de reformulações matemáticas para representar esse problema usando a notação de problemas Quadratic Unconstrained Binary Optimization (QUBO).
2. Reescreva o problema de otimização como um Hamiltoniano para o qual o estado fundamental corresponde à solução que minimiza a função de custo.
3. Crie um circuito quântico que preparará o estado fundamental desse Hamiltoniano por meio de um processo semelhante ao recozimento quântico.

Observação: Na metodologia QAOA, você deseja ter um operador (Hamiltoniano ) que represente a função de custo do nosso algoritmo híbrido, bem como um circuito parametrizado (Ansatz ) que represente os estados quânticos com soluções candidatas para o problema. Você pode fazer uma amostragem desses estados candidatos e, em seguida, avaliá-los usando a função de custo.

Gráfico → problema de otimização

A primeira etapa do mapeamento é uma mudança de notação. A seguir, expressamos o problema na notação QUBO:

onde $Q$ é uma matriz $n\times n$ de números reais, $n$ corresponde ao número de nós em seu gráfico, $x$ é o vetor de variáveis binárias apresentado acima e $x^T$ indica a transposição 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

Em seguida, você pode reformular o problema QUBO como um Hamiltoniano (aqui, uma matriz que representa a energia de um sistema):

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\}$ por meio de

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$ ), é possível obter uma formulação equivalente.

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

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

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

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

Lembre-se também de que as matrizes $Z$ estã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$.

Esse Hamiltoniano é chamado de função de custo Hamiltoniana. Ele tem a propriedade de que seu estado fundamental corresponde à solução que minimiza a função de custo $f(x)$. Portanto, para resolver o problema de otimização, agora é necessário preparar o estado fundamental de $H_C$ (ou um estado com alta sobreposição com ele) no computador quântico. Então, a amostragem desse estado produzirá, com alta probabilidade, 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)

Output:

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

Circuito quântico hamiltoniano

O Hamiltoniano $H_C$ contém a definição quântica de seu problema. Agora você pode criar um circuito quântico que ajudará a obter boas soluções do computador quântico. O QAOA é inspirado no recozimento quântico 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 depois direcionar o sistema para o estado fundamental do operador de custo no qual você está interessado. Isso é feito aplicando os operadores $\exp\{-i\gamma_k H_C\}$ e $\exp\{-i\beta_k H_m\}$ com os ângulos $\gamma_1,...,\gamma_p$ e $\beta_1,...,\beta_p~$.

O circuito quântico que você gera é parametrizado por $\gamma_i$ e $\beta_i$, de modo que você pode experimentar diferentes valores de $\gamma_i$ e $\beta_i$ e fazer uma amostra do estado resultante.

Neste caso, tentaremos um exemplo com 1 camada QAOA que contém 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")

Output:

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

Output:

circuit.parameters

Output:

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 em algoritmos quânticos, mas não é possível executá-los no hardware. Para poder ser executado em uma QPU, o circuito precisa ser submetido a uma série de operações que compõem a etapa de transpilação ou otimização do circuito do padrão.

A biblioteca Qiskit oferece uma série de passagens de transpilação que atendem a uma ampla gama de transformações de circuitos. Você precisa se certificar de que seu circuito seja otimizado para sua finalidade.

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.
• Desenrolamento das instruções no circuito quântico para as instruções nativas de hardware que o backend entende.
• Roteamento de quaisquer qubits no circuito que interagem com qubits físicos adjacentes uns aos outros.
• Supressão de erros com a adição de portas de um único qubit para suprimir o ruído com desacoplamento dinâmico.

Mais informações sobre a 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 ser executado em um dos dispositivos acessíveis pela nuvem usando o serviço Qiskit IBM® Runtime.

Observe que você pode testar seus programas localmente pelo "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)

Output:

  service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

3.4 Passo 3. Execute usando Qiskit primitives

No fluxo de trabalho do QAOA, os parâmetros ideais do QAOA são encontrados em um loop de otimização iterativo, que executa uma série de avaliações de circuitos e usa um otimizador clássico para encontrar os parâmetros ideais de $\beta_k$ e $\gamma_k$. Esse loop de execução é executado por meio das seguintes etapas:

1. Definir os parâmetros iniciais
2. Instanciar um novo site Session que contenha o loop de otimização e a primitiva usada para fazer a amostragem do circuito
3. Quando for encontrado um conjunto ideal de parâmetros, execute o circuito uma última vez para obter uma distribuição final que será usada na etapa de pós-processamento.

Defina o circuito com os 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]

Defina backend e primitiva de execução

Use as primitivas Qiskit Runtime para interagir com os backends IBM®. As duas primitivas são Sampler e Estimator, e a escolha da primitiva depende do tipo de medição que você deseja realizar no computador quântico. Para minimizar a função de custo $H_C$, utilize o Estimador, uma vez que a medida da função de custo é simplesmente a esperança de $\langle H_C \rangle$.

Executar o

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 iterativamente em uma sessão.

Você pode inserir a função de custo baseada no amostrador na rotina de minimização do SciPy para encontrar os parâmetros ideais.

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)

Output:

 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 foi capaz de reduzir o custo e encontrar parâmetros melhores para o circuito.

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

Output:

Depois de encontrar os parâmetros ideais para o circuito, você pode atribuir esses parâmetros e fazer uma amostragem da 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 bitstring que correspondem ao corte ideal do gráfico.

Observação: Isso significa preparar um estado quântico $\psi$ no computador e, em seguida, medi-lo. Uma medição colapsará o estado em um único estado de 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)

Output:

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)

Output:

{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-processamento, retornar resultado no formato clássico

A etapa de pós-processamento interpreta a saída da amostragem para retornar uma solução para seu problema original. Nesse caso, você está interessado no bitstring com a maior probabilidade, pois isso determina o corte ideal. As simetrias no problema permitem quatro soluções possíveis, e o processo de amostragem retornará uma delas com uma probabilidade um pouco maior, mas você pode ver na distribuição plotada abaixo que quatro das cadeias de bits são distintamente mais prováveis do que as 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)

Output:

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

Output:

Visualize o melhor corte

A partir da cadeia de bits ideal, você pode visualizar esse corte no gráfico 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)

Output:

E calcule o valor do corte. A solução não é ideal devido ao ruído (o valor de corte da solução ideal é 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)

Output:

The value of the cut is: 5

Isso conclui o tutorial de QAOA em pequena escala. Você aprenderá como adaptar o QAOA em uma escala de utilidade pública na "Parte 2: aumente a escala!" do tutorial do algoritmo de otimização aproximada Quantum.

# Check Qiskit version
import qiskit

qiskit.__version__

Output:

'2.0.2'