QAOA: полный авторский вариант 2021 года

Это полный авторский вариант упражнения из официальной тетрадки Qiskit Textbook по состоянию на 1 июля 2021 года, commit e03fac338f5a365daee994e1ec0d63ea353a376d. Он использует пятивершинный граф; позднейшая тетрадка с четырьмя вершинами в это условие не подмешивается.

Дан граф-бабочка V={0,1,2,3,4} с единичными рёбрами (0,1),(0,2),(1,2),(2,3),(3,4),(2,4). Требуется максимизировать число рёбер разреза. Обозначим

$$C(x)=\sum_{(j,k)\in E}[x_j\ne x_k],$$
$$H=\sum_{(j,k)\in E}\frac{I-Z_jZ_k}{2},$$
$$B=\sum_{j=0}^{4}X_j,\qquad |s\rangle=|+\rangle^{\otimes5}.$$

В разделе 5.3 автор задаёт $|\psi_1(\gamma,\beta)\rangle=e^{-i\beta B}e^{-i\gamma H}|s\rangle$, а в разделе 5.2 определяет $F_1=\langle\psi_1|H|\psi_1\rangle=2f_A+4f_B$. Тип A — рёбра (0,1),(3,4), тип B — остальные четыре; $f_{(j,k)}=\langle(I-Z_jZ_k)/2\rangle$. Унитарные операторы в локальном выводе записаны как $U_{jk}=e^{i\gamma(I-Z_jZ_k)/2}$ и $X_j(\beta)=e^{i\beta X_j}$, а дальние коммутирующие множители сокращаются. Поэтому f_A вычисляется по трём кубитам треугольника, f_B — по пяти вершинам объединённой окрестности ребра. Полный исходный вывод сохранён в происхождении и доступен по постоянной ссылке на источник.

Напечатанные автором выражения ниже — именно исходные выражения для проверки, а не заранее подтверждённые тождества. Для переноса по строкам введены только сокращения u,v,w:

$$\begin{aligned}f_A^{\mathrm{src}}={}&\tfrac12\sin4\gamma\sin4\beta\\&+\tfrac12\sin^22\beta\sin^22\gamma,\end{aligned}$$
$$f_B^{\mathrm{src}}=\tfrac12(1-u-v),$$
$$u=\sin^22\beta\sin^22\gamma\cos^24\gamma,$$
$$v=\tfrac14\sin4\beta\sin4\gamma(1+\cos^24\gamma),$$
$$F_1^{\mathrm{src}}=3-w(1+\cos^24\gamma),$$
$$w=\sin^22\beta\sin^22\gamma-\tfrac12\sin4\beta\sin4\gamma.$$

Пятикубитная схема задаётся полностью приведённой программой: равномерная подготовка, для каждого ребра CP(−2γ),P(γ) на обоих концах, затем RX(2β) на всех кубитах и измерение. Исходный поиск по сетке с шагом 0.1 сообщает γ=1.900,β=0.200 и максимум напечатанной функции 3.431. В сохранённом локальном моделировании было 10000 запусков, среднее 3.31 и лучшая найденная строка 01100 со стоимостью 4. В сохранённом опыте ibmq_santiago —2048 запусков, среднее 3.03 и лучшая строка 10001 со стоимостью 4. Все исходные графики состояния источника сохранены в условии, потому что пункт 0 требует их сравнить. Это архивные данные автора, не результаты новых запусков. Для их сравнения живой Santiago не требуется.

Полный авторский блок Problems:

