6 minutes

# Python - LP/SDP for Quantum Information with PICOS

## SDP interface and Solver

To write LP/SDP in Python I recommend the PICOS interface. 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 firendly developper community.

Regarding SDP solvers, in my experience, nothing compares to MOSEK. Sadly, this solver is proprietary (if price is of your concern, MOSEK offer a free one-year license for students). CVXOPT is a relatively good free-licensed alternative to MOSEK. However, in my experience, CXOPT have memory-leakage issue, bugs when solving 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 derministic strategies `p_d`

, i.e. `p`

is composed of 0s and 1s.
The local polytope can be written as a the set of its vertices `L=p_d`

.

Here a simple function that construct the local polytope

```
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.
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 value 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 value a two-qubit state can achieve. Such 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 expectedviolation of `2.8284...`

with 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)

1160 Words

13130-12-12 00:130