npa_hierarchy

Generates the NPA constraints.

npa_constraints

npa_constraints(
    assemblage: dict[tuple[int, int], Variable],
    k: int | str = 1,
    referee_dim: int = 1,
    no_signaling: bool = True,
) -> list[Constraint]

Generate the constraints specified by the NPA hierarchy up to a finite level.

¹

You can determine the level of the hierarchy by a positive integer or a string of a form like "1+ab+aab", which indicates that an intermediate level of the hierarchy should be used, where this example uses all products of 1 measurement, all products of one Alice and one Bob measurement, and all products of two Alice and one Bob measurement.

The commuting measurement assemblage operator must be given as a dictionary. The keys are tuples of Alice and Bob questions \(x, y\) and the values are cvxpy Variables which are matrices with entries:

\[ K_{xy}\Big(i + a \cdot dim_R, j + b \cdot dim_R \Big) = \langle i| \text{Tr}_{\mathcal{H}} \Big( \big( I_R \otimes A_a^x B_b^y \big) \sigma \Big) |j \rangle \]

Parameters:

• assemblage (dict[tuple[int, int], Variable]) –

  The commuting measurement assemblage operator.

• k (int | str, default: 1 ) –

  The level of the NPA hierarchy to use (default=1).

• referee_dim (int, default: 1 ) –

  The dimension of the referee's quantum system (default=1).

• no_signaling (bool, default: True ) –

  bool

Returns:

• list[Constraint] –

  A list of cvxpy constraints.

References

¹ Navascués, Miguel and Pironio, Stefano and Acín, Antonio. A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations. New Journal of Physics. vol. 10(7). (2008). doi:10.1088/1367-2630/10/7/073013.

Source code in toqito/state_opt/npa_hierarchy.py

