Funções de custo

Nesta lição, aprenderemos como avaliar uma função de custo:

• Primeiro, vamos conhecer as primitivas do Qiskit Runtime
• Definir uma função de custo $C(\vec\theta)$. Trata-se de uma função específica do problema que define o objetivo a ser minimizado (ou maximizado) pelo otimizador
• Definir uma estratégia de medição com as primitivas do Qiskit Runtime para equilibrar velocidade e precisão

 

Primitivas

Todos os sistemas físicos, sejam clássicos ou quânticos, podem existir em diferentes estados. Por exemplo, um carro em uma estrada pode ter determinada massa, posição, velocidade ou aceleração que caracterizam seu estado. Da mesma forma, sistemas quânticos também podem ter diferentes configurações ou estados, mas se diferenciam dos sistemas clássicos na maneira como lidamos com medições e evolução de estados. Isso leva a propriedades únicas, como superposição e emaranhamento, exclusivas da mecânica quântica. Assim como podemos descrever o estado de um carro usando propriedades físicas como velocidade ou aceleração, também podemos descrever o estado de um sistema quântico usando observáveis, que são objetos matemáticos.

Na mecânica quântica, os estados são representados por vetores coluna complexos normalizados, ou kets ($|\psi\rangle$), e os observáveis são operadores lineares hermitianos ($\hat{H}=\hat{H}^{\dagger}$) que atuam sobre os kets. Um autovetor ($|\lambda\rangle$) de um observável é chamado de autoestado. Ao medir um observável em um de seus autoestados ($|\lambda\rangle$), obtemos o autovalor correspondente ($\lambda$) como resultado.

Se você estiver se perguntando como medir um sistema quântico e o que é possível medir, o Qiskit oferece duas primitivas que podem ajudar:

• Sampler: Dado um estado quântico $|\psi\rangle$, essa primitiva obtém a probabilidade de cada possível estado da base computacional.
• Estimator: Dado um observável quântico $\hat{H}$ e um estado $|\psi\rangle$, essa primitiva calcula o valor esperado de $\hat{H}$.

A primitiva Sampler

A primitiva Sampler calcula a probabilidade de obter cada estado possível $|k\rangle$ da base computacional, dado um circuito quântico que prepara o estado $|\psi\rangle$. Ela calcula

$$p_k = |\langle k | \psi \rangle|^2 \quad \forall k \in \mathbb{Z}_2^n \equiv \{0,1,\cdots,2^n-1\},$$

onde $n$ é o número de qubits e $k$ é a representação inteira de qualquer string binária de saída possível $\{0,1\}^n$ (ou seja, inteiros na base $2$).

O Sampler do Qiskit Runtime executa o circuito várias vezes em um dispositivo quântico, realizando medições a cada execução e reconstruindo a distribuição de probabilidade a partir das bitstrings obtidas. Quanto mais execuções (ou shots) forem realizadas, mais precisos serão os resultados, mas isso exige mais tempo e recursos quânticos.

No entanto, como o número de saídas possíveis cresce exponencialmente com o número de qubits $n$ (ou seja, $2^n$), o número de shots também precisará crescer exponencialmente para capturar uma distribuição de probabilidade densa. Por isso, o Sampler é eficiente apenas para distribuições de probabilidade esparsas, nas quais o estado alvo $|\psi\rangle$ deve ser exprimível como uma combinação linear de estados da base computacional, com o número de termos crescendo no máximo polinomialmente com o número de qubits:

$$|\psi\rangle = \sum^{\text{Poly}(n)}_k w_k |k\rangle.$$

O Sampler também pode ser configurado para recuperar probabilidades de uma subseção do circuito, representando um subconjunto dos estados totais possíveis.

A primitiva Estimator

A primitiva Estimator calcula o valor esperado de um observável $\hat{H}$ para um estado quântico $|\psi\rangle$; onde as probabilidades dos observáveis podem ser expressas como $p_\lambda = |\langle\lambda|\psi\rangle|^2$, sendo $|\lambda\rangle$ os autoestados do observável $\hat{H}$. O valor esperado é então definido como a média de todos os resultados possíveis $\lambda$ (ou seja, os autovalores do observável) de uma medição do estado $|\psi\rangle$, ponderados pelas probabilidades correspondentes:

