{ "cells": [ { "cell_type": "markdown", "id": "ee388de5-2507-48ff-97a6-6777698d6256", "metadata": {}, "source": [ "---\n", "title: \"페르미온 격자 모델의 샘플 기반 Krylov 양자 대각화\"\n", "description: \"샘플 기반 양자 대각화 알고리즘을 사용하여 잡음이 있는 양자 하드웨어로 단일 불순물 앤더슨 모델을 시뮬레이션한다.\"\n", "---\n", "\n", "{/* cspell:ignore fontdict fontsize nocc SQKD DMRG textrm varepsilon vecs pqrs ijkl */}\n", "\n", "\n", "\n", "# 페르미온 격자 모델의 샘플 기반 Krylov 양자 대각화\n", "\n", "*사용량 추정치: Heron r2 프로세서에서 9초(참고: 이는 추정치일 뿐입니다. 런타임은 다를 수 있습니다.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "a1b2c3d4-e5f6-7890-abcd-ef1234567890", "metadata": {}, "source": [ "\n", "\n", "## 학습 성과\n", "\n", "* [SQD Qiskit 애드온](https://github.com/Qiskit/qiskit-addon-sqd) 을 사용하여 양자 처리 장치(QPU)에서 샘플링한 비트열을 바탕으로 격자 모델의 기저 상태 에너지를 근사하는 방법.\n", "* [ffsim을](https://github.com/qiskit-community/ffsim) 사용하여 페르미온 시뮬레이션을 위한 시간 발전 회로를 구축하는 방법.\n", "* 샘플 기반 크릴로프 대각화(SKQD) 알고리즘을 사용하여 후처리를 위해 여러 회로의 샘플을 결합하는 방법.\n", "\n", "\n", "\n", "## 전제조건\n", "\n", "* [화학 Hamiltonian의 샘플 기반 양자 대각화](/docs/tutorials/sample-based-quantum-diagonalization)\n", "* [격자 해밀토니안의 크릴로프 양자 대각화](/docs/tutorials/krylov-quantum-diagonalization)\n", "* [Qiskit primitives](/docs/guides/primitives)\n", "\n" ] }, { "cell_type": "markdown", "id": "dc5cc74e-06bf-45ac-a69b-81778138e08f", "metadata": {}, "source": [ "\n", "\n", "## 배경\n", "\n", "이 튜토리얼에서는 샘플 기반 양자 대각선화(SQD)를 사용하여 페르미온 격자 모델의 기저 상태 에너지를 추정하는 방법을 설명합니다. 특히 금속에 포함된 자성 불순물을 설명하는 데 사용되는 1차원 단일 불순물 앤더슨 모델(SIAM)을 연구합니다.\n", "\n", "이 튜토리얼은 관련 튜토리얼인 [화학 해밀턴의 샘플 기반 양자 대각선 화](/docs/tutorials/sample-based-quantum-diagonalization) 와 유사한 워크플로우를 따릅니다. 그러나 중요한 차이점은 양자 회로가 구축되는 방식에 있습니다. 다른 튜토리얼은 잠재적으로 수백만 개의 상호 작용 조건이 있는 화학 해밀턴에게 매력적인 휴리스틱 변형 안사츠를 사용합니다. 반면에 이 튜토리얼에서는 해밀턴의 시간 진화에 근사한 회로를 사용합니다. 이러한 회로는 깊이가 깊을 수 있으므로 격자 모델에 적용하는 데 이 접근 방식이 더 적합합니다. 이러한 회로에 의해 준비된 상태 벡터는 [크릴로프 부분 공간의](https://en.wikipedia.org/wiki/Krylov_subspace) 기초를 형성하며, 결과적으로 알고리즘은 적절한 가정 하에서 증명 가능하고 효율적으로 기저 상태로 수렴합니다.\n", "\n", "이 튜토리얼에서 사용된 접근 방식은 SQD와 [크릴로프 양자 대각선화(KQD)](https://arxiv.org/abs/2407.14431) 에 사용된 기술을 조합한 것으로 볼 수 있습니다. 이 결합된 접근 방식을 샘플 기반 크릴로프 양자 대각선화(SQKD)라고도 합니다. KQD 방법에 대한 튜토리얼은 [격자 해밀턴의 크릴로프 양자 대각선화를](/docs/tutorials/krylov-quantum-diagonalization) 참조하세요.\n", "\n", "이 튜토리얼은 [\"샘플 기반 크릴로프 대각선화를 위한 양자 중심 알고리즘\"](https://arxiv.org/abs/2501.09702) 을 기반으로 하며, 자세한 내용은 이 튜토리얼을 참조할 수 있습니다.\n", "\n", "\n", "\n", "### 단일 불순물 앤더슨 모델 (SIAM)\n", "\n", "1차원 SIAM 해밀턴은 세 개의 항의 합입니다:\n", "\n", "$$\n", "H = H_{\\textrm{imp}}+ H_\\textrm{bath} + H_\\textrm{hyb},\n", "$$\n", "\n", "여기서,\n", "\n", "$$\n", "\\begin{align*}\n", " H_\\textrm{imp} &= \\varepsilon \\left( \\hat{n}_{d\\uparrow} + \\hat{n}_{d\\downarrow} \\right) + U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}, \\\\\n", " H_\\textrm{bath} &= -t \\sum_{\\substack{\\mathbf{j} = 0\\\\ \\sigma\\in \\{\\uparrow, \\downarrow\\}}}^{L-1} \\left(\\hat{c}^\\dagger_{\\mathbf{j}, \\sigma}\\hat{c}_{\\mathbf{j}+1, \\sigma} + \\hat{c}^\\dagger_{\\mathbf{j}+1, \\sigma}\\hat{c}_{\\mathbf{j}, \\sigma} \\right), \\\\\n", " H_\\textrm{hyb} &= V\\sum_{\\sigma \\in \\{\\uparrow, \\downarrow \\}} \\left(\\hat{d}^\\dagger_\\sigma \\hat{c}_{0, \\sigma} + \\hat{c}^\\dagger_{0, \\sigma} \\hat{d}_{\\sigma} \\right).\n", "\\end{align*}\n", "$$\n", "\n", "여기서 $c^\\dagger_{\\mathbf{j},\\sigma}/c_{\\mathbf{j},\\sigma}$ 은 스핀이 있는 $\\mathbf{j}^{\\textrm{th}}$ 배스 사이트의 페르미온 생성/소멸 연산자 $\\sigma$, $\\hat{d}^\\dagger_{\\sigma}/\\hat{d}_{\\sigma}$ 은 불순물 모드의 생성/소멸 연산자이며, $\\hat{n}_{d\\sigma} = \\hat{d}^\\dagger_{\\sigma} \\hat{d}_{\\sigma}$, $t$, $U$, $V$ 은 호핑, 온사이트, 혼성화 상호작용을 나타내는 실수이고, $\\varepsilon$ 은 화학 전위를 명시하는 실수입니다.\n", "\n", "해밀턴은 일반적인 상호작용 전자 해밀턴의 특정 인스턴스라는 점에 유의하세요,\n", "\n", "$$\n", "\\begin{align*}\n", " H &= \\sum_{\\substack{p, q \\\\ \\sigma}} h_{pq} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{q\\sigma} + \\sum_{\\substack{p, q, r, s \\\\ \\sigma \\tau}} \\frac{h_{pqrs}}{2} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}^\\dagger_{q\\tau} \\hat{a}_{s\\tau} \\hat{a}_{r\\sigma} \\\\\n", " &= H_1 + H_2,\n", "\\end{align*}\n", "$$\n", "\n", "여기서 $H_1$ 은 페르미온 생성 및 소멸 연산자에서 2진법인 1진법으로 구성되고 $H_2$ 은 2진법인 2진법으로 구성됩니다. SIAM의 경우,\n", "\n", "$$\n", "H_2 = U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}\n", "$$\n", "\n", "와 $H_1$ 에는 해밀턴의 나머지 용어가 포함되어 있습니다. 해밀턴을 프로그래밍 방식으로 표현하기 위해 행렬 $h_{pq}$ 과 텐서 $h_{pqrs}$ 를 저장합니다.\n", "\n", "\n", "\n", "### 위치와 운동량 기반\n", "\n", "$H_\\textrm{bath}$ 의 대략적인 병진 대칭으로 인해 위치 기준(위에서 해밀턴이 지정된 궤도 기준)에서 지상 상태가 희박할 것으로 예상하지 않습니다. SQD의 성능은 기저 상태가 희소할 때, 즉 소수의 계산 기준 상태에만 상당한 가중치가 있는 경우에만 보장됩니다. 지상 상태의 희소성을 개선하기 위해 $H_\\textrm{bath}$ 이 대각선인 궤도 기준으로 시뮬레이션을 수행합니다. 우리는 이 기준을 *모멘텀 기준이라고* 부릅니다. $H_\\textrm{bath}$ 은 이차 페르미온 해밀턴이므로 궤도 회전에 의해 효율적으로 대각선화할 수 있습니다.\n", "\n", "\n", "\n", "### 해밀토니언에 의한 근사적 시간 진화\n", "\n", "해밀턴에 의한 시간 진화의 근사치를 구하기 위해 2 차 트로터-스즈키 분해를 사용합니다,\n", "\n", "$$\n", " e^{-i \\Delta t H} \\approx e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2}.\n", "$$\n", "\n", "[조던-위그너 변환에서](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation) $H_2$ 에 의한 시간 진화는 불순물 부위에서 스핀업 궤도와 스핀다운 궤도 사이의 단일 [CPhase](/docs/api/qiskit/qiskit.circuit.library.CPhaseGate) 게이트에 해당합니다. $H_1$ 은 이차 페르미온 해밀턴이므로 $H_1$ 에 의한 시간 진화는 궤도 회전에 해당합니다.\n", "\n", "크릴로프 기저는 $\\{ |\\psi_k\\rangle \\}_{k=0}^{D-1}$, 여기서 $D$ 은 크릴로프 부분 공간의 차원이며, 단일 트로터 단계의 반복 적용으로 형성됩니다\n", "\n", "$$\n", " |\\psi_k\\rangle \\approx \\left[e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2} \\right]^k\\ket{\\psi_0}.\n", "$$\n", "\n", "다음 SQD 기반 워크플로에서는 이 회로 세트에서 샘플링하고 결합된 비트스트링 세트를 SQD로 후처리합니다. 이 접근법은 관련 튜토리얼인 [화학 해밀턴의 샘플 기반 양자 대각선화에서](/docs/tutorials/sample-based-quantum-diagonalization) 사용된 것과는 대조적으로, 단일 휴리스틱 변형 회로에서 샘플을 도출했습니다.\n", "\n" ] }, { "cell_type": "markdown", "id": "df66d697-102c-4a2c-80f3-f67fdda05573", "metadata": {}, "source": [ "\n", "\n", "## 요구사항\n", "\n", "이 튜토리얼을 시작하기 전에 다음이 설치되어 있는지 확인하세요:\n", "\n", "* Qiskit SDK v1.0 또는 이후 버전, [시각화](/docs/api/qiskit/visualization) 지원 기능 포함\n", "* Qiskit Runtime v0.22 또는 이후 (`pip install qiskit-ibm-runtime`)\n", "* SQD Qiskit 애드온 v0.11 이상 (`pip install qiskit-addon-sqd`)\n", "* ffsim v0.0.72 이상 (`pip install ffsim`)\n", "\n" ] }, { "cell_type": "markdown", "id": "c3d4e5f6-a7b8-9012-cdef-123456789012", "metadata": {}, "source": [ "\n", "\n", "## 소규모 시뮬레이터 예시\n", "\n" ] }, { "cell_type": "markdown", "id": "8540487a-8033-49c2-9f30-022336105f64", "metadata": {}, "source": [ "\n", "\n", "### 1단계: 문제를 양자 회로에 매핑하기\n", "\n" ] }, { "cell_type": "markdown", "id": "4e6652bd-c97d-4f2f-98a7-d44857805cf1", "metadata": {}, "source": [ "먼저 위치 기준으로 SIAM 해밀턴을 생성합니다. 해밀턴은 행렬 $h_{pq}$ 과 텐서 $h_{pqrs}$ 로 표현됩니다. 그런 다음 이를 운동량 기준으로 회전합니다. 위치 기준에서는 첫 번째 위치에 불순물을 배치합니다. 그러나 모멘텀 기준으로 회전할 때는 다른 궤도와의 상호작용을 촉진하기 위해 불순물을 중앙 위치로 이동합니다.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "b5cb9c28-4721-4141-8665-96885038e210", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pyscf.fci\n", "\n", "\n", "def siam_hamiltonian(\n", " norb: int,\n", " hopping: float,\n", " onsite: float,\n", " hybridization: float,\n", " chemical_potential: float,\n", ") -> tuple[np.ndarray, np.ndarray]:\n", " \"\"\"Hamiltonian for the single-impurity Anderson model.\"\"\"\n", " # Place the impurity on the first site\n", " impurity_orb = 0\n", "\n", " # One body matrix elements in the \"position\" basis\n", " h1e = np.zeros((norb, norb))\n", " np.fill_diagonal(h1e[:, 1:], -hopping)\n", " np.fill_diagonal(h1e[1:, :], -hopping)\n", " h1e[impurity_orb, impurity_orb + 1] = -hybridization\n", " h1e[impurity_orb + 1, impurity_orb] = -hybridization\n", " h1e[impurity_orb, impurity_orb] = chemical_potential\n", "\n", " # Two body matrix elements in the \"position\" basis\n", " h2e = np.zeros((norb, norb, norb, norb))\n", " h2e[impurity_orb, impurity_orb, impurity_orb, impurity_orb] = onsite\n", "\n", " return h1e, h2e\n", "\n", "\n", "def momentum_basis(norb: int) -> np.ndarray:\n", " \"\"\"Get the orbital rotation to change from the position to the momentum basis.\"\"\"\n", " n_bath = norb - 1\n", "\n", " # Orbital rotation that diagonalizes the bath (non-interacting system)\n", " hopping_matrix = np.zeros((n_bath, n_bath))\n", " np.fill_diagonal(hopping_matrix[:, 1:], -1)\n", " np.fill_diagonal(hopping_matrix[1:, :], -1)\n", " _, vecs = np.linalg.eigh(hopping_matrix)\n", "\n", " # Expand to include impurity\n", " orbital_rotation = np.zeros((norb, norb))\n", " # Impurity is on the first site\n", " orbital_rotation[0, 0] = 1\n", " orbital_rotation[1:, 1:] = vecs\n", "\n", " # Move the impurity to the center\n", " new_index = n_bath // 2\n", " perm = np.r_[1 : (new_index + 1), 0, (new_index + 1) : norb]\n", " orbital_rotation = orbital_rotation[:, perm]\n", "\n", " return orbital_rotation\n", "\n", "\n", "def rotated(\n", " h1e: np.ndarray, h2e: np.ndarray, orbital_rotation: np.ndarray\n", ") -> tuple[np.ndarray, np.ndarray]:\n", " \"\"\"Rotate the orbital basis of a Hamiltonian.\"\"\"\n", " h1e_rotated = np.einsum(\n", " \"ab,Aa,Bb->AB\",\n", " h1e,\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " optimize=\"greedy\",\n", " )\n", " h2e_rotated = np.einsum(\n", " \"abcd,Aa,Bb,Cc,Dd->ABCD\",\n", " h2e,\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " optimize=\"greedy\",\n", " )\n", " return h1e_rotated, h2e_rotated\n", "\n", "\n", "# Total number of spatial orbitals, including the bath sites and the impurity\n", "# This should be an even number\n", "norb = 8\n", "\n", "# System is half-filled\n", "nelec = (norb // 2, norb // 2)\n", "# One orbital is the impurity, the rest are bath sites\n", "n_bath = norb - 1\n", "\n", "# Hamiltonian parameters\n", "hybridization = 1.0\n", "hopping = 1.0\n", "onsite = 10.0\n", "chemical_potential = -0.5 * onsite\n", "\n", "# Generate Hamiltonian in position basis\n", "h1e, h2e = siam_hamiltonian(\n", " norb=norb,\n", " hopping=hopping,\n", " onsite=onsite,\n", " hybridization=hybridization,\n", " chemical_potential=chemical_potential,\n", ")\n", "\n", "# Rotate to momentum basis\n", "orbital_rotation = momentum_basis(norb)\n", "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n", "# In the momentum basis, the impurity is placed in the center\n", "impurity_index = n_bath // 2\n", "\n", "# Use PySCF to compute the exact ground state energy\n", "reference_energy, _ = pyscf.fci.direct_spin1.kernel(h1e, h2e, norb, nelec)" ] }, { "cell_type": "markdown", "id": "e888edf2-7865-41a8-be57-6ddb72dd0cc7", "metadata": {}, "source": [ "다음으로 크릴로프 기저 상태를 생성하기 위해 회로를 생성합니다.\n", "각 스핀 종에 대해 초기 상태 $\\ket{\\psi_0}$ 는 페르미 레벨에 가장 가까운 세 전자의 가능한 모든 여기들을 $|00\\cdots 0011 \\cdots 11\\rangle$ 상태에서 시작하여 가장 가까운 4개의 빈 모드로 중첩하여 주어지며, 7개의 [XXPlusYY게이트를](/docs/api/qiskit/qiskit.circuit.library.XXPlusYYGate) 적용하여 실현됩니다.\n", "시간 진화 상태는 2차 트로터 단계를 연속적으로 적용하여 생성됩니다.\n", "\n", "이 모델과 회로 설계 방식에 대한 자세한 설명은 [\"샘플 기반 크릴로프 대각선화를 위한 양자 중심 알고리즘\"](https://arxiv.org/abs/2501.09702) 을 참조하세요.\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "0f729f86-1814-4d5a-ae65-3f9614e103b3", "metadata": {}, "outputs": [], "source": [ "from typing import Sequence\n", "\n", "import ffsim\n", "import scipy\n", "from qiskit import QuantumCircuit, QuantumRegister\n", "from qiskit.circuit import CircuitInstruction, Qubit\n", "from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate\n", "\n", "\n", "def prepare_initial_state(qubits: Sequence[Qubit], norb: int, nocc: int):\n", " \"\"\"Prepare initial state.\"\"\"\n", " assert norb >= 8\n", " x_gate = XGate()\n", " rot = XXPlusYYGate(0.5 * np.pi, -0.5 * np.pi)\n", " for i in range(nocc):\n", " yield CircuitInstruction(x_gate, [qubits[i]])\n", " yield CircuitInstruction(x_gate, [qubits[norb + i]])\n", " for i in range(3):\n", " for j in range(nocc - i - 1, nocc + i, 2):\n", " yield CircuitInstruction(rot, [qubits[j], qubits[j + 1]])\n", " yield CircuitInstruction(\n", " rot, [qubits[norb + j], qubits[norb + j + 1]]\n", " )\n", " yield CircuitInstruction(rot, [qubits[j + 1], qubits[j + 2]])\n", " yield CircuitInstruction(\n", " rot, [qubits[norb + j + 1], qubits[norb + j + 2]]\n", " )\n", "\n", "\n", "def trotter_step(\n", " qubits: Sequence[Qubit],\n", " time_step: float,\n", " one_body_evolution: np.ndarray,\n", " h2e: np.ndarray,\n", " impurity_index: int,\n", " norb: int,\n", "):\n", " \"\"\"A Trotter step.\"\"\"\n", " # Assume the two-body interaction is just the on-site interaction of the impurity\n", " onsite = h2e[\n", " impurity_index, impurity_index, impurity_index, impurity_index\n", " ]\n", " # Two-body evolution for half the time\n", " yield CircuitInstruction(\n", " CPhaseGate(-0.5 * time_step * onsite),\n", " [qubits[impurity_index], qubits[norb + impurity_index]],\n", " )\n", " # One-body evolution for the full time\n", " yield CircuitInstruction(\n", " ffsim.qiskit.OrbitalRotationJW(norb, one_body_evolution), qubits\n", " )\n", " # Two-body evolution for half the time\n", " yield CircuitInstruction(\n", " CPhaseGate(-0.5 * time_step * onsite),\n", " [qubits[impurity_index], qubits[norb + impurity_index]],\n", " )\n", "\n", "\n", "# Time step\n", "time_step = 0.2\n", "# Number of Krylov basis states\n", "krylov_dim = 8\n", "\n", "# Initialize circuit\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "\n", "# Generate initial state\n", "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n", " circuit.append(instruction)\n", "circuit.measure_all()\n", "\n", "# Create list of circuits, starting with the initial state circuit\n", "circuits = [circuit.copy()]\n", "\n", "# Add time evolution circuits to the list\n", "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n", "for i in range(krylov_dim - 1):\n", " # Remove measurements\n", " circuit.remove_final_measurements()\n", " # Append another Trotter step\n", " for instruction in trotter_step(\n", " qubits,\n", " time_step,\n", " one_body_evolution,\n", " h2e_momentum,\n", " impurity_index,\n", " norb,\n", " ):\n", " circuit.append(instruction)\n", " # Measure qubits\n", " circuit.measure_all()\n", " # Add a copy of the circuit to the list\n", " circuits.append(circuit.copy())" ] }, { "cell_type": "code", "execution_count": 3, "id": "9f2cc4d4-ecac-457a-bcae-558319668e1f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "circuits[0].draw(\"mpl\", scale=0.4, fold=-1)" ] }, { "cell_type": "code", "execution_count": 4, "id": "827976ec-4815-4707-80b1-e13fb2fef309", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "circuits[-1].draw(\"mpl\", scale=0.4, fold=-1)" ] }, { "cell_type": "markdown", "id": "b8ca6be4-61d9-47be-8099-8712c7ecc774", "metadata": {}, "source": [ "\n", "\n", "### 2단계: 양자 실행을 위한 문제 최적화\n", "\n" ] }, { "cell_type": "markdown", "id": "a3304e1b-9c7c-4212-8744-d1c62292eced", "metadata": {}, "source": [ "다음으로, 대상 하드웨어에 맞게 회로를 최적화합니다. 우선, 지정된 큐비트 수와 시간 진화 회로가 자연스럽게 분해되는 게이트 집합을 갖춘 일반적인 백엔드를 만들어 보겠습니다.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "2d2fdbff-1e22-45af-a2eb-c334e4328c59", "metadata": {}, "outputs": [], "source": [ "from qiskit.providers.fake_provider import GenericBackendV2\n", "\n", "backend = GenericBackendV2(\n", " 2 * norb, basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\"]\n", ")" ] }, { "cell_type": "markdown", "id": "ce3e6a52-7b99-49b6-8294-b005efa59cfc", "metadata": {}, "source": [ "이제 키스킷을 사용하여 대상 백엔드로 회로를 트랜스파일합니다.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "c8643533-9fec-40bf-a307-da8839b1e444", "metadata": {}, "outputs": [], "source": [ "from qiskit.transpiler import generate_preset_pass_manager\n", "\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend\n", ")\n", "isa_circuits = pass_manager.run(circuits)" ] }, { "cell_type": "markdown", "id": "6cfd3eea-e2d9-40a5-a449-1d3d790a5f2d", "metadata": {}, "source": [ "\n", "\n", "### 3단계: `Qiskit primitives`를 사용하여 실행합니다\n", "\n" ] }, { "cell_type": "markdown", "id": "ad48e3e6-1013-45a1-942c-63766fed5819", "metadata": {}, "source": [ "하드웨어 실행을 위해 회로를 최적화한 후, 대상 하드웨어에서 실행하고 접지 상태 에너지 추정을 위한 샘플을 수집할 준비가 되었습니다. 샘플러 프리미티브를 사용하여 각 회로에서 비트스트링을 샘플링한 후, 모든 결과를 단일 카운트 사전으로 결합하고 가장 일반적으로 샘플링된 상위 20개의 비트스트링을 플로팅합니다.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "80eee553-60d6-4258-88ab-d8d120418c36", "metadata": {}, "outputs": [], "source": [ "from qiskit.visualization import plot_histogram\n", "from qiskit.primitives import StatevectorSampler\n", "\n", "# Sample from the circuits\n", "sampler = StatevectorSampler()\n", "job = sampler.run(isa_circuits, shots=500)" ] }, { "cell_type": "code", "execution_count": 8, "id": "10af4663-7375-4b50-bae6-9f3d5106457b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qiskit.primitives import BitArray\n", "\n", "# Combine the shots from the individual Trotter circuits\n", "bit_array = BitArray.concatenate_shots(\n", " [result.data.meas for result in job.result()]\n", ")\n", "\n", "plot_histogram(bit_array.get_counts(), number_to_keep=20)" ] }, { "cell_type": "markdown", "id": "2aa74455-d16b-4ac3-a354-54a79d5c5759", "metadata": {}, "source": [ "\n", "\n", "### 4단계: 후처리 및 원하는 클래식 형식으로 결과 반환\n", "\n" ] }, { "cell_type": "markdown", "id": "d3f7713c-a4e6-407f-94de-0e5cfeb1134c", "metadata": {}, "source": [ "이제 `diagonalize_fermionic_hamiltonian` 함수를 사용하여 SQD 알고리즘을 실행합니다. 이 함수의 인수에 대한 설명은 [API 설명서를](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) 참조하세요.\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "7609d1e1-e8ef-48e1-a965-97927f403163", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -13.4222953188441\n", "\t\tSubspace dimension: 529\n", "\tSubsample 1\n", "\t\tEnergy: -13.42237556285828\n", "\t\tSubspace dimension: 784\n", "\tSubsample 2\n", "\t\tEnergy: -13.422045397387413\n", "\t\tSubspace dimension: 529\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -13.422379583305478\n", "\t\tSubspace dimension: 900\n", "\tSubsample 1\n", "\t\tEnergy: -13.422376197704326\n", "\t\tSubspace dimension: 841\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421162849295\n", "\t\tSubspace dimension: 1089\n", "Iteration 3\n", "\tSubsample 0\n", "\t\tEnergy: -13.422421164670345\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 1\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421205869572\n", "\t\tSubspace dimension: 1156\n", "Iteration 4\n", "\tSubsample 0\n", "\t\tEnergy: -13.422421494558726\n", "\t\tSubspace dimension: 1225\n", "\tSubsample 1\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n" ] } ], "source": [ "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", ")\n", "\n", "# List to capture intermediate results\n", "result_history = []\n", "\n", "\n", "def callback(results: list[SCIResult]):\n", " result_history.append(results)\n", " iteration = len(result_history)\n", " print(f\"Iteration {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\"\\tSubsample {i}\")\n", " print(f\"\\t\\tEnergy: {result.energy}\")\n", " print(\n", " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", "\n", "rng = np.random.default_rng(24)\n", "result = diagonalize_fermionic_hamiltonian(\n", " h1e_momentum,\n", " h2e_momentum,\n", " bit_array,\n", " samples_per_batch=100,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=3,\n", " max_iterations=5,\n", " symmetrize_spin=True,\n", " callback=callback,\n", " seed=rng,\n", ")" ] }, { "cell_type": "markdown", "id": "3dee9c61-fc42-48e9-8888-af8fc831cd5c", "metadata": {}, "source": [ "다음 코드 셀은 결과를 그래프로 표시합니다. 첫 번째 그래프는 구성 복원 반복 횟수에 따른 계산된 에너지를 나타내며, 두 번째 그래프는 최종 반복 후 각 공간 궤도의 평균 점유율을 보여줍니다. 이 문제가 매우 간단하기 때문에, 첫 번째 반복만으로도 정확한 에너지 값에 매우 근접하게 됩니다(y축의 눈금을 참고하세요).\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "b6879566-8bf5-4c28-bfb6-b2686692e3d3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reference energy: -13.42249\n", "SQD energy: -13.42242\n", "Absolute error: 0.00007\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "min_es = [\n", " min(result, key=lambda res: res.energy).energy\n", " for result in result_history\n", "]\n", "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n", "\n", "# Data for energies plot\n", "x1 = range(len(result_history))\n", "\n", "# Data for avg spatial orbital occupancy\n", "y2 = np.sum(result.orbital_occupancies, axis=0)\n", "x2 = range(len(y2))\n", "\n", "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "\n", "# Plot energies\n", "axs[0].plot(x1, min_es, label=\"energy\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].axhline(\n", " y=reference_energy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"reference energy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy\", fontdict={\"fontsize\": 12})\n", "axs[0].legend()\n", "\n", "# Plot orbital occupancy\n", "axs[1].bar(x2, y2, width=0.8)\n", "axs[1].set_xticks(x2)\n", "axs[1].set_xticklabels(x2)\n", "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n", "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n", "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n", "\n", "print(f\"Reference energy: {reference_energy:.5f}\")\n", "print(f\"SQD energy: {min_e:.5f}\")\n", "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0c9f4976-d770-426a-822e-e6756c2cfbe7", "metadata": {}, "source": [ "\n", "\n", "### 에너지 확인\n", "\n", "SQD가 반환하는 에너지는 실제 기저 상태 에너지의 상한값임을 보장합니다. SQD는 기저 상태를 근사하는 상태 벡터의 계수도 함께 반환하므로, 에너지 값을 확인할 수 있습니다. 다음 코드 셀에서 보여주는 것처럼, 1입자 및 2입자 축소 밀도 행렬을 사용하여 상태 벡터로부터 에너지를 계산할 수 있습니다.\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "e2b9de72-61cf-49d3-a1b5-f043e4b16956", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Recomputed energy: -13.42242\n" ] } ], "source": [ "rdm1 = result.sci_state.rdm(rank=1, spin_summed=True)\n", "rdm2 = result.sci_state.rdm(rank=2, spin_summed=True)\n", "\n", "energy = np.sum(h1e_momentum * rdm1) + 0.5 * np.sum(h2e_momentum * rdm2)\n", "\n", "print(f\"Recomputed energy: {energy:.5f}\")" ] }, { "cell_type": "markdown", "id": "5221928a-79ff-4f54-90b4-1fe4d7739aae", "metadata": {}, "source": [ "\n", "\n", "## 대규모 하드웨어 예시\n", "\n", "이제 실제 QPU에서 더 큰 규모의 예제를 실행해 보겠습니다.\n", "기준 에너지를 구하기 위해, 별도로 수행된 [DMRG](https://en.wikipedia.org/wiki/Density_matrix_renormalization_group) 계산 결과를 사용합니다.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "933037d8-847e-4986-80da-5ac8d677b2ff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using backend ibm_boston\n", "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -28.63965951544449\n", "\t\tSubspace dimension: 9801\n", "\tSubsample 1\n", "\t\tEnergy: -28.625588929202006\n", "\t\tSubspace dimension: 9409\n", "\tSubsample 2\n", "\t\tEnergy: -28.647371834135498\n", "\t\tSubspace dimension: 8281\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -28.67213260849567\n", "\t\tSubspace dimension: 29584\n", "\tSubsample 1\n", "\t\tEnergy: -28.670340686158816\n", "\t\tSubspace dimension: 27225\n", "\tSubsample 2\n", "\t\tEnergy: -28.669976379525988\n", "\t\tSubspace dimension: 31329\n", "Iteration 3\n", "\tSubsample 0\n", "\t\tEnergy: -28.68622875601382\n", "\t\tSubspace dimension: 36100\n", "\tSubsample 1\n", "\t\tEnergy: -28.698569623143126\n", "\t\tSubspace dimension: 34225\n", "\tSubsample 2\n", "\t\tEnergy: -28.694848533971882\n", "\t\tSubspace dimension: 33856\n", "Iteration 4\n", "\tSubsample 0\n", "\t\tEnergy: -28.69883392844593\n", "\t\tSubspace dimension: 42025\n", "\tSubsample 1\n", "\t\tEnergy: -28.701289495200996\n", "\t\tSubspace dimension: 38025\n", "\tSubsample 2\n", "\t\tEnergy: -28.699319594978245\n", "\t\tSubspace dimension: 45369\n", "Iteration 5\n", "\tSubsample 0\n", "\t\tEnergy: -28.701936886834154\n", "\t\tSubspace dimension: 51076\n", "\tSubsample 1\n", "\t\tEnergy: -28.702468711812013\n", "\t\tSubspace dimension: 53824\n", "\tSubsample 2\n", "\t\tEnergy: -28.702298147575938\n", "\t\tSubspace dimension: 52900\n", "Reference energy: -28.70660\n", "SQD energy: -28.70247\n", "Absolute error: 0.00413\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from qiskit_ibm_runtime import SamplerV2 as Sampler\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "\n", "# Model parameters\n", "norb = 20\n", "nelec = (norb // 2, norb // 2)\n", "n_bath = norb - 1\n", "hybridization = 1.0\n", "hopping = 1.0\n", "onsite = 10.0\n", "chemical_potential = -0.5 * onsite\n", "\n", "# Generate Hamiltonian and orbital rotation\n", "h1e, h2e = siam_hamiltonian(\n", " norb=norb,\n", " hopping=hopping,\n", " onsite=onsite,\n", " hybridization=hybridization,\n", " chemical_potential=chemical_potential,\n", ")\n", "orbital_rotation = momentum_basis(norb)\n", "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n", "impurity_index = n_bath // 2\n", "\n", "# Set reference energy to DMRG value computed separately\n", "reference_energy = -28.70659686\n", "\n", "# Algorithm parameters\n", "time_step = 0.2\n", "krylov_dim = 8\n", "\n", "# Construct circuits\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n", " circuit.append(instruction)\n", "circuit.measure_all()\n", "circuits = [circuit.copy()]\n", "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n", "for i in range(krylov_dim - 1):\n", " circuit.remove_final_measurements()\n", " for instruction in trotter_step(\n", " qubits,\n", " time_step,\n", " one_body_evolution,\n", " h2e_momentum,\n", " impurity_index,\n", " norb,\n", " ):\n", " circuit.append(instruction)\n", " circuit.measure_all()\n", " circuits.append(circuit.copy())\n", "\n", "# Initialize hardware backend\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " operational=True, simulator=False, min_num_qubits=127\n", ")\n", "print(f\"Using backend {backend.name}\")\n", "\n", "# Transpile to backend\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend\n", ")\n", "isa_circuits = pass_manager.run(circuits)\n", "\n", "# Sample from the circuits\n", "sampler = Sampler(backend)\n", "sampler.options.environment.job_tags = [\"TUT_SKQD\"]\n", "job = sampler.run(isa_circuits, shots=500)\n", "\n", "# Combine the shots from the individual Trotter circuits\n", "bit_array = BitArray.concatenate_shots(\n", " [result.data.meas for result in job.result()]\n", ")\n", "\n", "# Run configuration recovery and diagonalization\n", "result_history = []\n", "\n", "\n", "def callback(results: list[SCIResult]):\n", " result_history.append(results)\n", " iteration = len(result_history)\n", " print(f\"Iteration {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\"\\tSubsample {i}\")\n", " print(f\"\\t\\tEnergy: {result.energy}\")\n", " print(\n", " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", "\n", "rng = np.random.default_rng(24)\n", "result = diagonalize_fermionic_hamiltonian(\n", " h1e_momentum,\n", " h2e_momentum,\n", " bit_array,\n", " samples_per_batch=100,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=3,\n", " max_iterations=5,\n", " symmetrize_spin=True,\n", " callback=callback,\n", " seed=rng,\n", ")\n", "\n", "\n", "# Plot results\n", "min_es = [\n", " min(result, key=lambda res: res.energy).energy\n", " for result in result_history\n", "]\n", "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n", "x1 = range(len(result_history))\n", "y2 = np.sum(result.orbital_occupancies, axis=0)\n", "x2 = range(len(y2))\n", "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "axs[0].plot(x1, min_es, label=\"energy\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].axhline(\n", " y=reference_energy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"reference energy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy\", fontdict={\"fontsize\": 12})\n", "axs[0].legend()\n", "axs[1].bar(x2, y2, width=0.8)\n", "axs[1].set_xticks(x2)\n", "axs[1].set_xticklabels(x2)\n", "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n", "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n", "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n", "print(f\"Reference energy: {reference_energy:.5f}\")\n", "print(f\"SQD energy: {min_e:.5f}\")\n", "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "482ebea3-84b8-471b-bddc-b282e23162ad", "metadata": {}, "source": [ "\n", "\n", "## 다음 단계\n", "\n", "\n", " 이 글이 흥미로웠다면, 다음 자료도 참고해 보시기 바랍니다:\n", "\n", " * [화학 해밀토니안의 샘플 기반 양자 대각화](/docs/tutorials/sample-based-quantum-diagonalization) - 트로터 회로 대신 휴리스틱 변분 근사법을 사용한 관련 튜토리얼\n", " * [격자 해밀토니안의 크릴로프 양자 대각화(](/docs/tutorials/krylov-quantum-diagonalization) KQD) - KQD 방법에 대한 튜토리얼\n", " * [SQD 애드온 API 문서](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) - 함수 `diagonalize_fermionic_hamiltonian` 참조\n", " * [*샘플 기반 크릴로프 대각화법을 위한 양자 중심 알고리즘*](https://arxiv.org/abs/2501.09702) - 이 튜토리얼의 기초가 된 논문\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" }, "hours": 1.5, "qpuSeconds": 9 }, "nbformat": 4, "nbformat_minor": 4 }