Symmetric NPA relaxations¶
We initialize the RepLAB library.
[1]:
addpath([pwd, '/../../../external/replab']);
replab_init('verbose', 0);
replab_init: Initialization done.
We examine an upper bound on the quantum maximum of the CHSH inequality. We use the NPA hierarchy, and assume the moment matrix has been constructed in the symmetric subspace, as in arXiv:1808.09598.
We first declare the one and only variable, before constructing the moment matrix.
[2]:
y = sdpvar; % y is actually y_A0B0, but we save space by identifying y with it
C = [ 1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1];
A = [ 0 0 0 0 0
0 0 0 1 1
0 0 0 1 -1
0 1 1 0 0
0 1 -1 0 0];
X = C + A*y;
I_CHSH = 4*y; % it is <A0B0> + <A0B1> + <A1B0> - <A1B1>, and y_A1B1 = -y_A0B0.
optimize(X >= 0, -I_CHSH, sdpsettings('verbose', 0)); % sign change to maximize, and we don't show solver output
double(I_CHSH)
ans = 2.8284
To help with symmetrization, we’ll use an equivalent form of the constraint X, and check the result.
Now, we want to find a change of basis that brings X into a block-diagonal form. This is easier to check on the basis matrices C and A, i.e. we want to find U such that U*C*U' and U*A*U' are block-diagonal. Let us use RepLAB for that. First, remark that the moment matrix is invariant under the following signed permutations:
\(\vec{v} = \left( 1, A_0, A_1, B_0, B_1 \right) \rightarrow \left( 1, B_0, B_1, A_0, A_1 \right)\)
\(\vec{v} = \left( 1, A_0, A_1, B_0, B_1 \right) \rightarrow \left( 1, -A_0, -A_1, -B_0, -B_1 \right)\)
\(\vec{v} = \left( 1, A_0, A_1, B_0, B_1 \right) \rightarrow \left( 1, A_1, A_0, B_0, -B_1 \right)\)
We now show how to use RepLAB to find the change of basis for the symmetry group that includes only the first symmetry.
[3]:
g1 = [1 4 5 2 3]; % permutation of parties
g2 = [1 -2 -3 -4 -5]; % sign flip everywhere, note the signed permutation convention
g3 = [1 3 2 4 -5]; % additional symmetry
We build the symmetry group from those generators. For the permutation of parties only, the order (=size) of the group is 2. When adding the other symmetries, the order of the group should be 16.
[4]:
nElements = 5;G = replab.SignedPermutationGroup.of(g1, g2, g3);
G.order
ans =
16
If necessary, we can also expand the group elements (use G.elements.at(2) to get the second element for example).
[5]:
G.elements
ans =
{
[1,1] =
1 2 3 4 5
[1,2] =
1 -3 -2 -4 5
[1,3] =
1 2 -3 5 4
[1,4] =
1 3 -2 -5 4
[1,5] =
1 4 5 2 3
[1,6] =
1 -5 -4 -2 3
[1,7] =
1 4 -5 3 2
[1,8] =
1 5 -4 -3 2
[1,9] =
1 -5 4 3 -2
[1,10] =
1 -4 5 -3 -2
[1,11] =
1 5 4 2 -3
[1,12] =
1 -4 -5 -2 -3
[1,13] =
1 -3 2 5 -4
[1,14] =
1 -2 3 -5 -4
[1,15] =
1 3 2 4 -5
[1,16] =
1 -2 -3 -4 -5
}
Now, the representation we need is composed of signed permutation matrices; as we are considering a group of signed permutations, this is the group natural representation. Representations in RepLAB are described using the images of the generators. By calling G.sample, we get random elements from the group, and by calling rep.image(g) we get the image of a group element.
[6]:
rep = G.naturalRep
rep =
Orthogonal representation
dimension: 5
field: 'R'
group: Finite group of order 16
isUnitary: true
preimages{1}: [1, 4, 5, 2, 3]
images{1}: 5 x 5 double
preimages{2}: [1, -2, -3, -4, -5]
images{2}: 5 x 5 double
preimages{3}: [1, 3, 2, 4, -5]
images{3}: 5 x 5 double
[7]:
g = G.sample
g =
1 -5 4 3 -2
Now, we decompose the representation into irreducible components. For that, we call rep.decomposition. In the answer I(m)xR(d), we express that the component has \(m\) copies (multiplicity) of a Real representation of dimension d. We then play with the indices of the component, and copy.
Hint: explore rep.decomposition.component(i).irrep(j) for \((i,j) = (1,1), (2,1), (3,1)\).
[8]:
rep.decomposition
ans =
Orthogonal reducible representation
dimension: 5
field: 'R'
group: Finite group of order 16
isUnitary: true
parent: Orthogonal representation
basis(1,'double'): [1; 0; 0; 0; 0]
basis(2,'double'): [0; -0.5; 0.5; 1.1132e-17; -0.70711]
basis(3,'double'): [0; -0.5; -0.5; -0.70711; 0]
basis(4,'double'): [0; 0.60714; -0.36246; -0.17301; -0.68561]
basis(5,'double'): [0; -0.36246; -0.60714; 0.68561; -0.17301]
component(1): Isotypic component R(1) (trivial)
component(2): Isotypic component R(2) (nontrivial)
component(3): Isotypic component R(2) (nontrivial)
[9]:
rep.decomposition.component(1).irrep(1)
ans =
Orthogonal trivial irreducible real-type subrepresentation
dimension: 1
field: 'R'
group: Finite group of order 16
isUnitary: true
parent: Orthogonal representation
basis(1,'exact'): [1; 0; 0; 0; 0]
We now ask for the change of basis matrix, and verify that A and C block-diagonalize.
[10]:
U = rep.decomposition.basis
U'*C*U
U'*A*U
U =
1.0000 0 0 0 0
0 -0.5000 -0.5000 0.6071 -0.3625
0 0.5000 -0.5000 -0.3625 -0.6071
0 0.0000 -0.7071 -0.1730 0.6856
0 -0.7071 0 -0.6856 -0.1730
ans =
1.0000 0 0 0 0
0 1.0000 -0.0000 -0.0000 0.0000
0 -0.0000 1.0000 0.0000 0.0000
0 -0.0000 0.0000 1.0000 -0.0000
0 0.0000 0.0000 -0.0000 1.0000
ans =
0 0 0 0 0
0 1.4142 -0.0000 0 0.0000
0 -0.0000 1.4142 0.0000 -0.0000
0 0.0000 0.0000 -1.4142 0.0000
0 0.0000 -0.0000 0.0000 -1.4142
Once we have found the change of basis matrix, we transform our moment matrix X, and naturally it should have a block diagonal structure following the one of C and A – for that, it may be necessary to kill small off-block coefficients.
[11]:
X = (U'*C*U) + (U'*A*U)*y; % equivalent to the constraint above
optimize(X >= 0, -I_CHSH); % sign change to maximize
double(I_CHSH)
num. of constraints = 1
dim. of sdp var = 5, num. of sdp blk = 1 1th semidefinite block is actually diagonal
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
HKM 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
-------------------------------------------------------------------
0|0.000|0.000|8.0e-01|6.2e+00|5.0e+02| 5.000000e+01 0.000000e+00| 0:0:00| chol 1 1
1|1.000|1.000|6.5e-08|6.9e-02|4.3e+01| 3.984640e+01 4.400004e-01| 0:0:00| chol 1 1
2|0.933|1.000|1.5e-07|6.9e-03|2.9e+00| 3.700414e+00 8.665599e-01| 0:0:00| chol 1 1
3|1.000|0.949|3.5e-08|1.0e-03|2.8e-01| 3.024905e+00 2.752468e+00| 0:0:00| chol 1 1
4|0.987|0.988|1.1e-08|8.0e-05|3.4e-03| 2.830875e+00 2.827762e+00| 0:0:00| chol 1 1
5|0.989|0.989|3.5e-10|7.7e-06|3.8e-05| 2.828454e+00 2.828448e+00| 0:0:00| chol 1 1
6|0.989|0.989|8.9e-12|8.5e-08|4.2e-07| 2.828427e+00 2.828427e+00| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.00e-07
-------------------------------------------------------------------
number of iterations = 6
primal objective value = 2.82842742e+00
dual objective value = 2.82842735e+00
gap := trace(XZ) = 4.17e-07
relative gap = 6.26e-08
actual relative gap = 1.04e-08
rel. primal infeas (scaled problem) = 8.94e-12
rel. dual " " " = 8.49e-08
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual " " " = 0.00e+00
norm(X), norm(y), norm(Z) = 2.0e+00, 7.1e-01, 3.0e+00
norm(A), norm(b), norm(C) = 3.8e+00, 5.0e+00, 3.2e+00
Total CPU time (secs) = 0.11
CPU time per iteration = 0.02
termination code = 0
DIMACS: 8.9e-12 0.0e+00 1.4e-07 0.0e+00 1.0e-08 6.3e-08
-------------------------------------------------------------------
ans = 2.8284
As an exercice, recover an algebraic basis from U, guess the form of the semidefinite program when fully block-diagonal (i.e. for the full group of CHSH symmetries). Remark that now the problem can be solved by hand to recover the \(2\sqrt{2}\) bound, as the SDP matrix is fully diagonal, and thus corresponds to linear inequalities.