Examples file for LPAssumptions Mathematica package (https://github.com/rogercolbeck/LPAssumptions)

Example 1

Maximise y+ax (a≥1) subject to x≤1, x+y/2≤2, x,y≥0. For v={x,y}, writing this in the form Maximise c.v subject to M.v≤ b we have the following where bnew encodes the fact that both constraints are of “less-than-equal-to” form. Note that v≥ 0 is always implicitly assumed.

c={a,1};

M={{1,0},{1,1/2}};

b={1,2};
bnew=Thread[{b,-1}];

assump={a≥ 1};

Flip the sign of the objective function in order to turn maximization to minimization for LPAssumptions (see documentation for format). LPAssumptions lists the assumptions on the variable parameter (here, a) and the corresponding optimal point as a vector by default. If the OutputValues option is set to True, then it instead outputs the minimum value which is simply the dot product of the optimising vector and the objective. For our problem, this tells us that for a<2, the maximum value is attained at (0,4) and is equal to 4, while for a>2 it is attained at at (1,2) and equals 2+a.

out1=LPAssumptions[-c,M,bnew]

-out1[[1]][[-1]].c

-out1[[2]][[-1]].c

LPAssumptions[-c,M,bnew,OutputValues-> True]

Compare with Mathematica’s LinearProgramming for fixed values of a.

vec1=LinearProgramming[-c/.a->1/2,-M,-bnew]

vec2=LinearProgramming[-c/.a->3,-M,-bnew]

-vec1.c/.a->1/2

-out1[[1]][[-1]].c/.a->1/2

-vec2.c/.a->3

-out1[[2]][[-1]].c/.a->3

Example 2

Here we find the local weight of a distribution in the Bell scenario as a function of the detector efficiency eta and a noise parameter epsilon. The local weight of a distribution Π in a Bell scenario is the solution to the LP problem: Maximise (∑_i) (Subscript[x, i]), subject to (∑_i) (Subscript[x, i])(P^(L,i))≤ Π, (Subscript[x, i ])≥0  for all i, where (Subscript[{(P^(L,i))}, i]) is the set of all local deterministic distributions in that Bell scenario (see Refs. [1,2]). Here, we consider the bipartite Bell scenario with 2 inputs and 3 outputs (adding an additional “non detection” output to characterise detector inefficiency) per each party. The matrix mat contains as columns the 81 local deterministic distributions of this scenario. dist is the distribution whose local weight we wish to find as a function of eta and epsilon.  The bound we obtain in this example is the one found in Section IV D F of [3].

mat={{1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0},{1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0},{0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0},{0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0},{0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1},{0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1},{1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0},{1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0},{0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0},{0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0},{0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1},{0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1}};

Here is the distribution observed by noisy PR boxes with noise value eps (between 0 and 1/4) and detector efficiency eta (between 0 and 1)

dist[eps_,eta_]:={{(1/2) (1-2 eps) (eta^2),eps (eta^2),(1/2) (1-eta) eta,(1/2) (1-2 eps) (eta^2),eps (eta^2),(1/2) (1-eta) eta},{eps (eta^2),(1/2) (1-2 eps) (eta^2),(1/2) (1-eta) eta,eps (eta^2),(1/2) (1-2 eps) (eta^2),(1/2) (1-eta) eta},{(1/2) (1-eta) eta,(1/2) (1-eta) eta,(((1-eta))^2),(1/2) (1-eta) eta,(1/2) (1-eta) eta,(((1-eta))^2)},{(1/2) (1-2 eps) (eta^2),eps (eta^2),(1/2) (1-eta) eta,eps (eta^2),(1/2) (1-2 eps) (eta^2),(1/2) (1-eta) eta},{eps (eta^2),(1/2) (1-2 eps) (eta^2),(1/2) (1-eta) eta,(1/2) (1-2 eps) (eta^2),eps (eta^2),(1/2) (1-eta) eta},{(1/2) (1-eta) eta,(1/2) (1-eta) eta,(((1-eta))^2),(1/2) (1-eta) eta,(1/2) (1-eta) eta,(((1-eta))^2)}}

Note that our code only works when the variables are in the objective function, thus we take the dual of the above problem so that the variable distribution, dist now appears as the objective. The objective, obj of the primal problem (just a vector of ones) now becomes the constraint vector. Because all the constraints are of the form greater than or equal to we can use obj as the constraint vector without threading it with 1s.

obj=Table[1,{x,1,(3^4)}];

We restrict the region as below because outside this region we already know that the distribution is local (has local weight 1).

assump={eta>2/3,eta<1,eps>0,eps<1/8};

In this case we only care about the local weight (not the point where it occurs), so set 
OutputValues->True

out=LPAssumptions[Flatten[dist[eps,eta]],Transpose[mat],obj,assump,{OutputValues->True,PrintTemp->True}];

These are the values obtained:

For[i=1,i≤Dimensions[out][[1]],i++,Print[out[[i]][[-1]]]]

LPAssumptions in its current form does not provide the most compressed output, there are multiple assumptions under which the same optimum is obtained and these can be simplified further. For example, one can use Mathematica’s Reduce function to simply the assumptions

The assumptions under which the value is 1+2 eta+(-3+8 eps) (eta^2) are

out[[1]][[1]]

Reduce can be used to compress these to

Reduce[out[[1]][[1]]]

In all other cases the value is 1 and hence the distribution is local.  Therefore the local weight of the distribution is 1+2 eta+(-3+8 eps) (eta^2) for (2/3)<eta<1&&0<eps<((-2+3 eta)/(8 eta)) and 1 otherwise.

Compare with  Mathematica’s LinearProgramming for fixed values of eta and epsilon.

obj.LinearProgramming[-obj,-mat,-Flatten[dist[1/20,9/10]]]

Simplify[out[[1]][[-1]]]/.{eps->1/20,eta->9/10}

Example 3

This is the same as Example 1 but we replace a with (a^d) and take a≥1 , d≥0 as the assumptions.

c={(a^d),1};

M={{1,0},{1,1/2}};

b={1,2};

bnew=Thread[{b,-1}];

assump={a≥ 1,d≥ 1};

In this case, (as of v12.0) Mathematica is unable to resolve the following to True (although it is always true).

Simplify[(a^d)≥  1,assump]

For large enough a, Mathematica can resolve this:

Simplify[(2^d)≥ 1,assump]

If commented back in, the following line of code does not finish running because of the above problem. This can be spotted in the temporary messages which can be printed during run time by setting the PrintTemp option to True. With these we see that the code gets stuck making the same assumption again and again.

(* out1=LPAssumptions[-c,M,bnew,assump,PrintTemp->True] *)

By replacing (a^d) by t and assuming t ≥ 1 we recover Example 1, which is solved easily.

Furthermore, for fixed values of a or d, Simplify is able to resolve the assumptions and our code works as usual

a=2;

c={(a^d),1};

M={{1,0},{1,1/2}};

b={1,2};

bnew=Thread[{b,-1}];

assump={d≥1};

out1=LPAssumptions[-c,M,bnew,assump]

Compare with Mathematica’s LinearProgramming for fixed values of d.

-c.LinearProgramming[-c/.d-> 1,-M,-bnew]/.d-> 1

-out1[[1]][[2]].c/.d-> 1

-c.LinearProgramming[-c/.d-> 2,-M,-bnew]/.d-> 2

-out1[[2]][[2]].c/.d->2

Clear[a]

Setting d=1, we again retrieve Example 1.

d=1;

c={(a^d),1};

M={{1,0},{1,1/2}};

b={1,2};

bnew=Thread[{b,-1}];

assump={a≥ 1};

out1=LPAssumptions[-c,M,bnew,assump]

References

[1] M. Zukowski, D. Kaszlikowski, A. Baturo, and J.-Åke Larsson,  arXiv:quant-ph/9910058 (1999).
[2] T. Cope, and R. Colbeck, Physical Review A 100, 22114 (2019).
[3] S. Pirandola, U. L. Andersen, L. Banchi, M. Berta, D. Bunandar, R. Colbeck, D. Englund, T. Gehring, C. Lupo, C. Ottaviani, J. Pereira, M. Razavi, J. S. Shaari, M. Tomamichel, V. C. Usenko, G. Vallone, P. Villoresi and P. Wallden, arXiv:1906.01645 (2019).