def npa_constraints(
    assemblage: dict[tuple[int, int], cvxpy.Variable], k: int | str = 1, referee_dim: int = 1, no_signaling: bool = True
) -> list[cvxpy.constraints.constraint.Constraint]:
    r"""Generate the constraints specified by the NPA hierarchy up to a finite level.

    [@navascues2008convergent]

    You can determine the level of the hierarchy by a positive integer or a string
    of a form like "1+ab+aab", which indicates that an intermediate level of the hierarchy
    should be used, where this example uses all products of 1 measurement, all products of
    one Alice and one Bob measurement, and all products of two Alice and one Bob measurement.

    The commuting measurement assemblage operator must be given as a dictionary. The keys are
    tuples of Alice and Bob questions \(x, y\) and the values are cvxpy Variables which
    are matrices with entries:

    \[
        K_{xy}\Big(i + a \cdot dim_R, j + b \cdot dim_R \Big) =
        \langle i| \text{Tr}_{\mathcal{H}} \Big( \big(
            I_R \otimes A_a^x B_b^y \big) \sigma \Big) |j \rangle
    \]

    Args:
        assemblage: The commuting measurement assemblage operator.
        k: The level of the NPA hierarchy to use (default=1).
        referee_dim: The dimension of the referee's quantum system (default=1).
        no_signaling: bool

    Returns:
        A list of cvxpy constraints.

    """
    a_out, a_in, b_out, b_in = _get_nonlocal_game_params(assemblage, referee_dim)

    words = _gen_words(k, a_out, a_in, b_out, b_in)
    dim = len(words)

    if dim == 0:
        # Should not happen if IDENTITY_SYMBOL is always included
        raise ValueError("Generated word list is empty. Check _gen_words logic.")

    # Moment matrix (Gamma matrix in [@navascues2008convergent])
    # moment_matrix_R block corresponds to E[S_i^dagger S_j]
    moment_matrix_R = cvxpy.Variable((referee_dim * dim, referee_dim * dim), hermitian=True, name="R")

    # Referee's effective state rho_R = E[I]
    # This is the (0,0) block of moment_matrix_R since words[0] is Identity
    rho_R_referee = moment_matrix_R[0:referee_dim, 0:referee_dim]

    # Ensure rho_R_referee is a valid quantum state
    constraints = [
        cvxpy.trace(rho_R_referee) == 1,
        moment_matrix_R >> 0,
        # rho_R_referee >> 0 holds since it is a minor of moment_matrix_R
    ]

    # Store relations for (S_i^dagger S_j) -> block_index in moment_matrix_R
    # This helps enforce Γ(S_i^dagger S_j) = Γ(S_k^dagger S_l) if products are algebraically equal
    seen_reduced_products = {}

    for i in range(dim):
        for j in range(i, dim):  # Iterate over upper triangle including diagonal
            word_i_conj = tuple(reversed(words[i]))  # S_i^dagger

            # The product S_i^dagger S_j
            # For _reduce, ensure no IDENTITY_SYMBOL unless it's the only element.
            # If word_i_conj is (ID,), S_i_dagger_S_j is S_j. If word_j is (ID,), it's S_i_dagger.
            # If both are (ID,), product is (ID,).

            product_unreduced = []
            if word_i_conj != (IDENTITY_SYMBOL,):
                product_unreduced.extend(list(word_i_conj))
            if words[j] != (IDENTITY_SYMBOL,):
                product_unreduced.extend(list(words[j]))

            # This happens if both words[i] and words[j] were IDENTITY_SYMBOL
            if not product_unreduced:
                product_S_i_adj_S_j = (IDENTITY_SYMBOL,)
            else:
                product_S_i_adj_S_j = _reduce(tuple(product_unreduced))

            # Moment matrix (Gamma matrix in NPA paper [@navascues2008convergent] - arXiv:0803.4290)
            # This hierarchy can be generalized, e.g., to incorporate referee systems
            # as seen in extended nonlocal games (see, e.g., F. Speelman's thesis, [@speelman2016position]).
            current_block = moment_matrix_R[
                i * referee_dim : (i + 1) * referee_dim, j * referee_dim : (j + 1) * referee_dim
            ]

            if _is_zero(product_S_i_adj_S_j):  # Product is algebraically zero
                constraints.append(current_block == 0)
            elif _is_identity(product_S_i_adj_S_j):  # Product is identity operator
                # This occurs for (i,j) where S_i^dagger S_j = I. e.g. S_i = S_j and S_i is unitary (proj).
                # Or i=0, j=0 (I^dagger I = I).
                # This means current_block should be rho_R_referee if product_S_i_adj_S_j is I
                constraints.append(current_block == rho_R_referee)

            # Product is A_a^x B_b^y
            elif _is_meas(product_S_i_adj_S_j):
                alice_symbol, bob_symbol = product_S_i_adj_S_j
                constraints.append(
                    current_block
                    == assemblage[alice_symbol.question, bob_symbol.question][
                        alice_symbol.answer * referee_dim : (alice_symbol.answer + 1) * referee_dim,
                        bob_symbol.answer * referee_dim : (bob_symbol.answer + 1) * referee_dim,
                    ]
                )
            # Product is A_a^x or B_b^y (i.e., only one player involved)
            elif _is_meas_on_one_player(product_S_i_adj_S_j):  # Product is A_a^x or B_b^y
                symbol = product_S_i_adj_S_j[0]
                if symbol.player == "Alice":
                    # Sum over Bob's outcomes for a fixed Bob question (e.g., y=0)
                    # E[A_a^x] = sum_b K_x0(a,b)
                    sum_over_bob_outcomes = sum(
                        assemblage[symbol.question, 0][  # Assuming y=0 for Bob's marginal
                            symbol.answer * referee_dim : (symbol.answer + 1) * referee_dim,
                            b_ans * referee_dim : (b_ans + 1) * referee_dim,
                        ]
                        for b_ans in range(b_out)
                    )
                    constraints.append(current_block == sum_over_bob_outcomes)
                else:  # elif symbol.player == "Bob":
                    # Sum over Alice's outcomes for a fixed Alice question (e.g., x=0)
                    # E[B_b^y] = sum_a K_0y(a,b)
                    sum_over_alice_outcomes = sum(
                        assemblage[0, symbol.question][  # Assuming x=0 for Alice's marginal
                            a_ans * referee_dim : (a_ans + 1) * referee_dim,
                            symbol.answer * referee_dim : (symbol.answer + 1) * referee_dim,
                        ]
                        for a_ans in range(a_out)
                    )
                    constraints.append(current_block == sum_over_alice_outcomes)
            elif product_S_i_adj_S_j in seen_reduced_products:
                # This product S_k has been seen before as S_p^dagger S_q
                # So, Γ(S_i, S_j) = Γ(S_p, S_q)
                prev_i, prev_j = seen_reduced_products[product_S_i_adj_S_j]
                # Make sure to get the upper triangle element if current (i,j) is lower
                # The prev_i, prev_j should always refer to an upper triangle element by construction.
                previous_block = moment_matrix_R[
                    prev_i * referee_dim : (prev_i + 1) * referee_dim, prev_j * referee_dim : (prev_j + 1) * referee_dim
                ]
                constraints.append(current_block == previous_block)
            else:
                # First time seeing this operator product S_k
                seen_reduced_products[product_S_i_adj_S_j] = (i, j)

    # Constraints on the assemblage K_xy(a,b) itself --always apply all of these constraints!
    for x_alice_in in range(a_in):
        for y_bob_in in range(b_in):
            # Positivity: K_xy(a,b) >= 0 (operator PSD)
            for a_alice_out in range(a_out):
                for b_bob_out in range(b_out):
                    assemblage_block = assemblage[x_alice_in, y_bob_in][
                        a_alice_out * referee_dim : (a_alice_out + 1) * referee_dim,
                        b_bob_out * referee_dim : (b_bob_out + 1) * referee_dim,
                    ]
                    constraints.append(assemblage_block >> 0)

            # Normalization: Sum_{a,b} K_xy(a,b) = rho_R
            sum_over_outcomes_ab = sum(
                assemblage[x_alice_in, y_bob_in][
                    a * referee_dim : (a + 1) * referee_dim, b * referee_dim : (b + 1) * referee_dim
                ]
                for a in range(a_out)
                for b in range(b_out)
            )
            constraints.append(sum_over_outcomes_ab == rho_R_referee)
    if no_signaling:
        # No-signaling constraints on assemblage - ALWAYS APPLY
        # Bob's marginal rho_B(b|y) = Sum_a K_xy(a,b) must be independent of x
        for y_bob_in in range(b_in):
            for b_bob_out in range(b_out):
                sum_over_a_for_x0 = sum(
                    assemblage[0, y_bob_in][
                        a * referee_dim : (a + 1) * referee_dim, b_bob_out * referee_dim : (b_bob_out + 1) * referee_dim
                    ]
                    for a in range(a_out)
                )
                for x_alice_in in range(1, a_in):
                    sum_over_a_for_x_current = sum(
                        assemblage[x_alice_in, y_bob_in][
                            a * referee_dim : (a + 1) * referee_dim,
                            b_bob_out * referee_dim : (b_bob_out + 1) * referee_dim,
                        ]
                        for a in range(a_out)
                    )
                    constraints.append(sum_over_a_for_x0 == sum_over_a_for_x_current)

        # Alice's marginal rho_A(a|x) = Sum_b K_xy(a,b) must be independent of y
        for x_alice_in in range(a_in):
            for a_alice_out in range(a_out):  # For each Alice outcome a
                sum_over_b_for_y0 = sum(
                    assemblage[x_alice_in, 0][
                        a_alice_out * referee_dim : (a_alice_out + 1) * referee_dim,
                        b * referee_dim : (b + 1) * referee_dim,
                    ]
                    for b in range(b_out)
                )
                for y_bob_in in range(1, b_in):
                    sum_over_b_for_y_current = sum(
                        assemblage[x_alice_in, y_bob_in][
                            a_alice_out * referee_dim : (a_alice_out + 1) * referee_dim,
                            b * referee_dim : (b + 1) * referee_dim,
                        ]
                        for b in range(b_out)
                    )
                    constraints.append(sum_over_b_for_y0 == sum_over_b_for_y_current)

    return constraints

