Qiskit 예제 - 양자 위상 추정 (Quantum Phase Estimation, QPE)

 
이번 포스팅에서는 양자 위상 추정 방법에 대해 알아보도록 하겠습니다. 양자 위상 추정(Quantum Phase Estimation, QPE)은 고윳값 문제를 해결하는 알고리즘으로 양자 화학, 양자 기계 학습에서 활용되는 중요한 알고리즘입니다.
 
양자 위상 추정을 이용해 풀고 싶은 문제는 아래와 같은 고윳값 문제입니다
\begin{eqnarray} A {\pmb x} = \lambda {\pmb x} \end{eqnarray}
 
양자 위상 추정은 이 문제에 대해서, 주어진 연산자가 유니타리 연산자일 때 고윳값을 찾는 문제입니다.
 

유니타리 연산자란?

유니타리 연산자의 정의는 아래와 같습니다:
\begin{eqnarray} U^\dagger U =  I\end{eqnarray}
여기서 $U^\dagger$는 주어진 연산자 U의 Complex conjugate를 취하고, 그 결과를 Transpose 한 것을 말합니다. 
 
이러한 성질로 인해 유니타리 연산자의 경우, 고윳값의 절댓값이 1이 됩니다. 이는 아래와 같이 쉽게 증명할 수 있습니다:
\begin{eqnarray} U | u \rangle = \lambda | u \rangle  \Longrightarrow \langle u | U^\dagger U | u \rangle = \langle u | \lambda^* \lambda | u \rangle  \Longrightarrow |\lambda| = 1 \end{eqnarray} 
 
그럼, 이러한 결과에 따라 어떠한 유니타리 연산자에 대한 고윳값을 아래와 같이 위상에 대한 식으로 표현할 수 있습니다:
\begin{eqnarray} \lambda = e^{2\pi i \phi} \end{eqnarray}
여기서 위상 $\phi$는 $0 \le \phi \le 1$ 사이에 존재합니다. 때문에, 아래와 같이 이진 소숫법으로 아래와 같이 적을 수 있습니다:
\begin{eqnarray} \phi = 0 . \phi_1 \phi_2 \cdots \phi_n  \end{eqnarray} 
여기서 각각의 $\phi_i$는 0 또는 1로 표현이 됩니다.  (이진 소숫법에 대해서는 양자 푸리에 변환 포스팅을 참고해 주세요)
2024.12.30 - [Quantum Computing] - Qiskit을 이용한 양자 컴퓨팅 예제 (7) - 양자 푸리에 변환 (Quantum Fourier Transform)
 
결국, 위상 $\phi$를 알아낸다는 것은 고윳값 $\lambda$를 알아낸다는 것과 동일한 얘기가 되므로, QPE의 목적은 위상 $\phi$를 찾아내는 것이 될 것입니다.
 

양자 위상 추정 알고리즘

1) 초기 상태의 준비

양자 위상 추정 알고리즘을 이용하기 위해 n개의 측정 큐비트를 준비하고, $|0\rangle$ 상태로 초기화 합니다. 그럼 $| 0 \rangle^{\otimes n}$ 상태가 주어질 것입니다. 그리고 $U$에 대한 고유 상태 $| \psi \rangle$를 타겟 큐비트로 준비합니다. 따라서, 주어진 초기 상태는 아래와 같이 적을 수 있습니다:
\begin{eqnarray} | \psi_I \rangle = | 0 \rangle^{\otimes n} \otimes | \psi \rangle. \end{eqnarray} 
여기서 유니타리 연산자에 대한 고유 상태를 미리 알고 있어야 한다는 점에 유의하시기 바랍니다. 즉, 양자 위상 추정 알고리즘은 고윳값 문제에서 고유 벡터를 구하는 것이 아니라, 주어진 고유 벡터에 대해 고윳값을 계산하는 알고리즘 입니다. 
 
물리적으로는, 특정 고유 상태가 주어진 상황에서 이를 이용해 에너지나 다른 물리적 양을 계산하는 데 활용할 수 있습니다. 또한, 정확한 고유 상태를 준비할 수 없는 경우에도, 근사적으로 준비된 상태를 사용하여 해당 고윳값을 추정할 수 있습니다. 심지어 고유 상태를 알기 어려운 경우에는, 고유 상태를 구하는 다른 알고리즘과 결합하여 이를 구현할 수도 있습니다.
 

