6 February 2017

Есть две функции

PythonAlgorithmsMathematics
Привет

Есть две булевы функции n аргументов, одна — константная, другая — сбалансированная. На какую сам сядешь, на какую фронтендера посадишь? Вот только функции неизвестны, а вызвать их разрешается лишь один раз.

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

I’ve come to talk with you again


Давайте чуть более формально поставим задачу о двух функциях. Пусть дана булева функция f: \{0, 1\}^n \to \{0, 1\} и о ней априорно известно, что она или константна, то есть для любого из своих аргументов всегда возвращает либо 0, либо 1, или сбалансирована, то есть ровно на половине своих аргументов возвращает 0, а ровно на половине 1. Требуется определить, константна функция или сбалансирована. При этом считается, что время вызова функции несоизмеримо больше любых других операций, поэтому сложность алгоритма определяется количеством вызовов функции.

hard_choice

Пример:

  1. f(x_1, x_2) = x_1 \oplus x_2 сбалансирована:
    x_1
    x_2
    f
    0 0 0
    0 1 1
    1 0 1
    1 1 0

  2. f(x_1, x_2) = x_1 \vee x_2 \vee (\neg x_1 \wedge \neg x_2) константна:
    x_1
    x_2
    f
    0 0 1
    0 1 1
    1 0 1
    1 1 1

  3. f(x_1, x_2) = x_1 \wedge x_2 ни сбалансирована, ни константна:
    x_1
    x_2
    f
    0 0 0
    0 1 0
    1 0 0
    1 1 1

Задача, безусловно, искусственная и навряд ли когда-нибудь кому-нибудь встретится на практике, однако она — классический проводник в дивный новый мир квантовых вычислений, а нарушать традиции я не осмелюсь.

Классическое детерминированное решение


bruteforce

Давайте, для начала, решим задачу в классической модели вычислений. Для этого в худшем случае понадобится вызвать функцию на 2^{n-1}+1 аргументах: ровно на половине и ещё одном. Если все вычисленные значения одинаковы, то функция, очевидно, константна. Если же существуют хотя-бы два различных результата, функция сбалансирована. Сложность детерминированного алгоритма экспоненциальна и составляет O(2^{n-1} + 1).

Алгоритм на Python:

from itertools import product, starmap, tee

def pairwise(xs):
    a, b = tee(xs)
    next(b, None)
    return zip(a, b)

def is_constant(f, n):
    m = 2 ** (n - 1)

    for i, (x, y) in enumerate(pairwise(starmap(f, product({0, 1}, repeat=n)))):
        if i > m:
            break

        if x != y:
            return False

    return True

Классическое вероятностное решение


rasta

А что, если вместо половины аргументов мы проверим меньшее их число и вынесем вердикт? Точным ответ уже не будет, но с какой вероятностью мы ошибёмся? Допустим, мы вычислили функцию на k < 2^{n-1} + 1 аргументах. Если среди значений функции есть два различных, то всё просто — функция сбалансирована. Иначе, мы объявляем её константной с вероятностью p(k) < 1. Предположим, что мы ошиблись и функция на самом деле сбалансирована. Посчитаем вероятность ошибки 1 - p(k). Если мы выбирали аргументы равномерно, то вероятность того, что два подряд идущих значения функции одинаковы, равна 1/2, а вероятность же встретить k одинаковых подряд идущих значений равна 1/2^k. Таким образом:

1 - p(k) = \frac{1}{2^k},

p(k) = 1 - \frac{1}{2^k}.

Обратная функция:

k(p) = \log_2{\frac{1}{1 - p}}.

При фиксированном p сложность классического вероятностного алгоритма константна и равна O(\log_2{\frac{1}{1 - p})}. Чтобы быть уверенным в ответе на 99.99%, необходимо вызвать функцию всего 14 раз.

Алгоритм на Python:

import random
from itertools import product, starmap, tee


def pairwise(xs):
    a, b = tee(xs)
    next(b, None)
    return zip(a, b)


def is_constant(f, n, k=14):
    xs = list(product({0, 1}, repeat=n))
    random.shuffle(xs)
    xs = xs[:k]

    for x, y in pairwise(starmap(f, xs)):
        if x != y:
            return False

    return True

А если я скажу вам, что существует константное детерминированное решение со сложностью O(1), позволяющее вызвать функцию всего один раз?

Правда, прежде чем его рассмотреть, придётся отвлечься…

Because a vision softly creeping


Мифы


Пережде чем начать, хотелось бы обсудить несколько расхожих мифов, связанных с квантовыми вычислениями:

  1. Квантовые алгоритмы — это сложно.

    difficult

    Да, их сложно синтезировать, ибо требуется для этого математическое воображение и проницательность. Их сложно реализовывать на настоящих квантовых компьютерах: для этого необходимо превосходно знать физику и засиживаться каждый день допоздна в лаборатории на кафедре. Но что точно не требует никаких особых знаний и невероятного количества усердия, так это их понимание. Я утверждаю, что каждый может понимать квантовые алгоритмы: они опираются на крайне простую математику, доступную любому первокурснику. Всё, что от вас требуется — лишь немного времени на изучение.

  2. Уже существуют квантовые компьютеры на тысячи кубитов от D-Wave
    Нет, это не настоящие квантовые компьютеры.

  3. Не существует ни одного настоящего квантового компьютера
    Нет, существуют. В лабораторных условиях и располагают они лишь несколькими кубитами.

  4. Квантовые компьютеры позволят решать задачи, недоступные ранее
    Нет, множество задач, вычислимых в классической и квантовой моделях, совпадают. Квантовые вычисления позволяют лишь снизить сложность небольшого подмножества этих задач.

  5. На квантовом компьютере Crysis на максималках будет летать

    wat

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

  6. Квантовый компьютер — это чёрная коробочка с входом и выходом, заглянув в которую можно всё испортить

    blackbox

    Если вам 12 лет, такая аналогия сгодится. В любом другом случае она, как и все другие аналогии с коробками, котами, лупами и «связанными» нитками электронами, так активно продвигаемые во всех научно-популярных источниках, только лишь запутывает, создаёт иллюзию ложного понимания и скорее вредна, чем полезна. Откажитесь от этих аналогий.

