Krylov 시간 발전 회로(SKQD) 생성

중요

이 가이드의 개념은 현재 Python API에서만 사용할 수 있습니다. 동등한 기능은 향후 릴리스에서 C API로도 제공될 예정입니다.

샘플 기반 Krylov 양자 대각화(SKQD)는 부분 공간 대각화의 양자 중심 변형입니다. 단일 변분 앤사츠를 샘플링하는 방식(일반적인 샘플 기반 양자 대각화(SQD))이나 무작위 시간 발전 회로의 앙상블을 샘플링하는 방식(SqDRIFT)과 달리, Krylov 기저, 즉 다음 상태들의 집합을 샘플링합니다.

$$\lvert \psi_k \rangle = e^{-i k \, \Delta t \, H} \lvert \psi_0 \rangle, \qquad k = 0, 1, \dots, D - 1,$$

이는 기준 상태 $\lvert \psi_0 \rangle$를 고정된 시간 간격 $\Delta t$의 배수만큼 점점 늘려 가며 발전시켜 얻습니다. 이 $D$개의 회로에서 샘플링한 비트열은 이후 해밀토니안을 고전적으로 대각화하는 Krylov 부분 공간을 채웁니다.

이 가이드는 SKQD 논문에 따라, 비상호작용 배스에 결합된 자성 불순물인 단일 불순물 Anderson 모델(SIAM)에 대해 SKQD 알고리즘의 회로 구성 단계를 재현합니다.

참고

SKQD 샘플링 단계는 바닥 상태가 계산 기저에서 희소한 경우에만 동작합니다. 고전적 대각화가 관리 가능한 수의 샘플링된 비트열로 바닥 상태를 재구성할 수 있어야 하기 때문입니다. 아래에서 보듯이 이러한 희소성은 자동으로 주어지지 않으며, 문제를 표현하는 단일 입자 기저의 속성입니다. 기저를 올바르게 선택하는 것이 이 구성의 핵심입니다.

1. SIAM 해밀토니안을 페르미온 연산자로 구성하기

단일 불순물 Anderson 모델은 상호작용하는 불순물 오비탈 하나(인덱스 $0$)를 사슬 형태로 배열된 $L$개 사이트의 비상호작용 배스와 접촉시킵니다. 이 해밀토니안은 1체 부분 $H_1$과 2체 부분 $H_2$로 나뉩니다.

$$H = H_1 + H_2,$$

1체 부분은 불순물의 온사이트 에너지(화학 퍼텐셜 $\mu$), 불순물-배스 혼성화 $V$, 최근접 이웃 배스 호핑 $t$를 모읍니다.

$$H_1 = \mu \sum_\sigma a^\dagger_{0\sigma} a_{0\sigma} - V \sum_\sigma \left( a^\dagger_{0\sigma} a_{1\sigma} + \text{h.c.} \right) - t \sum_{\sigma} \sum_{i=1}^{L-1} \left( a^\dagger_{i\sigma} a_{i+1,\sigma} + \text{h.c.} \right),$$

2체 부분은 불순물에만 작용하는 온사이트 쿨롱 반발 $U$를 담고 있습니다.

$$H_2 = U \, a^\dagger_{0\uparrow} a_{0\uparrow} a^\dagger_{0\downarrow} a_{0\downarrow}.$$

이 가이드는 입자-정공 대칭점 $\mu = -U/2$에서 진행합니다. 이 두 부분은 1체 적분 h1e(배스 호핑, 불순물 혼성화, 불순물 온사이트 에너지)와 단일 2체 적분 h2e(불순물 쿨롱 항)에 직접 대응됩니다. 먼저 모델이 자연스럽게 정의되는 위치 기저에서 이들을 구성합니다.

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. 운동량 기저로 변경하기

SIAM의 배스 부분은 자유 페르미온 사슬이므로 운동량 기저로 변경하면 대각화됩니다. 이것이 SKQD의 결정적인 단계입니다. 배스는 병진 불변이므로 바닥 상태는 위치 기저에서는 희소하지 않지만, 배스 해밀토니안이 대각인 운동량 기저에서는 희소합니다. 이 변경을 수행하는 단일 입자 오비탈 회전은 문제에 특화된 선형 대수입니다. 이는 불순물 오비탈은 그대로 둔 채 배스 호핑 행렬을 대각화한 다음, 불순물을 중앙 사이트로 옮깁니다.

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