2) 측정 큐비트의 H 게이트와 제어 $U$ 연산자 적용하기

다음 단계로, 주어진 초기 상태에 대해서 $n$번째 큐비트에 H 게이트를 적용하고 제어 $U^{2^{n-1}}$ 연산자를 적용합니다. 
 
첫번째 측정 큐비트에 이를 적용하면 아래와 같이 쓰여집니다:
\begin{eqnarray} \frac{1}{\sqrt{2}}( | 0 \rangle + | 1\rangle ) \otimes | \psi \rangle = \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle + |1 \rangle | \psi \rangle \right) \end{eqnarray}
그 다음 여기에 컨트롤 큐비트가 $|0\rangle$일 때는 작동하지 않고, $|1\rangle$인 경우에만 작동하는 제어 U 연산자 타겟 큐비트에 적용합니다:
\begin{eqnarray} \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle +  | 1\rangle U | \psi \rangle \right) = \frac{1}{\sqrt{2}} \left( | 0\rangle | \psi \rangle + |1 \rangle e^{0.\phi_1 \phi_2 \phi_3 \cdots \phi_n} | \psi \rangle \right) = \frac{1}{\sqrt{2}} ( | 0 \rangle + e^{2\pi i 0. \phi_1 \phi_2 \cdots \phi_n} | 1 \rangle ) \otimes |\psi \rangle. \end{eqnarray}
여기서 $U=U^{2^0}$로 쓸 수 있습니다. 이를 통해 고유 상태의 위상 정보를 첫번째 레지스터에 인코딩 할 수 있습니다. 
 
두번째 측정 큐비트에는 H게이트와 제어 $U^{2^1}$ 연산자를 적용합니다:
\begin{eqnarray}  \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle +  | 1 \rangle U^2 | \psi \rangle \right) = \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle +  |1 \rangle e^{2 \pi i (2 \times 0.\phi_1 \phi_2 \cdots \phi_n )} | \psi \rangle \right). \end{eqnarray}
위 식에서 $\phi_1$은 이진 표기법에 의해 0 또는 1만을 가지므로, $e^{2 \pi i \phi_1} = 1$임을 알 수 있습니다. (이는 $\phi_1$뿐만 아니라 모든 $\phi_i$에 대해서 성립합니다.) 따라서, 주어진 식은
\begin{eqnarray} \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle + e^{2 \pi i 0.\phi_2 \phi_3 \cdots \phi_n} | 1 \rangle | \psi \rangle \right). \end{eqnarray}
이 될 것입니다. 
 
세번째 측정 큐비트에는 H 게이트와 제어 $U^{2^2}$ 연산자를 적용합니다:
\begin{eqnarray} \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle + | 1 \rangle U^{2^2} | \psi \rangle \right) = \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle + |1 \rangle e^{2\pi i (2^2 \times 0. \phi_1 \phi_2 \cdots \phi_n) } | \psi \rangle \right) = \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle + e^{2\pi i 0.\phi_3 \phi_4 \cdots \phi_n} | 1 \rangle |\psi \rangle \right). \end{eqnarray}
 
이러한 규칙을 일반화하면, k 번째 큐비트에 대해 아래와 같은 결과를 얻을 수 있을 것입니다:
\begin{eqnarray} \frac{1}{\sqrt{2}} \left( |0 \rangle | \psi \rangle + | 1\rangle U^{2^{k-1}} | \psi \rangle \right) = \frac{1}{\sqrt{2}} \left( | 0 \rangle | \psi \rangle + | 1\rangle e^{2 \pi i (2^{k-1}) 0.\phi_1 \phi_2 \cdots \phi_n } | \psi \rangle \right) = \frac{1}{\sqrt{2}} \left( |0 \rangle + e^{2\pi i 0. \phi_k \phi_{k+1} \cdots \phi_n} | 1\rangle \right) | \psi \rangle. \end{eqnarray}
 
