diff --git a/qmlcourseRU/_bibliography/references.bib b/qmlcourseRU/_bibliography/references.bib index c1b292c1c7..11222a9222 100644 --- a/qmlcourseRU/_bibliography/references.bib +++ b/qmlcourseRU/_bibliography/references.bib @@ -398,6 +398,21 @@ @article{diagonalQNewton publisher = {Springer} } +@article{embeddingqubo2019, + title = {Efficiently embedding QUBO problems on adiabatic quantum computers}, + author = {Date, Prasanna and Patton, Robert and Potok, Thomas}, + journal = {Quantum Information Processing}, + volume = {18}, + issue = {4}, + pages = {117}, + numpages = {0}, + year = {2019}, + month = {Mar}, + publisher = {Quantum Information Processing}, + doi = {10.1007/s11128-019-2236-3}, + url = {https://doi.org/10.1007/s11128-019-2236-3} +} + @misc{farhi2000quantum, title = {Quantum computation by adiabatic evolution}, author = {Farhi, Edward and Goldstone, Jeffrey and Gutmann, Sam and Sipser, Michael}, @@ -587,6 +602,21 @@ @article{jastrow1955many publisher = {APS} } +@article{kadowaki1998, + title = {Quantum annealing in the transverse Ising model}, + author = {Kadowaki, Tadashi and Nishimori, Hidetoshi}, + journal = {Phys. Rev. E}, + volume = {58}, + issue = {5}, + pages = {5355--5363}, + numpages = {0}, + year = {1998}, + month = {Nov}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevE.58.5355}, + url = {https://link.aps.org/doi/10.1103/PhysRevE.58.5355} +} + @article{lanczos1950iteration, title = {An iteration method for the solution of the eigenvalue problem of linear differential and integral operators}, author = {Lanczos, Cornelius}, @@ -747,6 +777,21 @@ @article{qsvmmap publisher = {Springer} } +@article{qubo3mlproblems2021, + title = {QUBO formulations for training machine learning models}, + author = {Date, Prasanna and Arthur, Davis and Pusey-Nazzaro, Lauren}, + journal = {Scientific Reports}, + volume = {11}, + issue = {1}, + pages = {10029}, + numpages = {0}, + year = {2021}, + month = {May}, + publisher = {Scientific Reports}, + doi = {10.1038/s41598-021-89461-4}, + url = {https://doi.org/10.1038/s41598-021-89461-4} +} + @article{seeley2012bravyi, title = {The Bravyi-Kitaev transformation for quantum computation of electronic structure}, author = {Seeley, Jacob T. and Richard, Martin J. and Love, Peter J.}, @@ -759,6 +804,21 @@ @article{seeley2012bravyi url = {https://arxiv.org/abs/1208.5986} } +@article{sherrington1975, + title = {Solvable Model of a Spin-Glass}, + author = {Sherrington, David and Kirkpatrick, Scott}, + journal = {Phys. Rev. Lett.}, + volume = {35}, + issue = {26}, + pages = {1792--1796}, + numpages = {0}, + year = {1975}, + month = {Dec}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevLett.35.1792}, + url = {https://link.aps.org/doi/10.1103/PhysRevLett.35.1792} +} + @article{sinchenko2019deep, title = {The deep learning and statistical physics applications to the problems of combinatorial optimization}, author = {Sinchenko, Semyon and Bazhanov, Dmitry}, @@ -904,33 +964,3 @@ @article{zachary1977 year = 1977, publisher = {University of New Mexico} } - -@article{kadowaki1998, - title = {Quantum annealing in the transverse Ising model}, - author = {Kadowaki, Tadashi and Nishimori, Hidetoshi}, - journal = {Phys. Rev. E}, - volume = {58}, - issue = {5}, - pages = {5355--5363}, - numpages = {0}, - year = {1998}, - month = {Nov}, - publisher = {American Physical Society}, - doi = {10.1103/PhysRevE.58.5355}, - url = {https://link.aps.org/doi/10.1103/PhysRevE.58.5355} -} - -@article{sherrington1975, - title = {Solvable Model of a Spin-Glass}, - author = {Sherrington, David and Kirkpatrick, Scott}, - journal = {Phys. Rev. Lett.}, - volume = {35}, - issue = {26}, - pages = {1792--1796}, - numpages = {0}, - year = {1975}, - month = {Dec}, - publisher = {American Physical Society}, - doi = {10.1103/PhysRevLett.35.1792}, - url = {https://link.aps.org/doi/10.1103/PhysRevLett.35.1792} -} diff --git a/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_1.png b/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_1.png new file mode 100644 index 0000000000..190b968d99 Binary files /dev/null and b/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_1.png differ diff --git a/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_2.png b/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_2.png new file mode 100644 index 0000000000..01a907fb90 Binary files /dev/null and b/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_2.png differ diff --git a/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_3.png b/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_3.png new file mode 100644 index 0000000000..9e123f4a14 Binary files /dev/null and b/qmlcourseRU/_static/problems2qml/qubo_linreg_svm_kmeans/fig_3.png differ diff --git a/qmlcourseRU/_toc.yml b/qmlcourseRU/_toc.yml index 4972d1229b..cfad6b545a 100644 --- a/qmlcourseRU/_toc.yml +++ b/qmlcourseRU/_toc.yml @@ -168,6 +168,8 @@ title: ⬜ О блоке - file: book/problems2qml/np2ising title: 🟩 NP как модели Изинга + - file: book/problems2qml/qubo_linreg_svm_kmeans + title: 🟧 QUBO-формулировки трех задач ML - file: book/problems2qml/jordanwigner title: 🟧 Преобразование Жордана-Вигнера - file: book/problems2qml/eigenvals diff --git a/qmlcourseRU/book/problems2qml/qubo_linreg_svm_kmeans.md b/qmlcourseRU/book/problems2qml/qubo_linreg_svm_kmeans.md new file mode 100644 index 0000000000..85be802147 --- /dev/null +++ b/qmlcourseRU/book/problems2qml/qubo_linreg_svm_kmeans.md @@ -0,0 +1,511 @@ +(qubo_linreg_svm_kmeans)= + +# QUBO-формулировки для линейной регрессии, SVM и метода k-средних + +Автор(ы): + +- [Бурдейный Дмитрий](https://github.com/dmburd) + + +## Описание лекции + +Здесь рассмотрим QUBO-формулировки трех задач ML (линейная регрессия, SVM, сбалансированная кластеризация методом k-средних) аналогично тому, как в [предыдущей лекции](np2ising) было продемонстрировано сведение к QUBO задач о максимальном разрезе в графе, о коммивояжере и о выделении сообществ в графе. Изложение полностью основывается на статье {cite}`qubo3mlproblems2021`. + +Некоторые обозначения, которые будут использоваться в дальнейшем: + +- $\mathbb{R}$ -- множество действительных чисел; +- $N$ -- количество объектов в обучающем наборе, $N \in \{1, 2, \dotsc\}$; +- $d$ -- количество признаков (*features*) для объектов в обучающем наборе, $d \in \{1, 2, \dotsc\}$; +- $X$ -- тренировочный набор, содержит по одному объекту в каждой из $N$ строк, каждая строка содержит $d$ значений признаков, $X \in \mathbb{R}^{N \times d}$; +- $Y$ -- набор истинных "ответов" (*ground truth*), соответствующих тренировочным объектам из $X$ ($Y \in \mathbb{R}^{N}$ в случае регрессии, $Y \in \{0,1\}^{N}$ в случае бинарной классификации); +- $\otimes$ -- произведение Кронекера (тензорное произведение); +- $\odot$ -- произведение Адамара (поэлементное произведение). + +Общая формулировка QUBO-задачи, которая используется в статье {cite}`qubo3mlproblems2021` и к которой всё сводится, выглядит так: + +$$ +z^T A z + z^T b + \to \min_{z \in \{0,1\}^M} +$$ (eqn:general_qubo_in_paper) + +где $M$ -- натуральное число, $z$ -- бинарный вектор решения, $A \in \mathbb{R}^{M \times M}$ -- QUBO-матрица с действительными элементами, $b \in \mathbb{R}^{M}$ -- QUBO-вектор. Как отмечалось в предыдущей лекции, при $z_i \in \{0,1\}$ выполняется равенство $z_i^2 = z_i$, так что линейные члены в {eq}`eqn:general_qubo_in_paper` можно включить в квадратичные, но этого делать не будем, т.к. для целей этой лекции и для лучшего понимания удобнее сохранить минимизируемую квадратичную форму именно в виде {eq}`eqn:general_qubo_in_paper`. + +## Линейная регрессия + +В задаче линейной регрессии предполагается, что зависимость истинных ответов от признаков тренировочных объектов приближенно линейная: + +$$ +y_i^{(pred)} = \braket{w, x_i} + b, \quad y_i^{(pred)} \approx y_i, +$$ + +где $y_i^{(pred)}$ -- предсказываемое значение, $y_i$ -- истинное значение (из разметки), $\braket{\cdot, \cdot}$ -- скалярное произведение. Удобно сразу избавиться от слагаемого $b$ (*bias*), добавив единицу к набору признаков. Тогда bias окажется включенным в веса $w$, а тренировочный набор будет иметь по $(d+1)$ признаков на объект: $X \in \mathbb{R}^{N \times (d+1)}$. Требуется найти веса $w$, при которых квадрат евклидовой нормы невязки минимален: + +$$ +E(w) = || Xw - Y ||^2 + \to \min_{w \in \mathbb{R}^{d+1}} +$$ (eqn:lin_reg_orig) + +```{figure} /_static/problems2qml/qubo_linreg_svm_kmeans/fig_1.png +:name: lin_reg_example +:width: 500px + +Иллюстрация к задаче линейной регрессии. +``` + + +Известно аналитическое решение задачи {eq}`eqn:lin_reg_orig`: + +$$ +w = \left( X^T X \right)^{-1} X^T Y +$$ + +Если $\left( X^T X \right)^{-1}$ не существует, нужно вычислить [псевдообратную](https://ru.wikipedia.org/wiki/Псевдообратная_матрица) матрицу. ВременнАя сложность решения задачи линейной регрессии равна $\mathcal{O}(N d^2)$, т.к. нужно $\mathcal{O}(N d^2)$ для вычисления матрицы $X^T X$ и $\mathcal{O}(d^3)$ для вычисления обратной к ней (предполагаем, что $N \gg d$). + + +### QUBO-формулировка + +Перепишем выражение {eq}`eqn:lin_reg_orig`: + +$$ +E(w) = w^T X^T X w - 2 w^T X^T Y + Y^T Y + \to \min_{w \in \mathbb{R}^{d+1}} +$$ (eqn:lin_reg_qubo_1) + +Наша цель -- найти вектор $w$, компоненты которого -- действительные числа. Но в QUBO-формулировке необходимо представить решение в виде вектора с бинарными компонентами. Как это сделать? Конечно, напрашивается идея использовать бинарное представление действительного числа $w_i$ (будем отдельно записывать в бинарном виде целую часть, отдельно -- дробную). Нужно помнить о том, что знак $w_i$ может быть как положительным, так и отрицательным. Формат представления придется выбрать фиксированным (т.е. с [фиксированной запятой](https://ru.wikipedia.org/wiki/Число_с_фиксированной_запятой), не с плавающей запятой). Пример: + +$$ +\begin{align} + \pm 110.101 &= \pm + \underbrace{ + \left( 0 \cdot 2^0 + 1 \cdot 2^1 + 1 \cdot 2^2 \right) + }_{\text{идем влево от разд. точки}} + \pm + \underbrace{ + \left( 1 \cdot 2^{-1} + 0 \cdot 2^{-2} + 1 \cdot 2^{-3} \right) + }_{\text{идем вправо от разд. точки}} \\ + &= \pm \left( 2 + 4 + \frac{1}{2} + \frac{1}{8} \right) \\ + &= \pm 6.625 +\end{align} +$$ + +Бинарные компоненты логично рассматривать как индикаторы наличия или отсутствия соответствующих степеней двойки в бинарном представлении каждого действительного числа. Введем вектор-столбец $P$, который отвечает за точность представления (_**P**recision vector_) и состоит из степеней двойки со знаками: + +$$ +P = \left[ -2^l, -2^{l-1}, \dotsc, -2^{-m+1}, -2^{-m}, 2^{-m}, 2^{-m+1}, \dotsc, 2^{l-1}, 2^l \right]^T, +$$ + +этот вектор отсортирован по возрастанию элементов. Отводится $l$ двоичных разрядов для целой части числа, $m$ разрядов для дробной. $l$ определяется максимальным по модулю действительным числом, которое хотим представлять; $m$ определяется желаемой точностью представления. Число элементов вектора $P$ равно $(m + 1 + l) \cdot2 = K$. + +Вводим вектор $\tilde{w}_i \in \{0,1\}^K$ такой, что + +$$ +\tilde{w}_i^T P += \sum_{k=1}^{K} p_k \tilde{w}_{i k} +\approx w_i \in \mathbb{R} +$$ + +Чтобы не было неоднозначности, нужно договориться о том, что для представления $w_i > 0$ используем только положительные элементы вектора $P$, а для $w_i < 0$ -- только отрицательные. + +Составляем бинарный вектор $\tilde{w} \in \{0,1\}^{K(d+1)}$ + +$$ +\tilde{w} = [ + \underbrace{ \tilde{w}_{11} \dots \tilde{w}_{1K} }_{\text{предст. } w_1} + \underbrace{ \tilde{w}_{21} \dots \tilde{w}_{2K} }_{\text{предст. } w_2} + \dots + \underbrace{ \tilde{w}_{(d+1)1} \dots \tilde{w}_{(d+1)K} }_{\text{предст. } w_{d+1}} + ]^T +$$ + +и матрицу точности $\mathcal{P}$ (_**P**recision matrix_), которая задает переход к бинарному представлению векторов: + +$$ +\mathcal{P} = I_{d+1} \otimes P^T, +$$ + +где $I_{d+1}$ -- единичная матрица размера $(d+1)$. Матрица $\mathcal{P}$ имеет размерность $(d+1) \times K(d+1)$. Исходный вектор весов можно записать как + +$$ +w = \mathcal{P} \tilde{w}, +$$ (eqn:w_by_tilde_w) + +где знак "=" на самом деле означает приближенное равенство (наше fixed-point представление имеет конечную точность). + +Всё готово для того, чтобы переписать исходную задачу {eq}`eqn:lin_reg_qubo_1` в QUBO-формулировке. Подставляем {eq}`eqn:w_by_tilde_w` в {eq}`eqn:lin_reg_qubo_1` и получаем + +$$ +E(\tilde{w}) = + \tilde{w}^T \mathcal{P}^T X^T X \mathcal{P} \tilde{w} + - 2 \tilde{w}^T \mathcal{P}^T X^T Y +\to \min_{\tilde{w} \in \{0,1\}^{(d+1)K}} +$$ (eqn:lin_reg_qubo_2) + +слагаемое $Y^T Y$ отброшено, т.к. это константа, никак не влияющая на решение задачи оптимизации без ограничений. + +### Оценка вычислительной сложности + +Для исходной задачи регрессии {eq}`eqn:lin_reg_orig` количество значений в датасете $X$ равно $\mathcal{O}(N d)$. Мы ввели $K$ бинарных переменных для каждого из $(d+1)$ весов. Значит, получилось $\mathcal{O}(K d)$ переменных в QUBO-формулировке {eq}`eqn:lin_reg_qubo_2`. Для решения задачи требуется $\mathcal{O}(K^2 d^2)$ кубитов (см. {cite}`embeddingqubo2019`), это пространственная сложность в рассматриваемом подходе. + +ВременнАя сложность в классической задаче $\mathcal{O}(N d^2)$. В случае QUBO-задачи для временнОй оценки нужно рассмотреть три части: + +1. Затраты времени для конвертации задачи регрессии в QUBO-формулировку. Здесь получаем $\mathcal{O}(N K^2 d^2)$ (проверьте это, оценив число умножений, необходимых для вычисления $\mathcal{P}^T X^T X \mathcal{P}$). +2. Время для реализации QUBO-задачи в квантовом "железе". Здесь потребуется $\mathcal{O}(K^2 d^2)$, если использовать алгоритм из {cite}`embeddingqubo2019`. +3. Время для выполнения квантового отжига. Существуют теоретические оценки времени для получения точного решения, но более практично рассматривать случай, когда можно просто довольствоваться достаточно высокой вероятностью (скажем, 99%) получения оптимального решения. Для современных квантовых компьютеров D-Wave с ограниченным числом кубитов на практике получается, что время отжига и число повторений можно считать константами. + +В итоге полная временнАя сложность решения QUBO-задачи на адиабатическом квантовом компьютере $\mathcal{O}(N K^2 d^2)$. Может показаться, что это хуже, чем временнАя сложность классического решения, если считать $K$ переменной. Но величина $K$ определяется только шириной диапазона числовых значений и желаемой точностью представления, $K$ не зависит от основных параметров задачи типа $N$ и $d$. Поэтому можно считать $K$ константой. Тогда число требуемых кубитов $\mathcal{O}(d^2)$, временнАя сложность $\mathcal{O}(N d^2)$, это эквивалентно классическому случаю. + + +## SVM + +Классический SVM подробно описан в соответствующей [лекции](classicsvm). Рассматривается тренировочный набор $X \in \mathbb{R}^{N \times d}$ и набор истинных меток $Y \in \{-1, +1\}^{N}$. Нужно решить задачу бинарной классификации, найдя веса $w \in \mathbb{R}^{d}$ и константу $b \in \mathbb{R}$, при которых классификатор $a(x_i) = \text{sign} \left( w^T x_i + b \right)$ допускает как можно меньше ошибок на обучающей выборке. + +```{figure} /_static/problems2qml/qubo_linreg_svm_kmeans/fig_2.png +:name: svm_example +:width: 500px + +Иллюстрация к задаче бинарной классификации с помощью SVM. +``` + +Двойственная задача {eq}`svmDual` в текущих обозначениях принимает вид + +$$ +\left\{ + \begin{aligned} + & \mathcal{L}(\lambda) = + \sum_{i = 1}^{N} \lambda_i + -\frac{1}{2} \sum_{i, j = 1}^{N} + \lambda_i \lambda_j y_i y_j + \langle x_i, x_j \rangle + \to \max_{\lambda}, \\ + & 0 \leq \lambda_i \leq C \quad \forall i \in \{1, 2, \dotsc, N\}, \\ + & \sum_{i = 1}^{N} \lambda_i y_i = 0. + \end{aligned} +\right. +$$ (eqn:svm_dual_1) + +### QUBO-формулировка + +Перепишем {eq}`eqn:svm_dual_1` как задачу минимизации: + +$$ +\left\{ + \begin{aligned} + & \mathcal{L}_{neg}(\lambda) = + \frac{1}{2} \sum_{i, j = 1}^{N} + \lambda_i \lambda_j y_i y_j + \langle x_i, x_j \rangle + -\sum_{i = 1}^{N} \lambda_i + \to \min_{\lambda}, \\ + & 0 \leq \lambda_i \leq C \quad \forall i \in \{1, 2, \dotsc, N\}, \\ + & \sum_{i = 1}^{N} \lambda_i y_i = 0. + \end{aligned} +\right. +$$ (eqn:svm_dual_2) + +Тренировочные объекты $x_i$, соответствующие $\lambda_i = 0$, называются *периферийными*, от них решение + +$$ +w = \sum_{i = 1}^{N} \lambda_i y_i x_i +$$ + +не зависит. Объекты, соответствующие $\lambda_i > 0$, называются *опорными*. + +Задача минимизации переписывается в виде + +$$ +\left\{ + \begin{aligned} + & \mathcal{L}_{neg}(\lambda) = + \frac{1}{2} \lambda^T \left( X X^T \odot Y Y^T \right) \lambda + -\lambda^T 1_{N} + \to \min_{\lambda}, \\ + & 0_{N} \leq \lambda \leq C_{N}, \\ + & \lambda^T Y = 0, + \end{aligned} +\right. +$$ (eqn:svm_dual_3) + +где $1_N$ -- вектор, состоящий из $N$ единиц (аналогичный смысл имеют $0_N$ и $C_N$), $\lambda \in \mathbb{R}^{N}$. + +Как в задаче линейной регрессии, будем представлять каждую $\lambda_i \in \mathbb{R}$ бинарным вектором $\tilde{\lambda}_i \in \{0,1\}^K$. Поскольку $\lambda_i \geq 0$, в *precision vector* $P$ нужно включить только положительные значения: + +$$ +P = \left[2^{-m}, 2^{-m+1}, \dotsc, 2^{l-1}, 2^l \right]^T +$$ + +Вектор $P$ содержит $m + 1 + l = K$ элементов (здесь $K$ никак не связано с $K$ из раздела про линейную регрессию, просто было решено не вводить новое обозначение для длины нового вектора $P$). При подходящем выборе $m$ и $l$ достаточно точно выполняется равенство + +$$ +\lambda_i \approx \sum_{k=1}^{K} p_k \tilde{\lambda}_{i k} +\quad \forall i \in \{1, 2, \dotsc, N\} +$$ (eqn:real_lambda_by_binary) + +Выполняем конкатенацию всех $\tilde{\lambda}_{i k}$ по вертикали + +$$ +\tilde{\lambda} = [ + \underbrace{ \tilde{\lambda}_{11} \dots \tilde{\lambda}_{1K} }_{\text{предст. } \lambda_1} + \underbrace{ \tilde{\lambda}_{21} \dots \tilde{\lambda}_{2K} }_{\text{предст. } \lambda_2} + \dots + \underbrace{ \tilde{\lambda}_{N1} \dots \tilde{\lambda}_{NK} }_{\text{предст. } \lambda_{N}} + ]^T +$$ + +и вводим матрицу $\mathcal{P}$ (*precision matrix*) + +$$ +\mathcal{P} = I_{N} \otimes P^T +$$ + +таким образом, что (приближенно) + +$$ +\lambda = \mathcal{P} \tilde{\lambda} +$$ + +Подставляя из этого выражения $\lambda$ в {eq}`eqn:svm_dual_3`, получаем: + +$$ +\left\{ + \begin{aligned} + & \mathcal{L}_{neg}(\tilde{\lambda}) = + \frac{1}{2} + \tilde{\lambda}^T \mathcal{P}^T \left( X X^T \odot Y Y^T \right) \mathcal{P} \tilde{\lambda} + -\tilde{\lambda}^T \mathcal{P}^T 1_{N} + \to \min_{\tilde{\lambda} \in \{0,1\}^{N K}} \\ + & \left( \mathcal{P} \tilde{\lambda} \right)^T Y = 0 + \end{aligned} +\right. +$$ (eqn:svm_qubo_1) + +Остается избавиться от ограничения в виде равенства в {eq}`eqn:svm_qubo_1`. Как обычно делается в таких случаях, вместо ограничения вводим соответствующий штраф за его нарушение: + +$$ +Penalty^{(hp)} = \frac{\gamma}{2} + \left( + \left( \mathcal{P} \tilde{\lambda} \right)^T Y + \right)^2 + = \frac{\gamma}{2} \left( + \left( \tilde{\lambda}^T \mathcal{P}^T \right) Y + \right)^2 + = \frac{\gamma}{2} \tilde{\lambda}^T \mathcal{P}^T \left( Y Y^T \right) \mathcal{P} \tilde{\lambda}, +$$ + +где $\gamma$ -- достаточно большая константа, *hp* означает *hyperplane* (гиперплоскость). + +Добавляем $Penalty^{(hp)}$ к $\mathcal{L}_{neg}(\tilde{\lambda})$ и получаем итоговую QUBO-формулировку: + +$$ + \frac{1}{2} + \tilde{\lambda}^T \mathcal{P}^T + \left( + X X^T \odot Y Y^T + \gamma Y Y^T + \right) + \mathcal{P} \tilde{\lambda} + -\tilde{\lambda}^T \mathcal{P}^T 1_{N} + \to \min_{\tilde{\lambda} \in \{0,1\}^{N K}} +$$ (eqn:svm_qubo_2) + + + +### Оценка вычислительной сложности + +Задача {eq}`eqn:svm_dual_3` содержит $\mathcal{O}(N d)$ значений в данных и $\mathcal{O}(N)$ параметров ($\lambda$). QUBO-формулировка {eq}`eqn:svm_qubo_2` содержит то же количество данных, а число параметров в $K$ раз больше, т.е. $\mathcal{O}(K N)$. Значит, потребуется $\mathcal{O}(N^2 K^2)$ кубитов. + +ВременнАя сложность классического SVM в типичных реализациях (например, [LIBSVM](https://en.wikipedia.org/wiki/LIBSVM)) равна $\mathcal{O}(N^3)$. Для оценки временнОй сложности QUBO рассматриваем три составляющие (как в задаче линейной регрессии): + +1. Затраты времени для конвертации в QUBO-формулировку. Из {eq}`eqn:svm_dual_2` и {eq}`eqn:real_lambda_by_binary` следует, что оценка времени $\mathcal{O}(N^2 K^2)$. +2. Для реализации QUBO-задачи в квантовом "железе" потребуется $\mathcal{O}(N^2 K^2)$. +3. Время для выполнения квантового отжига и число повторений можно считать константами (см. комментарии к тому же пункту в обсуждении линейной регрессии). + +В итоге временнАя сложность $\mathcal{O}(N^2 K^2)$. $K$ можно считать константой, т.к. она зависит только от диапазона и желаемой точности представления $\lambda$ и не зависит от параметров самой задачи классификации. Тогда получается временная сложность $\mathcal{O}(N^2)$, что гораздо лучше оценки $\mathcal{O}(N^3)$ в классическом случае. + + +## Сбалансированная кластеризация методом k-средних + +Кластеризация методом k-средних -- ML-задача обучения без учителя (*unsupervised*). Требуется распределить тренировочные объекты по $k$ кластерам так, чтобы суммарное отклонение тренировочных объектов, принадлежащих кластерам, от центроидов (центров масс) соответствующих кластеров было минимальным. Сбалансированная кластеризация методом k-средних -- частный случай, в котором каждый кластер содержит примерно одно и то же количество объектов $N/k$, как показано на {numref}`kmeans_clustering_example`. + +```{figure} /_static/problems2qml/qubo_linreg_svm_kmeans/fig_3.png +:name: kmeans_clustering_example +:width: 500px + +Иллюстрация к задаче сбалансированной кластеризации методом k-средних. +``` + +Нужно распределить $N$ объектов из тренировочного набора $X \in \mathbb{R}^{N \times d}$ по $k$ кластерам $\Phi = \{\phi_1, \dotsc, \phi_k\}$. Пусть $\mu_i$ -- центроид кластера $\phi_i$. В общем случае задача кластеризации методом k-средних формулируется так: + +$$ +\sum_{i=1}^{k} \frac{1}{2 |\phi_i|} + \sum_{x,y \in \phi_i} || x - y ||^2 + \to \min_{\Phi} +$$ (eqn:general_kmeans) + +Если размеры $|\phi_i|$ всех кластеров одинаковые, то формулировка переписывается так: + +$$ +\sum_{i=1}^{k} \sum_{x,y \in \phi_i} || x - y ||^2 + \to \min_{\Phi} +$$ (eqn:kmeans_equal_clusters) + +В прикладных задачах кластеризации размеры кластеров только приближенно равны друг другу, поэтому решение задачи {eq}`eqn:kmeans_equal_clusters` не является точным решением задачи {eq}`eqn:general_kmeans`. Для решения задачи кластеризации методом k-средних можно использовать, например, [алгоритм Ллойда](https://en.wikipedia.org/wiki/Lloyd%27s_algorithm). Существует модификация алгоритма Ллойда для случая сбалансированной кластеризации. + +### QUBO-формулировка + +Вводим матрицу $D \in \mathbb{R}^{N \times N}$, элементы которой равны квадратам попарных расстояний между тренировочными объектами: + +$$ +d_{i j} = || x_i - x_j ||^2 +$$ + +Также вводим бинарную матрицу $\tilde{W} \in \{0,1\}^{N \times k}$, каждый элемент которой $\tilde{w}_{ij} = 1$ в том и только в том случае, когда объект $x_i$ принадлежит кластеру $\phi_j$. Очевидно, есть два ограничения на $\tilde{W}$: + +1. Поскольку мы предполагаем, что все кластеры содержат примерно одно и тоже количество объектов, каждый столбец $\tilde{W}$ должен содержать примерно $N/k$ единиц. +2. Каждый объект принадлежит ровно одному кластеру, поэтому каждая строка $\tilde{W}$ должна содержать ровно одну единицу. + +Для получения QUBO-формулировки задачи нам потребуется избавиться от этих ограничений. Для этого, как обычно, введем в минимизируемую квадратичную форму штрафы за нарушение ограничений. Вернемся к этому через пару абзацев. + +Используя $D$ и $\tilde{W}$, мы можем переписать внутреннюю сумму в {eq}`eqn:kmeans_equal_clusters` в виде + +$$ +\sum_{x,y \in \phi_i} || x - y ||^2 = + {\tilde{w}_{j}^{'}}^T D \tilde{w}_{j}^{'}, +$$ + +где $\tilde{w}_{j}^{'}$ -- столбец номер j в $\tilde{W}$. Чтобы переписать в бинарных переменных полную (двойную) сумму в {eq}`eqn:kmeans_equal_clusters`, составим вектор-столбец из всех $N k$ элементов матрицы $\tilde{W}$: + +$$ +\tilde{w} = [ + \tilde{w}_{11} \dots \tilde{w}_{N1} + \tilde{w}_{12} \dots \tilde{w}_{N2} + \dots + \tilde{w}_{1k} \dots \tilde{w}_{Nk} + ]^T +$$ + +При условии, что ограничения на $\tilde{w}$ выполнены, запишем задачу {eq}`eqn:kmeans_equal_clusters` в эквивалентном виде + +$$ +\tilde{w}^T \left( I_k \otimes D \right) \tilde{w} + \to \min_{\tilde{w}} +$$ (eqn:sqr_form_without_constraints) + +Теперь разбираемся с ограничениями на $\tilde{w}$. + +Во-первых, каждый столбец $\tilde{W}$ должен содержать примерно $N/k$ единиц. Введем штраф, непосредственно отражающий это требование: + +$$ +Penalty_{j}^{(col)} = + \alpha + \left( + {\tilde{w}_{j}^{'}}^T \tilde{w}_{j}^{'} - \frac{N}{k} + \right)^2, +$$ + +где $\alpha$ -- достаточно большая константа. Раскрыв скобки в предыдущем выражении, можно убедиться, что + +$$ +Penalty_{j}^{(col)} + = + {\tilde{w}_{j}^{'}}^T + \alpha + \underbrace{ + \left( 1_N - \frac{2N}{k} I_N \right) + }_{\text{обозначим как } F} + \tilde{w}_{j}^{'} + + \text{const} +$$ + +Сумма всех штрафов для столбцов равна + +$$ +Penalty^{(col)} = \sum_j Penalty_{j}^{(col)} = + \tilde{w}^T + \left( I_k \otimes \alpha F \right) + \tilde{w} +$$ + +Во-вторых, каждая строка $\tilde{W}$ должна содержать ровно одну единицу. Соответствующий штраф + +$$ +Penalty_{i}^{(row)} = + \beta + \left( + \tilde{w}_{i}^T \tilde{w}_{i} - 1 + \right)^2, +$$ + +где $\beta$ -- достаточно большая константа. Раскрыв скобки в предыдущем выражении, получаем + +$$ +Penalty_{i}^{(row)} + = + \tilde{w}_{i}^T + \beta + \underbrace{ + \left( 1_k - 2 I_k \right) + }_{\text{обозначим как } G} + \tilde{w}_{i} + + \text{const} +$$ + +Чтобы найти сумму $Penalty^{(row)} = \sum_i Penalty_{i}^{(row)}$, преобразуем бинарный вектор $\tilde{w}$ в другой бинарный вектор $\tilde{v}$, получающийся из $\tilde{w}$ определенной перестановкой элементов: + +$$ +\tilde{v} = [ + \tilde{w}_{11} \dots \tilde{w}_{1k} + \tilde{w}_{21} \dots \tilde{w}_{2k} + \dots + \tilde{w}_{N1} \dots \tilde{w}_{Nk} + ]^T +$$ + +Переход от $\tilde{w}$ к $\tilde{v}$ можно представить как линейное преобразование + +$$ +\tilde{v} = Q \tilde{w} +$$ + +с некоторой матрицей $Q \in \{0,1\}^{Nk \times Nk}$ (матрица $Q$ в свою очередь получается из единичной матрицы $I_{Nk}$ определенной перестановкой элементов). + +Сумма штрафов для строк равна + +$$ +Penalty^{(row)} = \sum_i Penalty_{i}^{(row)} = + \tilde{v}^T + \left( I_N \otimes \beta G \right) + \tilde{v} + = + \tilde{w}^T Q^T + \left( I_N \otimes \beta G \right) + Q \tilde{w} +$$ + +Соберем вместе квадратичную форму {eq}`eqn:sqr_form_without_constraints` (записанную без учета ограничений) и штрафы $Penalty^{(col)}$, $Penalty^{(row)}$ в одно финальное выражение, которое и является QUBO-формулировкой задачи о сбалансированной кластеризации методом k-средних: + +$$ +\tilde{w}^T + \left( + I_k \otimes (D + \alpha F) + + Q^T (I_N \otimes \beta G) Q + \right) +\tilde{w} + \to \min_{\tilde{w}} +$$ (eqn:final_qubo_kmeans) + + +### Оценка вычислительной сложности + +Задача {eq}`eqn:kmeans_equal_clusters` содержит $\mathcal{O}(N d)$ значений в данных и $\mathcal{O}(N)$ переменных. В QUBO-формулировке вводим по $k$ бинарных переменных для каждой исходной. Получается $\mathcal{O}(N k)$ переменных, а значит, требуется $\mathcal{O}(N^2 k^2)$ кубитов. + +Известно, что классический алгоритм сбалансированной кластеризации методом k-средних сходится за время $\mathcal{O}(N^{3.5} k^{3.5})$ в худшем случае (см. ссылки в статье {cite}`qubo3mlproblems2021`). Для оценки временнОй сложности QUBO рассматриваем три составляющие: + +1. Затраты времени для конвертации в QUBO-формулировку. В выражении {eq}`eqn:final_qubo_kmeans` по вычислительной сложности доминирует слагаемое, содержащее $I_k \otimes D$. Соответствующая вычислительная сложность $\mathcal{O}(N^2 k d)$. +2. Для реализации QUBO-задачи в квантовом "железе" потребуется $\mathcal{O}(N^2 k^2)$. +3. Время для выполнения квантового отжига и число повторений можно считать константами (см. комментарии к тому же пункту в обсуждении линейной регрессии). + +В итоге получаем полную вычислительную сложность $\mathcal{O}(N^2 k (d+k))$. Это лучше, чем результат классического алгоритма в худшем случае. Но количество итераций в классическом алгоритме сильно зависит от "удачности" начального приближения для центроидов кластеров. Классический алгоритм может оказаться и быстрее квантового. + + +## Заключение + +В этой лекции были рассмотрены три важные задачи машинного обучения, которые можно переформулировать в виде QUBO (*Quadratic Unconstrained Binary Optimization*) для решения на квантовом аннилере путем сведения к задаче нахождения основного состояния квантовой системы. Общий подход заключается в том, что минимизируемый в классической формулировке функционал переписывается в виде квадратичной формы относительно бинарных переменных, а вместо условий-ограничений в финальную квадратичную форму QUBO-задачи вводятся штрафы за нарушение этих ограничений. Есть надежда на то, что при некотором количестве кубитов квантовые алгоритмы будут иметь преимущество перед классическими по времени выполнения.