Introduction to quantum science and technology - QUANT 400
\n",
"\n",
"
Lecturer:Prof. N. Macris
\n",
" \n",
"
Assistant: gian.gentinetta@epfl.ch"
]
},
{
"cell_type": "markdown",
"id": "46b0c9b1-3617-43e1-bfe5-2a76237204f9",
"metadata": {},
"source": [
"## Exercise 2 - Quantum computing"
]
},
{
"cell_type": "markdown",
"id": "caee92ab-871c-4b84-964f-4d1c2e75a463",
"metadata": {},
"source": [
"### Very brief introduction to Qiskit\n",
"\n",
"[Qiskit](https://www.ibm.com/quantum/qiskit) is an open-source software for working with quantum computers at the level of quantum circuits, operators, and primitives.\n",
"\n",
"See the [documentation](https://docs.quantum.ibm.com/) for more details.\n",
"\n",
"To install the required dependencies, run the cell below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b7c16b1c",
"metadata": {},
"outputs": [],
"source": [
"%pip install --upgrade pip\n",
"%pip install qiskit pylatexenc"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "93d8bb3f-10f2-49b8-a2df-ac17673694da",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"from qiskit import *\n",
"from qiskit.primitives import StatevectorEstimator, StatevectorSampler\n",
"from qiskit.quantum_info import Statevector, Pauli, SparsePauliOp"
]
},
{
"cell_type": "markdown",
"id": "e170b2eb-4fc0-4451-bbc3-4a0f37ac8b8b",
"metadata": {},
"source": [
"The basic element needed for your first program is the [QuantumCircuit](https://qiskit.org/documentation/stubs/qiskit.circuit.QuantumCircuit.html). We begin by creating a `QuantumCircuit` comprised of two qubits that start out in the $|00\\rangle$ state."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f27f574b-0af9-4a3a-9dce-82952baf59b0",
"metadata": {},
"outputs": [],
"source": [
"# Create a Quantum Circuit acting on a quantum register of two qubits\n",
"circ = QuantumCircuit(2)"
]
},
{
"cell_type": "markdown",
"id": "75e9aac3-8510-4a19-9514-dd9751d704ca",
"metadata": {},
"source": [
"Let's add gates that prepare a Bell-state (i.e. a maximally entangled two-qubit state):\n",
"\n",
"$$|\\psi\\rangle = \\left(|00\\rangle+|11\\rangle\\right)/\\sqrt{2}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "86ac0a79-d919-48d6-963a-e340a95614ef",
"metadata": {},
"outputs": [],
"source": [
"# Add a H gate on qubit 0, putting this qubit in superposition.\n",
"circ.h(0)\n",
"# Add a CX (CNOT) gate on control qubit 0 and target qubit 1, putting the qubits in a Bell state.\n",
"circ.cx(0, 1)"
]
},
{
"cell_type": "markdown",
"id": "9ecb37bd-b81a-49df-a022-a2d7e61bce26",
"metadata": {},
"source": [
"We can visualize the circuit"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "05d873c0-2667-4f90-b2a4-3c54af4d3b1f",
"metadata": {},
"outputs": [],
"source": [
"circ.draw(\"mpl\")"
]
},
{
"cell_type": "markdown",
"id": "f5c55df9-2ca5-4c5d-a242-f38f3ed29666",
"metadata": {},
"source": [
"For a full list of all gates that are available, check Qiskit's [circuit library](https://docs.quantum.ibm.com/api/qiskit/circuit_library)."
]
},
{
"cell_type": "markdown",
"id": "bb493210-e410-4b86-888d-15ec6d99e624",
"metadata": {},
"source": [
"We can extract the quantum state vector at the end of the circuit."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7cef7626-afea-4af0-8a37-33ed0b14816f",
"metadata": {},
"outputs": [],
"source": [
"Statevector(circ)"
]
},
{
"cell_type": "markdown",
"id": "2c10314f-aaaa-4a79-878d-d43a392a1a4a",
"metadata": {},
"source": [
"We can compute expectation values of observables using the `Estimator`. If you use the `StatevectorEstimator`, the expectation value is still calculated using the state vector representation of the quantum state."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d1767dc7-9157-480b-8fc8-895be74a6abe",
"metadata": {},
"outputs": [],
"source": [
"estimator = StatevectorEstimator()\n",
"observable = Pauli(\"XX\")\n",
"result = estimator.run([(circ, observable)]).result()[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "53b70be7-0565-4454-8289-fd75a3bdb847",
"metadata": {},
"outputs": [],
"source": [
"print(\"Value:\", result.data.evs)"
]
},
{
"cell_type": "markdown",
"id": "66094685",
"metadata": {},
"source": [
"You can compute the expectation values of multiple observables by providing a list of them."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2de15307-a8f1-40fd-9bc7-683ea1a841b9",
"metadata": {},
"outputs": [],
"source": [
"result = estimator.run([(circ, [Pauli(\"XX\"), Pauli(\"IY\")])]).result()[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ecba83cc-6019-4a2b-b0de-426231e7f486",
"metadata": {},
"outputs": [],
"source": [
"print(\"Value:\", result.data.evs)"
]
},
{
"cell_type": "markdown",
"id": "ae6a0a1d-38bb-4781-92c0-c393734ef1e4",
"metadata": {},
"source": [
"In practice on a real quantum computer, we can only measure the qubits in the computational basis (i.e., eigenstates of the Pauli-$\\hat{Z}$ operator) and obtain one of two possible outcomes for every qubit: 0, 1 corresponding to either the eigenvalue +1 (spin up) or eigenvalue -1 (spin down). The expectation values would then be computed by repeatedly running the circuit and measuring the qubits (possibly in different bases) and averaging over the measurement statistics. The more measurements (shots) you take, the more precise the expectation values can be estimated."
]
},
{
"cell_type": "markdown",
"id": "91e5510e-34b6-4b8e-8656-cc8dd8649bac",
"metadata": {},
"source": [
"We can obtain this raw measurement data using the `Sampler`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d25e4bfa-e6bc-4bc7-b45c-fe1556bd01a7",
"metadata": {},
"outputs": [],
"source": [
"circ.measure_all() # we first need to specify which qubits are measured\n",
"sampler = StatevectorSampler()\n",
"result = sampler.run([circ], shots=10000).result()[0]\n",
"counts = result.data.meas.get_counts()\n",
"counts"
]
},
{
"cell_type": "markdown",
"id": "ab3d981a",
"metadata": {},
"source": [
"Let's visualize the number of times each different bitstring was measured. Note that if you rerun the cell above, the outcome will change each time as the measurement process is inherently probabilistic."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e9781cd-f88d-4c00-aad1-033e1f9880b2",
"metadata": {},
"outputs": [],
"source": [
"from qiskit.visualization import plot_histogram\n",
"\n",
"plot_histogram(counts)"
]
},
{
"cell_type": "markdown",
"id": "f4b42d48-319d-41a5-bc32-5c5bef6e81f1",
"metadata": {},
"source": [
"### Time evolution of a quantum spin model"
]
},
{
"cell_type": "markdown",
"id": "04134078",
"metadata": {},
"source": [
"Let's simulate the time evolution of the Ising model given as \n",
"\n",
"$$\\hat H_{\\text{Ising}} = J \\sum_{i=1}^N \\hat Z_i \\hat Z_{i+1} - \\Gamma \\sum_{i=1}^N \\hat X_i$$\n",
"\n",
"using qiskit. We will use the 2nd order Trotter-Suzuki decomposition to approximate the total evolution unitary as\n",
"\n",
"$$ \\exp (-i \\hat{H} t) \\approx \\exp(-i \\frac{\\Delta_t}{2} \\hat H_\\mathrm{Z Z}) \\ldots \\exp(-i \\Delta_t \\hat H_\\mathrm{X}) \\exp(-i \\Delta_t \\hat H_\\mathrm{Z Z}) \\ldots \\exp(-i \\Delta_t \\hat H_\\mathrm{X}) \\exp(-i \\frac{\\Delta_t}{2} \\hat H_\\mathrm{Z Z}) $$\n",
"\n",
"In turn, each of these terms can be expressed as a product of commuting unitaries acting on one or two qubits only, e.g.,\n",
"\n",
"$$\\exp(-i \\Delta_t \\hat H_\\mathrm{X}) = \\prod_{i=1}^N \\exp(-i \\Delta_t \\Gamma \\hat{X}_i) = \\prod_{i=1}^N R_{X_i}(2 \\Delta_t \\Gamma ).$$\n",
"\n",
"$R_{X_i}$ denotes a rotation of qubit $i$ around the $X$ axis of the Bloch sphere. Note the extra factor of 2 since the definition of the $R_X$ gate in qiskit comes with a factor of 1/2."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "31fa1dd1-e96c-4d6e-8dbe-927373311dcd",
"metadata": {},
"outputs": [],
"source": [
"def time_evol_circuit(N, t, dt, J, gamma):\n",
" n_steps = round(t / dt) # number of Suzuki-Trotter steps\n",
" circ = QuantumCircuit(N)\n",
" for i in range(N - 1): # N-1 because we assume open boundary conditions\n",
" circ.rzz(2 * J * dt / 2, i, i + 1)\n",
" for step in range(n_steps - 1):\n",
" for i in range(N):\n",
" circ.rx(2 * gamma * dt, i)\n",
" for i in range(N - 1):\n",
" circ.rzz(2 * J * dt, i, i + 1)\n",
" for i in range(N):\n",
" circ.rx(2 * gamma * dt, i)\n",
" for i in range(N - 1):\n",
" circ.rzz(2 * J * dt / 2, i, i + 1)\n",
" return circ"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3fadac35-847c-4b31-a9a9-66a27aa47445",
"metadata": {},
"outputs": [],
"source": [
"circ = time_evol_circuit(4, 1.0, 0.1, 1.0, 0.5)\n",
"circ.draw(\"mpl\")"
]
},
{
"cell_type": "markdown",
"id": "f90b3c49-bcbc-4515-8bfc-735e51804a2e",
"metadata": {},
"source": [
"Starting from the $|0\\dots 0\\rangle$ initial state, let's measure the average magnetization along z-direction every 50 steps during the time evolution.\n",
"\n",
"$$ \\hat{O} = \\frac 1 N \\sum_i \\hat{Z}_i$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "19fb5d1d-70b8-492b-a969-8e5c17e12508",
"metadata": {},
"outputs": [],
"source": [
"N = 5 # number of qubits\n",
"t = 10.0 # total evolution time\n",
"dt = 1e-2 # Trotter time step size (Delta t)\n",
"J = 1.0 # Ising coupling strength\n",
"gamma = 0.8 # Ising transverse field\n",
"n_points = 50 # number of time stamps at which we want to measure the expectation value of the observable\n",
"ts = np.linspace(0, t, n_points + 1) # measure at these times\n",
"\n",
"# generate 50 different circuits corresponding to 50 different evolution times\n",
"circuits = [time_evol_circuit(N, t, dt, J, gamma) for t in ts]\n",
"\n",
"# define the observable we want to measure\n",
"obs_terms = [(\"Z\", [i], 1 / N) for i in range(N)]\n",
"observable = SparsePauliOp.from_sparse_list(obs_terms, num_qubits=N)\n",
"observable"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2bc0e159-c9a1-4dd4-a6b0-6bbbf8885418",
"metadata": {},
"outputs": [],
"source": [
"circuit_observable_pairs = [(circ, observable) for circ in circuits]\n",
"results = estimator.run(circuit_observable_pairs).result()\n",
"exp_values = [results[i].data.evs for i in range(n_points + 1)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "612b6f29",
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.plot(ts, exp_values)\n",
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(r\"$\\langle \\psi_t | \\hat{O} | \\psi_t \\rangle$\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "3fb9e495",
"metadata": {},
"source": [
"### Transpilation to native gates of the device"
]
},
{
"cell_type": "markdown",
"id": "eaede193",
"metadata": {},
"source": [
"When writing a quantum circuit you are free to use any quantum gate (unitary operator) that you like, alongside with a collection of non-gate operations such as qubit measurements and reset operations. However, when running a circuit on a real quantum device one no longer has this flexibility. Due to limitations in, for example, the physical interactions between qubits, difficulty in implementing multi-qubit gates, control electronics, etc., a quantum computing device can only natively support a handful of quantum gates and non-gate operations. In the present case of [IBM Q devices](https://quantum-computing.ibm.com/services/resources?tab=systems), the native gate set can be found by querying the devices themselves, and looking for the corresponding attribute in their configuration."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4757f690",
"metadata": {},
"outputs": [],
"source": [
"from qiskit.providers.fake_provider import GenericBackendV2\n",
"\n",
"fake_backend = GenericBackendV2(num_qubits=5, basis_gates=[\"u3\", \"cx\"])"
]
},
{
"cell_type": "markdown",
"id": "d4d84fcc",
"metadata": {},
"source": [
"Every quantum circuit run on an IBM Q device must be expressed using only these basis gates.\n",
"\n",
"Let's take the following circuit as an example"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2bc8b71a",
"metadata": {},
"outputs": [],
"source": [
"circ = time_evol_circuit(5, 0.1, 0.1, 1.0, 0.5)\n",
"circ.draw(output=\"mpl\")"
]
},
{
"cell_type": "markdown",
"id": "9ab2f210",
"metadata": {},
"source": [
"We have $R_{ZZ}$, and $R_X$ rotation gates, all of which are not in our devices basis gate set, and must be expanded. We can decompose the circuit to show what it would look like in the native gate set of the IBM Quantum devices."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4c202e52",
"metadata": {},
"outputs": [],
"source": [
"from qiskit import transpile\n",
"\n",
"circ_transpiled = transpile(circ, fake_backend)\n",
"circ_transpiled.draw(output=\"mpl\")"
]
},
{
"cell_type": "markdown",
"id": "f3cfe0a0",
"metadata": {},
"source": [
"A few things to highlight. First, the circuit has gotten longer with respect to the initial one. Second, although we had only 12 two-qubit gates in the original circuit, the fact that those were not in the basis set means that, when expanded, it requires more than a single cx gate to implement. All said, unrolling to the basis set of gates leads to an increase in the depth of a quantum circuit and the number of gates. The deeper a circuit, the more noise and the bigger the errors our quantum states will be subject to. Keep this in mind when designing circuits that are run on real quantum devices!"
]
},
{
"cell_type": "markdown",
"id": "d884df68",
"metadata": {},
"source": [
"In general, [transpilation](https://qiskit.org/documentation/apidoc/transpiler.html) is the process of rewriting a given input circuit to match the topology of a specific quantum device, and/or to optimize the circuit for execution on present day noisy quantum systems.\n",
"\n",
"Most circuits must undergo a series of transformations that make them compatible with a given target device, and optimize them to reduce the effects of noise on the resulting outcomes. Rewriting quantum circuits to match hardware constraints and optimizing for performance can be far from trivial. The flow of logic in the rewriting tool chain need not be linear, and can often have iterative sub-loops, conditional branches, and other complex behaviors. The details are not important for you at this stage, just be aware that transpilation is an important but potentially very complex process."
]
},
{
"cell_type": "markdown",
"id": "8819ab54",
"metadata": {},
"source": [
"### Problem 1: Noisy expectation values due to sampling"
]
},
{
"cell_type": "markdown",
"id": "caea2eb1",
"metadata": {},
"source": [
"As mentioned earlier, the measurement outcome of a digital quantum computation are always bitstrings. If we want to compute the expectation value of some observable, in practice, we need to measure reapeatedly, collect all occurences of each measured bit(string), and then compute weighted averages over those. Hence, the result will necessary be noisy. In this exercise you will study how the number of measurements (shots) affects the standard error on the estimated expectation value."
]
},
{
"cell_type": "markdown",
"id": "889cd1ec",
"metadata": {},
"source": [
"a) Define a function that given a circuit and a specified number of shots/measurements returns the expectation value of the average z-magnetization\n",
"\n",
"$$ \\hat{O} = \\frac 1 N \\sum_i \\hat{Z}_i$$\n",
"\n",
"You can only use the `Sampler()` and not the `Estimator()` for this exercise."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "70ed960e",
"metadata": {},
"outputs": [],
"source": [
"def measure_average_sz_magnetization(circ, n_shots):\n",
" circ.measure_all()\n",
" sampler = StatevectorSampler()\n",
" result = sampler.run([circ], shots=n_shots).result()[0]\n",
" counts = result.data.meas.get_counts()\n",
"\n",
" # TODO convert counts into an estimate of the average z-magnetization\n",
" return # TODO\n",
"\n",
"\n",
"# Test cases\n",
"test_circ = QuantumCircuit(2)\n",
"test_circ.x(0)\n",
"assert np.isclose(measure_average_sz_magnetization(test_circ, 1000), 0.0)\n",
"test_circ = QuantumCircuit(2)\n",
"test_circ.x(0)\n",
"test_circ.x(1)\n",
"assert np.isclose(measure_average_sz_magnetization(test_circ, 1000), -1.0)"
]
},
{
"cell_type": "markdown",
"id": "51339fc3",
"metadata": {},
"source": [
"💡Hint\n",
"\n",
"Given a observable that is diagonal in the $Z$-basis: $\\hat{O} = \\sum_\\mathbf{z} f(\\mathbf{z}) |\\mathbf{z}\\rangle\\langle \\mathbf{z}|$ with possible measurement outcomes $f(\\mathbf{z})$, the expectation value $\\langle \\psi |\\hat{O}|\\psi\\rangle$ can be estimated as\n",
"\n",
"$$ \\langle \\psi |\\hat{O}|\\psi\\rangle \\sim \\frac{1}{M_s} \\sum_{i=1}^{M_s} f(\\mathbf{z}^{(i)}) ,$$\n",
"\n",
"where $M_s$ is the number of measurements (shots) and $f(\\mathbf{z}^{(i)})$ is the measurement outcome for the ith measurement.\n",
"\n",
"As an example, if you measure the j-th qubit on the quantum device, $\\hat{Z}_j = \\sum_\\mathbf{z} (-1)^{z_j} |\\mathbf{z}\\rangle\\langle \\mathbf{z}|$, so $f(\\mathbf{z})=(-1)^{z_j}$.\n",
"\n",
"(Note that in Qiskit 0 corresponds to eigenvalue +1 (spin up) and 1 to eigenvalue -1 (spin down).)\n",
"\n",
"In case you are not so familiar with python, these two lines might be a good starting point:\n",
"```python\n",
"for (bitstring, count) in counts.items():\n",
" bits = np.array([int(bit) for bit in bitstring])\n",
" ...\n",
"```\n",
""
]
},
{
"cell_type": "markdown",
"id": "7029bd3a",
"metadata": {},
"source": [
"✨Solution\n",
"\n",
"```python\n",
"avg_mag = 0\n",
"for (bitstring, count) in counts.items():\n",
" bits = np.array([int(bit) for bit in bitstring])\n",
" szs = (-1)**(bits)\n",
" avg_mag += np.mean(szs) * count / n_shots\n",
"return avg_mag\n",
"```\n",
""
]
},
{
"cell_type": "markdown",
"id": "d80fe511",
"metadata": {},
"source": [
"b) Given the circuit below, compute the expectation value a few times (~20) and then calculate the standard deviation over all individual results you obtained. Set the number of shots to 1000."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1e2d3490",
"metadata": {},
"outputs": [],
"source": [
"circ = time_evol_circuit(N, 10.0, dt, J, gamma)\n",
"# TODO"
]
},
{
"cell_type": "markdown",
"id": "201c075d",
"metadata": {},
"source": [
"💡Hint\n",
"\n",
"You can use `np.std()` to calculate the standard deviation or do it by hand.\n",
""
]
},
{
"cell_type": "markdown",
"id": "3298d1c8",
"metadata": {},
"source": [
"c) Repeat the steps from above for different number of shots as specified below. Then plot the computed standard deviations against the number of shots."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e1d333bd",
"metadata": {},
"outputs": [],
"source": [
"n_shots_list = np.logspace(2, 5, 4, base=10, dtype=int)\n",
"stds = []\n",
"# TODO"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ddce04da",
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.loglog()\n",
"plt.plot(n_shots_list, stds)\n",
"plt.xlabel(\"number of shots\")\n",
"plt.ylabel(\"standard deviation\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "9aa84c88",
"metadata": {},
"source": [
"✨Solution\n",
"\n",
"```python\n",
"stds = []\n",
"for n_shots in n_shots_list:\n",
" result = []\n",
" for i in range(20):\n",
" result.append(measure_average_sz_magnetization(circ, n_shots))\n",
" stds.append(np.std(result))\n",
"```\n",
""
]
},
{
"cell_type": "markdown",
"id": "f03b58ca",
"metadata": {},
"source": [
"d) Can you figure out a simple mathematical relationship of how the standard deviation of the mean scales with the number of shots? What implications does this have on results obtained from a quantum computation?"
]
},
{
"cell_type": "markdown",
"id": "55269ae3",
"metadata": {},
"source": [
"💡Hint\n",
"\n",
"It might help to change the scales of the axes to a log-log one. You can simply use `plt.loglog()`. What does a straight line on a log-log plot correspond to? What is the approximate slope of the line? If you don't want to guess, you can use `np.polyfit(np.log(n_shots_list), np.log(stds), 1)` to fit a straight line to your curve.\n",
""
]
},
{
"cell_type": "markdown",
"id": "f4a184c7-96b9-40bb-8e43-e60ab26490ca",
"metadata": {},
"source": [
"### Problem 2: Digitial adiabatic state preperation "
]
},
{
"cell_type": "markdown",
"id": "45eca8ed-9848-42b0-9835-c081d9373ddf",
"metadata": {},
"source": [
"Let's turn to the problem of ground state preparation on quantum computers. One way of obtaining ground states is by __adiabatic time evolution__. Let's say that we know how to prepare the ground state of some Hamiltonian $H_i$ and our goal is to find the ground state of a different, more complicated Hamiltonian $H_f$. We then define a time-dependent Hamiltonian $H(t)$ that interpolates between the initial Hamiltonian $H_i$ and final Hamiltonian $H_f$, i.e.,\n",
"\n",
"$$H(t) = (1-u(t)) H_i + u(t) H_f ,$$\n",
"\n",
"where $u(t)$ is a smooth function of time which is 0 at $t=0$ and 1 at $t=T$. For example: $u(t) = t/T$.\n",
"\n",
"The [adiabatic theorem](https://en.wikipedia.org/wiki/Adiabatic_theorem) then guarantees that we stay in the instantaneous ground state of $H(t)$ throughout the time evolution as long as we evolve slowly enough and as long as there is a gap between the ground and first excited state.\n",
"\n",
"Let's apply the idea of adiabatic ground state preparation to our Ising model example. As an example imagine you want to prepare the ground state of $H_f = J \\sum_i X_i X_{i+1} - \\sum_i Z_i$ for some fixed coupling strength $J$. We know that the initial state on a quantum computer is $|0\\dots 0 \\rangle$ which is exactly the ground state of the transverse-field Hamiltonian, i.e., $H_i = - \\sum_i Z_i$. Hence, we need to adiabatically evolve with the following time-dependent Hamiltonian:\n",
"\n",
"$$H(t) = (1-u(t)) H_i + u(t) H_f = u(t) J \\sum_i X_i X_{i+1} - \\sum_i Z_i .$$\n",
"\n",
"So we adiabtically turn on the interaction term of the Ising model.\n",
"\n",
"We now want to simulate the time evolution of our initial $|0\\dots 0 \\rangle$ state with the Hamiltonian above on a quantum device. Let's figure out how to do this step by step. First, we compute the propagator $G_{\\Delta_t}(t)$ that evolves a quantum state at time $t$ by a small (infinitesimal) time increment $\\Delta_t$:\n",
"\n",
"$$ \\ket{\\Psi(t+\\Delta_t)} = G_{\\Delta_t}(t) \\ket{\\Psi(t)} = e^{-i \\Delta_t \\bar{H}_t} \\ket{\\Psi(t)}.$$\n",
"\n",
"a) Determine $\\bar{H}_t$ up to first order in the Magnus expansion for the time-dependent Ising Hamiltonian above using pen and paper (no need to write code). Assume a linear ramp, i.e., $u(t)= t/T$ and open boundary conditions."
]
},
{
"cell_type": "markdown",
"id": "6d87c50e",
"metadata": {},
"source": [
"💡Hint\n",
"Remember from class that the first term in the Magnus expansion is just the time-averaged Hamiltonian, i.e., \n",
"\n",
"$$\\bar{H}^{(1)}_t = \\frac{1}{\\Delta_t}\\int_{t}^{t+\\Delta_t} H(t') dt'$$\n",
""
]
},
{
"cell_type": "markdown",
"id": "e1d53a66",
"metadata": {},
"source": [
"b) Assuming that $\\Delta_t$ is small, use the first-order Trotter-Suzuki decomposition to split the propagator $G_{\\Delta_t}(t)$ into a product of local terms (again using pen and paper).\n",
"\n",
"How does the error of the obtained expression scale w.r.t. the exact time evolution?"
]
},
{
"cell_type": "markdown",
"id": "c2926ff7",
"metadata": {},
"source": [
"c) Now write down the propagator $G(t)$ for the time evolution from time $t=0$ to $t=n_{steps} \\Delta_t$ using the result above.\n",
"\n",
"$$ \\ket{\\Psi(t)} = G(t) \\ket{\\Psi(0)} = ... \\ket{\\Psi(0)}$$"
]
},
{
"cell_type": "markdown",
"id": "08f06c8f",
"metadata": {},
"source": [
"✨Solution\n",
"\n",
"$$G(t) = \\prod_{n=0}^{n_{steps}-1}\\left[\\prod_{i=1}^{N}\\exp(i\\Delta_t \\hat{Z}_i) \\prod_{i=1}^{N-1}\\exp(-i \\frac{\\Delta_t^2}{2T}(2n+1) \\hat{X}_i\\hat{X}_{i+1})\\right]$$\n",
""
]
},
{
"cell_type": "markdown",
"id": "59f42660",
"metadata": {},
"source": [
"d) Modify the function `time_evol_circuit` from above to the time-dependent case considered here.\n",
"\n",
"Note that as before $N$ is the total number of spins/qubits, $t$ is the time up to which we want to evolve, $(T,J)$ are parameters in the Hamiltonian, and `dt` is the step size $\\Delta_t$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4a36deed",
"metadata": {},
"outputs": [],
"source": [
"def adiabatic_time_evol_circuit(N, t, T, dt, J):\n",
" n_steps = round(t / dt) # number of Suzuki-Trotter steps\n",
" circ = QuantumCircuit(N)\n",
" for step in range(n_steps):\n",
" # TODO\n",
" return circ"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fadfe82f",
"metadata": {},
"outputs": [],
"source": [
"np.allclose(\n",
" -0.5282837412829096,\n",
" Statevector(adiabatic_time_evol_circuit(4, 1.0, 1.0, 0.1, 0.7))[0].real,\n",
") # should return True"
]
},
{
"cell_type": "markdown",
"id": "6026e1a8",
"metadata": {},
"source": [
"During the time evolution we want to measure the energy w.r.t. $H(t)$ to see if we actually stay in the instantaneous ground state at all times.\n",
"\n",
"e) Define a function that given a time $t$ returns the Hamiltonian $H(t)$ as a `SparsePauliOp`. There are different ways of doing this, any one is fine."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6b97e546",
"metadata": {},
"outputs": [],
"source": [
"def hamiltonian(N, t, T, J):\n",
" # TODO\n",
" return # TODO return the Hamiltonian as a SparsePauliOp"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15113d4b",
"metadata": {},
"outputs": [],
"source": [
"np.allclose(\n",
" np.array(\n",
" [\n",
" [-2.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.4 + 0.0j],\n",
" [0.0 + 0.0j, 0.0 + 0.0j, 0.4 + 0.0j, 0.0 + 0.0j],\n",
" [0.0 + 0.0j, 0.4 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],\n",
" [0.4 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 2.0 + 0.0j],\n",
" ]\n",
" ),\n",
" hamiltonian(2, 1.0, 2.0, 0.8).to_matrix(),\n",
") # should return True"
]
},
{
"cell_type": "markdown",
"id": "f13552b9",
"metadata": {},
"source": [
"💡Hint\n",
"\n",
"```python\n",
"def hamiltonian(N, t, T, J):\n",
" h_terms = [(\"Z\", [i], -1) for i in range(N)]\n",
" h_terms += # similar for the XX term\n",
" return SparsePauliOp.from_sparse_list(h_terms, num_qubits=N)\n",
"```\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "d7975976",
"metadata": {},
"source": [
"Now it's time to simulate the evolution and measure the energy of the Hamiltonian."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9ca539a5",
"metadata": {},
"outputs": [],
"source": [
"N = 5 # number of qubits\n",
"T = 1.0 # total evolution time\n",
"dt = 1e-2 # Trotter time step size (Delta t)\n",
"J = 0.5 # Ising coupling strength\n",
"n_points = 50 # number of time stamps at which we want to measure the expectation value of the observable\n",
"ts = np.linspace(0, T, n_points + 1) # measure at these times\n",
"measured_energies = [] # store measured energies here"
]
},
{
"cell_type": "markdown",
"id": "b53039c7",
"metadata": {},
"source": [
"f) Use the `StatevectorEstimator` to run the circuits and measure the corresponding Hamiltonian $H(t)$ at all times `ts`. Plot the expectation values as a function of time $t$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "51dee96e",
"metadata": {},
"outputs": [],
"source": [
"# TODO"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "99632779",
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.plot(ts, measured_energies)\n",
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(r\"$\\langle \\psi_t | H(t) | \\psi_t \\rangle $\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "f156a5fa",
"metadata": {},
"source": [
"g) To see whether we actually stay in the instantaneous ground state of $H(t)$ compute the true ground state at all times $t$ using exact diagonalization. Plot the results on top of the figure above.\n",
"\n",
"Note that you can convert a `SparsePauliOp` to a matrix format using the `to_matrix()` method."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1f4913ad",
"metadata": {},
"outputs": [],
"source": [
"exact_energies = [] # store the exact ground state energies here\n",
"\n",
"# TODO"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "762c8ccc",
"metadata": {},
"outputs": [],
"source": [
"np.allclose(-5.251939000256705, exact_energies[-1]) # should return True"
]
},
{
"cell_type": "markdown",
"id": "db38bdb1",
"metadata": {},
"source": [
"💡Hint\n",
"\n",
"Just import `from scipy.sparse.linalg import eigsh` and run `eigsh` like we did in the previous tutorial. \n",
""
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8de3db7c",
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.plot(ts, measured_energies)\n",
"plt.plot(ts, exact_energies, linestyle=\"--\")\n",
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(r\"$\\langle \\psi_t | H(t) | \\psi_t \\rangle $\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "f41934df",
"metadata": {},
"source": [
"h) Why do the two curves deviate at longer times? What parameters can you change to follow the instantaneous ground state all the way to the final time $t=T$ more accurately? Try it out!"
]
},
{
"cell_type": "markdown",
"id": "31d3b108",
"metadata": {},
"source": [
"💡Hint\n",
"\n",
"Recall the adiabatic theorem and it's conditions.\n",
""
]
},
{
"cell_type": "markdown",
"id": "a86d3f50",
"metadata": {},
"source": [
"i) 🚀Advanced: Critical slowing down\n",
"\n",
"One can show that the total required evolution time $T$ of an adiabatic protocol has to be larger than the inverse cubed energy gap between ground and first excited state at any point during the evolution:\n",
"\n",
"$$ T \\geq \\frac{C}{\\Delta^3}, \\qquad \\Delta = E_1 - E_0 ,$$\n",
"\n",
"where $C$ is some constant.\n",
"\n",
"Now imagine you want to prepare the ground state of the Ising model at the critical point $J/\\Gamma = 1$ for an infinitely-large spin chain using adiabatic evolution. How long would you have to evovle for? You can answer this question by pure reasoning, however, feel free to calculate (analytically or numerically) what happens to $T$ when you increase the number of spins $N$. \n",
"\n",
"Intuitively, why does adiabatic state preparation break down when it crosses a phase transition?"
]
},
{
"cell_type": "markdown",
"id": "0f089979",
"metadata": {},
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "env",
"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.13.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}