위에서의 규칙을 이용해 주어진 연산자에 대한 위상을 찾기 위해 아래와 같은 최종 상태를 구성할 수 있습니다:
\begin{eqnarray} |\psi_F \rangle &=& \frac{1}{\sqrt{2^n}} \left[ \left( |0\rangle + e^{2 \pi i 0.\phi_1 \phi_2 \cdots \phi_n}|1\rangle \right) \otimes \left( | 0 \rangle + e^{2 \pi i 0.\phi_2 \phi_3 \cdots \phi_n} | 1 \rangle \right) \otimes \cdots \otimes \left( |0 \rangle + e^{2 \pi i 0. \phi_n} |1 \rangle | \psi \rangle \right) \right] | \psi \rangle \end{eqnarray}
 
위 식은 모든 가능한 상태와 각각의 상태에 해당하는 위상을 갖는 식이 됩니다. 위에 식을 좀 더 간결하게 쓰면 아래와 같이 쓸 수 있습니다:
\begin{eqnarray} |\psi_F \rangle = \frac{1}{\sqrt{2^n}} \sum_{k=1}^{2^n-1} e^{2 \pi i k } | k\rangle | \psi \rangle,\end{eqnarray}
여기서 $|k\rangle$는 주어진 $k$에 대해 이진수로 표현된 상태를 의미하며, $e^{2\pi i k}$는 그 상태에 해당하는 위상을 의미합니다.
 

3) 푸리에 역변환 적용하기

마지막으로, 주어진 최종 상태에 대해 역 푸리에 변환을 적용하면 위상을 얻어낼 수 있습니다:
\begin{eqnarray} {\rm QFT}^{-1} | \psi_F \rangle = \frac{1}{2^n} \sum_{m=0}^{2^n-1} \left( \sum_{k=0}^{2^n-1}e^{2\pi i k (\phi - m/2^n)} \right) |m\rangle |\psi \rangle. \end{eqnarray} 
위에서 만약 $\phi$가 $m/2^n$과 비슷한 값을 가지지 않는다면, 그러한 상태들의 합은 위상차에 의해 상쇄될 것입니다. 따라서, $\phi \approx m/2^n$을 만족시키는 상태가 가장 많이 측정될 것입니다. 
 
다시 말해, 위에서 주어진 상태에 역 푸리에 변환을 적용하여 측정할 경우, 아래와 같이 상태가 가장 많이 측정될 것입니다:
\begin{eqnarray} {\rm QFT}^{-1} | \psi_F \rangle \approx |m \rangle |\psi \rangle, \end{eqnarray}
그리고 이 상태에서 $m$은 $\phi \approx m/2^n$을 만족하도록 주어질 것입니다. 측정된 상태 $| m \rangle$와 측정 큐비트의 수 $n$도 알기 때문에, 이 조건으로부터 $\phi$ 값 역시 찾아낼 수 있습니다.
 

Qiskit으로 구현한 QPE 알고리즘

이러한 QPE 알고리즘을 Qiskit으로 구현해 보겠습니다. 우리가 풀 문제는 아래와 같습니다:
\begin{eqnarray} U | \psi \rangle = \lambda | \psi \rangle \end{eqnarray}
여기서 $U$를 아래와 같이 주겠습니다:
\begin{eqnarray} U = \begin{pmatrix} 1 & 0 \\ 0 & e^{2\pi i \phi} \end{pmatrix}\end{eqnarray}
유니타리 연산자에서 위상 $\phi=1/5$로 설정하겠습니다.
 
주어진 유니타리 연산자에 대한 고유 벡터는 $| \psi \rangle = | 1 \rangle$이 됩니다. 따라서 위에 고윳값 문제는 쉽게 풀립니다: 
\begin{eqnarray} U | 1 \rangle = \begin{pmatrix} 1 & 0 \\ 0 & e^{2\pi i \phi} \end{pmatrix} \begin{pmatrix} 0 \\ 1 \end{pmatrix} = e^{2\pi i \phi} \begin{pmatrix} 0 \\ 1 \end{pmatrix} \end{eqnarray}
여기서 위상 $\phi = 1/5$이므로, 이 답을 Qiskit으로 풀어서 비교해 보겠습니다. 
 
우선 필요한 라이브러리들을 불러오고, 4개의 측정 큐비트와 1개의 타겟 큐비트를 준비하겠습니다: 

from qiskit import QuantumCircuit, transpile
from qiskit.circuit.library import QFT
from qiskit.visualization import plot_bloch_multivector, plot_histogram
from qiskit_aer import AerSimulator
import matplotlib.pyplot as plt
import numpy as np

