{ "cells": [ { "cell_type": "markdown", "id": "afaf178f", "metadata": {}, "source": [ "---\n", "title: \"Krylov-Quanten-Diagonalisierung von Gitter-Hamilton-Operatoren\"\n", "description: \"Lernen Sie die Krylov-Quanten-Diagonalisierung für Gitter-Hamilton-Operatoren anhand des Konzepts der Hamilton-Funktion kennen.\"\n", "---\n", "\n", "{/* cspell:ignore Fomichev GEVP Izmaylov KQDOLH Kang Nakatsukasa Toeplitz eigenbasis eigendecomposition eigvals eigvecs evals noncommuting orthogonalizer sharex thresholded */}\n", "\n", "\n", "\n", "# Krylov-Quanten-Diagonalisierung von Gitter-Hamilton-Operatoren\n", "\n", "Geschätzte *Laufzeit: 70 Minuten auf einem Heron- oder Nighthawk-Prozessor (HINWEIS: Dies ist lediglich eine Schätzung. (Die Laufzeit kann variieren.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "aa4ba954", "metadata": {}, "source": [ "\n", "\n", "## Lernziele\n", "\n", "* Wie lässt sich die Krylov-Quanten-Diagonalisierung (KQD) als das Erlernen einer endlichen Hamilton-Funktion interpretieren, die als Spektralfilter fungiert?\n", "* Wie man die projizierten Hamilton- und Überlappungsmatrizen mithilfe erweiterter Swap-Test-Messungen erstellt.\n", "* Wie lässt sich das daraus resultierende verallgemeinerte Eigenwertproblem (GEVP) lösen und eine Schätzung der Grundzustandsenergie für einen Gitter-Hamilton-Operator ermitteln?\n", "\n", "\n", "\n", "## Voraussetzungen\n", "\n", "* [Qiskit primitives](/docs/guides/primitives)\n", "* [Schätzung der Grundzustandsenergie der Heisenberg-Kette mit VQE](/docs/tutorials/spin-chain-vqe)\n", "* [Kurs zu Algorithmen zur Quantendiagonalisierung](/learning/courses/quantum-diagonalization-algorithms/krylov)\n", "\n", "\n", "\n", "## Hintergrund\n", "\n", "Dieses Tutorial zeigt, wie man den Krylov-Quanten-Diagonalisierungsalgorithmus (KQD) im Rahmen von Qiskit-Patterns implementiert.\n", "Zunächst lernen Sie die Theorie hinter dem Algorithmus kennen und sehen anschließend eine Demonstration seiner Ausführung auf einer QPU.\n", "\n", "Die Abschätzung der Eigenschaften von Vielteilchen-Hamilton-Operatoren im Niedrigenergiebereich ist eine zentrale Aufgabe der Quantensimulation.\n", "So stehen beispielsweise Grundzustandsenergien und niedrigliegende Anregungen in direktem Zusammenhang mit der chemischen Stabilität, der magnetischen Ordnung, Quantenphasenübergängen und dem Materialverhalten.\n", "Auf einem klassischen Computer wächst die Dimension des Hilbert-Raums exponentiell mit der Anzahl der Orbitalen oder Spins, sodass eine direkte Diagonalisierung schnell unpraktikabel wird.\n", "\n", "Für dieses Problem gibt es mehrere Ansätze im Bereich des Quantencomputings.\n", "Kurzfristige Variationsmethoden wie der [Variational Quantum Eigensolver (VQE)](/docs/tutorials/spin-chain-vqe) verwenden relativ flache parametrisierte Schaltkreise, erfordern jedoch eine nichtlineare klassische Optimierungsschleife mit zahlreichen Auswertungen der Quantenschaltkreise.\n", "Auf der anderen Seite bietet [die Quantenphasenschätzung (QPE)](https://arxiv.org/abs/quant-ph/0604193) einen direkteren Weg zur Eigenwertschätzung mit strengen Garantien, doch die Standard-QPE erfordert lange kohärente Schaltkreise und eignet sich vor allem für fehlertolerante Quantencomputer.\n", "KQD liegt zwischen diesen beiden Ansätzen: Es nutzt die Hamilton-Evolution in Echtzeit, wie bei Algorithmen auf Basis der Phasenschätzung, ersetzt jedoch die vollständige Phasenschätzung durch ein kompaktes projiziertes Eigenwertproblem, das klassisch gelöst werden kann.\n", "\n", "Betrachten wir einen $n$ -Qubit-Hamilton-Operator $H$ und einen Referenzzustand $\\lvert \\psi_{0}\\rangle$.\n", "Die KQD-Methode erstellt einen Krylov-Unterraum aus in Echtzeit entwickelten Zuständen,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\lvert \\psi_\\ell\\rangle = e^{-i\\ell\\Delta t H}\\lvert \\psi_{0}\\rangle,\n", "\\qquad\n", "\\ell = 0,1,\\ldots,r-1,\n", "\\end{equation*}\n", "$$\n", "\n", "wobei $r$ die Krylov-Dimension und $\\Delta t$ der Zeitschritt ist.\n", "Jeder Zustand im Krylov-Unterraum lässt sich dann als Linearkombination dieser Basiszustände darstellen,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\lvert \\psi(\\mathbf{c})\\rangle\n", "=\n", "\\frac{\\sum_{\\ell=0}^{r-1} c_\\ell \\lvert \\psi_\\ell\\rangle}\n", "{\\left\\|\\sum_{\\ell=0}^{r-1} c_\\ell \\lvert \\psi_\\ell\\rangle\\right\\|},\n", "\\end{equation*}\n", "$$\n", "\n", "wobei der Nenner den Zustand normiert.\n", "\n", "Mit einfacher Algebra lässt sich zeigen, dass die entsprechende Energie als Rayleigh-Quotient ausgedrückt wird:\n", "\n", "$$\n", "\\begin{equation*}\n", "E(\\mathbf{c})\n", "=\\langle \\psi(\\mathbf{c})|H|\\psi(\\mathbf{c}) \\rangle=\n", "\\frac{\n", "\\sum_{k,\\ell} c_k^* c_\\ell \\langle \\psi_k\\vert H\\vert\\psi_\\ell\\rangle\n", "}{\n", "\\sum_{k,\\ell} c_k^* c_\\ell \\langle \\psi_k\\vert\\psi_\\ell\\rangle\n", "}\n", "=\n", "\\frac{\n", "\\mathbf{c}^{\\dagger}\\mathcal{H}\\mathbf{c}\n", "}{\n", "\\mathbf{c}^{\\dagger}\\mathcal{S}\\mathbf{c}\n", "}.\n", "\\end{equation*}\n", "$$\n", "\n", "Hier sind die Matrizen $\\mathcal{S}$ und $\\mathcal{H}$,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{S}_{k\\ell}=\\langle \\psi_k\\vert\\psi_\\ell\\rangle,\n", "\\qquad\n", "\\mathcal{H}_{k\\ell}=\\langle \\psi_k\\vert H\\vert\\psi_\\ell\\rangle\n", "\\end{equation*}\n", "$$\n", "\n", "Definieren Sie die projizierte Überlappungsmatrix und die Hamilton-Matrix.\n", "Ihre Einträge werden anhand von Messungen an Quantenschaltungen geschätzt.\n", "\n", "Wir möchten den Koeffizienten $\\mathbf{c}$ ermitteln, der das Minimum $E(\\mathbf{c})$ liefert:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\min_{\\mathbf{c}\\neq \\mathbf{0}} E(\\mathbf{c}).\n", "\\end{equation*}\n", "$$\n", "\n", "Nach dem Rayleigh-Ritz-Satz entspricht diese Minimierung der Lösung des verallgemeinerten Eigenwertproblems (GEVP),\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{H}\\mathbf{c}=E \\mathcal{S}\\mathbf{c}.\n", "\\end{equation*}\n", "$$\n", "\n", "Beachten Sie, dass die Dimension $r$ klein genug sein kann, damit ein klassischer Computer das GEVP lösen kann.\n", "\n", "Es handelt sich um dasselbe Variationsprinzip, das auch bei der klassischen Unterraumdiagonalisierung zum Einsatz kommt, doch hier werden die Basiszustände durch die Quantenzeitentwicklung erzeugt.\n", "Im Vergleich zu VQE erfordert KQD in der Regel tiefere Schaltkreise, da es auf Echtzeit-Entwicklung basiert.\n", "Im Gegenzug vermeidet KQD die nichtlineare Parameteroptimierung und die iterative Ausführung auf Quantenhardware und verbessert sich systematisch mit zunehmender Größe des projizierten Unterraums.\n", "Der Algorithmus wurde in großem Maßstab auf bestehender Quantenhardware demonstriert [\\[2\\]](#references), und seine Leistungsfähigkeit lässt sich mit nachweisbaren Garantien analysieren [\\[1\\]](#references).\n", "\n" ] }, { "cell_type": "markdown", "id": "3285eedc", "metadata": {}, "source": [ "\n", "\n", "## Anforderungen\n", "\n", "Bevor Sie mit dieser Anleitung beginnen, stellen Sie bitte sicher, dass folgende Komponenten installiert sind:\n", "\n", "* Qiskit SDK v2.3 oder höher mit [Visualisierungsunterstützung](/docs/api/qiskit/visualization)\n", "* Qiskit Runtime v0.22 oder höher (`pip install qiskit-ibm-runtime`)\n", "* SciPy (`pip install scipy`)\n", "* Matplotlib (`pip install matplotlib`)\n", "* Pandas (`pip install pandas`)\n", "\n", "Für die Ausführung auf der Hardware sind `qiskit-ibm-runtime` sowie der Zugriff auf ein „ IBM Quantum® “-Konto erforderlich.\n", "\n" ] }, { "cell_type": "markdown", "id": "047d8dcd", "metadata": {}, "source": [ "\n", "\n", "## Einrichtung\n", "\n", "Die Setup-Zelle importiert die erforderlichen Module und definiert Hilfsfunktionen für den Workflow:\n", "\n", "1. den Heisenberg-Hamilton-Operator aufstellen;\n", "2. das GEVP mit Schwellenwert lösen;\n", "3. den trainierten Krylov-Filter auswerten;\n", "4. Filterwerte in spektrale Gewichte umwandeln;\n", "5. Stellen Sie die Referenz- und die gefilterten Energieverteilungen grafisch dar.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "7bb044cf", "metadata": {}, "outputs": [], "source": [ "from __future__ import annotations\n", "\n", "import warnings\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import scipy.linalg as la\n", "import matplotlib.pyplot as plt\n", "\n", "from qiskit import QuantumCircuit, transpile\n", "from qiskit.circuit import Parameter\n", "from qiskit.circuit.library import PauliEvolutionGate\n", "from qiskit.primitives import StatevectorEstimator\n", "from qiskit.quantum_info import Operator, SparsePauliOp\n", "from qiskit.synthesis import LieTrotter, SuzukiTrotter\n", "from qiskit.transpiler import PassManager, Layout\n", "from qiskit.transpiler.passes import CommutativeOptimization\n", "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", "from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2, Batch\n", "from qiskit_ibm_runtime.fake_provider import FakeMarrakesh\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "\n", "def make_heisenberg_hamiltonian(\n", " num_qubits: int,\n", " coupling: float = 1.0,\n", ") -> SparsePauliOp:\n", " \"\"\"Make a Heisenberg Hamiltonian for a 1D chain of qubits with nearest-neighbor interactions.\"\"\"\n", " terms: list[tuple[str, complex]] = []\n", "\n", " def append_term(q0: int, q1: int, pauli: str):\n", " label = [\"I\"] * num_qubits\n", " label[num_qubits - 1 - q0] = pauli[0]\n", " label[num_qubits - 1 - q1] = pauli[1]\n", " terms.append((\"\".join(label), coupling))\n", "\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " for q in range(num_qubits - 1):\n", " append_term(q, q + 1, pauli)\n", " return SparsePauliOp.from_list(terms).simplify()\n", "\n", "\n", "def _basis_state_transition_amplitude_sparse(\n", " hamiltonian: SparsePauliOp,\n", " bra_state: int,\n", " ket_state: int,\n", ") -> complex:\n", " \"\"\"Evaluate for computational-basis states.\"\"\"\n", " num_qubits = hamiltonian.num_qubits\n", " amplitude = 0.0 + 0.0j\n", " for pauli, coeff in zip(hamiltonian.paulis, hamiltonian.coeffs):\n", " new_state = ket_state\n", " phase = 1.0 + 0.0j\n", " for q in range(num_qubits):\n", " x = bool(pauli.x[q])\n", " z = bool(pauli.z[q])\n", " if not x and not z:\n", " continue\n", " bit = (new_state >> q) & 1\n", " if x and z:\n", " # Y|0> = i|1>, Y|1> = -i|0>\n", " phase *= 1j if bit == 0 else -1j\n", " new_state ^= 1 << q\n", " elif x:\n", " new_state ^= 1 << q\n", " else:\n", " # Z|0> = |0>, Z|1> = -|1>\n", " if bit:\n", " phase *= -1\n", " if new_state == bra_state:\n", " amplitude += coeff * phase\n", " return amplitude\n", "\n", "\n", "def basis_state_expectation_sparse(\n", " hamiltonian: SparsePauliOp,\n", " bitstring: str,\n", ") -> complex:\n", " \"\"\"Evaluate .\"\"\"\n", " state = int(bitstring, 2)\n", " return _basis_state_transition_amplitude_sparse(\n", " hamiltonian,\n", " bra_state=state,\n", " ket_state=state,\n", " )\n", "\n", "\n", "def diagonalize_single_1_subspace(\n", " hamiltonian: SparsePauliOp,\n", ") -> np.ndarray:\n", " \"\"\"Diagonalize the Hamiltonian projected onto the single-excitation subspace.\"\"\"\n", " num_qubits = hamiltonian.num_qubits\n", "\n", " # Integer basis states |...010...>, with the excitation at qubit k.\n", " basis = [1 << k for k in range(num_qubits)]\n", "\n", " h_single = np.empty((num_qubits, num_qubits), dtype=complex)\n", "\n", " for row, bra_state in enumerate(basis):\n", " for col, ket_state in enumerate(basis):\n", " h_single[row, col] = _basis_state_transition_amplitude_sparse(\n", " hamiltonian,\n", " bra_state=bra_state,\n", " ket_state=ket_state,\n", " )\n", " # Remove floating-point-level asymmetry.\n", " h_single = 0.5 * (h_single + h_single.conj().T)\n", " evals, _ = np.linalg.eigh(h_single)\n", " return np.real(evals)\n", "\n", "\n", "def simple_transpilation(circuit: QuantumCircuit) -> QuantumCircuit:\n", " \"\"\"Transpilation to simplify the circuit\"\"\"\n", " pm = PassManager(\n", " [\n", " CommutativeOptimization(),\n", " ]\n", " )\n", " circuit = transpile(circuit, optimization_level=3)\n", " circuit = pm.run(circuit)\n", " return circuit\n", "\n", "\n", "def summarize_circuit(circuit: QuantumCircuit) -> dict[str, int | str]:\n", " \"\"\"Summarize the circuit with depth, size, and 2-qubit gate information.\"\"\"\n", " two_qubit_total = sum(\n", " inst.operation.num_qubits == 2 for inst in circuit.data\n", " )\n", " two_qubit_depth = circuit.depth(lambda x: x[0].num_qubits == 2)\n", " return {\n", " \"depth\": circuit.depth(),\n", " \"size\": circuit.size(),\n", " \"2q gates\": two_qubit_total,\n", " \"2q depth\": two_qubit_depth,\n", " }\n", "\n", "\n", "def solve_thresholded_gevp(\n", " h_matrix: np.ndarray,\n", " s_matrix: np.ndarray,\n", " threshold: float = 1e-10,\n", ") -> tuple[float, np.ndarray, int]:\n", " \"\"\"Solve H c = E S c using canonical orthogonalization of S.\"\"\"\n", " s_vals, s_vecs = la.eigh(s_matrix)\n", "\n", " valid = s_vals > threshold\n", " if not np.any(valid):\n", " raise ValueError(\n", " \"All overlap eigenvalues were removed by thresholding.\"\n", " )\n", "\n", " keep = valid\n", "\n", " orthogonalizer = s_vecs[:, keep] @ np.diag(1.0 / np.sqrt(s_vals[keep]))\n", " h_orth = orthogonalizer.conj().T @ h_matrix @ orthogonalizer\n", " h_orth = 0.5 * (h_orth + h_orth.conj().T)\n", "\n", " eigvals, eigvecs = la.eigh(h_orth)\n", "\n", " coeffs = orthogonalizer @ eigvecs[:, 0]\n", " normalization = np.sqrt(np.real(coeffs.conj().T @ s_matrix @ coeffs))\n", " coeffs /= normalization\n", "\n", " return float(np.real(eigvals[0])), coeffs, int(np.sum(keep))" ] }, { "cell_type": "markdown", "id": "7a03a282", "metadata": {}, "source": [ "Im ersten Teil dieses Tutorials veranschaulichen wir die KQD-Methode anhand eines Simulators mit lokalem Zustandsvektor. Später nutzen wir ein echtes Quanten-Backend, um ein Problem im industriellen Maßstab anzugehen.\n", "\n", "Außerdem definieren wir ein fiktives Backend, um die backendspezifische Transpilierung zu veranschaulichen und die daraus resultierende Schaltung zu untersuchen.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8026556a", "metadata": {}, "outputs": [], "source": [ "try:\n", " service = QiskitRuntimeService()\n", "except Exception:\n", " QiskitRuntimeService.save_account(\n", " token=\"\", instance=\"\", overwrite=True\n", " )\n", " service = QiskitRuntimeService()\n", "\n", "backend = FakeMarrakesh()" ] }, { "cell_type": "markdown", "id": "c23d1b80", "metadata": {}, "source": [ "\n", "\n", "## Beispiel für einen Simulator im kleinen Maßstab\n", "\n" ] }, { "cell_type": "markdown", "id": "af6d069a", "metadata": {}, "source": [ "\n", "\n", "### Schritt 1: Klassische Eingaben einem Quantenproblem zuordnen\n", "\n", "\n", "\n", "#### Hamilton-Funktor und Referenzzustand\n", "\n", "In diesem Beispiel wird eine Heisenberg-Kette mit 12 Qubits und offenen Grenzen verwendet ( $n=12$ ),\n", "\n", "$$\n", "\\begin{equation*}\n", "H=\\sum_{i=0}^{n-2}\\left(X_iX_{i+1}+Y_iY_{i+1}+Z_iZ_{i+1}\\right),\n", "\\end{equation*}\n", "$$\n", "\n", "mit einem Produktzustand mit einer einzigen Anregung\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\psi_{0}\\rangle=|000001000000\\rangle\n", "\\end{equation*}\n", "$$\n", "\n", "als Referenzzustand.\n", "Da der oben definierte Heisenberg-Hamilton-Operator die Gesamtzahl der Anregungen bewahrt, verbleibt der Referenzzustand im Unterraum mit einer einzigen Anregung, dessen Dimension nur linear mit der Anzahl der Qubits wächst. Wir können daher die exakte Grundzustandsenergie effizient berechnen, indem wir den auf diesen Unterraum beschränkten Hamiltonoperator diagonalisieren, und sie ausschließlich als diagnostischen Maßstab für die KQD-Schätzung verwenden. Der KQD-Workflow selbst schätzt die projizierten Matrixelemente mithilfe von „ Qiskit primitives “ und löst das daraus resultierende projizierte Problem auf klassische Weise.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "e156e4c7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Hamiltonian:\n", "SparsePauliOp(['IIIIIIIIIIXX', 'IIIIIIIIIXXI', 'IIIIIIIIXXII', 'IIIIIIIXXIII', 'IIIIIIXXIIII', 'IIIIIXXIIIII', 'IIIIXXIIIIII', 'IIIXXIIIIIII', 'IIXXIIIIIIII', 'IXXIIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIIYYI', 'IIIIIIIIYYII', 'IIIIIIIYYIII', 'IIIIIIYYIIII', 'IIIIIYYIIIII', 'IIIIYYIIIIII', 'IIIYYIIIIIII', 'IIYYIIIIIIII', 'IYYIIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIIZZI', 'IIIIIIIIZZII', 'IIIIIIIZZIII', 'IIIIIIZZIIII', 'IIIIIZZIIIII', 'IIIIZZIIIIII', 'IIIZZIIIIIII', 'IIZZIIIIIIII', 'IZZIIIIIIIII', 'ZZIIIIIIIIII'],\n", " coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "Reference state: |000001000000>\n", "Reference energy: (7+0j)\n" ] } ], "source": [ "# Problem definition for the simulator example\n", "num_qubits = 12\n", "hamiltonian = make_heisenberg_hamiltonian(num_qubits=num_qubits, coupling=1.0)\n", "ref_bitstring = \"000001000000\"\n", "ref_energy = basis_state_expectation_sparse(hamiltonian, ref_bitstring)\n", "\n", "print(\"Hamiltonian:\")\n", "print(hamiltonian)\n", "print(f\"Reference state: |{ref_bitstring}>\")\n", "print(\"Reference energy: \", ref_energy)" ] }, { "cell_type": "markdown", "id": "7f874d2a", "metadata": {}, "source": [ "\n", "\n", "#### Legen Sie die Parameter für den Algorithmus fest\n", "\n", "Ausgehend von den Obergrenzen für die Hamilton-Norm schlägt Ref. [\\[1\\]](#references) heuristisch die Zeitschritt $\\Delta t$ als „ $\\pi/\\|H\\|$ “ vor.\n", "Da die Spektralnorm $\\|H\\|$ schwer zu berechnen ist, verwenden wir stattdessen ihre Obergrenze:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\|H\\| \\le \\sum_{i=0}^{n-2}\\underbrace{\\|X_i X_{i+1}+Y_{i} Y_{i+1}+Z_{i} Z_{i+1}\\|}_{4 \\times 4 \\text{ matrix, easy to calculate the norm}} = 3(n-1).\n", "\\end{equation*}\n", "$$\n", "\n", "Wir haben die Krylov-Dimension auf $r=10$ und die Anzahl der Trotter-Schritte pro Zeitschritt auf $5$ festgelegt: ein ausreichend großer Krylov-Raum, um das tiefliegende Spektrum aufzulösen und gleichzeitig die tiefste Schaltung ( $t_{\\max}=(r-1)\\Delta t$ ) rechnerisch überschaubar zu halten, sowie genügend Trotter-Schritte, um den Diskretisierungsfehler bei dieser tiefsten Schaltung gering zu halten.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "af9a86df", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dt in Krylov basis: 0.09519977738150888\n" ] } ], "source": [ "dt = np.pi / (3 * (num_qubits - 1))\n", "print(\"dt in Krylov basis: \", dt)\n", "\n", "krylov_dim = 10\n", "num_trotter_steps = 5" ] }, { "cell_type": "markdown", "id": "5f02e4f3", "metadata": {}, "source": [ "\n", "\n", "#### Schaltung aufbauen\n", "\n", "Hier bauen wir die Schaltungen auf, um die Matrixelemente $\\mathcal{H}_{k\\ell}$ und $\\mathcal{S}_{k\\ell}$ zu schätzen.\n", "Da alle Potenzen von $H$ kommutieren, gilt\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{H}_{k\\ell}\n", "=\n", "\\langle\\psi_{0}|H e^{-i(\\ell-k)\\Delta tH}|\\psi_{0}\\rangle=\\mathcal{H}_{0,\\ell-k},\n", "\\qquad\n", "\\mathcal{S}_{k\\ell}\n", "=\n", "\\langle\\psi_{0}|e^{-i(\\ell-k)\\Delta tH}|\\psi_{0}\\rangle=\\mathcal{S}_{0,\\ell-k}.\n", "\\end{equation*}\n", "$$\n", "\n", "Solche Matrizen, bei denen die Elemente auf diese Weise von den Indexunterschieden abhängen, werden als *Toeplitz-Matrizen* bezeichnet und lassen sich aus den Elementen der ersten Zeile rekonstruieren, die durch $d=\\ell-k$ indiziert sind.\n", "\n", "Hier stellen wir die Schaltung vor, die als **„Extended-Swap-Test“** bezeichnet wird und die\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\Phi_{0d}\\rangle=\n", "\\frac{|0\\rangle|\\psi_0\\rangle+|1\\rangle|\\psi_d\\rangle}{\\sqrt{2}},\n", "\\end{equation*}\n", "$$\n", "\n", "wobei $|\\psi_d\\rangle=e^{-id\\Delta tH}|\\psi_{0}\\rangle$ gilt.\n", "\n", "\n", "\n", "#### Referenzstatus\n", "\n", "Wir berechnen den Referenzzustand $|\\psi_0\\rangle$.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "fb2879c5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "qc_ref = QuantumCircuit(num_qubits)\n", "for i, b in enumerate(reversed(ref_bitstring)):\n", " if b == \"1\":\n", " qc_ref.x(i)\n", "display(qc_ref.draw(\"mpl\", scale=0.5))" ] }, { "cell_type": "markdown", "id": "4a96258c", "metadata": {}, "source": [ "\n", "\n", "#### Zeitliche Entwicklung\n", "\n", "Berechnen Sie den durch den Hamilton-Operator erzeugten Zeitentwicklungsoperator, der durch [eine](/docs/api/qiskit/qiskit.synthesis.LieTrotter) einfache Lie-Trotter-Approximation approximiert wird.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "017b8548", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t = Parameter(\"t\")\n", "\n", "evol_gate = PauliEvolutionGate(\n", " hamiltonian,\n", " time=t,\n", " synthesis=LieTrotter(reps=num_trotter_steps),\n", " label=\"U(t)\",\n", ")\n", "\n", "# Synthesize U(t) first and then control the synthesized circuit.\n", "# This makes the controlled structure visible in the circuit drawer.\n", "evolution_circuit = QuantumCircuit(num_qubits, name=\"U(t)\")\n", "evolution_circuit.append(evol_gate, range(num_qubits))\n", "evolution_circuit = simple_transpilation(evolution_circuit)\n", "\n", "# Make a controlled version of the evolution circuit.\n", "controlled_evolution_gate = evolution_circuit.to_gate(label=\"U(t)\").control(\n", " 1, label=\"C-U(t)\"\n", ")\n", "display(\n", " evolution_circuit.assign_parameters({t: 1.5}).draw(\n", " \"mpl\", scale=0.5, fold=-1\n", " )\n", ")" ] }, { "cell_type": "markdown", "id": "efa6076d", "metadata": {}, "source": [ "\n", "\n", "#### Erweiterte Swap-Testschaltung [\\[3\\]](#references)\n", "\n", "Die Schaltung legt zunächst den Referenzzustand im Systemregister fest, während sich die Ancilla weiterhin im Zustand „ $|0\\rangle$ “ befindet:\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle |0\\rangle^{\\otimes n}\n", "\\longrightarrow\n", "|0\\rangle |\\psi_{0}\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Wendet man dann ein Hadamard-Gatter auf die Ancilla an, entsteht eine kohärente Überlagerung zweier Verzweigungen:\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle |\\psi_0\\rangle\n", "\\longrightarrow\n", "\\frac{|0\\rangle + |1\\rangle}{\\sqrt{2}} |\\psi_0\\rangle\n", "=\n", "\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Schließlich gilt das Gatter für die kontrollierte Zeitentwicklung:\n", "\n", "$$\n", "\\begin{equation*}\n", "U(t=d\\Delta t) = e^{-id\\Delta t H}\n", "\\end{equation*}\n", "$$\n", "\n", "nur wenn sich die Ancilla im Zweig „ $|1\\rangle$ “ befindet. Dies bedeutet Folgendes:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle|\\psi_0\\rangle}{\\sqrt{2}}\n", "\\longrightarrow\n", "\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle U(d\\Delta t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Im nächsten Code-Block implementieren wir:\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\Psi(t)\\rangle = \\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}},\n", "\\end{equation*}\n", "$$\n", "\n", "die im Ausführungsschritt als „ $t=d\\Delta t$ “ für „ $d=1,\\cdots r-1$ “ zugewiesen wird.\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "1260f66a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ancilla = 0\n", "system_qubits = list(range(1, num_qubits + 1))\n", "\n", "extended_swap_test = QuantumCircuit(num_qubits + 1)\n", "\n", "# Append state preparation part\n", "extended_swap_test = extended_swap_test.compose(qc_ref, system_qubits)\n", "\n", "# Prepare the coherent branch label, (|0> + |1>) / sqrt(2).\n", "extended_swap_test.h(ancilla)\n", "\n", "# Apply U(t) only to the |1> branch of the ancilla.\n", "extended_swap_test.append(\n", " controlled_evolution_gate, [ancilla] + system_qubits\n", ")\n", "\n", "# Decompose once more for visualization so that control bullets are visible.\n", "display(extended_swap_test.draw(\"mpl\", fold=-1))" ] }, { "cell_type": "markdown", "id": "b05e2e7d", "metadata": {}, "source": [ "\n", "\n", "#### Beobachtbare Größen\n", "\n", "Für jedes hermitische System mit den Observablen $O$ legen wir hier die zu berechnenden Observablen wie folgt fest:\n", "\n", "$$\n", "\\begin{equation*}\n", "z=\\langle \\psi_0|O|\\psi_d\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Das liegt daran, dass $O=I$ das Überlappungselement $\\mathcal{S}_{0d}=\\langle\\psi_0|\\psi_d\\rangle$ liefert, während $O=H$ das Hamilton-Element $\\mathcal{H}_{0d}=\\langle\\psi_0|H|\\psi_d\\rangle$ liefert.\n", "\n", "Unter Verwendung von $X=|0\\rangle\\langle 1|+|1\\rangle\\langle 0|$ erhalten wir\n", "\n", "$$\n", "\\begin{align*}\n", "\\langle \\Psi_{0d}|X\\otimes O| \\Psi_{0d} \\rangle\n", "&=\n", "\\left(\\frac{\\langle0|\\langle\\psi_0| + \\langle1| \\langle\\psi_d|}{\\sqrt{2}}\\right)(X\\otimes O)\\left(\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle |\\psi_d\\rangle}{\\sqrt{2}}\\right) \\\\\n", "&=\n", "\\frac{1}{2}\n", "\\left(\n", "\\langle \\psi_0|O|\\psi_d\\rangle\n", "+\n", "\\langle \\psi_d|O|\\psi_0\\rangle\n", "\\right) \\\\\n", "&=\n", "\\frac{z+z^*}{2}\n", "=\n", "\\operatorname{Re} z.\n", "\\end{align*}\n", "$$\n", "\n", "Ebenso lässt sich mit „ $Y=-i|0\\rangle\\langle 1|+i|1\\rangle\\langle 0|$ “\n", "\n", "$$\n", "\\begin{align*}\n", "\\langle \\Psi_{0d}|Y\\otimes O| \\Psi_{0d} \\rangle\n", "&=\n", "\\left(\\frac{\\langle0|\\langle\\psi_0| + \\langle1| \\langle\\psi_d|}{\\sqrt{2}}\\right)(Y\\otimes O)\\left(\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle |\\psi_d\\rangle}{\\sqrt{2}}\\right) \\\\\n", "&=\n", "\\frac{1}{2}\n", "\\left(\n", "-i\\langle \\psi_0|O|\\psi_d\\rangle\n", "+\n", "i\\langle \\psi_d|O|\\psi_0\\rangle\n", "\\right) \\\\\n", "&=\n", "\\frac{-iz+iz^*}{2}\n", "=\n", "\\operatorname{Im} z.\n", "\\end{align*}\n", "$$\n", "\n", "Daher gilt:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\langle \\Phi_{0d}|X\\otimes O|\\Phi_{0d}\\rangle\n", "=\n", "\\operatorname{Re}\\langle\\psi_0|O|\\psi_d\\rangle ,\n", "\\quad\n", "\\langle \\Phi_{0d}|Y\\otimes O|\\Phi_{0d}\\rangle\n", "=\n", "\\operatorname{Im}\\langle\\psi_0|O|\\psi_d\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Schließlich müssen wir für jeden Zustands $|\\Psi_{0d}\\rangle$ Folgendes messen:\n", "\n", "$$\n", "\\begin{align*}\n", "X \\otimes I \\text{ for } \\operatorname{Re}\\mathcal{S}_{0d},\\\\\n", "Y \\otimes I \\text{ for } \\operatorname{Im}\\mathcal{S}_{0d},\\\\\n", "X \\otimes H \\text{ for } \\operatorname{Re}\\mathcal{H}_{0d},\\\\\n", "Y \\otimes H \\text{ for } \\operatorname{Im}\\mathcal{H}_{0d}.\\\\\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "4c563534", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Observable: Re S_0d\n", "SparsePauliOp(['IIIIIIIIIIIIX'],\n", " coeffs=[1.+0.j])\n", "\n", "Observable: Im S_0d\n", "SparsePauliOp(['IIIIIIIIIIIIY'],\n", " coeffs=[1.+0.j])\n", "\n", "Observable: Re H_0d\n", "SparsePauliOp(['IIIIIIIIIIXXX', 'IIIIIIIIIXXIX', 'IIIIIIIIXXIIX', 'IIIIIIIXXIIIX', 'IIIIIIXXIIIIX', 'IIIIIXXIIIIIX', 'IIIIXXIIIIIIX', 'IIIXXIIIIIIIX', 'IIXXIIIIIIIIX', 'IXXIIIIIIIIIX', 'XXIIIIIIIIIIX', 'IIIIIIIIIIYYX', 'IIIIIIIIIYYIX', 'IIIIIIIIYYIIX', 'IIIIIIIYYIIIX', 'IIIIIIYYIIIIX', 'IIIIIYYIIIIIX', 'IIIIYYIIIIIIX', 'IIIYYIIIIIIIX', 'IIYYIIIIIIIIX', 'IYYIIIIIIIIIX', 'YYIIIIIIIIIIX', 'IIIIIIIIIIZZX', 'IIIIIIIIIZZIX', 'IIIIIIIIZZIIX', 'IIIIIIIZZIIIX', 'IIIIIIZZIIIIX', 'IIIIIZZIIIIIX', 'IIIIZZIIIIIIX', 'IIIZZIIIIIIIX', 'IIZZIIIIIIIIX', 'IZZIIIIIIIIIX', 'ZZIIIIIIIIIIX'],\n", " coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n", "Observable: Im H_0d\n", "SparsePauliOp(['IIIIIIIIIIXXY', 'IIIIIIIIIXXIY', 'IIIIIIIIXXIIY', 'IIIIIIIXXIIIY', 'IIIIIIXXIIIIY', 'IIIIIXXIIIIIY', 'IIIIXXIIIIIIY', 'IIIXXIIIIIIIY', 'IIXXIIIIIIIIY', 'IXXIIIIIIIIIY', 'XXIIIIIIIIIIY', 'IIIIIIIIIIYYY', 'IIIIIIIIIYYIY', 'IIIIIIIIYYIIY', 'IIIIIIIYYIIIY', 'IIIIIIYYIIIIY', 'IIIIIYYIIIIIY', 'IIIIYYIIIIIIY', 'IIIYYIIIIIIIY', 'IIYYIIIIIIIIY', 'IYYIIIIIIIIIY', 'YYIIIIIIIIIIY', 'IIIIIIIIIIZZY', 'IIIIIIIIIZZIY', 'IIIIIIIIZZIIY', 'IIIIIIIZZIIIY', 'IIIIIIZZIIIIY', 'IIIIIZZIIIIIY', 'IIIIZZIIIIIIY', 'IIIZZIIIIIIIY', 'IIZZIIIIIIIIY', 'IZZIIIIIIIIIY', 'ZZIIIIIIIIIIY'],\n", " coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n" ] } ], "source": [ "n_qubits = hamiltonian.num_qubits\n", "\n", "observable_labels = [\n", " \"Re S_0d\",\n", " \"Im S_0d\",\n", " \"Re H_0d\",\n", " \"Im H_0d\",\n", "]\n", "\n", "# X ⊗ I and Y ⊗ I.\n", "# Qiskit's Pauli-label convention places qubit 0 on the rightmost character,\n", "# so the ancilla Pauli is appended to the right.\n", "obs_x_identity = SparsePauliOp(\"I\" * n_qubits + \"X\")\n", "obs_y_identity = SparsePauliOp(\"I\" * n_qubits + \"Y\")\n", "\n", "# X ⊗ H and Y ⊗ H.\n", "obs_x_hamiltonian = SparsePauliOp.from_list(\n", " [\n", " (label + \"X\", coeff)\n", " for label, coeff in zip(\n", " hamiltonian.paulis.to_labels(),\n", " hamiltonian.coeffs,\n", " )\n", " ]\n", ")\n", "\n", "obs_y_hamiltonian = SparsePauliOp.from_list(\n", " [\n", " (label + \"Y\", coeff)\n", " for label, coeff in zip(\n", " hamiltonian.paulis.to_labels(),\n", " hamiltonian.coeffs,\n", " )\n", " ]\n", ")\n", "\n", "observables = [\n", " obs_x_identity,\n", " obs_y_identity,\n", " obs_x_hamiltonian,\n", " obs_y_hamiltonian,\n", "]\n", "\n", "for obs, label in zip(observables, observable_labels):\n", " print(f\"Observable: {label}\")\n", " print(obs)\n", " print()" ] }, { "cell_type": "markdown", "id": "780486da", "metadata": {}, "source": [ "Bei der Abschätzung der Matrixelemente von „ $\\mathcal{H}$ “ ist die Anzahl der Pauli-Terme wesentlich größer als bei den Matrixelementen von „ $\\mathcal{S}$ “.\n", "\n", "Mithilfe einer Verschiebungstechnik [\\[4\\]](#references) können wir nun die Anzahl der gemessenen Hamilton-Terme reduzieren.\n", "Wir zerlegen den Hamilton-Operator wie folgt:\n", "\n", "$$\n", "\\begin{equation*}\n", "H = (H-T) + T,\n", "\\end{equation*}\n", "$$\n", "\n", "wobei $T$ so gewählt wird, dass der Referenzzustand sein Eigenzustand ist,\n", "\n", "$$\n", "\\begin{equation*}\n", "T|\\psi_0\\rangle = \\tau |\\psi_0\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Dann,\n", "\n", "$$\n", "\\begin{align*}\n", "\\mathcal{H}_{0d}\n", "&=\n", "\\langle\\psi_0|H|\\psi_d\\rangle \\\\\n", "&=\n", "\\langle\\psi_0|(H-T)|\\psi_d\\rangle\n", "+\n", "\\langle\\psi_0|T|\\psi_d\\rangle \\\\\n", "&=\n", "\\langle\\psi_0|(H-T)|\\psi_d\\rangle\n", "+\n", "\\tau \\langle\\psi_0|\\psi_d\\rangle \\\\\n", "&=\n", "\\widetilde{\\mathcal{H}}_{0d}\n", "+\n", "\\tau \\mathcal{S}_{0d}.\n", "\\end{align*}\n", "$$\n", "\n", "In dieser Datei gilt Folgendes:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\widetilde{\\mathcal{H}}_{0d}\n", "=\n", "\\langle\\psi_0|(H-T)|\\psi_d\\rangle\n", "\\end{equation*}\n", "$$\n", "\n", "ist das verschobene Matrixelement des Hamilton-Operators.\n", "Daher müssen wir lediglich $X\\otimes (H-T)$ und $Y\\otimes (H-T)$ messen.\n", "Der Beitrag von $T$ wird klassisch anhand des bereits gemessenen Überlappungsmatrixelements $\\mathcal{S}_{0d}$ rekonstruiert.\n", "\n", "In diesem Beispiel bietet sich der diagonale Teil des Heisenberg-Hamilton-Operators an,\n", "\n", "$$\n", "\\begin{equation*}\n", "T\n", "=\n", "\\sum_{i=0}^{n-2} Z_iZ_{i+1}.\n", "\\end{equation*}\n", "$$\n", "\n", "Da es sich bei dem Referenzzustand um einen Berechnungsbasiszustand handelt, ist er ein Eigenzustand jedes „ $Z_iZ_{i+1}$ “-Terms.\n", "\n", "Eine vorteilhaftere Wahl ist es jedoch, nicht nur die diagonalen $ZZ$ Terme einzubeziehen, sondern auch die $XX+YY$ Terme, die den Referenzzustand auslöschen.\n", "\n", "Für jedes Nachbarpaar gilt für den Operator $XX+YY$ :\n", "\n", "$$\n", "\\begin{equation*}\n", "(XX+YY)|00\\rangle = 0,\n", "\\qquad\n", "(XX+YY)|11\\rangle = 0,\n", "\\end{equation*}\n", "$$\n", "\n", "und\n", "\n", "$$\n", "\\begin{equation*}\n", "(XX+YY)|01\\rangle = 2|10\\rangle,\n", "\\qquad\n", "(XX+YY)|10\\rangle = 2|01\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Daher trägt der Term „ $XX+YY$ “ nur dann bei, wenn die beiden benachbarten Qubits unterschiedliche Belegungen in der Referenz-Bitfolge aufweisen.\n", "Sind die beiden Qubits entweder beide „ $0$ “ oder beide „ $1$ “, so löscht dieser Term den Referenzzustand aus und kann ebenfalls herausgenommen werden.\n", "\n", "Sei $|\\psi_{\\rm ref}\\rangle = |z_0\\cdots z_{n-1}\\rangle$, wobei $z_i \\in \\{0,1\\}$ gilt.\n", "Somit können wir wählen\n", "\n", "$$\n", "\\begin{equation*}\n", "T=\\sum_{i=0}^{n-2} Z_iZ_{i+1}\n", "+\n", "\\sum_{z_i={z_{i+1}}}\n", "\\left(\n", "X_iX_{i+1}\n", "+\n", "Y_iY_{i+1}\n", "\\right).\n", "\\end{equation*}\n", "$$\n", "\n", "Dieser Operator erfüllt weiterhin\n", "\n", "$$\n", "\\begin{equation*}\n", "T|\\psi_0\\rangle = \\tau |\\psi_0\\rangle ,\n", "\\end{equation*}\n", "$$\n", "\n", "denn die Terme von $ZZ$ wirken diagonal auf $|\\psi_0\\rangle$ ein, während die verschobenen Terme von $XX+YY$ den Wert Null ergeben.\n", "Der entsprechende Eigenwert wird daher ausschließlich durch die Terme „ $ZZ$ “ bestimmt,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\tau = \\sum_{i=0}^{n-2}\n", "(-1)^{z_i}\n", "(-1)^{z_{i+1}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Mit dieser Wahl lautet der verschobene Hamilton-Operator nun:\n", "\n", "$$\n", "\\begin{equation*}\n", "H-T = \\sum_{z_i\\ne z_{i+1}}\n", "\\left(\n", "X_iX_{i+1}\n", "+\n", "Y_iY_{i+1}\n", "\\right).\n", "\\end{equation*}\n", "$$\n", "\n", "Daher müssen nur die Kanten gemessen werden, an denen im Referenzzustand unterschiedliche Berufe ausgeübt werden.\n", "Alle Terme von $ZZ$ und alle inaktiven Terme von $XX+YY$ werden durch den Überlappungsbeitrag $\\tau \\mathcal{S}_{0d}$ rekonstruiert oder liefern per Konstruktion einen Beitrag von Null.\n", "\n", "Dadurch ergibt sich eine kleinere Observable als bei einer Verschiebung nur des diagonalen Teils.\n", "Insbesondere bleiben bei einem Referenzzustand auf Berechnungsbasis mit einer lokalisierten Anregung nur die an die Anregung angrenzenden Kanten in $H-T$ erhalten.\n", "Daher kann die Anzahl der Pauli-Terme in $X\\otimes(H-T)$ und $Y\\otimes(H-T)$ erheblich reduziert werden, während das rekonstruierte Matrixelement\n", "\n", "$$\n", "\\begin{equation*}\n", "\\widetilde{\\mathcal{H}}_{0d}\n", "+\n", "\\tau \\mathcal{S}_{0d}\n", "\\end{equation*}\n", "$$\n", "\n", "bleibt genau gleich.\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "52a9b116", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Observable: Re shifted H_0d\n", "SparsePauliOp(['IIIIIXXIIIIIX', 'IIIIIYYIIIIIX', 'IIIIXXIIIIIIX', 'IIIIYYIIIIIIX'],\n", " coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n", "Observable: Im shifted H_0d\n", "SparsePauliOp(['IIIIIXXIIIIIY', 'IIIIIYYIIIIIY', 'IIIIXXIIIIIIY', 'IIIIYYIIIIIIY'],\n", " coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n", "tau = 7.0\n" ] } ], "source": [ "def make_reduced_heisenberg_observables(\n", " ref_bitstring: str,\n", " coupling: float = 1.0,\n", ") -> tuple[SparsePauliOp, SparsePauliOp, float]:\n", " \"\"\"Build X⊗(H-T), Y⊗(H-T), and tau.\"\"\"\n", " n_qubits = len(ref_bitstring)\n", "\n", " shifted_terms: list[tuple[str, complex]] = []\n", " tau = 0.0\n", "\n", " def bit(q: int) -> str:\n", " return ref_bitstring[n_qubits - 1 - q]\n", "\n", " def append_term(q0: int, q1: int, pauli: str):\n", " label = [\"I\"] * n_qubits\n", " label[n_qubits - 1 - q0] = pauli[0]\n", " label[n_qubits - 1 - q1] = pauli[1]\n", " shifted_terms.append((\"\".join(label), coupling))\n", "\n", " for q in range(n_qubits - 1):\n", " same_occupation = bit(q) == bit(q + 1)\n", "\n", " # ZZ contribution to tau\n", " tau += coupling * (1.0 if same_occupation else -1.0)\n", "\n", " # XX + YY survives only for opposite occupations.\n", " if not same_occupation:\n", " append_term(q, q + 1, \"XX\")\n", " append_term(q, q + 1, \"YY\")\n", "\n", " if shifted_terms:\n", " obs_x_shifted_hamiltonian = SparsePauliOp.from_list(\n", " [(label + \"X\", coeff) for label, coeff in shifted_terms]\n", " )\n", " obs_y_shifted_hamiltonian = SparsePauliOp.from_list(\n", " [(label + \"Y\", coeff) for label, coeff in shifted_terms]\n", " )\n", " else:\n", " obs_x_shifted_hamiltonian = SparsePauliOp(\n", " \"I\" * n_qubits + \"X\", coeffs=[0.0]\n", " )\n", " obs_y_shifted_hamiltonian = SparsePauliOp(\n", " \"I\" * n_qubits + \"Y\", coeffs=[0.0]\n", " )\n", "\n", " return obs_x_shifted_hamiltonian, obs_y_shifted_hamiltonian, tau\n", "\n", "\n", "obs_x_shifted_hamiltonian, obs_y_shifted_hamiltonian, shift_tau = (\n", " make_reduced_heisenberg_observables(ref_bitstring)\n", ")\n", "\n", "print(\"Observable: Re shifted H_0d\")\n", "print(obs_x_shifted_hamiltonian)\n", "print()\n", "\n", "print(\"Observable: Im shifted H_0d\")\n", "print(obs_y_shifted_hamiltonian)\n", "print()\n", "\n", "print(\"tau =\", shift_tau)\n", "\n", "observables = [\n", " obs_x_identity,\n", " obs_y_identity,\n", " obs_x_shifted_hamiltonian,\n", " obs_y_shifted_hamiltonian,\n", "]\n", "observable_labels = [\n", " \"Re S_0d\",\n", " \"Im S_0d\",\n", " \"Re shifted H_0d\",\n", " \"Im shifted H_0d\",\n", "]" ] }, { "cell_type": "markdown", "id": "7c7daa55", "metadata": {}, "source": [ "\n", "\n", "### Schritt 2: Optimierung des Problems für die Ausführung auf Quantenhardware\n", "\n", "Nun wandeln wir die abstrakte Schaltung für den erweiterten Swap-Test in eine hardwareorientierte Vorlage um.\n", "Zuvor optimieren wir die Schaltung zunächst auf abstrakter Ebene weiter.\n", "\n", "\n", "\n", "#### Vergleich der Reihenfolge der Hamilton-Terme\n", "\n", "Zunächst vergleichen wir verschiedene Reihenfolgen der Pauli-Terme im Heisenberg-Hamilton-Operator für die Hamilton-Simulation.\n", "Der Hamilton-Operator selbst bleibt unverändert, doch die Reihenfolge beeinflusst, wie der Schaltkreis in Produktform generiert wird und inwieweit die Struktur im Schaltkreis parallelisiert werden kann.\n", "Beispielsweise listet die naive Reihenfolge zunächst alle „Nearest-Neighbor“- $XX$ -Terme auf, dann alle „ $YY$ “-Terme und schließlich alle „ $ZZ$ “-Terme.\n", "Dadurch werden benachbarte Kanten wie $(0,1)$ und $(1,2)$ nebeneinander platziert, sodass sie nicht parallel ausgeführt werden können.\n", "Bei der „even-then-odd“-Reihenfolge werden zunächst die disjunkten geraden Kanten und anschließend die ungeraden Kanten durchlaufen, wodurch parallele Zwei-Qubit-Schichten freigelegt werden.\n", "Die nach geraden und ungeraden Kanten gruppierte Reihenfolge geht noch einen Schritt weiter: Für jede Kante werden die lokalen Terme „ $XX$ “, „ $YY$ “ und „ $ZZ$ “ beibehalten, während weiterhin gerade Kanten vor ungeraden Kanten besucht werden.\n", "Wir gehen davon aus, dass die „even-then-odd“- und „even-odd“-Anordnungen mit Kantengruppierung die Schaltungstiefe verringern, indem sie parallele Zwei-Qubit-Schichten freilegen, und dass die Anordnung mit Kantengruppierung zudem den Trotter-Fehler reduziert, da die lokale Zwei-Qubit-Wechselwirkung an derselben Kante als kompakter Block behandelt wird.\n", "\n", "![Pauli-Entwicklungen mit unterschiedlicher Reihenfolge der Terme](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-quantum-diagonalization/ham_ordering.avif)\n", "\n", "Die Hamilton-Reihenfolge wirkt sich auch auf den Trotter-Fehler aus.\n", "Wenn wir nichtkommutierende Terme nebeneinander setzen, treten die Basisübergänge häufiger auf, was zu einem größeren Trotter-Fehler führt.\n", "Durch die Gruppierung von Termen, die dieselbe Pauli-Basistransformation erfordern, lassen sich überflüssige Basiswechsel vermeiden.\n", "\n", "Hier wird für den Vergleich die größte Evolutionszeit herangezogen, die in den Krylov-Schätzungen der ersten Zeile auftritt (\n", "$t_{\\max} = (r-1)\\Delta t$ ),\n", "wobei dieselbe Transpilationsbedingung verwendet wird.\n", "\n", "Um den Trotter-Fehler zu messen, verwenden wir **die Prozessuntreue** zwischen der Trotterisierten Schaltung und der exakten Hamilton-Entwicklung,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\text{Infidelity}(U,V)\n", "=\n", "1-F(U,V)\n", "=\n", "1 - \\frac{|\\operatorname{Tr}(U^\\dagger V)|^2}{d^2},\n", "\\end{equation*}\n", "$$\n", "\n", "wobei $d=2^n$ die Dimension des Hilbert-Raums ist.\n", "\n", "Diese Diagnose verwendet dichte Matrizen und eignet sich daher für dieses kleine 12-Qubit-Beispiel, ist jedoch nicht als skalierbare Subroutine gedacht.\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "d21eae64", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Comparison time: 0.8567979964335799\n", "\n", "naive:\n", "['IIIIIIIIIIXX', 'IIIIIIIIIXXI', 'IIIIIIIIXXII', 'IIIIIIIXXIII', 'IIIIIIXXIIII', 'IIIIIXXIIIII', 'IIIIXXIIIIII', 'IIIXXIIIIIII', 'IIXXIIIIIIII', 'IXXIIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIIYYI', 'IIIIIIIIYYII', 'IIIIIIIYYIII', 'IIIIIIYYIIII', 'IIIIIYYIIIII', 'IIIIYYIIIIII', 'IIIYYIIIIIII', 'IIYYIIIIIIII', 'IYYIIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIIZZI', 'IIIIIIIIZZII', 'IIIIIIIZZIII', 'IIIIIIZZIIII', 'IIIIIZZIIIII', 'IIIIZZIIIIII', 'IIIZZIIIIIII', 'IIZZIIIIIIII', 'IZZIIIIIIIII', 'ZZIIIIIIIIII']\n", "\n", "even-then-odd:\n", "['IIIIIIIIIIXX', 'IIIIIIIIXXII', 'IIIIIIXXIIII', 'IIIIXXIIIIII', 'IIXXIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIXXIII', 'IIIIIXXIIIII', 'IIIXXIIIIIII', 'IXXIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIYYII', 'IIIIIIYYIIII', 'IIIIYYIIIIII', 'IIYYIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIYYI', 'IIIIIIIYYIII', 'IIIIIYYIIIII', 'IIIYYIIIIIII', 'IYYIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIZZII', 'IIIIIIZZIIII', 'IIIIZZIIIIII', 'IIZZIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIZZI', 'IIIIIIIZZIII', 'IIIIIZZIIIII', 'IIIZZIIIIIII', 'IZZIIIIIIIII']\n", "\n", "even-odd edge-grouped:\n", "['IIIIIIIIIIXX', 'IIIIIIIIIIYY', 'IIIIIIIIIIZZ', 'IIIIIIIIXXII', 'IIIIIIIIYYII', 'IIIIIIIIZZII', 'IIIIIIXXIIII', 'IIIIIIYYIIII', 'IIIIIIZZIIII', 'IIIIXXIIIIII', 'IIIIYYIIIIII', 'IIIIZZIIIIII', 'IIXXIIIIIIII', 'IIYYIIIIIIII', 'IIZZIIIIIIII', 'XXIIIIIIIIII', 'YYIIIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIIIYYI', 'IIIIIIIIIZZI', 'IIIIIIIXXIII', 'IIIIIIIYYIII', 'IIIIIIIZZIII', 'IIIIIXXIIIII', 'IIIIIYYIIIII', 'IIIIIZZIIIII', 'IIIXXIIIIIII', 'IIIYYIIIIIII', 'IIIZZIIIIIII', 'IXXIIIIIIIII', 'IYYIIIIIIIII', 'IZZIIIIIIIII']\n", "\n", "Precomputing the exact evolution operator... Done\n" ] } ], "source": [ "# The comparison uses the largest time that appears in the first-row Krylov estimates.\n", "# Circuit depth does not depend on this numeric value, but the Trotter error does.\n", "comparison_time = (krylov_dim - 1) * dt\n", "\n", "\n", "def make_heisenberg_hamiltonian_ordered(\n", " num_qubits: int,\n", " ordering: str,\n", " coupling: float = 1.0,\n", ") -> SparsePauliOp:\n", " \"\"\"Return the same Heisenberg Hamiltonian with a specified term ordering.\"\"\"\n", " terms: list[tuple[str, complex]] = []\n", " even_edges = [(q, q + 1) for q in range(0, num_qubits - 1, 2)]\n", " odd_edges = [(q, q + 1) for q in range(1, num_qubits - 1, 2)]\n", "\n", " def append_term(q0: int, q1: int, pauli: str):\n", " label = [\"I\"] * num_qubits\n", " label[num_qubits - 1 - q0] = pauli[0]\n", " label[num_qubits - 1 - q1] = pauli[1]\n", " terms.append((\"\".join(label), coupling))\n", "\n", " if ordering == \"naive\":\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " for q in range(num_qubits - 1):\n", " append_term(q, q + 1, pauli)\n", " elif ordering == \"even-then-odd\":\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " for q0, q1 in even_edges + odd_edges:\n", " append_term(q0, q1, pauli)\n", " elif ordering == \"even-odd edge-grouped\":\n", " for q0, q1 in even_edges + odd_edges:\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " append_term(q0, q1, pauli)\n", " else:\n", " raise ValueError(f\"Unknown ordering: {ordering}\")\n", "\n", " return SparsePauliOp.from_list(terms).simplify()\n", "\n", "\n", "def build_numeric_evolution_circuit(\n", " hamiltonian: SparsePauliOp,\n", " synthesis,\n", " time_value: float,\n", " **synthesis_kwargs,\n", ") -> QuantumCircuit:\n", " \"\"\"Build a numeric circuit for exp(-i H t) with a chosen synthesis rule.\"\"\"\n", " evolution_gate = PauliEvolutionGate(\n", " hamiltonian,\n", " time=time_value,\n", " synthesis=synthesis(**synthesis_kwargs),\n", " )\n", " circuit = QuantumCircuit(hamiltonian.num_qubits)\n", " circuit.append(evolution_gate, range(hamiltonian.num_qubits))\n", " return circuit\n", "\n", "\n", "def process_infidelity(\n", " circuit: QuantumCircuit,\n", " exact_matrix: np.ndarray,\n", ") -> float:\n", " \"\"\"Return 1 - |Tr(U_circuit† U_exact) / d|².\"\"\"\n", " circuit_matrix = np.asarray(Operator(circuit).data)\n", " dim = circuit_matrix.shape[0]\n", " normalized_trace = np.vdot(circuit_matrix, exact_matrix) / dim\n", " fidelity = np.abs(normalized_trace) ** 2\n", " return float(np.clip(1.0 - fidelity, 0.0, 1.0))\n", "\n", "\n", "hamiltonians_by_ordering = {\n", " ordering: make_heisenberg_hamiltonian_ordered(\n", " num_qubits, ordering, coupling=1.0\n", " )\n", " for ordering in [\"naive\", \"even-then-odd\", \"even-odd edge-grouped\"]\n", "}\n", "\n", "print(f\"Comparison time: {comparison_time}\\n\")\n", "\n", "for order_name, ham_ordered in hamiltonians_by_ordering.items():\n", " print(f\"{order_name}:\")\n", " print([op for op, _ in ham_ordered.to_list()])\n", " print()\n", "\n", "print(\"Precomputing the exact evolution operator... \", end=\"\")\n", "exact_matrix = la.expm(-1j * comparison_time * hamiltonian.to_matrix())\n", "print(\"Done\")" ] }, { "cell_type": "markdown", "id": "d3fcb5e1", "metadata": {}, "source": [ "Zunächst belassen wir die Syntheseregel auf einen Trotter-Schritt erster Ordnung und variieren lediglich die Reihenfolge der Pauli-Terme.\n", "Das Ziel dieses Vergleichs besteht vor allem darin, zu ermitteln, um wie viel die Schaltungstiefe und die Kosten pro zwei Qubits reduziert werden können, wenn nicht überlappende Kanten zwischen nächsten Nachbarn dem Transpiler zugänglich gemacht werden.\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "21eaab38", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Circuit for naive ordering:\n", "['IIIIIIIIIIXX', 'IIIIIIIIIXXI', 'IIIIIIIIXXII', 'IIIIIIIXXIII', 'IIIIIIXXIIII', 'IIIIIXXIIIII', 'IIIIXXIIIIII', 'IIIXXIIIIIII', 'IIXXIIIIIIII', 'IXXIIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIIYYI', 'IIIIIIIIYYII', 'IIIIIIIYYIII', 'IIIIIIYYIIII', 'IIIIIYYIIIII', 'IIIIYYIIIIII', 'IIIYYIIIIIII', 'IIYYIIIIIIII', 'IYYIIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIIZZI', 'IIIIIIIIZZII', 'IIIIIIIZZIII', 'IIIIIIZZIIII', 'IIIIIZZIIIII', 'IIIIZZIIIIII', 'IIIZZIIIIIII', 'IIZZIIIIIIII', 'IZZIIIIIIIII', 'ZZIIIIIIIIII']\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Circuit for even-then-odd ordering:\n", "['IIIIIIIIIIXX', 'IIIIIIIIXXII', 'IIIIIIXXIIII', 'IIIIXXIIIIII', 'IIXXIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIXXIII', 'IIIIIXXIIIII', 'IIIXXIIIIIII', 'IXXIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIYYII', 'IIIIIIYYIIII', 'IIIIYYIIIIII', 'IIYYIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIYYI', 'IIIIIIIYYIII', 'IIIIIYYIIIII', 'IIIYYIIIIIII', 'IYYIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIZZII', 'IIIIIIZZIIII', 'IIIIZZIIIIII', 'IIZZIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIZZI', 'IIIIIIIZZIII', 'IIIIIZZIIIII', 'IIIZZIIIIIII', 'IZZIIIIIIIII']\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Circuit for even-odd edge-grouped ordering:\n", "['IIIIIIIIIIXX', 'IIIIIIIIIIYY', 'IIIIIIIIIIZZ', 'IIIIIIIIXXII', 'IIIIIIIIYYII', 'IIIIIIIIZZII', 'IIIIIIXXIIII', 'IIIIIIYYIIII', 'IIIIIIZZIIII', 'IIIIXXIIIIII', 'IIIIYYIIIIII', 'IIIIZZIIIIII', 'IIXXIIIIIIII', 'IIYYIIIIIIII', 'IIZZIIIIIIII', 'XXIIIIIIIIII', 'YYIIIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIIIYYI', 'IIIIIIIIIZZI', 'IIIIIIIXXIII', 'IIIIIIIYYIII', 'IIIIIIIZZIII', 'IIIIIXXIIIII', 'IIIIIYYIIIII', 'IIIIIZZIIIII', 'IIIXXIIIIIII', 'IIIYYIIIIIII', 'IIIZZIIIIIII', 'IXXIIIIIIIII', 'IYYIIIIIIIII', 'IZZIIIIIIIII']\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
orderingsynthesisinfidelitydepthsize2q gates2q depth
0naiveLieTrotter(reps=1)0.99991715333315
1naiveLieTrotter(reps=2)0.80599821666621
2naiveLieTrotter(reps=4)0.2713893313213233
3naiveLieTrotter(reps=8)0.0745905726426457
4even-then-oddLieTrotter(reps=1)0.999917633336
5even-then-oddLieTrotter(reps=2)0.80599812666612
6even-then-oddLieTrotter(reps=4)0.2713892413213224
7even-then-oddLieTrotter(reps=8)0.0745904826426448
8even-odd edge-groupedLieTrotter(reps=1)0.998432633336
9even-odd edge-groupedLieTrotter(reps=2)0.65384312666612
10even-odd edge-groupedLieTrotter(reps=4)0.1815292413213224
11even-odd edge-groupedLieTrotter(reps=8)0.0452444826426448
\n", "
" ], "text/plain": [ " ordering synthesis infidelity depth size \\\n", "0 naive LieTrotter(reps=1) 0.999917 15 33 \n", "1 naive LieTrotter(reps=2) 0.805998 21 66 \n", "2 naive LieTrotter(reps=4) 0.271389 33 132 \n", "3 naive LieTrotter(reps=8) 0.074590 57 264 \n", "4 even-then-odd LieTrotter(reps=1) 0.999917 6 33 \n", "5 even-then-odd LieTrotter(reps=2) 0.805998 12 66 \n", "6 even-then-odd LieTrotter(reps=4) 0.271389 24 132 \n", "7 even-then-odd LieTrotter(reps=8) 0.074590 48 264 \n", "8 even-odd edge-grouped LieTrotter(reps=1) 0.998432 6 33 \n", "9 even-odd edge-grouped LieTrotter(reps=2) 0.653843 12 66 \n", "10 even-odd edge-grouped LieTrotter(reps=4) 0.181529 24 132 \n", "11 even-odd edge-grouped LieTrotter(reps=8) 0.045244 48 264 \n", "\n", " 2q gates 2q depth \n", "0 33 15 \n", "1 66 21 \n", "2 132 33 \n", "3 264 57 \n", "4 33 6 \n", "5 66 12 \n", "6 132 24 \n", "7 264 48 \n", "8 33 6 \n", "9 66 12 \n", "10 132 24 \n", "11 264 48 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ordering_comparison_rows = []\n", "infidelity_reps = [1, 2, 4, 8]\n", "\n", "for ordering, ham_ordered in hamiltonians_by_ordering.items():\n", " for reps in infidelity_reps:\n", " circuit = build_numeric_evolution_circuit(\n", " ham_ordered,\n", " LieTrotter,\n", " comparison_time,\n", " reps=reps,\n", " )\n", " decomposed_circuit = simple_transpilation(circuit)\n", " infidelity = process_infidelity(decomposed_circuit, exact_matrix)\n", " ordering_comparison_rows.append(\n", " {\n", " \"ordering\": ordering,\n", " \"synthesis\": f\"LieTrotter(reps={reps})\",\n", " \"infidelity\": infidelity,\n", " **summarize_circuit(decomposed_circuit),\n", " }\n", " )\n", " if reps == 1:\n", " print(f\"Circuit for {ordering} ordering:\")\n", " print([op for op, _ in ham_ordered.to_list()])\n", " display(decomposed_circuit.draw(\"mpl\", fold=-1, scale=0.6))\n", "\n", "ordering_comparison_df = pd.DataFrame(ordering_comparison_rows)\n", "display(ordering_comparison_df)\n", "\n", "hamiltonian_for_synthesis = hamiltonians_by_ordering[\"even-odd edge-grouped\"]" ] }, { "cell_type": "markdown", "id": "1c15299c", "metadata": {}, "source": [ "Wir stellen fest `reps=1`, dass sowohl die „even-then-odd“- als auch die „even-odd“-Kantengruppierung die Tiefe von 15 auf 6 reduzieren, indem sie parallele Zwei-Qubit-Schichten freilegen, während die Anzahl der Zwei-Qubit-Gatter bei allen drei Anordnungen gleich bleibt.\n", "Da jedoch bei allen drei Ordnungen die Untreue nahe bei 1 liegt, variieren wir die Anzahl der Trotter-Wiederholungen, um die Ordnungen deutlicher voneinander abzugrenzen.\n", "Mit zunehmender Anzahl von Wiederholungen nimmt die Unzuverlässigkeit der Reihenfolge `even-odd edge-grouped` schneller ab als bei den beiden anderen und erreicht den Wert 0.045, während `reps=8` sie bei der naiven sowie der „zuerst gerade, dann ungerade“-Reihenfolge bei 0.075 liegt.\n", "\n", "\n", "\n", "#### Vergleich der Synthese von Produktformeln\n", "\n", "Als Nächstes untersuchen wir verschiedene erweiterte Einstellungen der Trotterisierung, wobei die feste Hamilton-Ordnung auf die nach geraden und ungeraden Kanten gruppierte Ordnung umgestellt wird.\n", "Wir betrachten den Lie-Trotter-Prozess erster Ordnung, den Suzuki-Trotter-Prozess zweiter Ordnung und den Suzuki-Trotter-Prozess vierter Ordnung.\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "7b6b7f2f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
synthesisrepsinfidelitydepthsize2q gates2q depth
0LieTrotter19.984324e-01633336
1SuzukiTrotter(order=2)19.733399e-01951519
3LieTrotter26.538427e-0112666612
4SuzukiTrotter(order=2)22.522533e-0115848415
6LieTrotter33.197242e-0118999918
7SuzukiTrotter(order=2)34.804050e-022111711721
9LieTrotter41.815291e-012413213224
10SuzukiTrotter(order=2)41.453103e-022715015027
12LieTrotter51.161770e-013016516530
2SuzukiTrotter(order=4)12.884402e-013318318333
13SuzukiTrotter(order=2)55.803455e-033318318333
5SuzukiTrotter(order=4)21.641162e-036334834863
8SuzukiTrotter(order=4)32.907076e-059351351393
11SuzukiTrotter(order=4)42.791061e-06123678678123
14SuzukiTrotter(order=4)54.736685e-07153843843153
\n", "
" ], "text/plain": [ " synthesis reps infidelity depth size 2q gates \\\n", "0 LieTrotter 1 9.984324e-01 6 33 33 \n", "1 SuzukiTrotter(order=2) 1 9.733399e-01 9 51 51 \n", "3 LieTrotter 2 6.538427e-01 12 66 66 \n", "4 SuzukiTrotter(order=2) 2 2.522533e-01 15 84 84 \n", "6 LieTrotter 3 3.197242e-01 18 99 99 \n", "7 SuzukiTrotter(order=2) 3 4.804050e-02 21 117 117 \n", "9 LieTrotter 4 1.815291e-01 24 132 132 \n", "10 SuzukiTrotter(order=2) 4 1.453103e-02 27 150 150 \n", "12 LieTrotter 5 1.161770e-01 30 165 165 \n", "2 SuzukiTrotter(order=4) 1 2.884402e-01 33 183 183 \n", "13 SuzukiTrotter(order=2) 5 5.803455e-03 33 183 183 \n", "5 SuzukiTrotter(order=4) 2 1.641162e-03 63 348 348 \n", "8 SuzukiTrotter(order=4) 3 2.907076e-05 93 513 513 \n", "11 SuzukiTrotter(order=4) 4 2.791061e-06 123 678 678 \n", "14 SuzukiTrotter(order=4) 5 4.736685e-07 153 843 843 \n", "\n", " 2q depth \n", "0 6 \n", "1 9 \n", "3 12 \n", "4 15 \n", "6 18 \n", "7 21 \n", "9 24 \n", "10 27 \n", "12 30 \n", "2 33 \n", "13 33 \n", "5 63 \n", "8 93 \n", "11 123 \n", "14 153 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "synthesis_comparison_rows = []\n", "\n", "for num_trotter_steps in [1, 2, 3, 4, 5]:\n", " synthesis_cases = [\n", " (\"LieTrotter\", LieTrotter, {\"reps\": num_trotter_steps}),\n", " (\n", " \"SuzukiTrotter(order=2)\",\n", " SuzukiTrotter,\n", " {\"order\": 2, \"reps\": num_trotter_steps},\n", " ),\n", " (\n", " \"SuzukiTrotter(order=4)\",\n", " SuzukiTrotter,\n", " {\"order\": 4, \"reps\": num_trotter_steps},\n", " ),\n", " ]\n", " for label, synthesis, kwargs in synthesis_cases:\n", " circuit = build_numeric_evolution_circuit(\n", " hamiltonian_for_synthesis,\n", " synthesis,\n", " comparison_time,\n", " **kwargs,\n", " )\n", " decomposed_circuit = simple_transpilation(circuit)\n", " synthesis_comparison_rows.append(\n", " {\n", " \"synthesis\": label,\n", " \"reps\": kwargs[\"reps\"],\n", " \"infidelity\": process_infidelity(\n", " decomposed_circuit, exact_matrix\n", " ),\n", " **summarize_circuit(decomposed_circuit),\n", " }\n", " )\n", "\n", "synthesis_comparison_df = pd.DataFrame(synthesis_comparison_rows)\n", "display(synthesis_comparison_df.sort_values([\"2q gates\", \"2q depth\"]))\n", "\n", "# For memory free\n", "exact_matrix = None" ] }, { "cell_type": "markdown", "id": "304f3357", "metadata": {}, "source": [ "Lie-Trotter liefert den flachsten Kreis, weist jedoch den größten Fehler auf, während der Suzuki-Trotter-Algorithmus vierter Ordnung zwar genauer ist, aber die Kreistiefe erhöht.\n", "Für den weiteren Verlauf des Tutorials entscheiden wir uns für die Suzuki-Trotter-Formel zweiter Ordnung, da sie eine Schaltung mit geringer Tiefe liefert und gleichzeitig den Trotter-Fehler im Vergleich zur Formel erster Ordnung erheblich reduziert.\n", "\n" ] }, { "cell_type": "markdown", "id": "2415fa62", "metadata": {}, "source": [ "\n", "\n", "#### Das Gatter für die kontrollierte zeitliche Entwicklung entfernen\n", "\n", "Im erweiterten Swap-Test erfordert die Steuerung der zeitlichen Entwicklung durch ein einzelnes Ancilla-Qubit, dass das Ancilla viele Gatter im gesamten System steuert.\n", "Dies kann zu einem erheblichen Routing-Aufwand führen und erfordert im schlimmsten Fall praktisch eine „All-to-One“-Konnektivität.\n", "Um dies zu vermeiden, ist eine weitere Optimierung möglich, indem das gesteuerte Zeitentwicklungsgatter durch eine steuerungsfreie Version ersetzt wird, wobei die Symmetrie des Hamilton-Operators ausgenutzt wird.\n", "Betrachten wir die folgende Schaltung.\n", "\n", "![Schaltungsoptimierung](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-quantum-diagonalization/circuit_opt.avif)\n", "\n", "Hier bereitet $B_{\\rm ref}$ den Referenzzustand vor: $B_{\\rm ref}|0^n\\rangle = |\\psi_{0}\\rangle$.\n", "\n", "Anstatt zunächst „ $\\frac{1}{\\sqrt{2}}(|0\\rangle+|1\\rangle)|\\psi_0\\rangle$ “ vorzubereiten und anschließend „ $U(t)$ “ nur auf den Zweig „ $|1\\rangle$ “ anzuwenden, bereitet der Schaltkreis die beiden Zweige direkt wie folgt vor:\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle|\\psi_0\\rangle\n", "\\quad \\text{and} \\quad\n", "|1\\rangle|\\psi_d\\rangle ,\n", "\\end{equation*}\n", "$$\n", "\n", "wo\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\psi_d\\rangle = U(d\\Delta t)|\\psi_0\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Die Schaltung wendet zunächst ein Hadamard-Gatter auf den Ancilla an und bereitet den Referenzzustand nur auf dem Zweig „ $|1\\rangle$ “ vor:\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle |0^n\\rangle\n", "\\longrightarrow\n", "\\frac{|0\\rangle |0^n\\rangle + |1\\rangle B_{\\rm ref}|0^n\\rangle}{\\sqrt{2}}\n", "=\n", "\\frac{|0\\rangle |0^n\\rangle + |1\\rangle |\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Anschließend wird der Operator der unkontrollierten Zeitentwicklung auf beide Zweige angewendet:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{|0\\rangle |0^n\\rangle + |1\\rangle |\\psi_0\\rangle}{\\sqrt{2}}\n", "\\longrightarrow\n", "\\frac{|0\\rangle U(t)|0^n\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Da der Hamiltonoperator die Anregungszahl erhält, lässt sich erkennen, dass $|0^n\\rangle$ sein Eigenzustand ist, und somit\n", "akumuliert der Evolutionsoperator unter dem Hamiltonoperator lediglich eine Phase:\n", "\n", "$$\n", "\\begin{equation*}\n", "U(t)|0^n\\rangle\n", "=\n", "e^{-iE_{\\rm vac}t}|0^n\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Dies bedeutet Folgendes:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{|0\\rangle U(t)|0^n\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}}\n", "=\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |0^n\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Anschließend wird „ $B_{\\rm ref}$ “ ausschließlich auf den Zweig „ $|0\\rangle$ “ angewendet.\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |\\psi_0\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Zu diesem Zeitpunkt weisen die beiden Zweige eine zusätzliche relative Phase auf. Um es zu entfernen, wenden wir das Ancilla-Phasengatter an\n", "\n", "$$\n", "\\begin{equation*}\n", "P(-E_{\\rm vac}t)\n", "=\n", "\\begin{pmatrix}\n", "1 & 0 \\\\\n", "0 & e^{-iE_{\\rm vac}t}\n", "\\end{pmatrix}.\n", "\\end{equation*}\n", "$$\n", "\n", "Dadurch wird der Zustand wie folgt transformiert:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |\\psi_0\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}}\n", "\\longrightarrow\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |\\psi_0\\rangle\n", "+\n", "e^{-iE_{\\rm vac}t}|1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "So haben wir schließlich, bis auf eine irrelevante globale Phase, Folgendes hergestellt:\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\Phi_{0d}\\rangle\n", "=\n", "\\frac{\n", "|0\\rangle|\\psi_0\\rangle\n", "+\n", "|1\\rangle U(t)|\\psi_0\\rangle\n", "}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Im folgenden Code-Block implementieren wir diese steuerungsfreie Schaltung.\n", "\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "790e5e0c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Vacuum energy E_vac = 11.0+0.0j\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "controlled_extended_swap_test = extended_swap_test\n", "\n", "# Reuse the ordering and synthesis rule selected by the comparison above.\n", "evol_gate_optimized = PauliEvolutionGate(\n", " hamiltonian_for_synthesis,\n", " time=t,\n", " synthesis=SuzukiTrotter(order=2, reps=num_trotter_steps),\n", ")\n", "\n", "uncontrolled_evolution = QuantumCircuit(num_qubits, name=\"U_ST2(t)\")\n", "uncontrolled_evolution.append(evol_gate_optimized, range(num_qubits))\n", "\n", "controlled_state_prep = QuantumCircuit(num_qubits + 1, name=\"C-Prep\")\n", "for q, bit in enumerate(reversed(ref_bitstring)):\n", " if bit == \"1\":\n", " controlled_state_prep.cx(ancilla, system_qubits[q])\n", "\n", "vacuum_bitstring = \"0\" * num_qubits\n", "vacuum_energy = basis_state_expectation_sparse(hamiltonian, \"0\" * num_qubits)\n", "\n", "optimized_extended_swap_test = QuantumCircuit(num_qubits + 1)\n", "optimized_extended_swap_test.h(ancilla)\n", "\n", "# Prepare |psi_ref> only on the |1> branch.\n", "optimized_extended_swap_test.compose(controlled_state_prep, inplace=True)\n", "optimized_extended_swap_test.barrier()\n", "\n", "# Apply the Trotterized time evolution without control.\n", "optimized_extended_swap_test.compose(\n", " uncontrolled_evolution,\n", " qubits=system_qubits,\n", " inplace=True,\n", ")\n", "optimized_extended_swap_test.barrier()\n", "\n", "# Map |0>|0...0> to |0>|psi_ref>, leaving the |1> branch unchanged.\n", "optimized_extended_swap_test.x(ancilla)\n", "optimized_extended_swap_test.compose(controlled_state_prep, inplace=True)\n", "optimized_extended_swap_test.x(ancilla)\n", "\n", "# Cancel the known vacuum phase so that the same X/Y observables can be used.\n", "optimized_extended_swap_test.p(-vacuum_energy * t, ancilla)\n", "optimized_extended_swap_test = simple_transpilation(\n", " optimized_extended_swap_test\n", ")\n", "\n", "print(f\"Vacuum energy E_vac = {vacuum_energy:.1f}\")\n", "display(\n", " optimized_extended_swap_test.assign_parameters({t: 1.0}).draw(\n", " \"mpl\", scale=0.5, fold=26\n", " )\n", ")" ] }, { "cell_type": "markdown", "id": "e223d1a1", "metadata": {}, "source": [ "\n", "\n", "#### Transpilation\n", "\n", "Nun transpilieren wir die gesteuerten und die ungesteuerten Schaltungen, damit sie auf der Hardware ausgeführt werden können.\n", "Vergleichen wir das Ergebnis der transpilierten Schaltungen.\n", "\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "6750197b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
labeldepthsize2q gates2q depth
0abstract controlled U(t)154572378646864580
1optimized non-controlled U(t)261171630757
\n", "
" ], "text/plain": [ " label depth size 2q gates 2q depth\n", "0 abstract controlled U(t) 15457 23786 4686 4580\n", "1 optimized non-controlled U(t) 261 1716 307 57" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Layout({\n", "18: ,\n", "5: ,\n", "6: ,\n", "7: ,\n", "8: ,\n", "9: ,\n", "10: ,\n", "11: ,\n", "12: ,\n", "13: ,\n", "14: ,\n", "15: ,\n", "19: \n", "})\n" ] } ], "source": [ "pass_manager = generate_preset_pass_manager(\n", " backend=backend,\n", " optimization_level=3,\n", ")\n", "\n", "isa_controlled_extended_swap_test = pass_manager.run(\n", " controlled_extended_swap_test\n", ")\n", "\n", "pass_manager = generate_preset_pass_manager(\n", " backend=backend, optimization_level=3, routing_method=\"none\"\n", ")\n", "isa_optimized_extended_swap_test = pass_manager.run(\n", " optimized_extended_swap_test\n", ")\n", "\n", "transpilation_result = [\n", " {\n", " \"label\": \"abstract controlled U(t)\",\n", " **summarize_circuit(isa_controlled_extended_swap_test),\n", " },\n", " {\n", " \"label\": \"optimized non-controlled U(t)\",\n", " **summarize_circuit(isa_optimized_extended_swap_test),\n", " },\n", "]\n", "\n", "display(pd.DataFrame(transpilation_result))\n", "\n", "\n", "def filter_qubits_from_layout(layout):\n", " q_layout = Layout(\n", " {\n", " physical: virtual\n", " for physical, virtual in layout.get_physical_bits().items()\n", " if virtual._register.name == \"q\"\n", " }\n", " )\n", " return q_layout\n", "\n", "\n", "print(\n", " filter_qubits_from_layout(\n", " isa_optimized_extended_swap_test.layout.initial_layout\n", " )\n", ")\n", "\n", "isa_observables = [\n", " op.apply_layout(isa_optimized_extended_swap_test.layout)\n", " for op in observables\n", "]" ] }, { "cell_type": "markdown", "id": "39aef2af", "metadata": {}, "source": [ "\n", "\n", "### Schritt 3: Mit „ Qiskit primitives “ ausführen\n", "\n", "Der nächste Schritt besteht darin, dieselbe parametrisierte Schaltung für mehrere Werte von $d$ zu betrachten. Für jeden Wert von $d$ schätzen wir vier Erwartungswerte: $X\\otimes I$, $Y\\otimes I$, $X\\otimes (H-T)$ und $Y\\otimes (H-T)$. Diese vier Zahlen werden dann zu den komplexen Elementen der ersten Zeile $\\mathcal{S}_{0d}$ und $\\mathcal{H}_{0d}$ zusammengefasst.\n", "\n", "Hier lassen sich $\\mathcal{S}_{00}=1$ und $\\mathcal{H}_{00}=\\langle\\psi_{0}|(H-T)|\\psi_{0}\\rangle=\\langle\\psi_{0}|H|\\psi_{0}\\rangle-\\tau$ klassisch berechnen, da $\\lvert\\psi_0\\rangle$ spärlich ist; daher überspringen wir den Fall $d=0$.\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "83a49027", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of Krylov basis states: r = 9\n", "Number of PUBs: 36\n", "Each d uses observables: ['Re S_0d', 'Im S_0d', 'Re shifted H_0d', 'Im shifted H_0d']\n" ] } ], "source": [ "pub_list = []\n", "d_values = list(range(1, krylov_dim))\n", "\n", "# Exact local statevector estimator.\n", "estimator = StatevectorEstimator()\n", "\n", "# We use the circuit before the transpilation for the local simulator,\n", "# but we will use the transpiled circuit for the real backend.\n", "for d in d_values:\n", " parameter_values = [d * dt]\n", "\n", " for ob in observables:\n", " pub_list.append(\n", " (\n", " optimized_extended_swap_test,\n", " ob,\n", " parameter_values,\n", " )\n", " )\n", "\n", "job = estimator.run(pub_list)\n", "\n", "# Local PrimitiveJob does not provide Runtime-style job inputs,\n", "# so preserve the inputs directly.\n", "inputs = pub_list\n", "result = job.result()\n", "\n", "print(f\"Number of Krylov basis states: r = {len(d_values)}\")\n", "print(f\"Number of PUBs: {len(pub_list)}\")\n", "print(f\"Each d uses observables: {observable_labels}\")" ] }, { "cell_type": "markdown", "id": "a0fd2bc9", "metadata": {}, "source": [ "\n", "\n", "### Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format\n", "\n", "Nach der Schätzung der projizierten Matrizen führen wir eine Regularisierung durch und lösen das GEVP\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{H}\\mathbf{c} = E\\mathcal{S}\\mathbf{c}.\n", "\\end{equation*}\n", "$$\n", "\n", "Der kleinste verallgemeinerte Eigenwert liefert die KQD-Schätzung der Grundzustandsenergie.\n", "\n" ] }, { "cell_type": "code", "execution_count": 16, "id": "1019edb3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated first row of S:\n", "[ 1. +0.j 0.758-0.596j 0.203-0.836j -0.291-0.636j -0.444-0.229j\n", " -0.276+0.053j -0.044+0.051j 0.006-0.121j -0.155-0.218j -0.345-0.101j]\n", "\n", "Estimated first row of H:\n", "[ 7. +0.j 4.842-4.76j 0.044-6.185j -3.791-3.653j -4.137+0.39j\n", " -1.495+2.654j 1.331+1.777j 1.85 -0.763j -0.032-2.278j -2.222-1.372j]\n", "\n", "Eigenvalues of the estimated overlap matrix S:\n", "[-0. 0. 0. 0. 0. 0. 0.01 0.355 3.526 6.109]\n", "Condition number above 1e-12: 1.904e+12\n", "\n" ] } ], "source": [ "h_shifted_row_est = np.zeros(krylov_dim, dtype=complex)\n", "s_row_est = np.zeros(krylov_dim, dtype=complex)\n", "\n", "h_shifted_row_est[0] = ref_energy - shift_tau\n", "s_row_est[0] = 1.0\n", "\n", "for idx, (pub_input, pub_result) in enumerate(zip(inputs, result)):\n", " d_index, obs_index = divmod(idx, len(observables))\n", " ev = np.asarray(pub_result.data.evs).reshape(-1)[0]\n", " std = np.asarray(pub_result.data.stds).reshape(-1)[0]\n", "\n", " if obs_index == 0:\n", " s_row_est[d_index + 1] = ev\n", " elif obs_index == 1:\n", " s_row_est[d_index + 1] += 1j * ev\n", " elif obs_index == 2:\n", " h_shifted_row_est[d_index + 1] = ev\n", " elif obs_index == 3:\n", " h_shifted_row_est[d_index + 1] += 1j * ev\n", "\n", "# H_0d = shifted_H_0d + tau * S_0d.\n", "h_row_est = h_shifted_row_est + shift_tau * s_row_est\n", "\n", "h_matrix_est = la.toeplitz(h_row_est.conj(), h_row_est)\n", "s_matrix_est = la.toeplitz(s_row_est.conj(), s_row_est)\n", "\n", "s_eigvals = la.eigvalsh(0.5 * (s_matrix_est + s_matrix_est.conj().T))\n", "positive_s_eigvals = s_eigvals[s_eigvals > 1e-12]\n", "s_condition_number = (\n", " positive_s_eigvals[-1] / positive_s_eigvals[0]\n", " if len(positive_s_eigvals) > 0\n", " else np.inf\n", ")\n", "\n", "\n", "with np.printoptions(precision=3, suppress=True):\n", " print(\"Estimated first row of S:\")\n", " print(s_row_est)\n", " print()\n", " print(\"Estimated first row of H:\")\n", " print(h_row_est)\n", " print()\n", " print(\"Eigenvalues of the estimated overlap matrix S:\")\n", " print(s_eigvals)\n", " print(f\"Condition number above 1e-12: {s_condition_number:.3e}\")\n", " print()" ] }, { "cell_type": "markdown", "id": "176d4283", "metadata": {}, "source": [ "Nun lösen wir das verallgemeinerte Eigenwertproblem unter Verwendung der Matrizen, die aus den Schaltungsschätzungen rekonstruiert wurden.\n", "Bei einer idealen Zustandsvektorberechnung mit exakter Echtzeitentwicklung sollte dies das exakte projizierte Ergebnis wiedergeben.\n", "In der Praxis können Abweichungen durch Trotterisierung, Stichprobenfehler und numerische Instabilität der Überlappungsmatrix entstehen.\n", "\n", "Wir beobachten, wie sich die Energie annähert, wenn wir die Dimension des Krylov-Unterraums erhöhen.\n", "\n" ] }, { "cell_type": "code", "execution_count": 17, "id": "e789b2e0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "exact ground state energy: 3.136296694843727\n", "Krylov ground state energy (dim=1, retained=1): 7.0\n", "Krylov ground state energy (dim=2, retained=2): 4.184510657551266\n", "Krylov ground state energy (dim=3, retained=3): 3.5539074630394136\n", "Krylov ground state energy (dim=4, retained=4): 3.3366270761341044\n", "Krylov ground state energy (dim=5, retained=5): 3.252017453225087\n", "Krylov ground state energy (dim=6, retained=6): 3.2300275138879186\n", "Krylov ground state energy (dim=7, retained=7): 3.2299154099085685\n", "Krylov ground state energy (dim=8, retained=7): 3.2298063744216776\n", "Krylov ground state energy (dim=9, retained=7): 3.2296778282872456\n", "Krylov ground state energy (dim=10, retained=8): 3.223647515867734\n" ] } ], "source": [ "exact_evals = diagonalize_single_1_subspace(hamiltonian)\n", "exact_ground = min(exact_evals)\n", "print(\"exact ground state energy: \", exact_ground)\n", "\n", "threshold = 1e-12\n", "energy_convergence = []\n", "\n", "for r in range(1, krylov_dim + 1):\n", " energy_est_kqd, coeffs_est, retained_est = solve_thresholded_gevp(\n", " h_matrix_est[:r, :r],\n", " s_matrix_est[:r, :r],\n", " threshold=threshold,\n", " )\n", " energy_convergence.append(energy_est_kqd)\n", " print(\n", " f\"Krylov ground state energy (dim={r}, retained={retained_est}): \",\n", " energy_est_kqd,\n", " )" ] }, { "cell_type": "code", "execution_count": 18, "id": "08d57f33", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def plot_energy_convergence(energy_convergence, exact_ground, krylov_dim):\n", " fig, ax = plt.subplots(figsize=(7, 4.5))\n", "\n", " ax.plot(\n", " range(1, krylov_dim + 1),\n", " energy_convergence,\n", " marker=\"o\",\n", " label=\"KQD estimate\",\n", " )\n", "\n", " ax.axhline(\n", " exact_ground,\n", " linestyle=\"--\",\n", " label=f\"Exact ground energy = {exact_ground:.6f}\",\n", " )\n", "\n", " ax.set_xlabel(\"Krylov dimension\")\n", " ax.set_ylabel(\"Ground-state energy\")\n", " ax.grid(True, alpha=0.3)\n", " ax.legend()\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "plot_energy_convergence(energy_convergence, exact_ground, krylov_dim)" ] }, { "cell_type": "markdown", "id": "23c4eb47", "metadata": {}, "source": [ "\n", "\n", "## Beispiel für groß angelegte Hardware\n", "\n", "Im vorangegangenen Abschnitt wurde ein 12-Qubit-Modell verwendet, damit die Zustandsvektor-Simulation als Diagnoseinstrument dienen konnte. Nun skalieren wir denselben KQD-Workflow auf eine Heisenberg-Kette mit 30 Qubits und bereiten die Arbeitslast für die Ausführung auf der „ IBM Quantum “-Hardware vor.\n", "\n", "\n", "\n", "### Schritte 1–4 in einem einzigen Codeblock zusammengefasst\n", "\n", "Hier führen wir nun all diese Details zu einem einzigen Workflow in größerem Maßstab zusammen, der anschließend auf unserer realen Quantenhardware ausgeführt wird. In diesem Abschnitt wenden wir realistische Einstellungen zur Fehlerminderung an, um die Zuverlässigkeit der Ergebnisse zu verbessern. Da die Matrixelemente, die verschiedenen Werten von $d$ entsprechen, parallel berechnet werden können, nutzen wir den Modus mode `Batch`, um sie effizient auszuführen.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e3e0d1f3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Backend: ibm_boston\n", "System qubits: 30\n", "Total circuit qubits: 156\n", "Krylov dimension: 7\n", "Time step: 0.036110\n", "Shift tau: 25.0\n", "ISA circuit depth: 219\n", "ISA two-qubit gates: 728\n", "ISA two-qubit depth: 52\n", "Job IDs by d: {1: 'd9fl4p4jeosc73fk4dg0', 2: 'd9fl4pineu4c739poecg', 3: 'd9fl4q2neu4c739poedg', 4: 'd9fl4qhhtsac739fjgn0', 5: 'd9fl4r4jeosc73fk4dhg', 6: 'd9fl4rkjeosc73fk4dj0'}\n", "Estimated first row of S:\n", "[ 1. +0.j 0.42055-0.65466j -0.17883-0.73125j -0.64523-0.31601j\n", " -0.70781+0.55574j -0.02467+0.70525j 0.59301+0.55602j]\n", "Estimated first row of H:\n", "[ 25. +0.j 10.3419 -16.5586j -5.017 -18.03933j\n", " -16.44797 -6.98173j -17.17081+14.84502j 0.62722+17.81196j\n", " 15.68886+12.74875j]\n", "exact ground state energy: 21.021912418526902\n", "Krylov ground state energy (dim=1, retained=1): 25.0\n", "Krylov ground state energy (dim=2, retained=2): 24.432409110686205\n", "Krylov ground state energy (dim=3, retained=3): 24.256856893278556\n", "Krylov ground state energy (dim=4, retained=4): 23.727409826799715\n", "Krylov ground state energy (dim=5, retained=4): 23.324720470780864\n", "Krylov ground state energy (dim=6, retained=5): 21.91957579005085\n", "Krylov ground state energy (dim=7, retained=5): 21.548331214122644\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# -------------------------Step 1-------------------------\n", "\n", "# Map the classical problem to quantum circuits and observables.\n", "# Problem and KQD parameters.\n", "large_num_qubits = 30\n", "large_krylov_dim = 7\n", "large_num_trotter_steps = 3\n", "large_dt = np.pi / (3 * (large_num_qubits - 1))\n", "large_t = Parameter(\"t_large\")\n", "\n", "# Use a single excitation near the center of the chain.\n", "large_excitation_qubit = large_num_qubits // 2\n", "large_ref_label = [\"0\"] * large_num_qubits\n", "large_ref_label[large_num_qubits - 1 - large_excitation_qubit] = \"1\"\n", "large_ref_bitstring = \"\".join(large_ref_label)\n", "\n", "# Use the ordering and product formula selected in the preceding section.\n", "large_hamiltonian = make_heisenberg_hamiltonian_ordered(\n", " large_num_qubits,\n", " ordering=\"even-odd edge-grouped\",\n", " coupling=1.0,\n", ")\n", "large_ref_energy = float(\n", " np.real(\n", " basis_state_expectation_sparse(large_hamiltonian, large_ref_bitstring)\n", " )\n", ")\n", "large_vacuum_energy = float(\n", " np.real(\n", " basis_state_expectation_sparse(\n", " large_hamiltonian, \"0\" * large_num_qubits\n", " )\n", " )\n", ")\n", "\n", "large_evolution_gate = PauliEvolutionGate(\n", " large_hamiltonian,\n", " time=large_t,\n", " synthesis=SuzukiTrotter(order=2, reps=large_num_trotter_steps),\n", ")\n", "large_uncontrolled_evolution = QuantumCircuit(\n", " large_num_qubits,\n", " name=\"U_ST2_large(t)\",\n", ")\n", "large_uncontrolled_evolution.append(\n", " large_evolution_gate,\n", " range(large_num_qubits),\n", ")\n", "\n", "# Build the control-free extended-swap-test circuit.\n", "large_ancilla = 0\n", "large_system_qubits = list(range(1, large_num_qubits + 1))\n", "\n", "large_controlled_state_prep = QuantumCircuit(\n", " large_num_qubits + 1,\n", " name=\"C-Prep-large\",\n", ")\n", "large_controlled_state_prep.cx(\n", " large_ancilla,\n", " large_system_qubits[large_excitation_qubit],\n", ")\n", "\n", "large_extended_swap_test = QuantumCircuit(large_num_qubits + 1)\n", "large_extended_swap_test.h(large_ancilla)\n", "large_extended_swap_test.compose(large_controlled_state_prep, inplace=True)\n", "large_extended_swap_test.compose(\n", " large_uncontrolled_evolution,\n", " qubits=large_system_qubits,\n", " inplace=True,\n", ")\n", "large_extended_swap_test.x(large_ancilla)\n", "large_extended_swap_test.compose(\n", " large_controlled_state_prep.inverse(),\n", " inplace=True,\n", ")\n", "large_extended_swap_test.x(large_ancilla)\n", "large_extended_swap_test.p(\n", " -large_vacuum_energy * large_t,\n", " large_ancilla,\n", ")\n", "\n", "# Reuse the Hamiltonian-shifting construction from the preceding section.\n", "(\n", " large_obs_x_shifted_hamiltonian,\n", " large_obs_y_shifted_hamiltonian,\n", " large_shift_tau,\n", ") = make_reduced_heisenberg_observables(large_ref_bitstring)\n", "\n", "large_observables = [\n", " SparsePauliOp(\"I\" * large_num_qubits + \"X\"),\n", " SparsePauliOp(\"I\" * large_num_qubits + \"Y\"),\n", " large_obs_x_shifted_hamiltonian,\n", " large_obs_y_shifted_hamiltonian,\n", "]\n", "large_observable_labels = [\n", " \"Re S_0d\",\n", " \"Im S_0d\",\n", " \"Re shifted H_0d\",\n", " \"Im shifted H_0d\",\n", "]\n", "\n", "# -------------------------Step 2-------------------------\n", "\n", "# Optimize the problem for quantum execution.\n", "# Select a real backend and transpile the parameterized circuit to ISA form.\n", "large_backend = service.backend(\"ibm_boston\")\n", "large_pass_manager = generate_preset_pass_manager(\n", " backend=large_backend, optimization_level=3, routing_method=\"none\"\n", ")\n", "large_isa_circuit = large_pass_manager.run(large_extended_swap_test)\n", "large_isa_observables = [\n", " observable.apply_layout(large_isa_circuit.layout)\n", " for observable in large_observables\n", "]\n", "\n", "large_two_qubit_gate_count = sum(\n", " instruction.operation.num_qubits == 2\n", " for instruction in large_isa_circuit.data\n", ")\n", "\n", "print(f\"Backend: {large_backend.name}\")\n", "print(f\"System qubits: {large_num_qubits}\")\n", "print(f\"Total circuit qubits: {large_isa_circuit.num_qubits}\")\n", "print(f\"Krylov dimension: {large_krylov_dim}\")\n", "print(f\"Time step: {large_dt:.6f}\")\n", "print(f\"Shift tau: {large_shift_tau}\")\n", "print(f\"ISA circuit depth: {large_isa_circuit.depth()}\")\n", "print(f\"ISA two-qubit gates: {large_two_qubit_gate_count}\")\n", "print(\n", " f\"ISA two-qubit depth: {large_isa_circuit.depth(lambda x: x[0].num_qubits == 2)}\"\n", ")\n", "\n", "# -------------------------Step 3-------------------------\n", "\n", "# Execute on quantum hardware with Qiskit Runtime primitives.\n", "# Submit one job per d, with all four observables in that job.\n", "large_d_values = list(range(1, large_krylov_dim))\n", "retrieve_batch_id = None\n", "large_jobs = []\n", "\n", "if retrieve_batch_id is None:\n", " large_estimator_options = {\n", " \"default_shots\": 8192,\n", " \"dynamical_decoupling\": {\n", " \"enable\": True,\n", " \"sequence_type\": \"XpXm\",\n", " },\n", " \"resilience\": {\n", " \"measure_mitigation\": True,\n", " \"measure_noise_learning\": {\n", " \"num_randomizations\": 32,\n", " \"shots_per_randomization\": 256,\n", " },\n", " \"layer_noise_learning\": {\n", " \"max_layers_to_learn\": 4,\n", " \"layer_pair_depths\": [0, 1, 2, 4, 16, 32],\n", " \"num_randomizations\": 32,\n", " \"shots_per_randomization\": 128,\n", " },\n", " \"zne_mitigation\": True,\n", " \"zne\": {\n", " \"amplifier\": \"pea\",\n", " \"noise_factors\": [1.0, 1.5, 2.0],\n", " \"extrapolator\": (\"exponential\", \"linear\"),\n", " },\n", " },\n", " \"twirling\": {\n", " \"enable_gates\": True,\n", " \"enable_measure\": True,\n", " \"num_randomizations\": 32,\n", " \"shots_per_randomization\": 256,\n", " \"strategy\": \"active-accum\",\n", " },\n", " }\n", "\n", " with Batch(backend=large_backend) as large_batch:\n", " large_batch_id = large_batch.session_id\n", " large_estimator = EstimatorV2(\n", " mode=large_batch,\n", " options=large_estimator_options,\n", " )\n", " # Krylov quantum diagonalization of lattice Hamiltonians -> TUT_KQDOLH.\n", " large_estimator.options.environment.job_tags = [\"TUT_KQDOLH\"]\n", "\n", " for d in large_d_values:\n", " parameter_values = [d * large_dt]\n", " pubs_for_d = [\n", " (large_isa_circuit, observable, parameter_values)\n", " for observable in large_isa_observables\n", " ]\n", " job = large_estimator.run(pubs_for_d)\n", " large_jobs.append(job)\n", " print(\n", " f\"Submitted d={d}: job_id={job.job_id()}, \"\n", " f\"PUBs={len(pubs_for_d)}\"\n", " )\n", "\n", " print(f\"Batch ID: {large_batch_id}\")\n", "else:\n", " large_batch_id = retrieve_batch_id\n", " large_jobs = service.jobs(\n", " session_id=large_batch_id,\n", " limit=None,\n", " descending=False,\n", " )\n", "\n", "large_job_ids_by_d = {\n", " d: job.job_id() for d, job in zip(large_d_values, large_jobs)\n", "}\n", "print(f\"Job IDs by d: {large_job_ids_by_d}\")\n", "\n", "# -------------------------Step 4-------------------------\n", "\n", "# Post-process the quantum results and solve the classical GEVP.\n", "# Reconstruct the first rows of S and the shifted Hamiltonian matrix.\n", "large_s_row_est = np.zeros(large_krylov_dim, dtype=complex)\n", "large_h_shifted_row_est = np.zeros(large_krylov_dim, dtype=complex)\n", "large_s_row_est[0] = 1.0\n", "large_h_shifted_row_est[0] = large_ref_energy - large_shift_tau\n", "\n", "for d, job in zip(large_d_values, large_jobs):\n", " job_result = job.result()\n", "\n", " if len(job_result) != len(large_observables):\n", " raise RuntimeError(\n", " f\"Expected {len(large_observables)} PUB results for d={d}, \"\n", " f\"but received {len(job_result)}.\"\n", " )\n", "\n", " expectation_values = [\n", " np.asarray(pub_result.data.evs).reshape(-1)[0]\n", " for pub_result in job_result\n", " ]\n", "\n", " large_s_row_est[d] = expectation_values[0] + 1j * expectation_values[1]\n", " large_h_shifted_row_est[d] = (\n", " expectation_values[2] + 1j * expectation_values[3]\n", " )\n", "\n", "# H_0d = shifted_H_0d + tau * S_0d.\n", "large_h_row_est = large_h_shifted_row_est + large_shift_tau * large_s_row_est\n", "\n", "large_s_matrix_est = la.toeplitz(\n", " large_s_row_est.conj(),\n", " large_s_row_est,\n", ")\n", "large_h_matrix_est = la.toeplitz(\n", " large_h_row_est.conj(),\n", " large_h_row_est,\n", ")\n", "\n", "large_exact_gnd = min(diagonalize_single_1_subspace(large_hamiltonian))\n", "\n", "large_energy_convergence = []\n", "\n", "with np.printoptions(precision=5, suppress=True):\n", " print(\"Estimated first row of S:\")\n", " print(large_s_row_est)\n", " print(\"Estimated first row of H:\")\n", " print(large_h_row_est)\n", " print(\"exact ground state energy: \", large_exact_gnd)\n", "\n", " for r in range(1, large_krylov_dim + 1):\n", " energy_est_kqd, coeffs_est, retained_est = solve_thresholded_gevp(\n", " large_h_matrix_est[:r, :r],\n", " large_s_matrix_est[:r, :r],\n", " threshold=5e-2,\n", " )\n", " large_energy_convergence.append(energy_est_kqd)\n", " print(\n", " f\"Krylov ground state energy (dim={r}, retained={retained_est}): \",\n", " energy_est_kqd,\n", " )\n", "\n", "plot_energy_convergence(\n", " large_energy_convergence, large_exact_gnd, large_krylov_dim\n", ")" ] }, { "cell_type": "markdown", "id": "app-md-intro", "metadata": {}, "source": [ "\n", "\n", "## Anhang: Die Sichtweise der Hamilton-Funktion (Spektralfilter)\n", "\n", "Der vorgestellte Hauptablauf veranschaulicht die praktische Anwendung von KQD: Erstellung einer Krylov-Basis aus in Echtzeit berechneten Zuständen, Schätzung der projizierten Matrizen $\\mathcal{H}$ und $\\mathcal{S}$ sowie Lösung des GEVP.\n", "In diesem Anhang wird *dieselbe* Berechnung aus einem ergänzenden Blickwinkel betrachtet, der erklärt, **warum** KQD funktioniert: aus der Perspektive *der* Hamilton-Funktion bzw. *des Spektralfilters* [\\[3\\]](#references), [\\[5\\]](#references).\n", "Es greift auf das 12-Qubit-Modell, die Zeitschritt- $\\Delta t$ und die oben bereits ermittelte Krylov-Lösung zurück; eine erneute Schaltungsausführung ist nicht erforderlich.\n", "\n", "\n", "\n", "### Der Referenzzustand als Energieverteilung\n", "\n", "Sei die Eigenwertzerlegung des Hamilton-Operators\n", "\n", "$$\n", "\\begin{equation*}\n", "H=\\sum_{m} E_m\\,|E_m\\rangle\\!\\langle E_m|,\n", "\\qquad E_0\\le E_1\\le\\cdots,\n", "\\end{equation*}\n", "$$\n", "\n", "mit den Energieeigenzuständen $|E_m\\rangle$.\n", "Jeder Referenzzustand lässt sich in dieser Eigenbasis entwickeln,\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\psi_{0}\\rangle=\\sum_m a_m\\,|E_m\\rangle,\n", "\\qquad a_m=\\langle E_m|\\psi_{0}\\rangle,\n", "\\end{equation*}\n", "$$\n", "\n", "Daher weist sie bei jedem Energie $E_m$ ein **Spektralgewicht** $p_m=|a_m|^2$ auf.\n", "Die Referenzenergie ist der Mittelwert dieser Verteilung, $\\langle H\\rangle_{0}=\\sum_m p_m E_m$.\n", "\n", "Die Eigenwertzerlegung eines allgemeinen Hamilton-Operators vom Typ „ $n$ “-Qubit ist exponentiell aufwendig, daher dient dieses Bild lediglich der Diagnose und ist *nicht* Bestandteil des Algorithmus.\n", "Hier können wir dies jedoch für *dasselbe* Problem „ $n=12$ “ kostengünstig berechnen: Der Heisenberg-Hamiltonoperator bewahrt die Gesamtanzahl der Anregungen, und die Referenzverteilung $|000001000000\\rangle$ weist eine einzige Anregung auf, sodass ihr gesamter Spektralinhalt im Unterraum der einzelnen Anregung liegt, dessen Dimension nur linear mit $n$ wächst.\n", "Wir verwenden daher erneut den exakten Block für eine einzelne Anregung (der oben bereits als Benchmark verwendet wurde) und lesen die Referenzverteilung $\\{(E_m, p_m)\\}$ innerhalb dieses Unterraums ab.\n", "\n" ] }, { "cell_type": "code", "execution_count": 20, "id": "app-code-decomp", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Single-excitation subspace dimension: 12\n", "Subspace ground-state energy: 3.136297\n", "Reference energy (sum p_m E_m): 7.000000\n", "Reference weight on subspace ground: 0.163827\n" ] } ], "source": [ "# Exact single-excitation-subspace decomposition of the reference state.\n", "# This reuses make_heisenberg_hamiltonian / _basis_state_transition_amplitude_sparse\n", "# and the n=12 `hamiltonian` and `ref_bitstring` defined in the small-scale example.\n", "single_excitation_states = [1 << k for k in range(num_qubits)]\n", "\n", "h_single = np.array(\n", " [\n", " [\n", " _basis_state_transition_amplitude_sparse(hamiltonian, bra, ket)\n", " for ket in single_excitation_states\n", " ]\n", " for bra in single_excitation_states\n", " ]\n", ")\n", "h_single = 0.5 * (h_single + h_single.conj().T)\n", "\n", "subspace_evals, subspace_evecs = la.eigh(h_single)\n", "subspace_evals = np.real(subspace_evals)\n", "\n", "# Reference-state coordinates inside the single-excitation subspace.\n", "ref_position = single_excitation_states.index(int(ref_bitstring, 2))\n", "ref_in_subspace = np.zeros(num_qubits, dtype=complex)\n", "ref_in_subspace[ref_position] = 1.0\n", "\n", "# Amplitudes and spectral weights of the reference in the energy eigenbasis.\n", "ref_eigen_amplitudes = subspace_evecs.conj().T @ ref_in_subspace\n", "ref_spectral_weights = np.abs(ref_eigen_amplitudes) ** 2\n", "\n", "print(f\"Single-excitation subspace dimension: {num_qubits}\")\n", "print(f\"Subspace ground-state energy: {subspace_evals[0]:.6f}\")\n", "print(\n", " f\"Reference energy (sum p_m E_m): {np.sum(ref_spectral_weights * subspace_evals):.6f}\"\n", ")\n", "print(f\"Reference weight on subspace ground: {ref_spectral_weights[0]:.6f}\")" ] }, { "cell_type": "markdown", "id": "app-md-filter", "metadata": {}, "source": [ "\n", "\n", "### KQD trainiert einen Filter, der diese Verteilung umformt\n", "\n", "Eine **Hamilton-Funktion** $f(H)$ wird mittels Spektralrechnung definiert,\n", "\n", "$$\n", "\\begin{equation*}\n", "f(H)=\\sum_m f(E_m)\\,|E_m\\rangle\\!\\langle E_m|,\n", "\\end{equation*}\n", "$$\n", "\n", "oder, mit anderen Worten, eine gewichtete Summe von Eigenprojektoren.\n", "Wird dies auf die Referenz angewendet, wird jede Spektralamplitude umgestaltet, $a_m \\to a_m f(E_m)$ :\n", "\n", "$$\n", "\\begin{equation*}\n", "f(H)\\,|\\psi_{0}\\rangle=\\sum_m a_m\\,f(E_m)\\,|E_m\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Wenn „ $f$ “ bei der niedrigsten Energie eine scharfe Spitze aufweisen würde (ansonsten „ $f(E_0)=1$ “ und „ $f(E_m)\\approx 0$ “), würde „ $f(H)$ “ als Grundzustands-Projektor fungieren, und die normierte Ausgabe wäre (nahezu) der Grundzustand.\n", "Daher ist ein guter **energiespektraler** Tiefpassfilter genau das, was wir brauchen.\n", "\n", "KQD schreibt „ $f$ “ nicht im Voraus vor. Stattdessen erweitert es den Filter auf der Basis der Echtzeitentwicklung,\n", "\n", "$$\n", "\\begin{equation*}\n", "f_{\\rm KQD}(E)=\\sum_{\\ell=0}^{r-1} c_\\ell\\, e^{-i\\ell\\Delta t E},\n", "\\end{equation*}\n", "$$\n", "\n", "eine *trigonometrische* Funktion der Energie, deren Koeffizienten $\\{c_\\ell\\}$ genau dem oben ermittelten GEVP-Eigenvektor entsprechen.\n", "Die Minimierung des Rayleigh-Quotienten $\\mathbf{c}^\\dagger\\mathcal{H}\\mathbf{c}/\\mathbf{c}^\\dagger\\mathcal{S}\\mathbf{c}$ entspricht daher dem *Erlernen* des Filters, der das Gewicht des angeregten Zustands der Referenz am besten unterdrückt.\n", "Ein größeres Krylov- $r$ s verleiht dem Filter mehr Freiheitsgrade und einen schärferen Peak bei der Grundzustandsenergie.\n", "\n", "Der folgende Helper wertet diesen trainierten Filter auf einer Energieachse aus; anschließend wenden wir ihn auf die oben ermittelte Referenzverteilung an.\n", "\n" ] }, { "cell_type": "code", "execution_count": 21, "id": "app-code-filter", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ground-state overlap (reference): 0.1638\n", "Ground-state overlap (filtered): 0.9638\n", "Mean energy (reference): 7.000000\n", "Mean energy (filtered): 3.164503\n" ] } ], "source": [ "def trigonometric_krylov_filter(\n", " coeffs: np.ndarray,\n", " energies: np.ndarray,\n", " time_step: float,\n", ") -> np.ndarray:\n", " \"\"\"Evaluate the learned Krylov filter f(E) = sum_l c_l exp(-i l dt E).\"\"\"\n", " values = np.zeros_like(energies, dtype=complex)\n", " for ell, coeff in enumerate(coeffs):\n", " values += coeff * np.exp(-1j * ell * time_step * energies)\n", " return values\n", "\n", "\n", "def filtered_spectral_weights(\n", " weights: np.ndarray,\n", " filter_values: np.ndarray,\n", ") -> np.ndarray:\n", " \"\"\"Reshape spectral weights by |f(E)|^2 and renormalize.\"\"\"\n", " reshaped = weights * np.abs(filter_values) ** 2\n", " return reshaped / np.sum(reshaped)\n", "\n", "\n", "# Recover the KQD coefficients from the already-estimated projected matrices.\n", "# The shift only moves H by tau * S, so it does not change the GEVP eigenvector;\n", "# we solve at the full Krylov dimension used in the small-scale example.\n", "_, kqd_coeffs, _ = solve_thresholded_gevp(\n", " h_matrix_est,\n", " s_matrix_est,\n", " threshold=1e-12,\n", ")\n", "\n", "filter_on_spectrum = trigonometric_krylov_filter(\n", " kqd_coeffs, subspace_evals, dt\n", ")\n", "filtered_weights = filtered_spectral_weights(\n", " ref_spectral_weights, filter_on_spectrum\n", ")\n", "\n", "print(f\"Ground-state overlap (reference): {ref_spectral_weights[0]:.4f}\")\n", "print(f\"Ground-state overlap (filtered): {filtered_weights[0]:.4f}\")\n", "print(\n", " f\"Mean energy (reference): {np.sum(ref_spectral_weights * subspace_evals):.6f}\"\n", ")\n", "print(\n", " f\"Mean energy (filtered): {np.sum(filtered_weights * subspace_evals):.6f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "app-md-plot", "metadata": {}, "source": [ "\n", "\n", "### Stellen Sie sich den Filter und seine Flexibilität vor\n", "\n", "Zunächst zeigen wir den trainierten Filter bei der oben verwendeten vollen Krylov-Dimension und verfolgen dann, wie er mit zunehmender Dimension $r$ an Schärfe gewinnt.\n", "\n", "Die Balken zeigen die Referenz-Spektralgewichte $p_m$ (vorher) und die gefilterten Gewichte $p_m|f_{\\rm KQD}(E_m)|^2$ (nachher) zusammen mit der gelernten Filterintensität $|f_{\\rm KQD}(E)|^2$ auf einer kontinuierlichen Energieachse.\n", "Der Filter konzentriert das Gewicht auf die niedrigste Energie des Unterraums mit Einzelanregung – dieselbe Energie, zu der die KQD-Schätzung im Beispiel mit kleinem Maßstab konvergierte.\n", "Beachten Sie, dass es sich hierbei um den Grundzustand *innerhalb des Sektors mit einer einzigen Anregung* handelt, der für diesen anregungserhaltenden Referenzzustand relevant ist, und nicht um den globalen Grundzustand.\n", "\n" ] }, { "cell_type": "code", "execution_count": 22, "id": "app-code-plot", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(8, 4))\n", "\n", "visible = (ref_spectral_weights > 1e-4) | (filtered_weights > 1e-4)\n", "ax.bar(\n", " subspace_evals[visible],\n", " ref_spectral_weights[visible],\n", " width=0.18,\n", " alpha=0.45,\n", " label=\"reference $p_m$\",\n", ")\n", "ax.bar(\n", " subspace_evals[visible],\n", " filtered_weights[visible],\n", " width=0.14,\n", " alpha=0.9,\n", " label=r\"filtered $p_m\\,|f_{\\rm KQD}(E_m)|^2$\",\n", ")\n", "\n", "energy_grid = np.linspace(\n", " subspace_evals.min() - 0.5, subspace_evals.max() + 0.5, 800\n", ")\n", "filter_intensity = (\n", " np.abs(trigonometric_krylov_filter(kqd_coeffs, energy_grid, dt)) ** 2\n", ")\n", "filter_intensity /= filter_intensity.max()\n", "ax.plot(\n", " energy_grid,\n", " filter_intensity,\n", " color=\"k\",\n", " linewidth=2,\n", " label=r\"$|f_{\\rm KQD}(E)|^2$ (normalized)\",\n", ")\n", "\n", "ax.axvline(\n", " subspace_evals[0],\n", " color=\"C3\",\n", " linestyle=\"--\",\n", " linewidth=1,\n", " label=\"subspace ground energy\",\n", ")\n", "ax.set_xlabel(\"Energy eigenvalue $E_m$\")\n", "ax.set_ylabel(\"Spectral weight\")\n", "ax.set_ylim(0, 1)\n", "ax.legend(loc=\"upper right\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "app-md-sweep", "metadata": {}, "source": [ "\n", "\n", "### Erhöhung der Krylov-Dimension: Flexibilität der gelernten Funktion\n", "\n", "Zur Erinnerung: Der trainierte Filter ist ein trigonometrisches Polynom in der Energie mit Koeffizienten der Form $r$,\n", "\n", "$$\n", "\\begin{equation*}\n", "f_{\\rm KQD}(E)=\\sum_{\\ell=0}^{r-1} c_\\ell\\, e^{-i\\ell\\Delta t E}.\n", "\\end{equation*}\n", "$$\n", "\n", "Die Krylov-Dimension $r$ entspricht genau der Anzahl der freien Koeffizienten und bestimmt somit die **Flexibilität** der Funktion.\n", "Ein kleiner „ $r$ “ kann nur einen breiten, sanft variierenden Filter erzeugen, der Gewicht in tiefer liegende angeregte Zustände abfließen lässt; wenn die „ $r$ “ zunimmt, kann der Filter einen schmaleren Peak bei der Zielenergie bilden und das verbleibende Gewicht der angeregten Zustände stärker unterdrücken.\n", "Dies ist das spektrale Filter-Äquivalent der im Beispiel mit kleinen Skalen beobachteten Energiekonvergenz: Mit zunehmendem Wert von $r$ kollabiert die gefilterte Verteilung auf den Grundzustand des Unterraums, und die geschätzte Energie nähert sich diesem Wert an.\n", "\n", "Wir verwenden die oben bereits geschätzten Projektionsmatrizen erneut und lösen einfach das GEVP bei jedem führenden „ $r\\times r$ “-Block, um anschließend den entsprechenden Filter zu berechnen und grafisch darzustellen.\n", "\n" ] }, { "cell_type": "code", "execution_count": 23, "id": "app-code-sweep", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " r retained ground overlap filtered energy\n", "-- -------- -------------- ---------------\n", " 2 2 0.4549 4.184510\n", " 4 4 0.8373 3.310168\n", " 6 6 0.9706 3.158100\n", " 8 7 0.9701 3.158725\n", "10 8 0.9638 3.164503\n", "\n", "Subspace ground-state energy (target): 3.136297\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Sweep the Krylov dimension using the leading r x r blocks of the estimated matrices.\n", "sweep_dims = [r for r in (2, 4, 6, 8, krylov_dim) if r <= krylov_dim]\n", "sweep_dims = sorted(set(sweep_dims))\n", "\n", "energy_grid = np.linspace(\n", " subspace_evals.min() - 0.5, subspace_evals.max() + 0.5, 800\n", ")\n", "\n", "sweep_cases = []\n", "print(\" r retained ground overlap filtered energy\")\n", "print(\"-- -------- -------------- ---------------\")\n", "for r in sweep_dims:\n", " _, coeffs_r, retained_r = solve_thresholded_gevp(\n", " h_matrix_est[:r, :r],\n", " s_matrix_est[:r, :r],\n", " threshold=1e-12,\n", " )\n", " filter_on_spectrum_r = trigonometric_krylov_filter(\n", " coeffs_r, subspace_evals, dt\n", " )\n", " filtered_weights_r = filtered_spectral_weights(\n", " ref_spectral_weights, filter_on_spectrum_r\n", " )\n", " filtered_energy_r = float(np.sum(filtered_weights_r * subspace_evals))\n", "\n", " sweep_cases.append((r, coeffs_r, filtered_weights_r))\n", " print(\n", " f\"{r:2d} {retained_r:8d} {filtered_weights_r[0]:14.4f} {filtered_energy_r:15.6f}\"\n", " )\n", "\n", "print(f\"\\nSubspace ground-state energy (target): {subspace_evals[0]:.6f}\")\n", "\n", "# One panel per Krylov dimension: filtered spectrum (bars) + filter intensity (curve).\n", "fig, axes = plt.subplots(\n", " len(sweep_cases),\n", " 1,\n", " figsize=(8, 2.1 * len(sweep_cases)),\n", " sharex=True,\n", ")\n", "axes = np.atleast_1d(axes)\n", "\n", "for idx, (ax, (r, coeffs_r, filtered_weights_r)) in enumerate(\n", " zip(axes, sweep_cases)\n", "):\n", " visible = (ref_spectral_weights > 1e-4) | (filtered_weights_r > 1e-4)\n", " ax.bar(\n", " subspace_evals[visible],\n", " ref_spectral_weights[visible],\n", " width=0.18,\n", " alpha=0.35,\n", " color=\"C0\",\n", " label=\"reference $p_m$\" if idx == 0 else None,\n", " )\n", " ax.bar(\n", " subspace_evals[visible],\n", " filtered_weights_r[visible],\n", " width=0.14,\n", " alpha=0.9,\n", " color=\"C1\",\n", " label=\"filtered weights\" if idx == 0 else None,\n", " )\n", "\n", " filter_intensity_r = (\n", " np.abs(trigonometric_krylov_filter(coeffs_r, energy_grid, dt)) ** 2\n", " )\n", " filter_intensity_r /= filter_intensity_r.max()\n", " ax.plot(energy_grid, filter_intensity_r, color=\"k\", linewidth=2)\n", "\n", " ax.axvline(subspace_evals[0], color=\"C3\", linestyle=\"--\", linewidth=1)\n", " ax.set_ylim(0, 1)\n", " ax.set_ylabel(\"weight\")\n", " ax.legend(loc=\"upper right\", title=f\"$r={r}$\")\n", "\n", "axes[-1].set_xlabel(\"Energy eigenvalue $E_m$\")\n", "fig.suptitle(\n", " r\"KQD-learned filter $|f_{\\rm KQD}(E)|^2$ sharpening with Krylov dimension $r$\",\n", " y=1.0,\n", ")\n", "fig.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4a6b7ea8", "metadata": {}, "source": [ "\n", "\n", "## Nächste Schritte\n", "\n", "Wenn Sie diese Arbeit interessant fanden, könnten Sie auch an folgenden Materialien interessiert sein:\n", "\n", "\n", " * [Probenbasierte Krylov-Quantendiagonalisierung eines fermionischen Gittermodells](/docs/tutorials/sample-based-krylov-quantum-diagonalization)\n", " * [Schätzung der Grundzustandsenergie der Heisenberg-Kette mit VQE](/docs/tutorials/spin-chain-vqe)\n", " * [Kurs zu Algorithmen zur Quantendiagonalisierung](/learning/courses/quantum-diagonalization-algorithms/krylov)\n", "\n", "\n", "\n", "\n", "## Referenzen\n", "\n", "\\[1] E. N. Epperly, L. Lin und Y. Nakatsukasa, *Eine Theorie zur Diagonalisierung von Quantenunterräumen*, *SIAM Journal on Matrix Analysis and Applications* **43**, 1263–1290 (2022).\n", "\n", "\\[2] N. Yoshioka, M., Amico, W., Kirby u. a., *Diagonalisierung großer Vielteilchen-Hamilton-Operatoren auf einem Quantenprozessor*, [arXiv:2407.14431](https://arxiv.org/abs/2407.14431) (2024).\n", "\n", "\\[3] R. M. Parrish und P. L. McMahon, *Diagonalisierung von Quantenfiltern: Quanten-Eigenwertzerlegung ohne vollständige Quantphasenschätzung*, *Physical Review Letters* **122**, 230401 (2019).\n", "\n", "\\[4] G. Lee, S. Choi, J. Huh und AF Izmaylov, *Effiziente Strategien zur Reduzierung des Abtastfehlers bei der Quanten-Krylov-Unterraum-Diagonalisierung*, *Digital Discovery* **4**, 954-969 (2025).\n", "\n", "\\[5] G. Lee, M. Kang, J. Hong, S. Fomichev und J. Hm, *„Filtered Quantum Phase Estimation“*, [arXiv:2510.04294](https://arxiv.org/abs/2510.04294) (2025).\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "hours": 2, "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" }, "qpuSeconds": 4200 }, "nbformat": 4, "nbformat_minor": 5 }