Advanced

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

Performance considerations

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(solve_dual=False)
print("The primal formulation was solved in", time()-time0, "seconds.")
time0 = time()
sdprelax.solve(solve_dual=True)
print("The dual formulation was solved in", time()-time0, "seconds.")

The primal formulation was solved in 20.820358276367188 seconds.
The dual formulation was solved in 0.8410844802856445 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._calculate_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 (see this example from RepLAB). A more in-depth example with code detailing this will be added to the Examples section in the future.