변분 워크플로와 달리, 여기서는 해밀토니안 적분을 운동량 기저로 회전시킵니다. 그 결과 전체 회로(해밀토니안과 상태 모두)가 바닥 상태가 희소한 기저에 놓이게 되며, 샘플링된 비트열이 그만큼 유용한 정보를 담게 됩니다. 1체 및 2체 적분을 회전시키는 것은 전자 구조 텐서의 표준적인 기저 변환입니다.

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

회전은 불순물 오비탈을 고정한 채(배스에만 작용) 이루어지므로 온사이트 상호작용은 변하지 않습니다. h2e_momentum은 h2e와 마찬가지로 단일 수-수 항으로 남습니다. 유일한 변화는 불순물이 중앙 인덱스로 옮겨졌다는 점입니다. 이는 각각 0이 아닌 항목을 하나씩 갖는 적분 텐서에서 바로 확인할 수 있습니다.

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

이러한 국소성 덕분에 나중에 상호작용을 저렴하게 합성할 수 있습니다.

전자 적분 생성자 from_1body_tril_spin_sym()를 사용해 1체 적분을 FermionOperator로 조립합니다. 이 생성자는 하삼각 순서로 패킹된 적분을 기대합니다. 2체 부분은 단일 수-수 항이므로 2체 적분 처리 과정을 거치지 않고 cre()/ann()으로 직접 구성합니다. 즉 불순물 모드 $p$에 대한 온사이트 쿨롱 연산자 $U \, n_{p\uparrow} n_{p\downarrow}$입니다(블록 스핀 배치에서 $\uparrow$ 모드는 인덱스 p, $\downarrow$ 모드는 인덱스 p + norb). 두 생성자 모두 norb개의 공간 오비탈을 스핀 두 배로 늘려 2 * norb개의 페르미온 모드로 만듭니다.

1체 부분 $H_1$과 2체 부분 $H_2$를 미리 합치지 않고 별도의 연산자로 유지하세요. 이는 1절 방정식의 항 $H_1$과 $H_2$에 대응합니다. 각 부분을 (두 기저 모두에서) 한 번만 만들어 두면 아래에서 재사용할 수 있습니다. 이들을 합하면 두 기저 어느 쪽에서든 정확한 대각화에 쓸 전체 해밀토니안이 됩니다. 개별적으로는 $H_2$(운동량 기저)가 4단계에서 Trotter 스텝이 발전시키는 연산자입니다.

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

기저 변경은 유니터리 닮음 변환이므로 스펙트럼은 변하지 않습니다. $(\text{norb}, \text{nelec})$ 섹터에서의 최저 고윳값을 정확한 대각화와 비교해 연산자가 올바른지 확인하세요. FermionOperator는 (네이티브 완전 배치 상호작용(FCI) 행렬-벡터 커널을 기반으로 하는) SciPy LinearOperator를 제공하므로, scipy.sparse.linalg.eigsh()에 바로 전달할 수 있습니다.

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

기저 변경의 이점은 희소성입니다. 이를 구체적으로 보기 위해, 바닥 상태 가중치의 99%를 담는 데 필요한 계산 기저 행렬식의 수를 세어 운동량 기저와 위치 기저를 비교해 보겠습니다. 에너지는 기저와 무관하지만 행렬식의 수는 그렇지 않습니다.

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

두 연산자는 같은 에너지를 갖지만, 위치 기저의 바닥 상태는 수천 개의 행렬식에 퍼져 있는 반면 운동량 기저에서는 약 20개에 집중됩니다. 이것이 SKQD 샘플링이 위치 기저가 아니라 운동량 기저에서 동작하는 이유입니다.

3. 기준 상태 준비하기

SKQD 기준 상태 $\lvert \psi_0 \rangle$는 페르미 준위에 가장 가까운 전자들을 가장 가까운 빈 운동량 모드로 들뜨게 한 모든 들뜸의 중첩입니다. SKQD 논문은 이를 XXPlusYYGate들을 손으로 짠 네트워크로 구성합니다. 이 상태는 강하게 상관된 것처럼 보이지만(0이 아닌 계산 기저 진폭이 수백 개입니다), 수 보존 2모드 회전의 네트워크는 단일 입자(페르미온 가우시안) 연산이므로 그것이 만드는 상태는 단일 Slater 행렬식입니다. 즉 네트워크가 구현하는 단일 입자 회전이 주어지면, 하나의 PrepareSlaterDeterminant 게이트로 이를 준비할 수 있습니다.