Зачем?


Зачем прикладному математику (программисту) уметь разбираться в квантовых алгоритмах на прикладном уровне? Тут всё просто, я готов предложить вам два довода:

  1. Для саморазвития. Почему бы и нет?
  2. Они уже идут. Они, квантовые компьютеры. Они уже рядом. Моргнуть не успеете, как парочка появится в серверной вашей компании, а ещё через несколько лет и в виде сопроцессора в вашем ноутбуке. И бежать уже будет некуда. Придётся под них программировать, вызывать квантовые сопрограммы. А без понимания это делать сложно, согласитесь.

Left its seeds while I was sleeping


Самой базовой составляющей квантовых вычислений является квантовая система. Квантовая система — физическая система, все действия которой сравнимы по величине с постоянной Планка. Это определение и тот факт, что квантовые системы подчиняются законам матричной механики — все знания, которые нам потребуются из физики. Дальше — только математика.

loldontinterrupt

Как и любая другая физическая система, квантовая система может находиться в определённом состоянии. Все возможные состояния квантовой системы образуют гильбертово пространство \mathcal{H} над полем комплексных чисел. Надеюсь, с концепцией комплексных чисел читатель знаком — их понимание необходимо везде в дальнейшем. Если нет, советую познакомиться и возвращаться. Гильбертово пространство же является полным нормированным метрическим линейным пространством с нормой \Vert x \Vert = \sqrt{(x, x)}, где (x, y) — скалярное произведение. По порядку с конца:

  1. Линейное (векторное) пространство — множество элементов X с введёнными на нём операциями сложения элементов x + y и умножения x \cdot \lambda на элемент поля K (в нашем случае поля комплексных чисел). Операции эти должны быть замкнуты (результат должен принадлежать множеству X) и должны выполняться 8 аксиом. Посмотреть полный их список, а так же познакомиться подробнее с линейными пространствами рекомендую здесь.

  2. В метрическом пространстве X для любых элементов x, y \in X определено расстояние \rho(x, y), которое удовлетворяет требованиям (аксиомам метрического
    пространства):

    • \rho(x, y) \geqslant 0, при этом \rho(x, y) = 0 тогда и только тогда, когда x и y совпадают;
    • \rho(x, y) = \rho(y, x);
    • \rho(x, y) \leqslant \rho(x, z) + \rho(z, y) — неравенство треугольника.

  3. В нормированном пространстве X для любого элемента x \in X существует вещественное число \Vert x \Vert \in R, называемое его нормой и удовлетворяюещее, опять же, трём аксиомам:

    • \Vert x \Vert \geqslant 0, если же \Vert x \Vert = 0, то x = 0 — нулевой элемент;
    • \Vert \lambda \cdot x \Vert = \vert \lambda \vert \cdot \Vert x \Vert;
    • \Vert x + y \Vert \leqslant \Vert x \Vert + \Vert y \Vert.


Введём в полученном пространстве скалярное произведение, удовлетворяющее обычным требованиям скалярного произведения, введём норму как показано выше и получим пространство Гильберта.

hilbertalmost

Обсудим ещё понятие сопряжённого пространства. Пространством, сопряжённым к \mathcal{H} называется пространство \mathcal{H^*} линейных операторов над \mathcal{H}. Что такое линейный оператор? Можно думать о нём как об обобщении линейной функции: для линейноного оператора A: \mathcal{H} \to Y должно выполняться:

A (\lambda_1 x_1 + \lambda_2 x_2) = \lambda_1 A x_1 + \lambda_2 A x_2,

где x_1, x_2 \in \mathcal{H}, \lambda_1, \lambda_2 \in K. (На самом деле, ещё и его норма должна быть ограничена на единичной гиперсфере, но, дабы избежать ещё десятка громоздких определений, ограничимся таким интуитивным представлением.)

В квантовой информатике по историческим причинам используются обозначения Дирака. Они могут казаться необоснованно громоздкими и вычурными, но они — стандарт, которого стоит придерживаться. В этих обозначениях элемент нашего гильбертова пространства, описывающий состояние системы, называется кет-вектором и обозначается

\left| \psi \right\rangle \in \mathcal{H}.

Бра-вектором называется элемент сопряжённого пространства

\langle\left \phi \right| \in \mathcal{H^*},

такой что

(\langle\left \phi \right|) \left| \psi \right\rangle = (\left| \phi \right\rangle, \left| \psi \right\rangle) = \langle\left \phi | \psi \right\rangle.

То есть, это линейный оператор, применение которого к нашему вектору состояния аналогично скалярному произведению на соответствующий ему элемент «оригинального» гильбертова пространства. Для удобства записи при приминении бра-вектора к кет-вектору две вертикальные черты сливаются в одну, что показано в выражении выше.

Важно то, что вектора, отличающиеся лишь умножением на некоторую ненулевую константу, отвечают одному и тому же физическому состоянию, поэтому часто рассматривают не все возможные состояния, а лишь нормированные, то есть такое подмножество \mathcal{H}, что

\langle\left \psi | \psi \right\rangle = 1

— норма каждого элемента равна единице. Все такие вектора живут на единичной гиперсфере.

Если мы выделим в нашем гильбертовом пространстве состояний некоторый базис \{\left| e_i \right\rangle\}_{i=1}^m \subset \mathcal{H}, то любой кет-вектор сможем записать в виде вектора комплексных чисел — коэффициентов его разложения по этому базису:

