Advanced

In this section we will consider some advanced aspects related to the package.

Solving the dual versus the primal formulation of the SDP

For semidefinite programs that appear often in causal compatibility problems, using the dual formulation speeds up the solve time and significantly lowers RAM usage.

Consider the following example, where we use the MOSEK Fusion API to solve the primal version of a program, and then the dual:

[1]:

from inflation import InflationProblem, InflationSDP
from time import time
import numpy as np

qtriangle = InflationProblem(dag={"rho_AB": ["A", "B"],
                                  "rho_BC": ["B", "C"],
                                  "rho_AC": ["A", "C"]},
                             outcomes_per_party=[2, 2, 2],
                             settings_per_party=[1, 1, 1],
                             inflation_level_per_source=[2, 2, 2])
sdprelax = InflationSDP(qtriangle, verbose=0)
sdprelax.generate_relaxation('npa2')

P_W = np.zeros((2, 2, 2, 1, 1, 1))
for a, b, c in np.ndindex((2, 2, 2)):
    if a + b + c == 1:
        P_W[a, b, c, 0, 0, 0] = 1 / 3

sdprelax.set_distribution(P_W)

time0 = time()
sdprelax.solve(dualise=False)
print("The primal formulation was solved in", time() - time0, "seconds.")
time0 = time()
sdprelax.solve(dualise=True)
print("The dual formulation was solved in", time() - time0, "seconds.")

SDP solver took 3.70 seconds
The primal formulation was solved in 3.823288917541504 seconds.
SDP solver took 0.46 seconds
The dual formulation was solved in 0.5328986644744873 seconds.

Notice that there is an order of magnitude difference between the primal and dual formulations of the same problem. This is not true for all problems, but for the semidefinite programming relaxations generated for causal compatibility, almost always the dual formulation is more efficient. This should be taken into account when attempting to solve a relaxation. In what follows, we recompile some useful information for different interfaces.

CVXPY. If you export the problem to CVXPY, the behaviour depends on the solver you choose to use. When choosing MOSEK, note that CVXPY dualises by default all continuous problems. There is no automatic dualisation option. There is no option to specify whether to solve the primal or dual problem. Thus if you wanted to solve the primal with MOSEK, you would need to write the dual formulation manually, which when dualised would solve the primal (it is not expected that the user will need to do this!).
PICOS 2.4. Picos supports dualisation with the dualise=True options flag. See this issue for more details.
YALMIP. Like CVXPY, YALMIP automatically dualises problems, however there is a flag, dualize, in sdpsettings to disable this feature if so desired.
MOSEK Fusion API. Our implementation of the semidefinite programming relaxation supports both the primal and dual formulations, as seen in the example above. This is done manually, as MOSEK Fusion API does not have functionality to change from the primal to the dual formulations.

Large scale problems

For solving large scale semidefinite programs, it is recommended to use the MOSEK Fusion API, as going through interfaces for conic problems, such as PICOS or CVXPY, usually has an overhead in the pre-processing state (for example, there can be a higher RAM usage in the preprocessing stage than when solving the problem, which can lead to out-of-memory errors). There does not seem to be such an overhead when using YALMIP. For using YALMIP, the user can export the problem to .dat-s format using InflationSDP.write_to_file(), and load it in MATLAB using YALMIP’s loadsdpafile.

For large problems, it is recommended to try using a first-order SDP solver, such as SCS, if using second-order SDP solvers, such as MOSEK, is too slow or too memory-consuming. To use SCS the problem needs to be exported to the user’s interface of choice and have SCS installed.

It is also worth considering using symmetries to block-diagonalise the semidefinite program. This can be done with RepLAB in MATLAB. Symmetries arising from inflation can be calculated with InflationSDP._discover_inflation_symmetries(), and they are encoded as permutations of the list of generating monomials which leave the SDP invariant. This then can be used in RepLAB to block-diagonalise the problem, such as in the following example.

Using RepLAB in MATLAB to block-diagonalise an inflation SDP

In the following example, we generate the SDP relaxation for a problem, in this case, optimisation of the Mermin inequality over the triangle scenario with quantum sources, and show how to use RepLAB to block diagonalise the SDP and solve it in MATLAB.

