Приклади та застосування

У цьому уроці ми розглянемо кілька прикладів варіаційних алгоритмів і способи їх застосування:

• Як написати власний варіаційний алгоритм
• Як застосувати варіаційний алгоритм для пошуку мінімальних власних значень
• Як використовувати варіаційні алгоритми для розв'язання прикладних задач

Зверни увагу: фреймворк Qiskit patterns можна застосувати до всіх задач, які ми тут розглядаємо. Втім, щоб уникнути повторень, ми явно виділятимемо кроки цього фреймворку лише в одному прикладі, що виконується на реальному обладнанні.

Формулювання задач

Уяви, що ми хочемо скористатися варіаційним алгоритмом для знаходження власного значення такого спостережуваного:

$$\hat{O}_1 = 2 II - 2 XX + 3 YY - 3 ZZ,$$

Це спостережуване має такі власні значення:

$$\left\{ \begin{array}{c} \lambda_0 = -6 \\ \lambda_1 = 4 \\ \lambda_2 = 4 \\ \lambda_3 = 6 \end{array} \right\}$$

І власні стани:

$$\left\{ \begin{array}{c} |\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)\\ |\phi_1\rangle = \frac{1}{\sqrt{2}}(|00\rangle - |11\rangle)\\ |\phi_2\rangle = \frac{1}{\sqrt{2}}(|01\rangle - |10\rangle)\\ |\phi_3\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |10\rangle) \end{array} \right\}$$

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

from qiskit.quantum_info import SparsePauliOp

observable_1 = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

Власний VQE

Спочатку розглянемо, як вручну побудувати екземпляр VQE для знаходження найменшого власного значення $\hat{O}_1$. Це охопить різноманітні техніки, які ми вивчали впродовж цього курсу.