\left| \psi \right\rangle = \sum_{i=1}^m{\langle\left e_i | \psi \right\rangle \left| e_i \right\rangle},

при этом матричная механика говорит нам, что квадраты модулей коэффициентов разложения \vert \langle\left e_i | \psi \right\rangle \vert^2 физически означают вероятности нахождения квантовой системы в соответствующем базисном состоянии при измерении в этом базисе.

Вот оно — первое и основное свойство квантовых систем, которое так часто мусолят в популярных статьях: если измерить систему в некотором базисе, то она перейдёт в одно из базисных состояний, потеряет информацию и не сможет вернуться назад. Вот только при прочтении складывается ощущение, что всё происходит абсолютно случайно и повлиять на это никак нельзя, тогда как на самом деле вероятности перехода известны заранее и, более того, зависят от измерительного базиса. Если бы всё было так случайно, как нам представляют, детерминированные квантовые алгоритмы были бы невозможны.

Если мы можем представить элемент Гильбертова пространства вектором при некотором фиксированном базисе, то линейный оператор над этим пространством мы можем представить матрицей.

Действительно,

A \left| \psi \right\rangle = \sum_{i=1}^m{\langle\left e_i | \psi \right\rangle A \left| e_i \right\rangle}

эквивалентно

A \left| \psi \right\rangle = \hat{A} \hat{\psi},

где \hat{A} получается поочерёдным применением оператора A к базисным элементам \left| e_i \right\rangle и выписыванием получившихся элементов по строкам, а \hat{\psi} — разложение \left| \psi \right\rangle в этом же базисе.

Пусть оператор A представляется матрицей

A = \begin{pmatrix}
    a_{11} & a_{12} & a_{13} & \dots  & a_{1m} \\
    a_{21} & a_{22} & a_{23} & \dots  & a_{2m} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    a_{m1} & a_{m2} & a_{m3} & \dots  & a_{mm}
\end{pmatrix}.

Элементы матрицы — комплексные числа. Давайте возьмём каждый элемент и заменим его комплексно сопряжённым (комплексно сопряжённым к x + iy называется число x - iy) и, заодно, транспонируем всю матрицу:

A^\dagger = \overline{A}^T = \begin{pmatrix}
    \overline{a_{11}} & \overline{a_{21}} & \overline{a_{31}} & \dots  & \overline{a_{m1}} \\
    \overline{a_{12}} & \overline{a_{22}} & \overline{a_{32}} & \dots  & \overline{a_{m2}} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    \overline{a_{1m}} & \overline{a_{2m}} & \overline{a_{3m}} & \dots  & \overline{a_{mm}}
\end{pmatrix}.

Такая матрица A^\dagger называется эрмитово-сопряжённой к A. Если A A^\dagger = A^\dagger A = I, где I — единичный оператор (с единичной матрицей), то соответствующий оператор называется унитарным.

Второе правило, которое диктует нам матричная механика: на квантовую систему могут действовать только унитарные операторы. Почему? Потому что такие преобразования обратимы во времени и не теряют информацию. Действительно, если

U \left| \psi_0 \right\rangle = \left| \psi_1 \right\rangle,

то можно применить обратное преобразование

U^\dagger \left| \psi_1 \right\rangle = U^\dagger U \left| \psi_0 \right\rangle = \left| \psi_0 \right\rangle

и получить исходное состояние системы.

Напоследок самое важное: тензорное произведение. Тензорным произведением двух гильбертовых пространств \mathcal{H}_1 и \mathcal{H}_2 называется гильбертово пространство, обозначаемое \mathcal{H}_1 \otimes \mathcal{H}_2. Я не буду приводить формальное определение, лишь отмечу важные нам свойства:

  1. Размерность получившегося пространства равна произведению размерностей исходных пространств:
    \dim{\mathcal{H}_1 \otimes \mathcal{H}_2$} = \dim{\mathcal{H}_1} \cdot \dim{\mathcal{H}_2};
  2. Если \{\left| e_i \right\rangle\}_{i=1}^m — базис \mathcal{H}_1, а \{\left| f_i \right\rangle\}_{i=1}^n — базис \mathcal{H}_2, то \{\left| e_i \otimes f_j \right\rangle\}_{i=1, j=1}^{m, n} — порождающий базис \mathcal{H}_1 \otimes \mathcal{H}_2.

Тензорным же произведенем операторов A из \mathcal{H}_1^* и B из \mathcal{H}_2^* (оператор A представлен матрицей, показанной выше) называется оператор A \otimes B из \big[ \mathcal{H}_1 \otimes \mathcal{H}_2 \big]^* с матрицей

A \otimes B = \begin{pmatrix}
    a_{11} \cdot B & a_{12} \cdot B & a_{13} \cdot B & \dots  & a_{1m} \cdot B \\
    a_{21} \cdot B & a_{22} \cdot B & a_{23} \cdot B & \dots  & a_{2m} \cdot B \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    a_{m1} \cdot B & a_{m2} \cdot B & a_{m3} \cdot B & \dots  & a_{mm} \cdot B
\end{pmatrix}.

Такое произведение ещё называется произведением Кронекера: мы умножаем вторую матрицу на каждый элемент первой матрицы и из получившихся блоков составляем блочную матрицу. Если размерность A равнялась n_1 \times n_2, а размерность B равнялась m_1 \times m_2, то размерность матрицы, полученной в ходе тензорного произведения, будет равна n_1 \cdot m_1 \times n_2 \cdot m_2.

Пример:

A = \begin{bmatrix}
    1 & 2 \\
    3 & 0
\end{bmatrix},\;
B = \begin{bmatrix}
    1 & 1 \\
    1 & 2
