{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "453ab30f-a0c6-448b-8579-710968b0f288", "metadata": {}, "source": [ "---\n", "title: Instances and extensions\n", "description: This lesson shows the complete workflow for some variational approaches and describes how they can be modified.\n", "---\n", "\n", "{/* cspell: ignore disp SSVQE's forall */}\n", "\n", "# Instances and extensions\n", "\n", "This chapter will cover several quantum variational algorithms, including\n", "\n", "* [Variational Quantum Eigensolver (VQE)](https://arxiv.org/abs/1304.3061)\n", "* [Subspace Search VQE (SSVQE)](https://arxiv.org/abs/1810.09434)\n", "* [Variational Quantum Deflation (VQD)](https://arxiv.org/abs/1805.08138)\n", "* [Quantum Sampling Regression (QSR)](https://arxiv.org/pdf/2012.02338)\n", "\n", "By using these algorithms, we will learn about several design ideas that can be incorporated into custom variational algorithms, such as weights, penalties, over-sampling, and under-sampling. We encourage you to experiment with these concepts and share your findings with the community.\n", "\n", "The Qiskit patterns framework applies to all these algorithms - but we will explicitly call out the steps only in the first example.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "341eacfa-f6fe-4727-a9de-5babe59455f6", "metadata": {}, "source": [ "## Variational Quantum Eigensolver (VQE)\n", "\n", "[VQE](https://arxiv.org/abs/1304.3061) is one of the most widely used variational quantum algorithms, setting up a template for other algorithms to build upon.\n", "\n", "![A diagram showing how VQE uses the reference state and ansatz to estimate a cost function, and then iterate using variational parameters.](https://quantum.cloud.ibm.com/learning/images/courses/variational-algorithm-design/instances-and-extensions/instances-vqe.svg)\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "9453754c-03c5-44cf-accc-42ad434d0645", "metadata": {}, "source": [ "### Step 1: Map classical inputs to a quantum problem\n", "\n", "### Theoretical layout\n", "\n", "VQE's layout is simple:\n", "\n", "* Prepare reference operators $U_R$\n", " * We start from the state $|0\\rangle$ and go to the reference state $|\\rho\\rangle$\n", "* Apply the variational form $U_V(\\vec\\theta_{i,j})$ to create an ansatz $U_A(\\vec\\theta_{i,j})$\n", " * We go from the state $|\\rho\\rangle$ to $U_V(\\vec\\theta_{i,j})|\\rho\\rangle = |\\psi(\\vec\\theta_{i,j})\\rangle$\n", "* Bootstrap at $i=0$ if we have a similar problem (typically found via classical simulation or sampling)\n", " * Each optimizer will be bootstrapped differently, resulting in an initial set of parameter vectors $\\Theta_0 := \\\\{ {\\vec\\theta_{0,j} | j \\in \\mathcal{J}_\\text{opt}^0} \\\\}$ (for example, from an initial point $\\vec\\theta_0$).\n", "* Evaluate the cost function $C(\\vec\\theta_{i,j}) := \\langle \\psi(\\vec{\\theta}) | \\hat{H} | \\psi(\\vec{\\theta})\\rangle$ for all prepared states on a quantum computer.\n", "* Use a classical optimizer to select the next set of parameters $\\Theta_{i+1}$.\n", "* Repeat the process until convergence is reached.\n", "\n", "This is a simple classical optimization loop where we evaluate the cost function. Some optimizers may require multiple evaluations to calculate a gradient, determine the next iteration, or assess convergence.\n", "\n", "Here's the example for the following observable:\n", "\n", "$$\n", "\\hat{O}_1 = 2 II - 2 XX + 3 YY - 3 ZZ,\n", "$$\n", "\n" ] }, { "cell_type": "markdown", "id": "ee353736-d996-4c96-bac5-15b144da9a74", "metadata": {}, "source": [ "### Implementation\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "db20526e-48d4-44eb-9d2b-6c5d2ce5777e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qiskit import QuantumCircuit\n", "from qiskit.quantum_info import SparsePauliOp\n", "from qiskit.circuit.library import TwoLocal\n", "import numpy as np\n", "\n", "theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()\n", "observable = SparsePauliOp.from_list([(\"II\", 2), (\"XX\", -2), (\"YY\", 3), (\"ZZ\", -3)])\n", "\n", "reference_circuit = QuantumCircuit(2)\n", "reference_circuit.x(0)\n", "\n", "variational_form = TwoLocal(\n", " 2,\n", " rotation_blocks=[\"rz\", \"ry\"],\n", " entanglement_blocks=\"cx\",\n", " entanglement=\"linear\",\n", " reps=1,\n", ")\n", "\n", "ansatz = reference_circuit.compose(variational_form)\n", "\n", "ansatz.decompose().draw(\"mpl\")" ] }, { "cell_type": "code", "execution_count": 2, "id": "88797da9-1519-4ebb-8e5c-447fa9ccf662", "metadata": {}, "outputs": [], "source": [ "def cost_func_vqe(parameters, ansatz, hamiltonian, estimator):\n", " \"\"\"Return estimate of energy from estimator\n", "\n", " Parameters:\n", " params (ndarray): Array of ansatz parameters\n", " ansatz (QuantumCircuit): Parameterized ansatz circuit\n", " hamiltonian (SparsePauliOp): Operator representation of Hamiltonian\n", " estimator (Estimator): Estimator primitive instance\n", "\n", " Returns:\n", " float: Energy estimate\n", " \"\"\"\n", "\n", " estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])\n", " estimator_result = estimator_job.result()[0]\n", "\n", " cost = estimator_result.data.evs[0]\n", " return cost" ] }, { "cell_type": "code", "execution_count": 3, "id": "08859a29-ee3c-416b-92fd-f60176fe6a80", "metadata": {}, "outputs": [], "source": [ "from qiskit.primitives import StatevectorEstimator\n", "\n", "estimator = StatevectorEstimator()" ] }, { "attachments": {}, "cell_type": "markdown", "id": "ec815a0e-7f38-4152-b979-b92953d204f0", "metadata": {}, "source": [ "We can use this cost function to calculate optimal parameters\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "bfbffa18-b281-4b76-a57e-700c55ee2daa", "metadata": {}, "outputs": [ { "data": { "text/plain": [ " message: Optimization terminated successfully.\n", " success: True\n", " status: 1\n", " fun: -5.999999982445723\n", " x: [ 1.741e+00 9.606e-01 1.571e+00 2.115e-05 1.899e+00\n", " 1.243e+00 6.063e-01 6.063e-01]\n", " nfev: 136\n", " maxcv: 0.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# SciPy minimizer routine\n", "from scipy.optimize import minimize\n", "\n", "x0 = np.ones(8)\n", "\n", "result = minimize(\n", " cost_func_vqe, x0, args=(ansatz, observable, estimator), method=\"COBYLA\"\n", ")\n", "\n", "result" ] }, { "cell_type": "markdown", "id": "009eb817-583c-4240-9951-34f76adc92ea", "metadata": {}, "source": [ "### Step 2: Optimize problem for quantum execution\n", "\n", "We will select the least-busy backend, and import the necessary components from `qiskit_ibm_runtime`.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c06c9393-b4f6-4e4c-92b1-4ec062852e2c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "from qiskit_ibm_runtime import SamplerV2 as Sampler\n", "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n", "from qiskit_ibm_runtime import Session, EstimatorOptions\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "\n", "\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(operational=True, simulator=False)\n", "print(backend)" ] }, { "cell_type": "markdown", "id": "94d86bc0-b0fe-466c-a6e2-582f240ae014", "metadata": {}, "source": [ "We will transpile the circuit using the preset pass manager with optimization level 3, and we will apply the corresponding layout to the observable.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "09c82d1d-9a14-4601-baef-8080a2758201", "metadata": {}, "outputs": [], "source": [ "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", "\n", "pm = generate_preset_pass_manager(backend=backend, optimization_level=3)\n", "isa_ansatz = pm.run(ansatz)\n", "isa_observable = observable.apply_layout(layout=isa_ansatz.layout)" ] }, { "cell_type": "markdown", "id": "07772b29-e880-476c-a4bc-07f7591ad969", "metadata": {}, "source": [ "### Step 3: Execute using Qiskit Runtime primitives\n", "\n", "We are now ready to run our calculation on IBM Quantum® hardware. Because the cost function minimization is highly iterative, we will start a Runtime session. This way, we will only have to wait in a queue once. Once the job begins running, each iteration with updates to parameters will run immediately.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c785cd0b-2a26-4414-a320-9746d016083e", "metadata": {}, "outputs": [], "source": [ "x0 = np.ones(8)\n", "\n", "estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)\n", "\n", "with Session(backend=backend) as session:\n", " estimator = Estimator(mode=session, options=estimator_options)\n", "\n", " result = minimize(\n", " cost_func_vqe,\n", " x0,\n", " args=(isa_ansatz, isa_observable, estimator),\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 200, \"disp\": True},\n", " )\n", "session.close()\n", "print(result)" ] }, { "cell_type": "markdown", "id": "34ad36bb-e83a-4cfc-be9b-f5ecde1122ab", "metadata": {}, "source": [ "### Step 4: Post-process, return result in classical format\n", "\n", "We can see that the minimization routine successfully terminated, meaning we reached the default tolerance of the COBYLA classical optimizer. If we require a more precise result, we can specify a smaller tolerance. This may indeed be the case, since the result was several percent off compared to the result obtained by the simulator above.\n", "\n", "The value of x obtained is the current best guess for the parameters that minimize the cost function. If iterating to obtain a higher precision, those values should be used in place of the x0 initially used (a vector of ones).\n", "\n", "Finally, we note that the function was evaluated 96 times in the process of optimization. That might be different from the number of optimization steps, since some optimizers require multiple function evaluations in a single step, such as when estimating a gradient.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "015a8beb-ec4e-4e02-a937-6588c6a81b64", "metadata": {}, "source": [ "## Subspace Search VQE (SSVQE)\n", "\n", "[SSVQE](https://arxiv.org/abs/1810.09434) is a variant of VQE that allows obtaining the first $k$ eigenvalues of an observable $\\hat{H}$ with eigenvalues $\\{\\lambda_0, \\lambda_1,...,\\lambda_{N-1}\\}$, where $N\\geq k$. Without loss of generality, we assume that $\\lambda_0<\\lambda_1<...<\\lambda_{N-1}$. SSVQE introduces a new idea by adding weights to help prioritize optimizing for the term with the largest weight.\n", "\n", "![A diagram showing how subspace-search VQE uses the components of variational algorithm.](https://quantum.cloud.ibm.com/learning/images/courses/variational-algorithm-design/instances-and-extensions/instances-ssvqe.svg)\n", "\n", "To implement this algorithm, we need $k$ mutually orthogonal reference states $\\{ |\\rho_j\\rangle \\}_{j=0}^{k-1}$, meaning $\\langle \\rho_j | \\rho_l \\rangle = \\delta_{jl}$ for $j,lw_l$, and $U_V(\\vec{\\theta})$ is the user-defined variational form.\n", "\n", "The SSVQE algorithm relies on the fact that eigenstates corresponding to different eigenvalues are mutually orthogonal. Specifically, the inner product of $U_V(\\vec{\\theta})|\\rho_j\\rangle$ and $U_V(\\vec{\\theta})|\\rho_l\\rangle$ can be expressed as:\n", "\n", "$$\n", "\\begin{aligned}\n", "\\langle \\rho_j | U_{V}^{\\dagger}(\\vec{\\theta})U_{V}(\\vec{\\theta})|\\rho_l \\rangle\n", "\n", "& = \\langle \\rho_j | I |\\rho_l \\rangle \\\\[1mm]\n", "\n", "& = \\langle \\rho_j | \\rho_l \\rangle \\\\[1mm]\n", "\n", "& = \\delta_{jl}\n", "\n", "\\end{aligned}\n", "$$\n", "\n", "The first equality holds because $U_{V}(\\vec{\\theta})$ is a quantum operator and is therefore unitary. The last equality holds because of the orthogonality of the reference states $|\\rho_j\\rangle$. The fact that orthogonality is preserved through unitary transformations is deeply related to the principle of conservation of information, as expressed in quantum information science. Under this view, non-unitary transformations represent processes where information is either lost or injected.\n", "\n", "Weights $w_j$ help ensure that all the states are eigenstates. If the weights are sufficiently different, the term with the largest weight (that is, $w_0$) will be given priority during optimization over the others. As a result, the resulting state $U_{V}(\\vec{\\theta})|\\rho_0 \\rangle$ will become the eigenstate corresponding to $\\lambda_0$. Because $\\{ U_{V}(\\vec{\\theta})|\\rho_j\\rangle \\}_{j=0}^{k-1}$ are mutually orthogonal, the remaining states will be orthogonal to it and, therefore, contained in the subspace corresponding to the eigenvalues $\\{\\lambda_1,...,\\lambda_{N-1}\\}$.\n", "\n", "Applying the same argument to the rest of the terms, the next priority would then be the term with weight $w_1$, so $U_{V}(\\vec{\\theta})|\\rho_1 \\rangle$ would be the eigenstate corresponding to $\\lambda_1$, and the other terms would be contained in the eigenspace of $\\{\\lambda_2,...,\\lambda_{N-1}\\}$.\n", "\n", "By reasoning inductively, we deduce that $U_{V}(\\vec{\\theta})|\\rho_j \\rangle$ will be an approximate eigenstate of $\\lambda_j$ for $0\\leq j < k.$\n", "\n", "### Theoretical layout\n", "\n", "SSVQE's can be summarized as follows:\n", "\n", "* Prepare several reference states by applying a unitary U\\_R to k different computational basis states\n", " * This algorithm requires the usage of $k$ mutually orthogonal reference states $\\{ |\\rho_j\\rangle \\}_{j=0}^{k-1}$, such that $\\langle \\rho_j | \\rho_l \\rangle = \\delta_{jl}$ for $j,l \\lambda_0$, instead of minimizing the cost function $C_0(\\vec{\\theta}) := \\langle \\psi(\\vec{\\theta}) | \\hat{H} | \\psi(\\vec{\\theta})\\rangle$, we optimize:\n", "\n", "$$\n", "C_1(\\vec{\\theta}) :=\n", "C_0(\\vec{\\theta})+ \\beta_0 |\\langle \\psi(\\vec{\\theta})| \\psi(\\vec{\\theta^0})\\rangle |^2\n", "$$\n", "\n", "The positive value $\\beta_0$ should ideally be greater than $\\lambda_1-\\lambda_0$.\n", "\n", "This introduces a new cost function that can be viewed as a constrained problem, where we minimize $C_\\text{VQE}(\\vec{\\theta}) = \\langle \\psi(\\vec{\\theta}) | \\hat{H} | \\psi(\\vec{\\theta})\\rangle$ subject to the constraint that the state must be orthogonal to the previously obtained $|\\psi(\\vec{\\theta^0})\\rangle$, with $\\beta_0$ acting as a penalty term if the constraint is not satisfied.\n", "\n", "Alternatively, this new problem can be interpreted as running VQE on the new observable:\n", "\n", "$$\n", "\\hat{H_1} := \\hat{H} + \\beta_0 |\\psi(\\vec{\\theta^0})\\rangle \\langle \\psi(\\vec{\\theta^0})|\n", "\\quad \\Rightarrow \\quad\n", "C_1(\\vec{\\theta}) = \\langle \\psi(\\vec{\\theta}) | \\hat{H_1} | \\psi(\\vec{\\theta})\\rangle,\n", "$$\n", "\n", "Assuming that the solution to the new problem is $|\\psi(\\vec{\\theta^1})\\rangle$, the expected value of $\\hat{H}$ (not $\\hat{H_1}$) should be $ \\langle \\psi(\\vec{\\theta^1}) | \\hat{H} | \\psi(\\vec{\\theta^1})\\rangle = \\lambda_1$.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "0b713dad-29ff-4693-9e72-df9b6de3afd0", "metadata": {}, "source": [ "To obtain the third eigenvalue $\\lambda_2$, the cost function to optimize is:\n", "\n", "$$\n", "C_2(\\vec{\\theta}) :=\n", "C_1(\\vec{\\theta}) + \\beta_1 |\\langle \\psi(\\vec{\\theta})| \\psi(\\vec{\\theta^1})\\rangle |^2\n", "$$\n", "\n", "where $\\beta_1$ is a positive constant large enough to enforce orthogonality of the solution state to both $|\\psi(\\vec{\\theta^0})\\rangle$ and $|\\psi(\\vec{\\theta^1})\\rangle$. This penalizes states in the search space that do not meet this requirement, effectively restricting the search space. Thus, the optimal solution of the new problem should be the eigenstate corresponding to $\\lambda_2$.\n", "\n", "Like the previous case, this new problem can also be interpreted as VQE with the observable:\n", "\n", "$$\n", "\\hat{H_2} := \\hat{H_1} + \\beta_1 |\\psi(\\vec{\\theta^1})\\rangle \\langle \\psi(\\vec{\\theta^1})|\n", "\\quad \\Rightarrow \\quad\n", "C_2(\\vec{\\theta}) = \\langle \\psi(\\vec{\\theta}) | \\hat{H_2} | \\psi(\\vec{\\theta})\\rangle.\n", "$$\n", "\n", "If the solution to this new problem is $|\\psi(\\vec{\\theta^2})\\rangle$, the expected value of $\\hat{H}$ (not $\\hat{H_2}$) should be $ \\langle \\psi(\\vec{\\theta^2}) | \\hat{H} | \\psi(\\vec{\\theta^2})\\rangle = \\lambda_2$.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "a17f468d-37bf-4d7a-8465-ab6cda804d0a", "metadata": {}, "source": [ "Analogously, to obtain the $k$-th eigenvalue $\\lambda_{k-1}$, you would minimize the cost function:\n", "\n", "$$\n", "C_{k-1}(\\vec{\\theta}) :=\n", "C_{k-2}(\\vec{\\theta}) + \\beta_{k-2} |\\langle \\psi(\\vec{\\theta})| \\psi(\\vec{\\theta^{k-2}})\\rangle |^2,\n", "$$\n", "\n", "Remember that we defined $\\vec{\\theta^j}$ such that $\\langle \\psi(\\vec{\\theta^j}) | \\hat{H} | \\psi(\\vec{\\theta^j})\\rangle = \\lambda_j, \\forall j" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qiskit.circuit.library import TwoLocal\n", "\n", "ansatz = TwoLocal(2, rotation_blocks=[\"ry\", \"rz\"], entanglement_blocks=\"cz\", reps=1)\n", "\n", "ansatz.decompose().draw(\"mpl\")" ] }, { "cell_type": "markdown", "id": "20934503-046f-4daa-b99c-43c99060b46d", "metadata": {}, "source": [ "First, we'll set up a function that calculates the state fidelity -- a percentage of overlap between two states that we'll use as a penalty for VQD:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "188d306a-48f2-4c36-9d90-2a243269f847", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "\n", "def calculate_overlaps(ansatz, prev_circuits, parameters, sampler):\n", " def create_fidelity_circuit(circuit_1, circuit_2):\n", " \"\"\"\n", " Constructs the list of fidelity circuits to be evaluated.\n", " These circuits represent the state overlap between pairs of input circuits,\n", " and their construction depends on the fidelity method implementations.\n", " \"\"\"\n", "\n", " if len(circuit_1.clbits) > 0:\n", " circuit_1.remove_final_measurements()\n", " if len(circuit_2.clbits) > 0:\n", " circuit_2.remove_final_measurements()\n", "\n", " circuit = circuit_1.compose(circuit_2.inverse())\n", " circuit.measure_all()\n", " return circuit\n", "\n", " overlaps = []\n", "\n", " for prev_circuit in prev_circuits:\n", " fidelity_circuit = create_fidelity_circuit(ansatz, prev_circuit)\n", " sampler_job = sampler.run([(fidelity_circuit, parameters)])\n", " meas_data = sampler_job.result()[0].data.meas\n", "\n", " counts_0 = meas_data.get_int_counts().get(0, 0)\n", " shots = meas_data.num_shots\n", " overlap = counts_0 / shots\n", " overlaps.append(overlap)\n", "\n", " return np.array(overlaps)" ] }, { "cell_type": "markdown", "id": "36345bbd-ae25-41d3-b3f4-04b93d822018", "metadata": {}, "source": [ "It's time to write VQD's cost function. As before when we calculated only the ground state, we will determine the lowest energy state using the Estimator primitive. However, as described above, we will now add a penalty term to ensure orthogonality of higher-energy states. That is, for each new excited state, a penalty is added for any overlap between the current variational state and the lower-energy eigenstates already found.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5367c47a-5e76-45ca-b9c6-fc711db2e64a", "metadata": {}, "outputs": [], "source": [ "def cost_func_vqd(\n", " parameters, ansatz, prev_states, step, betas, estimator, sampler, hamiltonian\n", "):\n", " estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])\n", "\n", " total_cost = 0\n", "\n", " if step > 1:\n", " overlaps = calculate_overlaps(ansatz, prev_states, parameters, sampler)\n", " total_cost = np.sum(\n", " [np.real(betas[state] * overlap) for state, overlap in enumerate(overlaps)]\n", " )\n", "\n", " estimator_result = estimator_job.result()[0]\n", "\n", " value = estimator_result.data.evs[0] + total_cost\n", "\n", " return value" ] }, { "cell_type": "markdown", "id": "6be376ac-cfef-4151-841c-d93fcbeb6c76", "metadata": {}, "source": [ "Note especially that the cost function above refers to the `calculate_overlaps` function, which actually creates a new quantum circuit. If we want to run on real hardware, that new circuit must also be transpiled, hopefully in an optimal way, to run on the backend we select. Note that transpilation has been built in to the `calculate_overlaps` or `cost_func_vqd` functions. Feel free to try modifying the code yourself to build in this additional (and conditional) transpilation - but this will also be done for you in the next lesson.\n", "\n", "In this lesson, we will run the VQD algorithm using the Statevector Sampler and Statevector Estimator:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "825e6e01-b0b0-4cf5-a05d-eeac33349fdd", "metadata": {}, "outputs": [], "source": [ "from qiskit.primitives import StatevectorEstimator as Estimator\n", "\n", "sampler = Sampler()\n", "estimator = Estimator()" ] }, { "cell_type": "markdown", "id": "08a92ff7-b0c0-467d-a322-cdcc7afb1add", "metadata": {}, "source": [ "We will introduce an observable to be estimated. In the next lesson we will add some physical context to this, like the excited state of a molecule. It may be helpful to think of this observable as the Hamiltonian of a system that can have excited states, even though this observable has not been chosen to match any particular molecule or atom.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "aa8ee41c-a0a7-4445-90b0-8554a31c3ce6", "metadata": {}, "outputs": [], "source": [ "from qiskit.quantum_info import SparsePauliOp\n", "\n", "observable = SparsePauliOp.from_list([(\"II\", 2), (\"XX\", -2), (\"YY\", 3), (\"ZZ\", -3)])" ] }, { "cell_type": "markdown", "id": "dabaa1a5-11df-46fc-92e9-0e1272b06a20", "metadata": {}, "source": [ "Here, we set the total number of states we wish to calculate (ground state and excited states, k), and the penalties (betas) for overlap between statevectors that should be orthogonal. The consequences of choosing betas to be too high or too low will be explored a bit in the next lesson. For now, we will simply use those provided below. We will start by using all zeros as our parameters. In your own calculations, you may want to use more clever starting parameters based on your knowledge of the system or on previous calculations.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "798fc5cf-c961-405d-a224-cf8189589190", "metadata": {}, "outputs": [], "source": [ "k = 3\n", "betas = [33, 33, 33]\n", "x0 = np.zeros(8)" ] }, { "cell_type": "markdown", "id": "f180bdeb-4fb8-4dc9-8f88-ddc60874c1db", "metadata": {}, "source": [ "We can now run the calculation:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dabcc978-0f86-41ab-92c5-9768645ce923", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " message: Optimization terminated successfully.\n", " success: True\n", " status: 1\n", " fun: -5.999999979545955\n", " x: [-5.150e-01 -5.452e-02 -1.571e+00 -2.853e-05 2.671e-01\n", " -2.672e-01 -8.509e-01 -8.510e-01]\n", " nfev: 131\n", " maxcv: 0.0\n", " message: Optimization terminated successfully.\n", " success: True\n", " status: 1\n", " fun: 4.024550284767612\n", " x: [-3.745e-01 1.041e+00 8.637e-01 1.202e+00 -8.847e-02\n", " 1.181e-02 7.611e-01 -3.006e-01]\n", " nfev: 110\n", " maxcv: 0.0\n", " message: Optimization terminated successfully.\n", " success: True\n", " status: 1\n", " fun: 5.608925562838559\n", " x: [-2.670e-01 1.280e+00 1.070e+00 -8.031e-01 -1.524e-01\n", " -6.956e-02 7.018e-01 1.514e+00]\n", " nfev: 90\n", " maxcv: 0.0\n" ] } ], "source": [ "from scipy.optimize import minimize\n", "\n", "prev_states = []\n", "prev_opt_parameters = []\n", "eigenvalues = []\n", "\n", "for step in range(1, k + 1):\n", " if step > 1:\n", " prev_states.append(ansatz.assign_parameters(prev_opt_parameters))\n", "\n", " result = minimize(\n", " cost_func_vqd,\n", " x0,\n", " args=(ansatz, prev_states, step, betas, estimator, sampler, observable),\n", " method=\"COBYLA\",\n", " options={\n", " \"maxiter\": 200,\n", " },\n", " )\n", " print(result)\n", "\n", " prev_opt_parameters = result.x\n", " eigenvalues.append(result.fun)" ] }, { "cell_type": "markdown", "id": "8d605a38-9711-4295-8245-201d3fb55e67", "metadata": {}, "source": [ "The values we obtained from the cost function are approximately -6.00, 4.02, and 5.61. The important thing about these results is that the function values are increasing. If we had obtained a first excited state that is lower in energy than our initial, unconstrained calculation of the ground state, that would have indicated an error somewhere in our code.\n", "\n", "The values of x are the parameters that yielded a statevector corresponding to each of these costs (energies).\n", "\n", "Finally, we note that all three minimizations converged to within the default tolerance of the classical optimizer (here COBYLA). They required 131, 110, and 90 function evaluations, respectively.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "c4b4d9b8-d730-41cb-b979-0c6e02faa442", "metadata": {}, "source": [ "## Quantum Sampling Regression (QSR)\n", "\n", "One of the main issues with VQE is the multiple calls to a quantum computer that are required to obtain the parameters for each step, for example, $k$, $k-1$, and so on. This is especially problematic when access to quantum devices is queued. While a [`Session`](/docs/guides/run-jobs-session) can be used to group multiple iterative calls, an alternative approach is to use sampling. By utilizing more classical resources, we can complete the full optimization process in a single call. This is where [Quantum Sampling Regression](https://arxiv.org/pdf/2012.02338) comes into play. Since access to quantum computers is still a low-offer/high-demand commodity, we find this trade-off to be both possible and convenient for many current studies. This approach harnesses all available classical capabilities while still capturing many of the inner workings and intrinsic properties of quantum computations that do not appear in simulation.\n", "\n", "![A diagram showing how QSR uses the components of a variational algorithm.](https://quantum.cloud.ibm.com/learning/images/courses/variational-algorithm-design/instances-and-extensions/instances-qsr.svg)\n", "\n", "The idea behind QSR is that the cost function $C(\\theta) := \\langle \\psi(\\theta) | \\hat{H} | \\psi(\\theta)\\rangle$ can be expressed as a Fourier series in the following manner:\n", "\n", "$$\n", "\\begin{aligned}\n", "C(\\vec{\\theta})\n", "\n", "& := \\langle \\psi(\\theta) | \\hat{H} | \\psi(\\theta)\\rangle \\\\[1mm]\n", "\n", "& := a_0 + \\sum_{k=1}^S[a_k\\cos(k\\theta)+ b_k\\sin(k\\theta)] \\\\[1mm]\n", "\n", "\\end{aligned}\n", "$$\n", "\n", "Depending on the periodicity and bandwidth of the original function, the set $S$ may be finite or infinite. For the purposes of this discussion, we will assume it is infinite. The next step is to sample the cost function $C(\\theta)$ multiple times in order to obtain the Fourier coefficients $\\{a_0, a_k, b_k\\}_{k=1}^S$. Specifically, since we have $2S+1$ unknowns, we will need to sample the cost function $2S+1$ times.\n", "\n", "If we then sample the cost function for $2S+1$ parameter values $\\{\\theta_1,...,\\theta_{2S+1}\\}$, we can obtain the following system:\n", "\n", "$$\n", "\\begin{pmatrix} 1 & \\cos(\\theta_1) & \\sin(\\theta_1) & \\cos(2\\theta_1) & ... & \\sin(S\\theta_1) \\\\\n", "1 & \\cos(\\theta_2) & \\sin(\\theta_2) & \\cos(2\\theta_2) & \\cdots & \\sin(S\\theta_2)\\\\\n", "\\vdots & \\vdots & \\vdots & \\vdots & \\ddots & \\vdots\\\\\n", "1 & \\cos(\\theta_{2S+1}) & \\sin(\\theta_{2S+1}) & \\cos(2\\theta_{2S+1}) & \\cdots & \\sin(S\\theta_{2S+1})\n", "\\end{pmatrix} \\begin{pmatrix} a_0 \\\\ a_1 \\\\ b_1 \\\\ a_2 \\\\ \\vdots \\\\ b_S \\end{pmatrix} = \\begin{pmatrix} C(\\theta_1) \\\\ C(\\theta_2) \\\\ \\vdots \\\\ C(\\theta_{2S+1}) \\end{pmatrix},\n", "$$\n", "\n", "that we'll rewrite as\n", "\n", "$$\n", "Fa=c.\n", "$$\n", "\n", "In practice, this system is generally not consistent because the cost function values $c$ are not exact. Therefore, it is usually a good idea to normalize them by multiplying them by $F^\\dagger$ on the left, which results in:\n", "\n", "$$\n", "F^\\dagger Fa = F^\\dagger c.\n", "$$\n", "\n", "This new system is always consistent, and its solution is a least-squares solution to the original problem. If we have $k$ parameters instead of just one, and each parameter $\\theta^i$ has its own $S_i$ for $i \\in {1,...,k}$, then the total number of samples required is:\n", "\n", "$$\n", "T=\\prod_{i=1}^k(2S_i+1)\\leq \\prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n,\n", "$$\n", "\n", "where $S_{\\max} = \\max_i(S_i)$. Furthermore, adjusting $S_{\\max}$ as a tunable parameter (instead of inferring it) opens up new possibilities, such as:\n", "\n", "* **Over-sampling** to improve accuracy.\n", "* **Under-sampling** to boost performance by reducing runtime overhead or eliminating local minima.\n", "\n" ] }, { "cell_type": "markdown", "id": "cd5399d1-8166-42e5-938c-fb2e2d600cb3", "metadata": {}, "source": [ "### Theoretical layout\n", "\n", "QSR's layout can be summarized as follows:\n", "\n", "* Prepare reference operators $U_R$.\n", " * We'll go from the state $|0\\rangle$ to the reference state $|\\rho\\rangle$\n", "* Apply the variational form $U_V(\\vec\\theta_{i,j})$ to create an ansatz $U_A(\\vec\\theta_{i,j})$.\n", " * Determine the bandwidth associated with each parameter in the ansatz. An upper bound is sufficient.\n", "* Bootstrap at $i=0$ if we have a similar problem (typically found by classical simulation or sampling).\n", "* Sample the cost function $C(\\vec\\theta) := a_0 + \\sum_{k=1}^S[a_k\\cos(k\\theta)+ b_k\\sin(k\\theta)]$ at least $T$ times.\n", " * $T=\\prod_{i=1}^k(2S_i+1)\\leq \\prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$\n", " * Decide whether to over-sample/under-sample to balance speed vs accuracy by adjusting $T$.\n", "* Compute the Fourier coefficients from the samples (that is, solve the normalized linear system of equations).\n", "* Solve for the global minimum of the resulting regression function on a classical machine.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "02bd8e30-adf8-4261-8b35-77f6709dce0d", "metadata": {}, "source": [ "## Summary\n", "\n", "With this lesson, you learned about multiple variational instances available:\n", "\n", "* General layout\n", "* Introducing weights and penalties to adjust a cost function\n", "* Exploring under-sampling vs over-sampling to trade-off speed vs accuracy\n", "\n", "These ideas can be adapted to form a custom variational algorithm that fits your problem. We encourage you to share your results with the community. The next lesson will explore how to use a variational algorithm to solve an application.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" } }, "nbformat": 4, "nbformat_minor": 2 }