First, we generate the SDP relaxation and write it to file. We generate a particularly large SDP (moment matrix with around 1700 columns) to showcase the advantages of block diagonalisation.

[2]:

from inflation import InflationProblem, InflationSDP
from scipy.io import savemat

scenario = InflationProblem(dag={"rho_AB": ["A", "B"],
                                 "rho_BC": ["B", "C"],
                                 "rho_AC": ["A", "C"]},
                           outcomes_per_party=(2, 2, 2),
                           settings_per_party=(2, 2, 2),
                           inflation_level_per_source=(3, 3, 3))
sdp = InflationSDP(scenario, commuting=False, verbose=1)
sdp.generate_relaxation("npa2")

mmnts = sdp.measurements
A0, B0, C0, A1, B1, C1 = (1 - 2*mmnts[party][0][setting][0]  # As correlators
                          for setting in range(2) for party in range(3))
sdp.set_objective(A1*B0*C0 + A0*B1*C0 + A0*B0*C1 - A1*B1*C1)

sdp.write_to_file('inflation_sdp.mat')
print("Writing the inflation symmetries to inflation_symmetries.mat")
savemat('inflation_symmetries.mat',
        {'inflation_symmetries': sdp.inflation_symmetries + 1} # because matlab indexing starts at 1
        )

Number of columns in the moment matrix: 1729

Calculating moment matrix: 100%|██████████| 1495585/1495585 [04:03<00:00, 6131.07it/s]
Applying symmetries      : 100%|██████████| 215/215 [00:00<00:00, 770.54it/s]

Number of variables after symmetrization: 6476

Initializing monomials   : 100%|██████████| 6476/6476 [00:08<00:00, 762.65it/s]

Writing the SDP program to inflation_sdp.mat
Writing the inflation symmetries to inflation_symmetries.mat

Within a MATLAB script, after installing RepLAB, we need to initialise it by running replab_init:

run "X:\...\replab-develop\replab_init.m"

Next, we load the SDP and symmetries from file, and create a Yalmip sdpvar for the moment matrix:

load('inflation_sdp.mat');
load('inflation_symmetries.mat');

% Convert symmetries to cell for compatibility with RepLAB functions
symmetries = num2cell(double(inflation_symmetries), 2);

% Build momentmatrix as a Yalmip sdpvar
IndexMatrix = double(momentmatrix);  % int32 -> double
vars = sdpvar(1, max(IndexMatrix(:)));
Gamma = vars(IndexMatrix);

Computing the block structure with RepLAB is as easy as running:

GammaBlocks = replab.CommutantVar.fromSymSdpMatrix(Gamma, symmetries);

which generates a block diagonal moment matrix (in about 1.5 minutes):

>> disp(GammaBlocks)
Commutant variable 1729x1729 (23 blocks, 6476 scalar variables)
           dim: 1729
   dimensions1: [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 8]
         field: 'real'
    matrixType: 'symmetric'
multiplicities: [70, 3, 3, 3, 14, 14, 14, 12, 12, 12, 12, 12, 12, 68, 68, 68, 4, 4, 4, 60, 60, 60, 36]
         types: 'RRRRRRRRRRRRRRRRRRRRRRR'
             U: 1729 x 1729 double
        blocks: [70, 3, 3, 3, 14, 14, 14, 12, 12, 12, 12, 12, 12, 68, 68, 68, 4, 4, 4, 60, 60, 60, 36]

Now we solve the SDP with Yalmip:

constraints = [GammaBlocks(1,1) == 1, GammaBlocks >= 0];
mermin = [vars(objective(:, 1))] * objective(:, 2);
optimize(constraints, mermin);

While solving the SDP without block diagonalisation takes about 4.5 minutes, the block diagonalised SDP takes about 10 seconds, a 25x speedup!

While this is a significant speedup, it should be noted that this is notable only for problems with a large number of symmetries, such as the previous example which was an inflation with 3 copies of each source. For problems with moment matrices of similar dimension but with less symmetries, the speedup is present but not as significant. For example, the same problem but with (2,2,2) copies per source and NPA level 2 union local level 1 generates a moment matrix of 873 columns. Solving the original SDP on one computer takes around 2.3 minutes, while solving the block diagonalised SDP takes around 1.7 minutes, giving a speedup of around 1.3x. While still relevant, it is not as significant as the previous example.