bell_npa_constraints

bell_npa_constraints(
    p_var: Variable, desc: list[int], k: int | str = 1
) -> list[Constraint]

Generate NPA hierarchy constraints for Bell inequalities ¹.

The constraints are based on the positivity of a moment matrix constructed from measurement operators. This function generates constraints for a CVXPY variable representing probabilities or correlations in the Collins-Gisin notation. ²

The level of the hierarchy k can be an integer (standard NPA level) or a string specifying intermediate levels (e.g., "1+ab", "2+aab").

The input p_var is a CVXPY variable representing the probabilities in the Collins-Gisin (CG) notation. It should have dimensions \(((oa-1) \times ma+1, (ob-1) \times mb+1)\), where \(oa, ob\) are the number of outputs and \(ma, mb\) are the number of inputs for Alice and Bob, respectively, as specified in desc = [\(oa\), \(ob\), \(ma\), \(mb\)]. The entries of p_var correspond to: - p_var[0, 0]: The overall probability (should be 1). - p_var[i, 0] (for \(i > 0\)): Marginal probabilities/correlations for Alice. - p_var[0, j] (for \(j > 0\)): Marginal probabilities/correlations for Bob. - p_var[i, j] (for \(i > 0, j > 0\)): Joint probabilities/correlations for Alice and Bob.

The mapping from indices \((i, j)\) to specific operators depends on the ordering defined by desc. Specifically, if \(i = (oa-1) \times x + a + 1\) and \(j = (ob-1) \times y + b + 1\)