각 XXPlusYYGate는 모드 쌍에 대한 Givens 회전으로 동작하므로, 2큐비트 게이트를 직접 내보내는 대신 단일 입자 관점에서(게이트당 Givens 행렬 하나) 그 회전을 바로 만들 수 있습니다. 이 게이트들은 페르미 준위에서 확장되는 층으로 퍼져 나가는 최근접 이웃 회전의 벽돌 쌓기 패턴을 이루며, 각 페르미 준위 근처 전자를 그 바로 위의 빈 모드들에 퍼뜨립니다.

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

참고

위의 회전은 해밀토니안과 맞추어 운동량 기저로 표현됩니다. 페르미 준위 근처의 전자는 에너지가 가장 높은 점유 운동량 모드이며, 네트워크는 이들을 에너지가 가장 낮은 빈 모드로 들뜨게 합니다. 이 때문에 기준 상태는 이전 단계의 기저 변경 이전이 아니라 이후에 구성해야 합니다.

Slater 행렬식 준비는 스핀 섹터별로 수행됩니다. 각 섹터는 고유한 (전자, 오비탈) 형태를 갖기 때문입니다. 하나의 PrepareSlaterDeterminant 게이트(점유, 즉 가장 낮은 nocc개 모드를 채운 상태와 단일 입자 회전의 조합)가 한 섹터를 준비합니다. 여기서는 두 스핀 섹터가 같은 기준 점유와 회전을 공유하므로, 같은 게이트를 알파 모드 0..norb에 한 번, 베타 모드 norb..2*norb에 한 번 적용합니다.

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

이 선언적 게이트가 손으로 작성한 XXPlusYYGate 네트워크를 대체하며, Transpiler가 Givens 회전 합성을 생성합니다.

4. Krylov 회로 조립하기

각 Krylov 회로는 기준 상태 준비 후에 $e^{-i k \, \Delta t \, H}$에 따른 발전이 이어지는 구조입니다. 해밀토니안을 1체 및 2체 부분 $H = H_1 + H_2$로 나누고, 이 분할에 따라 발전을 Trotter화합니다. 이는 회로 깊이에 결정적인 선택입니다. 1체 발전 $e^{-i t H_1}$은 자유 페르미온(페르미온 가우시안) 연산이므로, 항별로 Pauli-Trotter화해야 하는 것이 아니라 단일 입자 기저 회전(유니터리 $e^{-i t \, \mathbf{h}_1}$을 갖는 OrbitalRotation)입니다. 2체 발전 $e^{-i t H_2}$는 단일 온사이트 상호작용 항에만 작용합니다. 이렇게 스텝을 구성하면 Transpiler가 Givens 회전의 얕은 벽돌 쌓기 패턴과 하나의 2큐비트 상호작용 게이트를 내보낼 수 있으며, 이는 SKQD 논문이 손으로 구성하는 얕은 구조입니다.

2체 부분은 이미 가지고 있습니다. 2단계에서 만든 two_body_momentum 연산자(단일 온사이트 항)입니다. 1체 발전 $e^{-i t H_1}$은 스핀에 무관하므로, 유니터리 $e^{-i t \, \mathbf{h}_1}$을 갖는 섹터별 OrbitalRotation 하나를 각 스핀 섹터에 한 번씩 적용하며, 이는 기준 상태 준비와 같은 병렬 구조입니다.

각 Krylov 회로는 총 시간 $k \, \Delta t$ 동안 발전합니다. 이는 크기 $\Delta t$인 $k$개 스텝의 2차 Trotter 곱으로 구현되므로, Krylov 차원이 커져도 스텝당 오차는 고정됩니다. 각 스텝은 (저렴한 단일 항) 2체 발전의 두 반 스텝 사이에 전체 1체 회전(스핀 섹터당 하나)을 끼워 넣습니다. Krylov 차원 $D$는 이러한 회로의 개수이며, dim은 0(기준 상태만)부터 $D - 1$까지의 범위를 갖습니다.

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

참고

여기서 만든 회로에는 아직 측정이 없습니다. 같은 페르미온 회로가 6단계의 정확한 시뮬레이션에도 쓰이는데, 측정은 유니터리 연산이 아니기 때문입니다. 5단계에서 FermionicCircuit.measure_all()로 측정을 추가하며, 이 메서드는 모드 i를 고전 비트 i로 읽습니다.

트랜spile된 회로가 아니라 페르미온 회로를 측정하는 것이 좋습니다. 페르미온 측정은 읽을 모드를 이름으로 지정하므로, 레이아웃과 라우팅이 그 모드가 위치하게 된 큐비트와 함께 이를 따라가고 샘플링된 비트열은 모드 기준으로 인덱싱됩니다. 반대로 트랜spile된 회로를 측정하면 final_index_layout()을 직접 되돌려야 하며, 잊으면 조용히 순서가 바뀐 결과가 반환됩니다. 이는 프리셋이 사용하는 점유 기저 Jordan-Wigner 인코딩, 즉 각 모드가 큐비트 하나로 표현되는 방식에 의존합니다.

