← Materiały bloku: Kwantowe wyżarzanie kombinatorycznych problemów optymalizacyjnych  ·  podgląd notatnika 07_Twierdzenie_adiabatyczne_i_wyzarzanie_kwantowe.ipynb  ·  pobierz plik .ipynb

Twierdzenie adiabatyczne

Dla potrzeb kursu, sformułujemy intuicyjnie twierdzenie adiabatyczne. Pełne i formalne zapisy mechaniki kwantowej pominiemy, zostawiając je do dogłębnej analizy zainteresowanym kursantom.

Twierdzenie adiabatyczne (w mechanice kwantowej) Rozważmy system kwantowy zmieniający się w czasie. Jeżeli na początku ewolucji znajduje się on w stanie podstawowym, to przy pewnych technicznych założeniach, o ile zmiany systemu następują wystarczająco wolno, uzyskany na końcu system kwantowy jest również w stanie podstawowym. W skrócie, powolna ewolucja systemu przekształca stan podstawowy na stan podstawowy.

Sformułowanie wyżarzania kwantowego

W wyżarzaniu kwantowym wykorzystuje się to twierdzenie, aby znaleźć stan podstawowy (minimum globalne) hamiltonianu Isinga, który koduje problem optymalizacyjny. Proces wygląda tak:

1. Startujemy od prostego hamiltonianu $H_{initial}$​, którego stan podstawowy jest znany (np. prosta superpozycja).

2. Stopniowo przekształcamy go w hamiltonian problemu $H_{problem}$​, którego minimum chcemy znaleźć.

3. Jeśli proces jest adiabatyczny (czyli wystarczająco powolny), to układ przejdzie z podstawowego stanu $H_{initial}$​ do podstawowego stanu $H_{problem}$​, dając rozwiązanie.

tutaj wyżarzanie kwantowe zostanie zaimplementowane następująco:

$$ \mathcal{H}(s) = A(s)H_{initial} + B(s)H_{Ising} =A(s) \sum_{i}^{N} \sigma_i^x + B(s) \left( \sum_{j < i} J_{ij} \sigma_i^z \sigma_j^z + \sum_{i=1}^{N} h_i \sigma_i^z \right) $$

gdzie $s = t/\tau$, $\tau$ to całkowity czas wyżarzania a funkcje $A(s)$ i $B(s)$ to dowolne ciągłe funkcje spełniające warunki $A(0) \gg B(0)$, oraz $B(1) \gg A(1)$ oraz $A(1) \approx 0$.

Ewolucja układu

Ewolucja układu jest opisywana przez równanie Schrödingera

$$ i \frac{\partial}{\partial t} \ket{\Psi (t)} = \mathcal{H}\ket{\Psi (t)} $$

gdzie $t$ jest czasem, $\ket{\Psi (t)}$ jest wektorem stanu a $\mathcal{H}$ jest Hamiltonianem opisującym układ.

Jako że mechanika kwantowa jest z natury probabilistyczna i na końcu ewolucji otrzymujemy pewien wektor prawdopodobieństw, interesuje nas wartość oczekiwana pomiaru energi $<\mathcal{H}(t)>$. Jest to średnia ważona wszystkich możliwych pomiarów klasycznego stanu (tj. wektoru spinów).

# Funkcje pomocniczne które będziemy używać do robienia symulacji wyżarzania kwantowego
from setup import setup_paths
setup_paths()

import os

import qutip as qt
import pandas as pd
import numpy as np

from typing import Optional
from math import pi
from qutip.measurement import measure, measurement_statistics


def create_spin_operators(N: int):
    #  Konstruujemy odpowiednie operatory działające na odpowiednim qubicie
    sx_list = []
    sz_list = []
    for i in range(N):
        # sigma_x/z na i-tym miejscu, na reszcze identyczności
        op_list_x = [qt.sigmax() if k==i else qt.qeye(2) for k in range(N)]
        op_list_z = [qt.sigmaz() if k==i else qt.qeye(2) for k in range(N)]
        # z listy tworzymy jedna dużą macierz (operator) przez użycie iloczynu tensorowego
        sx_list.append(qt.tensor(op_list_x))
        sz_list.append(qt.tensor(op_list_z))
    return sx_list, sz_list


