# Optimal randomness generation from entangled quantum states

Posted on 05 August 2015

# Introduction¶

This notebook is the computational appendix to the paper Optimal randomness generation from entangled quantum states (arXiv:1505.03837). We detail how to obtain the numerical results.

To run this notebook, the dependencies are NumPy, MatPlotLib, QuTip, SymPy, and Ncpol2sdpa>=1.7. Furthermore, to solve the semidefinite programmes (SDPs), a member of the SDPA family of solvers must be in the path.

As a first step, we import all necessary functions:

In :
from math import atan
import matplotlib.pyplot as plt
from numpy import *
from qutip import *
from ncpol2sdpa import *


## Maximum guessing probability with the elegant inequality and a POVM¶

First we set the Bell scenario: Alice has three dichotomic measurements and one four-output POVM. Bob has four dichotomic measurements. We replicate the elegant Bell inequality  in the Collins-Gisin notation , leaving blank spaces for Alice's POVM:

In :
A_configuration = [2, 2, 2, 4]
B_configuration = [2, 2, 2, 2]
I = [[ 0, -1.5, 0.5, 0.5,  0.5],
[ 0,    1,   1,  -1,   -1],
[ 0,    1,  -1,   1,   -1],
[ 0,    1,  -1,  -1,    1],
[ 0,    0,   0,   0,    0],
[ 0,    0,   0,   0,    0],
[ 0,    0,   0,   0,    0]]


Next we define the algebra of the measurement operators:

In :
A = generate_measurements(A_configuration, 'A')
B = generate_measurements(B_configuration, 'B')
substitutions = projective_measurement_constraints(A, B)


We are looking for the maximum guessing probability, so the constraint will be the Bell inequality which also includes the POVM tetrahedron. The tetrahedron coefficients have a $k$ parameter, but the numerical calculations are insensitive to it. Since we want the Bell inequality to match the maximum violation, we add a pair of inequality constraints.

In :
k = 1.0
bell_inequality = - define_objective_with_I(I, A, B) - \
k*(A*B + A*B + A*B +
(1-A-A-A)*B)
moments = [-sqrt(3) + bell_inequality]


The objective is the maximum guessing probability. The SDPA solver can only minimize, so we flip the sign. It is sufficient to study one output, as the results are symmetric in the POVM.

In :
objective = - A


Here we give an example of a level 2+AAB+ABB relaxation. Higher level relaxations are straightforward generalization of this. Obtaining the relaxation can take a few minutes.

In :
level = 2
AAB = [Ak*Ai*Bj for Ak in flatten(A) for Ai in flatten(A) for Bj in flatten(B)]
ABB = [Ak*Bi*Bj for Ak in flatten(A) for Bi in flatten(B) for Bj in flatten(B)]
sdpRelaxation = SdpRelaxation(flatten([A, B]))
sdpRelaxation.get_relaxation(level, objective=objective,
momentequalities=moments, substitutions=substitutions,
extramonomials=flatten([AAB, ABB]))


Solving the SDP takes another few minutes.

In :
sdpRelaxation.solve()
print(sdpRelaxation.primal, sdpRelaxation.dual)

(-0.2514585418483301, -0.2508329684907338)


The regular-precision solver will not give an accurate result, it will yield a large gap between the primal and dual. To solve the problem with higher precision, we export it to a file.

In :
sdpRelaxation.write_to_file("elegant_povm.dat-s")


Oddly, solving this with the arbitrary-precision solver SDPA_GMP or the quad-double-precision solver SDPA_QD will lead to divergent results. This is due to the constraint that involves $\sqrt{3}$. It is written to the file at regular double precision, which leads to the divergent behaviour with the higher precision solvers. To rectify the problem, the file can be manually edited to include $\sqrt{3}$ at a higher precision. We used six extra digits, that is, $1.7320508075688772$. With this setting, we get $0.25$ up to $10^{-8}$ precision at this level of relaxation using the quad-double-precision solver, with a relative gap between the primal and dual on the order of $10^{-57}$.

## Noise sensitivity of the analytical result¶