5. 큐비트 회로로 트랜spile하기

FermionicCircuit을 트랜spile하면 페르미온 게이트가 큐비트 게이트로 매핑됩니다. 여기에는 Jordan-Wigner jordan_wigner() 매핑이 사용됩니다. OrbitalRotation 게이트는 Givens 회전의 벽돌 쌓기 패턴으로 합성되고, 단일 항 Evolution은 하나의 2큐비트 상호작용 게이트로 낮춰집니다.

여기의 모든 페르미온 게이트에는 기본 합성이 있으므로, 단계를 직접 조립하지 않고 바로 사용할 수 있는 generate_preset_jw_pass_manager()를 사용할 수 있습니다. 키워드 인수는 큐비트 단계로 전달되며, 합성된 깊이를 최적화 없이 있는 그대로 보려면 optimization_level=0을 전달하세요. 그러면 Krylov 회로 전체 집합의 생성은 Krylov 차원에 대한 루프가 됩니다. 트랜spile되지 않은 FermionicCircuit도 트랜spile된 큐비트 회로와 함께 보관하세요. 페르미온 회로는 6단계의 정확한(상태 벡터) 시뮬레이션에 쓰이므로 측정이 없는 상태로 유지해야 합니다. measure_all(inplace=False)는 트랜spile할 측정된 복사본을 반환하고 원본은 그대로 두며, 백엔드가 실행하는 것은 이 복사본입니다.

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         measured = fermionic_circuit.measure_all(inplace=False)
...         circuits.append(pass_manager.run(measured))
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

Krylov 차원은 회로의 2큐비트 깊이에 그대로 나타납니다. 기준 상태 준비는 고정 비용이며, Krylov 거듭제곱이 하나 늘 때마다 시간 발전 연산자의 2차 Trotter 스텝이 하나씩 추가됩니다.

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

참고

해밀토니안을 나누면 이러한 깊이가 적당한 수준으로 유지됩니다. 1체 발전은 진정한 오비탈 회전이므로 인접한 큐비트에 대한 XXPlusYYGate 벽돌 쌓기 패턴으로 합성됩니다. 온사이트 상호작용의 각 반 스텝은 불순물의 두 스핀 모드를 결합하는 단일 RZZGate이며, 이는 2차 곱이 회전 양쪽에 끼워 넣는 저렴한 항입니다. 대신 전체 해밀토니안을 하나의 Evolution 게이트에 넘기면 1체 부분을 항별로 Pauli-Trotter화하게 되어, (장거리) 불순물-배스 결합에 걸쳐 긴 Jordan-Wigner Z 문자열이 전파되고 깊이가 자릿수 하나만큼 부풀어 오릅니다.

첫 번째 회로는 기준 행렬식만 준비하고, 마지막 회로는 가장 많은 Trotter 스텝(여기서는 네 개)을 포함합니다.

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. Krylov 회로에서 비트열 샘플링하기

하드웨어에서는 트랜spile된 회로를 실행하고 측정합니다. 여기서는 그 단계의 잡음 없는(상태 벡터) 버전으로, SKQD가 수집할 비트열을 보여 줍니다. 트랜spile되지 않은 각 FermionicCircuit을 ffsim.apply_unitary()로 시뮬레이션하고, 결과 상태 벡터를 ffsim.sample_state_vector()로 샘플링합니다.

PrepareSlaterDeterminant 게이트는 시뮬레이션에서 검증 후 회전 방식으로 동작합니다(입력 상태가 기준 행렬식을 점유하는지 확인한 다음 회전을 적용합니다). 따라서 입력은 회전이 없는 단순한 점유 행렬식, 즉 스핀별로 가장 낮은 nocc개 모드를 채운 상태입니다. 게이트가 reference_rotation을 직접 적용합니다.

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

다섯 개의 Krylov 회로 모두에서 얻은 샘플을 합치면 지지 집합은 매우 작아서, (4, 4) 섹터의 $\binom{8}{4}^2 = 4900$개 행렬식 중 서로 다른 비트열이 수백 개에 불과합니다. 2단계에서 확인한 운동량 기저의 희소성에서 직접 비롯된 이러한 집중 덕분에 이후의 고전적 대각화를 다룰 수 있게 됩니다.