\end{bmatrix},

A \otimes B = \begin{bmatrix}
    1 \cdot \begin{pmatrix}
    1 & 1 \\
    1 & 2
\end{pmatrix} & 2 \cdot \begin{pmatrix}
    1 & 1 \\
    1 & 2
\end{pmatrix} \\
    3 \cdot \begin{pmatrix}
    1 & 1 \\
    1 & 2
\end{pmatrix} & 0 \cdot \begin{pmatrix}
    1 & 1 \\
    1 & 2
\end{pmatrix}
\end{bmatrix} = \begin{bmatrix}
    1 & 1 & 2 & 2 \\
    1 & 2 & 2 & 4 \\
    3 & 3 & 0 & 0 \\
    3 & 6 & 0 & 0 \\
\end{bmatrix}.


A = \begin{bmatrix}
    1 \\
    0 \\
\end{bmatrix},\;
B = \begin{bmatrix}
    1 \\
    2
\end{bmatrix},

A \otimes B = \begin{bmatrix}
    1 \\
    2 \\
    0 \\
    0 \\
\end{bmatrix}.

Третье важное свойство квантовых систем: две квантовые системы могут находиться в состоянии суперпозиции, при этом новое пространство состояний представляет собой тензорное произведение исходных пространств, а состояние новой системы будет являться тензорным произведением состояний исходных систем. Так, суперпозицией систем в состояниях \left| \psi \right\rangle \in \mathcal{H}_1 и \left| \phi \right\rangle \in \mathcal{H}_2 будет новая система в состоянии \left| \psi \right\rangle \otimes \left| \phi \right\rangle, описываемая гильбертовым пространством \mathcal{H}_1 \otimes \mathcal{H}_2.

And the vision that was planted in my brain


Вот и вся математика, что нам понадобится. На всякий случай, резюмируя:

  1. При фиксированном базисе квантовую систему можно описать комплексным вектором, а эволюцию этой системы — унитарной комплексной матрицей;
  2. Квантовую систему можно измерить в каком-либо базисе и она перейдёт в одно из базисных состояний в соответствии с заранее определёнными вероятностями.

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

woooho

Кубит


Давайте рассмотрим некоторую квантовую систему, описываемую двумерным гильбертовым пространством \mathcal{H}^2 и выделим в ней некоторый базис, вектора которого обозначим как \left| 0 \right\rangle и \left| 1 \right\rangle. Внутри скобок записывается индекс базисного вектора в двоичной системе счисления, начиная с нуля без дополнительных символов. Такие обозначения окажутся очень удобными. Таким образом,

\left| 0 \right\rangle = \begin{bmatrix}
    1 \\
    0 
\end{bmatrix},\;
\left| 1 \right\rangle = \begin{bmatrix}
    0 \\
    1 
\end{bmatrix},

и произвольный вектор \left| \psi \right\rangle \in \mathcal{H}^2 можно выразить следующим образом:

\left| \psi \right\rangle = \alpha \left| 0 \right\rangle + \beta \left| 1 \right\rangle,

где \alpha и \beta — некоторые комплексные числа, такие что \vert \alpha \vert^2 + \vert \beta \vert^2 = 1 (вспомните интерпретацию коэффициентов разложения и условие нормировки из предыдущего параграфа). Так вот, такая простая квантовая система называется кубитом (quantum bit, qbit). Кубит — аналог классического бита в квантовой модели вычислений. Нетрудно видеть, что пространство всевозможных состояний одного кубита \mathcal{H}^2 представляет из себя трёхмерную сферу, называемую сферой Блоха, где \left| 0 \right\rangle соответствует нижнему полюсу, а \left| 1 \right\rangle — верхнему.

sphere


Регистр


Один кубит, как и один бит — слишком скучно, поэтому сразу рассмотрим суперпозицию нескольких кубитов. Такая суперпозиция называется квантовым регистром (quantum register, qregister) из n кубитов. Например, квантовый регистр из 2 кубитов будет описываться пространством \mathcal{H}^4 и иметь 4 базисных состояния:

\left| 00 \right\rangle = \left| 0 \right\rangle \otimes \left| 0 \right\rangle = \begin{bmatrix}
    1 \\
    0 \\
    0 \\
    0
\end{bmatrix},

\left| 01 \right\rangle = \left| 0 \right\rangle \otimes \left| 1 \right\rangle = \begin{bmatrix}
    0 \\
    1 \\
    0 \\
    0
\end{bmatrix},

\left| 10 \right\rangle = \left| 1 \right\rangle \otimes \left| 0 \right\rangle = \begin{bmatrix}
    0 \\
    0 \\
    1 \\
    0
\end{bmatrix},

\left| 11 \right\rangle = \left| 1 \right\rangle \otimes \left| 1 \right\rangle = \begin{bmatrix}
    0 \\
    0 \\
    0 \\
    1
\end{bmatrix}.

Соответственно, любое состояние такого регистра \left| \phi \right\rangle \in \mathcal{H}^4 можно представить в виде

\left| \phi \right\rangle = \alpha_1 \left| 00 \right\rangle + \alpha_2 \left| 01 \right\rangle + \alpha_3 \left| 10 \right\rangle + \alpha_4 \left| 11 \right\rangle,

где \vert  \alpha_1 \vert^2 + \vert  \alpha_2 \vert^2 + \vert  \alpha_3 \vert^2 + \vert  \alpha_4 \vert^2 = 1. Заметьте, что в наших обозначениях базисный вектор с единицей на k-ом месте обозначется числом k, записанным в двоичном виде.

Далее аналогично. Квантовый регистр из m кубитов будет описываться 2^m-мерным гильбертовым пространством \mathcal{H}^{2^m}, иметь 2^m базисных состояний, формируемых аналогичным образом. Не откладывая надолго, научимся эмулировать квантовые регистры:

import numpy as np


class QRegister:
    def __init__(self, n_qbits, init):
        self._n = n_qbits
        assert len(init) == self._n

        self._data = np.zeros((2 ** self._n), dtype=np.complex64)
        self._data[int('0b' + init, 2)] = 1

3 строчки кода для создания квантового регистра — совсем не сложно, согласитесь. Использовать можно таким образом:

a = QRegister(1, '0')    # |0>
b = QRegister(1, '1')    # |1>
c = QRegister(3, '010')  # |010>

Квантовый алгоритм включает в себя:

  1. Инициализацию квантового регистра;
  2. Набор унитарных преобразований над ним;
  3. Измерение результата.

Измерение


С первым пунктом разобрались и научились его эмулировать, давайте теперь научимся эмулировать последний: измерение. Как вы помните, квадраты коэффициентов вектора состояния физически означают вероятности перехода в это состояние. В соответствии с этим реализуем новый метод в классе QRegister:

def measure(self):
    probs = np.real(self._data) ** 2 + np.imag(self._data) ** 2
    states = np.arange(2 ** self._n)
    mstate = np.random.choice(states, size=1, p=probs)[0]
    return f'{mstate:>0{self._n}b}'

Генерируем вероятности probs выбора одного из 2^n состояний states и случайно выбираем его с помощью np.random.choice. Остаётся только вернуть бинарную строку с соответствующим количеством дополняющих нулей. Очевидно, что для базисных состояний ответ всегда будет одинаков и равен самому этому состоянию. Проверим:

>>> QRegister(1, '0').measure()
'0'
>>> QRegister(2, '10').measure()
'10'
>>> QRegister(8, '01001101').measure()
'01001101'

Почти всё готово, чтобы решить нашу задачу! Осталось лишь научиться влиять на квантовые регистры. Мы уже знаем, что сделать это можно унитарными преобразованиями. В квантовой информатике унитарное преобразование называется гейтом (quantum gate, qgate, gate).

Гейты


В рамках данной статьи рассмотрим лишь малое число самых основных гейтов, которые нам пригодятся. На самом деле, их намного больше.

Единичный


Единичный гейт — самый простой, какой только можно рассмотреть. Его матрица выглядит следующим образом:

I = \begin{bmatrix}
    1 & 0\\
    0 & 1
\end{bmatrix}.

Он никак не изменяет кубит, на который действует:

I \left| 0 \right\rangle = \left| 0 \right\rangle,\; I \left| 1 \right\rangle = \left| 1 \right\rangle, \; I \left| \psi \right\rangle = \left| \psi \right\rangle,

однако не стоит считать его бесполезным — он нам понадобится, и не раз.

Гейт Адамара


H = \frac{1}{\sqrt{2}}
\begin{bmatrix}
1 & 1\\
1 & -1
\end{bmatrix}

Не составляет труда проверить, что матрица унитарна:

H H^\dagger = \frac{1}{\sqrt{2}}
\begin{bmatrix}
1 & 1\\
1 & -1
\end{bmatrix} \cdot \frac{1}{\sqrt{2}}
\begin{bmatrix}
1 & 1\\
1 & -1
\end{bmatrix} = \frac{1}{2} \begin{bmatrix}
2 & 0\\
0 & 2
\end{bmatrix} = I.

Рассмотрим действие гейта Адамара на базисные кубиты:

H \left| 0 \right\rangle = \frac{1}{\sqrt{2}}
\begin{bmatrix}
1 & 1\\
1 & -1
\end{bmatrix}
\begin{bmatrix}
1\\
0
\end{bmatrix} = \frac{1}{\sqrt{2}}
\begin{bmatrix}
1\\
1
\end{bmatrix} = \frac{1}{\sqrt{2}} (\left| 0 \right\rangle + \left| 1 \right\rangle),

H \left| 1 \right\rangle = \frac{1}{\sqrt{2}}
\begin{bmatrix}
1 & 1\\
1 & -1
\end{bmatrix}
\begin{bmatrix}
0\\
1
\end{bmatrix} = \frac{1}{\sqrt{2}}
\begin{bmatrix}
1\\
-1
\end{bmatrix} = \frac{1}{\sqrt{2}} (\left| 0 \right\rangle - \left| 1 \right\rangle).

Или в общем виде [1]:

H \left| x \right\rangle = \frac{1}{\sqrt{2}} \sum_{y \in \lbrace 0,1 \rbrace} (-1)^{x \cdot y} \left| y \right\rangle.

Как можно заметить, гейт Адамара любое базисное состояние переводит в равновероятное — при измерении с равной вероятностью можно получить любой результат.

Гейты Паули


Три крайне важных гейта, которым соответствуют матрицы, введённые Вольфгангом Паули:

X =
\begin{bmatrix}
0 & 1\\
1 & 0
\end{bmatrix},\;
Y =
\begin{bmatrix}
0 & -i\\
i & 0
\end{bmatrix},\;
Z =
\begin{bmatrix}
1 & 0\\
0 & -1
\end{bmatrix}.

Гейт X ещё называется NOT-гейтом:

X \left| 0 \right\rangle = \left| 1 \right\rangle, \; X \left| 1 \right\rangle = \left| 0 \right\rangle,

а геометрически его применение эквивалентно повороту на сфере Блоха на \pi радиан вокруг оси x.

X gate

Гейты Y и Z аналогичны гейту X за тем исключением, что вращение производится вокруг соответствующей оси.

Существует теорема, согласно которой с помощью гейтов I, X, Y и Z можно выразить любой однокубитный гейт. Например:

H = \frac{1}{\sqrt{2}}\big(X + Z\big),