def cost_func_vqe(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

    return cost

from qiskit.circuit.library.n_local import n_local
from qiskit import QuantumCircuit

import numpy as np

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

variational_form = n_local(
    num_qubits=2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)

raw_ansatz = reference_circuit.compose(variational_form)
raw_ansatz.decompose().draw("mpl")

Розпочнемо налагодження на локальних симуляторах.

from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler

estimator = Estimator()
sampler = Sampler()

Тепер задамо початковий набір параметрів:

x0 = np.ones(raw_ansatz.num_parameters)
print(x0)

[1. 1. 1. 1. 1. 1. 1. 1.]

Мінімізуємо цю функцію вартості для обчислення оптимальних параметрів

# SciPy minimizer routine
from scipy.optimize import minimize
import time

start_time = time.time()

result = minimize(
    cost_func_vqe,
    x0,
    args=(raw_ansatz, observable_1, estimator),
    method="COBYLA",
    options={"maxiter": 1000, "disp": True},
)

end_time = time.time()
execution_time = end_time - start_time

Return from COBYLA because the trust region radius reaches its lower bound.
Number of function values = 103   Least value of F = -5.999999998357189
The corresponding X is:
[2.27483579e+00 8.37593091e-01 1.57080508e+00 5.82932911e-06
 2.49973063e+00 6.41884255e-01 6.33686904e-01 6.33688223e-01]

result

message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -5.999999998357189
       x: [ 2.275e+00  8.376e-01  1.571e+00  5.829e-06  2.500e+00
            6.419e-01  6.337e-01  6.337e-01]
    nfev: 103
   maxcv: 0.0

Оскільки ця навчальна задача використовує лише два Qubit-и, ми можемо перевірити результат за допомогою вбудованого розв'язувача лінійної алгебри NumPy.

from numpy.linalg import eigvalsh

solution_eigenvalue = min(eigvalsh(observable_1.to_matrix()))

print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")

print(
    f"Percent error: {100*abs((result.fun - solution_eigenvalue)/solution_eigenvalue):.2e}"
)

Number of iterations: 103
Time (s): 0.4394676685333252
Percent error: 2.74e-08

Як бачиш, результат надзвичайно близький до ідеального.

Експерименти для підвищення швидкості та точності

Додаємо референсний стан

У попередньому прикладі ми не використовували жодного референсного оператора $U_R$. Тепер подумаймо, як можна отримати ідеальний власний стан $\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$. Розглянемо наступну схему.

from qiskit import QuantumCircuit

ideal_qc = QuantumCircuit(2)
ideal_qc.h(0)
ideal_qc.cx(0, 1)

ideal_qc.draw("mpl")

Швидко перевіримо, що ця схема дає нам потрібний стан.

from qiskit.quantum_info import Statevector

Statevector(ideal_qc)

Statevector([0.70710678+0.j, 0.        +0.j, 0.        +0.j,
             0.70710678+0.j],
            dims=(2, 2))

Тепер, коли ми бачили, як виглядає схема, що підготовлює розв'язок, логічно використати Gate Адамара як референсну схему, тоді повний ансац набуде вигляду:

reference = QuantumCircuit(2)
reference.h(0)
reference.cx(0, 1)
# Include barrier to separate reference from variational form
reference.barrier()

ref_ansatz = variational_form.decompose().compose(reference, front=True)

ref_ansatz.draw("mpl")

Для цієї нової схеми ідеальний розв'язок можна отримати при всіх параметрах, рівних $0$, що підтверджує розумність вибору референсної схеми.

Тепер порівняємо кількість обчислень функції вартості, ітерацій оптимізатора та витрачений час із результатами попередньої спроби.

import time

start_time = time.time()

ref_result = minimize(
    cost_func_vqe, x0, args=(ref_ansatz, observable_1, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

Використовуємо оптимальні параметри для обчислення мінімального власного значення:

experimental_min_eigenvalue_ref = cost_func_vqe(
    ref_result.x, ref_ansatz, observable_1, estimator
)
print(experimental_min_eigenvalue_ref)

-5.999999996759607

print("ADDED REFERENCE STATE:")
print(f"""Number of iterations: {ref_result.nfev}""")
print(f"""Time (s): {execution_time}""")
print(
    f"Percent error: {100*abs((experimental_min_eigenvalue_ref - solution_eigenvalue)/solution_eigenvalue):.2e}"
)

ADDED REFERENCE STATE:
Number of iterations: 127
Time (s): 0.5620882511138916
Percent error: 5.40e-08

Залежно від конкретної системи, це може або не може покращити швидкість чи точність у цьому дуже маленькому прикладі. Головне те, що використання фізично вмотивованих референсних станів стає дедалі важливішим для підвищення швидкості та точності в міру масштабування задач.

Змінюємо початкову точку

Тепер, коли ми побачили ефект додавання референсного стану, розглянемо, що відбувається при виборі різних початкових точок $\vec{\theta_0}$. Зокрема, ми використаємо $\vec{\theta_0}=(0,0,0,0,6,0,0,0)$ та $\vec{\theta_0}=(6,6,6,6,6,6,6,6,6)$.

Пам'ятай, що, як обговорювалось при введенні референсного стану, ідеальний розв'язок знаходиться при всіх параметрах, рівних $0$, тому перша початкова точка повинна давати менше обчислень.

import time

start_time = time.time()

x0 = [0, 0, 0, 0, 6, 0, 0, 0]

x0_1_result = minimize(
    cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

print("INITIAL POINT 1:")
print(f"""Number of iterations: {x0_1_result.nfev}""")
print(f"""Time (s): {execution_time}""")

INITIAL POINT 1:
Number of iterations: 108
Time (s): 0.4492197036743164

Змінюємо початкову точку на $\vec{\theta_0}=(6,6,6,6,6,6,6,6,6)$:

import time

start_time = time.time()

x0 = 6 * np.ones(raw_ansatz.num_parameters)

x0_2_result = minimize(
    cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

print("INITIAL POINT 2:")
print(f"""Number of iterations: {x0_2_result.nfev}""")
print(f"""Time (s): {execution_time}""")

INITIAL POINT 2:
Number of iterations: 107
Time (s): 0.40889453887939453

Експериментуючи з різними початковими точками, ти можеш досягти збіжності швидше та з меншою кількістю обчислень функції.

Експерименти з різними оптимізаторами

Оптимізатор можна змінити через аргумент method функції minimize з SciPy; більше варіантів можна знайти тут. Спочатку ми використовували обмежений мінімізатор (COBYLA). У цьому прикладі спробуємо необмежений мінімізатор (BFGS).

import time

start_time = time.time()

result = minimize(
    cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="BFGS"
)

end_time = time.time()
execution_time = end_time - start_time

print("CHANGED TO BFGS OPTIMIZER:")
print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")

CHANGED TO BFGS OPTIMIZER:
Number of iterations: 117
Time (s): 0.31656408309936523

Приклад VQD

Тут ми явно реалізуємо фреймворк Qiskit patterns.

Крок 1: Відображення класичних вхідних даних на квантову задачу

Тепер замість пошуку лише найменшого власного значення наших спостережуваних ми шукатимемо всі $4$ (де $k=4$).

Пригадай, що функції вартості VQD мають вигляд:

$$C_{l}(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle + \sum_{j=0}^{l-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2 \quad \forall l\in\{1,\cdots,k\}=\{1,\cdots,4\}$$

Це особливо важливо, оскільки вектор $\vec{\beta}=(\beta_0,\cdots,\beta_{k-1})$ (у цьому випадку $(\beta_0, \beta_1, \beta_2, \beta_3)$) має передаватись як аргумент при визначенні об'єкта VQD.

Крім того, в реалізації VQD у Qiskit замість розгляду ефективних спостережуваних, описаних у попередньому зошиті, точності $|\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$ обчислюються безпосередньо за допомогою алгоритму ComputeUncompute, який використовує примітив Sampler для вибірки ймовірності отримання $|0\rangle$ для схеми $U_A^\dagger(\vec{\theta})U_A(\vec{\theta^j})$. Це працює саме тому, що ця ймовірність дорівнює

$$\begin{aligned} p_0 & = |\langle 0|U_A^\dagger(\vec{\theta})U_A(\vec{\theta^j})|0\rangle|^2 \\[1mm] & = |\big(\langle 0|U_A^\dagger(\vec{\theta})\big)\big(U_A(\vec{\theta^j})|0\rangle\big)|^2 \\[1mm] & = |\langle \psi(\vec{\theta}) |\psi(\vec{\theta^j}) \rangle|^2 \\[1mm] \end{aligned}$$

ansatz = n_local(
    num_qubits=2,
    rotation_blocks=["ry", "rz"],
    entanglement_blocks="cz",
    #    entanglement="linear",
    reps=1,
)

ansatz.decompose().draw("mpl")

Почнемо з розгляду такого спостережуваного:

$$\hat{O}_2 := 2 II - 3 XX + 2 YY - 4 ZZ$$

Це спостережуване має такі власні значення:

$$\left\{ \begin{array}{c} \lambda_0 = -7 \\ \lambda_1 = 3\\ \lambda_2 = 5 \\ \lambda_3 = 7 \end{array} \right\}$$

І власні стани:

$$\left\{ \begin{array}{c} |\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)\\ |\phi_1\rangle = \frac{1}{\sqrt{2}}(|00\rangle - |11\rangle)\\ |\phi_2\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |10\rangle)\\ |\phi_3\rangle = \frac{1}{\sqrt{2}}(|01\rangle - |10\rangle) \end{array} \right\}$$

from qiskit.quantum_info import SparsePauliOp

observable_2 = SparsePauliOp.from_list([("II", 2), ("XX", -3), ("YY", 2), ("ZZ", -4)])

Тепер скористаємося такою функцією для обчислення штрафу за перекриття. Зверни увагу, що це все ще частина відображення задачі на квантові схеми. Однак, як обговорювалось у попередньому уроці, ця функція обчислює перекриття між поточною варіаційною схемою та оптимізованою схемою з попереднього стану з нижчою енергією/вартістю. Нова схема, що генерується, також має бути транспільована для запуску на реальному обладнанні. Ми вже бачили цю функцію раніше, коли використовували симулятор. Тут нам треба вже враховувати транспіляцію та пов'язану оптимізацію для використання реального Backend-у, звідси рядки навколо if realbackend == 1. Це трохи змішує крок 2, але ми явно виділимо крок 2 пізніше.

import numpy as np
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

def calculate_overlaps(
    ansatz, prev_circuits, parameters, sampler, realbackend, backend
):
    def create_fidelity_circuit(circuit_1, circuit_2):
        if len(circuit_1.clbits) > 0:
            circuit_1.remove_final_measurements()
        if len(circuit_2.clbits) > 0:
            circuit_2.remove_final_measurements()

        circuit = circuit_1.compose(circuit_2.inverse())
        circuit.measure_all()
        return circuit

    overlaps = []

    for prev_circuit in prev_circuits:
        fidelity_circuit = create_fidelity_circuit(ansatz, prev_circuit)
        if realbackend == 1:
            pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
            fidelity_circuit = pm.run(fidelity_circuit)
        sampler_job = sampler.run([(fidelity_circuit, parameters)])
        meas_data = sampler_job.result()[0].data.meas

        counts_0 = meas_data.get_int_counts().get(0, 0)
        shots = meas_data.num_shots
        overlap = counts_0 / shots
        overlaps.append(overlap)

    return np.array(overlaps)

Тепер додамо функцію вартості VQD. Зверни увагу, що порівняно з попереднім уроком тут з'явилися два додаткові аргументи (realbackend та backend) для підтримки транспіляції при використанні реальних Backend-ів.

def cost_func_vqd(
    parameters,
    ansatz,
    prev_states,
    step,
    betas,
    estimator,
    sampler,
    hamiltonian,
    realbackend,
    backend,
):
    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])

    total_cost = 0

    if step > 1:
        overlaps = calculate_overlaps(
            ansatz, prev_states, parameters, sampler, realbackend, backend
        )
        total_cost = np.sum(
            [np.real(betas[state] * overlap) for state, overlap in enumerate(overlaps)]
        )

    estimator_result = estimator_job.result()[0]

    value = estimator_result.data.evs[0] + total_cost

    return value

Знову ж таки, для налагодження використаємо симулятори, а потім перейдемо до реального обладнання.

from qiskit.primitives import StatevectorSampler
from qiskit.primitives import StatevectorEstimator

sampler = StatevectorSampler(default_shots=4092)
estimator = StatevectorEstimator()

Тут вводимо кількість станів, які ми хочемо обчислити, штрафи та початковий набір параметрів x0.

from qiskit.quantum_info import SparsePauliOp

k = 4
betas = [50, 60, 40]
x0 = np.ones(8)

Тепер протестуємо алгоритм за допомогою симуляторів:

from scipy.optimize import minimize

prev_states = []
prev_opt_parameters = []
eigenvalues = []

realbackend = 0

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(
            ansatz,
            prev_states,
            step,
            betas,
            estimator,
            sampler,
            observable_2,
            realbackend,
            None,
        ),
        method="COBYLA",
        options={"maxiter": 200, "tol": 0.000001},
    )
    print(result)

    prev_opt_parameters = result.x
    eigenvalues.append(result.fun)

message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -6.9999999999996
       x: [ 1.571e+00  1.571e+00  2.519e+00  2.100e+00  1.242e+00
            6.935e-01  2.298e+00  1.991e+00]
    nfev: 151
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 3.698974255258432
       x: [ 1.269e+00  1.109e+00  1.080e+00  1.200e+00  1.094e+00
            1.163e+00  9.752e-01  9.519e-01]
    nfev: 103
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 4.731320121938101
       x: [ 1.533e+00  2.451e+00  2.526e+00  2.406e+00  1.968e+00
            2.105e+00  8.537e-01  8.442e-01]
    nfev: 110
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 7.008239313655201
       x: [ 4.150e+00  2.120e+00  3.495e+00  7.262e-01  1.953e+00
           -1.982e-01  3.263e-01  2.563e+00]
    nfev: 126
   maxcv: 0.0

eigenvalues

[np.float64(-6.9999999999996),
 np.float64(3.698974255258432),
 np.float64(4.731320121938101),
 np.float64(7.008239313655201)]

Ці результати досить близькі до очікуваних, за винятком похибки апроксимації та глобальної фази. Ми могли б відрегулювати допуск на класичному оптимізаторі та/або штрафи за перекриття вектора стану, щоб отримати точніші значення.

solution_eigenvalues = [-7, 3, 5, 7]

for index, experimental_eigenvalue in enumerate(eigenvalues):
    solution_eigenvalue = solution_eigenvalues[index]

    print(
        f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
    )

Percent error: 5.71e-14
Percent error: 2.33e-01
Percent error: 5.37e-02
Percent error: 1.18e-03

Зміна бета-значень

Як згадувалось у попередньому уроці, значення $\vec{\beta}$ мають бути більшими за різницю між власними значеннями. Подивімося, що відбувається, коли ця умова не виконується для $\hat{O}_2$

$$\hat{O}_2 = 2 II - 3 XX + 2 YY - 4 ZZ$$

з власними значеннями

$$\left\{ \begin{array}{c} \lambda_0 = -7 \\ \lambda_1 = 3\\ \lambda_2 = 5 \\ \lambda_3 = 7 \end{array} \right\}$$

from qiskit.quantum_info import SparsePauliOp

k = 4
betas = np.ones(3)
x0 = np.zeros(8)

from scipy.optimize import minimize

prev_states = []
prev_opt_parameters = []
eigenvalues = []

realbackend = 0

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(
            ansatz,
            prev_states,
            step,
            betas,
            estimator,
            sampler,
            observable_2,
            realbackend,
            None,
        ),
        method="COBYLA",
        options={"tol": 0.01, "maxiter": 200},
    )
    print(result)

    prev_opt_parameters = result.x
    eigenvalues.append(result.fun)

message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -6.999916534745094
       x: [ 1.568e+00 -1.569e+00  1.385e-01  1.398e-01 -7.972e-01
            7.835e-01 -2.375e-01  4.539e-02]
    nfev: 125
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -1.515139929812874
       x: [-5.317e-04 -2.514e-03  1.016e+00  9.998e-01  3.890e-04
            1.772e-04  1.568e-04  8.497e-04]
    nfev: 35
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -0.509948114293115
       x: [-3.796e-03  8.853e-03  3.015e-04  9.997e-01  6.271e-04
           -2.554e-03  1.017e-04  2.766e-04]
    nfev: 37
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 0.4914672235935682
       x: [-7.178e-03 -8.652e-03  1.125e+00 -5.428e-02 -1.586e-03
            2.031e-03 -3.462e-03  5.734e-03]
    nfev: 35
   maxcv: 0.0

solution_eigenvalues = [-7, 3, 5, 7]

for index, experimental_eigenvalue in enumerate(eigenvalues):
    solution_eigenvalue = solution_eigenvalues[index]

    print(
        f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
    )

Percent error: 1.19e-05
Percent error: 1.51e+00
Percent error: 1.10e+00
Percent error: 9.30e-01

Цього разу оптимізатор повертає той самий стан $|\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$ як пропонований розв'язок для всіх власних станів — що явно невірно. Це трапляється через те, що значення бета виявились занадто малими, щоб штрафувати мінімальний власний стан у послідовних функціях вартості. Через це він не виключався з ефективного простору пошуку в наступних ітераціях алгоритму, і щоразу обирався як найкращий можливий розв'язок.

Рекомендуємо поекспериментувати зі значеннями $\vec{\beta}$ та переконатися, що вони більші за різницю між власними значеннями.

Крок 2: Оптимізація задачі для квантового виконання

Щоб запустити це на реальному залізі, потрібно оптимізувати квантові схеми під обраний квантовий комп'ютер. У нашому випадку ми просто скористаємося найменш завантаженим Backend.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
# Or use a specific backend
# backend = service.backend("ibm_brisbane")
print(backend)

<IBMBackend('ibm_brisbane')>

Ми виконаємо транспіляцію нашої схеми за допомогою готового менеджера проходів з рівнем оптимізації 3.

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = observable_2.apply_layout(layout=isa_ansatz.layout)

Крок 3: Виконання за допомогою примітивів Qiskit

Подбавши про те, щоб скинути значення бета до достатньо великих, тепер можемо запустити обчислення на реальному квантовому залізі.

# Estimated compute resource usage: 25 minutes. Benchmarked at 24 min, 30 sec on an Eagle r3 processor on 5-30-24

k = 2
betas = [30, 50, 80]
x0 = np.zeros(8)

real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []

realbackend = 1

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

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

    for step in range(1, k + 1):
        if step > 1:
            real_prev_states.append(isa_ansatz.assign_parameters(prev_opt_parameters))

        result = minimize(
            cost_func_vqd,
            x0,
            args=(
                isa_ansatz,
                real_prev_states,
                step,
                betas,
                estimator,
                sampler,
                isa_observable,
                realbackend,
                backend,
            ),
            method="COBYLA",
            options={"maxiter": 200},
        )
        print(result)

        real_prev_opt_parameters = result.x
        real_eigenvalues.append(result.fun)

session.close()
print(real_eigenvalues)

Крок 4: Постобробка, повернення результату в класичному форматі

Результат за структурою схожий на те, що обговорювалось у попередніх уроках і прикладах. Проте в наведених вище результатах є дещо проблематичне, з чого можна зробити застережний висновок щодо збуджених станів. Щоб обмежити обчислювальний час у цьому навчальному прикладі, ми встановили максимальну кількість ітерацій для класичного оптимізатора, яка, можливо, була занадто малою: 200 ітерацій. Попереднє обчислення на симуляторі не змогло зійтися за 200 ітерацій. Тут наше — зійшлося... але з якою точністю? Ми не задали допуск, за якого COBYLA вважав би себе «збіжним». Поглянувши на значення функції та порівнявши з попередніми запусками, бачимо, що COBYLA був далеко від збіжності до потрібної нам точності.

Є ще одна проблема: енергія першого збудженого стану виявляється нижчою за енергію основного стану! Спробуй пояснити, як таке може трапитися. Підказка: це пов'язано з точкою збіжності, яку ми щойно розглянули. Детальне пояснення цієї поведінки наведено нижче — після застосування VQD до молекули H2.

Квантова хімія: розрахунок основного та збудженого енергетичних станів

Наша мета — мінімізувати математичне сподівання спостережуваної величини, що відповідає енергії (гамільтоніан $\hat{\mathcal{H}}$):

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

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

H2_op = SparsePauliOp.from_list(
    [
        ("II", -1.052373245772859),
        ("IZ", 0.39793742484318045),
        ("ZI", -0.39793742484318045),
        ("ZZ", -0.01128010425623538),
        ("XX", 0.18093119978423156),
    ]
)

chem_ansatz = efficient_su2(H2_op.num_qubits)

chem_ansatz.decompose().draw("mpl")

from qiskit import QuantumCircuit

def cost_func_vqe(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 = np.ones(chem_ansatz.num_parameters)

Ми можемо мінімізувати цю функцію вартості для знаходження оптимальних параметрів, і спочатку перевіримо код на локальному симуляторі.

from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler

estimator = Estimator()
sampler = Sampler()

# SciPy minimizer routine
from scipy.optimize import minimize
import time

start_time = time.time()

result = minimize(
    cost_func_vqe, x0, args=(chem_ansatz, H2_op, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

result

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.857275029048451
       x: [ 7.326e-01  1.354e+00 ...  1.040e+00  1.508e+00]
    nfev: 242
   maxcv: 0.0

Мінімальне значення функції вартості (-1.857...) — це енергія основного стану молекули H2 в одиницях хартрі.

Збуджені стани

Ми також можемо скористатися VQD, щоб розв'язати задачу для $k=2$ загальних станів (основного стану та першого збудженого стану).

from qiskit.quantum_info import SparsePauliOp
import numpy as np

k = 2
betas = [33, 33]
# x0 = np.zeros(ansatz.num_parameters)
x0 = [
    1.164e00,
    -2.438e-01,
    9.358e-04,
    6.745e-02,
    1.990e00,
    9.810e-02,
    6.154e-01,
    5.454e-01,
]

Додамо обчислення перекриття:

from scipy.optimize import minimize

prev_states = []
prev_opt_parameters = []
eigenvalues = []

realbackend = 0

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(
            ansatz,
            prev_states,
            step,
            betas,
            estimator,
            sampler,
            H2_op,
            realbackend,
            None,
        ),
        method="COBYLA",
        options={"tol": 0.001, "maxiter": 2000},
    )
    print(result)

    prev_opt_parameters = result.x
    eigenvalues.append(result.fun)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.8572671093941977
       x: [ 1.164e+00 -2.437e-01  2.118e-03  6.448e-02  1.990e+00
            9.870e-02  6.167e-01  5.476e-01]
    nfev: 58
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.0322873777662176
       x: [ 3.205e+00  1.502e+00  1.699e+00 -1.107e-02  3.086e+00
            1.530e+00  4.445e-02  7.013e-02]
    nfev: 99
   maxcv: 0.0

eigenvalues

[-1.8572671093941977, -1.0322873777662176]

Реальне залізо і фінальне застереження

Щоб запустити це на реальному залізі, потрібно оптимізувати квантові схеми під обраний квантовий комп'ютер. У нашому випадку ми просто скористаємося найменш завантаженим Backend.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)

Для транспіляції скористаємося готовим менеджером проходів і максимально оптимізуємо схему, використовуючи рівень оптимізації 3.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = H2_op.apply_layout(layout=isa_ansatz.layout)

Оскільки VQD є дуже ітеративним, ми виконаємо всі кроки всередині Runtime-сесії, щоб завдання стояли в черзі лише на початку, а не між кожним оновленням параметрів. В усьому іншому синтаксис функції вартості та Estimator залишається незмінним.

x0 = [
    1.306e00,
    -2.284e-01,
    6.913e-02,
    -2.530e-02,
    1.849e00,
    7.433e-02,
    6.366e-01,
    5.600e-01,
]

# Estimated hardware usage: 20 min benchmarked on an Eagle r3 processor on 5-30-24

real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []

realbackend = 1

estimator_options = EstimatorOptions(resilience_level=1, default_shots=4096)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session)
    sampler = Sampler(mode=session)

    for step in range(1, k + 1):
        if step > 1:
            real_prev_states.append(
                isa_ansatz.assign_parameters(real_prev_opt_parameters)
            )

        result = minimize(
            cost_func_vqd,
            x0,
            args=(
                isa_ansatz,
                real_prev_states,
                step,
                betas,
                estimator,
                sampler,
                isa_observable,
                realbackend,
                backend,
            ),
            method="COBYLA",
            options={"tol": 0.001, "maxiter": 300},
        )
        print(result)

        real_prev_opt_parameters = result.x
        real_eigenvalues.append(result.fun)

session.close()
print(real_eigenvalues)

Отримана енергія основного стану (-1.83 хартрі) не надто далека від правильного значення (-1.85 хартрі). Однак енергія збудженого стану досить сильно відхиляється. Це схоже на помилкову поведінку, яку ми бачили раніше в цьому уроці. Reported енергія збудженого стану майже збігається з енергією основного стану. У попередньому випадку ми навіть спостерігали енергію збудженого стану, яка була нижчою за reported енергію основного стану.

Варіаційний розрахунок не може дати енергію, нижчу за справжню енергію основного стану. У попередньому випадку отримана нами енергія основного стану була далека від справжньої. Оскільки ми не знайшли справжньої енергії основного стану, ніякого протиріччя немає. У цьому ж випадку енергія основного стану була досить близькою до правильного значення, і тим не менш енергія збудженого стану виявляється дивно близькою до цього самого значення.

Щоб краще зрозуміти, як це трапилося, пригадай: збуджений стан ми знаходимо, вимагаючи ортогональності варіаційного стану до основного стану (за допомогою схем перекриття та штрафних доданків). Якщо нам не вдається точно отримати енергію основного стану (або ми помиляємося на кілька відсотків), то ми також не отримуємо точного вектора основного стану! Тому, вимагаючи ортогональності збудженого стану до першого знайденого стану, ми накладали умову ортогональності не до справжнього основного стану, а до якогось його наближення (іноді дуже грубого). Відповідно, збуджений стан не був змушений бути ортогональним до справжнього основного стану, і наші оцінки енергій збуджених станів виявились дуже близькими до енергії основного стану.

Це завжди буде проблемою у VQD. Але в принципі це можна виправити, збільшивши максимальну кількість ітерацій для класичного оптимізатора, встановивши менший допуск для класичного оптимізатора, а також, можливо, спробувавши інший анзац, якщо нам систематично не вдається знайти справжній основний стан. Як ми бачили, може також знадобитися змінити штрафні коефіцієнти перекриття (бета). Але це, по суті, окрема проблема. Жоден штраф за перекриття не допоможе тобі уникнути справжнього основного стану, якщо ти не знайшов достатньо точної оцінки цього самого основного стану для схеми перекриття.

Оптимізація: Max-Cut

Задача максимального розрізу (Max-Cut) — це задача комбінаторної оптимізації, яка полягає у поділі вершин графа на дві непересічних множини таким чином, щоб кількість ребер між цими двома множинами була максимальною. Формально: дано неорієнтований граф $G=(V,E)$, де $V$ — множина вершин, а $E$ — множина ребер; задача Max-Cut вимагає розбити вершини на дві непересічних підмножини, $S$ і $T$, так щоб кількість ребер, один кінець яких належить $S$, а інший — $T$, була максимальною.

Max-Cut можна застосовувати для розв'язання різноманітних задач, таких як кластеризація, проектування мереж та фазові переходи. Почнімо зі створення графа задачі:

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

Цю задачу можна сформулювати як задачу бінарної оптимізації. Для кожного вузла $0 \leq i < n$, де $n$ — кількість вузлів графа (у нашому випадку $n=4$), розглянемо бінарну змінну $x_i$. Ця змінна матиме значення $1$, якщо вузол $i$ належить до групи, яку ми позначимо як $1$, і $0$ — якщо він у іншій групі, яку позначимо як $0$. Також позначимо $w_{ij}$ (елемент $(i,j)$ матриці суміжності $w$) — вагу ребра між вузлами $i$ та $j$. Оскільки граф неорієнтований, $w_{ij}=w_{ji}$. Тоді нашу задачу можна сформулювати як максимізацію такої функції вартості:

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

Щоб розв'язати цю задачу за допомогою квантового комп'ютера, виразимо функцію вартості як математичне сподівання деякого оператора. Однак оператори, які Qiskit підтримує нативно, складаються з операторів Паулі, що мають власні значення $1$ та $-1$, а не $0$ та $1$. Тому зробимо таку заміну змінних:

Де $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Для зручного доступу до ваг ребер можна використовувати матрицю суміжності $w$. Це допоможе отримати нашу функцію вартості:

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

Звідси випливає:

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

Отже, нова функція вартості, яку ми хочемо максимізувати, має вигляд:

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

Крім того, квантовий комп'ютер природно шукає мінімуми (зазвичай найнижчий рівень енергії), а не максимуми — тому замість максимізації $C(\vec{z})$ ми будемо мінімізувати:

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

Тепер, маючи функцію вартості для мінімізації, змінні якої можуть набувати значень $-1$ та $1$, проведемо таку аналогію з оператором Паулі $Z$:

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

Іншими словами, змінна $z_i$ буде еквівалентна Gate $Z$, що діє на Qubit $i$. Крім того:

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

Тоді оператор, який ми розглядатимемо, матиме вигляд:

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

до якого потім потрібно додати незалежний доданок:

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

Оператор є лінійною комбінацією членів з операторами Z на вузлах, з'єднаних ребром (нагадаємо, що нульовий Qubit знаходиться крайнім праворуч): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Після побудови оператора схему-анзац для алгоритму QAOA можна легко створити за допомогою Circuit QAOAAnsatz з бібліотеки схем Qiskit.

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

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

max_ansatz = QAOAAnsatz(max_hamiltonian, reps=2)
# Draw
max_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

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

from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler

estimator = Estimator()
sampler = Sampler()

Тепер задамо початковий набір випадкових параметрів:

import numpy as np

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

[6.0252949  0.58448176 2.15785731 1.13646074]

Для мінімізації функції вартості можна використати будь-який класичний оптимізатор. На реальній квантовій системі зазвичай краще працюють оптимізатори, розроблені для нерівних ландшафтів функцій вартості. Тут ми використовуємо процедуру COBYLA з SciPy через функцію minimize.

Оскільки ми ітеративно виконуємо багато викликів до Runtime, використовуємо сесію для виконання всіх викликів у межах одного блоку. Більш того, у QAOA розв'язок закодовано у вихідному розподілі схеми-анзацу, прив'язаного до оптимальних параметрів мінімізації. Тому нам знадобиться примітив Sampler, який ми ініціалізуємо з тією самою сесією. Запустимо процедуру мінімізації:

result = minimize(
    cost_func, x0, args=(max_ansatz, max_hamiltonian, estimator), method="COBYLA"
)

print(result)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -2.585287311689236
       x: [ 7.332e+00  3.904e-01  2.045e+00  1.028e+00]
    nfev: 80
   maxcv: 0.0

Вектор розв'язку з кутами параметрів (x), підставлений у схему-анзац, дає розбиття графа, яке ми шукали.

eigenvalue = cost_func(result.x, max_ansatz, max_hamiltonian, estimator)
print(f"""Eigenvalue: {eigenvalue}""")
print(f"""Max-Cut Objective: {eigenvalue + offset}""")

Eigenvalue: -2.585287311689236
Max-Cut Objective: -5.085287311689235

from qiskit.result import QuasiDistribution
from qiskit.primitives import StatevectorSampler

sampler = StatevectorSampler()

# Assign solution parameters to ansatz
qc = max_ansatz.assign_parameters(result.x)

# Add measurements to our circuit
qc.measure_all()

# Sample ansatz at optimal parameters
# samp_dist = sampler.run(qc).result().quasi_dists[0]

shots = 1024
job = sampler.run([qc], shots=shots)

qc.decompose().draw("mpl")

data_pub = job.result()[0].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()}
)
probabilities = quasi_dist

# Close the session since we are now done with it
# session.close()

from qiskit.visualization import plot_distribution

plot_distribution(counts)

binary_string = max(counts.items(), key=lambda kv: kv[1])[0]
x = np.asarray([int(y) for y in reversed(list(binary_string))])

colors = ["r" if x[i] == 0 else "c" for i in range(n)]
mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color=colors
)

Підсумок

У цьому уроці ти навчився:

• Як писати власний варіаційний алгоритм
• Як застосовувати варіаційний алгоритм для пошуку мінімальних власних значень
• Як використовувати варіаційні алгоритми для розв'язання прикладних задач