def create_problem_hamiltonian(N: int, sz_list: list, J: dict, h: Optional[dict] = None):
    #  Zakłądamy że qubity są ponumerowane od 0 do N-1
    linear = sum([h[i] * sz_list[i] for i in range(N)]) if h else 0
    H_problem = linear + sum([J[(i, j)] * sz_list[i] * sz_list[j] for (i, j) in J.keys()])
    return H_problem


def load_dwave_schedule(path: os.PathLike):
    df = pd.read_excel(path, usecols=["s", "A(s) (GHz)", "B(s) (GHz)"])
    # Przejście z GHz na naturalne jednostki oraz aplikacja skalaru 1/2
    df["A(s)"] = df["A(s) (GHz)"].map(lambda x: x*pi)
    df["B(s)"] = df["B(s) (GHz)"].map(lambda x: x*pi)
    return df


# Kod do symulacji wyżarzania kwantowego. W tej symulacji ewoluujemy równanie Schrödingera 
def create_full_hamiltonian(J: dict, h: dict,  time: int, A: Optional[callable] = None, B: Optional[callable] = None):
    if A is None:
        def A(t, T):
            return 1 - t/T

    if B is None:
        def B(t, T):
            return t/T
        
    N = len(h)

    # Tworzymy odpowiednie operatory sigma_x i sigma_z 
    sx_list, sz_list = create_spin_operators(N)
    
    # Tworzymy  Hamiltonian początkowy H_initial 
    H_initial =  -1 * sum(sx_list)

    # Tworzymy hamiltonian probelmu
    H_problem = create_problem_hamiltonian(N, sz_list, J, h)

    # Pełen hamiltonian H(t)
    # H(t) = A(t)*H_initial + B(t)*H_problem
    H_t = qt.QobjEvo([[H_initial, A], [H_problem, B]], args={"T": time})

    # Przygotuj stan podstawowy dla początkowego hamiltonianiu
    psi0 = H_initial.groundstate()[1]

    # Do symulacji potrzebujemy H_t i psi0, H_problem potrzebujemy do wykonania pomiaru
    return H_t, psi0, H_problem


def simulate(hamiltonian, psi0, H_problem, time: int, num_steps: int):
    
    t_list = np.linspace(0, time, num_steps, endpoint=True)
    result = qt.sesolve(hamiltonian, psi0, t_list, e_ops=H_problem, options={"store_states": True})
    states = result.states
    expectatnions = result.expect
    
    return states, expectatnions
    
    

# Ewolucja wartości oczekiwanej w czasie
# Ewolucja prawdopodobienstwa pomiaru stanu w czasie

import matplotlib.pyplot as plt

from qutip.measurement import measure, measurement_statistics
from funkcje_pomocnicze import read_instance_dict, K8


J, h = read_instance_dict(K8.path)

# Czas wyżarzania
T = 50
num_steps = 50

H_t, psi0, H_problem = create_full_hamiltonian(J, h, T)


states, expect = simulate(H_t, psi0, H_problem, T, num_steps)
final_state = states[-1]

# Obrazek wartość oczekiwana
x = range(num_steps)
y1 = expect[0]


plt.figure(figsize=(8, 6))

plt.plot(x, y1, label=r"$\langle \mathcal{H}(t) \rangle$")

plt.axhline(y=-9.5, color='black', linestyle='--')

plt.title("Ewolucja wartości oczekiwanej pomiaru energi w czasie")
plt.xlabel("Czas")
plt.ylabel(r"Wartość oczekiwana energii $\langle \mathcal{H}(t) \rangle$")

plt.text(1, -9.25, "Stan podstawowy", fontsize=10, color='black')
plt.grid(True)
plt.legend()
plt.show()

y2_1 = [measurement_statistics(state, H_problem)[2] for state in states]
y2_2 = [data[:4] for data in y2_1]
y2 = [[x * 100 for x in data] for data in y2_2]