откуда видно, что гейт Адамара геометрически означает вращение на \pi радиан вокруг оси \frac{x + z}{\sqrt{2}}.

Реализуем все рассмотренные гейты на Python. Для этого создадим ещё один класс:

class QGate:
    def __init__(self, matrix):
        self._data = np.array(matrix, dtype=np.complex64)

        assert len(self._data.shape) == 2
        assert self._data.shape[0] == self._data.shape[1]

        self._n = np.log2(self._data.shape[0])

        assert self._n.is_integer()

        self._n = int(self._n)


А в класс QRegister добавим операцию применения гейта:

def apply(self, gate):
    assert isinstance(gate, QGate)
    assert self._n == gate._n
    self._data = gate._data @ self._data


И создадим уже известные нам гейты:

I = QGate([[1, 0], [0, 1]])
H = QGate(np.array([[1, 1], [1, -1]]) / np.sqrt(2))
X = QGate([[0, 1], [1, 0]])
Y = QGate([[0, -1j], [1j, 0]])
Z = QGate([[1, 0], [0, -1]])

Орёл или решка?



coin

Давайте для примера рассмотрим простейший квантовый алгоритм: он будет генерировать случайный бит — ноль или один, орла или решку. Это будет самая честная монета во вселенной — результат станет известен лишь при измерении, а природа случайности зашита в само основание вселенной и повлиять на неё невозможно никоим образом.

Для алгоримта нам понадобится всего один кубит. Пусть в начальный момент времени он находится в состоянии \left| 0 \right\rangle:

\left| 0 \right\rangle = \begin{bmatrix}
1\\
0
\end{bmatrix}.

Теперь применим к нему гейт Адамара H чтобы получить состояние

H \left| 0 \right\rangle = \frac{1}{\sqrt{2}}
\begin{bmatrix}
1\\
1
\end{bmatrix}.

Если мы сейчас измерим полученную систему, с вероятностью \Big\vert \frac{1}{\sqrt{2}} \Big\vert^2 = \frac{1}{2} она окажется в состоянии \left| 0 \right\rangle и с точно такой же вероятностью в состоянии \left| 1 \right\rangle. Останется только записать полученный результат.

Проверим алгоритм, эмулируя его на нашем классическом компьютере:

from quantum import QRegister, H


def quantum_randbit():
    a = QRegister(1, '0')
    a.apply(H)

    return a.measure()


for i in range(32):
    print(quantum_randbit(), end='')
print()

Результаты:

➜ python example-randbit.py
11110011101010111010011100000111

➜ python example-randbit.py
01110000111100011000101010100011

➜ python example-randbit.py
11111110011000001101010000100000

Весь приведённый выше алгоритм может быть записан парой формул:

\left| q_0 \right\rangle = \left| 0 \right\rangle\\
\left| q_1 \right\rangle = H \left| q_0 \right\rangle\\
r = measure(\left| q_1 \right\rangle)

Однако, с такой записью не очень удобно работать: структура списка — последовательности действий, хорошо подходящая для классических алгоритмов, не применима в квантовом случае: здесь у нас нет ни циклов, ни условий, лишь течение состояния вперёд (а иногда и назад) во времени. Поэтому для описания алгоритмов в квантовой информатике широко используются квантовые схемы. Вот схема приведённого выше алгоритма:
simple

Слева всегда обозначено исходное состояние системы. В прямоугольниках указываются унитарные преобразования, производимые над этим состоянием, а в конце на всех или на нескольких кубитах распологается значок измерительного прибора — операция измерения. Существует ещё «синтаксический сахар» для некоторых многокубитных преобразований в виде точек, разветвлений и кружочков. И всё. Если вы умеете отличать квадрат от треугольника и круга, вы без проблем поймёте любую схему квантового алгоритма.

holes

Больше кубитов богу кубитов


А что, если мы работаем не с одним кубитом, а с целым их регистром? И, допустим, хотим применить гейт лишь к одному кубиту? На помощь приходят свойства тензорного произведения. По определению тензорного произведения оператора A: V \to W и B: X \to Y:

(A \otimes B)(v \otimes x) = Av \otimes Bx,

где v \in V, x \in X. Тогда:

(A \otimes I)(v \otimes x) = Av \otimes Ix = Av \otimes x,

то есть, всё равно — применить ли гейт к одному кубиту, а затем связать его со вторым и получить квантовый регистр или же применить ко всему регистру оператор A \otimes I. Аналогично, если мы хотим применить оператор A только ко второму кубиту, мы можем применить ко всему регистру оператор I \otimes A. Это значит, что такая схема:

example0

Полностью аналогична такой:

example1

Только единичные гейты для удобства опускают.

А если мы, наоборот, хотим применить гейт сразу к нескольким кубитам? Опять же, из определения тензорного произведения, для этого мы можем применить к ним этот гейт, тензорно умноженный сам на себя необходимое количество раз:

Hv \otimes Hx = (H \otimes H) (v \otimes x) = H^{\otimes2} (v \otimes x).

A^{\otimes n} означает тензорное возведение в степень. Кстати, \left| \psi \right \rangle^{\otimes n} для краткости записывают как \left| \psi^n \right \rangle. Таким образом,

H^{\otimes n} \left| 0^n \right \rangle

означает: к каждому кубиту квантового регистра из n нулевых кубитов применить преобразование Адамара.

Добавим тензорное произведение и возведение в степень в наш QGate:

def __matmul__(self, other):
    return QGate(np.kron(self._data, other._data))

def __pow__(self, n, modulo=None):
    x = self._data.copy()

    for _ in range(n - 1):
        x = np.kron(x, self._data)

    return QGate(x)

Квантовый оракул


