Model QUBO¶
Quadratic Unconstrained Binary Optimization (QUBO) to model do rozwiązywania problemów optymalizacji kombinatorycznej. Jak sugeruje nazwa, dotyczy on modeli optymalizacyjnych (O), w których występują zależności kwadratowe (Q) między zmiennymi binarnymi (B), brak ograniczeń (U). Co zaskakujące, warunki QUBO spełnia duża klasa problemów, takich jak MAX-CUT, problem komiwojażera czy harmonogramowanie zadań.
Model QUBO - definicja¶
Dla danej symetrycznej macierzy $Q$, należy znaleźć binarny wektor $\bm{x}^*$, taki że:
$$ \bm{x}^* = \arg \min_x \bm{x}Q\bm{x}^T = \sum_i \sum_j Q_{ij} x_i x_j $$
Często przekształca się macierz $Q$ do formy górnotrójkątnej - dla wszystkich $i$, $j$, gdzie $j > i$, zastępuje się element powyżej przekątnej $Q_{ij}$ przez $Q_{ij} + Q_{ji}$. Następnie wszystkie elementy poniżej przekątnej $Q_{ji}$ dla $j < i$ zamienia się na 0.
Warto zauważyć, że dla zmiennych binarnych zachodzi $x^2 = x$ ($0^2=0$, $1^2=1$). Pozwala to podzielić QUBO na część liniową znajdującą się na diagonali oraz część kwadratową:
$$ \bm{x}Q\bm{x}^T = \sum_i Q_{ii} x_i + \sum_{i \neq j} Q_{ij} x_i x_j $$
Rozwiązywanie QUBO - przykład¶
Podążając za przykładem spróbujemy lepiej zrozumieć, w jaki sposób można przekształcić problemy w formę QUBO. Na początek rozważmy następujący problem optymalizacyjny - chcemy zminimalizować wartość wyrażenia: $$ y = -5x_1 -3x_2 -8x_3 - 6x_4 + 4x_1x_2 + 8x_1x_3 + 2x_2x_3 + 10x_3x_4 $$
Obserwacje:
$y$ składa się z części liniowej $-5x_1 -3x_2 -8x_3 - 6x_4$ oraz części kwadratowej $4x_1x_2 + 8x_1x_3 + 2x_2x_3 + 10x_3x_4$.
Ponieważ dla zmiennych binarnych $x_i = x_i^2$, w części liniowej można zastąpić każdą zmienną jej kwadratem:
$$ -5x_1^2 -3x_2^2 -8x_3^2 - 6x_4^2 $$
co daje postać homogeniczną: $$ y=-5x_1^2 -3x_2^2 -8x_3^2 - 6x_4^2+4x_1x_2 + 8x_1x_3 + 2x_2x_3 + 10x_3x_4 $$
- Następnie możemy przekształcić model do następującej postaci macierzowej:
$$ y = [x_1 \, x_2 \, x_3 \, x_4] \begin{bmatrix} -5 & 4 & 8 & 0 \\ 0 & -3 & 2 & 0 \\ 0 & 0 & -8 & 10 \\ 0 & 0 & 0 & -6 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} $$
lub krócej:
$$
y = \arg \min_x \bm{x}Q\bm{x}^T
$$
Należy zauważyć, że stosujemy tu konwencję zachowania $Q$ jako macierzy górnotrójkątnej.
- Poza wymogiem binarności zmiennych decyzyjnych, QUBO jest modelem nieograniczonym, gdzie wszystkie dane problemu są zawarte w macierzy $Q$. Fakt ten sprawia, że model QUBO jest atrakcyjną formą dla wielu problemów optymalizacyjnych.
# Metota naiwnego wyczerpującego przeszukiwania
from dimod import ExactSolver
import numpy as np
Q = np.array([[-5, 4, 8, 0], [0, -3, 2, 0], [0, 0, -8, 10], [0, 0, 0, -6]])
sampleset = ExactSolver().sample_qubo(Q)
solution_vector = list(sampleset.first.sample.values())
solution = sampleset.first.energy
print("rozwiązanie:", end = " ")
for idx, i in enumerate(solution_vector):
print(f"x_{idx + 1} = {i},", end = " ")
print(f"\ny_min = {solution}")
rozwiązanie: x_1 = 1, x_2 = 0, x_3 = 0, x_4 = 1, y_min = -11.0
Podobieństwo między QUBO a modelem Isinga¶
Nie sposób nie zauważyć, że wyrażenia opisujące QUBO i model Isinga są bliźniaczo podobne
$$ H(\textbf{s}) = \sum_{(i,j)} J_{ij} s_i s_j + \sum_i h_i s_i \quad \text{oraz} \quad \bm{x}Q\bm{x}^T = \sum_{i \neq j} Q_{ij} x_i x_j + \sum_i Q_{ii} x_i $$ gdzie zakładamy, że $\textbf{J}$ i $\textbf{Q}$ są macierzami górnotrójkątnymi, $i = 1, \ldots, N$, a różnicą są wartości $s_i = \pm 1$, $x_i\in\lbrace 0, 1\rbrace$ .
W rzeczywistości rozwiązanie QUBO i znajdowanie stanów podstawowych modelu Isinga to równoważne problemy! Można je stosunkowo łatwo przekształcić między sobą za pomocą następującego podstawienia:
$$
x_i = \frac{s_i + 1}{2}
$$
$$
s_i = 2x_i - 1
$$
Podstawiając te zmienne do odpowiednich wzorów, otrzymujemy:
$$
Q_{ij} = 4 J_{ij}
$$
$$
Q_{ii} = 2h_i - 2\sum_{\mathcal{N}(i)} J_{ij}
$$
$$ \text{offset} = \sum_{(i,j)}J_{ij} - \sum_{i}h_i $$
oraz
$$
J_{ij} = \frac{1}{4} Q_{ij}
$$
$$
h_{i} = \frac{1}{2}Q_{ii} + \frac{1}{4} \sum_{j \neq i} Q_{ij}
$$
$$ \text{offset} = \frac{1}{2}\sum_{i}Q_{ii} + \frac{1}{4}\sum_{i < j}Q_{ij} $$
Należy pamiętać, że transformacja generuje pewną dodatkową stałą, więc uzyskane wartości funkcji celu nie mogą być porównywane bezpośrednio.
Zamiana pomiędzy modelami¶
import numpy as np
rng = np.random.default_rng(seed=42)
# Przykładowy model Isinga
N = 15
h = {i: 0.5 for i in range(N)}
J = {(i, i+1): -1 for i in range(N-1)}
# Ising to qubo
Q = {}
# część liniowa
for i in h.keys():
if i != 0 and i != N-1:
Q[(i,i)] = 2*h[i] - 2*(J[(i,i+1)] + J[(i-1, i)])
elif i == 0:
Q[(i,i)] = 2*h[i] - 2*(J[(i,i+1)])
else:
Q[(i,i)] = 2*h[i] - 2*(J[(i-1,i)])
# część kwadratowa
for (i, j) in J.keys():
Q[(i,j)] = 4 * J[(i,j)]
offset = sum(list(J.values())) - sum(list(h.values()))
state_ising = [rng.choice([-1, 1]) for i in range(N)]
state_qubo = [(s+1)/2 for s in state_ising]
energy_ising = sum([coupling * state_ising[i] * state_ising[j] for (i, j), coupling in J.items()]) + \
sum(bias * state_ising[i] for i, bias in h.items())
energy_qubo = sum([x * state_qubo[i] * state_qubo[j] for (i, j), x in Q.items()])
print("offset: ", offset)
print("energia Isinga: ", energy_ising)
print("energia QUBO: ", energy_qubo)
print("energia QUBO + offset", energy_qubo + offset)
offset: -21.5 energia Isinga: 1.5 energia QUBO: 23.0 energia QUBO + offset 1.5
Na szczęście istnieje wiele narzędzi pozwalających automatycznie przechodzic pomiędzy modelami. Poniżej zostanie zaprezentowane jak to zrobić za pomocą paczki dimod i zawartej w niej klasy BinaryQadraticModel
# Zamiana między modelami.
import numpy as np
from dimod import BinaryQuadraticModel
rng = np.random.default_rng(seed=42)
# Przykładowy model Isinga
N = 15
h = {i: 0.5 for i in range(N)}
J = {(i, i+1): -1 for i in range(N-1)}
# Utworzenie BQM
bqm = BinaryQuadraticModel(h, J, vartype="SPIN")
# zamiana modelu Isinga na QUBO
# Warto pamiętać, że zamiana daje nam zarówno model w nowej formie jak i tzw. offset czyli stałą którą musimy brac pod uwagę przy zamienianiu modeli
QUBO, offset = bqm.to_qubo()
state_ising = [rng.choice([-1, 1]) for i in range(N)]
state_qubo = [(s+1)/2 for s in state_ising]
energy_ising = sum([coupling * state_ising[i] * state_ising[j] for (i, j), coupling in J.items()]) + \
sum(bias * state_ising[i] for i, bias in h.items())
energy_qubo = sum([x * state_qubo[i] * state_qubo[j] for (i, j), x in QUBO.items()])
print("offset: ", offset)
print("energia Isinga: ", energy_ising)
print("energia QUBO: ", energy_qubo)
print("energia QUBO + offset", energy_qubo + offset)
# Tak samo można zamienić QUBO na Isinga
# ćwiczenie dla grupy
offset: -21.5 energia Isinga: 1.5 energia QUBO: 23.0 energia QUBO + offset 1.5
Przykład - dyskretny problem plecakowy¶
Dyskretny problem plecakowy (ang. 0-1 Knapsack Problem) należy do klasycznych problemów optymalizacji kombinatorycznej. Dane są:
- zbiór przedmiotów $i = 1,2,\dots,n$,
- wartość każdego przedmiotu $v_i$,
- waga każdego przedmiotu $w_i$,
- pojemność plecaka $W$.
Celem jest maksymalizacja łącznej wartości wybranych przedmiotów przy zachowaniu ograniczenia pojemności plecaka:
Formalna definicja¶
Dla każdego przedmiotu podejmujemy decyzję binarną:
$$ x_i \in \{0,1\}, $$
gdzie:
- $x_i = 1$ oznacza, że przedmiot został wybrany,
- $x_i = 0$ oznacza, że przedmiot nie został wybrany.
$$ \max \sum_{i=1}^{n} v_i x_i $$
przy warunku:
$$ \sum_{i=1}^{n} w_i x_i \leq W. $$
Jest to więc problem wyboru takiego podzbioru przedmiotów, który daje jak największą wartość, ale nie przekracza dopuszczalnej wagi.
Przejście do formulacji QUBO¶
Model QUBO (Quadratic Unconstrained Binary Optimization) wymaga zapisania problemu w postaci funkcji kwadratowej zmiennych binarnych bez jawnych ograniczeń. Oznacza to, że ograniczenie pojemności trzeba włączyć do funkcji celu w postaci składnika karnego.
Standardowa postać QUBO ma postać:
$$ \min \; x^T Q x, \qquad x \in \{0,1\}^m. $$
Ponieważ problem plecakowy jest sformułowany jako zadanie maksymalizacyjne z ograniczeniem, przepisujemy go do postaci minimalizacyjnej i dodajemy karę za naruszenie ograniczenia.
Wprowadzenie zmiennych pomocniczych¶
Aby zamienić nierówność
$$ \sum_{i=1}^{n} w_i x_i \leq W $$
na równość, wprowadzamy zmienną luzu $ s $, taką że:
$$ \sum_{i=1}^{n} w_i x_i + s = W, \qquad s \geq 0. $$
Ponieważ w modelu QUBO wszystkie zmienne muszą być binarne, zmienną $ s $ kodujemy binarnie:
$$ s = \sum_{k=0}^{K} 2^k y_k, \qquad y_k \in \{0,1\}. $$
Liczba bitów $K+1$ powinna być wystarczająca do reprezentacji wartości od $0$ do $W$.
Po podstawieniu otrzymujemy warunek:
$$ \sum_{i=1}^{n} w_i x_i + \sum_{k=0}^{K} 2^k y_k = W. $$
Przekształcając:
$$ \sum_{i=1}^{n} w_i x_i + \sum_{k=0}^{K} 2^k y_k - W = 0. $$¶
Funkcja celu z karą¶
Chcemy maksymalizować wartość przedmiotów, więc w wersji minimalizacyjnej używamy składnika:
$$ - \sum_{i=1}^{n} v_i x_i. $$
Aby wymusić spełnienie ograniczenia, dodajemy karę kwadratową:
$$ A \left( \sum_{i=1}^{n} w_i x_i + \sum_{k=0}^{K} 2^k y_k - W \right)^2, $$
gdzie $A > 0$ jest współczynnikiem kary (często też używa się oznaczenia $P$).
Ostatecznie otrzymujemy funkcję QUBO:
$$ \min \; A \left( \sum_{i=1}^{n} w_i x_i + \sum_{k=0}^{K} 2^k y_k - W \right)^2 - B \sum_{i=1}^{n} v_i x_i, $$
gdzie:
- $A$ odpowiada za egzekwowanie ograniczenia pojemności,
- $B$ wzmacnia wpływ maksymalizacji wartości.
Zazwyczaj wybiera się $A$ istotnie większe od $B$ $(A \gg B)$, aby rozwiązania niezgodne z ograniczeniem były silnie penalizowane.
items = [
{"name": "A", "weight": 4, "value": 8},
{"name": "B", "weight": 5, "value": 10},
{"name": "C", "weight": 2, "value": 5},
{"name": "D", "weight": 3, "value": 6},
]
W = 9 # pojemność plecaka
weights = [item["weight"] for item in items]
values = [item["value"] for item in items]
n = len(items)
# Współczynniki modelu
A = 20 # kara za niespełnienie ograniczenia
B = 1
Bibliografia¶
- Lucas, A. (2014). Ising formulations of many NP problems. Frontiers in Physics, Volume 2
- Glover, F., Kochenberger, G., Hennig, R. et al. (2022). Quantum bridge analytics I: a tutorial on formulating and using QUBO models. Annals of Operations Research 314, 141–183