ThreeStarFNNInequalities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computational appendix of [Nat. Commun. 14, 2153 (2023)](https://www.nature.com/articles/s41467-023-37842-w)\n",
    "---\n",
    "\n",
    "In this notebook we perform a step-by-step implementation of the hybrid classical-nonsignaling inflation of the three-branch scenario presented in [Nat. Commun. 14, 2153 (2023)](https://www.nature.com/articles/s41467-023-37842-w) [[arXiv:2212.09765](https://www.arxiv.org/abs/2212.09765)] and obtain the full network nonlocality witnesses presented there, namely Equation (4).\n",
    "\n",
    "Authors: Alejandro Pozas-Kerstjens\n",
    "\n",
    "Requires: [mosek](https://www.mosek.com/) >= 9.1.13 & <= 9.2.28, [numpy](https://numpy.org/), [picos](https://picos-api.gitlab.io/picos/) <= 2.2.55, [qutip](http://qutip.org/), [sympy](https://www.sympy.org/)\n",
    "\n",
    "Last updated: 21 Dec, 2022\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this notebook we follow the first part in the Materials and Methods of the manuscript in order to implement the linear program arising from application of the inflation technique to detect whether a four-partite distribution admits a realisation in the three-branch star network when one of the sources is assumed to be classical and the other ones distribute general no-signaling resources. We will implement, step by step, the linear program defined there through direct implementation of Equations (6-10), solving the problem, extracting the associated certificate of infeasibility, and deriving from it the full network nonlocality witness in Equation (4) in the manuscript.\n",
    "\n",
    "We begin by importing the necessary libraries and by defining several objects that will be useful throughout the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import qutip as qt\n",
    "import sympy as sp\n",
    "import picos as pic\n",
    "\n",
    "from sympy import Symbol as S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "oA=2; mA=2; oB=2;\n",
    "\n",
    "l = r'\\langle{}'; r = r'\\rangle{}';    # For pretty visuals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define the family of probability distributions that we consider. These are generated by measuring a (potentially noisy) maximally entangled state between each of the branch parties and the central party, the latter performing a two-outcome measurement described by the projectors $|\\text{GHZ}\\rangle\\langle\\text{GHZ}|$ and $1-|\\text{GHZ}\\rangle\\langle\\text{GHZ}|$, and the formers being able to perform the same two binary-outcome measurements defined on the X-Z plane, characterized by angles $\\theta_0$ and $\\theta_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prob(v, θ0, θ1):\n",
    "    ϕplus = qt.ket2dm(qt.bell_state('00'))\n",
    "    ϕplus.dims = [[4], [4]]\n",
    "    ghzplus = qt.ket2dm((qt.basis([2, 2, 2], [0, 0, 0])\n",
    "                         + qt.basis([2, 2, 2], [1, 1, 1]))\n",
    "                        / np.sqrt(2))\n",
    "    ghzminus = qt.ket2dm((qt.basis([2, 2, 2], [0, 0, 0])\n",
    "                          - qt.basis([2, 2, 2], [1, 1, 1]))\n",
    "                         / np.sqrt(2))\n",
    "    id2 = qt.qeye(2); id4 = qt.qeye(4);\n",
    "\n",
    "    ρ = qt.tensor([v*ϕplus+(1-v)*id4/4,\n",
    "                   v*ϕplus+(1-v)*id4/4,\n",
    "                   v*ϕplus+(1-v)*id4/4])\n",
    "    ρ.dims = [[2, 2, 2, 2, 2, 2], [2, 2, 2, 2, 2, 2]]\n",
    "    ρ = ρ.permute([0,2,4,1,3,5])\n",
    "\n",
    "    A = [np.sin(θ0)*qt.sigmax() + np.cos(θ0)*qt.sigmaz(),\n",
    "         np.sin(θ1)*qt.sigmax() + np.cos(θ1)*qt.sigmaz()]\n",
    "    B = [ghzplus, qt.qeye([2, 2, 2]) - ghzplus]\n",
    "\n",
    "    p = np.zeros((oA,oA,oA,oB,mA,mA,mA))\n",
    "    for a1,a2,a3,b,x1,x2,x3 in np.ndindex(oA,oA,oA,oB,mA,mA,mA):\n",
    "        p[a1,a2,a3,b,x1,x2,x3] = np.real(ρ.overlap(\n",
    "            qt.tensor([id2+(-1)**a1*A[x1], id2+(-1)**a2*A[x2],\n",
    "                       id2+(-1)**a3*A[x3], B[b]]))/8)\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define the function that will set up the linear program corresponding to compatibility with the inflation in Fig. 4 in the manuscript, where the information sent from the source connecting the central party with $A^{(1)}$ is assumed to be encoded in a classical system, and thus is copied and sent to multiple copies of $A^{(1)}$. We will do so using the library [PICOS](https://picos-api.gitlab.io/picos/), which will allow us to define the probability distribution $p_{\\text{inf}}$ subject to the normalization, no-signaling and marginalization constraints, as well as to symmetry constraints derived from the fact that there exist multiple copies of the party $A^{(1)}$. Within PICOS we will also extract the certificate of infeasibility in case the linear program does not admit any solution, and later transform this into a Bell-like inequality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_problem(prob_dist, solve=False):\n",
    "    oA = prob_dist.shape[0]\n",
    "    oB = prob_dist.shape[3]\n",
    "    mA = prob_dist.shape[4]\n",
    "    probabilities = pic.RealVariable('pinf',\n",
    "                                     shape=oA*oA*oA*oA*oA*oB*mA*mA*mA*mA*mA,\n",
    "                                     lower=0) # Directly impose positivity (Eq. 6) here\n",
    "    pinf = np.reshape(probabilities,\n",
    "                      (oA,oA,oA,oA,oA,oB,mA,mA,mA,mA,mA))\n",
    "\n",
    "    # Normalisation (Eq. 7)\n",
    "    norm_const = [probsum == 1\n",
    "                  for probsum in np.sum(pinf, (0, 1, 2, 3, 4, 5)\n",
    "                                       ).flatten().tolist()]\n",
    "\n",
    "    # No-signaling for party A11 (Eq. 8)\n",
    "    nsA_const = [[] for _ in range(mA-1)]\n",
    "    nsA_0_terms = np.sum(pinf[:,:,:,:,:,:,0,:,:,:,:], 0)\n",
    "    for x in range(1, mA):\n",
    "        nsA_x_terms = np.sum(pinf[:,:,:,:,:,:,x,:,:,:,:], 0)\n",
    "        nsA_const[x-1] = [term == 0 for term in\n",
    "                          (nsA_0_terms - nsA_x_terms).flatten().tolist()]\n",
    "    nsA1_const = sum(nsA_const, [])\n",
    "    # We do not need to impose no-signaling on A12 and A13\n",
    "    # because they will be enforced when setting the \n",
    "    # inflation symmetries\n",
    "\n",
    "    # No-signaling for party A2 (Eq. 8)\n",
    "    nsA_const = [[] for _ in range(mA-1)]\n",
    "    nsA_0_terms = np.sum(pinf[:,:,:,:,:,:,:,:,:,0,:], 3)\n",
    "    for x in range(1, mA):\n",
    "        nsA_x_terms = np.sum(pinf[:,:,:,:,:,:,:,:,:,x,:], 3)\n",
    "        nsA_const[x-1] = [term == 0 for term in\n",
    "                          (nsA_0_terms - nsA_x_terms).flatten().tolist()]\n",
    "    nsA2_const = sum(nsA_const, [])\n",
    "    \n",
    "    # No-signaling for party A3 (Eq. 8)\n",
    "    nsA_const = [[] for _ in range(mA-1)]\n",
    "    nsA_0_terms = np.sum(pinf[:,:,:,:,:,:,:,:,:,:,0], 4)\n",
    "    for x in range(1, mA):\n",
    "        nsA_x_terms = np.sum(pinf[:,:,:,:,:,:,:,:,:,:,x], 4)\n",
    "        nsA_const[x-1] = [term == 0 for term in\n",
    "                          (nsA_0_terms - nsA_x_terms).flatten().tolist()]\n",
    "    nsA3_const = sum(nsA_const, [])\n",
    "\n",
    "    # Inflation (Eq. 9). A11 <-> A12\n",
    "    inf_const = [term == 0 for term in\n",
    "                 (pinf-np.transpose(pinf, [1,0,2,3,4,5,7,6,8,9,10])\n",
    "                  ).flatten().tolist()]\n",
    "    # Inflation (Eq. 9). A11 <-> A13\n",
    "    inf_const += [term == 0 for term in\n",
    "                  (pinf-np.transpose(pinf, [2,1,0,3,4,5,8,7,6,9,10])\n",
    "                   ).flatten().tolist()]\n",
    "    # Inflation (Eq. 9). A12 <-> A13\n",
    "    inf_const += [term == 0 for term in\n",
    "                  (pinf-np.transpose(pinf, [0,2,1,3,4,5,6,8,7,9,10])\n",
    "                   ).flatten().tolist()]\n",
    "    # Inflation (Eq. 9). A11 -> A12 -> A13 -> A11\n",
    "    inf_const += [term == 0 for term in\n",
    "                  (pinf-np.transpose(pinf, [1,2,0,3,4,5,7,8,6,9,10])\n",
    "                   ).flatten().tolist()]\n",
    "    # Inflation (Eq. 9). A11 -> A13 -> A12 -> A11\n",
    "    inf_const += [term == 0 for term in\n",
    "                  (pinf-np.transpose(pinf, [2,0,1,3,4,5,8,6,7,9,10])\n",
    "                   ).flatten().tolist()]\n",
    "\n",
    "    # Matching original probabilities (Eq. 10)\n",
    "    marginal_const = [term == 0 for term in\n",
    "                      (np.sum(pinf[:,:,:,:,:,:,:,0,0,:,:], (1, 2))\n",
    "                       -prob_dist).flatten().tolist()]\n",
    "    \n",
    "    # Set up the problem\n",
    "    problem = pic.Problem()\n",
    "    for const in (norm_const + nsA1_const + nsA2_const\n",
    "                  + nsA3_const + inf_const + marginal_const):\n",
    "        problem.add_constraint(const)\n",
    "\n",
    "    # Solve if requested\n",
    "    if solve:\n",
    "        try:\n",
    "            problem.solve(solver='mosek')\n",
    "            print('Optimal solution found')\n",
    "            return [np.reshape(probabilities.value,\n",
    "                               (oA,oA,oA,oA,oA,oB,mA,mA,mA,mA,mA)),\n",
    "                    [None, None]]\n",
    "        except pic.SolutionFailure:\n",
    "            print('The problem is infeasible. Getting certificate...')\n",
    "            problem.solve(solver='mosek', primals=False)            \n",
    "            print('Primal or dual infeasibility certificate found')\n",
    "            ynorm     = np.array([const.dual for const in norm_const])\n",
    "            ymarginal = np.array([const.dual for const in marginal_const])\n",
    "            return [None, [ynorm, ymarginal]]\n",
    "    else:\n",
    "        return problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The problem is infeasible. Getting certificate...\n",
      "Primal or dual infeasibility certificate found\n"
     ]
    }
   ],
   "source": [
    "[pinf, [ynorm, ymarginal]] = set_problem(prob(0.882, -1.865, -0.4146),\n",
    "                                         solve=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The certificate of infeasibility is a vector $\\mathbf{y}$ that satisfies $\\mathbf{y}\\cdot A = \\mathbf{0}$ and $\\mathbf{y}\\cdot\\mathbf{b} > 0$, which signals the impossibility of finding a vector $\\mathbf{x}$ that satisfies $A\\cdot\\mathbf{x} \\geq \\mathbf{b}$. Importantly, since in the problem we are dealing with the matrix $A$ is constant regardless the distribution that we want to study, the quantity $\\mathbf{y}\\cdot\\mathbf{b} > 0$ can be understood as a Bell-like inequality that applies not only to the distribution under scrutiny but also to arbitrary ones. In the following, we write the inequality in terms of arbitrary probabilities, to read it as a Bell-like inequality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build symbols for the probabilities\n",
    "symp = np.zeros((oA,oA,oA,oB,mA,mA,mA), dtype=object)\n",
    "for a1,a2,a3,b,x1,x2,x3 in np.ndindex(oA,oA,oA,oB,mA,mA,mA):\n",
    "    symp[a1,a2,a3,b,x1,x2,x3] = sp.symbols(\n",
    "        f\"p({a1}{a2}{a3}{b}|{x1}{x2}{x3})\")\n",
    "\n",
    "# Function to write probabilities in terms of expectation values\n",
    "def prob_to_expect(symbol):\n",
    "    a1, a2, a3, b, _, x1, x2, x3 = str(symbol)[2:-1]\n",
    "    a1 = int(a1); a2 = int(a2); a3 = int(a3); b = int(b)\n",
    "    \n",
    "    expec = 1/16 * (\n",
    "        1\n",
    "        + (-1)**b * S(l + \"B\" + r)\n",
    "        + (-1)**a1 * S(l + \"A^{(1)}_\" + x1 + r)\n",
    "        + (-1)**a2 * S(l + \"A^{(2)}_\" + x2 + r)\n",
    "        + (-1)**a3 * S(l + \"A^{(3)}_\" + x3 + r)\n",
    "        + (-1)**(a1 + b) * S(l + \"A^{(1)}_\" + x1 + \"B\" + r)\n",
    "        + (-1)**(a2 + b) * S(l + \"A^{(2)}_\" + x2 + \"B\" + r)\n",
    "        + (-1)**(a3 + b) * S(l + \"A^{(3)}_\" + x3 + \"B\" + r)\n",
    "        + (-1)**(a1 + a2) * S(l + \"A^{(1)}_\" + x1 + \"A^{(2)}_\" + x2 + r)\n",
    "        + (-1)**(a1 + a3) * S(l + \"A^{(1)}_\" + x1 + \"A^{(3)}_\" + x3 + r)\n",
    "        + (-1)**(a2 + a3) * S(l + \"A^{(2)}_\" + x2 + \"A^{(3)}_\" + x3 + r)\n",
    "        + (-1)**(a1 + a2 + b) * S(l + \"A^{(1)}_\" + x1 + \"A^{(2)}_\" + x2 + \"B\" + r)\n",
    "        + (-1)**(a1 + a3 + b) * S(l + \"A^{(1)}_\" + x1 + \"A^{(3)}_\" + x3 + \"B\" + r)\n",
    "        + (-1)**(a2 + a3 + b) * S(l + \"A^{(2)}_\" + x2 + \"A^{(3)}_\" + x3 + \"B\" + r)\n",
    "        + (-1)**(a1 + a2 + a3) * S(l + \"A^{(1)}_\" + x1 + \"A^{(2)}_\" + x2 + \"A^{(3)}_\" + x3 + r)\n",
    "        + (-1)**(a1 + a2 + a3 + b) * S(l + \"A^{(1)}_\" + x1 + \"A^{(2)}_\" + x2 + \"A^{(3)}_\" + x3 + \"B\" + r)\n",
    "    )\n",
    "    return expec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 1.0 \\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_0B\\rangle{} - 1.0 \\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_0\\rangle{} - 1.0 \\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_1B\\rangle{} - 1.0 \\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_1\\rangle{} - 1.0 \\langle{}A^{(1)}_0A^{(3)}_0B\\rangle{} - 1.0 \\langle{}A^{(1)}_0A^{(3)}_0\\rangle{} - 1.0 \\langle{}A^{(1)}_0A^{(3)}_1B\\rangle{} - 1.0 \\langle{}A^{(1)}_0A^{(3)}_1\\rangle{} - 1.0 \\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_0B\\rangle{} - 1.0 \\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_0\\rangle{} + \\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_1B\\rangle{} + \\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_1\\rangle{} - 1.0 \\langle{}A^{(1)}_1A^{(3)}_0B\\rangle{} - 1.0 \\langle{}A^{(1)}_1A^{(3)}_0\\rangle{} + \\langle{}A^{(1)}_1A^{(3)}_1B\\rangle{} + \\langle{}A^{(1)}_1A^{(3)}_1\\rangle{} - 2.0 \\langle{}A^{(2)}_0B\\rangle{} - 2.0 \\langle{}A^{(2)}_0\\rangle{} - 2.0 \\langle{}B\\rangle{} - 2.0$"
      ],
      "text/plain": [
       "-1.0*\\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_0B\\rangle{} - 1.0*\\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_0\\rangle{} - 1.0*\\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_1B\\rangle{} - 1.0*\\langle{}A^{(1)}_0A^{(2)}_0A^{(3)}_1\\rangle{} - 1.0*\\langle{}A^{(1)}_0A^{(3)}_0B\\rangle{} - 1.0*\\langle{}A^{(1)}_0A^{(3)}_0\\rangle{} - 1.0*\\langle{}A^{(1)}_0A^{(3)}_1B\\rangle{} - 1.0*\\langle{}A^{(1)}_0A^{(3)}_1\\rangle{} - 1.0*\\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_0B\\rangle{} - 1.0*\\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_0\\rangle{} + \\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_1B\\rangle{} + \\langle{}A^{(1)}_1A^{(2)}_0A^{(3)}_1\\rangle{} - 1.0*\\langle{}A^{(1)}_1A^{(3)}_0B\\rangle{} - 1.0*\\langle{}A^{(1)}_1A^{(3)}_0\\rangle{} + \\langle{}A^{(1)}_1A^{(3)}_1B\\rangle{} + \\langle{}A^{(1)}_1A^{(3)}_1\\rangle{} - 2.0*\\langle{}A^{(2)}_0B\\rangle{} - 2.0*\\langle{}A^{(2)}_0\\rangle{} - 2.0*\\langle{}B\\rangle{} - 2.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Extract certificate in terms of probabilities\n",
    "I1 = ymarginal @ symp.flatten() + sum(ynorm)\n",
    "\n",
    "# Write in terms of correlators\n",
    "16*I1.subs({prob: prob_to_expect(prob) for prob in I1.free_symbols})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This quantity, $\\mathcal{I}_1$, is negative for all distributions that can be generated in the three-branch star network when the source connecting $A^{(1)}$ to the central party is classical, and thus obtaining a positive value signals that the distribution that achieves this value cannot be produced using a classical source there. Because the star network considered is symmetric under cyclic permutations of the branch parties, $\\mathcal{I}_2$ and $\\mathcal{I}_3$, corresponding to the situations where the classical sources are the one connecting, respectively, $A^{(2)}$ and $A^{(3)}$ with $B$, can be obtained by performing the permutations $A^{(1)}\\rightarrow A^{(2)}\\rightarrow A^{(3)} \\rightarrow A^{(1)}$ and $A^{(1)}\\rightarrow A^{(3)}\\rightarrow A^{(2)} \\rightarrow A^{(1)}$. A simultaneous violation of $\\mathcal{I}_1\\leq 0$, $\\mathcal{I}_2\\leq 0$ and $\\mathcal{I}_3\\leq 0$ by a probability distribution implies that the distribution is fully network nonlocal."
   ]
  }
 ],
 "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}