Каждой бинарной функции f : \{0, 1\}^n \to \{0, 1\} соответствует единственный многокубитный гейт размерности n + 1, который обозначается U_f и называется квантовым оракулом для этой функции. Он содержит в себе всю информацию о функции f и «позволяет» одновременно вызвать её на всех её возможных аргументах.

Почему его размерность n + 1? Всё дело в том, что согласно ещё одному фундаментальному свойству квантовых вычислений, в них не теряется информация и они обратимы во времени. Если мы в классической модели вычислений вызовем функцию f(x, y) = x \oplus y и получим результат, то по нему одному мы не сможем сказать, с какими аргументами функция была вызвана и, соответственно, не сможем обратить вычисления во времени.

diagram0


Один из способов избежать этого — запомнить аргументы, с которыми была вызвана функция.

diagram1


То же самое происходит и в квантовой модели вычислений, только там оно встроено в саму их природу — без сохранения полной информации невозможно построить унитарное преобразование.

Квантовый оракул U_f определяется как унитарное преобразование, переводящее состояние \left| x \right \rangle \left| y \right \rangle в состояние \left| x \right \rangle \left| y \oplus f(x) \right \rangle, где n кубитов, обозначенных x, несут в себе информацию об аргументе функции и остаются неизменными, а единственный кубит y — её результат. \oplus обозначает сложение по модулю 2.

Разберёмся на примере. Допустим, мы хотим построить оракул для функции одного аргумента f(x) = x. Значит, соответствующий оператор U_f будет двухкубитным и может быть описан квадратной матрицей размерности 2^2 = 4. В соответствии с описанным выше правилом посмотрим, в какие состояния должны перейти все возможные базисные состояния регистра:

\left| 00 \right \rangle  = \left| 0 \right \rangle \left| 0 \right \rangle \to \left| 0 \right \rangle \left| 0 \oplus f(0) \right \rangle = \left| 00 \right \rangle

\left| 01 \right \rangle  = \left| 0 \right \rangle \left| 1 \right \rangle \to \left| 0 \right \rangle \left| 0 \oplus f(1) \right \rangle = \left| 01 \right \rangle

\left| 10 \right \rangle  = \left| 1 \right \rangle \left| 0 \right \rangle \to \left| 1 \right \rangle \left| 1 \oplus f(0) \right \rangle = \left| 11 \right \rangle

\left| 11 \right \rangle  = \left| 1 \right \rangle \left| 1 \right \rangle \to \left| 1 \right \rangle \left| 1 \oplus f(1) \right \rangle = \left| 10 \right \rangle

Оператор с единичной матрицей переводит любое состояние в себя же, так как при умножении матрицы на i-ую строку мы берём i-ый же компонент вектора и ставим его на то же самое место. Если в матрице на i-ой строке все элементы кроме j-го будут нулями, а он единицей, то на j-ое место результирующего вектора мы поставим элемент, который стоял на i-ом месте исходного вектора. В примере выше нам нужно оставить на месте {00}_2 = 0 и {01}_2 = 1 элементы вектора состояния и поменять местами {10}_2 = 3 и {11}_2 = 4 элементы. Тогда U_f будет выглядеть следующим образом:

U_f = \begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & 1 & 0\\
\end{bmatrix}.

Проверим:

\begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & 1 & 0\\
\end{bmatrix}\begin{bmatrix}
x_{00} \\
x_{01} \\
x_{10} \\
x_{11} \\
\end{bmatrix} = \begin{bmatrix}
x_{00} \\
x_{01} \\
x_{11} \\
x_{10} \\
\end{bmatrix}.

Аналогично для любой другой функции произвольного числа аргументов. Реализуем оракл:

def U(f, n):
    m = n + 1

    U = np.zeros((2**m, 2**m), dtype=np.complex64)

    def bin2int(xs):
        r = 0
        for i, x in enumerate(reversed(xs)):
            r += x * 2 ** i
        return r

    for xs in product({0, 1}, repeat=m):
        x = xs[:~0]
        y = xs[~0]

        z = y ^ f(*x)

        instate = bin2int(xs)
        outstate = bin2int(list(x) + [z])
        U[instate, outstate] = 1

    return QGate(U)

Still remains


Ну вот мы и рассмотрели весь базис, который необходим для решения поставленной задачи и даже успели написать небольшой фреймворк на питоне, который нам поможет это решение проверить. Алгоритм, позволяющий за один вызов функции определить, константна функция или сбалансирована, именуется алгоритмом Дойча — Йожи. Версия алгоритма для функции одной переменной была разработанна в 1985 году Дэвидом Дойчем, а затем обобщёна на случай нескольких переменных в 1992 году с помощью Ричарда Йожи. Вот схема этого алгоритма:

dj

Если результатом измерения станет 0, то функция константна, иначе — сбалансирована. Сразу реализуем алгоритм:

from quantum import QRegister, H, I, U


def is_constant(f, n):
    q = QRegister(n + 1, '0' * n + '1')
    q.apply(H ** (n + 1))
    q.apply(U(f, n))
    q.apply(H ** n @ I)

    if q.measure()[:~0] == '0' * n:
        return True
    else:
        return False

И проверим:

def f1(x):
    return x


def f2(x):
    return 1


def f3(x, y):
    return x ^ y


def f4(x, y, z):
    return 0


print('f(x) = x is {}'.format('constant' if is_constant(f1, 1) else 'balansed'))
print('f(x) = 1 is {}'.format('constant' if is_constant(f2, 1) else 'balansed'))
print('f(x, y) = x ^ y is {}'.format('constant' if is_constant(f3, 2) else 'balansed'))
print('f(x, y, z) = 0 is {}'.format('constant' if is_constant(f4, 3) else 'balansed'))

Результат:

f(x) = x is balansed
f(x) = 1 is constant
f(x, y) = x ^ y is balansed
f(x, y, z) = 0 is constant

