source/ja/nonlocal_correlations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 【課題】量子相関を調べる\n",
    "\n",
    "第一回の実習ではCHSH不等式の破れを調べるために、2つの量子ビットの相関関数$C^{i} \\, (i=0,1,2,3)$という量を量子コンピュータを使って計算しました。この課題では、この量をもう少し細かく調べてみましょう。\n",
    "\n",
    "```{contents} 目次\n",
    "---\n",
    "local: true\n",
    "---\n",
    "```\n",
    "\n",
    "$\\newcommand{\\ket}[1]{|#1\\rangle}$\n",
    "$\\newcommand{\\bra}[1]{\\langle#1|}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QCシミュレータの使い方\n",
    "\n",
    "実習で見たように、QCで実現される量子状態は、量子力学の公理に基づいて理論的に計算・予測できます。そこで用いられる数学的操作も単なる足し算や掛け算（線形演算）なので、実はQCの量子状態は（古典）計算機で比較的簡単に計算できます。当然のことですが、QCは何も魔法のブラックボックスというわけではありません。\n",
    "\n",
    "ただし、古典計算機で量子状態を再現するためには、特殊な場合を除いて、量子ビット数の指数関数的な量のメモリが必要になります。これも前半で見たように、$n$量子ビットあれば、系の自由度（degrees of freedom / dof: 実数自由パラメータの数）は$2^{n+1} - 2$ですので、例えば各自由度を64ビット（＝8バイト）の浮動小数点で表現するとしたら、必要なメモリは(-2を無視して)\n",
    "\n",
    "$$\n",
    "2^3\\, \\mathrm{(bytes / dof)} \\times 2^{n+1}\\, \\mathrm{(dof)} = 2^{n+4}\\, \\mathrm{(bytes)}\n",
    "$$\n",
    "\n",
    "なので、$n=16$で1 MiB、$n=26$で1 GiB、$n=36$で1 TiBです。現在の計算機では、ハイエンドワークステーションでRAMが$\\mathcal{O}(1)$ TiB、スパコン「富岳」で5 PB (~2<sup>52</sup> bytes)なのに対し、QCではすでに$n=65$のものが存在するので、既に古典計算機でまともにシミュレートできない機械が存在していることになります。\n",
    "\n",
    "しかし、逆に言うと、$n \\sim 30$程度までの回路であれば、ある程度のスペックを持った計算機で厳密にシミュレートできるということが言えます。じっさい世の中には[数多くの](https://quantiki.org/wiki/list-qc-simulators)シミュレータが存在します。Qiskitにも様々な高機能シミュレータが同梱されています。\n",
    "\n",
    "シミュレーションはローカル（手元のPythonを動かしているコンピュータ）で実行できるので、ジョブを投げて結果を待つ時間が省けます。この課題ではたくさんの細かい量子計算をするので、実機を使わず、`qasm_simulator`というQiskitに含まれるシミュレータを利用します。\n",
    "\n",
    "Qiskitのシミュレータには`Aer`というオブジェクトからアクセスします。`Aer`は実習で登場した`IBMQ`と同様の構造をしており、複数のシミュレータをバックエンドとして管理しています。そのうちの`qasm_simulator`を取り出します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# まずは全てインポート\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize, Bounds\n",
    "from qiskit import QuantumCircuit, Aer, execute\n",
    "from qiskit.visualization import plot_histogram\n",
    "\n",
    "print('notebook ready')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulator = Aer.get_backend('qasm_simulator')\n",
    "print(simulator.name())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "実習の内容を再現してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "circuits = []\n",
    "for ic in range(4):\n",
    "    circuit = QuantumCircuit(2, name='circuit{}'.format(ic))\n",
    "    circuit.h(0)\n",
    "    circuit.cx(0, 1)\n",
    "    circuits.append(circuit)\n",
    "\n",
    "circuits[0].ry(-np.pi / 4., 1)\n",
    "circuits[1].ry(-3. * np.pi / 4., 1)\n",
    "circuits[2].ry(-np.pi / 4., 1)\n",
    "circuits[3].ry(-3. * np.pi / 4., 1)\n",
    "\n",
    "circuits[2].ry(-np.pi / 2., 0)\n",
    "circuits[3].ry(-np.pi / 2., 0)\n",
    "\n",
    "for circuit in circuits:\n",
    "    circuit.measure_all()\n",
    "    \n",
    "shots = 10000\n",
    "    \n",
    "sim_job = execute(circuits, backend=simulator, shots=shots)\n",
    "\n",
    "sim_result = sim_job.result()\n",
    "\n",
    "C = []\n",
    "for circuit in circuits:\n",
    "    c = sim_result.get_counts(circuit)\n",
    "    ax = plt.figure().add_subplot()\n",
    "    plot_histogram(c, ax=ax)\n",
    "    \n",
    "    C.append((c['00'] + c['11'] - c['01'] - c['10']) / shots)\n",
    "    \n",
    "S = C[0] - C[1] + C[2] + C[3]\n",
    "print('S =', S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上のように、`qasm_simulator`は実機と同様に`execute`関数で実行でき、ヒストグラムデータを返します。実機ではショット数に制限がありますが、シミュレータにはありません。ただしショット数が多いほど、当然実行に時間がかかります。といってもこの程度の回路であれば常識的なショット数ならほぼ瞬間的にジョブの実行が終わるので、上の例では実習で使った`job_monitor()`関数を使用していません。また、シミュレータにはノイズがない[^simulator_noise]ので、$S$の計算結果が統計誤差の範囲内で理論値と一致していることが見て取れます。\n",
    "\n",
    "[^simulator_noise]: 標準設定において。実機の振る舞いに従うよう、あえてノイズを加えるような設定も存在します。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 測定基底の変換\n",
    "\n",
    "さて、おさらいをすると、上の$C^{0,1,2,3}$を計算する4つの回路は以下のようなものでした。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "for circuit in circuits:\n",
    "    ax = plt.figure().add_subplot()\n",
    "    circuit.draw('mpl', ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ベル状態を作るところまですべての回路で共通で、その後それぞれ異なる角度の$R_y$ゲートをかけています。実習では深く立ち入りませんでしたが、この$R_y$ゲートはベル状態を$\\{\\ket{0}, \\ket{1}\\}$とは異なる基底で測定するために置かれています。どういうことか、以下で説明します。\n",
    "\n",
    "ここまで「測定」とはレジスタの量子状態$\\sum_{j=0}^{2^n-1} c_j \\ket{j}$からビット列$j$を確率$|c_j|^2$で得る行為である、と説明してきました。しかし、本来「測定」はもっと一般的な概念です。それを理解するために、まず、量子力学的には計算基底状態$\\ket{j}$は何ら特別な状態ではないということを理解しておきましょう。\n",
    "\n",
    "例えば1量子ビットにおいて、計算基底状態は$\\ket{0}$と$\\ket{1}$ですが、以下のような状態$\\ket{\\theta}$と$\\ket{\\theta + \\pi}$を考えます。\n",
    "\n",
    "```{math}\n",
    ":label: theta_ket_def\n",
    "\\ket{\\theta} := R_y(\\theta)\\ket{0} = \\cos\\frac{\\theta}{2}\\ket{0} + \\sin\\frac{\\theta}{2}\\ket{1} \\\\\n",
    "\\ket{\\theta + \\pi} := R_y(\\theta)\\ket{1} = -\\sin\\frac{\\theta}{2}\\ket{0} + \\cos\\frac{\\theta}{2}\\ket{1}\n",
    "```\n",
    "\n",
    "すると、\n",
    "\n",
    "$$\n",
    "\\ket{0} = \\cos\\frac{\\theta}{2}\\ket{\\theta} - \\sin\\frac{\\theta}{2}\\ket{\\theta + \\pi} \\\\\n",
    "\\ket{1} = \\sin\\frac{\\theta}{2}\\ket{\\theta} + \\cos\\frac{\\theta}{2}\\ket{\\theta + \\pi},\n",
    "$$\n",
    "\n",
    "つまり、計算基底状態が$\\ket{\\theta}$と$\\ket{\\theta + \\pi}$の重ね合わせとして表現できます。量子ビットの任意の状態は$\\ket{0}$と$\\ket{1}$の重ね合わせで表現できるので、$\\ket{\\theta}$と$\\ket{\\theta + \\pi}$の重ね合わせでも表現できるということになります。そのようなときは$\\ket{\\theta}$と$\\ket{\\theta + \\pi}$を基底として状態を表しているわけです。\n",
    "\n",
    "一般に、量子力学的には、2つの異なる状態$\\ket{a} \\nsim \\ket{b}$を考えれば、それらの重ね合わせで量子ビットの任意の状態が表現できます。そして、$\\ket{a}$と$\\ket{b}$が直交する[^orthogonal]ときに状態$\\ket{\\psi}$が\n",
    "\n",
    "$$\n",
    "\\ket{\\psi} = \\alpha \\ket{a} + \\beta \\ket{b}\n",
    "$$\n",
    "\n",
    "と表現されるならば、「基底$\\ket{a}$と$\\ket{b}$についての測定」という操作を考えることができます。$\\ket{\\psi}$に対してそのような測定をすると、状態$\\ket{a}$が確率$|\\alpha|^2$で、状態$\\ket{b}$が確率$|\\beta|^2$で得られます。\n",
    "\n",
    "量子計算においても、アルゴリズムの一部として、計算の結果実現した状態を特定の基底で測定するということが多々あります。ところが、ここで若干問題があります。量子コンピュータは実装上、計算基底でしか測定ができないのです。量子力学の理論的には特別でない$\\ket{0}$と$\\ket{1}$ですが、量子コンピュータという実態にとっては具体的な対応物があるのです。\n",
    "\n",
    "そこで、量子計算では、状態を任意の基底で測定することを諦め、反対に状態を変化させてしまいます。例えば、本当は上の$\\ket{\\theta}$と$\\ket{\\theta + \\pi}$という基底で量子ビットの状態$\\ket{\\psi} = \\alpha \\ket{\\theta} + \\beta \\ket{\\theta + \\pi}$を測定したいとします。しかし計算基底でしか測定ができないので、代わりに$R_y(-\\theta)$を$\\ket{\\psi}$にかけます。すると式{eq}`theta_ket_def`から\n",
    "\n",
    "$$\n",
    "R_y(-\\theta)\\ket{\\theta} = \\ket{0} \\\\\n",
    "R_y(-\\theta)\\ket{\\theta + \\pi} = \\ket{1}\n",
    "$$\n",
    "\n",
    "なので、\n",
    "\n",
    "$$\n",
    "R_y(-\\theta)\\ket{\\psi} = \\alpha \\ket{0} + \\beta \\ket{1}\n",
    "$$\n",
    "\n",
    "が得られます。この$R_y(-\\theta)\\ket{\\psi}$を計算基底で測定した結果は、$\\ket{\\psi}$を$\\ket{\\theta}, \\ket{\\theta + \\pi}$基底で測定した結果と等価です。\n",
    "\n",
    "このように、測定を行いたい基底（ここでは$\\ket{\\theta}, \\ket{\\theta + \\pi}$）を$\\ket{0}, \\ket{1}$から得るための変換ゲート（$R_y(\\theta)$）の逆変換を測定したい状態にかけることで、計算基底での測定で求める結果を得ることができます。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 観測量の期待値とその計算法\n",
    "\n",
    "課題の説明に入る前に、さらに話が込み入ってきますが、量子計算でも多出する（ワークブックでは特に{doc}`vqe`以降）概念である「観測量の期待値」について説明します。\n",
    "\n",
    "観測量とはそのまま「観測できる量」のことで、量子状態から取り出せる（古典的）情報のこととも言えます。例えば、何かしらの粒子の運動を量子力学的に記述した場合、その粒子の位置や運動量などが観測量です。\n",
    "\n",
    "ケットで表される量子状態に対して、観測量はケットに作用する「エルミート演算子」で表現されます。細かい定義は参考文献に譲り、ここで必要な最小限のことだけ述べると、エルミート演算子は**対角化可能で、実数の固有値を持つ**という性質を持っています。つまり、$A$が$N$次元の量子状態の空間のエルミート演算子であれば、\n",
    "\n",
    "$$\n",
    "A \\ket{\\phi_j} = a_j \\ket{\\phi_j}\n",
    "$$\n",
    "\n",
    "が成り立つような状態$\\ket{\\phi_j}$と実数$a_j$の組が$N$個存在し、各$\\ket{\\phi_j} \\, (j=0,\\dots,N-1)$は互いに直交します。このとき、$\\ket{\\phi_j}$を固有ベクトル、$a_j$を固有値と呼びます。そして、この固有値$a_j$が、演算子$A$で表される観測量の値に対応します[^continuous_observables]。\n",
    "\n",
    "さて、$A$が$n$ビット量子レジスタに対する演算子であるとします。レジスタの状態$\\ket{\\psi}$が$A$の固有ベクトルで\n",
    "\n",
    "$$\n",
    "\\ket{\\psi} = \\sum_{j=0}^{2^n-1} \\gamma_j \\ket{\\phi_j}\n",
    "$$\n",
    "\n",
    "と分解されるとき、この状態を固有ベクトル基底$\\{\\ket{\\phi_j}\\}$で測定すると、状態$\\ket{\\phi_j}$が確率$|\\gamma_j|^2$で得られます。そして、状態が$\\ket{\\phi_j}$であれば観測量$A$の値は$a_j$です。\n",
    "\n",
    "そのような測定を多数回繰り返して$A$の値の期待値を求めることを考えます。期待値の定義は「確率変数のすべての値に確率の重みを付けた加重平均」です。$A$が確率分布$\\{|\\gamma|^2_j\\}$に従って値$\\{a_j\\}$を取るので、\n",
    "\n",
    "$$\n",
    "\\bra{\\psi} A \\ket{\\psi} = \\sum_{j=0}^{2^n-1} a_j |\\gamma_j|^2\n",
    "$$\n",
    "\n",
    "となります。ここから、量子コンピュータにおいて観測量の期待値を計算する方法を見出すことができます。具体的には、\n",
    "\n",
    "1. 観測量をエルミート演算子で表現する\n",
    "1. 演算子を対角化し、固有値と対応する固有ベクトルを求める\n",
    "1. 固有ベクトルを基底として、レジスタの状態を測定する\n",
    "1. 測定から得られた確率分布を重みとして固有値の平均値を取る\n",
    "\n",
    "です。3の測定の際には、上のセクションで説明した測定基底の変換を利用します。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CHSH不等式を見直す\n",
    "\n",
    "実習の中で、\n",
    "\n",
    "$$\n",
    "C = P_{00} - P_{01} - P_{10} + P_{11}\n",
    "$$\n",
    "\n",
    "という量を計算しました。ここで、$P_{lm}$は2つの量子ビットでそれぞれ$l, m \\, (=0,1)$が得られる確率でした。実はこの量$C$は、2ビットレジスタにおけるある観測量の期待値として捉えることができます。\n",
    "\n",
    "まず、1つの量子ビットに対して、固有値が$\\pm 1$であるような観測量$\\sigma^{\\theta}$を考えます。そのような$\\sigma^{\\theta}$は無数に存在し、固有ベクトルで区別できます。$\\theta$は固有ベクトルを決める何らかのパラメータです[^specifying_eigenvectors]。例えば$\\sigma^0$という観測量を、計算基底を固有ベクトルとして、\n",
    "\n",
    "$$\n",
    "\\sigma^0 \\ket{0} = \\ket{0} \\\\\n",
    "\\sigma^0 \\ket{1} = -\\ket{1}\n",
    "$$\n",
    "\n",
    "で定義できる、という具合です。\n",
    "\n",
    "次に、2つの量子ビットA, Bからなるレジスタを考え、$\\sigma^{\\kappa}$をAの、$\\sigma^{\\lambda}$をBの観測量とします。また、それぞれの演算子の固有ベクトルを\n",
    "\n",
    "$$\n",
    "\\sigma^{\\kappa} \\ket{\\kappa_{\\pm}} = \\pm \\ket{\\kappa_{\\pm}} \\\\\n",
    "\\sigma^{\\lambda} \\ket{\\lambda_{\\pm}} = \\pm \\ket{\\lambda_{\\pm}}\n",
    "$$\n",
    "\n",
    "で定義します。これらの固有ベクトルを使って、レジスタの状態$\\ket{\\psi}$を\n",
    "\n",
    "$$\n",
    "\\ket{\\psi} = c_{++} \\ket{\\lambda_+}_B \\ket{\\kappa_+}_A + c_{+-} \\ket{\\lambda_+}_B \\ket{\\kappa_-}_A + c_{-+} \\ket{\\lambda_-}_B \\ket{\\kappa_+}_A + c_{--} \\ket{\\lambda_-}_B \\ket{\\kappa_-}_A\n",
    "$$\n",
    "\n",
    "と分解します。ケットはAが右、Bが左になるよう並べました。すると、積$\\sigma^{\\lambda}_B \\sigma^{\\kappa}_A$の$\\ket{\\psi}$に関する期待値は、\n",
    "\n",
    "$$\n",
    "\\bra{\\psi} \\sigma^{\\lambda}_B \\sigma^{\\kappa}_A \\ket{\\psi} = |c_{++}|^2 - |c_{+-}|^2 - |c_{-+}|^2 + |c_{--}|^2\n",
    "$$\n",
    "\n",
    "です。\n",
    "\n",
    "最後に、同じ結果を計算基底での測定で表すために、$\\ket{\\psi}$に対して基底変換を施します。$\\{\\ket{\\kappa_{\\pm}}, \\ket{\\lambda_{\\pm}}\\}$が何らかのパラメータ付きゲート$R(\\theta)$を通して計算基底と\n",
    "\n",
    "$$\n",
    "\\ket{\\kappa_+} = R(\\kappa) \\ket{0} \\\\\n",
    "\\ket{\\kappa_-} = R(\\kappa) \\ket{1} \\\\\n",
    "\\ket{\\lambda_+} = R(\\lambda) \\ket{0} \\\\\n",
    "\\ket{\\lambda_-} = R(\\lambda) \\ket{1}\n",
    "$$\n",
    "\n",
    "で結びついているなら、状態$\\ket{\\psi'} = R^{-1}_B(\\lambda) R^{-1}_A(\\kappa) \\ket{\\psi}$を計算基底で測定したとき、\n",
    "\n",
    "$$\n",
    "P_{00} = |c_{++}|^2 \\\\\n",
    "P_{01} = |c_{+-}|^2 \\\\\n",
    "P_{10} = |c_{-+}|^2 \\\\\n",
    "P_{11} = |c_{--}|^2\n",
    "$$\n",
    "\n",
    "が成り立ちます。確認のためはっきりさせておくと、左辺は$\\ket{\\psi'}$を計算基底で測定し、ビット列00, 01, 10, 11を得る確率です。つまり、最初の$C$は\n",
    "\n",
    "$$\n",
    "C = \\bra{\\psi'} \\sigma^0_B \\sigma^0_A \\ket{\\psi'} = \\bra{\\psi} \\sigma^{\\lambda}_B \\sigma^{\\kappa}_A \\ket{\\psi}\n",
    "$$\n",
    "\n",
    "を表していたのでした。\n",
    "\n",
    "これを踏まえて、CHSH不等式の左辺は結局何を計算していたのか、見直してみましょう。ここでベル状態を\n",
    "\n",
    "$$\n",
    "\\ket{\\Psi} = \\frac{1}{\\sqrt{2}} \\left(\\ket{00} + \\ket{11}\\right)\n",
    "$$\n",
    "\n",
    "とおき、第0ビットをA、第1ビットをBとします。$R_y(\\pi/4)\\ket{0}, R_y(\\pi/2)\\ket{0}, R_y(3\\pi/4)\\ket{0}$が固有値$+1$の固有ベクトルとなるような演算子をそれぞれ$\\sigma^{\\pi/4}, \\sigma^{\\pi/2}, \\sigma^{3\\pi/4}$とすると、\n",
    "\n",
    "$$\n",
    "S = C^0 - C^1 + C^2 + C^3 = \\bra{\\Psi} \\sigma^{\\pi/4}_B \\sigma^0_A \\ket{\\Psi} - \\bra{\\Psi} \\sigma^{3\\pi/4}_B \\sigma^0_A \\ket{\\Psi} + \\bra{\\Psi} \\sigma^{\\pi/4}_B \\sigma^{\\pi/2}_A \\ket{\\Psi} + \\bra{\\Psi} \\sigma^{3\\pi/4}_B \\sigma^{\\pi/2}_A \\ket{\\Psi}\n",
    "$$\n",
    "\n",
    "がわかります。\n",
    "\n",
    "観測量$\\sigma^{\\theta}$を用いてCHSH不等式をより正確に表現すると、\n",
    "\n",
    "> 4つのパラメータ$\\kappa, \\lambda, \\mu, \\nu$を用いて  \n",
    "> $S(\\kappa, \\lambda, \\mu, \\nu) = \\langle \\sigma^{\\kappa}\\sigma^{\\lambda} \\rangle - \\langle \\sigma^{\\kappa}\\sigma^{\\nu} \\rangle + \\langle \\sigma^{\\mu}\\sigma^{\\lambda} \\rangle + \\langle \\sigma^{\\mu}\\sigma^{\\nu} \\rangle$  \n",
    "> という量を定義すると、エンタングルメントのない古典力学において$|S| \\leq 2$である\n",
    "\n",
    "となります。\n",
    "\n",
    "実習で用いた$\\sigma^{\\theta}$のパラメータは、$|S|$の値を最大化するものでした。次のセルで数値的にこのことを確かめています。\n",
    "\n",
    "```python tags=[\"remove-output\"]\n",
    "def quantum_S(x):\n",
    "    \"\"\"\n",
    "    Three-parameter function to be minimized. Returns -|S| = -|<sigma^k sigma^l> - <sigma^k sigma^n> + <sigma^m sigma^l> + <sigma^m sigma^n>| with k fixed to 0.\n",
    "    \"\"\"\n",
    "    \n",
    "    return -np.abs(np.cos(-x[0]) - np.cos(-x[2]) + np.cos(x[1] - x[0]) + np.cos(x[1] - x[2]))\n",
    "\n",
    "# Initial values\n",
    "x0 = np.array([np.pi / 2., 0., np.pi / 2.])\n",
    "# Bounds on all parameters ([0, pi])\n",
    "bounds = Bounds([0.] * 3, [np.pi] * 3)\n",
    "# Minimize using scipy.optimize.minimize\n",
    "res = minimize(quantum_S, x0, method='trust-constr', bounds=bounds)\n",
    "print('argmax |S|: kappa=0, lambda={:.2f}pi, mu={:.2f}pi, nu={:.2f}pi'.format(*(res.x / np.pi)))\n",
    "```\n",
    "\n",
    "[^orthogonal]: ここまで状態ケットの内積に関する説明をしていませんが、線形代数に慣れている方は「量子力学の状態ケットはヒルベルト空間の要素である」と理解してください。慣れていない方は、量子ビットの直行する状態とは$\\ket{0}$と$\\ket{1}$に同じゲートをかけて得られる状態のことを言うと覚えておいてください。\n",
    "[^continuous_observables]: 上で粒子の位置や運動量という例を挙げたので、$a_j \\, (j=0,\\dots,N-1)$と離散的な固有値が観測量の値だというと混乱するかもしれません。位置や運動量といった連続的な値を持つ観測量は、無数に固有値を持つエルミート演算子で表されます。そのような演算子のかかる状態空間は無限次元です。\n",
    "[^specifying_eigenvectors]: 量子ビットの一般の状態は実パラメータ2つで決まるので、$\\sigma^{\\theta, \\phi}$などと書いたほうがより明示的ですが、ここでの議論では結局1パラメータしか使わないので、「何か一般の（次元を指定しない）パラメータ」として$\\theta$と置いています。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 問題: ベル状態と可分状態の違い\n",
    "\n",
    "### 一般の$\\sigma$演算子の期待値\n",
    "\n",
    "上のように$R_y(\\theta)\\ket{0}$が固有値$+1$の固有ベクトルとなるような演算子を$\\sigma^{\\theta}$として、$\\bra{\\Psi} \\sigma^{\\phi}_B \\sigma^{\\theta}_A \\ket{\\Psi}$を計算してみましょう。まず基底変換を具体的に書き下します。\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\ket{\\Psi'} & = R_{y,B}(-\\phi) R_{y,A}(-\\theta) \\frac{1}{\\sqrt{2}} \\left( \\ket{00} + \\ket{11} \\right) \\\\\n",
    "& = \\frac{1}{\\sqrt{2}} \\left\\{ \\left[ \\cos\\left(\\frac{\\phi}{2}\\right)\\ket{0} - \\sin\\left(\\frac{\\phi}{2}\\right)\\ket{1} \\right] \\left[ \\cos\\left(\\frac{\\theta}{2}\\right)\\ket{0} - \\sin\\left(\\frac{\\theta}{2}\\right)\\ket{1} \\right] + \\left[ \\sin\\left(\\frac{\\phi}{2}\\right)\\ket{0} + \\cos\\left(\\frac{\\phi}{2}\\right)\\ket{1} \\right] \\left[ \\sin\\left(\\frac{\\theta}{2}\\right)\\ket{0} + \\cos\\left(\\frac{\\theta}{2}\\right)\\ket{1} \\right] \\right\\} \\\\\n",
    "& = \\frac{1}{\\sqrt{2}} \\left[ \\cos\\left(\\frac{\\theta - \\phi}{2}\\right)\\ket{00} - \\sin\\left(\\frac{\\theta - \\phi}{2}\\right)\\ket{01} + \\sin\\left(\\frac{\\theta - \\phi}{2}\\right)\\ket{10} + \\cos\\left(\\frac{\\theta - \\phi}{2}\\right)\\ket{11} \\right].\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "したがって、\n",
    "\n",
    "```{math}\n",
    ":label: quantum_correlation\n",
    "\\begin{align}\n",
    "\\bra{\\Psi} \\sigma^{\\phi}_B \\sigma^{\\theta}_A \\ket{\\Psi} & = \\bra{\\Psi'} \\sigma^{0}_B \\sigma^{0}_A \\ket{\\Psi'} \\\\\n",
    "& = \\cos^2\\left(\\frac{\\theta - \\phi}{2}\\right) - \\sin^2\\left(\\frac{\\theta - \\phi}{2}\\right) \\\\\n",
    "& = \\cos(\\theta - \\phi)\n",
    "\\end{align}\n",
    "```\n",
    "\n",
    "となります。\n",
    "\n",
    "### 問題\n",
    "\n",
    "上の計算結果を量子回路でも確認してみましょう。実習のように2ビット量子レジスタをベル状態にし、2つの量子ビットに適当な$R_y$ゲートをかけ、期待値$C$を$R_y$ゲートのパラメータの値の差の関数としてプロットします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception",
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "# Define theta and phi values\n",
    "thetas = [0., np.pi / 4., np.pi / 2., 3. * np.pi / 4., np.pi]\n",
    "theta_labels = ['0', r'$\\pi/4$', r'$\\pi/2$', r'$3\\pi/4$', r'$\\pi$']\n",
    "phis = np.linspace(0., np.pi, 16, endpoint=True)\n",
    "\n",
    "# Construct a circuit for each (theta, phi) pair\n",
    "circuits = []\n",
    "for itheta, theta in enumerate(thetas):\n",
    "    for iphi, phi in enumerate(phis):\n",
    "        circuit = QuantumCircuit(2, name='circuit_{}_{}'.format(itheta, iphi))\n",
    "\n",
    "        ##################\n",
    "        ### EDIT BELOW ###\n",
    "        ##################\n",
    "\n",
    "        #circuit.?\n",
    "        \n",
    "        ##################\n",
    "        ### EDIT ABOVE ###\n",
    "        ##################\n",
    "        \n",
    "        circuit.measure_all()\n",
    "\n",
    "        circuits.append(circuit)\n",
    "\n",
    "# Execute the circuit in qasm_simulator and retrieve the results\n",
    "simulator = Aer.get_backend('qasm_simulator')\n",
    "shots = 10000\n",
    "sim_job = execute(circuits, backend=simulator, shots=shots)\n",
    "result = sim_job.result()\n",
    "\n",
    "# Plot C versus (theta - phi) for each theta\n",
    "icirc = 0\n",
    "for itheta, theta in enumerate(thetas):\n",
    "    x = theta - phis\n",
    "    y = np.zeros_like(x)\n",
    "\n",
    "    for iphi, phi in enumerate(phis):\n",
    "        c = result.get_counts(circuits[icirc])\n",
    "\n",
    "        ##################\n",
    "        ### EDIT BELOW ###\n",
    "        ##################\n",
    "\n",
    "        #y[iphi] = ?\n",
    "        \n",
    "        ##################\n",
    "        ### EDIT ABOVE ###\n",
    "        ##################\n",
    "\n",
    "        icirc += 1\n",
    "\n",
    "    plt.plot(x, y, 'o', label=theta_labels[itheta])\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(r'$\\theta - \\phi$')\n",
    "plt.ylabel(r'$\\langle \\Psi | \\sigma^{\\phi}_B \\sigma^{\\theta}_A | \\Psi \\rangle$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "式{eq}`quantum_correlation`では、2体系の期待値がパラメータ$\\theta$と$\\phi$だけを含む関数の積として記述できないことがはっきり示されています。これはベル状態の持つエンタングルメントの現れで、2体系をバラバラに考えることができないということを言っています。\n",
    "\n",
    "それでは、$C$や$S$をエンタングルメントのない状態（可分状態）に対して計算してみたら、何が得られるでしょうか。例えば、ベル状態$1/\\sqrt{2}(\\ket{00} + \\ket{11})$の代わりに、「確率1/2で$\\ket{00}$、確率1/2で$\\ket{11}$」という状態を考えます。\n",
    "\n",
    "まずは$|S(\\kappa, \\lambda, \\mu, \\nu)|$の最大値を求めます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "def classical_S(params):\n",
    "    \"\"\"\n",
    "    Four-parameter function to be minimized. Returns -|S| = -|<sigma^k sigma^l> - <sigma^k sigma^n> + <sigma^m sigma^l> + <sigma^m sigma^n>|.\n",
    "    \n",
    "    Args:\n",
    "        params (np.ndarray): Values of kappa, lambda, mu, nu\n",
    "        \n",
    "    Returns:\n",
    "        float: -|S|\n",
    "    \"\"\"\n",
    "    \n",
    "    k, l, m, n = params\n",
    "    \n",
    "    S = 0.\n",
    "    \n",
    "    # S from |00>\n",
    "    ##################\n",
    "    ### EDIT BELOW ###\n",
    "    ##################\n",
    "    #sksl_00 = \n",
    "    #sksn_00 = \n",
    "    #smsl_00 = \n",
    "    #smsn_00 = \n",
    "    #S += (sksl_00 - sksn_00 + smsl_00 + smsn_00) * 0.5\n",
    "    ##################\n",
    "    ### EDIT ABOVE ###\n",
    "    ##################\n",
    "    \n",
    "    # S from |11>\n",
    "    ##################\n",
    "    ### EDIT BELOW ###\n",
    "    ##################\n",
    "    #sksl_11 = \n",
    "    #sksn_11 = \n",
    "    #smsl_11 = \n",
    "    #smsn_11 = \n",
    "    #S += (sksl_11 - sksn_11 + smsl_11 + smsn_11) * 0.5\n",
    "    ##################\n",
    "    ### EDIT ABOVE ###\n",
    "    ##################\n",
    "    \n",
    "    return -np.abs(S)\n",
    "\n",
    "# Initial values\n",
    "x0 = np.array([0., np.pi / 2., 0., np.pi / 2.])\n",
    "# Bounds on all parameters ([0, pi])\n",
    "bounds = Bounds([0.] * 4, [np.pi] * 4)\n",
    "# Minimize using scipy.optimize.minimize\n",
    "res = minimize(classical_S, x0, method='trust-constr', bounds=bounds)\n",
    "print('argmax |S|: kappa={:.3f}pi, lambda={:.3f}pi, mu={:.3f}pi, nu={:.3f}pi'.format(*(res.x / np.pi)))\n",
    "print('max |S|:', -res.fun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "source": [
    "この状態を量子回路で実装して、上の実験を繰り返してみましょう。ただし、実際に確率1/2で異なる状態が現れるような「混合状態」は量子回路では表現しにくいので、$\\ket{00}$と$\\ket{11}$を初期状態とする回路を一つずつ用意して、それらから得る$C$の平均値を求めることにします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "# Construct a circuit for each (theta, phi) pair\n",
    "circuits = []\n",
    "for itheta, theta in enumerate(thetas):\n",
    "    for iphi, phi in enumerate(phis):\n",
    "        circuit_00 = QuantumCircuit(2, name='circuit_{}_{}_00'.format(itheta, iphi))\n",
    "        circuit_11 = QuantumCircuit(2, name='circuit_{}_{}_11'.format(itheta, iphi))\n",
    "\n",
    "        ##################\n",
    "        ### EDIT BELOW ###\n",
    "        ##################\n",
    "\n",
    "        #circuit_00.?\n",
    "        #circuit_11.?\n",
    "        \n",
    "        ##################\n",
    "        ### EDIT ABOVE ###\n",
    "        ##################\n",
    "        \n",
    "        circuit_00.measure_all()\n",
    "        circuit_11.measure_all()\n",
    "\n",
    "        circuits.append(circuit_00)\n",
    "        circuits.append(circuit_11)\n",
    "\n",
    "# Execute the circuit in qasm_simulator and retrieve the results\n",
    "shots = 10000\n",
    "sim_job = execute(circuits, backend=simulator, shots=shots)\n",
    "result = sim_job.result()\n",
    "\n",
    "# Plot C versus (theta - phi) for each theta\n",
    "icirc = 0\n",
    "for itheta, theta in enumerate(thetas):\n",
    "    x = theta - phis\n",
    "    y = np.zeros_like(x)\n",
    "\n",
    "    for iphi, phi in enumerate(phis):\n",
    "        c_00 = result.get_counts(circuits[icirc])\n",
    "        c_11 = result.get_counts(circuits[icirc + 1])\n",
    "\n",
    "        ##################\n",
    "        ### EDIT BELOW ###\n",
    "        ##################\n",
    "\n",
    "        #y[iphi] = ?\n",
    "        \n",
    "        ##################\n",
    "        ### EDIT ABOVE ###\n",
    "        ##################\n",
    "\n",
    "        icirc += 2\n",
    "\n",
    "    plt.plot(x, y, 'o', label=theta_labels[itheta])\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(r'$\\theta - \\phi$')\n",
    "plt.ylabel(r'$\\langle \\Psi | \\sigma^{\\theta}_A \\sigma^{\\phi}_B | \\Psi \\rangle$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "raises-exception",
     "remove-output"
    ]
   },
   "source": [
    "ベル状態からの結果と比べてみて、何が言えるでしょうか。\n",
    "\n",
    "**提出するもの**\n",
    "\n",
    "- 完成した回路のコード（EDIT BELOW / EDIT ABOVEの間を埋める）とシミュレーション結果によるプロット\n",
    "- ベル状態と可分状態の混合状態とでの2体相関の違いに関する考察"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}