Алгоритм оценки фазы

Этот алгоритм лежит в основе многих других, включая алгоритм Шора. Задача, которую он решает, заключается в следующем.

Пусть дан унитарный оператор \(U\) в виде черного ящика-оракула, и его собственный вектор \(|\psi\rangle)\) в виде состояния, которое невозможно скопировать из-за законов квантовой механики. Поскольку \(U\) является унитарным, все его собственные значения и собственные вектора имеют вид \[ U|\psi\rangle = e^{2\pi i \theta}|\psi\rangle. \] Необходимо вычислить \(n\)-битовую оценку значения \(\theta\).

Алгоритм оценки фазы

Общая идея

Если бы мы могли записать фазу \(\theta\) в \(n\) разных кубитов в базисе Фурье, то для решения было бы достаточно применить обратное преобразование Фурье и измерить регистр. Даже если фаза не целая, как в нашем случае, мы получим какую-то её оценку (может быть показано, что вероятность измерения выше \(4/\pi^2 \approx 40%\)).

Фазовая отдача (phase kickback)

Что получится, если мы выполним операцию \(CU\), поместив в управляющий кубит состояние \(\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)\), а в управляемый — \(|\psi\rangle\)?

Получим

\begin{align*} CU\Big(\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle) \otimes |\psi\rangle\Big) &= \frac{1}{\sqrt{2}}\Big(|0\rangle \otimes |\psi\rangle + |1\rangle \otimes e^{2\pi i \theta}|\psi\rangle\Big) = \\ &= \frac{1}{\sqrt{2}}(|0\rangle + e^{2\pi i \theta} |1\rangle) \otimes |\psi\rangle. \end{align*}

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

Повышение точности

Повторяя этот процесс с разными кубитами \(j\) и разными степенями \(CU^{2^j}\), получим следующее состояние:

\begin{align*} \vert \psi' \rangle &= \frac{1}{2^{n/2}} \big(|0\rangle + e^{2\pi i\theta 2^{n-1}}|1\rangle\big) \otimes \cdots \otimes \big(|0\rangle + e^{2\pi i\theta 2^1}|1\rangle\big) \otimes \big(|0\rangle + e^{2\pi i\theta 2^0}|1\rangle\big) \otimes |\psi\rangle \\ &= \frac{1}{2^{n/2}} \sum_{k=0}^{2^n-1} e^{2\pi i \theta k}|k\rangle \otimes |\psi\rangle. \end{align*}

После обратного квантового преобразования Фурье получаем \[ \vert \psi'' \rangle = \frac{1}{2^{n/2}} \sum_{x=0}^{2^n-1} \sum_{k=0}^{2^n-1} e^{-\frac{2\pi ik}{2^n} (x - 2^n\theta)}|x\rangle \otimes |\psi\rangle. \]

Измерив первый регистр, с большой вероятностью получим \(x = 2^n \theta\). Для большей точности потребуется больше кубитов в этом регистре.

Код

Реализуем описанное выше в виде кода.

import cirq
import numpy as np

def make_phase_estimation(n, U, u_bit):
    qubits = cirq.LineQubit.range(n)
    phase_estimator = cirq.Circuit(cirq.H.on_each(*qubits))
    for i, bit in enumerate(qubits):
        phase_estimator.append(cirq.ControlledGate(U).on(bit, u_bit) ** (2**(n - i - 1)))
    phase_estimator.append(cirq.qft(*qubits[::-1], without_reverse=True, inverse=True))
    phase_estimator.append(cirq.measure(*qubits, key='theta'))
    return phase_estimator

Если операция \(U\) такова, что \(2^n \theta\) — целое, то измерения дадут точное значение:

n1 = 3
u1 = cirq.NamedQubit('u')
pe1 = make_phase_estimation(n1, cirq.T, u1)
pe1.insert(0, cirq.X(u1))
print(pe1)

sim = cirq.Simulator()
res1 = sim.run(pe1, repetitions = 10)
print(res1.measurements['theta'])
theta1 = np.sum(2 ** np.arange(n1) * res1.measurements['theta'], axis=1) / 2**n1
print(theta1)

Здесь мы всегда получаем \(\theta = 0.125\) потому что мы находим фазу для оператора \(T = Z^{1/4}\) и собственного значения \(|1\rangle\) (для того, чтобы начать с этого состояния мы и вставили \(X\) в схему): \[ T|1\rangle = e^{2\pi i/8}|1\rangle. \]

Если же \(2^n \theta\) нельзя представить точно в \(n\) битах, то мы получим какое-то приближение:

n2 = 3
u2 = cirq.NamedQubit('u')
pe2 = make_phase_estimation(n2, cirq.Z**(1/3), u2)
pe2.insert(0, cirq.X(u2))
print(pe2)

res2 = sim.run(pe2, repetitions = 10)
print(res2.measurements['theta'])
theta2 = np.sum(2 ** np.arange(n2) * res2.measurements['theta'], axis=1) / 2**n2
print(theta2)
print(res2.histogram(key = 'theta'))

И чем больше мы используем кубитов, тем точнее результат:

n3 = 4
u3 = cirq.NamedQubit('u')
pe3 = make_phase_estimation(n3, cirq.Z**(1/3), u3)
pe3.insert(0, cirq.X(u3))
print(pe3)

res3 = sim.run(pe3, repetitions = 1000)
print(res3.histogram(key = 'theta'))

n4 = 5
u4 = cirq.NamedQubit('u')
pe4 = make_phase_estimation(n4, cirq.Z**(1/3), u4)
pe4.insert(0, cirq.X(u4))
print(pe4)

res4 = sim.run(pe4, repetitions = 1000)
print(res4.histogram(key = 'theta'))

Стандартная функция histogram собирает биты в обратном порядке к тому, что мы ожидаем: \(4 = 100_2\), \(12 = 1100_2\), \(20 = 10100_2\). Получаем такие приближения: \(001, 0011, 00101, \ldots\), или \(\theta \approx 0.125, 0.1875, 0.15625, \ldots \) при искомом значении \(\theta = 1/6 = 0.1(6) = 0.1666\ldots = 0.0(01)_2 = 0.00101010101\ldots_2\).