plt.figure(figsize=(8, 6))
plt.plot(x, y2)
plt.title("Ewolucja prawdopodobienstwa pomiaru poszczególnych stanow")
plt.xlabel("Czas")
plt.ylabel("Prawdopodobienstwo pomiaru stanu [%]")
plt.grid(True)
plt.legend([r"Stan podstawowy $H_{\text{problem}}$", r"1. stan wzbudzony $H_{\text{problem}}$", r"2. stan wzbudzony $H_{\text{problem}}$", r"3. stan wzbudzony $H_{\text{problem}}$"])
plt.show()

Tak jak się spodziewaliśmy, wraz z czasem prawdopodobienstwo pomiaru stanu podstawowego rośnie. Pomiar jest probabilistyczny, więc rozpatrujemy wartość oczekiwaną pomiaru.

Demonstracja adiabatyczności

Przeprowadzimy wiele symulacji, każdorazowo wykonując wyżarzanie kwantowe przez podany czas

# rysunek dla demonstracji adiabatyczności

import matplotlib.pyplot as plt
from funkcje_pomocnicze import read_instance_dict, K8
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset


J, h = read_instance_dict(K8.path)
num_steps = 10
T = 50
gaps = []
for t in range(1, T+1):
    H_t, psi0, H_problem = create_full_hamiltonian(J, h, t)
    states, expect = simulate(H_t, psi0, H_problem, t, num_steps)
    energy = expect[0][-1]
    gap = (K8.best_energy - energy)/(2*K8.best_energy) * 100
    gaps.append(gap)

x = range(1, T+1)
fig, ax = plt.subplots(figsize=(8, 6))

plt.plot(x, gaps)
plt.axhline(y=0, color='black', linestyle='--')

# Define the region of interest for the zoomed-in plot
x1, x2 = 30, 40  # x-limits for inset
# Compute corresponding y-limits for the sine function
y1, y2 = gaps[x1], gaps[x2]
# To have a little margin around the data points, adjust y-limits
delta_y = 0.1
y1 -= delta_y
y2 += delta_y

axins = zoomed_inset_axes(ax, zoom=4, loc="right")
axins.plot(x, gaps)

# Limit the view to the region of interest
axins.set_xlim(x1, x2)
axins.set_ylim(0, 1)



# Draw lines indicating the area of the inset on the main plot
mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")

ax.set_title("Luka energetyczna w zależności od całkowitego czasu wyżarzania")
ax.set_xlabel(r"Całkowity czas wyżarzania $\tau$")
ax.set_ylabel("Oczekiwana luka energetyczna [%]")
ax.grid(True)

plt.show()

# Demonstracja luki energetycznej w zależności od całkowitego czasu wyżarzania

import matplotlib.pyplot as plt
from funkcje_pomocnicze import read_instance_dict, K8
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset


J, h = read_instance_dict(K8.path)
num_steps = 25
times = [1, 10, 25, 50]
energies = {}
for t in times:
    H_t, psi0, H_problem = create_full_hamiltonian(J, h, t)
    states, expect = simulate(H_t, psi0, H_problem, t, num_steps)
    energies[t] = expect

x = np.linspace(0, 1, num=num_steps, endpoint=True)
fig, ax = plt.subplots(figsize=(8, 6))

for t in times:
    plt.plot(x, energies[t][0], label=rf"$\tau = ${t}")

plt.text(0, -9.25, "Stan podstawowy", fontsize=10, color='black')
plt.axhline(y=-9.5, color='black', linestyle='--')
plt.legend()
plt.title(r"Ewolucja wartości oczekiwanej pomiaru energi dla różnych wartości $\tau$")
plt.xlabel(r"$s = t / \tau$")
plt.ylabel(r"Wartość oczekiwana energii $\langle \mathcal{H}(t) \rangle$")
plt.grid()
plt.show()

Bibliografia

• Albash, T., and Lidar, D. A. (2018). Adiabatic quantum computation. Reviews of Modern Physics, 90(1), 015002.
• Tanaka, S., Tamura, R., and Chakrabarti, B. K. (2017). Quantum spin glasses, annealing and computation. Cambridge University Press.