0. QAOA выдаёт битовую строку. Является ли она оптимальным решением для этого графа? Сравните экспериментальные результаты сверхпроводящего чипа с результатами локального моделирования QASM.
1. В разделе 5.2 аналитически вычислена F₁. Проверьте шаги и вычислите f_A(γ,β) и f_B(γ,β).
2. В реализации Qiskit дано точное выражение F₁.
- Напишите процедуру оценки среднего F₁(γ,β) по выборке измерений; подсказка — cost_function_C(x,G) из раздела 5.4 и обработка данных разделов 5a/5b.
- Численно оптимизируйте параметры выборочной оценки, например SPSA из примера VQE. Получаются ли те же γ*,β*?
3. Пробная схема раздела 5.3 имеет p=1 и согласована с оборудованием.
- Процедурой из пункта 2 оцените F_p при p=2,3. Что ожидается на реальном оборудовании?
- Обобщите семейство пробных состояний, например на форму Ref.4 с учётом оборудования. Ref.4: Kandala et al., Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature549 (2017),242.
4. Рассмотрите пример MAX3-SAT, как в разделе примеров, и измените cost_function_C(c,G) из раздела 5.4, использованную для F_p. Запустите QAOA для этой задачи с формой, учитывающей оборудование, и разберите результат. В исходном теоретическом примере приведено $c_1(x)=(1-x_1)(1-x_3)x_{132}$ для 133 переменных; этот пример сохранён как написан, а конкретная новая задача выбирается в рамках открытого пункта автора.

Все пять пунктов, включая пятый пункт 4, отсутствующий в поздней версии, сохранены вместе. Исходная программа содержит аппаратный вызов только как контекст архивного опыта; новое выполнение на Santiago в этих пунктах не требуется.

import numpy as np
import networkx as nx  # tool to handle general Graphs 
import matplotlib.pyplot as plt 
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

from qiskit import Aer, IBMQ
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, transpile, assemble
from qiskit.providers.ibmq import least_busy
from qiskit.tools.monitor import job_monitor
from qiskit.visualization import plot_histogram

# Generating the butterfly graph with 5 nodes 
n = 5
V = np.arange(0,n,1)
E =[(0,1,1.0),(0,2,1.0),(1,2,1.0),(3,2,1.0),(3,4,1.0),(4,2,1.0)] 

G = nx.Graph()
G.add_nodes_from(V)
G.add_weighted_edges_from(E)

# Generate plot of the Graph
colors       = ['r' for node in G.nodes()]
default_axes = plt.axes(frameon=True)
pos          = nx.spring_layout(G)

nx.draw_networkx(G, node_color=colors, node_size=600, alpha=1, ax=default_axes, pos=pos)

# Evaluate the function
step_size   = 0.1;

a_gamma         = np.arange(0, np.pi, step_size)
a_beta          = np.arange(0, np.pi, step_size)
a_gamma, a_beta = np.meshgrid(a_gamma,a_beta)

F1 = 3-(np.sin(2*a_beta)**2*np.sin(2*a_gamma)**2-0.5*np.sin(4*a_beta)*np.sin(4*a_gamma))*(1+np.cos(4*a_gamma)**2)

# Grid search for the minimizing variables
result = np.where(F1 == np.amax(F1))
a      = list(zip(result[0],result[1]))[0]

gamma  = a[0]*step_size;
beta   = a[1]*step_size;

# Plot the expetation value F1
fig = plt.figure()
ax  = fig.gca(projection='3d')

surf = ax.plot_surface(a_gamma, a_beta, F1, cmap=cm.coolwarm, linewidth=0, antialiased=True)

ax.set_zlim(1,4)
ax.zaxis.set_major_locator(LinearLocator(3))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

plt.show()

#The smallest parameters and the expectation can be extracted
print('\n --- OPTIMAL PARAMETERS --- \n')
print('The maximal expectation value is:  M1 = %.03f' % np.amax(F1))
print('This is attained for gamma = %.03f and beta = %.03f' % (gamma,beta))

# prepare the quantum and classical resisters
QAOA = QuantumCircuit(len(V), len(V))

# apply the layer of Hadamard gates to all qubits
QAOA.h(range(len(V)))
QAOA.barrier()

# apply the Ising type gates with angle gamma along the edges in E
for edge in E:
    k = edge[0]
    l = edge[1]
    QAOA.cp(-2*gamma, k, l)
    QAOA.p(gamma, k)
    QAOA.p(gamma, l)
    
