## SDP interface and Solver

PICOS is the Python interface I recommend to write/solve LP/SDP problem. It is under the GPLv3 free-license, probably supports your favorite solver, has crucial features for quantum information and, most importantly, has an attentive and friendly developer community.

Regarding SDP solvers, I had best performance with MOSEK. Sadly, this solver is proprietary (if price is of your concern, MOSEK offers a free one-year license for students). CVXOPT is a relatively good free-licensed alternative to MOSEK. However, in my experience, CVXOPT might lead to some issues (memory-leakage, unstable with complex SDP), and an important optimization time compare to MOSEK.

## LP: Local correlations, a membership problem

For this section we need the following import

import picos as pcs
import numpy as np
from itertools import product


Consider a scenario where a state is shared between N parties, each having m measurement choices (inputs) with d possible results (outputs). The correlations between the parties admit a local hidden-variable model iff the correlations belong to a local polytope [1,2].
To know if correlations of a given scenario are local, it is sufficient to check if these correlations belong, i.e. are “inside”, the local polytope.

Correlations are in a vector of probability $p = {p(a_1 \dots a_N|x_1 \dots x_N)}$ for all inputs x_i and output a_i. The vertices of the local polytope are the deterministic strategies p_d, i.e. p is composed of 0s and 1s. The local polytope can be written as the set of its vertices L=p_d.
Here a simple function that construct such a set

def local_polytope(N,m,d):
"""Local Polytope generation.

Generate :math:d^{m.N} vertices of the Local Polytpe for a system with N partite, each having m measurement/input and d output.

Parameters
----------
N : int
number of partite
m : int
number of input of one partite
d : int
number of output of one partite

Returns
-------
vertices : array_like
Polytope. Each row is a vertex (corresponding to a deterministic strategy), each column is the :math:n^{th} element to that vertice

"""
D = np.zeros((d**m,m), dtype=int)
for _ in range(d**m):
D[_][:] = np.array(np.unravel_index(_,(d,)*m))
vertices = np.zeros(((d**(m*N),)+(m,)*N+(d,)*N))
c = 0
for _ in product(range(d**m), repeat=N):
for x in product(range(m), repeat=N):
vertices[(c,)+x+tuple([D[_[i]][x[i]] for i in range(N)])]=1
c += 1
shape = np.prod(np.delete(vertices.shape[:],0,0))
polytope = vertices.reshape(vertices.shape[0],shape)
return polytope


A correlation vector p is described by a local model iff p can be written as a convex sum of the local-polytope vertices. I.e. iff p is inside the local polytope. This is a linear programming problem: find a vector v such that v*L=p, with sum(v)=1 (convexity) and v>0 (each element of v positive). If v exists then p is compatible with a local model.

Using PICOS we write this linear program

def is_local(correlation,polytope):
pb = pcs.Problem() # Create a picos Problem instance
# The two constants of our problem
L = pcs.Constant("L", polytope)     # local polytope
p = pcs.Constant("p", correlation)  # correlation vector
# The variable v (the coefficient of the convex sum)
dim = L.shape[0]                    # number of vertices
v = pcs.RealVariable("v", dim)      # v is a real vector with dim element
# Constraints
pb.add_constraint(v >= 0)           # all elements of v are positive
pb.add_constraint(sum(v) == 1)      # sum of elements of v is 1
pb.add_constraint(L.T*v == p)       # v is the coefficient of the convex sum of all L vertices
# Objective
pb.set_objective("min", 1 | v)      # dummy objective

# Solving the LP
pb.solve(solver="cvxopt",primals=False)
if pb.status == "optimal":
return True                     # if an optimal solution exists, p belongs to L
else:
return False                    # otherwise p is non-local


To test our LP implementation let’s now create a function that returns a local correlation vector p, by construction.

def local_correlation(polytope):
nb_vertices = len(polytope) # Number of vertices composing the polytope
coeff = np.random.random(nb_vertices) # Random coefficients, one per vertices
coeff /= np.sum(coeff) # Normalised coefficients
p = np.sum([c*v for (c,v) in zip(coeff,polytope)],axis=0) # p is a convex sum of the polytopes vertices
return p


We verify our implementation with two examples, one local, one non-local

L = local_polytope(2,2,2)
p_l = local_correlation(L)
p_nl = L[1]+1e-3           # A non local correlation -> we take a vertice and we add 0.01 to each of its elements
print(is_local(p_l,L))     # Output True
print(is_local(p_nl,L))    # Output False


## (Complex) SDP: Maximize a Bell inequality

For this section, import and constant definition needed

