Алгоритм HHL

Автор(ы):

• Решетова Карина

• Токарев Игорь

Сегодня пришла пора поговорить о знаменитом алгоритме Харроу, Хассидима и Ллойда, более известном как HHL-алгоритме, способном решать системы линейных уравнений.

Очень надеюсь, что к данному занятию у вас уже есть представление об алгоритме фазовой оценки (QPE), использующем обратное квантовое преобразование Фурье, на котором и базируется HHL. Глубокое понимание всех тонкостей этого алгоритма потребует от вас уверенного владения математическим аппаратом. За детальным описанием вы всегда можете обратиться к статьям [DHM+18], [HHL09], и [HZL+17]. Приготовьтесь потратить время и умственные ресурсы, если алгоритм вас зацепит и вы решите в нем как следует покопаться. Мы же поможем вам заинтересоваться, рассмотрим основные принципы и небольшой пример.

Note

Именно HHL-алгоритм произвел настоящую революцию в области квантового машинного обучения. Ведь решение систем линейных уравнений так или иначе находится “под капотом” почти любого известного алгоритма машинного обучения. И действительно:

• классические линейная и логистическая регрессия сводятся именно к этой задаче;

• задача SVM может быть переформулирована в терминах решений систем линейных уравнений;

• задача нахождения обартной матрицы (часто используется в глубоком обучении) внутри обычно решается через решение линейной системы;

И это только малая часть примеров!

Так что знакомство с QML не будет полным без ознакомления с этим прекрасным, но очень сложным алгоритмом!

Fig. 64 Сет Ллойд, профессор MIT и один из создателей HHL-алгоритма

Задача

Представим обычную систему линейных уравнений:

\[\begin{split} \left\{\begin{array}{l} a_{11} x_{1}+a_{12} x_{2} = b_{1} \\ a_{21} x_{1}+a_{22} x_{2} = b_{2} \end{array}\right. \end{split}\]

Что в операторной форме можно переписать как:

\[ \large A\vec{x} = \vec{b}, \]

где \(A\) – эрмитова матрица.

Мы будем решать задачу на квантовом компьютере, то нам нужно перейти к квантовым состояниям:

\[ \large A|x\rangle = |b\rangle \]

Чтобы найти искомый вектор \(|x\rangle\), все что нам по сути нужно сделать – это найти обратный к \(A\) оператор (обозначаемый \(A^{-1}\)), который находится из равенства:

\[ \large AA^{-1} = A^{-1}A = I \]

Распишем применительно к нашей задаче поиска вектора \(|b\rangle\):

\[\begin{split} \large A|x\rangle = |b\rangle \\ \large A^{-1}A|x\rangle = A^{-1}|b\rangle \\ \large |x\rangle = A^{-1}|b\rangle \\ \end{split}\]

Оказывается, что и это все можно провернуть с помощью известных квантовых преобразований. Принципиальный вид нашей схемы представлен следующим образом:

Fig. 65 Квантовая схема, реализующая алгоритм HHL

В нижний регистр загружается вектор \(|b\rangle\), средний и нижний регистры участвуют в фазовой оценке, а верхний дополнительный кубит нужен для так называемого вращения, обусловленного собственными значениями. Давайте разбираться.

Реализация HHL

Для начала мы должны подготовить наши регистры по всем квантовым законам: \(|b\rangle\) и \(|x\rangle\) должны быть пронормированы, а оператор \(A\) должен быть эрмитовым. Надеемся, что Вы помните про ортонормированный базис, сферу Блоха, комплексное представление векторов \(|0\rangle\) и \(|1\rangle\)… Если нет, то обратитесь к предыдущим разделам курса.

Мы будем использовать оператор \(U = e^{iAt}\), и нужно, чтобы он был обратим –- для этого \(А\) должна быть эрмитовой.

Note

В случае, когда \(A\) не является эрмитовой, нужно перейти к эрмитовой матрице \(C\):

\[\begin{split} C = \begin{pmatrix} 0 & A\\ A^{\dagger} & 0 \end{pmatrix} \end{split}\]

И рассматривается задача \(C \vec{y} = \left(\begin{array}{l} \;\vec{b}\; \\ \;0\; \end{array}\right)\) для того, чтобы найти решение \(y = \left(\begin{array}{l} \;0\; \\ \;\vec{x}\; \end{array}\right)\)

Вспомним, что эрмитову матрицу \(A\) можно представить в виде суммы собственных векторов, умноженных на собственные значения, т.е. в виде спектрального разложения:

\[\begin{split} \begin{aligned} &A=\sum_{j=0}^{N-1} \lambda_{j}\left|u_{j}\right\rangle\left\langle u_{j}\right| \\ &A^{-1}=\sum_{j=0}^{N-1} \lambda_{j}^{-1}\left|u_{j}\right\rangle\left\langle u_{j}\right| \end{aligned} \end{split}\]

Тогда вектор\(|b\rangle\) можно представить через собственные векторы \(A\):

\[ |b\rangle = \sum_{j=0}^{N-1} b_{j} | u_j \rangle \]

Чтобы понять, почему это тоже ключевой момент, давайте вспомним, что значит собственный вектор и собственное значение матрицы.

Note

Собственным вектором \(|u\rangle\) оператора \(A\) называется такой ненулевой вектор, для которого выполняется:

\[\large A|u\rangle = \lambda |u\rangle\]

\(\lambda\) – собственное значение оператора \(A\).

Таким образом, искомый вектор \(|x\rangle\) – не что иное, как:

\[ |x\rangle = A^{-1} |b \rangle = \sum_{j=0}^{N-1} \lambda_{j}^{-1}b_j | u_{j}\rangle \]

Итак, фазовая оценка. Мы применяем к кубитам второго регистра матрицы Адамара, тем самым приводим их в суперпозицию. Следом запускаем оператор \(U\):

\[ U = e^{iAt} = \sum_{j=0}^{N-1} e^{i \lambda_j t} \left| u_j \rangle \langle u_j \right| \]

Для того, чтобы узнать собственное значение оператора \(U\), получения фазы (Quantum Phase Estimation – QPE), результатом которого получится следующее состояние:

\[ QPE(U, |0\rangle|u\rangle) = \]

\[ \frac{1}{2^{m/2}}( |0\rangle + e^{2\pi i 2^{m-1} \psi} |1\rangle) \otimes (|0\rangle + e^{2\pi i 2^{m-2} \psi} |1\rangle) \otimes \dots \otimes (|0\rangle + e^{2\pi i 2^{0} \psi} |1\rangle) \otimes |u\rangle = \]

\[ \frac{1}{2^{m/2}}\sum_{j=0}^{2^{m-1}} e^{2\pi i \psi j} |j\rangle |u\rangle = |\psi_{u}\rangle |u\rangle \]

Параметр \(t\) это нормировочная константа в случае \(U = e^{iAt}\):

\[ e^{2\pi i \psi} = e^{i \lambda_j t} \]

\[ \psi = \frac{\lambda_j t}{2\pi} \]

Параметр \(t\) подбирается с учетом того, что на выходе алгоритма QPE собственные значения \(\lambda_j\) нормализуются к виду \(0 \leq \lambda_j \leq 1\) и обычно мы располагаем ограниченным числом кубитов, которое можно использовать для аппроксимации.

Алгоритм обратного квантового Фурье переводит фазу в конкретный вектор.

Принципиальная схема QPE выглядит следующим образом:

Fig. 66 Схема алгоритма QPE

Итак, мы подготовились, вспомнили много хорошего, теперь пошагово распишем наш алгоритм.

Стартуем мы со следующим состоянием:

\[ \large |0\rangle_{a}|0\rangle_{r}|b\rangle_{m} \]

Т.е. наше состояние будет храниться в трех регистрах, в каждом из которых содержится столько кубитов, сколько нужно для решения задачи.

1. Применение QPE с использованием преобразования \(e^{iAt}\), после чего мы получим собственное значение оператора \(A\) во втором регистре:

   \[ |0\rangle_{a}|0\rangle_{r}|b\rangle_{m} \rightarrow \sum_{j=0}^{N-1}b_j|0\rangle_{a}|\lambda_j\rangle_r|u_j\rangle_m \]

2. Поворачиваем первый кубит (с индексом \(a\)), используя специальный оператор вращения \(R\):

   \[ R|0\rangle_{a} = \sum_{j=0}^{N-1}\left(\sqrt{1-\frac{C^{2}}{\lambda_{j}^{2}}}|0\rangle_{a} + \frac{C}{\lambda_{j}}|1\rangle_{a}\right), \]

   где \(C\) – константа, которая должна быть меньше минимального из лямбда: \(|C| < \lambda_{min}\) [Почему?].

   Переводим первый кубит \(|0\rangle_a\):

   \[ \sum_{j=0}^{N-1}b_j|0\rangle_{a}|\lambda_j\rangle_r|u_j\rangle_m \rightarrow \]
   \[ \sum_{j=0}^{N-1}\left(\sqrt{1-\frac{C^{2}}{\lambda_{j}^{2}}}|0\rangle+\frac{C}{\lambda_j}|1\rangle\right)b_j\left|\lambda_{j}\right\rangle_{n}\left|u_{j}\right\rangle_{m} \]

3. Применяем \(QPE^{\dagger}\) (т.е. обратное получение фазы) и получаем следующее состояние:

   \[ \sum_{j=0}^{N-1}\left(\sqrt{1-\frac{C^{2}}{\lambda_{j}^{2}}}|0\rangle+\frac{C}{\lambda_j}|1\rangle\right)b_{j}|0\rangle_{n}\left|u_{j}\right\rangle_{m} \]

   В конце мы измеряем верхний кубит и если получаем единицу, то знаем, что в нижнем регистре хранится искомый \(|x\rangle\) с учетом нормировки:

   \[ |x\rangle \approx \sum_{j=0}^{N-1}C(\frac{b_j}{\lambda_j})|u_j\rangle \]

Пример

Рассмотрим небольшой, но удобный пример. Удобный в том отношении, что, вообще говоря, алгоритм HHL имеет определенное приближение. Если собственные значения не представимы в бинарной форме, то о 100% точности говорить не приходится. Мы также опустим ряд вопросов, связанных с подбором параметров и количества кубитов второго регистра. Главное сейчас – понять, что происходит, и для этого наша матрица \(А\) эрмитова, все условия подобраны, а преобразования точны. Стоит помнить, что знать заранее значение собственных векторов и собственных значений нам совершенно не обязательно – это нужно лишь для наглядности.

Итак, пусть задача выглядит так :

\[\begin{split} \begin{aligned} &A=\left(\begin{array}{ll} 1 & \frac{3}{5} \\ \frac{3}{5} & 1 \end{array}\right) \\ &|b\rangle=|1\rangle=\left[\begin{array}{l} 0 \\ 1 \end{array}\right] \end{aligned} \end{split}\]

\[\begin{split} \left\{\begin{array}{l} x_{1}+\frac{3}{5} x_{2}=0 \\ \frac{3}{5} x_{1}+x_{2}=1 \end{array}\right. \end{split}\]

Собственные значения и соответствующие собственные векторы:

\[\begin{split} \large \begin{aligned} &\lambda_{0}=\frac{2}{5},\left|u_{0}\right\rangle=\left[\begin{array}{c} {-}1 \\ \;\;1 \end{array}\right] \\ &\lambda_{1}=\frac{8}{5},\left|u_{1}\right\rangle=\left[\begin{array}{l} \;\;1\; \\ \;\;1\; \end{array}\right] \end{aligned} \end{split}\]

Зададим параметр \(t\) и проанализируем фазу:

\[ \ t = 2\pi \frac{5}{16} \]

\[ e^{2\pi i \psi} = e^{i \lambda_j t} \]

\[ \psi = \frac{\lambda_j t}{2\pi} \]

\[ \frac{\lambda_0 t}{2\pi} = \frac{1}{8}\text{,} \;\;\;\;\;\;\;\; \frac{\lambda_1 t}{2\pi} = \frac{1}{2} \]

Как мы видим, для перевода угла \(\psi\) в векторную форму, нам понадобятся три кубита. После преобразования QPE мы имеем следующее состояние:

\[ \begin{aligned} &{QPE(|0\rangle_{a}|0\rangle_{r}|b\rangle_{m})=\sum_{j=0}^{1} \frac{1}{\sqrt{2}} \left |0\rangle_{a} |\lambda_{j} \rangle_{r} |u_{j}\right\rangle} = \frac{1}{\sqrt{2}}\left(|0\rangle_{a} |001\rangle_r \left|u_{0}\right\rangle + |0\rangle_{a} |100\rangle_r \left|u_{1}\right\rangle\right) \end{aligned} \]

Подберем константу \(C\) (как мы помним, она должна быть меньше наименьшего из собственных чисел) и произведем вращение:

\[\begin{split} \begin{aligned} &\frac{1}{\sqrt{2}} \left(\sqrt{1-\frac{(1 / 16)^{2}}{(1 / 8)^{2}}}|0\rangle+\frac{1 / 16}{1 / 8}|1\rangle\right)|001\rangle_r \left|u_{0}\right\rangle_m +\\ &+\frac{1}{\sqrt{2}}\left(\sqrt{1-\frac{(1 / 16)^{2}}{(1 / 2)^{2}}}|0\rangle_a+\frac{1 / 16}{1 / 2}|1\rangle_a\right)|100\rangle_r |u_{1}\rangle_m \end{aligned} \end{split}\]

В конце мы производим измерение верхнего кубита (с индексом \(a\)) и при получении единицы можем быть уверены, что нижний регистр содержит искомое решение с учетом нормировки:

\[\begin{split} \begin{aligned} \frac{1}{2\sqrt{2}}|1\rangle_a|000\rangle_r\left|u_{0}\right\rangle_m + \frac{1}{8\sqrt{2}}|1\rangle_a|000\rangle_r\left|u_{1}\right\rangle_m \end{aligned}\\ \end{split}\]

\[ |x\rangle \approx \frac{1}{2 \sqrt{2}}/(\sqrt{\frac{17}{128}})|u_0\rangle + \frac{1}{8 \sqrt{2}}/(\sqrt{\frac{17}{128}})|u_1\rangle \]