notebooks/5.3_quantum_approximate_optimazation_algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5-3. Quantum Approximate Optimazation Algorithm (QAOA): 量子近似最適化アルゴリズム"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 概要\n",
    "\n",
    "この節では、NISQアルゴリズムの一つと考えられる Quantum Approximate Optimazation Algorithm (QAOA; 量子近似最適化アルゴリズム)を学ぶ。QAOAは、量子アニーリングと同様に、組み合わせ最適化問題の解を求めるためのアルゴリズムである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 問題設定\n",
    "QAOAでは、$z = z_{1}z_{2}\\cdots z_{n} \\: (z_i =0,1)$ という$n$桁のビット列$z$に関して、コスト関数$C(z) = \\sum_\\alpha C_\\alpha(z)$が最小になるような$z$を探す問題を考える。$C_\\alpha(z)$はビット列$z$を引数にとる何らかの関数で、ここでは特に、イジングモデル的な$C_\\alpha(z) = z_i\\cdot z_j$といった項を考えれば良い。\n",
    "\n",
    "この最小化問題を解くために、$n$**ビットの量子系を用いる**。そして、$\\beta = (\\beta^{(1)}, \\cdots \\beta^{(p)}), \\gamma = (\\gamma^{(1)}, \\cdots \\gamma^{(p)})$ をパラメータとして、次のような量子状態を考える。\n",
    "\n",
    "\\begin{align}\n",
    "&|s\\rangle = |+\\rangle^{\\otimes n} = \\frac{1}{2^{n/2}} \\sum_{z=0}^{2^{n}-1} |z\\rangle, \\\\\n",
    "&|\\beta, \\gamma \\rangle = U_X(\\beta^{(p)}) U_C(\\gamma^{(p)}) \\cdots U_X(\\beta^{(1)}) U_C(\\gamma^{(1)}) |s\\rangle. \n",
    "\\end{align}\n",
    "\n",
    "ここで $|+\\rangle=\\frac{1}{\\sqrt{2}}(|0\\rangle+|1\\rangle)$は$X$演算子の固有状態$X|+\\rangle=|+\\rangle$であり、 $U_C(\\gamma), U_X(\\beta)$ は次のように定義される。\n",
    "\n",
    "$$\n",
    "U_C(\\gamma^{(i)}) = e^{-i\\gamma^{(i)} C(Z)} = \\prod_{\\alpha} e^{-i\\gamma^{(i)} C_{\\alpha}(Z)}, \\\\\n",
    "U_X(\\beta^{(i)}) = e^{-i\\beta^{(i)} \\sum_{j=1}^n X_j} = \\prod_{j =1}^n e^{-i\\beta^{(i)} X_j}.\n",
    "$$\n",
    "\n",
    "状態$|\\beta, \\gamma \\rangle$やこれらの演算子の意味を説明するには量子アニーリングの知識が必要になる。とりあえず、QAOAを使うだけならこういうものだと受け入れて使ってしまえば良い。    \n",
    "（なお、$C(Z)$というのはビット列を引数にとる関数$C(z)$の入力にパウリ演算子$Z_1\\cdots Z_n$と代入したものであることに注意。）\n",
    "\n",
    "そして、$F(\\beta, \\gamma) = \\langle{\\bf \\gamma, \\,\\beta}|C(Z)|{\\bf \\gamma, \\,\\beta}\\rangle$ を最小にするような$\\beta,\\gamma$を探索することで、元々の最適化問題の答えを探そうとするのが、QAOAアルゴリズムである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QAOAアルゴリズムの手順\n",
    "具体的なQAOAアルゴリズムの手順は以下の通りである。\n",
    "\n",
    "1. 量子コンピュータ上で重ね合わせ状態$|s\\rangle = |+\\rangle^{\\otimes n}$を作る。\n",
    "2. パラメータ$\\beta, \\gamma$に応じて、量子状態に$U_C(\\gamma^{(i)}),U_X(\\beta^{(i)})$をかけていき、状態$|\\beta, \\gamma \\rangle$を得る。\n",
    "3. 量子コンピュータを用いて $\\langle \\beta, \\gamma |C(Z)|\\beta, \\gamma \\rangle$ を測定する。\n",
    "4. 古典コンピュータで、$\\langle \\beta, \\gamma |C(Z)|\\beta, \\gamma \\rangle$ がより小さくなるようにパラメータ $\\beta, \\gamma$ をアップデートする。\n",
    "5. 1〜4を繰り返し、最適な $\\beta^*, \\gamma^*$ を得る。\n",
    "6. 状態 $|\\beta^*, \\gamma^* \\rangle$ に対して、$z$方向の射影測定を複数回実行し、**得られた（良さそうな）測定結果** $z_1\\cdots z_n$ **を元々の最適化問題の解として採用**する。（注：測定結果 $z_1\\cdots z_n$ は古典ビット）\n",
    "\n",
    "少々ややこしいので、具体例を実装しながら確認していこう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 実装：Maxcut問題をQAOAで解く\n",
    "\n",
    "ここでは具体例としてMaxcut問題という問題を取り上げる。   \n",
    "[Maxcut問題](https://ja.wikipedia.org/wiki/カット_(グラフ理論))は、$n$個の頂点を持つグラフ（例えば下図）を２つに分割する時に、分割される辺の数の最大値を求める問題である。\n",
    "\n",
    "![maxcut-example](figs/5/maxcut_example.png)\n",
    "(図の出典：Wikipedia [カット_(グラフ理論)](https://ja.wikipedia.org/wiki/カット_(グラフ理論)))\n",
    "\n",
    "この問題をQAOAで扱えるような最適化問題に帰着させるには、以下のようにする。\n",
    "頂点を２つのグループに分けた時、片方のグループに属する頂点に+1、もう一方のグループに-1を付与するとすれば、コスト関数\n",
    "\n",
    "$$\n",
    "C(z) = -\\frac{1}{2} \\sum_{\\text{辺で繋がっている頂点}i,j} ( 1 - z_i z_j)\n",
    "$$\n",
    "\n",
    "は (グループ分けによって分割される辺の数) $\\times (-1)$ を表す。\n",
    "ゆえに、$C(z)$を最小化するようなビット列$z=z_1\\cdots z_n$を見つければ、分割する辺の数を最大化するような頂点の分け方を見つけたことになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下では、長方形(頂点が４つの図形)のmaxcut問題を解いてみよう。\n",
    "![maxcut](figs/5/maxcut_4nodes.png)\n",
    "\n",
    "この場合、$C(Z)$は\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "C(Z) &=& -\\frac{1}{2}(1-Z_{0}Z_{1})-\\frac{1}{2}(1-Z_{1}Z_{2})-\\frac{1}{2}(1-Z_{2}Z_{3})-\\frac{1}{2}(1-Z_{3}Z_{1})\\\\\n",
    "&=&\\frac{1}{2}(Z_{0}Z_{1}+Z_{1}Z_{2}+Z_{2}Z_{3}+Z_{3}Z_{1}) - 2\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "となる。第二項は定数だから、以下では\n",
    "\n",
    "$$\n",
    "C(Z) = \\frac{1}{2}(Z_{0}Z_{1}+Z_{1}Z_{2}+Z_{2}Z_{3}+Z_{3}Z_{1}) \n",
    "$$\n",
    "\n",
    "とおく。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $p=1$の場合\n",
    "まずは$p=1$の場合の実装をやってみよう。この時、$|\\beta, \\gamma \\rangle = U_X(\\beta^{(1)}) U_C(\\gamma^{(1)}) |s\\rangle$ である。\n",
    "\n",
    "$U_C(\\gamma^{(1)}) = \\prod_{i=0}^3 e^{-i\\gamma^{(1)} Z_i Z_{i+1} }$ を実装するには、4-2節でも用いた関係\n",
    "\n",
    "$$\n",
    "e^{-i \\delta Z_i Z_{i+1}}　= \\operatorname{CNOT}_{i,i+1} \\cdot e^{-i\\delta Z_{i+1}} \\cdot \\operatorname{CNOT}_{i,i+1}.\n",
    "$$\n",
    "\n",
    "を使えば良い。（行列に直して計算すると、合っていることがわかる。）\n",
    "\n",
    "以上を踏まえて、$|\\beta, \\gamma \\rangle$ を構成して $\\langle \\beta, \\gamma | C(Z) |\\beta, \\gamma \\rangle$ を測定し、それを最小化する工程を実装するとこのようになる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Google Colaboratoryの場合・Qulacsがインストールされていないlocal環境の場合のみ実行してください\n",
    "!pip install qulacs\n",
    "\n",
    "## Google Colaboratory / (Linux or Mac)のjupyter notebook 環境の場合にのみ実行してください。\n",
    "## Qulacsのエラーが正常に出力されるようになります。\n",
    "!pip3 install wurlitzer\n",
    "%load_ext wurlitzer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.999999999499185\n",
      "[1.17809152 0.39269362]\n"
     ]
    }
   ],
   "source": [
    "#必要なライブラリをインポートする\n",
    "from qulacs import QuantumState, QuantumCircuit, Observable, PauliOperator\n",
    "from qulacs.gate import H, CNOT, RX, RZ\n",
    "from scipy.optimize import minimize\n",
    "import numpy as np\n",
    "\n",
    "## 頂点の数\n",
    "n = 4 \n",
    "\n",
    "## C(Z)をqulacs.Observableとして定義\n",
    "cost_observable = Observable(n) \n",
    "for i in range(n):\n",
    "    cost_observable.add_operator( PauliOperator(\"Z {:} Z {:}\".format(i, (i+1)%n), 0.5) )\n",
    "\n",
    "# circuit に U_C(gamma) を加える関数\n",
    "def add_U_C(circuit, gamma):\n",
    "    for i in range(n):\n",
    "        j = (i+1) % n\n",
    "        circuit.add_CNOT_gate(i, j)\n",
    "        circuit.add_gate(RZ(j, -2*gamma)) ## qulacsでは RZ(theta)=e^{i*theta/2*Z}\n",
    "        circuit.add_CNOT_gate(i, j)\n",
    "    return circuit\n",
    "    \n",
    "# circuit に U_X(beta) を加える関数\n",
    "def add_U_X(circuit, beta):\n",
    "    for i in range(n):\n",
    "        circuit.add_gate(RX(i, -2*beta))\n",
    "    return circuit\n",
    "\n",
    "# p=1 の |beta, gamma> を作って <beta, gamma| C(Z) |beta, gamma> を返す関数\n",
    "# x = [beta, gamma]\n",
    "def QAOA_output_onelayer(x): \n",
    "    beta, gamma = x\n",
    "\n",
    "    circuit = QuantumCircuit(n)\n",
    "    ## 重ね合わせを作るため、アダマールゲートをかける\n",
    "    for i in range(n):\n",
    "        circuit.add_H_gate(i)\n",
    "    ## U_C, U_Xをかける\n",
    "    circuit =  add_U_C(circuit, gamma)\n",
    "    circuit =  add_U_X(circuit, beta)\n",
    "\n",
    "    ## |beta, gamma>を作る\n",
    "    state = QuantumState(n)\n",
    "    state.set_zero_state()                     \n",
    "    circuit.update_quantum_state(state)\n",
    "    return cost_observable.get_expectation_value(state)                 \n",
    "\n",
    "## 初期値\n",
    "x0 = np.array( [0.1, 0.1 ])\n",
    "\n",
    "## scipy.minimize を用いて最小化\n",
    "result = minimize(QAOA_output_onelayer, x0, options={'maxiter':500}, method='powell')\n",
    "print(result.fun) # 最適化後の値\n",
    "print(result.x) # 最適化後の(beta, gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最小値-1と、その時の$\\beta^{(1)}, \\gamma^{(1)}$の値が得られた。   \n",
    "この最適な状態 $|\\beta, \\gamma\\rangle$ を$z$方向に射影測定した時にどのような値が得られるか、具体的に見てみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.01562503 0.01562568 0.01562568 0.0781236  0.01562568 0.26562503\n",
      " 0.0781236  0.01562568 0.01562568 0.0781236  0.26562503 0.01562568\n",
      " 0.0781236  0.01562568 0.01562568 0.01562503]\n"
     ]
    }
   ],
   "source": [
    "# 最適なbeta, gammaを使って |beta, gamma> をつくる\n",
    "beta_opt, gamma_opt = result.x\n",
    "\n",
    "circuit = QuantumCircuit(n)\n",
    "## 重ね合わせを作るため、アダマールゲートをかける\n",
    "for i in range(n):\n",
    "    circuit.add_H_gate(i)\n",
    "## U_C, U_Xをかける\n",
    "circuit =  add_U_C(circuit, gamma_opt)\n",
    "circuit =  add_U_X(circuit, beta_opt)\n",
    "\n",
    "## |beta, gamma>を作る\n",
    "state = QuantumState(n)\n",
    "state.set_zero_state()                     \n",
    "circuit.update_quantum_state(state)\n",
    "\n",
    "## z方向に観測した時の確率分布を求める. (状態ベクトルの各成分の絶対値の二乗=観測確率)\n",
    "probs = np.abs(state.get_vector())**2\n",
    "print(probs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFACAYAAAD07atFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAF5dJREFUeJzt3Xu0pXdZH/DvYyIBgUBixhgodbh5SVUgDBglskAkhEuJFkRoC2mLjbcsQGGVeGmBLliNIGgRFxgkJVSLN2CRkhQb04BGBZlAyAXEpHRYEAIZyiUBgRLy9I/9xnUyzJlzyd77d+acz2etvfbe7/vu93ne/Z5z5jvvtbo7AAAs1zeNbgAAYCcSwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGOHJ0A+tx3HHH9e7du0e3AQCwpssvv/wz3b1rrekOixC2e/fu7N27d3QbAABrqqqPrWc6uyMBAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABjgs7h0JO8Husy9cSp195zxxKXU4vPj5g+WzJQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABggIWFsKq6T1VdWlUfqqprquq50/AXV9X1VXXF9HjConoAANiqFnkD71uSPL+7319Vd09yeVVdPI37je7+9QXWBgDY0hYWwrr7hiQ3TK9vrqoPJ7n3ouoBABxOlnJMWFXtTvKQJO+dBp1VVVdW1XlVdcwqnzmzqvZW1d79+/cvo00AgKVZeAirqrsleUuS53X3TUlem+T+SR6c2ZayVx7sc919bnfv6e49u3btWnSbAABLtdAQVlXfnFkA+/3ufmuSdPenu/vr3X1rktcnefgiewAA2IoWeXZkJXlDkg9396tWDD9hxWQ/nuTqRfUAALBVLfLsyEckeWaSq6rqimnYLyd5RlU9OEkn2ZfkpxfYAwDAlrTIsyMvS1IHGXXRomoCABwuXDEfAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGCAhYWwqrpPVV1aVR+qqmuq6rnT8GOr6uKqunZ6PmZRPQAAbFWL3BJ2S5Lnd/eJSU5O8vNVdWKSs5Nc0t0PTHLJ9B4AYEdZWAjr7hu6+/3T65uTfDjJvZOcnuT8abLzk/zYonoAANiqlnJMWFXtTvKQJO9Ncnx33zCN+lSS41f5zJlVtbeq9u7fv38ZbQIALM3CQ1hV3S3JW5I8r7tvWjmuuztJH+xz3X1ud+/p7j27du1adJsAAEu10BBWVd+cWQD7/e5+6zT401V1wjT+hCQ3LrIHAICtaJFnR1aSNyT5cHe/asWoC5KcMb0+I8nbF9UDAMBWdeQC5/2IJM9MclVVXTEN++Uk5yT5o6p6dpKPJXnaAnsAANiSFhbCuvuyJLXK6Mcsqi4AwOHAFfMBAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGOHKtCarqzkmelOSHk9wryZeTXJ3kwu6+ZrHtAQBsT4cMYVX1kswC2LuSvDfJjUnunOQ7k5wzBbTnd/eVC+4TAGBbWWtL2N9094tWGfeqqvq2JP94zj0BAGx7hwxh3X3hgcOmrV936u6buvvGzLaOAQCwAWseE7ZSVf1UkqcmOaKq9nb3Ly2mLQCA7e2QZ0dW1ZMPGPSj3X1adz82yRMW1xYAwPa21iUqvq+q3l5VD57eX1lVv1tVr0/izEgAgE1a65iwl1XVtyf5j1VVSf59krsnuYszIgEANm89x4R9KcnzkjwwyblJ9iZ5+SKbAgDY7tY6JuylSd6S5B1JHt3dT05yRZKLqupZS+gPAGBbWuuYsCd196lJHpPkWUnS3RckOTXJMQvuDQBg21prd+TVVXVukrskefdtA7v7liT/eZGNAQBsZ2sdmP8vq+r7knytu/92ST0BAGx7ax0Tdkp3X7VaAKuqo6vqexfTGgDA9rXW7sinVNXLk7wzyeVJ9md2A+8HJHl0ku9I8vyFdggAsA2ttTvyF6rq2CRPSfITSU5I8uUkH07yO9192eJbBADYfta8Tlh3f7aqzuvu1y+jIQCAnWCtS1Tc5tqqekVVfc9CuwEA2CHWG8IelOTvkryhqt5TVWdW1dGH+kBVnVdVN1bV1SuGvbiqrq+qK6aHm4ADADvSukJYd9/c3a/v7h9K8sIkL0pyQ1WdX1UPWOVjb0xy2kGG/0Z3P3h6XLSprgEADnPrCmFVdURVPbmq3pbkN5O8Msn9kvz3JAcNUt3950k+O69GAQC2k/XcwDtJrk1yaZJXdPdfrRj+J1X1yA3WPGu67+TeJM/v7s9t8PMAAIe99R4T9qzufvbKAFZVj0iS7n7OBuq9Nsn9kzw4yQ2ZbVE7qOm4s71VtXf//v0bKAEAsPWtN4S9+iDDfmujxbr709399e6+Ncnrkzz8ENOe2917unvPrl27NloKAGBLO+TuyKr6wSQ/lGRXVf3iilFHJzlio8Wq6oTuvmF6++NJrj7U9AAA29Vax4TdKcndpunuvmL4TUmeeqgPVtWbkzwqyXFV9YnMzqh8VFU9OEkn2ZfkpzfVNQDAYW6t2xa9O8m7q+qN3f2xjcy4u59xkMFv2Mg8AAC2q7V2R/5mdz8vyWuqqg8c391PXlhnAADb2Fq7I//r9Pzri24EAGAnWWt35OXT87uX0w4AwM6w1u7IqzI7iP6guvv7594RAMAOsNbuyCctpQsAgB1mrd2RGzojEgCA9TnkFfOr6rLp+eaquunA5+W0CACw/ay1JeyU6fnuh5oOAICNWeuYsH9QVSclOSWzA/Uv6+4PLKwrAIBtbl038K6q/5Dk/CTfmuS4JG+sql9dZGMAANvZereE/YskD+ruryRJVZ2T5IokL11UYwAA29m6toQl+WSSO694f1SS6+ffDgDAzrDWxVp/K7NjwL6Q5Jqqunh6/9gkf7P49gAAtqe1dkfunZ4vT/K2FcPftZBuAAB2iLUuUXH+shoBANhJ1nVgflU9MMl/SnJiVhwb1t33W1BfAADb2noPzP8vSV6b5JYkj07ypiS/t6imAAC2u/WGsLt09yVJqrs/1t0vTvLExbUFALC9rfc6YV+tqm9Kcm1VnZXZ5Snutri2AAC2t/VuCXtukm9J8pwkD03yzCRnLKopAIDtbl1bwrr7fUkybQ17TnffvNCuAAC2ufXeO3JPVV2V5MokV1XVB6vqoYttDQBg+1rvMWHnJfm57v6LJKmqUzI7Y/L7F9UYAMB2tt5jwr5+WwBLku6+LLPLVQAAsAlr3TvypOnlu6vqd5K8ObN7R/5k3LoIAGDT1tod+coD3r9oxeuecy8AADvGWveOfPSyGgEA2EnWe3bkParqVVW1d3q8sqrusejmAAC2q/UemH9ekpuTPG163JTZ2ZEAAGzCei9Rcf/ufsqK9y+pqisW0RAAwE6w3i1hX56uDZYkqapHJPnyYloCANj+1rsl7GeSvGnFcWCfi3tHAgBs2pohbLpf5Hd194Oq6ugk6e6bFt4ZAMA2tubuyO6+Ncm/m17fJIABANxx6z0m7M+q6gVVdZ+qOva2x0I7AwDYxtZ7TNhPZnaF/J87YPj95tsOAMDOsN4QdmJmAeyUzMLYXyR53aKaAgDY7tYbws7P7AKtr57e//Np2NMW0RQAwHa33hD2vd194or3l1bVhxbREADATrDeA/PfX1Un3/amqn4gyd7FtAQAsP2tN4Q9NMlfVdW+qtqX5K+TPKyqrqqqKw/2gao6r6purKqrVww7tqourqprp+dj7vASAAAchta7O/K0Tcz7jUlek+RNK4adneSS7j6nqs6e3r9wE/MGADisrSuEdffHNjrj7v7zqtp9wODTkzxqen1+kndFCAMAdqD17o6cl+O7+4bp9aeSHL/ahFV1ZlXtraq9+/fvX053AABLsuwQ9g+6uzO75thq48/t7j3dvWfXrl1L7AwAYPGWHcI+XVUnJMn0fOOS6wMAbAnLDmEXJDljen1GkrcvuT4AwJawsBBWVW/O7FIW31VVn6iqZyc5J8ljq+raJD86vQcA2HHWe4mKDevuZ6wy6jGLqgkAcLgYdmA+AMBOJoQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMcOToBuA2u8++cOE19p3zxIXXOFwt4/tPrIPV+P7H8zeIZbMlDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGAAIQwAYAAhDABgACEMAGCAI0cUrap9SW5O8vUkt3T3nhF9AACMMiSETR7d3Z8ZWB8AYBi7IwEABhgVwjrJ/6yqy6vqzINNUFVnVtXeqtq7f//+JbcHALBYo0LYKd19UpLHJ/n5qnrkgRN097ndvae79+zatWv5HQIALNCQENbd10/PNyZ5W5KHj+gDAGCUpYewqrprVd39ttdJTk1y9bL7AAAYacTZkccneVtV3Vb/v3X3Owf0AQAwzNJDWHd/NMmDll0XAGArcYkKAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABhDAAgAGEMACAAYQwAIABjhzdwFax++wLl1Jn3zlP3JL1l9XDoeqPttOXf7TRvwOj6+90vv/xf4N2Qv21elg2W8IAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYYEsKq6rSq+khVXVdVZ4/oAQBgpKWHsKo6IslvJ3l8khOTPKOqTlx2HwAAI43YEvbwJNd190e7+/8l+YMkpw/oAwBgmBEh7N5JPr7i/SemYQAAO0Z193ILVj01yWnd/VPT+2cm+YHuPuuA6c5Mcub09ruSfGSpja7PcUk+o/6Orb8VelBfffV3bv2t0MNOr7+a7+juXWtNdOQyOjnA9Unus+L9P5qG3U53n5vk3GU1tRlVtbe796i/M+tvhR7UV1/9nVt/K/Sw0+vfUSN2R74vyQOr6r5VdackT09ywYA+AACGWfqWsO6+parOSvKnSY5Icl53X7PsPgAARhqxOzLdfVGSi0bUnrPRu0vVH290D+qrr/7OrZ+M72Gn179Dln5gPgAAblsEADCEEAYAMIAQdoCD3ddyOpPzvdOwP5zO6kxVHTW9v24av3vFfH5pGv6RqnrcMutX1bdW1aVV9cWqes0Cl/+RVfX+qrpluv7byvmcUVXXTo8zBtR/Z1V9vqresZHlP0QPZ03vu6qOWzFtVdWrp3FXVtVJC/oONlN/U9/BBut/d1X9dVV9tapesNZ8llz/vKq6saqu3uDyf8PnqurYqrp4WpcXV9Ux0/BFrP951d/s+t9I/UWs/3nVn+f6/4mquqaqbq2qPQdMf9C/9XNe/s3UX/jy1yH+ramqh1bVVVNvr66qWnL9l1XVx6vqixtZ/qXrbo/pkdnZmv87yf2S3CnJBzO7v+UfJXn6NM3rkvzs9Prnkrxuev30JH84vT5x+uxRSe47zfOIJda/a5JTkvxMktcscPl3J/n+JG9K8tQV8zk2yUen52Om18csq/407jFJ/mmSd8zpZ+AhU719SY5bMf0TkvyPJJXk5CTvXdB3sKH6m/0ONlH/25I8LMnLkrxgrfksq/407pFJTkpy9QZ/Br7hc0lenuTs6fXZSX5tEet/XvXv4O/ARurPdf3Pq/4C1v/3ZHbR8Hcl2bNi+EH/1i9g+TdUf4nLv+q/NUn+ZvqZrOln9PFLrn9ykhOSfHEjy7/shy1ht7fafS1/JMmfTNOcn+THptenT+8zjX/MlPZPT/IH3f3V7v4/Sa6b5r2U+t39pe6+LMlXNrb4G6vf3fu6+8oktx4wn8clubi7P9vdn0tycZLTllg/3X1JkpvXt9hr99DdH+jufQeZ/vQkb+qZ9yS5Z1WdkDl/B5uov9nvYEP1u/vG7n5fkq+tZz5LrJ/u/vMkn11HzfV8buXv2oG/g/Nc//Oqv+nfgY3UX8D6n1f9ua7/7v5wdx/sri2r/a2f6/Jvov5Sln+1f2umn8Gju/s93d2Z/Uf5xw78/KLqT+Pe0903rKfmSELY7a12X8vPd/ctBwy73fTT+C8k+dZDzGdZ9Tdro/U3Op9l1b8jNtr7atPP+zuY1/Tzrj/v+WzVe8sev+IP+qeSHD+9nvf6n1f9eVut/mqWtfyjLWv9b7T+aPfOrJfbbJW+thwhDGADpv/ZD7u2j/pj68M8CWG3t9p9Le9ZVUceMOx200/j75Hk/x5iPsuqv1kbrb/R+Syr/h2x0d5Xm37e38G8pp93/XnPZ97LMy+fvm033/R84zR83ut/XvXnbbX6q1nW8o+2rPW/0fqjXZ9ZL7fZKn1tOULY7a12X8tLk9x29t0ZSd4+vb5gep9p/P+a/pd2QZKn1+zsxfsmeWBmBykuq/5mbbT+av40yalVdcx0FtOp07Bl1b8jNnpv0wuSPKtmTk7yhWm3yby/g43W36x53dt1s/PZqveWXfm7duDv4DzX/7zqz9tq9Vcz7/W40frLstrf+mX9HG/235qFmn4Gb6qqk6fjpJ+VrbPOtpbeAmcHbKVHZmcb/V1mZ5n8yjTsfpn9YF+X5I+THDUNv/P0/rpp/P1WzOdXpnl8JOs8K2TO9fdldnDjFzPbH7/eM3M2Uv9h07y/lNkWuGtWzOffTNNfl+RfL2j5D1X/L5LsT/LlaZrH3cEenjPN55Ykn0zyu9PwSvLb07RX5fZn7szzO9hM/U19Bxus/+3T8JuSfH56ffRq81ly/TcnuSGzg7Y/keTZ66z/DZ/L7FjLS5Jcm+TPkhy7wPU/r/qbXf8bqb+I9T+v+vNc/z8+vf5qkk8n+dMV0x/0b/2cl38z9Ze1/PtykH9rkuxJcvXU22sy3aFnifVfPr2/dXp+8XrXwTIfblsEADCA3ZEAAAMIYQAAAwhhAAADCGEAAAMIYQAAAwhhwLZTVc+rqm+Z13QAi+ASFcC2U1X7Mrtm1mfmMR3AItgSBhzWququVXVhVX2wqq6uqhcluVeSS6vq0mma11bV3qq6pqpeMg17zkGmO7Wq/rqq3l9Vf1xVd5uGn1NVH6qqK6vq18csKbDd2BIGHNaq6ilJTuvufzu9v0eSD2bFFq6qOra7P1tVR2R25fXndPeVK7eEVdVxSd6a2VXHv1RVL0xyVGZXpP+rJN/d3V1V9+zuzy99QYFtx5Yw4HB3VZLHVtWvVdUPd/cXDjLN06rq/Uk+kOSfJDnxINOcPA3/y6q6IrN7FH5Hki8k+UqSN1TVP0vy94tYCGDnOXJ0AwB3RHf/XVWdlNm9+l5aVZesHD/d2PgFSR7W3Z+rqjdmdt/VA1WSi7v7Gd8wourhSR6T2Y3kz0ryI/NdCmAnsiUMOKxV1b2S/H13/16SVyQ5KcnNSe4+TXJ0Zjd5/0JVHZ/k8Ss+vnK69yR5RFU9YJrvXavqO6fjwu7R3Rcl+YUkD1r0MgE7gy1hwOHu+5K8oqpuTfK1JD+b5AeTvLOqPtndj66qDyT52yQfT/KXKz577gHT/askb66qo6bxv5pZUHt7Vd05s61lv7iUpQK2PQfmAwAMYHckAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAA/x9CODL0v4JenQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# プロットする\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "## z方向に射影測定した時に得られる可能性があるビット列\n",
    "z_basis = [format(i,\"b\").zfill(n) for i in range(probs.size)]\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.xlabel(\"states\")\n",
    "plt.ylabel(\"probability(%)\")\n",
    "plt.bar(z_basis, probs*100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "つまり、$z$方向の射影測定を行うと、`0101` か `1010` が測定される確率が高いことが分かった。これらのビット列は頂点1と頂点3、頂点2と頂点4が同じグループになるということを意味するから、以下のような分割を表している。\n",
    "\n",
    "![maxcut-p-1](figs/5/maxcut_result_p1.png)\n",
    "\n",
    "この時、図形を分割する曲線が横切る辺の数は4本であり、この図形を分割する時に通る辺の数の最大値である。\n",
    "\n",
    "よって、最適化した$|\\beta, \\gamma\\rangle$に射影測定を行い、ある程度多数の測定結果を行って測定確率が高いビット列を採用すれば、もともと解きたかった最適化問題$C(z)$の解が得られた、ということになる。   \n",
    "一応、これでめでたしめでたし、と言えるのだが、最適化したコスト関数$\\langle \\beta, \\gamma | C(Z) |\\beta, \\gamma \\rangle$ の値は−1だったことを思い出してほしい。$\\langle 0101 | C(Z) |0101 \\rangle = \\langle 1010 | C(Z) |1010 \\rangle = -2$なので、コスト関数については正しい値が得られていない！ これはvariationalな状態$|\\beta, \\gamma \\rangle$ が十分な表現能力を持たず、真の解$|0101\\rangle, |1010\\rangle$ を表現できなかったことに由来する。\n",
    "\n",
    "そこで、回路をより複雑にした$p=2$の場合に結果がどう変わるか見てみよう。\n",
    "\n",
    "※ ちなみに、$|\\beta, \\gamma\\rangle$に100回の測定を行ってビット列$z$を100通り得て、それぞれについて$C(z)$を古典コンピュータで計算してみて最も良かったものを採用する、といった戦略を用いれば、このような問題は生じないかもしれない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $p=2$の場合"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.9999996223652372\n",
      "[1.01150163 1.11872823 0.45213228 0.55937865]\n",
      "[1.10913271e-15 9.96406007e-09 9.96406007e-09 2.72762241e-08\n",
      " 9.96406007e-09 4.99999906e-01 2.72762241e-08 9.96406007e-09\n",
      " 9.96406007e-09 2.72762241e-08 4.99999906e-01 9.96406007e-09\n",
      " 2.72762241e-08 9.96406007e-09 9.96406007e-09 1.10913271e-15]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFACAYAAAD07atFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAF+NJREFUeJzt3X2wbXdZH/DvYyKgQCAhMUYQAwXFKPJ2oVHQGYhgVErSgghjJba0GWUYRGUkVlu0o9MAiogwaoCUS215EWSSAsXGGFCUtxte8sJbIg0jMZAohAREauTpH3vdenK55559zt17/87L5zNzZu/1+jxr733O/d611l6rujsAAKzW14xuAABgLxLCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAY4dnQD8zjxxBP71FNPHd0GAMCGLr/88r/p7pM2mm9HhLBTTz01Bw4cGN0GAMCGquqT88zncCQAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEu9WGtVXZfk1iT/mOS27t5XVSckeV2SU5Ncl+TJ3f25ZfYBALDdrGJP2KO7+8HdvW8aPi/Jpd19/ySXTsMAAHvKiMORZyXZPz3fn+TsAT0AAAy17HtHdpL/XVWd5Pe6+4IkJ3f3DdP0Tyc5+XALVtW5Sc5Nknvf+95LbhPGO/W8tyy9xnXn//DSa7AzreLzl/gMwlrLDmGP6u7rq+obklxSVR9dO7G7ewpoX2UKbBckyb59+w47DwDATrXUw5Hdff30eGOSNyV5RJLPVNUpSTI93rjMHgAAtqOlhbCqunNV3fXg8ySPS3JVkouTnDPNdk6Si5bVAwDAdrXMw5EnJ3lTVR2s8z+6+21V9b4kr6+qpyf5ZJInL7EHAIBtaWkhrLs/keRBhxn/t0nOWFZdAICdwBXzAQAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABlh6CKuqY6rqA1X15mn4PlX1nqq6tqpeV1V3WHYPAADbzSr2hP10ko+sGX5+kt/s7vsl+VySp6+gBwCAbWWpIayq7pXkh5O8YhquJI9J8oZplv1Jzl5mDwAA29Gy94S9OMnPJ/nKNHyPJDd3923T8KeS3HPJPQAAbDtLC2FV9fgkN3b35Vtc/tyqOlBVB2666aYFdwcAMNYy94Q9MskTquq6JK/N7DDkbyW5e1UdO81zryTXH27h7r6gu/d1976TTjppiW0CAKze0kJYd/9Cd9+ru09N8pQkf9LdP5bksiRPmmY7J8lFy+oBAGC7GnGdsOcm+dmqujazc8ReOaAHAIChjt14lqPX3W9P8vbp+SeSPGIVdQEAtitXzAcAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGEAIAwAYQAgDABhACAMAGODYjWaoqjsleXyS703yTUm+lOSqJG/p7quX2x4AwO50xBBWVb+SWQB7e5L3JLkxyZ2SfGuS86eA9nPdfcWS+wQA2FU22hP23u5+3jrTXlRV35Dk3gvuCQBg1ztiCOvutxw6btr7dYfuvqW7b8xs7xgAAJuw4Tlha1XVv0vypCTHVNWB7v6F5bQFALC7HfHbkVX1hENGfX93n9ndj03yQ8trCwBgd9voEhUPrKqLqurB0/AVVfWKqnp5Et+MBADYoo3OCfu1qvrGJP+5qirJf0xy1yRf5xuRAABbN885YV9M8uwk909yQZIDSV6wzKYAAHa7jc4J+9Ukb0zy5iSP7u4nJPlgkrdW1dM2WPZOVfXeqvpQVV09XXMsVXWfqnpPVV1bVa+rqjssaFsAAHaMjc4Je3x3Py7JGUmeliTdfXGSxyU5foNlv5zkMd39oCQPTnJmVZ2e5PlJfrO775fkc0mefhT9AwDsSBuFsKuq6oIkr07yjoMju/u27v6tIy3YM1+YBr92+ukkj0nyhmn8/iRnb6VxAICdbKMT8/91VT0wyT9090c3u/KqOibJ5Unul+RlSf4yyc3dfds0y6eS3HOdZc9Ncm6S3PveLsoPAOwuG50T9qjuvnK9AFZVx1XVd663fHf/Y3c/OMm9kjwiyQPmbay7L+jufd2976STTpp3MQCAHWGjb0c+sapekORtme3RuimzG3jfL8mjk3xLkp/bqEh331xVlyX57iR3r6pjp71h90py/VH0DwCwI210OPJnquqEJE9M8iNJTknypSQfSfJ73f3O9ZatqpMyO4x5c1V9XZLHZnZS/mWZ3frotUnOSXLRIjYEAGAn2fA6Yd392aq6sLtfvsl1n5Jk/3Re2NckeX13v7mqPpzktdPlLz6Q5JWb7hoAYIeb9wbe11TVG5Nc2N0fmWeB6Yr6DznM+E9kdn4YAMCetdElKg56UJKPJ3llVb27qs6tquOW2BcAwK42Vwjr7lu7++Xd/T1JnpvkeUluqKr9VXW/pXYIALALzRXCquqYqnpCVb0pyYuT/EaS+yb5n0neusT+AAB2pbnPCcvsW40v7O6/WDP+DVX1fYtvCwBgd5s3hD3t0MtRVNUju/vPu/tZS+gLAGBXm/fE/JccZtxvL7IRAIC95Ih7wqrqu5N8T5KTqupn10w6Lskxy2wMAGA32+hw5B2S3GWa765rxt+S2VXvAQDYgo1uW/SOJO+oqld19ydX1BMAwK630eHIF3f3s5O8tKr60Ond/YSldQYAsIttdDjyv02Pv77sRgAA9pKNDkdePj2+YzXtAADsDRsdjrwyyVcdhjyou79r4R0BAOwBGx2OfPxKugAA2GM2OhzpG5EAAEtwxCvmV9U7p8dbq+qWQx9X0yIAwO6z0Z6wR02Pdz3SfAAAbM68N/BOVT00yaMyO1H/nd39gaV1BQCwy811A++q+k9J9ie5R5ITk7yqqn5pmY0BAOxm8+4J+7EkD+ruv0+Sqjo/yQeT/OqyGgMA2M3m2hOW5K+T3GnN8B2TXL/4dgAA9oaNLtb625mdA/b5JFdX1SXT8GOTvHf57QEA7E4bHY48MD1enuRNa8a/fSndAADsERtdomL/qhoBANhL5joxv6run+S/JDkta84N6+77LqkvAIBdbd4T8/9rkt9JcluSRyd5dZLfX1ZTAAC73bwh7Ou6+9Ik1d2f7O5fTvLDy2sLAGB3m/c6YV+uqq9Jck1VPTOzy1PcZXltAQDsbvPuCfvpJF+f5FlJHpbkx5Ocs6ymAAB2u7n2hHX3+5Jk2hv2rO6+daldAQDscvPeO3JfVV2Z5IokV1bVh6rqYcttDQBg95r3nLALkzyju/8sSarqUZl9Y/K7ltUYAMBuNu85Yf94MIAlSXe/M7PLVQAAsAUb3TvyodPTd1TV7yV5TWb3jvzRuHURAMCWbXQ48jcOGX7emue94F4AAPaMje4d+ehVNQIAsJfM++3Iu1XVi6rqwPTzG1V1t2U3BwCwW817Yv6FSW5N8uTp55bMvh0JAMAWzHuJin/W3U9cM/wrVfXBZTQEALAXzLsn7EvTtcGSJFX1yCRfWk5LAAC737x7wn4yyavXnAf2ubh3JADAlm0Ywqb7RX5bdz+oqo5Lku6+ZemdAQDsYhsejuzuryT5+en5LQIYAMDRm/ecsD+uqudU1TdX1QkHf460wDTvZVX14aq6uqp+ehp/QlVdUlXXTI/HH/VWAADsMPOeE/ajmV0h/xmHjL/vEZa5LcnPdff7q+quSS6vqkuS/ESSS7v7/Ko6L8l5SZ67ubYBAHa2efeEnZbkZUk+lOSDSX47yXccaYHuvqG73z89vzXJR5LcM8lZSfZPs+1Pcvbm2wYA2NnmDWH7k3x7kpdkFsBOyz8FqQ1V1alJHpLkPUlO7u4bpkmfTnLyOsuce/AK/TfddNO8pQAAdoR5D0d+Z3eftmb4sqr68DwLVtVdkrwxybO7+5aq+v/Turur6rA3Au/uC5JckCT79u1zs3AAYFeZd0/Y+6vq9IMDVfXPkxzYaKGq+trMAth/7+4/nEZ/pqpOmaafkuTGzbUMALDzzRvCHpbkL6rquqq6Lsm7kjy8qq6sqisOt0DNdnm9MslHuvtFayZdnH+60Os5SS7aUucAADvYvIcjz9zCuh+Z5MeTXLnmPpP/Icn5SV5fVU9P8snMbggOALCnzBXCuvuTm11xd78zSa0z+YzNrg8AYDeZ93AkAAALJIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMIIQBAAwghAEADCCEAQAMsLQQVlUXVtWNVXXVmnEnVNUlVXXN9Hj8suoDAGxny9wT9qokZx4y7rwkl3b3/ZNcOg0DAOw5Swth3f2nST57yOizkuyfnu9Pcvay6gMAbGerPifs5O6+YXr+6SQnrzdjVZ1bVQeq6sBNN920mu4AAFZk2In53d1J+gjTL+jufd2976STTlphZwAAy7fqEPaZqjolSabHG1dcHwBgW1h1CLs4yTnT83OSXLTi+gAA28IyL1HxmiTvSvJtVfWpqnp6kvOTPLaqrkny/dMwAMCec+yyVtzdT11n0hnLqgkAsFO4Yj4AwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQhgAwABCGADAAEIYAMAAQ0JYVZ1ZVR+rqmur6rwRPQAAjLTyEFZVxyR5WZIfTHJakqdW1Wmr7gMAYKQRe8IekeTa7v5Ed//fJK9NctaAPgAAhhkRwu6Z5K/WDH9qGgcAsGccO7qB9VTVuUnOnQa/UFUfG9nPOk5M8jfqq79Teqjnj62/BOrvsPoL/gzuuO3fZfW3Qw+j66/nW+aZaUQIuz7JN68Zvtc07na6+4IkF6yqqa2oqgPdvU999fdqD+qrr/7erb8dehhd/2iNOBz5viT3r6r7VNUdkjwlycUD+gAAGGble8K6+7aqemaSP0pyTJILu/vqVfcBADDSkHPCuvutSd46ovaCjT5cqv7erp+M70F99dXfu/WT8T2Mrn9UqrtH9wAAsOe4bREAwABCGADAAELYIQ53X8vpm5zvmca9bvpWZ6rqjtPwtdP0U9es5xem8R+rqh9YZf2qukdVXVZVX6iqly5x+7+vqt5fVbdV1ZMOWc85VXXN9HPOgPpvq6qbq+rNm9n+I/TwzGm4q+rENfNWVb1kmnZFVT10Sa/BVupv6TXYZP0HVNW7qurLVfWcjdaz4voXVtWNVXXVJmp/1TJVdUJVXTK9j5dU1fHT+GW894uqfzSf/830sIz3f1H1N/3+H6H+j1TV1VX1larad8j8h/1bv+Dt30r9pW9/HeHfmqp6WFVdOfX2kqqqFdf/tar6q6r6wma2f+W628/0k9m3Nf8yyX2T3CHJhzK7v+Xrkzxlmud3k/zU9PwZSX53ev6UJK+bnp82LXvHJPeZ1nnMCuvfOcmjkvxkkpcucftPTfJdSV6d5Elr1nNCkk9Mj8dPz49fVf1p2hlJ/kWSNy/oM/CQqd51SU5cM/8PJflfSSrJ6Unes6TXYFP1t/oabKH+NyR5eJJfS/KcjdazqvrTtO9L8tAkV21i+79qmSQvSHLe9Py8JM9fxnu/qPpH8/nfQg8Lff8XVX+r7/8R6n97km9L8vYk+9aMP+zf+iVs/6bqr3D71/23Jsl7p89lTZ/TH1xx/dOTnJLkC5v9HVjljz1ht7fefS0fk+QN0zz7k5w9PT9rGs40/Ywp7Z+V5LXd/eXu/j9Jrp3WvZL63f3F7n5nkr/f3OZvrn53X9fdVyT5yiHr+YEkl3T3Z7v7c0kuSXLmCuunuy9Ncut8m71xD939ge6+7jDzn5Xk1T3z7iR3r6pTsuDXYAv1t/oabKp+d9/Y3e9L8g/zrGeF9dPdf5rks3PU3GiZtb9nh/7+LfK9X1T9o/n8b6qHJbz/i6q/pfd/veW6+yPdfbi7tqz3t36h27+F+ivZ/vX+rZk+h8d197u7uzP7j/LZhy6/rPrTtHd39w3z1BxJCLu99e5reXN333bIuNvNP03/fJJ7HGE9q6q/VZutv9n1rKr+0dhs7+vNv+jXYFHzL7r+otezHe8te/KaP+afTnLy9HzR7/2i6i/Dej2sZ1WvwWir+gxstv5o98ysl4O2S1/bjhAGMKfpf/XDruszuv526GF0fVgkIez21ruv5d2r6thDxt1u/mn63ZL87RHWs6r6W7XZ+ptdz6rqH43N9r7e/It+DRY1/6LrL3o9i96eRfjMwcN80+ON0/hFv/eLqr8M6/WwnlW9BqOt6jOw2fqjXZ9ZLwdtl762HSHs9ta7r+VlSQ5+++6cJBdNzy+ehjNN/5Ppf2kXJ3lKzb69eJ8k98/sJMVV1d+qzdZfzx8leVxVHT99i+lx07hV1T8am7236cVJnlYzpyf5/HTYZNGvwWbrb9Wi7u261fVsx3vLrv09O/T3b5Hv/aLqL8N6Paxn0e/jZuuvynp/61f1Od7qvzVLNX0Ob6mq06fzpJ+W7fOebS+9Db4dsJ1+MvvG0ccz+5bJL07j7pvZB/vaJH+Q5I7T+DtNw9dO0++7Zj2/OK3jY5nzWyELrn9dZic3fiGz4/HzfjNnM/UfPq37i5ntgbt6zXr+7TT/tUn+zZK2/0j1/yzJTUm+NM3zA0fZw7Om9dyW5K+TvGIaX0leNs17ZW7/zZ1FvgZbqb+l12CT9b9xGn9Lkpun58ett54V139NkhsyO2n7U0mePkftr1oms/MsL01yTZI/TnLCEt/7RdU/ms//ZnpYxvu/qPqbfv+PUP9fTs+/nOQzSf5ozfyH/Vu/4O3fSv1Vbf91Ocy/NUn2Jblq6u2lme7Qs8L6L5iGvzI9/vK878Eqf9y2CABgAIcjAQAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwoBdp6qeXVVfv6j5AJbBJSqAXaeqrsvsull/s4j5AJbBnjBgR6uqO1fVW6rqQ1V1VVU9L8k3Jbmsqi6b5vmdqjpQVVdX1a9M4551mPkeV1Xvqqr3V9UfVNVdpvHnV9WHq+qKqvr1MVsK7Db2hAE7WlU9McmZ3f3vp+G7JflQ1uzhqqoTuvuzVXVMZldef1Z3X7F2T1hVnZjkDzO76vgXq+q5Se6Y2VXp/yLJA7q7q+ru3X3zyjcU2HXsCQN2uiuTPLaqnl9V39vdnz/MPE+uqvcn+UCS70hy2mHmOX0a/+dV9cHM7lH4LUk+n+Tvk7yyqv5Vkr9bxkYAe8+xoxsAOBrd/fGqemhm9+r71aq6dO306cbGz0ny8O7+XFW9KrP7rh6qklzS3U/9qglVj0hyRmY3kn9mkscsdiuAvcieMGBHq6pvSvJ33f37SV6Y5KFJbk1y12mW4zK7yfvnq+rkJD+4ZvG18707ySOr6n7Teu9cVd86nRd2t+5+a5KfSfKgZW8TsDfYEwbsdA9M8sKq+kqSf0jyU0m+O8nbquqvu/vRVfWBJB9N8ldJ/nzNshccMt9PJHlNVd1xmv5LmQW1i6rqTpntLfvZlWwVsOs5MR8AYACHIwEABhDCAAAGEMIAAAYQwgAABhDCAAAGEMIAAAYQwgAABvh/vlgnOeLX9S0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#必要なライブラリをインポートする\n",
    "from qulacs import QuantumState, QuantumCircuit, Observable, PauliOperator\n",
    "from qulacs.gate import H, CNOT, RX, RZ\n",
    "from scipy.optimize import minimize\n",
    "import numpy as np\n",
    "\n",
    "## 頂点の数\n",
    "n = 4 \n",
    "\n",
    "## C(Z)をqulacs.Observableとして定義\n",
    "cost_observable = Observable(n) \n",
    "for i in range(n):\n",
    "    cost_observable.add_operator( PauliOperator(\"Z {:} Z {:}\".format(i, (i+1)%n), 0.5) )\n",
    "\n",
    "# circuit に U_C(gamma) を加える関数\n",
    "def add_U_C(circuit, gamma):\n",
    "    for i in range(n):\n",
    "        j = (i+1) % n\n",
    "        circuit.add_CNOT_gate(i, j)\n",
    "        circuit.add_gate(RZ(j, -2*gamma)) ## qulacsでは RZ(theta)=e^{i*theta/2*Z}\n",
    "        circuit.add_CNOT_gate(i, j)\n",
    "    return circuit\n",
    "    \n",
    "# circuit に U_X(beta) を加える関数\n",
    "def add_U_X(circuit, beta):\n",
    "    for i in range(n):\n",
    "        circuit.add_gate(RX(i, -2*beta))\n",
    "    return circuit\n",
    "\n",
    "# p=2 の |beta, gamma> を作って <beta, gamma| C(Z) |beta, gamma> を返す関数\n",
    "# x = [beta0, beta1, gamma0, gamma1]\n",
    "def QAOA_output_twolayer(x): \n",
    "    beta0, beta1, gamma0, gamma1 = x\n",
    "\n",
    "    circuit = QuantumCircuit(n)\n",
    "    ## 重ね合わせを作るため、アダマールゲートをかける\n",
    "    for i in range(n):\n",
    "        circuit.add_H_gate(i)\n",
    "    ## U_C, U_Xをかける\n",
    "    circuit =  add_U_C(circuit, gamma0)\n",
    "    circuit =  add_U_X(circuit, beta0)\n",
    "    circuit =  add_U_C(circuit, gamma1)\n",
    "    circuit =  add_U_X(circuit, beta1)\n",
    "\n",
    "    ## |beta, gamma>を作る\n",
    "    state = QuantumState(n)\n",
    "    state.set_zero_state()                     \n",
    "    circuit.update_quantum_state(state)\n",
    "    return cost_observable.get_expectation_value(state)                 \n",
    "\n",
    "## 初期値\n",
    "x0 = np.array( [0.1, 0.1, 0.2, 0.3 ])\n",
    "\n",
    "## scipy.minimize を用いて最小化\n",
    "result = minimize(QAOA_output_twolayer, x0, options={'maxiter':500}, method='powell')\n",
    "print(result.fun) # 最適化後の値\n",
    "print(result.x) # 最適化後の[beta0, beta1, gamma0, gamma1]\n",
    "\n",
    "## 最適化後の状態を測定した時の確率分布を調べる\n",
    "beta0, beta1, gamma0, gamma1 = result.x\n",
    "\n",
    "circuit = QuantumCircuit(n)\n",
    "## 重ね合わせを作るため、アダマールゲートをかける\n",
    "for i in range(n):\n",
    "    circuit.add_H_gate(i)\n",
    "## U_C, U_Xをかける\n",
    "circuit =  add_U_C(circuit, gamma0)\n",
    "circuit =  add_U_X(circuit, beta0)\n",
    "circuit =  add_U_C(circuit, gamma1)\n",
    "circuit =  add_U_X(circuit, beta1)\n",
    "\n",
    "## |beta, gamma>を作る\n",
    "state = QuantumState(n)\n",
    "state.set_zero_state()                     \n",
    "circuit.update_quantum_state(state)\n",
    "\n",
    "## 状態ベクトルの各成分の絶対値の二乗=観測確率\n",
    "probs = np.abs(state.get_vector())**2\n",
    "print(probs)\n",
    "\n",
    "## z方向に射影測定した時に得られる可能性があるビット列\n",
    "z_basis = [format(i,\"b\").zfill(n) for i in range(probs.size)]\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.xlabel(\"states\")\n",
    "plt.ylabel(\"probability(%)\")\n",
    "plt.bar(z_basis, probs*100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$p=1$の時に比べ、圧倒的に大きい確率で真の解$|0101\\rangle, |1010\\rangle$が得られる確率が高いことが分かる。また、コスト関数の値も正しく−2に近づいている。\n",
    "\n",
    "このように、QAOAを用いる際には、変分量子回路の複雑さ$p$の大きさにも注意しながら実装する必要があるだろう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 参考文献\n",
    "\n",
    "[1] E. Farhi, J. Goldstone, and S. Gutmann, “A Quantum Approximate Optimization Algorithm”, [arXiv:1411.4028](https://arxiv.org/abs/1411.4028) (2014).\n",
    "\n",
    "[2] Eddie Farhi: A Quantum Approximate Optimization Algorithm, https://www.youtube.com/watch?v=J8y0VhnISi8"
   ]
  }
 ],
 "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.7.5"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}