$$\langle\hat{H}\rangle_\psi := \sum_\lambda p_\lambda \lambda = \langle \psi | \hat{H} | \psi \rangle$$

Porém, calcular o valor esperado de um observável nem sempre é possível, pois frequentemente não conhecemos sua eigenbasis. O Estimator do Qiskit Runtime utiliza um processo algébrico sofisticado para estimar o valor esperado em um dispositivo quântico real, decompondo o observável em uma combinação de outros observáveis cuja eigenbasis já conhecemos.

Em termos mais simples, o Estimator decompõe qualquer observável que não sabe como medir em observáveis mais simples e mensuráveis, chamados de operadores de Pauli. Qualquer operador pode ser expresso como uma combinação de $4^n$ operadores de Pauli.

$$\hat{P}_k := \sigma_{k_{n-1}}\otimes \cdots \otimes \sigma_{k_0} \quad \forall k \in \mathbb{Z}_4^n \equiv \{0,1,\cdots,4^n-1\}, \\$$

de forma que

$$\hat{H} = \sum^{4^n-1}_{k=0} w_k \hat{P}_k$$

onde $n$ é o número de qubits, $k \equiv k_{n-1} \cdots k_0$ para $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (ou seja, inteiros na base $4$), e $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

Após realizar essa decomposição, o Estimator deriva um novo circuito $V_k|\psi\rangle$ para cada observável $\hat{P}_k$ (a partir do circuito original), a fim de efetivamente diagonalizar o observável de Pauli na base computacional e medi-lo. Podemos medir observáveis de Pauli com facilidade porque conhecemos $V_k$ antecipadamente, o que geralmente não acontece para outros observáveis.

Para cada $\hat{P}_{k}$, o Estimator executa o circuito correspondente em um dispositivo quântico várias vezes, mede o estado de saída na base computacional e calcula a probabilidade $p_{kj}$ de obter cada saída possível $j$. Em seguida, ele busca o autovalor $\lambda_{kj}$ de $P_k$ correspondente a cada saída $j$, multiplica por $w_k$ e soma todos os resultados para obter o valor esperado do observável $\hat{H}$ para o estado $|\psi\rangle$ dado.

$$\langle\hat{H}\rangle_\psi = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj},$$

Como calcular o valor esperado de $4^n$ Paulis é impraticável (ou seja, cresce exponencialmente), o Estimator só é eficiente quando uma grande quantidade de $w_k$ é zero (ou seja, decomposição de Pauli esparsa em vez de densa). Formalmente, dizemos que, para que essa computação seja eficientemente resolúvel, o número de termos não nulos deve crescer no máximo polinomialmente com o número de qubits $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

O leitor pode notar a suposição implícita de que a amostragem de probabilidade também precisa ser eficiente, como explicado para o Sampler, o que significa

$$\langle\hat{H}\rangle_\psi = \sum_{k}^{\text{Poly}(n)} w_k \sum_{j}^{\text{Poly}(n)}p_{kj} \lambda_{kj}.$$

Exemplo guiado para calcular valores esperados

Vamos assumir o estado de um único qubit $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ e o observável

$$\begin{aligned} \hat{H} & = \begin{pmatrix} -1 & 2 \\ 2 & 1 \\ \end{pmatrix}\\[1mm] & = 2X - Z \end{aligned}$$

com o seguinte valor esperado teórico $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Como não sabemos como medir esse observável, não podemos calcular seu valor esperado diretamente e precisamos reescrevê-lo como $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. É possível mostrar que isso resulta no mesmo valor, notando que $\langle+|X|+\rangle = 1$ e $\langle+|Z|+\rangle = 0$.

Vamos ver como calcular $\langle X \rangle_+$ e $\langle Z \rangle_+$ diretamente. Como $X$ e $Z$ não comutam (ou seja, não compartilham a mesma eigenbasis), eles não podem ser medidos simultaneamente, portanto precisamos dos circuitos auxiliares:

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

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

# The following code will work for any other initial single-qubit state and observable
original_circuit = QuantumCircuit(1)
original_circuit.h(0)

H = SparsePauliOp(["X", "Z"], [2, -1])

aux_circuits = []
for pauli in H.paulis:
    aux_circ = original_circuit.copy()
    aux_circ.barrier()
    if str(pauli) == "X":
        aux_circ.h(0)
    elif str(pauli) == "Y":
        aux_circ.sdg(0)
        aux_circ.h(0)
    else:
        aux_circ.id(0)
    aux_circ.measure_all()
    aux_circuits.append(aux_circ)

