update

Author: mhjensen

doc/Programs/VQEcodes/.ipynb_checkpoints/vqecode44-checkpoint.ipynb

@@ -0,0 +1,295 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "ceab64a3",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)\n",
+    "doconce format html vqecode44.do.txt  -->\n",
+    "<!-- dom:TITLE: VQE Simulation for a Two-Qubit Hamiltonian -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5d72e03a",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "# VQE Simulation for a Two-Qubit Hamiltonian\n",
+    "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo, Norway\n",
+    "\n",
+    "Date: **January 21, 2026**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "deebcb04",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "## Simple Hamiltonian\n",
+    "We consider the given 4×4 two-qubit Hamiltonian"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a942189d",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "$$\n",
+    "H = \\begin{pmatrix}\n",
+    "e_1 + V_z & 0      & 0      & V_x \\\\\n",
+    "0         & e_2 - V_z & V_x    & 0   \\\\\n",
+    "0         & H_x    & e_3 - V_z & 0   \\\\\n",
+    "H_x       & 0      & 0      & e_4 + V_z\n",
+    "\\end{pmatrix}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7baed6f3",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "Any two-qubit Hamiltonian can be expressed as a linear combination of tensor products of Pauli matrices (including the identity) .  In general we write\n",
+    "$$H = \\sum_{i,j\\in{I,X,Y,Z}} c_{ij},\\sigma_i\\otimes\\sigma_j,$$\n",
+    "where the coefficients \\(c_{ij} = \\tfrac14\\Tr[(\\sigma_i\\otimes\\sigma_j)H]\\) are real for a Hermitian H.  We implement this decomposition in code:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "02c45757",
+   "metadata": {
+    "collapsed": false,
+    "editable": true
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "# Define 2x2 Pauli matrices\n",
+    "I = np.array([[1,0],[0,1]])\n",
+    "X = np.array([[0,1],[1,0]])\n",
+    "Y = np.array([[0,-1j],[1j,0]])\n",
+    "Z = np.array([[1,0],[0,-1]])\n",
+    "\n",
+    "# Example values for constants (can be changed)\n",
+    "e1, e2, e3, e4 = 0.5, 1.0, 1.5, 2.0\n",
+    "Vx, Vz, Hx = 0.2, 0.3, 0.4\n",
+    "\n",
+    "# Construct the Hamiltonian matrix H\n",
+    "H = np.zeros((4,4), dtype=complex)\n",
+    "H[0,0] = e1 + Vz\n",
+    "H[1,1] = e2 - Vz\n",
+    "H[2,2] = e3 - Vz\n",
+    "H[3,3] = e4 + Vz\n",
+    "H[0,3] = Vx\n",
+    "H[3,0] = Hx\n",
+    "H[1,2] = Vx\n",
+    "H[2,1] = Hx\n",
+    "\n",
+    "# Decompose H into Pauli basis coefficients c_{ij}\n",
+    "paulis = {'I': I, 'X': X, 'Y': Y, 'Z': Z}\n",
+    "coeffs = {}\n",
+    "for i, Ai in paulis.items():\n",
+    "    for j, Bj in paulis.items():\n",
+    "        P = np.kron(Ai, Bj)\n",
+    "        c = np.trace(P.dot(H)) / 4.0\n",
+    "        if abs(c) > 1e-8:  # keep non-zero terms\n",
+    "            coeffs[(i,j)] = np.real_if_close(c)\n",
+    "\n",
+    "# Display non-zero Pauli terms of H\n",
+    "for (i,j), c in coeffs.items():\n",
+    "    print(f\"{c:+.4f} * {i}⊗{j}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eea60454",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "The code prints the nonzero terms of the Hamiltonian in the Pauli basis.  For example, with the above values one finds. yielding"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "697434ba",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "$$\n",
+    "H = 2.5\\,I\\otimes I +0.5\\,X\\otimes X-0.5\\,I\\otimes Z-1.0\\,Z\\otimes I+0.3\\,Z\\otimes Z.  T\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0681554e",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "This matches the matrix and confirms the decomposition in Pauli products."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "722b47ea",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "## Ansatz and Expectation Values\n",
+    "\n",
+    "We use a simple hardware-efficient ansatz: apply parameterized single-qubit rotations (here R_y(\\theta) and R_z(\\phi) gates) on each qubit, followed by an entangling CNOT gate .  Concretely, starting from $|00\\rangle$, we apply\n",
+    "$U(\\theta_0,\\theta_1,\\phi_0,\\phi_1) = \\bigl(R_z(\\phi_0)R_y(\\theta_0)\\otimes R_z(\\phi_1)R_y(\\theta_1)\\bigr)\\,\\text{CNOT}_{0\\to1}$.\n",
+    "This prepares a trial state $|\\psi(\\boldsymbol\\theta)\\rangle$. The energy expectation is\n",
+    "$\\langle E\\rangle = \\langle\\psi|H|\\psi\\rangle = \\sum_{(i,j)} c_{ij}\\,\\langle\\psi|\\sigma_i\\otimes\\sigma_j|\\psi\\rangle$.  In other words, we compute the expectation of each Pauli term and sum them .  The code below constructs the ansatz state and computes \\langle H\\rangle:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "6e6f9366",
+   "metadata": {
+    "collapsed": false,
+    "editable": true
+   },
+   "outputs": [],
+   "source": [
+    "# Define single-qubit rotation gates\n",
+    "def RY(theta):\n",
+    "    return np.array([[np.cos(theta/2), -np.sin(theta/2)],\n",
+    "                     [np.sin(theta/2),  np.cos(theta/2)]])\n",
+    "def RZ(theta):\n",
+    "    return np.array([[np.exp(-1j*theta/2), 0],\n",
+    "                     [0, np.exp(1j*theta/2)]])\n",
+    "\n",
+    "# 4x4 CNOT (control=qubit0, target=qubit1)\n",
+    "CNOT = np.array([[1,0,0,0],\n",
+    "                 [0,1,0,0],\n",
+    "                 [0,0,0,1],\n",
+    "                 [0,0,1,0]])\n",
+    "\n",
+    "# Initial state |00> as a vector\n",
+    "psi0 = np.kron(np.array([1,0]), np.array([1,0]))\n",
+    "\n",
+    "# Function to prepare ansatz state and compute energy\n",
+    "def energy_of_params(params):\n",
+    "    theta0, theta1, phi0, phi1 = params\n",
+    "    # Apply RY rotations on each qubit\n",
+    "    U1 = np.kron(RY(theta0), RY(theta1))\n",
+    "    psi = U1.dot(psi0)\n",
+    "    # Apply RZ rotations on each qubit\n",
+    "    U2 = np.kron(RZ(phi0), RZ(phi1))\n",
+    "    psi = U2.dot(psi)\n",
+    "    # Apply entangling CNOT gate\n",
+    "    psi = CNOT.dot(psi)\n",
+    "    # Compute ⟨ψ|H|ψ⟩\n",
+    "    return np.real(np.vdot(psi, H.dot(psi)))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8d0b3b41",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "Here energy$\\_$of$\\_$params(params) returns the variational energy $\\langle\\psi|H|\\psi\\rangle$.\n",
+    "This function automatically sums over the Pauli terms because we directly use the matrix H; one could\n",
+    "equivalently $sum c_{ij}\\,\\langle\\psi|\\sigma_i\\otimes\\sigma_j|\\psi\\rangle$ term-by-term.\n",
+    "By the variational principle, this expectation is always ≥ the true ground-state energy."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "adae4dd0",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "## Classical Optimization\n",
+    "\n",
+    "We now perform a classical minimization over the parameters $\\{\\theta_0,\\theta_1,\\phi_0,\\phi_1\\}$ to find the lowest energy.  We can use, for example, the Nelder–Mead routine from SciPy or a simple gradient descent.  For illustration we use SciPy’s minimize (Nelder–Mead) on our 4-parameter ansatz:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "48935bc7",
+   "metadata": {
+    "collapsed": false,
+    "editable": true
+   },
+   "outputs": [],
+   "source": [
+    "from scipy.optimize import minimize\n",
+    "\n",
+    "# Random initial guess for parameters\n",
+    "init_params = np.random.rand(4) * 2*np.pi\n",
+    "\n",
+    "# Minimize the energy\n",
+    "res = minimize(energy_of_params, init_params, method='Nelder-Mead')\n",
+    "theta_opt = res.x\n",
+    "E_vqe = res.fun\n",
+    "print(\"Optimized parameters:\", theta_opt)\n",
+    "print(\"VQE ground state energy (approx):\", E_vqe)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3b27191a",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "After optimization, res.fun holds the minimal energy found.  We also compare to the exact ground-state energy (the lowest eigenvalue of H) computed by direct diagonalization :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "6aa2f677",
+   "metadata": {
+    "collapsed": false,
+    "editable": true
+   },
+   "outputs": [],
+   "source": [
+    "# Exact diagonalization for verification\n",
+    "eigvals = np.linalg.eigvalsh(H)\n",
+    "E_exact = np.min(eigvals)\n",
+    "print(\"Exact ground state energy (lowest eigenvalue):\", E_exact)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bb7610ee",
+   "metadata": {
+    "editable": true
+   },
+   "source": [
+    "For the sample parameter values above, the output might look like:\n",
+    "\n",
+    "In summary: the code defines the Hamiltonian and decomposes it into Pauli tensor terms , constructs a parameterized ansatz (RY, RZ, CNOT) , computes the expectation values of each term , and uses a classical optimizer to minimize the energy. The returned E_vqe is the approximate ground-state energy, which matches the exact eigenvalue for this two-qubit example (as expected for a sufficiently expressive ansatz and a well-behaved optimization)  ."
+   ]
+  }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}