This notebook is the computational appendix of arXiv:1612.08551. We demonstrate how to use a modified version of the NPA hierarchy to test whether a set of observed correlations is local. We also give the details how to reproduce the numerical results shown in the paper. Furthermore, we show how to determine the robustness to noise of a set of given nonlocal correlations and to extract a Bell inequality from the dual of the SDP. To improve readability of this notebook, we placed the supporting functions to a separate file; please download this in the same folder as the notebook if you would like to evaluate it. The following dependencies must also be available: at least one SDP solver (SDPA as an executable in the path or Mosek with its Python interface installed; cvxopt as a solver is not recommended) together with the Ncpol2sdpa and qutip packages. First, we import everything we will need:
In [ ]:
from __future__ import print_function, division
from math import sqrt
from qutip import sigmax, sigmaz
from ncpol2sdpa import flatten, SdpRelaxation, generate_variables
from time import time
from sympy import S
from local_tools import generate_commuting_measurements, get_W_reduced, \
get_GHZ_reduced, get_moment_constraints, get_fullmeasurement
Then we define the scenario we are considering. In full generality we define it with the three parameters $(N,m,d)$, corresponding to the case of $N$ parties, $m$ measurement choices with $d$ outcomes each. For instance
In [ ]:
N, m, d = 7, 2, 2
configuration = [d for _ in range(m)]
measurements = [generate_commuting_measurements(configuration, chr(65+i))
for i in range(N)]
We generated the symbolic variables defining the measurements. They are treated as commuting variables, since to test for locality, we use the NPA hierarchy with commuting measurements. Given that we will always work in the correlator space, we define the substitution rule $\mathcal({M}_{x_i}^{(i)})^2 = \mathbb{1}$ for any $i = 1,...,N$ and $x_i = 0,..,m-1$.
In [ ]:
substitutions = {M**2: S.One for M in flatten(measurements)}
In [ ]:
W_state = get_W_reduced(N)
W_operators = [[sigmax(), sigmaz()] for _ in range(N)]
Then we generate the moments, i.e. the values of the correlators:
In [ ]:
time0 = time()
moments = get_moment_constraints(N, W_state, measurements, W_operators)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
We construct the hierarchy at level two with corresponding substitutions and moments assignments:
In [ ]:
time0 = time()
sdp = SdpRelaxation(flatten(measurements))
sdp.get_relaxation(2, substitutions=substitutions, momentsubstitutions=moments)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Finally, we solve the SDP feasibility problem:
In [ ]:
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("SDP status is " + sdp.status)
In [ ]:
GHZ_state = get_GHZ_reduced(4)
GHZ_operators = [[sigmax(), 1/sqrt(2)*(sigmaz()+sigmax())] for _ in range(N)]
For the $GHZ$ state, we need to add the values of the two full-body correlators $\langle M_0^{(1)}M_0^{(2)}...M_0^{(N)} \rangle$ and $\langle M_1^{(1)}M_0^{(2)}...M_0^{(N)} \rangle$.
In [ ]:
time0 = time()
moments = get_moment_constraints(N, GHZ_state, measurements, GHZ_operators)
extra = get_fullmeasurement([0 for _ in range(N)], measurements)
extramonomials = [extra]
moments[extra] = 1
extra = get_fullmeasurement(flatten([1,[0 for _ in range(N-1)]]), measurements)
extramonomials.append(extra)
moments[extra] = (1/sqrt(2))
print("Constraints were generated in " + str(time()-time0) + " seconds.")
We construct the hierarchy at the hybrid level $\lbrace \mathcal{O}_2,M_0^{(1)}M_0^{(2)}...M_0^{(N)},M_1^{(1)}M_0^{(2)}...M_0^{(N)} \rbrace$ with corresponding substitutions and moments assignments:
In [ ]:
time0 = time()
sdp = SdpRelaxation(flatten(measurements), verbose=1, parallel=True)
sdp.get_relaxation(2, substitutions=substitutions,
momentsubstitutions=moments, extramonomials=extramonomials)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
We solve the SDP feasibility problem
In [ ]:
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("SDP status is " + sdp.status)
We show how to estimate the noise robustness of some observed nonlocal correlations numerically. Given a set of correlations generated by a state $\rho$ and some measurements $M_{x_i}^{(i)}$, we estimate the minimal amount of noise $\lambda$ such that the correlations arise from $(1 - \lambda )\rho + \lambda \frac{\mathbb{1}}{2^N}$. If we perform measurements on a Bloch sphere, the noise affects the values of the correlators the same way. That is, if we denote $v_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4}$ the value of $\langle M_{j_1}^{(i_1)}M_{j_2}^{(i_2)}M_{j_3}^{(i_3)}M_{j_4}^{(i_4)} \rangle$, the noisy state will result in the corresponding value of $(1 - \lambda ) v_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4} $ We cast this problem as an SDP by introducing $\lambda$ as the objective function of the minimization problem and by replacing the moments substitutions with the modified values depending on $\lambda$. First, we define $\lambda$ as an additional symbolic variable:
In [ ]:
lambda_ = generate_variables("\lambda")[0]
In [ ]:
time0 = time()
moments = get_moment_constraints(N, W_state, measurements, W_operators, lambda_)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
Then the relaxation:
In [ ]:
sdp = SdpRelaxation(flatten(measurements), parameters=[lambda_],
verbose=1, parallel=True)
sdp.get_relaxation(2, objective=lambda_, substitutions=substitutions,
momentsubstitutions=moments)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
and finally we solve the SDP:
In [ ]:
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("lambda_min is " + str(sdp.primal))
In [ ]:
time0 = time()
moments = get_moment_constraints(N, GHZ_state, measurements, GHZ_operators, lambda_)
extra = get_fullmeasurement([0 for _ in range(N)], measurements)
extramonomials = [extra]
moments[extra] = (1 - lambda_)
extra = get_fullmeasurement(flatten([1, [0 for _ in range(N-1)]]), measurements)
extramonomials.append(extra)
moments[extra] = (1/sqrt(2))*(1 - lambda_)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
The relaxation is generated as:
In [ ]:
sdp = SdpRelaxation(flatten(measurements), parameters=[lambda_],verbose=1, parallel=True)
sdp.get_relaxation(2, objective=lambda_, substitutions=substitutions,
momentsubstitutions=moments, extramonomials=extramonomials)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Solving it, we get:
In [ ]:
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("lambda_min is " + str(sdp.primal))
From the dual of the noise robustness SDP, one can extract a Bell inequality that detects the observed nonlocal correlations. The derived inequality will be in the form $\sum \alpha_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4} \langle M_{j_1}^{(i_1)}M_{j_2}^{(i_2)}M_{j_3}^{(i_3)}M_{j_4}^{(i_4)} \rangle + 1 \geq 0$. We will present the code we used for the case of the $GHZ$ and $m = 2$ measurement choices. First, we solve the noise robustness SDP:
In [ ]:
time0 = time()
moments = get_moment_constraints(N, GHZ_state,measurements, GHZ_operators,
lambda_, 2)
extra = get_fullmeasurement([0 for _ in range(N)], measurements)
extramonomials = [extra]
moments[extra] = (1 - lambda_)
extra = get_fullmeasurement(flatten([1, [0 for _ in range(N-1)]]), measurements)
extramonomials.append(extra)
moments[extra] = (1/sqrt(2))*(1 - lambda_)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
In [ ]:
sdp = SdpRelaxation(flatten(measurements), parameters=[lambda_], verbose=1, parallel=True)
sdp.get_relaxation(2, objective=lambda_, substitutions=substitutions,
momentsubstitutions=moments, extramonomials=extramonomials)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("lambda_min is " + str(sdp.primal))
We only assign the values of up to the two-body correlators. This will ensure that we get the same expression for the Bell inequality as the one presented in the manuscript. Now, we proceed to generate a second relaxation, but we do not substitute the monomials corresponding to the observed correlations.
In [ ]:
time0 = time()
sdp2 = SdpRelaxation(flatten(measurements), verbose=1, parallel=True)
sdp2.get_relaxation(2, substitutions=substitutions, extramonomials=extramonomials)
print("Second SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Then we assign to the second SDP the same dual variables as for the one solved before.
In [ ]:
sdp2.status = sdp.status
sdp2.y_mat = [sdp.y_mat[1]]
We extract the value of the dual variable for each monomial appearing in the moment matrix. For the function "extract_dual_value" to work properly we needed to generate a moment matrix in which the values of the correlators where not substituted yet.
In [ ]:
time0 = time()
bound = sdp2.extract_dual_value(1)
ineq = 0
for monomial, index in sdp2.monomial_index.items():
ineq += round(2*sdp2.extract_dual_value(monomial),4)*monomial
print("Dual was generated in " + str(time()-time0) + " seconds")
Lastly, we normalize the inequality so to be in the form presented before (i.e. with the classical bound being equal to $1$ )
In [ ]:
print(ineq/bound)