Работает. Замечательно. А почему он работает? Докажем корректность алгоритма и заодно посмотрим на принцип его работы.

Рассмотрим действие гейта H на квантовый регистр в произвольном базисном состоянии [1]:

H^{\otimes n} \left|x_1, x_2, \ldots, x_n\right\rangle =

(по свойству тензорного произведения)

= H \left| x_1 \right\rangle \otimes H \left| x_2 \right\rangle \otimes \dots \otimes H \left| x_n \right\rangle =

(применим к каждому отдельному кубиту)

=\frac{1}{\sqrt{2^n}} \Big(\sum \limits_{y_1 \in \lbrace 0,1 \rbrace} (-1)^{x_1 \cdot y_1} \left| y_1 \right\rangle \otimes \sum \limits_{y_2 \in \lbrace 0,1 \rbrace} (-1)^{x_2 \cdot y_2} \left| y_2 \right\rangle \otimes \ldots \otimes \sum \limits_{y_n \in \lbrace 0,1 \rbrace} (-1)^{x_n \cdot y_n} \left| y_n \right\rangle \Big) =

(вынесем -1 за знак тензорного произведения)

=\frac{1}{\sqrt{2^n}} \sum \limits_{(y_1, y_2, \ldots, y_n) \in \lbrace 0,1 \rbrace^n} (-1)^{x_1 \cdot y_1 \oplus x_2 \cdot y_2 \oplus \dots \oplus x_n \cdot y_n} \left| {y_1 y_2 \ldots y_n} \right\rangle =

(переобозначим для компактности)

=\frac{1}{\sqrt{2^n}} \sum \limits_{y \in \lbrace 0,1 \rbrace^n} (-1)^{(x, y)} \left| y \right\rangle,

где \oplus — сложение по модулю 2, а (x, y) — скалярное произведение по модулю 2.

В начальный момент времени система находится в состоянии \left| {q_0} \right\rangle \left| {q_1} \right\rangle = \left| 0 \right\rangle^{\otimes n} \left| 1 \right\rangle и под действием преобразования Адамара перейдёт в состояние

H^{\otimes n} \left| 0 \right\rangle^{\otimes n} (H \left| 1 \right\rangle) = \frac{1}{\sqrt{2^{n+1}}}\displaystyle\sum_{x=0}^{2^n} \left| x \right\rangle \big(\left| 0 \right\rangle - \left|1 \right\rangle \big).

Оракул U_f переведёт систему в состояние

\frac{1}{\sqrt{2^{n+1}}}\displaystyle\sum_{x=0}^{2^n} \left| x \right\rangle \Big(\left| {f(x)} \right\rangle \oplus \big(\left| 0 \right\rangle - \ \left| 1 \right\rangle \big) \Big) = \frac{1}{\sqrt{2^{n+1}}}\displaystyle\sum_{x=0}^{2^n} \left| x \right\rangle \Big( \left| {f(x)} \right\rangle - \left| {f(x) \oplus 1} \right\rangle \Big).

Заметьте, что если на некотором фиксированном аргументе x функция f(x) примет значение 0, то коэффициент перед \left| x \right\rangle не изменится, иначе — поменяет знак. Исходя из этого наблюдения можно переписать текущее состояние как

\frac{1}{\sqrt{2^{n+1}}}\displaystyle\sum_{x=0}^{2^n} (-1)^{f(x)} \left| x \right\rangle \big( \left| 0 \right\rangle - \left| 1 \right\rangle \big).

На данном этапе отбросим неиспользуемый далее последний кубит и применим повторное преобразование Адамара к первым n кубитам, чтобы получить состояние

H^{\otimes n} \Big( \frac{1}{\sqrt{2^n}}\displaystyle\sum_{x=0}^{2^n} (-1)^{f(x)} \left| x \right\rangle \Big) = \frac{1}{2^n}\displaystyle\sum_{y=0}^{2^n} \Big( \displaystyle\sum_{x=0}^{2^n} (-1)^{f(x) + x \cdot y} \Big) \left| y \right\rangle.

А так как x \cdot 0^n = 0, вероятность наблюдения исхода \left| {0^n} \right\rangle равна

\frac{1}{2^n}\displaystyle\sum_{x=0}^{2^n} (-1)^{f(x)} = 
\begin{cases}
   1 &\text{, если $f(x)$ константна}\\
   0 &\text{, если $f(x)$ сбалансирована}
\end{cases}.

Вот и всё. Мы рассмотрели достаточно упрощённую модель квантовых вычислений, написали мини-фреймворк на питоне для их эмуляции на классическом компьютере и разобрали один из самых простых «настоящих» квантовых алгоритмов — алгоритм Дойча-Йожи. Вот котик для тех, кто не поленился, прочитал статью и разобрался в теме:

.

→ С кодом можно ознакомиться на GitHub

Спасибо за прочтение! Ставьте лайки, подписывайтесь на профиль, оставляйте комментарии, занимайтесь полезными делами. Дополнения приветствуются.

Литература


  • [1] М. Вялый, «Квантовые алгоритмы: возможности и
    ограничения»
  • [2] А. С. Холево, «Введение в квантовую теорию информации»
  • [3] M. A. Nielsen, «Quantum Information Theory»

P.S: Формулы на хабре — отстой. Спасибо Roman Parpalak за сервис upmath.me, без него пост с таким количеством формул был бы невозможен.
Tags:квантовый компьютерквантовые алгоритмыквантовые вычисленияpythonnumpy
Hubs: Python Algorithms Mathematics
+109
44.2k 381
Comments 61
Popular right now
Python для работы с данными
December 7, 202031,500 ₽Нетология
Python QA Engineer
December 21, 202060,000 ₽OTUS
Python для анализа данных
December 31, 202018,990 ₽Level UP