docs/source/optimal_povm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XdRhxg_msiR5"
   },
   "source": [
    "# Optimal POVM\n",
    "\n",
    "One can open this tutorial in Google Colab (is recommended)\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/LuchnikovI/QGOpt/blob/master/docs/source/optimal_povm.ipynb)\n",
    "\n",
    "In the following tutorial, we show how to perform optimization over the manifold of POVMs by using the QGOpt library. It is known that measurements of a qubit induced by tetrahedral POVM allow reconstructing an unknown qubit state with a minimal variance if there is no prior information about a qubit state. Let us check this fact numerically using optimization over the manifold of POVMs.\n",
    "In the beginning, let us import some libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "BE_lX6SwsgaF"
   },
   "outputs": [],
   "source": [
    "import tensorflow as tf  # tf 2.x\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "try:\n",
    "    import QGOpt as qgo\n",
    "except ImportError:\n",
    "    !pip install git+https://github.com/LuchnikovI/QGOpt@Dev\n",
    "    import QGOpt as qgo\n",
    "\n",
    "# Fix random seed to make results reproducable.\n",
    "tf.random.set_seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_7f8yOFGvrlJ"
   },
   "source": [
    "## 1. Prior information about a quantum state\n",
    "\n",
    "We represent a prior probability distribution over a quantum state approximately, by using a set of samples from a prior distribution. Since tetrahedral POVM is optimal when there is no prior information about a state, we consider uniform distribution across the Bloch ball."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "NmDWObAKNMR5"
   },
   "outputs": [],
   "source": [
    "#-------------------------------------------------------------------------#\n",
    "num_of_samples = 10000  # number of samples representing prior information\n",
    "#-------------------------------------------------------------------------#\n",
    "\n",
    "# Pauli matrices\n",
    "sigma_x = tf.constant([[0, 1], [1, 0]], dtype=tf.complex128)\n",
    "sigma_y = tf.constant([[0 + 0j, -1j], [1j, 0 + 0j]], dtype=tf.complex128)\n",
    "sigma_z = tf.constant([[1, 0], [0, -1]], dtype=tf.complex128)\n",
    "\n",
    "# All Pauli matrices in one tensor of shape (3, 2, 2)\n",
    "sigma = tf.concat([sigma_x[tf.newaxis],\n",
    "                   sigma_y[tf.newaxis],\n",
    "                   sigma_z[tf.newaxis]], axis=0)\n",
    "\n",
    "# Set of points distributed uniformly across Bloch ball\n",
    "x = tf.random.normal((num_of_samples, 3), dtype=tf.float64)\n",
    "x = x / tf.linalg.norm(x, axis=-1, keepdims=True)\n",
    "x = tf.cast(x, dtype=tf.complex128)\n",
    "u = tf.random.uniform((num_of_samples, 1), maxval=1, dtype=tf.float64)\n",
    "u = u ** (1 / 3)\n",
    "u = tf.cast(u, dtype=tf.complex128)\n",
    "x = x * u\n",
    "\n",
    "# Set of density matrices distributed uniformly across Bloch ball\n",
    "# (prior information)\n",
    "rho = 0.5 * (tf.eye(2, dtype=tf.complex128) + tf.tensordot(x, sigma, axes=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "McVrCrB62_dk"
   },
   "source": [
    "## 2. Search for the optimal POVM with given prior information about a state\n",
    "\n",
    "Here we search for the optimal POVM via minimizing the variance of a posterior distribution over density matrices. First, we define an example of the POVMs manifold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = qgo.manifolds.POVM()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The manifolds of POVMs is represented through the quadratic parametrization $M_i = A_iA_i^\\dagger$ with an equivalence relation $A_i\\sim A_iQ_i$, where $Q_i$ is an arbitrary unitary matrix. Here, we initialize a variable that represents the parametrization of each element of POVM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# randon initial parametrization of POVM\n",
    "A = m.random((4, 2, 2), dtype=tf.complex128)\n",
    "# real representtion of A\n",
    "A = qgo.manifolds.complex_to_real(A)\n",
    "# tf.Variable to be tuned\n",
    "A = tf.Variable(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we initialize Riemannian Adam optimizer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = 0.03\n",
    "opt = qgo.optimizers.RAdam(m, lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we ran the part of code that calculates forward pass, gradients, and optimization step several times until convergence to the optimal point is reached:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "xSk4-4x3NjVO"
   },
   "outputs": [],
   "source": [
    "for i in range(1000):\n",
    "    with tf.GradientTape() as tape:\n",
    "\n",
    "        # Complex representation of A\n",
    "        Ac = qgo.manifolds.real_to_complex(A)\n",
    "        # POVM from its parametrization\n",
    "        povm = Ac @ tf.linalg.adjoint(Ac)\n",
    "        # Inverce POVM (is needed to map a probability distribution to a density matrix)\n",
    "        povm_inv = tf.linalg.inv(tf.reshape(povm, (4, 4)))\n",
    "        # Matrix T maps probability vector to four real parameters representing\n",
    "        # a quantum state (equivalent to inverse POVM)\n",
    "        T = tf.concat([tf.math.real(povm_inv[0, tf.newaxis]),\n",
    "                        tf.math.real(povm_inv[3, tf.newaxis]),\n",
    "                        tf.math.real(povm_inv[2, tf.newaxis]),\n",
    "                        tf.math.imag(povm_inv[2, tf.newaxis])], axis=0)\n",
    "\n",
    "        # POVM maps a quantum state to a probability vector\n",
    "        p = tf.tensordot(rho, povm, axes=[[2], [1]])\n",
    "        p = tf.transpose(p, (0, 2, 1, 3))\n",
    "        p = tf.math.real(tf.linalg.trace(p))\n",
    "\n",
    "        # Covariance matrix of a reconstructed density matrix\n",
    "        cov = -p[:, tf.newaxis] * p[..., tf.newaxis]\n",
    "        cov = cov + tf.linalg.diag(p ** 2)\n",
    "        cov = cov + tf.linalg.diag(p * (1 - p))\n",
    "        cov = tf.tensordot(T, cov, [[1], [1]])\n",
    "        cov = tf.tensordot(cov, T, [[2], [1]])\n",
    "        cov = tf.transpose(cov, (1, 0, 2))\n",
    "\n",
    "        # Covariance matrix avaraged over prior distribution\n",
    "        av_cov = tf.reduce_mean(cov, axis=0)\n",
    "\n",
    "        # loss function (log volume of Covariance matrix)\n",
    "        loss = tf.reduce_sum(tf.math.log(tf.linalg.svd(av_cov)[0][:-1]))\n",
    "    grad = tape.gradient(loss, A)  # gradient\n",
    "    opt.apply_gradients(zip([grad], [A]))  # minimization step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "paAU_lip8ex5"
   },
   "source": [
    "## 3. Verification\n",
    "\n",
    "Here we check the resulting optimal POVM. For tetrahedral POVM one has the following relation ${\\rm Tr}\\left(M^\\alpha M^\\beta\\right) = \\frac{2\\delta_{\\alpha\\beta} + 1}{12}$. One can see, that this relation is almost true for a resulting POVM. The small error appears due to the approximate Monte-Carlo averaging of a covariance matric by using a set of samples from the prior distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 360
    },
    "colab_type": "code",
    "id": "8rPatx71bYPs",
    "outputId": "8b7042e3-e85a-48a9-9120-39bdb07b1fe3"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(\n",
      "[[0.24927765 0.08337808 0.08333673 0.08328467]\n",
      " [0.08337808 0.24939711 0.08328523 0.08333633]\n",
      " [0.08333673 0.08328523 0.25029829 0.08337795]\n",
      " [0.08328467 0.08333633 0.08337795 0.25102899]], shape=(4, 4), dtype=float64)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAECCAYAAADesWqHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAACItJREFUeJzt3c+LXfUdxvHn6TgmGgVpzSJMhuhCBHFhzJCN0EWwNA0BuzQLV8KshAjduPUfcNdNwNAWRBHiQkLaEEpEQjXmBzEYR2sQimOEpIpoWpqY8OlibktiB+ZMer73zJ3n/YKBucnlznOSvHPm3pnkuKoEIMtPhh4AYPwIHwhE+EAgwgcCET4QiPCBQBMdvu3dtj+1fdH2S0Pv6ZPtg7Yv2/5o6C0t2J61fdz2gu0LtvcPvakvtjfa/sD2h6Nje3noTT/mSf06vu0pSX+V9AtJi5JOSdpXVR8POqwntn8u6aqkP1TV40Pv6ZvtLZK2VNVZ2/dLOiPp1+vh98+2JW2qqqu2pyWdkLS/qt4feNp/TfIZf6eki1X1eVVdl/SGpGcG3tSbqnpX0jdD72ilqr6qqrOj97+XtCBpZthV/aglV0c3p0dva+oMO8nhz0j64pbbi1onf3DS2H5I0nZJJ4dd0h/bU7bPSbos6VhVraljm+TwvcyPram/VbEy2/dJOiTpxar6bug9famqm1X1hKStknbaXlNP1yY5/EVJs7fc3irp0kBbcAdGz38PSXqtqt4aek8LVfWtpHck7R54ym0mOfxTkh6x/bDtuyU9K+ntgTeho9ELYK9KWqiqV4be0yfbm20/MHr/HklPS/pk2FW3m9jwq+qGpBckHdXSC0NvVtWFYVf1x/brkt6T9KjtRdvPD72pZ09Jek7SLtvnRm97hh7Vky2Sjts+r6UT1LGqOjzwpttM7JfzANy5iT3jA7hzhA8EInwgEOEDgQgfCDTx4dueH3pDSxzfZFurxzfx4Utak7+wPeL4JtuaPL71ED6AVWryDTwP/nSqts3e1fvjLufK1ze1+WdTY/lY//HZ+U1j+1g/6JqmtWFsH2/cOL5+/Uv/0PW6ttw/YLtNkzq3zd6lv/xp/f4L2b0zO4aeACzrZP250/34VB8IRPhAIMIHAhE+EIjwgUCEDwQifCAQ4QOBCB8IRPhAIMIHAhE+EIjwgUCEDwQifCAQ4QOBCB8IRPhAIMIHAhE+EIjwgUCEDwTqFL7t3bY/tX3R9kutRwFoa8XwbU9J+q2kX0l6TNI+24+1HgagnS5n/J2SLlbV51V1XdIbkp5pOwtAS13Cn5H0xS23F0c/BmBCdQl/uetw/c8F92zP2z5t+/SVr2/+/8sANNMl/EVJs7fc3irp0o/vVFUHqmququbGfRFLAKvTJfxTkh6x/bDtuyU9K+nttrMAtLTi1XKr6obtFyQdlTQl6WBVXWi+DEAznS6TXVVHJB1pvAXAmPCde0AgwgcCET4QiPCBQIQPBCJ8IBDhA4EIHwhE+EAgwgcCET4QiPCBQIQPBCJ8IBDhA4EIHwhE+EAgwgcCET4QiPCBQIQPBCJ8IFCn/157tT47v0l7Z3a0eOg14fCXZ4ae0NR6/r3DEs74QCDCBwIRPhCI8IFAhA8EInwgEOEDgQgfCET4QCDCBwIRPhCI8IFAhA8EInwgEOEDgQgfCET4QCDCBwIRPhCI8IFAhA8EInwgEOEDgQgfCLRi+LYP2r5s+6NxDALQXpcz/u8k7W68A8AYrRh+Vb0r6ZsxbAEwJjzHBwL1dtFM2/OS5iVpo+7t62EBNNDbGb+qDlTVXFXNTWtDXw8LoAE+1QcCdfly3uuS3pP0qO1F28+3nwWgpRWf41fVvnEMATA+fKoPBCJ8IBDhA4EIHwhE+EAgwgcCET4QiPCBQIQPBCJ8IBDhA4EIHwhE+EAgwgcCET4QiPCBQIQPBCJ8IBDhA4EIHwhE+EAgwgcC9XYJrSR7Z3YMPaGpI1+eHXpCU3tmnhx6wuA44wOBCB8IRPhAIMIHAhE+EIjwgUCEDwQifCAQ4QOBCB8IRPhAIMIHAhE+EIjwgUCEDwQifCAQ4QOBCB8IRPhAIMIHAhE+EIjwgUCEDwQifCDQiuHbnrV93PaC7Qu2949jGIB2ulxJ54ak31TVWdv3Szpj+1hVfdx4G4BGVjzjV9VXVXV29P73khYkzbQeBqCdVT3Ht/2QpO2STrYYA2A8Ol800/Z9kg5JerGqvlvm5+clzUvSRt3b20AA/et0xrc9raXoX6uqt5a7T1UdqKq5qpqb1oY+NwLoWZdX9S3pVUkLVfVK+0kAWutyxn9K0nOSdtk+N3rb03gXgIZWfI5fVSckeQxbAIwJ37kHBCJ8IBDhA4EIHwhE+EAgwgcCET4QiPCBQIQPBCJ8IBDhA4EIHwhE+EAgwgcCET4QiPCBQIQPBCJ8IBDhA4EIHwhE+EAgwgcCdb6EFnLsmXly6AlNHb10bugJzez85T873Y8zPhCI8IFAhA8EInwgEOEDgQgfCET4QCDCBwIRPhCI8IFAhA8EInwgEOEDgQgfCET4QCDCBwIRPhCI8IFAhA8EInwgEOEDgQgfCET4QCDCBwKtGL7tjbY/sP2h7Qu2Xx7HMADtdLmSzjVJu6rqqu1pSSds/7Gq3m+8DUAjK4ZfVSXp6ujm9OitWo4C0Fan5/i2p2yfk3RZ0rGqOrnMfeZtn7Z9+gdd63sngB51Cr+qblbVE5K2Stpp+/Fl7nOgquaqam5aG/reCaBHq3pVv6q+lfSOpN1N1gAYiy6v6m+2/cDo/XskPS3pk9bDALTT5VX9LZJ+b3tKS39RvFlVh9vOAtBSl1f1z0vaPoYtAMaE79wDAhE+EIjwgUCEDwQifCAQ4QOBCB8IRPhAIMIHAhE+EIjwgUCEDwQifCAQ4QOBCB8IRPhAIMIHAhE+EIjwgUCEDwQifCAQ4QOBvHRNzJ4f1L4i6W+9P/DyHpT09zF9rCFwfJNt3Me3rao2r3SnJuGPk+3TVTU39I5WOL7JtlaPj0/1gUCEDwRaD+EfGHpAYxzfZFuTxzfxz/EBrN56OOMDWCXCBwIRPhCI8IFAhA8E+jefFK6oHytQyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cross = tf.tensordot(povm, povm, [[2], [1]])\n",
    "cross = tf.transpose(cross, (0, 2, 1, 3))\n",
    "cross = tf.linalg.trace(cross)\n",
    "cross = tf.math.real(cross)\n",
    "plt.matshow(cross)\n",
    "print(cross)"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "optimal_povm.ipynb",
   "provenance": []
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}