import numpy as np
import picos as pcs

# Paulis
Z = np.array([[1.,0.],[0.,-1.]])
X = np.array([[0.,1.],[1.,0.]])
# Kroenecker product
K = np.kron


Consider a scenario where Alice and Bob share a two-qubit state rho. They each have two measurements : A0,A1 for Alice, B0,B1 for Bob. These measurements are of the form $M_i = \cos(\theta)\sigma_x + (-1)^i \sin(\theta)\sigma_z$ with i the measurement choice, 0 or 1, and θ the angle between measurement 0 and 1. Alice set θ=a while Bob set θ=b.
A bell operator can be seen as a linear combination of the four correlators AxBy. For example the well-known CHSH operator S=A0(B0+B1)+A1(B0-B1). Here is a function to generate such operators (return CHSH by default)

def bell_op(a=np.pi/4,b=np.pi/4,coeff=[1,1,1,-1]):
""" Return Bell operator, for specific angle between measurements (for each partite resp.)

Parameters
----------
a : float
Angle between observabe, Alice side.
b : float
Angle between observabe, Bob side.
coeff : array
Numpy array of four elements, for A0B0, A1B0, A0B1, A1B1.

Return
------
bell : ndarray
Bell operator
"""

O = lambda x,a : np.cos(a)*X + (-1)**x*np.sin(a)*Z
if coeff == [1,1,1,-1]:
#CHSH case
bell = K(O(0,a),O(0,b)+O(1,b)) + K(O(1,a),O(0,b)-O(1,b))
else:
bell = coeff[0]*K(O(0,a),O(0,b))+coeff[1]*K(O(1,a),O(0,b))+coeff[2]*K(O(0,a),O(1,b))+coeff[3]*K(O(1,a),O(1,b))
return bell


The CHSH score achieved by rho is given by the linear operation $\text{Tr}[S\rho]$. Since all two-qubits state are semi-positive definite and the CHSH value is linear (in rho) we can use semi-definite programming (SDP) to find the maximum CHSH score a two-qubit state can achieve. Such an SDP reads $\max_\rho \qquad \text{Tr}[S\rho]$ $\text{s.t.} \qquad \rho \succeq 0$ $\qquad \quad \text{Tr}[\rho]=1$ with rho a 4x4 Hermitian matrix.

Using PICOS, we can write this SDP as follow

def max_bell(operator):
""" Compute the max(Tr[operator*rho]) over all rho = two-qubits state

Parameters
----------
operator : np.ndarray
Two-qubits operator
"""
sdp = pcs.Problem()
#SDP "constant", the bell operator
bell_operator = pcs.Constant(operator)
#SDP variable (quantum state)
rho = pcs.HermitianVariable('rho',4)
#Constraints of the SDP :
#semi positive definite
#Trace 1
# Objective of the SDP: Tr[Bell*rho]
obj = pcs.trace(bell_operator*rho).real
sdp.set_objective('max',obj)

# Some print stuff
print("Solving: ")
print(sdp)
sol = sdp.solve(solver='cvxopt',verbosity=0)
print(f"Solved with {sol.claimedStatus} status")
print(f"Optimal violation: {obj}")
print(f"State reaching this violation:\n{rho}")


Let’s test our implementation

S = bell_op()
max_bell(S)


Running this output

Solving max(Tr[CHSH*rho] over all two-qubit states
Solving:
Complex Semidefinite Program
maximize Re(tr([4×4]·rho))
over
4×4 hermitian variable rho
subject to
rho ≽ 0
tr(rho) = 1
Solved with optimal status
Optimal violation: 2.8284271145023476
State reaching this violation:
[ 7.32e-02-j0.00e+00  1.77e-01-j1.85e-17  1.77e-01+j5.72e-18 -7.32e-02+j2.73e-17]
[ 1.77e-01+j1.85e-17  4.27e-01-j0.00e+00  4.27e-01+j5.60e-17 -1.77e-01+j4.77e-17]
[ 1.77e-01-j5.72e-18  4.27e-01-j5.60e-17  4.27e-01-j0.00e+00 -1.77e-01+j7.19e-17]
[-7.32e-02-j2.73e-17 -1.77e-01-j4.77e-17 -1.77e-01-j7.19e-17  7.32e-02-j0.00e+00]


We obtain the expected violation of 2.8284... which corresponds to 2√2 achieved by Bell states [2].

## References

[1] Pironio S., J. Phys. A: Math. theor. 47 (424020)
[2] Brunner N. et. al., Rev. Mod. Phys. 86, (419)