• p_var[i, 0] corresponds to the expectation of Alice's operator \(A_{a|x}\) (using \(0\) to \(oa-2\) for \(a\)).
• p_var[0, j] corresponds to the expectation of Bob's operator \(B_{b|y}\) (using \(0\) to \(ob-2\) for \(b\)).
• p_var[i, j] corresponds to the expectation of the product \(A_{a|x} B_{b|y}\).

Parameters:

• p_var (Variable) –

  A CVXPY Variable representing probabilities/correlations in Collins-Gisin notation. Shape: \(((oa-1) \times ma+1, (ob-1) \times mb+1)\).

• desc (list[int]) –

  A list [\(oa\), \(ob\), \(ma\), \(mb\)] specifying outputs and inputs for Alice and Bob.

• k (int | str, default: 1 ) –

  The level of the NPA hierarchy (integer or string like "1+ab"). Default is 1.

Returns:

• list[Constraint] –

  A list of CVXPY constraints.

Raises:

• ValueError –

  If internal identity mapping fails.

Examples:

Consider the CHSH inequality scenario with desc = [2, 2, 2, 2]. We want to generate the NPA level 1 constraints.

import cvxpy
import numpy as np
from toqito.state_opt.npa_hierarchy import bell_npa_constraints
desc = [2, 2, 2, 2]
oa, ob, ma, mb = desc
p_var_dim = ((oa - 1) * ma + 1, (ob - 1) * mb + 1)
p_var = cvxpy.Variable(p_var_dim, name="p_cg")
constraints = bell_npa_constraints(p_var, desc, k=1)
print(len(constraints))
print(constraints[0])

14
Gamma + Promote(-0.0, (5, 5)) >> 0

We can also use intermediate levels, like "1+ab":

constraints_1ab = bell_npa_constraints(p_var, desc, k="1+ab")
print(len(constraints_1ab))
print(constraints_1ab[0])

34
Gamma + Promote(-0.0, (9, 9)) >> 0

For the CGLMP inequality with dim=3, desc = [3, 3, 2, 2], level 1:

import cvxpy
import numpy as np
from toqito.state_opt.npa_hierarchy import bell_npa_constraints
desc_cglmp = [3, 3, 2, 2]
oa_c, ob_c, ma_c, mb_c = desc_cglmp
p_var_dim_c = ((oa_c - 1) * ma_c + 1, (ob_c - 1) * mb_c + 1)
p_var_c = cvxpy.Variable(p_var_dim_c, name="p_cglmp")
constraints_c = bell_npa_constraints(p_var_c, desc_cglmp, k=1)
print(len(constraints_c))
print(constraints_c[0])

38
Gamma + Promote(-0.0, (9, 9)) >> 0

References

¹ Navascués, Miguel and Pironio, Stefano and Acín, Antonio. A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations. New Journal of Physics. vol. 10(7). (2008). doi:10.1088/1367-2630/10/7/073013.
² Collins, Daniel and Gisin, Nicolas. A relevant two qubit Bell inequality inequivalent to the CHSH inequality. Journal of Physics A: Mathematical and General. vol. 37(5). (2004). doi:10.1088/0305-4470/37/5/021.

Source code in toqito/state_opt/npa_hierarchy.py

