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

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

Пусть дан унитарный оператор \(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)

0: ───H───@───────────#3──────────────M('theta')───
          │           │               │
1: ───H───┼───@───────#2──────────────M────────────
          │   │       │               │
2: ───H───┼───┼───@───qft[norev]^-1───M────────────
          │   │   │
u: ───X───Z───S───T────────────────────────────────
[[1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]]
[0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125]

Здесь мы всегда получаем \(\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'))

0: ───H───@──────────────────────────────#3──────────────M('theta')───
          │                              │               │
1: ───H───┼──────────@───────────────────#2──────────────M────────────
          │          │                   │               │
2: ───H───┼──────────┼─────────@─────────qft[norev]^-1───M────────────
          │          │         │
u: ───X───Z^(-2/3)───Z^(2/3)───Z^(1/3)────────────────────────────────
[[0 0 1]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [1 0 0]
 [0 1 0]
 [1 0 0]
 [0 0 1]
 [1 0 0]]
[0.5   0.125 0.125 0.125 0.125 0.125 0.25  0.125 0.5   0.125]
Counter({4: 7, 1: 2, 2: 1})

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

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'))

0: ───H───@────────────────────────────────────────#4──────────────M('theta')───
          │                                        │               │
1: ───H───┼─────────@──────────────────────────────#3──────────────M────────────
          │         │                              │               │
2: ───H───┼─────────┼──────────@───────────────────#2──────────────M────────────
          │         │          │                   │               │
3: ───H───┼─────────┼──────────┼─────────@─────────qft[norev]^-1───M────────────
          │         │          │         │
u: ───X───Z^(2/3)───Z^(-2/3)───Z^(2/3)───Z^(1/3)────────────────────────────────
Counter({12: 679, 4: 172, 2: 54, 8: 31, 0: 12, 15: 10, 10: 10, 6: 8, 14: 4, 1: 4, 11: 4, 7: 3, 13: 3, 3: 2, 5: 2, 9: 2})

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'))

0: ───H───@───────────────────────────────────────────────────#5──────────────M('theta')───
          │                                                   │               │
1: ───H───┼──────────@────────────────────────────────────────#4──────────────M────────────
          │          │                                        │               │
2: ───H───┼──────────┼─────────@──────────────────────────────#3──────────────M────────────
          │          │         │                              │               │
3: ───H───┼──────────┼─────────┼──────────@───────────────────#2──────────────M────────────
          │          │         │          │                   │               │
4: ───H───┼──────────┼─────────┼──────────┼─────────@─────────qft[norev]^-1───M────────────
          │          │         │          │         │
u: ───X───Z^(-2/3)───Z^(2/3)───Z^(-2/3)───Z^(2/3)───Z^(1/3)────────────────────────────────
Counter({20: 685, 12: 169, 4: 45, 28: 22, 24: 18, 2: 13, 8: 7, 18: 7, 16: 5, 10: 4, 6: 3, 31: 3, 0: 3, 13: 2, 1: 2, 26: 2, 15: 1, 9: 1, 29: 1, 7: 1, 19: 1, 23: 1, 22: 1, 11: 1, 17: 1, 30: 1})

Стандартная функция 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\).