We analyze the noise sensitivity using the full correlations , therefore we need a function that generates a behaviour for a given state and a measurement. For convenience, we define a function to calculate the dot product $\vec{\sigma}\vec{v}$.

In :
def dp(v, sigma):
result = v*sigma
for i in range(1, 3):
result += v[i]*sigma[i]
return result


Next we define the joint probability as the function of $p$ the parameter in a Werner state. The measurements correspond to the ones described in the manuscript, with the POVM being the tetrahedron.

In :
def joint_probabilities(p):
psi = (tensor(basis(2,0),basis(2,0)) + tensor(basis(2,1),basis(2,1))).unit()
rho = Qobj(p*(psi*psi.dag()).data+(1-p)*qeye(4).data/4, dims=[[2,2], [2,2]])
A_00 = (qeye(2)+sigmaz())/2
A_10 = (qeye(2)+sigmax())/2
A_20 = (qeye(2)+sigmay())/2

B_00 = (qeye(2)+(sigmaz()+sigmax())/sqrt(2))/2
B_10 = (qeye(2)+(sigmaz()-sigmax())/sqrt(2))/2
B_20 = (qeye(2)+(sigmaz()-sigmay())/sqrt(2))/2
B_30 = (qeye(2)+(sigmaz()+sigmay())/sqrt(2))/2
B_40 = (qeye(2)+(sigmax()-sigmay())/sqrt(2))/2
B_50 = (qeye(2)+(sigmax()+sigmay())/sqrt(2))/2

b = [array([ 1, -1, -1])/sqrt(3),
array([-1, -1,  1])/sqrt(3),
array([ 1,  1,  1])/sqrt(3),
array([-1,  1, -1])/sqrt(3)]
B_6= [(qeye(2)+dp(bj, [sigmax(), sigmay(), sigmaz()]))/4 for bj in b]

probabilities = [expect(tensor(A_00, qeye(2)), rho),
expect(tensor(A_10, qeye(2)), rho),
expect(tensor(A_20, qeye(2)), rho),

expect(tensor(qeye(2), B_00), rho),
expect(tensor(qeye(2), B_10), rho),
expect(tensor(qeye(2), B_20), rho),
expect(tensor(qeye(2), B_30), rho),
expect(tensor(qeye(2), B_40), rho),
expect(tensor(qeye(2), B_50), rho),
expect(tensor(qeye(2), B_6), rho),
expect(tensor(qeye(2), B_6), rho),
expect(tensor(qeye(2), B_6), rho),

expect(tensor(A_00, B_00), rho),
expect(tensor(A_00, B_10), rho),
expect(tensor(A_00, B_20), rho),
expect(tensor(A_00, B_30), rho),
expect(tensor(A_00, B_40), rho),
expect(tensor(A_00, B_50), rho),
expect(tensor(A_00, B_6), rho),
expect(tensor(A_00, B_6), rho),
expect(tensor(A_00, B_6), rho),

expect(tensor(A_10, B_00), rho),
expect(tensor(A_10, B_10), rho),
expect(tensor(A_10, B_20), rho),
expect(tensor(A_10, B_30), rho),
expect(tensor(A_10, B_40), rho),
expect(tensor(A_10, B_50), rho),
expect(tensor(A_10, B_6), rho),
expect(tensor(A_10, B_6), rho),
expect(tensor(A_10, B_6), rho),

expect(tensor(A_20, B_00), rho),
expect(tensor(A_20, B_10), rho),
expect(tensor(A_20, B_20), rho),
expect(tensor(A_20, B_30), rho),
expect(tensor(A_20, B_40), rho),
expect(tensor(A_20, B_50), rho),
expect(tensor(A_20, B_6), rho),
expect(tensor(A_20, B_6), rho),
expect(tensor(A_20, B_6), rho)]

return probabilities


We formally define the scenario and the projector algebra again. The difference is that in the full-probability picture, the maximum guessing probability may not be extremal. Hence we construct four separate algebras, corresponding to the four outputs.

In :
A_configuration = [2, 2, 2]
B_configuration = [2, 2, 2, 2, 2, 2, 4]