qiskit.visualization.plot_histogram()으로 카운트를 그려 보면 희소성을 확인할 수 있습니다. 기준 행렬식을 비롯해 소수의 구성이 지배적입니다.

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. 샘플링된 부분공간에서 대각화하기

샘플링된 비트열은 SKQD의 고전 연산 부분에 대한 입력이 됩니다. 해밀토니안을 샘플링된 구성들이 펼치는 부분공간에 사영한 다음 그 안에서 대각화합니다. 이 작업은 qiskit-addon-sqd 패키지에 맡기며, 이 패키지의 diagonalize_fermionic_hamiltonian()이 샘플 기반 양자 대각화 루프 전체를 실행합니다. 이 함수는 샘플링된 구성들의 배치로 부분공간을 만들고, 그 안에서 해밀토니안을 대각화하며, 구성 복구(configuration recovery)를 사용해 부분공간을 반복적으로 다듬습니다. 구성 복구란 알려진 입자 수 대칭성을 위반하는 구성을 바로잡기 위해 점유 상태를 뒤집는 것입니다. 구성 복구는 하드웨어 잡음이 샘플링된 비트열을 잘못된 입자 수 섹터로 오염시키는 문제를 되돌리도록 설계되었습니다. 여기서의 샘플링은 잡음이 없으므로 그 측면에서는 복구할 것이 없습니다.

솔버는 샘플링된 카운트를 BitArray로 받으며, 이는 6단계에서 모은 counts로 직접 만듭니다. num_bits는 명시적으로 전달하세요. 전달하지 않으면 너비가 가장 넓은 샘플링된 비트열에서 추론되는데, 마지막 모드를 점유한 샘플이 하나도 없으면 이 너비가 2 * norb보다 좁아지고, 그러면 솔버가 배열을 거부합니다. 솔버는 나머지 구성에서 쓴 것과 같은 운동량 기저에서 대각화하므로, 2단계의 h1e_momentum과 h2e_momentum 적분을 다시 정리하지 않고 그대로 받습니다. 복구된 에너지는 단 몇백 개의 샘플링된 구성만으로도 정확한 기준값과 밀리 하트리보다 훨씬 작은 오차로 일치합니다. SQD 에너지는 실제 바닥 상태에 대한 변분 상한이므로, 항상 reference_energy(2단계에서 정의된 값이 아직 유효합니다)와 같거나 그보다 높습니다.

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts), num_bits=2 * norb)
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

실행 과정을 시각화합니다. 반복에 따른 에너지의 수렴과, 복구된 바닥 상태에서 각 공간 오비탈의 평균 점유도를 그립니다. 반복별 에너지는 해당 반복의 배치들 중 가장 낮은 값이고, 점유도는 두 스핀 섹터에 걸친 orbital_occupancies를 합한 것입니다. 보정할 잡음이 없는데도 에너지는 반복할수록 감소합니다. 각 반복은 풀에서 부분 샘플링한 samples_per_batch개의 구성으로 이루어진 num_batches개의 배치 안에서 대각화하며, 비트열 이월(carryover)에 의한 반복적 개선으로 어떤 구성이 그 배치에 들어가는지가 꾸준히 좋아지기 때문에 복구된 에너지가 개선됩니다.

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

이것으로 SKQD 루프가 마무리됩니다. 1~5단계의 운동량 기저 회로가 희소한 샘플(6단계)을 만들고, 위의 샘플 기반 대각화가 그 샘플을 다시 바닥 상태 에너지로 바꿔 줍니다.

다음 단계

이 가이드는 잡음 없는 상태 벡터 시뮬레이터에서 SKQD 파이프라인 전체를 실행했습니다. 하드웨어에서는 6단계의 상태 벡터 샘플링 대신 트랜스파일된 circuits를 실행하고 측정하며, 그 결과 카운트는 7단계의 대각화에 그대로 입력됩니다. 회로 실행에 대한 도움말은 Qiskit 문서를, 잡음이 있는 샘플에서의 동작을 포함해 구성 복구가 가장 큰 역할을 하는 부분공간 대각화 후처리에 대한 정보는 SQD 애드온 튜토리얼을 참고하세요.

ffsim.apply_unitary()와 ffsim.sample_state_vector()가 이 패키지의 페르미온 회로에서 왜 기본적으로 동작하는지, 그리고 6단계에서 사용한 고정 입자 수 섹터가 2단계에서 확립한 운동량 기저의 희소성을 왜 샘플 효율적으로 만드는지 이해하려면 ffsim 관계 가이드를 읽어 보세요.