examples/QAOA_portfolio_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fcdf8dbb-4ccc-41af-9fb1-64362589db06",
   "metadata": {},
   "source": [
    "# Portfolio optimization\n",
    "\n",
    "This example shows how to generate the QAOA with customized initial state and mixers for portfolio optimization\n",
    "\n",
    "It contains the example of QAOA with Dicke state as the initial state and ring-XY mixer (https://arxiv.org/abs/2305.03857)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d0d23bd-ed12-4b6e-b44f-046de43641f7",
   "metadata": {},
   "source": [
    "The portfolio problem formulation is as below\n",
    "\\begin{aligned}\n",
    "\\min_{\\mathbf{x} \\in \\{0,1\\}^N} f(\\mathbf{x}) & = q \\mathbf{x}^T \\mathbf{W} \\mathbf{x} - \\boldsymbol{\\mu}^T \\mathbf{x}, \\\\\n",
    "\\mathop{\\mathrm{s.t.}} \\quad & \\mathbf{1}^T \\mathbf{x} = K    \n",
    "\\end{aligned}\n",
    "N is the number of assets, $\\boldsymbol{\\mu}$ is a vector of expected returns, $\\mathbf{W}$ is the covariance matrix between assets, q controls the trade-off between return and risk.\n",
    "The constraint is the Hamming-weight constraint, where K denotes how many assets we can buy (equivalent to the budget constraint under the assumption of every asset being equally priced)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed815229-c0a0-4d8d-8886-c06f5e5a79c4",
   "metadata": {},
   "source": [
    "# Generate (scaled) portfolio optimization problem with Qiskit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c3de09e7-5315-4170-8263-6041a7819fac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qokit.portfolio_optimization import get_problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35694327-b15e-400b-98b5-2387473232ed",
   "metadata": {},
   "source": [
    "For weighted combinatorial optimization, a global scaling on the objective function is beneficial for the parameter optimziaiton. \n",
    "The scaling rule is from Eq. (8) of https://arxiv.org/abs/2305.15201"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "57544f68-6538-4dda-a47b-254252caa54a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "po_problem = get_problem(N=6,K=3,q=0.5,seed=1,pre=1)\n",
    "means_in_spins = np.array([po_problem['means'][i] - po_problem['q'] * np.sum(po_problem['cov'][i, :]) for i in range(len(po_problem['means']))])\n",
    "scale = 1 / np.sqrt(np.mean((( po_problem['q']*po_problem['cov'])**2).flatten())+np.mean((means_in_spins**2).flatten()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e49f9152-e49c-4cbd-a2b3-20894a927a9a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "po_problem = get_problem(N=6,K=3,q=0.5,seed=1,pre=scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a800d2c6-d7a9-449a-befc-e1620676847e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "po_problem2 = get_problem(N=6,K=3,q=0.5,seed=1,pre='rule')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "60942736-20e4-4519-a584-65f613b7bddf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# confirm that the scaling rule in the function matches the one above\n",
    "assert np.allclose(po_problem['cov'], po_problem2['cov'])\n",
    "assert np.allclose(po_problem['means'], po_problem2['means'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dcaa734-b9bd-4d23-be18-c7cbd4892525",
   "metadata": {},
   "source": [
    "# Get QAOA objective and evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c4cb1e0a-f0eb-48e4-a1dc-8d342462c156",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qokit.portfolio_optimization import portfolio_brute_force, get_sk_ini\n",
    "from qokit.qaoa_objective_portfolio import get_qaoa_portfolio_objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1fa46296-89cb-40e7-aaf8-09b708eda12d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "p = 1\n",
    "qaoa_obj = get_qaoa_portfolio_objective(po_problem=po_problem,p=p,ini='dicke',mixer='trotter_ring',T=1,simulator='python')\n",
    "best_portfolio = portfolio_brute_force(po_problem,return_bitstring=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4bb70cb8-db79-45b3-b808-96b7d6647b68",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "x0 = get_sk_ini(p=p) \n",
    "# Alternative: random initial point# x0 = np.random.rand(2*p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "238479eb-db75-4e19-a159-36a2bed9c2ae",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energy = 0.004782805426149568, Approximation ratio = 0.767014699611637\n"
     ]
    }
   ],
   "source": [
    "po_energy = qaoa_obj(x0).real\n",
    "po_ar = (po_energy-best_portfolio[1])/(best_portfolio[0]-best_portfolio[1])\n",
    "print(f\"energy = {po_energy}, Approximation ratio = {po_ar}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "800fbb5c-29fd-4f29-971c-99802daab5cf",
   "metadata": {},
   "source": [
    "# Optimize QAOA parameters\n",
    "\n",
    "Note that we are using NLopt for optimization as it supports better-performing BOBYQA optimizer. \n",
    "\n",
    "Run `pip install nlopt` to install this dependency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9d730369-3d9a-4e0e-bd2b-a81d92baa7b7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "import nlopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "55a761e4-2122-4b00-bc1f-39730d6a3dd7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def minimize_nlopt(f, x0, rhobeg=None, p=None):\n",
    "    def nlopt_wrapper(x, grad):\n",
    "        if grad.size > 0:\n",
    "            sys.exit(\"Shouldn't be calling a gradient!\")\n",
    "        return f(x).real\n",
    "\n",
    "    opt = nlopt.opt(nlopt.LN_BOBYQA, 2 * p)\n",
    "    opt.set_min_objective(nlopt_wrapper)\n",
    "\n",
    "    opt.set_xtol_rel(1e-8)\n",
    "    opt.set_ftol_rel(1e-8)\n",
    "    opt.set_initial_step(rhobeg)\n",
    "    xstar = opt.optimize(x0)\n",
    "    minf = opt.last_optimum_value()\n",
    "\n",
    "    return xstar, minf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "86b7ff7f-43f4-4242-95b9-82784f0632a5",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energy = 0.004135567715727546, Approximation ratio = 0.8050399835274252\n"
     ]
    }
   ],
   "source": [
    "_, opt_energy = minimize_nlopt(qaoa_obj, x0, p=1, rhobeg=0.01/1)\n",
    "opt_ar = (opt_energy-best_portfolio[1])/(best_portfolio[0]-best_portfolio[1])\n",
    "print(f\"energy = {opt_energy}, Approximation ratio = {opt_ar}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "90c5c6af-003c-4724-a943-c747bd24f3c6",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energy = 0.004135567763346608, Approximation ratio = 0.8050399807298003\n"
     ]
    }
   ],
   "source": [
    "res = minimize(qaoa_obj, x0, method='COBYLA', options={'rhobeg':0.001})\n",
    "print(f\"energy = {res.fun}, Approximation ratio = {(res.fun-best_portfolio[1])/(best_portfolio[0]-best_portfolio[1])}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7af5e53e-2703-4f4b-9554-efe8196a14a0",
   "metadata": {},
   "source": [
    "# Get qiskit circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0c6e7780-0392-4a9f-b0a7-4b23c0af84dc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qokit.qaoa_circuit_portfolio import get_qaoa_circuit, get_parameterized_qaoa_circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d5fec9f4-09e0-409f-a613-07f23194c94d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "qc = get_qaoa_circuit(po_problem, gammas=x0[:p]/2, betas=x0[p:]/2, depth=p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30f21f20-f6ec-4fd8-86ee-3f3ba7155e78",
   "metadata": {},
   "source": [
    "# Verify the qiskit circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07154688-aa49-42e2-b1a0-b3c7c400537d",
   "metadata": {},
   "source": [
    "It is important to verify the qiskit circuit simulation with fur-backend simulations for hardware experiments. For portfolio optimization, there is a 1/2 scaling of (gamma, beta) in the qiskit circuit. \n",
    "Also see: `tests/test_portfolio_optimization/test_portfolio_qokitandqiskit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bc26bf47-72ce-4944-9a9f-d217e4c77894",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit import execute, Aer\n",
    "from qokit.utils import reverse_array_index_bit_order\n",
    "from qokit.qaoa_circuit_portfolio import get_energy_expectation_sv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0c4f35dc-9402-48ce-89ef-c5a80883aaf8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "backend = Aer.get_backend(\"statevector_simulator\")\n",
    "qc = get_qaoa_circuit(po_problem, gammas=x0[:p]/2, betas=x0[p:]/2, depth=p)\n",
    "result = execute(qc, backend).result()\n",
    "sv1 = reverse_array_index_bit_order(result.get_statevector())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8581ac1f-2675-4ccb-b311-ee7e658b2782",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_25329/3897112670.py:2: DeprecationWarning: The method ``qiskit.circuit.quantumcircuit.QuantumCircuit.bind_parameters()`` is deprecated as of qiskit 0.45.0. It will be removed no earlier than 3 months after the release date. Use assign_parameters() instead\n",
      "  qc2 = parameterized_qc.bind_parameters(np.hstack([x0[p:]/2, x0[:p]/2]))\n"
     ]
    }
   ],
   "source": [
    "parameterized_qc = get_parameterized_qaoa_circuit(po_problem, depth=p)\n",
    "qc2 = parameterized_qc.bind_parameters(np.hstack([x0[p:]/2, x0[:p]/2]))\n",
    "result = execute(qc2, backend).result()\n",
    "sv2 = reverse_array_index_bit_order(result.get_statevector())\n",
    "assert np.allclose(sv1, sv2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "2ceed19a-9391-4eed-95b9-dbfc3b8d6a9d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "qaoa_obj_qiskit = get_qaoa_portfolio_objective(po_problem=po_problem, p=p, simulator=\"qiskit\", ini=\"dicke\", mixer=\"trotter_ring\", T=1)\n",
    "qaoa_obj_auto = get_qaoa_portfolio_objective(po_problem=po_problem, p=p, simulator=\"auto\", ini=\"dicke\", mixer=\"trotter_ring\", T=1)\n",
    "assert np.allclose(qaoa_obj_auto(x0), qaoa_obj_qiskit(x0))\n",
    "assert np.allclose(get_energy_expectation_sv(po_problem,sv1), qaoa_obj_qiskit(x0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64997624-2c3f-4b24-85b6-0a6589b2f034",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qokit",
   "language": "python",
   "name": "qokit"
  },
  "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}