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
    sdp.add_constraint(rho >> 0)
        #Trace 1
    sdp.add_constraint(pcs.trace(rho) == 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)