original_circuit.draw("mpl")

# Auxiliary circuit for X
aux_circuits[0].draw("mpl")

# Auxiliary circuit for Z
aux_circuits[1].draw("mpl")

Agora podemos realizar o cálculo manualmente usando o Sampler e verificar os resultados com o Estimator:

from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.result import QuasiDistribution
import numpy as np

## SAMPLER
shots = 10000
sampler = StatevectorSampler()
job = sampler.run(aux_circuits, shots=shots)

# Run the sampler job and step through results
expvals = []
for index, pauli in enumerate(H.paulis):
    data_pub = job.result()[index].data
    bitstrings = data_pub.meas.get_bitstrings()
    counts = data_pub.meas.get_counts()
    quasi_dist = QuasiDistribution(
        {outcome: freq / shots for outcome, freq in counts.items()}
    )

    # Use the probabilities and known eigenvalues of Pauli operators to estimate the expectation value.
    val = 0

    if str(pauli) == "X":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Y":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Z":
        val += 1 * quasi_dist.get(0, 0)
        val += -1 * quasi_dist.get(1, 0)

    expvals.append(val)

# Print expectation values

print("Sampler results:")
for pauli, expval in zip(H.paulis, expvals):
    print(f"  >> Expected value of {str(pauli)}: {expval:.5f}")

total_expval = np.sum(H.coeffs * expvals).real
print(f"  >> Total expected value: {total_expval:.5f}")

# Use estimator for comparison
observables = [
    *H.paulis,
    H,
]  # Note: run for individual Paulis as well as full observable H

estimator = StatevectorEstimator()
job = estimator.run([(original_circuit, observables)])
estimator_expvals = job.result()[0].data.evs

# Print results
print("Estimator results:")
for obs, expval in zip(observables, estimator_expvals):
    if obs is not H:
        print(f"  >> Expected value of {str(obs)}: {expval:.5f}")
    else:
        print(f"  >> Total expected value: {expval:.5f}")

Sampler results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00420
  >> Total expected value: 1.99580
Estimator results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00000
  >> Total expected value: 2.00000

Rigor matemático (opcional)

Expressando $|\psi\rangle$ na base dos autoestados de $\hat{H}$, $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, segue-se:

$$\begin{aligned} \langle \psi | \hat{H} | \psi \rangle & = \bigg(\sum_{\lambda'}a^*_{\lambda'} \langle \lambda'|\bigg) \hat{H} \bigg(\sum_{\lambda} a_\lambda | \lambda\rangle\bigg)\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \langle \lambda'|\hat{H}| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \langle \lambda'| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \cdot \delta_{\lambda, \lambda'}\\[1mm] & = \sum_\lambda |a_\lambda|^2 \lambda\\[1mm] & = \sum_\lambda p_\lambda \lambda\\[1mm] \end{aligned}$$

Como não conhecemos os autovalores nem os autoestados do observável alvo $\hat{H}$, primeiro precisamos considerar sua diagonalização. Dado que $\hat{H}$ é Hermitiano, existe uma transformação unitária $V$ tal que $\hat{H}=V^\dagger \Lambda V,$ onde $\Lambda$ é a matriz diagonal de autovalores, de modo que $\langle j | \Lambda | k \rangle = 0$ se $j\neq k$ e $\langle j | \Lambda | j \rangle = \lambda_j$.

Isso implica que o valor esperado pode ser reescrito como:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \langle\psi|V^\dagger \Lambda V|\psi\rangle\\[1mm] & = \langle\psi|V^\dagger \bigg(\sum_{j=0}^{2^n-1} |j\rangle \langle j|\bigg) \Lambda \bigg(\sum_{k=0}^{2^n-1} |k\rangle \langle k|\bigg) V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |k\rangle \langle k| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |j\rangle \langle j| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}|\langle j| V|\psi\rangle|^2 \lambda_j\\[1mm] \end{aligned}$$

Dado que, se um sistema está no estado $|\phi\rangle = V |\psi\rangle$, a probabilidade de medir $| j\rangle$ é $p_j = |\langle j|\phi \rangle|^2$, o valor esperado acima pode ser expresso como:

$$\langle\psi|\hat{H}|\psi\rangle = \sum_{j=0}^{2^n-1} p_j \lambda_j.$$

É muito importante observar que as probabilidades são obtidas do estado $V |\psi\rangle$ e não de $|\psi\rangle$. É por isso que a matriz $V$ é absolutamente necessária. Você pode estar se perguntando como obter a matriz $V$ e os autovalores $\Lambda$. Se você já tivesse os autovalores, não haveria necessidade de usar um computador quântico, uma vez que o objetivo dos algoritmos variacionais é justamente encontrar esses autovalores de $\hat{H}$.

Felizmente, existe uma alternativa: qualquer matriz $2^n \times 2^n$ pode ser escrita como uma combinação linear de $4^n$ produtos tensoriais de $n$ matrizes de Pauli e identidades, todas hermitianas e unitárias com $V$ e $\Lambda$ conhecidos. É isso que o Estimator do Runtime faz internamente, decompondo qualquer objeto Operator em um SparsePauliOp.

Aqui estão os Operadores que podem ser usados:

$$\begin{array}{c|c|c|c} \text{Operator} & \sigma & V & \Lambda \\[1mm] \hline I & \sigma_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & V_0 = I & \Lambda_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\[4mm] X & \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & V_1 = H =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} & \Lambda_1 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Y & \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & V_2 = HS^\dagger =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\cdot \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -i \\ 1 & i \end{pmatrix}\quad & \Lambda_2 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Z & \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & V_3 = I & \Lambda_3 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{array}$$

Então vamos reescrever $\hat{H}$ em termos dos Paulis e identidades:

$$\hat{H} = \sum_{k_{n-1}=0}^3... \sum_{k_0=0}^3 w_{k_{n-1}...k_0} \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0} = \sum_{k=0}^{4^n-1} w_k \hat{P}_k,$$

onde $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ para $k_{n-1},...,k_0\in \{0,1,2,3\}$ (ou seja, base $4$) e $\hat{P}_{k} := \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0}$:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}|\langle j| V_k|\psi\rangle|^2 \langle j| \Lambda_k |j\rangle \\[1mm] & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj}, \\[1mm] \end{aligned}$$

onde $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ e $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, de forma que: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Funções de custo

Em geral, funções de custo são usadas para descrever o objetivo de um problema e avaliar o quão bem um estado candidato está se saindo em relação a esse objetivo. Essa definição pode ser aplicada a diversos exemplos em química, aprendizado de máquina, finanças, otimização, entre outros.

Vamos considerar um exemplo simples de encontrar o estado fundamental de um sistema. Nosso objetivo é minimizar o valor esperado do observável que representa a energia (Hamiltoniano $\hat{\mathcal{H}}$):

$$\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle$$

Podemos usar o Estimator para calcular o valor esperado e passar esse valor a um otimizador para minimizá-lo. Se a otimização for bem-sucedida, ela retornará um conjunto de valores ótimos de parâmetros $\vec\theta^*$, a partir dos quais poderemos construir o estado solução proposto $|\psi(\vec\theta^*)\rangle$ e calcular o valor esperado observado como $C(\vec\theta^*)$.

Observe que só conseguiremos minimizar a função de custo para o conjunto limitado de estados que estamos considerando. Isso nos leva a duas possibilidades distintas:

• Nosso ansatz não define o estado solução ao longo do espaço de busca: Se for esse o caso, o otimizador nunca encontrará a solução, e precisamos experimentar outros ansatzes que possam representar nosso espaço de busca com mais precisão.
• Nosso otimizador não consegue encontrar essa solução válida: A otimização pode ser definida de forma global ou local. Exploraremos o que isso significa na seção seguinte.

No geral, realizaremos um loop de otimização clássica, mas dependendo da avaliação da função de custo em um computador quântico. Nessa perspectiva, pode-se pensar na otimização como um processo puramente clássico em que chamamos um oráculo quântico de caixa-preta cada vez que o otimizador precisa avaliar a função de custo.

