12 - hamiltonian_vqe_qaoa.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.opflow import X, Y, Z, I, MatrixEvolution, PauliTrotterEvolution, PauliExpectation, StateFn\n",
    "from qiskit.algorithms import VQE, QAOA\n",
    "from qiskit.algorithms.optimizers import COBYLA\n",
    "from qiskit.utils import QuantumInstance, algorithm_globals\n",
    "from qiskit import Aer, QuantumCircuit\n",
    "from qiskit.circuit.library import RealAmplitudes, EfficientSU2, ZFeatureMap\n",
    "from qiskit.circuit import Parameter\n",
    "from qiskit.opflow import CircuitSampler\n",
    "import numpy as np\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hamiltonian time evolution simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\hat{H} = \\frac{1}{2} \\left( \\hat{I}\\otimes \\hat{I} +  \\hat{\\sigma}_x \\otimes \\hat{\\sigma}_x + \\hat{\\sigma}_y \\otimes \\hat{\\sigma}_y + \\hat{\\sigma}_z \\otimes \\hat{\\sigma}_z \\right) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "operator_H = ((I ^ I) + (X ^ X) + (Y ^ Y) + (Z ^ Z))*1/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5 * II\n",
      "+ 0.5 * XX\n",
      "+ 0.5 * YY\n",
      "+ 0.5 * ZZ\n"
     ]
    }
   ],
   "source": [
    "print(operator_H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 1.+0.j 0.+0.j]\n",
      " [0.+0.j 1.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 1.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "H_matrix = operator_H.to_matrix()\n",
    "\n",
    "print(H_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = Parameter('t')\n",
    "\n",
    "time_H = (t * operator_H).exp_i()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "e^(-i*1.0*t * (\n",
      "  0.5 * II\n",
      "  + 0.5 * XX\n",
      "  + 0.5 * YY\n",
      "  + 0.5 * ZZ\n",
      "))\n"
     ]
    }
   ],
   "source": [
    "print(time_H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.87758256-0.47942554j 0.        +0.j         0.        +0.j\n",
      "  0.        +0.j        ]\n",
      " [0.        +0.j         0.87758256+0.j         0.        -0.47942554j\n",
      "  0.        +0.j        ]\n",
      " [0.        +0.j         0.        -0.47942554j 0.87758256+0.j\n",
      "  0.        +0.j        ]\n",
      " [0.        +0.j         0.        +0.j         0.        +0.j\n",
      "  0.87758256-0.47942554j]]\n",
      "     ┌──────────────┐\n",
      "q_0: ┤0             ├\n",
      "     │  Hamiltonian │\n",
      "q_1: ┤1             ├\n",
      "     └──────────────┘\n"
     ]
    }
   ],
   "source": [
    "evolution_time = 0.5\n",
    "\n",
    "mat_evo = MatrixEvolution()\n",
    "\n",
    "result_mat_exp = mat_evo.convert(time_H)\n",
    "result_mat_exp = result_mat_exp.bind_parameters({t: evolution_time})\n",
    "\n",
    "print(result_mat_exp.to_matrix())\n",
    "print(result_mat_exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CircuitStateFn(\n",
      "     ┌──────────────────────┐\n",
      "q_0: ┤0                     ├\n",
      "     │  Initialize(0,0,1,0) │\n",
      "q_1: ┤1                     ├\n",
      "     └──────────────────────┘\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "observable = X ^ X\n",
    "\n",
    "observable_measurement = StateFn(observable).adjoint()\n",
    "\n",
    "eigenvalues, eigenstates = np.linalg.eigh(operator_H.to_matrix())\n",
    "initial_state = StateFn(eigenstates[0])\n",
    "\n",
    "print(initial_state.to_circuit_op())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  OperatorMeasurement(XX),\n",
      "  e^(-i*1.0*t * (\n",
      "    0.5 * II\n",
      "    + 0.5 * XX\n",
      "    + 0.5 * YY\n",
      "    + 0.5 * ZZ\n",
      "  )),\n",
      "  VectorStateFn(Statevector([0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],\n",
      "              dims=(2, 2)))\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "evo_and_measure = observable_measurement @ time_H @ initial_state\n",
    "\n",
    "print(evo_and_measure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  OperatorMeasurement(XX),\n",
      "       ┌──────────────────────────────────────┐\n",
      "  q_0: ┤0                                     ├\n",
      "       │  exp(-it (II + XX + YY + ZZ))(1.0*t) │\n",
      "  q_1: ┤1                                     ├\n",
      "       └──────────────────────────────────────┘,\n",
      "  VectorStateFn(Statevector([0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],\n",
      "              dims=(2, 2)))\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "pauli_trot = PauliTrotterEvolution(trotter_mode='trotter', reps = 1)\n",
    "\n",
    "trotterized_op = pauli_trot.convert(evo_and_measure)\n",
    "\n",
    "print(trotterized_op)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observable at time [5.5, 0.75]: [ 0.+0.j -0.+0.j]\n"
     ]
    }
   ],
   "source": [
    "# Run Hamiltoninan in ideal device\n",
    "evo_time_points = [5.5, 0.75]\n",
    "\n",
    "hamiltonian_trotter_expectations = trotterized_op.bind_parameters({t: evo_time_points})\n",
    "\n",
    "print(f\"Observable at time {evo_time_points}: {np.round(hamiltonian_trotter_expectations.eval(), 3)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energies:  [1.    0.996]\n"
     ]
    }
   ],
   "source": [
    "# Run Hamiltoninan in simulatoror device\n",
    "\n",
    "sampler = CircuitSampler(backend=Aer.get_backend(\"qasm_simulator\"))\n",
    "\n",
    "sampled_trotter_exp_op = sampler.convert(hamiltonian_trotter_expectations)\n",
    "sampled_trotter_energies = sampled_trotter_exp_op.eval()\n",
    "\n",
    "print(\"Energies: \", np.round(np.real(sampled_trotter_energies), 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# VQE for Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = Aer.get_backend('statevector_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "quantum_instance = QuantumInstance(backend, \n",
    "                                   shots = 1024, \n",
    "                                   seed_simulator = algorithm_globals.random_seed, \n",
    "                                   seed_transpiler = algorithm_globals.random_seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ansatz = RealAmplitudes(num_qubits = 1, reps = 4)\n",
    "ansatz = ZFeatureMap(1, reps = 4)\n",
    "#ansatz = EfficientSU2(num_qubits = 1, reps = 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = COBYLA(maxiter = 120, tol = 0.001)\n",
    "\n",
    "vqe = VQE(ansatz=ansatz, quantum_instance=quantum_instance, optimizer=optimizer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum energy eigenvalue: (0.1819510439424697+0j)\n"
     ]
    }
   ],
   "source": [
    "minimum_energy = vqe.compute_minimum_eigenvalue(operator_H).eigenvalue\n",
    "\n",
    "print(\"Minimum energy eigenvalue:\", minimum_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QAOA for Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenstate:\n",
      "\n",
      " [-0.48661084-0.11493429j -0.48661084-0.11493429j -0.48661084-0.11493429j\n",
      " -0.48661084-0.11493429j]\n",
      "\n",
      "Minimum energy eigenvalue:\n",
      "\n",
      " (1.0000000000000004+0j)\n"
     ]
    }
   ],
   "source": [
    "qaoa = QAOA(optimizer, quantum_instance=Aer.get_backend('statevector_simulator'))\n",
    "\n",
    "result = qaoa.compute_minimum_eigenvalue(operator_H)\n",
    "\n",
    "minimum_energy = result.eigenstate\n",
    "\n",
    "print(\"Eigenstate:\\n\\n\", minimum_energy)\n",
    "print(\"\\nMinimum energy eigenvalue:\\n\\n\",  result.eigenvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}