num_qubits = 4  # 계수 비트 개수 (precision)
ancilla_qubit = 1  # 타겟 큐빗

# 전체 큐비트: 계수 비트 + 타겟 큐빗
total_qubits = num_qubits + ancilla_qubit

 
그 다음 회로를 구성하겠습니다. 4개의 큐비트는 H 게이트를 적용할 것이며, 타겟 큐비트는 우리가 알고 있는 고유 벡터 $|1\rangle$ 상태로 만들어줍니다:

# 양자 회로 초기화
qc = QuantumCircuit(total_qubits)

# 계수 비트 초기화 (|0⟩ -> |+⟩ 상태로 변환)
qc.h(range(num_qubits))

# 타겟 큐빗 초기화: |1⟩
qc.x(num_qubits)

 
그 다음 Controlled-U 연산자를 구성하여, 위에서 설명한 알고리즘에 맞게 $U^{2^k}$를 k번째 큐비트에 적용해 주겠습니다 (지금 풀고 있는 예제에서 주어진 제어 U 연산자는 Controlled Phase gate라고 불리며, cp 명령어를 통해 불러올 수 있습니다):

# Controlled-U 연산 정의
def controlled_unitary_with_phase(qc, control, target, phase, power):
    """Controlled U gate with phase adjusted by 2^j."""
    qc.cp(phase * 2 * np.pi * (2**power), control, target)

# Controlled-U 연산 추가
for qubit in range(num_qubits):
    controlled_unitary_with_phase(qc, qubit, num_qubits, 1/5, qubit)

 
마지막으로 역 푸리에 변환을 적용하고 측정하겠습니다:

# 역 QFT 적용
qc.append(QFT(num_qubits, do_swaps=True).inverse(), range(num_qubits))

# 계수 비트 측정
qc.measure_all()

 
이렇게 주어진 회로를 1024번 시뮬레이션하고, 이에 대한 결과를 히스토그램으로 출력하겠습니다:

# 시뮬레이션 실행
simulator = AerSimulator()
compiled_circuit = transpile(qc, simulator)
result = simulator.run(compiled_circuit, shots=1024).result()
counts = result.get_counts()

# 결과 출력
print("측정 결과:", counts)

fig1 = plot_histogram(counts)
fig1.savefig("histogram.png", dpi=300)  #png로 저장

 
양자회로도 역시 함께 출력하겠습니다:

# 회로 시각화
fig2 = qc.draw(output="mpl")  # Matplotlib 객체 반환
fig2.savefig("quantum_circuit_qpe.png", dpi=300)  # PNG로 저장

 
먼저 히스토그램 실행 결과를 보도록 하겠습니다:

Qiskit 코드 실행으로 얻은 히스토그램.

결과 해석:

보시는 것처럼 $|10011\rangle$ 상태가 가장 많이 나오는 것을 알 수 있습니다. 이 순서는 $|q_4 q_3 q_2 q_1 q_0 \rangle $ 순서이므로, 가장 앞에 놓인 상태는 타겟 큐비트(고유 벡터)인 $|1 \rangle$ 상태를 의미합니다. 따라서 측정 큐비트에서 가장 많이 측정된 상태는 $|0011\rangle$을 의미하며, 이는 십진법으로 표현했을 때 $|3 \rangle$이 됩니다. 
 
앞서 설명한 조건에서 위상 $\phi$는
\begin{eqnarray} \phi \approx \frac{m}{2^n} \end{eqnarray}
의 조건을 가지므로, 우리는 $\phi$를 아래와 같이 계산할 수 있습니다:
\begin{eqnarray} \phi = \frac{3}{2^4} = \frac{3}{16} = 0.1875 \end{eqnarray}
이는 실제 직접 풀어서 얻은 0.2와 유사합니다. 큐비트 수를 늘려가면서 풀다 보면, 보다 0.2에 가까운 값을 얻을 수 있을 것입니다! 
(주어진 코드에서 측정 큐비트의 수(num_qubits)를 늘려보면 쉽게 확인할 수 있습니다. 측정 큐비트가 6개만 되어도, $\phi$ 값은 약 0.203 정도로 실제값과 매우 비슷한 결과가 나옵니다!)
 

양자 회로도

위에 예제에 대한 양자 회로도 출력 결과는 아래와 같습니다: 

Quantum Phase Estimation에 대한 양자 회로도