def cost_func_vqe(params, circuit, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (circuit, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    return cost

from qiskit.circuit.library import TwoLocal

observable = SparsePauliOp.from_list([("XX", 1), ("YY", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)
ansatz = reference_circuit.compose(variational_form)

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
ansatz.decompose().draw("mpl")

Primeiro, executaremos isso usando um simulador: o StatevectorEstimator. Isso geralmente é recomendado para depuração, mas em seguida faremos o cálculo em hardware quântico real. Cada vez mais, os problemas de interesse deixam de ser classicamente simuláveis sem instalações de supercomputação de ponta.

estimator = StatevectorEstimator()
cost = cost_func_vqe(theta_list, ansatz, observable, estimator)
print(cost)

[-0.58744589]

Agora procederemos com a execução em um computador quântico real. Observe as mudanças de sintaxe. As etapas envolvendo o pass_manager serão discutidas com mais detalhes no próximo exemplo. Uma etapa de particular importância em algoritmos variacionais é o uso de uma sessão do Qiskit Runtime. Iniciar uma sessão permite executar várias iterações de um algoritmo variacional sem precisar aguardar em uma nova fila cada vez que os parâmetros são atualizados. Isso é importante quando os tempos de fila são longos e/ou muitas iterações são necessárias. Apenas parceiros da IBM Quantum® Network podem usar sessões Runtime. Se você não tiver acesso a sessões, pode reduzir o número de iterações enviadas de uma vez e salvar os parâmetros mais recentes para uso em execuções futuras. Se você enviar iterações demais ou encontrar tempos de fila muito longos, poderá receber o código de erro 1217, que se refere a longos atrasos entre envios de jobs.

# Estimated usage: < 1 min. Benchmarked at 7 seconds on an Eagle processor
# Load necessary packages:

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Select the least busy backend:

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)
# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(theta_list, isa_ansatz, isa_observable, estimator)

session.close()
print(cost)

Observe que os valores obtidos nos dois cálculos acima são bem similares. Técnicas para melhorar os resultados serão discutidas mais adiante.

Exemplo de mapeamento para sistemas não físicos

O problema de corte máximo (Max-Cut) é um problema de otimização combinatória que consiste em dividir os vértices de um grafo em dois conjuntos disjuntos, de modo que o número de arestas entre os dois conjuntos seja maximizado. Mais formalmente, dado um grafo não direcionado $G=(V,E)$, onde $V$ é o conjunto de vértices e $E$ é o conjunto de arestas, o problema Max-Cut pede para particionar os vértices em dois subconjuntos disjuntos, $S$ e $T$, de forma que o número de arestas com um extremo em $S$ e o outro em $T$ seja maximizado.

Podemos aplicar Max-Cut para resolver diversos problemas, incluindo: agrupamento (clustering), design de redes, transições de fase, entre outros. Começaremos criando um grafo para o problema:

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

Este problema pode ser expresso como um problema de otimização binária. Para cada nó $0 \leq i < n$, onde $n$ é o número de nós do grafo (neste caso $n=4$), consideraremos a variável binária $x_i$. Essa variável terá o valor $1$ se o nó $i$ estiver em um dos grupos que chamaremos de $1$, e $0$ se estiver no outro grupo, que chamaremos de $0$. Também denotaremos como $w_{ij}$ (elemento $(i,j)$ da matriz de adjacência $w$) o peso da aresta que vai do nó $i$ ao nó $j$. Como o grafo é não direcionado, $w_{ij}=w_{ji}$. Assim, podemos formular nosso problema como a maximização da seguinte função de custo:

$$\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}$$

Para resolver esse problema com um computador quântico, vamos expressar a função de custo como o valor esperado de um observável. Contudo, os observáveis aceitos nativamente pelo Qiskit são compostos de operadores de Pauli, que têm autovalores $1$ e $-1$ em vez de $0$ e $1$. Por isso, faremos a seguinte mudança de variável:

Onde $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Podemos usar a matriz de adjacência $w$ para acessar convenientemente os pesos de todas as arestas. Isso será usado para obter nossa função de custo:

$$z_i = 1-2x_i \rightarrow x_i = \frac{1-z_i}{2}$$

Isso implica que:

$$\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}$$

Portanto, a nova função de custo que queremos maximizar é:

$$\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}$$

Além disso, a tendência natural de um computador quântico é encontrar mínimos (geralmente a menor energia) em vez de máximos. Por isso, em vez de maximizar $C(\vec{z})$, vamos minimizar:

$$-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

Agora que temos uma função de custo a minimizar cujas variáveis podem assumir os valores $-1$ e $1$, podemos fazer a seguinte analogia com o operador de Pauli $Z$:

$$z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}$$

