BellPolytopes.jl

This package addresses the membership problem for local polytopes: it constructs Bell inequalities and local models in multipartite Bell scenarios with arbitrary settings.

The original article for which it was written can be found here:

Improved local models and new Bell inequalities via Frank-Wolfe algorithms.

Installation

The most recent release is available via the julia package manager, e.g., with

using Pkg
Pkg.add("BellPolytopes")

or the main branch:

Pkg.add(url="https://github.com/ZIB-IOL/BellPolytopes.jl", rev="main")

Getting started

Let's say we want to characterise the nonlocality threshold obtained with the two-qubit maximally entangled state and measurements whose Bloch vectors form an icosahedron. Using BellPolytopes.jl, here is what the code looks like.

julia> using BellPolytopes, Ket, LinearAlgebra

julia> rho = state_phiplus(Float64) # two-qubit maximally entangled state
4×4 Hermitian{Float64, Matrix{Float64}}:
 0.5  0.0  0.0  0.5
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.5  0.0  0.0  0.5

julia> φ = (1 + √5) / 2;

julia> v = [0 1 φ; 0 1 -φ; 1 φ 0; 1 -φ 0; φ 0 1; φ 0 -1] / sqrt(2 + φ) # Bloch vectors forming an icosahedron
6×3 Matrix{Float64}:
 0.0        0.525731   0.850651
 0.0        0.525731  -0.850651
 0.525731   0.850651   0.0
 0.525731  -0.850651   0.0
 0.850651   0.0        0.525731
 0.850651   0.0       -0.525731

julia> mes = povm_dichotomic.(bloch_operator.(eachrow(v)));

julia> p = tensor_correlation(rho, mes, 2; marg=false)
6×6 Matrix{Float64}:
  0.447214  -1.0       -0.447214   0.447214   0.447214  -0.447214
 -1.0        0.447214  -0.447214   0.447214  -0.447214   0.447214
 -0.447214  -0.447214  -0.447214   1.0        0.447214   0.447214
  0.447214   0.447214   1.0       -0.447214   0.447214   0.447214
  0.447214  -0.447214   0.447214   0.447214   1.0        0.447214
 -0.447214   0.447214   0.447214   0.447214   0.447214   1.0

julia> lower_bound, upper_bound, local_model, bell_inequality = nonlocality_threshold(p);

julia> println([lower_bound, upper_bound])
[0.778, 0.779]

julia> final_iterate = sum(weight * atom for (weight, atom) in local_model);

julia> norm(final_iterate - lower_bound * p) < 1e-3 # checking local model
true

julia> local_bound_correlation(bell_inequality)[1] / dot(bell_inequality, p) # checking the Bell inequality
0.7785490499446976

Under the hood

The computation is based on an efficient variant of the Frank-Wolfe algorithm to iteratively find the local point closest to the input correlation tensor. See this review for an introduction to the method and the package FrankWolfe.jl for the implementation on which this package relies.

In a nutshell, each step gets closer to the objective point:

• either by moving towards a good vertex of the local polytope,
• or by astutely combining the vertices (or atoms) already found and stored in the active set.

julia> res = bell_frank_wolfe(p; v0=0.8, verbose=3, callback_interval=10^2, mode_last=-1);
   #Inputs: 6
 Symmetric: true
 Dimension: 21
Visibility: 0.8
   Iteration        Primal      Dual gap    Time (sec)       #It/sec    #Atoms       #LMO
         100    7.5008e-03    1.2306e-01    7.0849e-03    1.4114e+04        15         23
         200    1.8452e-03    4.2474e-02    9.2252e-03    2.1680e+04        14         27
         300    1.8093e-03    2.2514e-07    1.2580e-02    2.3847e+04        14         36
        Last    1.8093e-03    7.6190e-08    1.2739e-02    2.3786e+04        14         37
        Last    1.8093e-03    7.6190e-08    1.3419e-02    2.2580e+04        14         38
v_c ≤ 0.778392

Exact certificates after a run

dyadic_certificate rounds the active-set weights to a common binary denominator, reconstructs the local mixture with exact integer arithmetic, and returns a rational visibility lower bound. It accepts the solver result, its active set (res[5]), or a saved BellPolytopes.ActiveSetStorage.

For a full exact correlation tensor, the certificate concerns the finite scenario. The number of parties and the marginal convention are inferred from the active set:

p_exact = [1//1 1//1; 1//1 -1//1]
res = bell_frank_wolfe(Float64.(p_exact); v0=0.49)
cert = dyadic_certificate(res, p_exact; v0=49//100)
cert.visibility_lower  # Rational{BigInt}

For a bipartite Gram target, round the measurement directions and include their exact convex-hull shrinking factors in the certificate:

res = bell_frank_wolfe(v*v'; v0=0.69)
cert = dyadic_certificate(res; measurements=v, v0=69//100)

For other states, including multipartite states and marginals, provide target=vertices -> exact_correlations(vertices). This function receives a tuple of the rounded rational measurement matrices and must return the full exact correlation tensor, with identities last on each axis. The geometry form compensates marginal blocks before applying white-noise shrinking. Its finite_visibility refers to that compensated finite target; visibility_lower also includes shrinking. Exact noise tensors and certified local radii are supported by the finite-target form through o and radius.

An optional q supplies consecutive orbit labels for a uniform symmetry average. The caller must guarantee that this averaging preserves locality; no orbit verification or symmetry discovery is performed. Without q, the raw stored strategies are used. Floating-point targets, visibility values, and shrinking estimates are not accepted as exact certificates. The docstrings describe rounding precision, reusable squared shrinking bounds (shr2), and the backend=:integer alternative to exact integer products computed by BLAS.

Going further

More examples can be found in the corresponding folder of the package. They include the construction of a Bell inequality with a higher tolerance to noise as CHSH as well as multipartite and high-dimensional instances.