Meta-VQE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Meta-VQE algorithm (demo) <img src=\"img/ML_logo.png\" width=\"20%\" align=\"right\">\n",
    "_______________________________________________________\n",
    "\n",
    "_Authors:_ Alba Cervera-Lierta, Jakob S. Kottmann, Alán Aspuru-Guzik.\n",
    "\n",
    "This notebook presents a demo code to reproduce the results of the article \"The Meta-Variational Quantum Eigensolver (Meta-VQE): Learning energy profiles of parameterized Hamiltonians for quantum simulation\" (http://arxiv.org/abs/2009.13545).\n",
    "\n",
    "This demo simulates a meta-VQE for the 1D XXZ model with external field. It also shows the comparison with a standard VQE, opt-meta-VQE, GA-VQE and opt-GA-VQE algorithms detailed in the main article. \n",
    "\n",
    "For computational reasons, the notebook is set to simulate a $n=4$ spin chain and a meta-VQE with $2$ encoding layer and $1$ processing layers. Results shown in the main article can be obtained by adjusting the number of qubits and layers. Article results can be found in \"data\" (link) folder in this repository.\n",
    "\n",
    "The code is written using [_Tequila_](https://github.com/aspuru-guzik-group/tequila) language and [_qulacs_](https://github.com/qulacs/qulacs) simulator as a default backend.\n",
    "\n",
    "__Table of contents:__\n",
    "1. [Model and PQC ansatz](#Model-and-PQC-ansatz)\n",
    "    1. [Hamiltonian: XXZ with external field](#Hamiltonian:-XXZ-with-external-field)\n",
    "    1. [Parametrized Quantum Circuit](#Parametrized-Quantum-Circuit)\n",
    "1. [Main code: meta-VQE algorithm](#Main-code:-meta-VQE-algorithm)\n",
    "    1. [Meta-VQE](#Meta-VQE)\n",
    "        1. [Training](#Training)\n",
    "        1. [Test](#Test)\n",
    "    1. [Comparison with other algorithms](#Comparison-with-other-algorithms)\n",
    "        1. [Standard VQE](#Standard-VQE)\n",
    "        1. [Globally-averaged VQE (GA-VQE)](#Globally-averaged-VQE-(GA-VQE))\n",
    "        1. [VQE with meta-VQE initialization (opt-meta-VQE)](#VQE-with-meta-VQE-initialization-(opt-meta-VQE))\n",
    "        1. [VQE with GA-VQE initialization (opt-GA-VQE)](#VQE-with-GA-VQE-initialization-(opt-GA-VQE))\n",
    "    1. [Plots](#Plots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import all we need\n",
    "# tested with tq git revision: ab252adc4090465ff47f182de32d69b84d7fc13f\n",
    "import tequila as tq\n",
    "import numpy as np\n",
    "from numpy import random\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Seed generator from alphabetic string\n",
    "def seed_gen(string):\n",
    "    num_str = [ord(letter) - 96 for letter in string]\n",
    "    # seed must be between 0 and 2**32 - 1\n",
    "    return(np.mod(int(''.join(map(str,num_str))),2**32 - 1))\n",
    "\n",
    "seeds = ['matterlab', 'alba', 'jakob', 'alan']\n",
    "\n",
    "s = 0 # default seed, used in the main article results\n",
    "random.seed(seed_gen(seeds[s])) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model and PQC ansatz\n",
    "\n",
    "First, we will define the Hamiltonian and circuit ansatz used in this work. These functions can be easily modified to accomodate other hamiltonians or ansatzes.\n",
    "\n",
    "### Hamiltonian: XXZ with external field\n",
    "\n",
    "The spin Hamiltonian used is the 1D antiferromagnetic XXZ spin chain with external field $\\lambda$ and periodic boundary conditions:\n",
    "$$ \\mathcal{H}=\\sum_{i=1}^{n}\\sigma_{i}^{x}\\sigma_{i+1}^{x}+\\sigma_{i}^{y}\\sigma_{i+1}^{y}+\\Delta\\sigma_{i}^{z}\\sigma_{i+1}^{z} + \\lambda \\sum_{i=1}^{n}\\sigma_{i}^{z}.$$\n",
    "To check the accuracy of the results obtained, we also compute the exact ground state energy of this Hamiltonian by diagonalizing it; we called that function `exact`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def HXXZ(num_qubits,delta,lam):\n",
    "    Lx = int(np.sqrt(num_qubits))\n",
    "    ham = tq.paulis.X(0)*tq.paulis.X(1)\n",
    "    ham += tq.paulis.Y(0)*tq.paulis.Y(1)\n",
    "    ham += tq.paulis.X(0)*tq.paulis.X(1+Lx)\n",
    "    ham += tq.paulis.Y(0)*tq.paulis.Y(1+Lx)\n",
    "    ham += delta*tq.paulis.Z(0)*tq.paulis.Z(1)\n",
    "    ham += delta*tq.paulis.Z(0)*tq.paulis.Z(1+Lx)\n",
    "    for i in range(num_qubits):\n",
    "        if i+1 % Lx != 0 and i != 0:\n",
    "            ham += tq.paulis.X(i)*tq.paulis.X(i+1)\n",
    "            ham += tq.paulis.Y(i)*tq.paulis.Y(i+1)\n",
    "            ham += delta*tq.paulis.Z(i)*tq.paulis.Z(i+1)\n",
    "        if i+Lx < num_qubits and i != 0:\n",
    "            ham += tq.paulis.X(i)*tq.paulis.X(i+Lx)\n",
    "            ham += tq.paulis.Y(i)*tq.paulis.Y(i+Lx)\n",
    "            ham += delta*tq.paulis.Z(i)*tq.paulis.Z(i+Lx)\n",
    "        ham += lam*tq.paulis.Z(i)\n",
    "    return(ham)\n",
    "def exact(num_qubits,delta,lam):\n",
    "    ham_matrix = HXXZ(num_qubits,delta,lam).to_matrix()\n",
    "    energ = np.linalg.eigvals(ham_matrix)\n",
    "    return(min(energ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+2.0000X(0)X(1)+2.0000Y(0)Y(1)+1.0000X(0)X(4)+1.0000Y(0)Y(4)+1.0000X(0)X(3)+1.0000Y(0)Y(3)+1.0000X(1)X(2)+1.0000Y(1)Y(2)+1.0000X(1)X(4)+1.0000Y(1)Y(4)+1.0000X(2)X(3)+1.0000Y(2)Y(3)+1.0000X(2)X(5)+1.0000Y(2)Y(5)+1.0000X(3)X(4)+1.0000Y(3)Y(4)+1.0000X(3)X(6)+1.0000Y(3)Y(6)+1.0000X(4)X(5)+1.0000Y(4)Y(5)+1.0000X(4)X(7)+1.0000Y(4)Y(7)+1.0000X(5)X(6)+1.0000Y(5)Y(6)+1.0000X(5)X(8)+1.0000Y(5)Y(8)+1.0000X(6)X(7)+1.0000Y(6)Y(7)+1.0000X(7)X(8)+1.0000Y(7)Y(8)+1.0000X(8)X(9)+1.0000Y(8)Y(9)\n"
     ]
    }
   ],
   "source": [
    "print(HXXZ(9,0,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parametrized Quantum Circuit\n",
    "\n",
    "The meta-VQE algorithm needs a Parameterized Quantum Circuit (PQC) both for the encoding and processsing parts. We use a unique ansatz for both of them, where the parameters appear in the single-qubit rotational gates. The ansatz is divided in layers, each of them including single-qubit gates on each qubit plus CNOT gates between neighbouring pairs. \n",
    "\n",
    "First, we define the general single-qubit unitary gater `Ugen`,\n",
    "$$ U_{gen} = R_{z}(\\phi)R_{y}(\\theta)R_{z}(\\varphi). $$\n",
    "\n",
    "The `PQC` depends on the number of qubits, layers, Hamiltonian argument and the optimization variable name. If the argument is zero (which will be the case for the processing layers, where there is no encoding), then it only uses `Ugen` gates with two optimization angles. On the contrary, if the arguments are different from zero (corresponding to the encoding layers), the `Ugen` gate angles will be defined as\n",
    "$$ R_{z}(w_{il}^{(1)}\\Delta + \\phi_{il}^{(1)}) R_{y}(w_{il}^{(2)}\\Delta + \\phi_{il}^{(2)}), $$\n",
    "where $i$ and $l$ correspond to the qubit and layer respectively.\n",
    "Then, `PQC` will output a quantum circuit with $4L$ parameters for $L$ encoding layers, and $2L$ parameters for $L$ processing layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# General 1-qubit gate\n",
    "# q0 = target qubit\n",
    "def Ugen(th, phi, lam, q0):\n",
    "    ugate = tq.gates.Rz(target=q0, angle=phi) + tq.gates.Ry(target=q0, angle=th) + tq.gates.Rz(target=q0, angle=lam)\n",
    "    return ugate\n",
    "\n",
    "def PQC(num_qubits, layers, arg, var_name):   \n",
    "    U = tq.QCircuit()\n",
    "    arg_tmp = [arg, arg]\n",
    "        \n",
    "    for l in range(layers):\n",
    "        # Layer single-qubit gates\n",
    "        for q in range(num_qubits):\n",
    "            theta = tq.Variable(name=\"{}_{}{}\".format(var_name, l, q))\n",
    "            phi = tq.Variable(name=\"{}2_{}{}\".format(var_name, l, q))\n",
    "            if arg_tmp[0] == 0 and arg_tmp[1] == 0:\n",
    "                U += Ugen(theta, phi, 0, q)  # for no encoding, reduce the total number of variables\n",
    "            else:\n",
    "                wth = tq.Variable(name=\"{}w_{}{}\".format(var_name, l, q))\n",
    "                wph = tq.Variable(name=\"{}2w_{}{}\".format(var_name, l, q))\n",
    "                U += Ugen(wth * arg_tmp[0] + theta, wph * arg_tmp[1] + phi, 0, q)\n",
    "        # Layer alternating CNOTs\n",
    "        ent0 = np.mod(l, 2)\n",
    "        for ent in range(ent0, num_qubits - 1, 2):\n",
    "            U += tq.gates.CNOT(control=ent, target=ent + 1)\n",
    "    return (U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main code: meta-VQE algorithm\n",
    "\n",
    "First, we define the main characteristics of the meta-VQE algorithm and the tests that we would like to run to compare its performace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PQC ansatz: same for encoding and processing layers\n",
    "ans_enc = PQC # ansatz variational circuit\n",
    "ans_var = PQC # ansatz encoding\n",
    "\n",
    "n_qub = 9 # qubits\n",
    "n_lay = [2,1] # layers [encoding, processing]\n",
    "\n",
    "# Hamiltonian argument limits\n",
    "arg_max = -1.1\n",
    "arg_min = 1.1\n",
    "lamXXZ = 0.75 # we set lambda, so the code only works with a single meta-parameter\n",
    "\n",
    "# Training and test points\n",
    "n_train = 20 #20 \n",
    "n_test = 100 \n",
    "equispaced_train = True # if False, points are distributed at random\n",
    "equispaced_test = False\n",
    "\n",
    "# Optimizer options\n",
    "methods = 'BFGS'\n",
    "grad_methods = '2-point'\n",
    "backend = 'qulacs'\n",
    "lr = 0.01\n",
    "# make sure scipy version is >= 1.5\n",
    "mthd_opt = {'finite_diff_rel_step': 0.0001}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Meta-VQE\n",
    "\n",
    "As discussed in the main article, the cost function used to train the meta-VQE is the sum of the Hamiltonian expectation values w.r.t the parameterized quantum circuit initialized in the $|0\\rangle$ state:\n",
    "$$ \\mathcal{L}oss(\\vec{\\Phi},\\vec{\\Theta}) = \\sum_{i=1}^{M}\\langle\\psi_{i}\\lvert H(\\Delta_{i})|\\psi_{i}\\rangle, $$\n",
    "where\n",
    "$$ |\\psi_{i}\\rangle \\equiv |\\psi(\\Delta_{i},\\vec{\\Phi},\\vec{\\Theta})\\rangle = \\mathcal{U}(\\vec{\\Theta})\\mathcal{S}(\\Delta_{i},\\vec{\\Phi})|0\\rangle^{\\otimes n},$$\n",
    "and $\\Delta_{i}$ the Hamiltonian $\\Delta$ parameters from the training points, $\\mathcal{U}$ and $\\mathcal{S}$ the PQC corresponding to processing and encoding layers respectiverly (in these simulations, they have the same circuit structure, see above), and $\\vec{\\Phi},\\vec{\\Theta}$ are the optimization parameters aranged in a vector form.\n",
    "\n",
    "#### Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the training points\n",
    "if equispaced_train == True:\n",
    "    arg_train = [arg_min + i*(arg_max-arg_min)/(n_train-1) for i in range(n_train)]\n",
    "else:\n",
    "    arg_train = [random.uniform(arg_min, arg_max) for i in range(n_train)]\n",
    "    arg_train.sort()\n",
    "\n",
    "t0_metaVQE = time.time()\n",
    "# Construct the objective (cost function)\n",
    "Obj = tq.Objective()\n",
    "for i in range(n_train):\n",
    "    ham = HXXZ(n_qub, arg_train[i], lamXXZ)\n",
    "    # meta-VQE circuit: \n",
    "    # encoding variables share the prefix \"th_en\" and processing variables, the prefix \"th\"\n",
    "    total_U = ans_enc(n_qub, n_lay[0], arg_train[i], \"th_en\") + ans_var(n_qub, n_lay[1], 0, \"th\")\n",
    "\n",
    "    Obj += tq.ExpectationValue(H=ham, U=total_U)\n",
    "\n",
    "variables = Obj.extract_variables()\n",
    "variables = sorted(variables, key=lambda x: x.name)\n",
    "\n",
    "# Random initialization of variables\n",
    "th0 = {key: random.uniform(0, np.pi) for key in variables}\n",
    "\n",
    "initial_values = th0\n",
    "metaVQE = tq.minimize(objective=Obj, adaptive = True, lr=lr, method_options=mthd_opt, method=methods, gradient=grad_methods, samples=None,\n",
    "                      initial_values=initial_values, backend=backend, noise=None, device=None, silent=True)\n",
    "\n",
    "t1_metaVQE = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable `metaVQE` contains the result of the optimization, i.e. the final values of the loss function and optimization angles $\\vec{\\Phi}_{opt},\\vec{\\Theta}_{opt}$, together with optimization history. With these values, we can check what is the result of the training points and use them for the test points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t2_metaVQE = time.time()\n",
    "# Results training\n",
    "x_train = arg_train\n",
    "y_train_ex = []\n",
    "y_train_metaVQE = []\n",
    "error_train_metaVQE = []\n",
    "\n",
    "for i in range(n_train):\n",
    "    ham = HXXZ(n_qub, arg_train[i], lamXXZ)\n",
    "\n",
    "    # construct the corresponding PQC with the same variable names\n",
    "    total_U = ans_enc(n_qub, n_lay[0], arg_train[i], \"th_en\") + ans_var(n_qub, n_lay[1], 0, \"th\")\n",
    "\n",
    "    exp_val = tq.ExpectationValue(H=ham, U=total_U)\n",
    "    # Use the angles obtained in the training part\n",
    "    res = tq.simulate(exp_val, variables=metaVQE.angles) \n",
    "    # Exact result to compare\n",
    "    res_ex = exact(n_qub,arg_train[i], lamXXZ)\n",
    "\n",
    "    # print(arg_train[i], res_ex, res, abs(res-res_ex))\n",
    "\n",
    "    y_train_ex.append(res_ex)\n",
    "    y_train_metaVQE.append(res)\n",
    "    error_train_metaVQE.append(abs(res-res_ex))\n",
    "t3_metaVQE = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the testing points\n",
    "if equispaced_test == True:\n",
    "    arg_test = [arg_min + i*(arg_max-arg_min)/(n_test-1) for i in range(n_test)]\n",
    "    arg_test.sort()\n",
    "else:\n",
    "    arg_test = [random.uniform(arg_min, arg_max) for i in range(n_test)]\n",
    "    arg_test.sort()\n",
    "\n",
    "t4_metaVQE = time.time()\n",
    "x_test= arg_test\n",
    "y_test_ex = []\n",
    "y_test_metaVQE = []\n",
    "error_test_metaVQE = []\n",
    "\n",
    "Obj = tq.Objective()\n",
    "for i in range(n_test):\n",
    "    ham = HXXZ(n_qub, arg_test[i], lamXXZ)\n",
    "    # construct the corresponding PQC with the same variable names\n",
    "    total_U = ans_enc(n_qub, n_lay[0], arg_test[i], \"th_en\") + ans_var(n_qub, n_lay[1], 0, \"th\")\n",
    "\n",
    "    exp_val = tq.ExpectationValue(H=ham, U=total_U)\n",
    "    # Use the angles obtained in the training part\n",
    "    res = tq.simulate(exp_val, variables=metaVQE.angles)\n",
    "    # Exact result to compare\n",
    "    res_ex = exact(n_qub, arg_test[i], lamXXZ)\n",
    "\n",
    "    # print(arg_test[i], res_ex, res, abs(res - res_ex))\n",
    "\n",
    "    y_test_metaVQE.append(res)\n",
    "    y_test_ex.append(res_ex.real)\n",
    "    error_test_metaVQE.append(abs(res-res_ex.real))\n",
    "t5_metaVQE = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLOTS\n",
    "\n",
    "# g.s. energy\n",
    "plt.plot(x_test, y_test_ex, color=\"black\", ls=\":\", label=\"exact\")\n",
    "plt.plot(x_test, y_test_metaVQE,  ls=\"-\", label=\"metaVQE test\")\n",
    "plt.scatter(x_train, y_train_metaVQE, color=\"red\", marker=\"x\", label=\"metaVQE train\")\n",
    "plt.xlabel(\"Δ\")\n",
    "plt.ylabel('g.s. energy')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Absolute error\n",
    "plt.plot(x_test, error_test_metaVQE, ls=\"-\", label=\"metaVQE test\")\n",
    "plt.scatter(x_train, error_train_metaVQE, marker=\"x\", color=\"red\", label=\"metaVQE train\")\n",
    "plt.xlabel(\"Δ\")\n",
    "plt.ylabel(\"Absolute error\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison with other algorithms\n",
    "\n",
    "#### Standard VQE\n",
    "\n",
    "The first comparison test correponds with the well-known VQE algorithm. The standard VQE optimizes the values of the variables of the PQC for each Hamiltonian, i.e. for each value of $\\Delta$. Then, in contrast with meta-VQE, we have to run individual optimizations for each Hamiltonian point instead of a single big optimization. \n",
    "\n",
    "The other main difference between VQE and meta-VQE is the encoding part of the circuit. Standard VQE algorithms do not encode the variables of the Hamiltonian into the PQC. Then, to compare with the above meta-VQE results, we fix the circuit depth, i.e. the total number of layers is the same, but all single-qubit rotations contain only two optimization variables (as if all layers were processing layers).\n",
    "\n",
    "Due to the computational cost of running individual optimization for each value of $\\Delta$, we compare meta-VQE and VQE algorithms using the equispaced training points generated for meta-VQE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0_VQE = time.time()\n",
    "\n",
    "y_stVQE = []\n",
    "error_stVQE = []\n",
    "\n",
    "for i in range(n_train):\n",
    "    ham = HXXZ(n_qub, arg_train[i], lamXXZ)\n",
    "\n",
    "    # VQE circuit layers same as meta-VQE layers (encoding + variational parts)\n",
    "    total_U = ans_var(n_qub, n_lay[0]+n_lay[1], 0, \"th\")\n",
    "\n",
    "    Obj = tq.ExpectationValue(H=ham, U=total_U)\n",
    "\n",
    "    variables = Obj.extract_variables()\n",
    "    variables = sorted(variables, key=lambda x: x.name)\n",
    "\n",
    "    # Random initialization of variables\n",
    "    # new initialization because the circuit is different (no encoding part)\n",
    "    th0_bis = {key: random.uniform(0, np.pi) for key in variables}\n",
    "\n",
    "    initial_values = th0_bis\n",
    "    standardVQE = tq.minimize(objective=Obj, adaptive = True, lr=lr, method_options=mthd_opt, method=methods, gradient=grad_methods, samples=None,\n",
    "                              initial_values=initial_values, backend=backend, noise=None,\n",
    "                              device=None, silent=True)\n",
    "    res = standardVQE.energy\n",
    "    # res_ex = Exact results same as meta-VQE training points\n",
    "\n",
    "    # print(arg_train[i], y_train_ex[i], res, abs(res - y_train_ex[i]))\n",
    "\n",
    "    y_stVQE.append(res)\n",
    "    error_stVQE.append(abs(res-y_train_ex[i]))\n",
    "    \n",
    "t1_VQE = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Globally-averaged VQE (GA-VQE)\n",
    "\n",
    "To check the power of the encoding strategy of the meta-VQE, we run a meta-VQE-like algorithm without encoding the Hamiltonian variables. We call this algorithm globally-averaged VQE, since it has the same PQC strucutre of a standard VQE but the optimization is done globally, by using the meta-VQE cost function defined above.\n",
    "\n",
    "We use the same points used in the meta-VQE to train and test the GA-VQE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0_GAVQE = time.time()\n",
    "\n",
    "Obj = tq.Objective()\n",
    "for i in range(n_train):\n",
    "    ham = HXXZ(n_qub, arg_train[i], lamXXZ)\n",
    "\n",
    "    # same number of layers as VQE and meta-VQE (encoding + variational parts)\n",
    "    total_U = ans_var(n_qub, n_lay[0] + n_lay[1], 0, \"th\")\n",
    "\n",
    "    Obj += tq.ExpectationValue(H=ham, U=total_U)\n",
    "\n",
    "# same variables as a starndard VQE, same random initiaization\n",
    "initial_values = th0_bis\n",
    "GA_VQE = tq.minimize(objective=Obj, adaptive = True, lr=lr, method_options=mthd_opt, method=methods, gradient=grad_methods, samples=None,\n",
    "                            initial_values=initial_values, backend=backend, noise=None,\n",
    "                            device=None, silent=True)\n",
    "\n",
    "y_test_GAVQE = []\n",
    "error_test_GAVQE = []\n",
    "\n",
    "for i in range(n_test):\n",
    "    ham = HXXZ(n_qub, arg_test[i], lamXXZ)\n",
    "\n",
    "    # same number of layers as VQE and meta-VQE (encoding + variational parts)\n",
    "    total_U = ans_var(n_qub, n_lay[0] + n_lay[1], 0, \"th\")\n",
    "\n",
    "    exp_val = tq.ExpectationValue(H=ham, U=total_U)\n",
    "    res = tq.simulate(exp_val, variables=GA_VQE.angles)\n",
    "    # res_ex = same as meta-VQE test\n",
    "\n",
    "    # print(arg_test[i], y_test_ex[i], res, abs(res-y_test_ex[i]))\n",
    "\n",
    "    y_test_GAVQE.append(res)\n",
    "    error_test_GAVQE.append(abs(res-y_test_ex[i]))\n",
    "    \n",
    "t1_GAVQE = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### VQE with meta-VQE initialization (opt-meta-VQE)\n",
    "\n",
    "Another potential advantage of running a meta-VQE algorithm is using the optimized variables of the later as initialization for the former. It is well known the barren plateau problem of PQC: if the minimization algorithm starts at a random point, both the gradient and variance of optimization variables tend to zero exponentially, leading a local minima result. There are several proposals to circunvent this problem. Ours is to start in a region close to the ground state solution by using the result of a meta-VQE as initialization point for a standard VQE (opt-meta-VQE).\n",
    "\n",
    "Since the standard VQE do not contain encoding layers (with more optimization variables than the processing layers), we need to use a meta-VQE-like circuit to run the opt-meta-VQE (which uses the variables of the meta-VQE as starting point for the optimization). The use of encoding layes for a standard VQE had no sense, since the variables are initialized at random and, therefore, the two variables per rotational gate of the encoding layers can be redifined as a single variable. However, in this case, we are not starting at random but in a particular point suggested by the meta-VQE solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0_optmetaVQE = time.time()\n",
    "\n",
    "y_optmetaVQE = []\n",
    "error_optmetaVQE = []\n",
    "\n",
    "for i in range(n_train):\n",
    "    ham = HXXZ(n_qub, arg_train[i], lamXXZ)\n",
    "\n",
    "    # VQE with encoding layers (a meta-VQE PQC)\n",
    "    total_U = ans_enc(n_qub, n_lay[0], arg_train[i], \"th_en\") + ans_var(n_qub, n_lay[1], 0, \"th\")\n",
    "\n",
    "    Obj = tq.ExpectationValue(H=ham, U=total_U)\n",
    "\n",
    "    # Initial values are the angles obtained after running the meta-VQE\n",
    "    initial_values = metaVQE.angles\n",
    "\n",
    "    optmetaVQE = tq.minimize(objective=Obj, adaptive = True, lr=lr, method_options=mthd_opt, method=methods, gradient=grad_methods,\n",
    "                                           samples=None,\n",
    "                                           initial_values=initial_values, backend=backend, noise=None, device=None,\n",
    "                                           silent=True)\n",
    "    res = optmetaVQE.energy\n",
    "    # res_ex: Exact results same as meta-VQE training points\n",
    "\n",
    "    # print(arg_train[i], y_train_ex[i], res, abs(res - y_train_ex[i]))\n",
    "\n",
    "    y_optmetaVQE.append(res)\n",
    "    error_optmetaVQE.append(abs(res - y_train_ex[i]))\n",
    "    \n",
    "t1_optmetaVQE = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### VQE with GA-VQE initialization (opt-GA-VQE)\n",
    "\n",
    "Finally, last test consist on using the result of the GA-VQE as the starting point for a standard VQE, similarly to the previous test. This time, the number of variables and the PQC is the same as the standard VQE (with no encoding layers)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0_optGAVQE = time.time()\n",
    "\n",
    "y_optGAVQE = []\n",
    "error_optGAVQE= []\n",
    "\n",
    "for i in range(n_train):\n",
    "    ham = HXXZ(n_qub, arg_train[i], lamXXZ)\n",
    "        \n",
    "    # VQE number circuit layers same as meta-VQE layers (encoding + variational parts)\n",
    "    total_U = ans_var(n_qub, n_lay[0] + n_lay[1], 0, \"th\")\n",
    "\n",
    "    Obj = tq.ExpectationValue(H=ham, U=total_U)\n",
    "\n",
    "    # Initial values are the ones obtained in the GA-VQE optimization\n",
    "    initial_values = GA_VQE.angles\n",
    "\n",
    "    optGAVQE = tq.minimize(objective=Obj, adaptive = True, lr=lr, method_options=mthd_opt, method=methods, gradient=grad_methods,\n",
    "                                                 samples=None,\n",
    "                                                 initial_values=initial_values, backend=backend, noise=None,\n",
    "                                                 device=None, silent=True)\n",
    "    res = optGAVQE.energy\n",
    "    # res_ex: Exact results same as meta-VQE training points\n",
    "\n",
    "    #print(arg_train[i], y_train_ex[i], res, abs(res - y_train_ex[i]))\n",
    "\n",
    "    y_optGAVQE.append(res)\n",
    "    error_optGAVQE.append(abs(res - y_train_ex[i]))\n",
    "    \n",
    "t1_optGAVQE = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x_test, y_test_ex, color=\"black\", ls=\":\", label=\"exact\")\n",
    "plt.plot(x_test, y_test_metaVQE,  ls=\"-\", label=\"metaVQE test\")\n",
    "plt.scatter(x_train, y_train_metaVQE, color=\"red\", marker=\"x\", label=\"metaVQE train\")\n",
    "plt.plot(x_train, y_stVQE,  marker=\".\", label=\"VQE\")\n",
    "plt.plot(x_test, y_test_GAVQE,  ls=\"--\", label=\"GA-VQE\")\n",
    "plt.plot(x_train, y_optmetaVQE, marker=\"*\", label=\"opt-meta-VQE\")\n",
    "plt.plot(x_train, y_optGAVQE, marker=\"v\", label=\"opt-GA-VQE\")\n",
    "\n",
    "plt.xlabel(\"Δ\")\n",
    "plt.ylabel('g.s. energy')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "# Absolute error\n",
    "\n",
    "plt.plot(x_test, error_test_metaVQE, ls=\"-\", label=\"metaVQE test\")\n",
    "plt.scatter(x_train, error_train_metaVQE, marker=\"x\", color=\"red\", label=\"metaVQE train\")\n",
    "plt.plot(x_train, error_stVQE, marker=\".\", label=\"VQE\")\n",
    "plt.plot(x_test, error_test_GAVQE, ls=\"--\", label=\"GA-VQE\")\n",
    "plt.plot(x_train, error_optmetaVQE, marker=\"*\", label=\"opt-meta-VQE\")\n",
    "plt.plot(x_train, error_optGAVQE, marker=\"v\", label=\"opt-GA-VQE\")\n",
    "\n",
    "plt.xlabel(\"Δ\")\n",
    "plt.ylabel(\"Absolute error\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}