Em outras palavras, a variável $z_i$ será equivalente a uma porta $Z$ atuando no qubit $i$. Além disso:

$$Z_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i$$

Assim, o observável que vamos considerar é:

$$\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j$$

ao qual teremos de adicionar o termo independente depois:

$$\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

O operador é uma combinação linear de termos com operadores Z nos nós conectados por uma aresta (lembrando que o qubit 0 está mais à direita): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Uma vez construído o operador, o ansatz para o algoritmo QAOA pode ser facilmente construído usando o circuito QAOAAnsatz da biblioteca de circuitos do Qiskit.

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp.from_list(
    [("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)

ansatz = QAOAAnsatz(hamiltonian, reps=2)
# Draw
ansatz.decompose(reps=3).draw("mpl")

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")

Offset: -2.5

Com o Runtime Estimator recebendo diretamente um Hamiltoniano e um ansatz parametrizado e retornando a energia necessária, a função de custo para uma instância QAOA é bastante simples:

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

import numpy as np

x0 = 2 * np.pi * np.random.rand(ansatz.num_parameters)

estimator = StatevectorEstimator()
cost = cost_func_vqe(x0, ansatz, hamiltonian, estimator)
print(cost)

1.473098768180865

# Estimated usage: < 1 min, benchmarked at 6 seconds on ibm_osaka, 5-23-24
# Load some necessary packages:

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator

# Select the least busy backend:

backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)

# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_hamiltonian = hamiltonian.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(x0, isa_ansatz, isa_hamiltonian, estimator)

# Close session after done
session.close()
print(cost)

1.1120776913677988

Vamos revisitar este exemplo em Aplicações para explorar como usar um otimizador para percorrer o espaço de busca. De maneira geral, isso inclui:

• Usar um otimizador para encontrar os parâmetros ótimos
• Vincular os parâmetros ótimos ao ansatz para encontrar os autovalores
• Traduzir os autovalores para a definição do nosso problema

Estratégia de medição: velocidade versus precisão

Como mencionado, estamos usando um computador quântico ruidoso como um oráculo de caixa-preta, onde o ruído pode tornar os valores obtidos não determinísticos, gerando flutuações aleatórias que, por sua vez, prejudicarão — ou até impedirão completamente — a convergência de certos otimizadores para uma solução proposta. Este é um problema geral que precisamos enfrentar à medida que exploramos incrementalmente a utilidade quântica e avançamos em direção à vantagem quântica:

Podemos usar as opções de supressão e mitigação de erros das Primitivas do Qiskit Runtime para lidar com o ruído e maximizar a utilidade dos computadores quânticos atuais.

Supressão de Erros

A supressão de erros refere-se a técnicas usadas para otimizar e transformar um circuito durante a compilação, a fim de minimizar erros. Essa é uma técnica básica de tratamento de erros que geralmente resulta em algum overhead de pré-processamento clássico no tempo total de execução. Esse overhead inclui a transpilação de circuitos para execução em hardware quântico por meio de:

• Expressão do circuito usando as portas nativas disponíveis no sistema quântico
• Mapeamento dos qubits virtuais para qubits físicos
• Adição de SWAPs com base nos requisitos de conectividade
• Otimização de portas de 1Q e 2Q
• Adição de desacoplamento dinâmico a qubits ociosos para prevenir os efeitos da decoerência.

As Primitivas permitem o uso de técnicas de supressão de erros por meio da configuração da opção optimization_level e da seleção de opções avançadas de transpilação. Em um curso posterior, aprofundaremos diferentes métodos de construção de circuitos para melhorar os resultados, mas para a maioria dos casos, recomendamos definir optimization_level=3.

Vamos visualizar o valor de aumentar a otimização no processo de transpilação observando um circuito de exemplo com um comportamento ideal simples.

from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")

qc = QuantumCircuit(2)
qc.x(1)
qc.h(0)
qc.cp(theta, 0, 1)
qc.h(0)
observables = SparsePauliOp.from_list([("ZZ", 1)])

qc.draw("mpl")

O circuito acima pode produzir valores esperados senoidais do observável fornecido, desde que inseramos fases abrangendo um intervalo adequado, como $[0,2\pi]$.

## Setup phases
import numpy as np

phases = np.linspace(0, 2 * np.pi, 50)

# phases need to be expressed as a list of lists in order to work
individual_phases = [[phase] for phase in phases]

Podemos usar um simulador para demonstrar a utilidade de uma transpilação otimizada. Voltaremos abaixo ao uso de hardware real para demonstrar a utilidade da mitigação de erros. Usaremos o QiskitRuntimeService para obter um backend real (neste caso, ibm_brisbane) e o AerSimulator para simular esse backend, incluindo seu comportamento de ruído.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

# get a real backend from the runtime service
service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# generate a simulator that mimics the real quantum system with the latest calibration results
backend_sim = AerSimulator.from_backend(backend)

Agora podemos usar um pass manager para transpilar o circuito para a "arquitetura de conjunto de instruções" (ISA) do backend. Este é um novo requisito no Qiskit Runtime: todos os circuitos enviados a um backend devem estar em conformidade com as restrições do alvo do backend, ou seja, devem ser escritos em termos da ISA do backend — isto é, o conjunto de instruções que o dispositivo consegue entender e executar. Essas restrições de alvo são definidas por fatores como as portas de base nativas do dispositivo, sua conectividade de qubits e — quando relevante — suas especificações de pulso e de temporização de outras instruções.

Observe que, no caso presente, faremos isso duas vezes: uma com optimization_level = 0 e outra com ele definido como 3. Em cada caso, usaremos a Primitiva Estimator para estimar os valores esperados do observável em diferentes valores de fase.

# Import estimator and specify that we are using the simulated backend:

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend_sim)

circuit = qc

# Use a pass manager to transpile the circuit and observable for the backend being simulated.
# Start with no optimization:

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=0)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

