Python으로 큐비트 시뮬레이션하기 (고전적으로)
양자 컴퓨터를 직접 다뤄 보면 어떨지 궁금해 본 적이 있습니까? 누구나 지하실에 초저온 양자 프로세서를 둘 수는 없지만(아직은), 그다음으로 좋은 방법이 있습니다. 바로 평범한 Python 코드로 양자 컴퓨터를 시뮬레이션하는 것입니다. 물론 우스울 정도로 비효율적이지만, 놀라울 만큼 많은 것을 배울 수 있습니다! 이제 간단한 양자 시뮬레이터를 처음부터 만들어 보겠습니다. 확장성은 (전혀) 좋지 않겠지만, 큐비트와 중첩을 직접 다뤄 보고 이것들을 소프트웨어로 시뮬레이션하면 어떤 모습인지 체감할 수 있을 것입니다.
양자 정보의 기본 단위는 큐비트입니다. 0 또는 1 중 하나로 확실히 고정되는 고전 비트와 달리, 큐비트는 두 상태가 동시에 결합된 중첩 상태로 존재할 수 있습니다. 흔히 |ψ>로 표기하는 단일 큐비트의 상태는 수학적으로 다음과 같이 나타냅니다: α|0> + β|1>. 여기서 |0>과 |1>은 두 가지 기본 기저 상태(고전적인 0과 1에 해당)를 나타내며, α와 β는 진폭이라고 부르는 복소수입니다. 이 진폭은 임의의 값이 아닙니다. 진폭의 절댓값을 제곱한 |α|²와 |β|²는 각각 큐비트를 측정했을 때 |0> 상태 또는 |1> 상태로 관측될 확률을 나타냅니다. 핵심(Core) 규칙은 이 확률들의 합이 반드시 1이 되어야 한다는 것입니다: |α|² + |β|² = 1.
큐비트가 여러 개가 되면 상황은 기하급수적으로 복잡해집니다. n개의 큐비트로 이루어진 시스템에는 2^n가지의 가능한 고전적 구성이 있습니다(예를 들어 큐비트가 두 개라면 00, 01, 10, 11). 이 n-큐비트 시스템의 양자 상태를 나타내려면 2^n개의 기저 상태 각각에 대한 복소 진폭이 필요합니다. 예를 들어 3-큐비트 시스템(n=3)은 기저 상태 |000>, |001>, …, |111>을 갖는 8차원 복소 벡터 공간에 존재합니다. 시스템 상태에 대한 완전한 기술은 이 2^n개의 복소 진폭을 담은 리스트, 즉 상태 벡터입니다.
이러한 기하급수적 증가는 양자 컴퓨팅이 지닌 힘의 핵심이자, 고전적으로 시뮬레이션하기가 그토록 어려운 이유이기도 합니다. 큐비트 수를 두 배로 늘리면 고전 상태의 수는 제곱이 되므로, 시뮬레이션에 필요한 메모리와 처리 능력이 기하급수적으로 늘어납니다. 우리 시뮬레이터는 이 상태 벡터를 효과적으로 관리하여 실제 양자 컴퓨터의 하드웨어에서 수행될 작업을 시뮬레이션합니다.
QuantumSystem 클래스로 캡슐화된 Python 구현을 살펴봅시다.
import math
import random
import cmath
from typing import List, Tuple, Sequence
SQRT2 = math.sqrt(2.0)
class QuantumSystem:
def __init__(self, num_qubits: int):
if num_qubits < 1:
raise ValueError("Number of qubits must be at least 1.")
self.num_qubits: int = num_qubits
self.num_states: int = 1 << num_qubits
self.amplitudes: List[complex] = [complex(0.0)] * self.num_states
self.amplitudes[0] = complex(1.0)
def _validate_qubit_index(self, qubit_index: int):
if not (0 <= qubit_index < self.num_qubits):
raise ValueError(
f"Invalid qubit index: {qubit_index}. "
f"Must be between 0 and {self.num_qubits - 1}."
)
def _validate_qubit_indices(self, *indices: int):
seen = set()
for index in indices:
self._validate_qubit_index(index)
if index in seen:
raise ValueError(f"Duplicate qubit index specified: {index}")
seen.add(index)
def measure(self) -> Tuple[int, ...]:
probabilities: List[float] = [abs(amp) ** 2 for amp in self.amplitudes]
measured_state_index: int = random.choices(
population=range(self.num_states), weights=probabilities, k=1
)[0]
self.amplitudes = [complex(0.0)] * self.num_states
self.amplitudes[measured_state_index] = complex(1.0)
measured_bits = tuple(
(measured_state_index >> bit_pos) & 1
for bit_pos in range(self.num_qubits)
)
return measured_bits
def get_probabilities(self) -> List[float]:
return [abs(amp) ** 2 for amp in self.amplitudes]
def apply_h(self, target_qubit: int):
self._validate_qubit_index(target_qubit)
new_amplitudes = self.amplitudes[:]
target_mask = 1 << target_qubit
for i in range(self.num_states // 2):
lower_mask = target_mask - 1
upper_mask = ~lower_mask
i0 = (i & lower_mask) | ((i << 1) & upper_mask)
i1 = i0 | target_mask
amp0 = self.amplitudes[i0]
amp1 = self.amplitudes[i1]
new_amp0 = (amp0 + amp1) / SQRT2
new_amp1 = (amp0 - amp1) / SQRT2
new_amplitudes[i0] = new_amp0
new_amplitudes[i1] = new_amp1
self.amplitudes = new_amplitudes
def apply_cnot(self, control_qubit: int, target_qubit: int):
self._validate_qubit_indices(control_qubit, target_qubit)
control_mask = 1 << control_qubit
target_mask = 1 << target_qubit
new_amplitudes = self.amplitudes[:]
for i0 in range(self.num_states):
if (i0 >> control_qubit) & 1:
i1 = i0 ^ target_mask
if i0 < i1:
new_amplitudes[i0], new_amplitudes[i1] = \
self.amplitudes[i1], self.amplitudes[i0]
self.amplitudes = new_amplitudes
def apply_t(self, target_qubit: int):
self._validate_qubit_index(target_qubit)
phase_shift = cmath.exp(1j * cmath.pi / 4.0)
target_mask = 1 << target_qubit
for i in range(self.num_states):
if (i >> target_qubit) & 1:
self.amplitudes[i] *= phase_shift
def apply_original_pi_over_eight(self, target_qubit: int):
self._validate_qubit_index(target_qubit)
phase_zero = cmath.exp(-1j * cmath.pi / 8.0)
phase_one = cmath.exp(1j * cmath.pi / 8.0)
target_mask = 1 << target_qubit
for i in range(self.num_states):
if (i >> target_qubit) & 1:
self.amplitudes[i] *= phase_one
else:
self.amplitudes[i] *= phase_zero
def __repr__(self) -> str:
state_strs = []
for i, amp in enumerate(self.amplitudes):
if not cmath.isclose(amp, 0.0):
basis_state = format(i, f'0{self.num_qubits}b')
state_strs.append(f"{amp:.3f}|{basis_state}>")
if not state_strs:
return "QuantumSystem(num_qubits={}, state=Zero Vector)".format(
self.num_qubits
)
return " + ".join(state_strs)
def __str__(self) -> str:
return self.__repr__()
system = QuantumSystem(3)처럼 QuantumSystem 인스턴스를 생성하면 __init__ 메서드가 상태 벡터를 설정합니다. 이 메서드는 큐비트 수가 유효한지 확인하고, 이를 self.num_qubits에 저장한 뒤, 전체 상태 수(self.num_states = 2**num_qubits)를 계산합니다. 핵심(Core) 속성인 self.amplitudes는 self.num_states개의 복소수 0으로 이루어진 리스트로 초기화됩니다. 중요한 점은 맨 첫 번째 상태(인덱스 0, |00...0>에 해당)의 진폭을 complex(1.0)으로 설정한다는 것인데, 이는 양자 계산의 표준 시작 상태를 나타냅니다.
현재 상태를 확인할 수 있도록 __repr__ 메서드는 편리한 문자열 형식을 제공하여, 0이 아닌 진폭을 이진수로 쓴 해당 기저 상태와 함께 표시합니다(예: (0.707+0.000j)|00> + (0.707+0.000j)|11>).
양자 계산은 큐비트에 양자 게이트를 적용하는 방식으로 진행됩니다. 이 게이트들은 고전 논리 게이트(NOT, AND, OR 등)에 대응하는 양자적 개념입니다. 수학적으로 이 게이트들은 상태 벡터에 적용되는 유니터리 변환에 해당합니다. 게이트를 적용하면 진폭이 변경되며, 이를 통해 중첩이나 얽힘이 생성될 수 있습니다. 우리 클래스는 몇 가지 기본 게이트를 구현합니다.
apply_h에 구현된 아다마르 게이트는 중첩 상태를 만드는 데 필수적입니다. |0> 상태의 큐비트에 적용하면 (|0> + |1>) / sqrt(2)를 만들고, |1>에 적용하면 (|0> - |1>) / sqrt(2)를 만듭니다. 다중 큐비트 시스템에서 특정 target_qubit에 적용하면, 해당 큐비트의 위치에서만 차이가 나는 모든 기저 상태 쌍에 영향을 줍니다. 코드는 num_states / 2개의 쌍을 순회하면서 이를 효율적으로 처리합니다. 각 쌍 (i0, i1)(i0은 대상 큐비트가 0이고 i1은 1인 경우)에 대해, 아다마르 변환 규칙인 amp0' = (amp0 + amp1) / sqrt(2)와 amp1' = (amp0 - amp1) / sqrt(2)를 사용하여 새로운 진폭 amp0'와 amp1'를 계산합니다. 모든 연산이 잠재적으로 모든 진폭에 간접적으로 영향을 미친다는 점에 주목하십시오. 이는 양자 상태의 전일적인 특성을 반영합니다.
apply_cnot에 구현된 제어 NOT(CNOT) 게이트는 얽힘을 만드는 데 핵심적인 2큐비트 게이트입니다. 얽힘이란 큐비트들이 서로 떨어져 있어도 상관관계를 유지하는 기묘한 연결을 말합니다. 이 게이트는 control_qubit이 |1> 상태일 때, 그리고 그때에만 target_qubit의 상태를 뒤집습니다. 상태 벡터 표현에서 이는 제어 비트가 1이면서 대상 비트 위치에서만 차이가 나는 기저 상태 쌍의 진폭을 서로 맞바꾸는 것을 의미합니다. 구현은 상태 인덱스(i0)를 순회합니다. i0의 제어 비트가 1이면, 대상 비트가 뒤집힌 대응 인덱스 i1을 찾습니다(i1 = i0 ^ target_mask). 이중 교환을 피하기 위해 i0 < i1인 경우에만 self.amplitudes[i0]와 self.amplitudes[i1] 사이의 진폭 교환을 수행하므로, 관련된 각 쌍은 정확히 한 번만 처리됩니다.
위상 게이트는 측정 확률을 바꾸지 않고 진폭의 복소 위상을 변경합니다. apply_t 메서드는 표준 T 게이트를 구현하며, 이 게이트는 target_qubit이 1인 모든 기저 상태의 진폭에 exp(i*pi/4)의 위상 이동을 적용합니다. apply_original_pi_over_eight 메서드는 초기 설계에서 비롯된 약간 다른 게이트를 구현하며, 큐비트가 0이면 exp(-i*pi/8)를, 1이면 exp(i*pi/8)를 적용합니다. 두 메서드 모두 진폭을 순회하면서, 상태 인덱스 i가 target_qubit과 관련된 조건을 만족하면 적절한 복소 위상 인자를 곱하는 방식으로 동작합니다.
마지막으로, 고전적인 결과는 어떻게 추출할까요? 이는 measure 메서드에 구현된 측정을 통해 이루어집니다. 양자 측정은 확률적입니다. measure 메서드는 먼저 각 기저 상태 i의 진폭의 절댓값을 제곱하여(abs(amp)**2) 해당 상태가 측정될 확률을 계산합니다. 그런 다음 Python의 random.choices를 사용하여, 계산된 확률을 가중치로 삼아 가능한 상태 인덱스(range(self.num_states)) 중에서 가중 무작위 선택을 수행합니다. 이는 시스템이 붕괴하여 도달하는 단일 고전 상태의 인덱스 measured_state_index를 반환합니다. 그런 다음 시뮬레이션은 이 붕괴를 강제합니다. 상태 벡터 self.amplitudes는 측정된 상태의 진폭만 1.0으로 설정되고 나머지는 모두 0으로 재설정됩니다. 이 메서드는 measured_state_index에서 도출한 고전 비트(0과 1)의 튜플로 결과를 반환합니다.
예제를 통해 이러한 구성 요소가 실제로 동작하는 모습을 볼 수 있습니다. 먼저 3큐비트 시스템을 만들고 두 번째 큐비트에 아다마르 게이트를 적용해 봅시다.
print("--- Example 1: Original Sequence ---")
system1 = QuantumSystem(3)
print(f"Initial state: {system1}")
system1.apply_h(1)
print(f"After H(1): {system1}")
system1.apply_h(2)
print(f"After H(2): {system1}")
system1.apply_original_pi_over_eight(1)
print(f"After Pi/8(1): {system1}")
system1.apply_cnot(control_qubit=1, target_qubit=0)
print(f"After CNOT(1,0): {system1}")
print(f"Probabilities: {[f'{p:.3f}' for p in system1.get_probabilities()]}")
measurement_result = system1.measure()
print(f"Measured state: {measurement_result}")
print(f"State after measurement: {system1}")
다음으로 벨 상태를 준비해 봅시다.
print("--- Example 2: Bell State Preparation ---")
bell_system = QuantumSystem(2)
print(f"Initial state: {bell_system}")
bell_system.apply_h(0)
print(f"After H(0): {bell_system}")
bell_system.apply_cnot(control_qubit=0, target_qubit=1)
print(f"After CNOT(0,1): {bell_system}")
print("Measuring Bell state 10 times:")
counts = {}
for _ in range(10):
bell_system = QuantumSystem(2)
bell_system.apply_h(0)
bell_system.apply_cnot(control_qubit=0, target_qubit=1)
result = bell_system.measure()
counts[result] = counts.get(result, 0) + 1
print(f"Measurement counts: {counts}")
이제 그린버거-혼-차일링거 상태를 준비해 봅시다.
print("--- Example 3: GHZ State Preparation ---")
ghz_system = QuantumSystem(3)
print(f"Initial state: {ghz_system}")
ghz_system.apply_h(0)
print(f"After H(0): {ghz_system}")
ghz_system.apply_cnot(control_qubit=0, target_qubit=1)
print(f"After CNOT(0,1): {ghz_system}")
ghz_system.apply_cnot(control_qubit=0, target_qubit=2)
print(f"After CNOT(0,2): {ghz_system}")
print("Measuring GHZ state 10 times:")
counts_ghz = {}
for _ in range(10):
ghz_system = QuantumSystem(3)
ghz_system.apply_h(0)
ghz_system.apply_cnot(control_qubit=0, target_qubit=1)
ghz_system.apply_cnot(control_qubit=0, target_qubit=2)
result = ghz_system.measure()
counts_ghz[result] = counts_ghz.get(result, 0) + 1
print(f"Measurement counts: {counts_ghz}")
이 예제들은 벨 상태((|00> + |11>)/sqrt(2))와 그린버거–혼–차일링거(GHZ) 상태((|000> + |111>)/sqrt(2)) 같은 잘 알려진 얽힘 상태를 준비하는 과정을 보여 줍니다. (시뮬레이션을 위해 상태를 재설정한 뒤) 측정을 여러 번 실행하면 얽힘에 내재된 상관관계가 분명히 드러납니다. 벨 상태에서는 결과가 항상 (0, 0) 또는 (1, 1)이고, GHZ 상태에서는 (0, 0, 0) 또는 (1, 1, 1)이며, 각각 거의 같은 확률로 나타납니다.
상태 벡터를 사용하여 양자 시스템을 고전적으로 시뮬레이션하는 것은 사실상 교육용 도구로서만 유용합니다. 이 방식을 쓰면 필요한 자원이 지수적으로 늘어난다는 사실(2^n개의 복소수, 모든 상태에 영향을 미치는 연산)을 직접 마주하게 됩니다. 따라서 작은 양자 시스템조차 시뮬레이션하기에 실용적인 방법이 아닙니다. 이 접근 방식은 큐비트 수가 많아지면 금세 실행할 수 없게 되지만(50~60개의 큐비트를 시뮬레이션하는 것만으로도 현재 슈퍼컴퓨터의 한계에 이릅니다), 양자 컴퓨터가 무엇을 해낼 수 있을지 짐작하게 해 줍니다. 물론 언젠가 양자 컴퓨터를 실제로 만들 수 있다면 말입니다.