P_0_A = generate_measurements(A_configuration, 'P_0_A')
P_0_B = generate_measurements(B_configuration, 'P_0_B')
substitutions1 = projective_measurement_constraints(P_0_A, P_0_B)
P_1_A = generate_measurements(A_configuration, 'P_1_A')
P_1_B = generate_measurements(B_configuration, 'P_1_B')
substitutions2 = projective_measurement_constraints(P_1_A, P_1_B)
P_2_A = generate_measurements(A_configuration, 'P_2_A')
P_2_B = generate_measurements(B_configuration, 'P_2_B')
substitutions3 = projective_measurement_constraints(P_2_A, P_2_B)
P_3_A = generate_measurements(A_configuration, 'P_3_A')
P_3_B = generate_measurements(B_configuration, 'P_3_B')
substitutions4 = projective_measurement_constraints(P_3_A, P_3_B)
substitutions = dict(substitutions1.items() + substitutions2.items() +
substitutions3.items() + substitutions4.items())


The guessing probability in this case is the sum of the POVM outcomes. Since we work in the Collins-Gisin notation, the last outcome of a measurement is not included in the algebra. This can easily be constructed by deducting the rest of the operators in a measurement from the identity. This, however, is more complicated here, because each algebra generates a moment matrix, and the top left is not the expectation value with the identity operator. So instead, we will have to include a special variable in the objective function, the top left element of the last moment matrix.

In :
guessing_probability = -(P_0_B + P_1_B + P_2_B +
(-P_3_B-P_3_B-P_3_B))
extraobjexpr = "-3[0,0]"


We do a level 2 relaxation.

In :
level = 2
sdpRelaxation = SdpRelaxation([flatten([P_0_A, P_0_B]),
flatten([P_1_A, P_1_B]),
flatten([P_2_A, P_2_B]),
flatten([P_3_A, P_3_B])],
normalized=False)


The following function obtains the probabilities and defines the constraint $\sum_a p_a=p$. It also normalizes the top-left elements of the four moment matrices. On its first call it obtains the relaxation, and in subsequent calls, it simply replaces the constraints and solves the SDP.

In :
def get_maximum_guessing_probability(p, first=False):
probabilities = joint_probabilities(p)

moments = []
k=0
#Marginals first
for i in range(3):
moments.append(P_0_A[i]+P_1_A[i]+P_2_A[i]+P_3_A[i]-probabilities[k])
k += 1
for j in range(6):
moments.append(P_0_B[j]+P_1_B[j]+P_2_B[j]+P_3_B[j]-probabilities[k])
k += 1
for j in range(3):
moments.append(P_0_B[j]+P_1_B[j]+P_2_B[j]+P_3_B[j]-probabilities[k])
k += 1