noisy_exp_values = []
pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
noisy_exp_values = cost[0]

# Repeat above steps, but now with optimization = 3:

exp_values_with_opt_es = []
pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=3)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
exp_values_with_opt_es = cost[0]

Por fim, podemos plotar os resultados e observar que a precisão do cálculo foi razoavelmente boa mesmo sem otimização, mas melhorou definitivamente ao aumentar a otimização para o nível 3. Observe que, em circuitos mais profundos e complexos, a diferença entre os níveis de otimização 0 e 3 tende a ser mais significativa. Este é um circuito muito simples usado como modelo de brinquedo.

import matplotlib.pyplot as plt

plt.plot(phases, noisy_exp_values, "o", label="opt=0")
plt.plot(phases, exp_values_with_opt_es, "o", label="opt=3")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Mitigação de Erros

A mitigação de erros refere-se a técnicas que permitem aos usuários reduzir erros em circuitos modelando o ruído do dispositivo no momento da execução. Tipicamente, isso resulta em overhead de pré-processamento quântico relacionado ao treinamento do modelo e overhead de pós-processamento clássico para mitigar erros nos resultados brutos usando o modelo gerado.

A opção resilience_level das Primitivas do Qiskit Runtime especifica a quantidade de resiliência a ser construída contra erros. Níveis mais altos geram resultados mais precisos à custa de tempos de processamento mais longos, devido ao overhead de amostragem quântica. Os níveis de resiliência podem ser usados para configurar o equilíbrio entre custo e precisão ao aplicar mitigação de erros na sua consulta de Primitiva.

Ao implementar qualquer técnica de mitigação de erros, esperamos que o viés em nossos resultados seja reduzido em relação ao viés anterior não mitigado. Em alguns casos, o viés pode até desaparecer. No entanto, isso tem um custo. À medida que reduzimos o viés em nossas estimativas, a variabilidade estatística aumenta (ou seja, a variância), o que podemos compensar aumentando ainda mais o número de shots por circuito no processo de amostragem. Isso introduz um overhead além do necessário para reduzir o viés, portanto não é feito por padrão. Podemos facilmente optar por esse comportamento ajustando o número de shots por circuito em options.executions.shots, conforme mostrado no exemplo abaixo.

Neste curso, exploraremos esses modelos de mitigação de erros em alto nível para ilustrar a mitigação que as Primitivas do Qiskit Runtime podem realizar, sem exigir detalhes completos de implementação.

Extinção de erros de leitura com twirling (T-REx)

A extinção de erros de leitura com twirling (T-REx) usa uma técnica conhecida como Pauli twirling para reduzir o ruído introduzido durante o processo de medição quântica. Essa técnica não assume nenhuma forma específica de ruído, o que a torna bastante geral e eficaz.

Fluxo de trabalho geral:

1. Adquirir dados para o estado zero com inversões de bits aleatórias (Pauli X antes da medição)
2. Adquirir dados para o estado desejado (ruidoso) com inversões de bits aleatórias (Pauli X antes da medição)
3. Calcular a função especial para cada conjunto de dados e dividir.

 

Podemos configurar isso com options.resilience_level = 1, demonstrado no exemplo abaixo.

Extrapolação para zero ruído

A extrapolação para zero ruído (ZNE) funciona amplificando primeiro o ruído no circuito que está preparando o estado quântico desejado, obtendo medições para vários níveis diferentes de ruído e usando essas medições para inferir o resultado sem ruído.

Fluxo de trabalho geral:

1. Amplificar o ruído do circuito para vários fatores de ruído
2. Executar cada circuito com ruído amplificado
3. Extrapolar de volta para o limite de zero ruído

 

Podemos configurar isso com options.resilience_level = 2. Podemos otimizar ainda mais explorando uma variedade de noise_factors, noise_amplifiers e extrapolators, mas isso está fora do escopo deste curso. Encorajamos você a experimentar essas opções conforme descritas aqui.

Cada método vem com seu próprio overhead associado: um compromisso entre o número de computações quânticas necessárias (tempo) e a precisão dos nossos resultados:

$$\begin{array}{c|c|c|c} \text{Methods} & R=1 \text{, T-REx} & R=2 \text{, ZNE} \\[1mm] \hline \text{Assumptions} & \text{None} & \text{Ability to scale noise} \\[1mm] \text{Qubit overhead} & 1 & 1 \\[1mm] \text{Sampling overhead} & 2 & N_{\text{noise-factors}} \\[1mm] \text{Bias} & 0 & \mathcal{O}(\lambda^{N_{\text{noise-factors}}}) \\[1mm] \end{array}$$

Usando as opções de mitigação e supressão do Qiskit Runtime

Veja como calcular um valor esperado usando mitigação e supressão de erros no Qiskit Runtime. Podemos usar exatamente o mesmo circuito e observável de antes, mas desta vez mantendo o nível de otimização fixo em 2 e ajustando a resiliência ou as técnicas de mitigação de erros utilizadas. Este processo de mitigação de erros ocorre várias vezes ao longo de um loop de otimização.

Realizamos esta parte em hardware real, pois a mitigação de erros não está disponível em simuladores.

# Estimated usage: 8 minutes, benchmarked on an Eagle processor, 5-23-24

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import (
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)

# We select the least busy backend

# Select the least busy backend
# backend = service.least_busy(
#    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
# )

# Or use a specific backend
backend = service.backend("ibm_brisbane")

# Initialize some variables to save the results from different runs:

exp_values_with_em0_es = []
exp_values_with_em1_es = []
exp_values_with_em2_es = []

# Use a pass manager to optimize the circuit and observables for the backend chosen:

pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

# Open a session and run with no error mitigation:

estimator_options = EstimatorOptions(resilience_level=0, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em0_es = cost[0]

# Open a session and run with resilience = 1:

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em1_es = cost[0]

# Open a session and run with resilience = 2:

estimator_options = EstimatorOptions(resilience_level=2, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em2_es = cost[0]

Como antes, podemos plotar os valores esperados resultantes em função do ângulo de fase para os três níveis de mitigação de erros usados. Com algum esforço, é possível perceber que a mitigação de erros melhora ligeiramente os resultados. Novamente, esse efeito é muito mais pronunciado em circuitos mais profundos e complexos.

import matplotlib.pyplot as plt

plt.plot(phases, exp_values_with_em0_es, "o", label="unmitigated")
plt.plot(phases, exp_values_with_em1_es, "o", label="resil = 1")
plt.plot(phases, exp_values_with_em2_es, "o", label="resil = 2")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Resumo

Nesta lição, você aprendeu a criar uma função de custo:

• Criar uma função de custo
• Como usar as Primitivas do Qiskit Runtime para mitigar e suprimir ruído
• Como definir uma estratégia de medição para otimizar velocidade versus precisão

Aqui está nossa carga de trabalho variacional em alto nível:

Nossa função de custo é executada a cada iteração do loop de otimização. A próxima lição explorará como o otimizador clássico usa a avaliação da nossa função de custo para selecionar novos parâmetros.

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

1.1.0
0.23.0