def bell_npa_constraints(
    p_var: cvxpy.Variable,
    desc: list[int],
    k: int | str = 1,
) -> list[cvxpy.constraints.constraint.Constraint]:
    r"""Generate NPA hierarchy constraints for Bell inequalities [@navascues2008convergent].

    The constraints are based on the positivity of a moment matrix constructed from measurement
    operators. This function generates constraints for a CVXPY variable representing probabilities
    or correlations in the Collins-Gisin notation. [@collins2004relevant]

    The level of the hierarchy ``k`` can be an integer (standard NPA level) or a string specifying
    intermediate levels (e.g., "1+ab", "2+aab").

    The input ``p_var`` is a CVXPY variable representing the probabilities in the Collins-Gisin (CG)
    notation. It should have dimensions \(((oa-1) \times ma+1, (ob-1) \times mb+1)\),
    where \(oa, ob\) are the number of outputs and \(ma, mb\) are the number of inputs for Alice
    and Bob, respectively, as specified in ``desc`` = [\(oa\), \(ob\), \(ma\), \(mb\)].
    The entries of ``p_var`` correspond to:
    - ``p_var[0, 0]``: The overall probability (should be 1).
    - ``p_var[i, 0]`` (for \(i > 0\)): Marginal probabilities/correlations for Alice.
    - ``p_var[0, j]`` (for \(j > 0\)): Marginal probabilities/correlations for Bob.
    - ``p_var[i, j]`` (for \(i > 0, j > 0\)): Joint probabilities/correlations for Alice and Bob.

    The mapping from indices \((i, j)\) to specific operators depends on the ordering defined by ``desc``.
    Specifically, if \(i = (oa-1) \times x + a + 1\) and \(j = (ob-1) \times y + b + 1\)

    - ``p_var[i, 0]`` corresponds to the expectation of Alice's operator
                      \(A_{a|x}\) (using \(0\) to \(oa-2\) for \(a\)).
    - ``p_var[0, j]`` corresponds to the expectation of Bob's operator
                      \(B_{b|y}\) (using \(0\) to \(ob-2\) for \(b\)).
    - ``p_var[i, j]`` corresponds to the expectation of the product \(A_{a|x} B_{b|y}\).

    Args:
        p_var: A CVXPY Variable representing probabilities/correlations in Collins-Gisin notation.
                  Shape: \(((oa-1) \times ma+1, (ob-1) \times mb+1)\).
        desc: A list [\(oa\), \(ob\), \(ma\), \(mb\)]
                    specifying outputs and inputs for Alice and Bob.
        k: The level of the NPA hierarchy (integer or string like "1+ab"). Default is 1.


    Returns:
        A list of CVXPY constraints.

    Raises:
        ValueError: If internal identity mapping fails.

    Examples:
        Consider the CHSH inequality scenario with ``desc = [2, 2, 2, 2]``. We want to generate the NPA level 1
        constraints.

        ```python exec="1" source="above" result="text" session="npa_example"
        import cvxpy
        import numpy as np
        from toqito.state_opt.npa_hierarchy import bell_npa_constraints
        desc = [2, 2, 2, 2]
        oa, ob, ma, mb = desc
        p_var_dim = ((oa - 1) * ma + 1, (ob - 1) * mb + 1)
        p_var = cvxpy.Variable(p_var_dim, name="p_cg")
        constraints = bell_npa_constraints(p_var, desc, k=1)
        print(len(constraints))
        print(constraints[0])
        ```

        We can also use intermediate levels, like "1+ab":

        ```python exec="1" source="above" result="text" session="npa_example"
        constraints_1ab = bell_npa_constraints(p_var, desc, k="1+ab")
        print(len(constraints_1ab))
        print(constraints_1ab[0])
        ```

        For the CGLMP inequality with ``dim=3``, ``desc = [3, 3, 2, 2]``, level 1:

        ```python exec="1" source="above" result="text"
        import cvxpy
        import numpy as np
        from toqito.state_opt.npa_hierarchy import bell_npa_constraints
        desc_cglmp = [3, 3, 2, 2]
        oa_c, ob_c, ma_c, mb_c = desc_cglmp
        p_var_dim_c = ((oa_c - 1) * ma_c + 1, (ob_c - 1) * mb_c + 1)
        p_var_c = cvxpy.Variable(p_var_dim_c, name="p_cglmp")
        constraints_c = bell_npa_constraints(p_var_c, desc_cglmp, k=1)
        print(len(constraints_c))
        print(constraints_c[0])
        ```


    """
    oa, ob, ma, mb = desc
    words = _gen_words(k, oa, ma, ob, mb)
    dim = len(words)
    gamma = cvxpy.Variable((dim, dim), name="Gamma", symmetric=True)
    constraints = [gamma >> 0]
    p_flat = p_var.flatten(order="F")
    seen_constraints = {}

    constraints.append(gamma[0, 0] == p_var[0, 0])

    seen_constraints[()] = (0, 0)

    for i in range(dim):
        for j in range(i, dim):
            if i == 0 and j == 0:
                continue
            word_i = words[i]
            word_j = words[j]
            word_i_conj = tuple(reversed(word_i))
            combined_word = _reduce(word_i_conj + word_j)

            if not combined_word:
                constraints.append(gamma[i, j] == 0)
                continue

            constraint_key = combined_word
            if constraint_key in seen_constraints:
                prev_i, prev_j = seen_constraints[constraint_key]
                constraints.append(gamma[i, j] == gamma[prev_i, prev_j])
                continue

            p_cg_index = _word_to_p_cg_index(combined_word, oa, ob, ma, mb)
            if p_cg_index is not None:
                constraints.append(gamma[i, j] == p_flat[p_cg_index])
                seen_constraints[constraint_key] = (i, j)
            else:
                seen_constraints[constraint_key] = (i, j)

    return constraints