#Joint probabilities
for i in range(3):
for j in range(6):
moments.append(P_0_A[i]*P_0_B[j]+P_1_A[i]*P_1_B[j]+P_2_A[i]*P_2_B[j]+P_3_A[i]*P_3_B[j]-probabilities[k])
k += 1
for j in range(3):
moments.append(P_0_A[i]*P_0_B[j]+P_1_A[i]*P_1_B[j]+P_2_A[i]*P_2_B[j]+P_3_A[i]*P_3_B[j]-probabilities[k])
k += 1
# We need to normalize the top-left elements of the moment matrices
moments.append(["0[0,0]+1[0,0]+2[0,0]+3[0,0]-1")
if first:
sdpRelaxation.get_relaxation(level, objective=guessing_probability,
momentequalities=moments,
substitutions=substitutions,
extraobjexpr=extraobjexpr)
else:
sdpRelaxation.process_constraints(momentequalities=moments,)

sdpRelaxation.solve()
return -sdpRelaxation.primal


We call this function over a range of values to obtain the results. This can take over half an hour to calculate.

In :
noise = linspace(0.85, 1.0, 30)
sensitivity = [get_maximum_guessing_probability(noise, first = True)]
for p in noise[1:]:
sensitivity.append(get_maximum_guessing_probability(p))
print(sensitivity)

[1.0000004872175416, 1.0000005240736838, 1.0000005470789861, 1.000000499589071, 1.0000004379654095, 1.0000005968134946, 1.0000005408672348, 0.990881754645857, 0.9760158851876728, 0.9608173204563474, 0.9452661731202696, 0.9293556282111481, 0.913061447027384, 0.8963600608723868, 0.8792227510613924, 0.8616165289443828, 0.8435031830118433, 0.824835204595386, 0.8055536364857269, 0.7855802884285172, 0.7648133581539607, 0.7431125290346845, 0.7202829011964192, 0.696045828145155, 0.6699809384876207, 0.6414145008962466, 0.6091843638554746, 0.5709367331409689, 0.520258797053067, 0.3016279816545112]

In :
%matplotlib inline
fig, ax = plt.subplots()
plt.plot(noise, [-log2(s) for s in sensitivity], linewidth=2)

Out:
[<matplotlib.lines.Line2D at 0x7fcce711f050>] ## Noise sensitivity of the [(2, 2, 2, 4) (2, 2, 2, 2)] scenario¶

The case is very similar to the previous one. First we define the joint probability:

In :
def joint_probabilities(p):
psi = (tensor(basis(2,0),basis(2,0)) + tensor(basis(2,1),basis(2,1))).unit()
rho = Qobj(p*(psi*psi.dag()).data+(1-p)*qeye(4).data/4, dims=[[2,2], [2,2]])
A_00 = (qeye(2)+sigmaz())/2
A_10 = (qeye(2)+sigmax())/2
A_20 = (qeye(2)+sigmay())/2
sigma = [sigmax(), sigmay(), sigmaz()]
b = [array([ 1, -1, -1])/sqrt(3),
array([-1, -1,  1])/sqrt(3),
array([ 1,  1,  1])/sqrt(3),
array([-1,  1, -1])/sqrt(3)]
A_3 = [(qeye(2)+dp(bj, sigma))/4 for bj in b]

B_00 = (qeye(2)+dp([ 1, -1,  1], sigma)/sqrt(3))/2
B_10 = (qeye(2)+dp([-1,  1,  1], sigma)/sqrt(3))/2
B_20 = (qeye(2)+dp([ 1,  1, -1], sigma)/sqrt(3))/2
B_30 = (qeye(2)+dp([-1, -1, -1], sigma)/sqrt(3))/2

probabilities = [expect(tensor(A_00, qeye(2)), rho),
expect(tensor(A_10, qeye(2)), rho),
expect(tensor(A_20, qeye(2)), rho),
expect(tensor(A_3, qeye(2)), rho),
expect(tensor(A_3, qeye(2)), rho),
expect(tensor(A_3, qeye(2)), rho),

expect(tensor(qeye(2), B_00), rho),
expect(tensor(qeye(2), B_10), rho),
expect(tensor(qeye(2), B_20), rho),
expect(tensor(qeye(2), B_30), rho),

expect(tensor(A_00, B_00), rho),
expect(tensor(A_00, B_10), rho),
expect(tensor(A_00, B_20), rho),
expect(tensor(A_00, B_30), rho),

expect(tensor(A_10, B_00), rho),
expect(tensor(A_10, B_10), rho),
expect(tensor(A_10, B_20), rho),
expect(tensor(A_10, B_30), rho),

expect(tensor(A_20, B_00), rho),
expect(tensor(A_20, B_10), rho),
expect(tensor(A_20, B_20), rho),
expect(tensor(A_20, B_30), rho),

expect(tensor(A_3, B_00), rho),
expect(tensor(A_3, B_10), rho),
expect(tensor(A_3, B_20), rho),
expect(tensor(A_3, B_30), rho),

expect(tensor(A_3, B_00), rho),
expect(tensor(A_3, B_10), rho),
expect(tensor(A_3, B_20), rho),
expect(tensor(A_3, B_30), rho),

expect(tensor(A_3, B_00), rho),
expect(tensor(A_3, B_10), rho),
expect(tensor(A_3, B_20), rho),
expect(tensor(A_3, B_30), rho)]

return probabilities


Next we define the algebra, the objective function, and we formally initialize the level-2 relaxation.

In :
A_configuration = [2, 2, 2, 4]
B_configuration = [2, 2, 2, 2]
P_0_A = generate_measurements(A_configuration, 'P_0_A')
P_0_B = generate_measurements(B_configuration, 'P_0_B')
substitutions1 = projective_measurement_constraints(P_0_A, P_0_B)
P_1_A = generate_measurements(A_configuration, 'P_1_A')
P_1_B = generate_measurements(B_configuration, 'P_1_B')
substitutions2 = projective_measurement_constraints(P_1_A, P_1_B)
P_2_A = generate_measurements(A_configuration, 'P_2_A')
P_2_B = generate_measurements(B_configuration, 'P_2_B')
substitutions3 = projective_measurement_constraints(P_2_A, P_2_B)
P_3_A = generate_measurements(A_configuration, 'P_3_A')
P_3_B = generate_measurements(B_configuration, 'P_3_B')
substitutions4 = projective_measurement_constraints(P_3_A, P_3_B)
substitutions = dict(substitutions1.items() + substitutions2.items() +
substitutions3.items() + substitutions4.items())
guessing_probability = -(P_0_A + P_1_A + P_2_A + (-P_3_A-P_3_A-P_3_A))
extraobjexpr = "-3[0,0]"
level = 2
sdpRelaxation = SdpRelaxation([flatten([P_0_A, P_0_B]),
flatten([P_1_A, P_1_B]),
flatten([P_2_A, P_2_B]),
flatten([P_3_A, P_3_B])],
normalized=False)


We need to define the function to get the maximum guessing probability again.

In :
def get_maximum_guessing_probability(p, first=False):
probabilities = joint_probabilities(p)

moments = []
k = 0
#Marginals first
for i in range(3):
moments.append(P_0_A[i]+P_1_A[i]+P_2_A[i]+P_3_A[i]-probabilities[k])
k += 1
for j in range(3):
moments.append(P_0_A[j]+P_1_A[j]+P_2_A[j]+P_3_A[j]-probabilities[k])
k += 1
for j in range(4):
moments.append(P_0_B[j]+P_1_B[j]+P_2_B[j]+P_3_B[j]-probabilities[k])
k += 1

#Joint probabilities
for i in range(3):
for j in range(4):
moments.append(P_0_A[i]*P_0_B[j]+P_1_A[i]*P_1_B[j]+P_2_A[i]*P_2_B[j]+P_3_A[i]*P_3_B[j]-probabilities[k])
k += 1
for i in range(3):
for j in range(4):
moments.append(P_0_A[i]*P_0_B[j]+P_1_A[i]*P_1_B[j]+P_2_A[i]*P_2_B[j]+P_3_A[i]*P_3_B[j]-probabilities[k])
k += 1
# We need to normalize the top-left elements of the moment matrices
moments.extend("0[0,0]+1[0,0]+2[0,0]+3[0,0]-1")
if first:
sdpRelaxation.get_relaxation(level, objective=guessing_probability,
momentequalities=moments,
substitutions=substitutions,
extraobjexpr=extraobjexpr)
else:
sdpRelaxation.process_constraints(momentequalities=moments)

sdpRelaxation.solve()
return -sdpRelaxation.primal


The rest of the solution is identical.

In :
noise = linspace(0.85, 1.0, 30)
sensitivity = [get_maximum_guessing_probability(noise, first = True)]
for p in noise[1:]:
sensitivity.append(get_maximum_guessing_probability(p))
print(sensitivity)

[1.0000006077543446, 1.0000006496529297, 1.0000006452327943, 1.000000554308008, 1.0000005385866668, 1.000000576607919, 1.0000006427110755, 1.000000189347059, 1.0000000971849392, 1.0000001969583598, 0.9968305657655212, 0.9738482865458454, 0.9500254170240093, 0.9257106571332189, 0.900979776661172, 0.8758534765779566, 0.8503319167826222, 0.824396612304085, 0.7980154925968846, 0.7711340838143288, 0.7436695481474971, 0.7155067003731853, 0.6864836031710995, 0.6563755399276776, 0.6248556844800475, 0.5914286714525884, 0.5551494202256076, 0.5127876722677261, 0.4522162956410998, 0.2568499549706953]

In :
plt.plot(noise, [-log2(s) for s in sensitivity], linewidth=2)

Out:
[<matplotlib.lines.Line2D at 0x7fcce96bcf50>] ## Global randomness in the [(2, 2, 2, 2, 2, 2, 2, 4) (2, 2, 2, 2, 2, 2, 2, 4)] scenario¶

For the global randomness, we use a specific Bell inequality, defined by the following matrix in the Collins-Gisin picture:

In :
delta = 0.0676946837698
k = 1.0
I= [[      0,          -0.5,      -0.5,      -0.5, -(1+delta)/2, -(1-delta)/2,-(-1+delta)/2,-(-1-delta)/2-k,    0,    0,     0],
[   -0.5,             1,         0,         0,            1,            1,           -1,             -1,    0,    0,     0],
[   -0.5,             0,         1,         0,        delta,       -delta,            0,              0,    0,    0,     0],
[   -0.5,             0,         0,         1,            0,            0,        delta,         -delta,    0,    0,     0],
[  -(1+delta)/2,  delta,         0,         1,            0,            0,            0,              0,   -k,    0,     0],
[  -(1-delta)/2, -delta,         0,         1,            0,            0,            0,              0,    0,   -k,     0],
[  -(-1+delta)/2,     0,     delta,        -1,            0,            0,            0,              0,    0,    0,    -k],
[  -(-1-delta)/2-k,   0,    -delta,        -1,            0,            0,            0,              0,    k,    k,     k],
[      0,             0,         0,         0,           -k,            0,            0,              k,    0,    0,     0],
[      0,             0,         0,         0,            0,           -k,            0,              k,    0,    0,     0],
[      0,             0,         0,         0,            0,            0,           -k,              k,    0,    0,     0]]


We obtained the $\delta$ value by nonconvex optimization. Defining the operators is straightforward:

In :
A_configuration = [2, 2, 2, 2, 2, 2, 2, 4]
B_configuration = [2, 2, 2, 2, 2, 2, 2, 4]
A = generate_measurements(A_configuration, 'A')
B = generate_measurements(B_configuration, 'B')
AB = [Ai*Bj for Ai in flatten(A) for Bj in flatten(B)]
substitutions = projective_measurement_constraints(A, B)


We impose the maximum violation as a constraint and optimize for the guessing probability. The guessing probability is symmetric for all joint output of the POVM, hence it is enough to look at the first pair.

In :
level = 1
equality = -define_objective_with_I(I, A, B) - ((3.0+8.0*sqrt(1+delta**2)-11)/4+2)
sdpRelaxation = SdpRelaxation(flatten([A, B]), verbose=2)
sdpRelaxation.get_relaxation(level, substitutions=substitutions, objective=-A*B,
equalities=[equality], extramonomials=AB)
sdpRelaxation.solve()
print(sdpRelaxation.primal, sdpRelaxation.dual)

Calculating block structure...
Estimated number of SDP variables: 7380
Generating moment matrix...
Reduced number of SDP variables: 5412
Processing 2/2 constraints...
(-0.20333663051476122, -0.19042314701932123)


The gap between the primal and dual solution is large. Solving the above SDP with SeDuMi provides a smaller gap with the value reported in the manuscript. Exporting the SDP is done with the following command:

In :
sdpRelaxation.write_to_file("global_randomness.dat-s")


 Gisin, N. Bell inequalities: Many questions, a few answers. Quantum Reality, Relativistic Causality, and Closing the Epistemic Circle, Springer Netherlands, 2009, 73, 125--138.

 Collins, D. & Gisin, N. A relevant two qubit Bell inequality inequivalent to the CHSH inequality. Journal of Physics A: Mathematical and General, IOP Publishing, 2004, 37, 1775.

 Nieto-Silleras, O.; Pironio, S. & Silman, J. Using complete measurement statistics for optimal device-independent randomness evaluation. New Journal of Physics, 2014, 16, 013035.

Share on:    