Get linear equality constraints working.

[eliewolfe]

This means both automatically EXTRACTING the implicit equality constraints from the problem set, as well as accepting further equality constraints from the user. When I attempt to implement the implicit constraints I am getting "infeasible" without set_values. This is bad. Either the implicit constraints I am extracting are wrong, or they are being processed incorrectly. My suspicion is the latter.

Code to reproduce the example:

from causalinflation.InflationProblem import InflationProblem
from causalinflation.quantum.InflationSDP import InflationSDP
import numpy as np
prob = InflationProblem(dag={'U_AB': ['A', 'B'],
                             'A': ['B']},
                        outcomes_per_party=(2, 2),
                        settings_per_party=(3, 1),
                        inflation_level_per_source=(1,),
                        order=('A', 'B'),
                        verbose=2)
sdp = InflationSDP(prob, commuting=False, verbose=1)
sdp.generate_relaxation('local1')
for equality_dict in sdp.moment_linear_equalities:
    print(f"Equality dict: {equality_dict}")
import sympy
objective = sympy.Symbol('pAB(00|00)')
objective += sympy.Symbol('pA(1|0)') - sympy.Symbol('pAB(10|00)')  # sympy.Symbol('pAB(11|00)')
objective += sympy.Symbol('pAB(00|10)')
objective += sympy.Symbol('pAB(10|10)')
objective += sympy.Symbol('pA(0|2)') - sympy.Symbol('pAB(00|20)')  # sympy.Symbol('pAB(01|20)')
sdp.set_objective(objective, direction='max')
sdp.solve()
print(f"Objective value: {sdp.objective_value}")
for k, v in sdp.solution_object['x'].items():
    if sdp.compound_monomial_from_name_dict[k].knowable_q:
        print(f"{k} := {v}")

[ecboghiu]

When solving the primal instead of the dual, it seems to work giving output:

Solving the model...
Objective value: 2.207106776722439
pA(1|2) := 0.29281414078103335
pAB(00|20) := 0.2071900961706035
pAB(00|10) := 0.42678101644235916
pAB(10|10) := 0.4267740350186997
pAB(00|00) := 0.42678081525412187
pA(0|1) := 0.5000011234347465
pA(0|0) := 0.5000048718121373
pAB(10|20) := 0.1464078444377339
pA(1|0) := 0.49999512818786274
pAB(10|00) := 0.07321998122896768
pA(1|1) := 0.49999887656525355
pA(0|2) := 0.7071858592189666

This means for sure that there is a bug with how the dual problem handles the constraints. I'll look into it tomorrow.

[ecboghiu]

I also checked and all the equalities are correctly satisfied with the primal formulation, you can use the following code for that:

for equality_dict in sdp.moment_linear_equalities:
    summ = 0
    for k, v in equality_dict.items():
        if k != '1':
            summ += sdp.solution_object['x'][k] * v
        else:
            summ += v
    if not np.allclose(summ, 0):
        print(f"Equality dict: {equality_dict} Numerical evaluation: {summ}")

[ecboghiu]

Turns out there was a sign wrong :'( Now it should be working! https://github.com/ecboghiu/inflation/commit/19a4c384721e9a3797ea2e00f586c1572f221541 Reopen if other issues arise.

[eliewolfe]

REOPENING ISSUE. The first example of a non-network scenario in examples.ipynb, namely

prob = InflationProblem(dag={'U_AB': ['A', 'B'],
                             'U_AC': ['A', 'C'],
                             'U_AD': ['A', 'D'],
                             'C': ['D'],
                             'A': ['B', 'C', 'D']},
                        outcomes_per_party=(2, 2, 2, 2),
                        settings_per_party=(1, 1, 1, 1),
                        inflation_level_per_source=(1, 1, 1),
                        order=('A', 'B', 'C', 'D'),
                        verbose=2)
sdp = InflationSDP(prob, commuting=False, verbose=1)
sdp.generate_relaxation('local1')
dist = np.zeros((2, 2, 2, 2, 1, 1, 1, 1), dtype=float)
dist[0,0,0,0] = 1/4
dist[0,1,0,1] = 1/4
dist[1,0,0,0] = 1/4
dist[1,1,1,0] = 1/4
sdp.set_distribution(dist)
sdp.solve()

is leading to a KeyError: '0' in line 286 of solveSDP_MosekFUSION. I suppose this is a problem caused by zero-probability known_value monomials appearing inside the list of equality constraints? I can try to patch, but probably best if @ecboghiu takes lead.

[eliewolfe]

@ecboghiu I was thinking a bit more about this, and I'd like to petition you to set up solve_SDPMosekFusion without any CONSTANT_KEY. How about you simply DEFINE a new mask matrix (all zeroes) and associated with some internal key you call CONSTANT_KEY. Then, we dump all the numerical values in known_moments into it by taking the weighted sum of their mask matrices, weighted by their values. The advantage to doing this without assuming that 1 will be found in known_moments is that it makes it much much easier to adapted the code for supports testing, when the 1 monomial is found in lower_bounds instead of known_moments. Let me be clear: by using an internal CONSTANT_KEY in solve_SDPMosekFusion we get the following benefits:

1. Safer handling of zero-valued constants throughout.
2. No need to make constant_term a _keywordargument as was necessitated thus far in the for_supports branch.
3. No need to patch known_moments to contain a fake monomial assigned the value 1 with a blank mask matrix. Again, this is currently needed and implemented ad hoc in the for_supports branch. So, would you mind please doing it this way?

[ecboghiu]

Sure, this is sort of what I had in mind in this comment

https://github.com/ecboghiu/inflation/issues/65#issuecomment-1256014727

The only change would be changing CONSTANT_IDX from '1' to something that will never be the name of a monomial, such as 'CONSTANT', and then the other changes should be minimal.

[eliewolfe]

Appears to be fixed pursuant to Issue #72 . Closing for now; note that we still need to write tests for the equality constraints, as per issue #69 .