# then apply the single qubit X rotations with angle beta to all qubits
QAOA.barrier()
QAOA.rx(2*beta, range(len(V)))

# Finally measure the result in the computational basis
QAOA.barrier()
QAOA.measure(range(len(V)),range(len(V)))

### draw the circuit for comparison
QAOA.draw()

# Compute the value of the cost function
def cost_function_C(x,G):
    
    x = x[::-1]
    E = G.edges()
    if( len(x) != len(G.nodes())):
        return np.nan
        
    C = 0;
    for index in E:
        e1 = index[0]
        e2 = index[1]
        
        w      = G[e1][e2]['weight']
        C = C + w*x[e1]*(1-x[e2]) + w*x[e2]*(1-x[e1])
        
    return C

# run on local simulator
backend = Aer.get_backend("aer_simulator")
shots = 10000

TQAOA = transpile(QAOA, backend)
qobj = assemble(TQAOA, shots=shots)
QAOA_results = backend.run(qobj).result()

plot_histogram(QAOA_results.get_counts(),figsize = (8,6),bar_labels = False)

# Evaluate the data from the simulator
counts = QAOA_results.get_counts()

avr_C       = 0
max_C       = [0,0]
hist        = {}

for k in range(len(G.edges())+1):
    hist[str(k)] = hist.get(str(k),0)

for sample in list(counts.keys()):

    # use sampled bit string x to compute C(x)
    x         = [int(num) for num in list(sample)]
    tmp_eng   = cost_function_C(x,G)
    
    # compute the expectation value and energy distribution
    avr_C     = avr_C    + counts[sample]*tmp_eng
    hist[str(round(tmp_eng))] = hist.get(str(round(tmp_eng)),0) + counts[sample]
    
    # save best bit string
    if( max_C[1] < tmp_eng):
        max_C[0] = sample
        max_C[1] = tmp_eng
                
M1_sampled   = avr_C/shots

print('\n --- SIMULATION RESULTS ---\n')
print('The sampled mean value is M1_sampled = %.02f while the true value is M1 = %.02f \n' % (M1_sampled,np.amax(F1)))
print('The approximate solution is x* = %s with C(x*) = %d \n' % (max_C[0],max_C[1]))
print('The cost function is distributed as: \n')
plot_histogram(hist,figsize = (8,6),bar_labels = False)

provider = IBMQ.load_account()
backend = provider.get_backend('ibmq_santiago')
shots = 2048

TQAOA = transpile(QAOA, backend)
job_exp = backend.run(TQAOA, shots=shots)
job_monitor(job_exp)

exp_results = job_exp.result()
plot_histogram(exp_results.get_counts(),figsize = (10,8),bar_labels = False)

# Evaluate the data from the experiment
counts = exp_results.get_counts()

avr_C       = 0
max_C       = [0,0]
hist        = {}

for k in range(len(G.edges())+1):
    hist[str(k)] = hist.get(str(k),0)

for sample in list(counts.keys()):

    # use sampled bit string x to compute C(x)
    x         = [int(num) for num in list(sample)]
    tmp_eng   = cost_function_C(x,G)
    
    # compute the expectation value and energy distribution
    avr_C     = avr_C    + counts[sample]*tmp_eng
    hist[str(round(tmp_eng))] = hist.get(str(round(tmp_eng)),0) + counts[sample]
    
    # save best bit string
    if( max_C[1] < tmp_eng):
        max_C[0] = sample
        max_C[1] = tmp_eng
                
M1_sampled   = avr_C/shots

print('\n --- EXPERIMENTAL RESULTS ---\n')
print('The sampled mean value is M1_sampled = %.02f while the true value is M1 = %.02f \n' % (M1_sampled,np.amax(F1)))
print('The approximate solution is x* = %s with C(x*) = %d \n' % (max_C[0],max_C[1]))
print('The cost function is distributed as: \n')
plot_histogram(hist,figsize = (8,6),bar_labels = False)