has_symmetric_extension

Determine whether there exists a symmetric extension for a given quantum state.

has_symmetric_extension

has_symmetric_extension(rho: ndarray, level: int = 2, dim: ndarray | int | None = None, ppt: bool = True, tol: float = 0.0001) -> bool

Determine whether there exists a symmetric extension for a given quantum state.

For more information, see ¹.

Determining whether an operator possesses a symmetric extension at some level level can be used as a check to determine if the operator is entangled or not.

This function was adapted from QETLAB.

Examples:

2-qubit symmetric extension:

In ², it was shown that a 2-qubit state \(\rho_{AB}\) has a symmetric extension if and only if

\[ \text{Tr}(\rho_B^2) \geq \text{Tr}(\rho_{AB}^2) - 4 \sqrt{\text{det}(\rho_{AB})}. \]

This closed-form equation is much quicker to check than running the semidefinite program.

import numpy as np
from toqito.state_props import has_symmetric_extension
from toqito.matrix_ops import partial_trace
rho = np.array([[1, 0, 0, -1], [0, 1, 1/2, 0], [0, 1/2, 1, 0], [-1, 0, 0, 1]])
# Show the closed-form equation holds
print(
np.trace(np.linalg.matrix_power(partial_trace(rho, 1), 2))
>= np.trace(rho**2) - 4 * np.sqrt(np.linalg.det(rho)))

True

# Now show that the `has_symmetric_extension` function recognizes this case.
print(has_symmetric_extension(rho))

True

Higher qubit systems:

Consider a density operator corresponding to one of the Bell states.

\[ \rho = \frac{1}{2} \begin{pmatrix} 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 \end{pmatrix} \]

To make this state over more than just two qubits, let's construct the following state

\[ \sigma = \rho \otimes \rho. \]

As the state \(\sigma\) is entangled, there should not exist a symmetric extension at some level. We see this being the case for a relatively low level of the hierarchy.

import numpy as np
from toqito.states import bell
from toqito.state_props import has_symmetric_extension
rho = bell(0) @ bell(0).conj().T
sigma = np.kron(rho, rho)
print(has_symmetric_extension(sigma))

True

Raises:

• ValueError –

  If dimension does not evenly divide matrix length.

Parameters:

• rho (ndarray) –

  A matrix or vector.

• level (int, default: 2 ) –

  Level of the hierarchy to compute.

• dim (ndarray | int | None, default: None ) –

  The default has both subsystems of equal dimension.

• ppt (bool, default: True ) –

  If True, this enforces that the symmetric extension must be PPT.

• tol (float, default: 0.0001 ) –

  Tolerance when determining whether a symmetric extension exists.

Returns:

• bool –

  True if mat has a symmetric extension; False otherwise.

References

¹ Doherty, A. and Parrilo, Pablo and Spedalieri, Federico. Distinguishing Separable and Entangled States. Physical Review Letters. vol. 88(18). (2002). doi:10.1103/physrevlett.88.187904.
² Chen, Jianxin and Ji, Zhengfeng and Kribs, David and Lütkenhaus, Norbert and Zeng, Bei. Symmetric extension of two-qubit states. Physical Review A. vol. 90(3). (2014). doi:10.1103/physreva.90.032318.

Source code in toqito/state_props/has_symmetric_extension.py

def has_symmetric_extension(
    rho: np.ndarray,
    level: int = 2,
    dim: np.ndarray | int | None = None,
    ppt: bool = True,
    tol: float = 1e-4,
) -> bool:
    r"""Determine whether there exists a symmetric extension for a given quantum state.

    For more information, see [@Doherty_2002_Distinguishing].

    Determining whether an operator possesses a symmetric extension at some level `level`
    can be used as a check to determine if the operator is entangled or not.

    This function was adapted from QETLAB.

    Examples:
        2-qubit symmetric extension:

        In [@Chen_2014_Symmetric], it was shown that a 2-qubit state \(\rho_{AB}\) has a
        symmetric extension if and only if

        \[
            \text{Tr}(\rho_B^2) \geq \text{Tr}(\rho_{AB}^2) - 4 \sqrt{\text{det}(\rho_{AB})}.
        \]

        This closed-form equation is much quicker to check than running the semidefinite program.

        ```python exec="1" source="above" session="has_symmetric_example"
        import numpy as np
        from toqito.state_props import has_symmetric_extension
        from toqito.matrix_ops import partial_trace
        rho = np.array([[1, 0, 0, -1], [0, 1, 1/2, 0], [0, 1/2, 1, 0], [-1, 0, 0, 1]])
        # Show the closed-form equation holds
        print(
        np.trace(np.linalg.matrix_power(partial_trace(rho, 1), 2))
        >= np.trace(rho**2) - 4 * np.sqrt(np.linalg.det(rho)))
        ```

        ```python exec="1" source="above" session="has_symmetric_example"
        # Now show that the `has_symmetric_extension` function recognizes this case.
        print(has_symmetric_extension(rho))
        ```

        Higher qubit systems:

        Consider a density operator corresponding to one of the Bell states.

        \[
            \rho = \frac{1}{2} \begin{pmatrix}
                                1 & 0 & 0 & 1 \\
                                0 & 0 & 0 & 0 \\
                                0 & 0 & 0 & 0 \\
                                1 & 0 & 0 & 1
                               \end{pmatrix}
        \]

        To make this state over more than just two qubits, let's construct the following state

        \[
            \sigma = \rho \otimes \rho.
        \]

        As the state \(\sigma\) is entangled, there should not exist a symmetric extension at some
        level. We see this being the case for a relatively low level of the hierarchy.

        ```python exec="1" source="above"
        import numpy as np
        from toqito.states import bell
        from toqito.state_props import has_symmetric_extension
        rho = bell(0) @ bell(0).conj().T
        sigma = np.kron(rho, rho)
        print(has_symmetric_extension(sigma))
        ```

    Raises:
        ValueError: If dimension does not evenly divide matrix length.

    Args:
        rho: A matrix or vector.
        level: Level of the hierarchy to compute.
        dim: The default has both subsystems of equal dimension.
        ppt: If `True`, this enforces that the symmetric extension must be PPT.
        tol: Tolerance when determining whether a symmetric extension exists.

    Returns:
        `True` if `mat` has a symmetric extension; `False` otherwise.

    """
    len_mat = rho.shape[1]

    # Set default dimension if none was provided.
    if dim is None:
        dim_val = int(np.round(np.sqrt(len_mat)))
    elif isinstance(dim, int):
        dim_val = dim
    else:
        dim_val = None

    # Allow the user to enter in a single integer for dimension.
    if dim_val is not None:
        dim_arr = np.array([dim_val, len_mat / dim_val])
        if np.abs(dim_arr[1] - np.round(dim_arr[1])) >= 2 * len_mat * np.finfo(float).eps:
            raise ValueError("If `dim` is a scalar, it must evenly divide the length of the matrix.")
        dim_arr[1] = int(np.round(dim_arr[1]))
    else:
        dim_arr = np.array(dim)

    dim_arr = np.int_(dim_arr)

    dim_x, dim_y = int(dim_arr[0]), int(dim_arr[1])
    # In certain situations, we don't need semidefinite programming.
    if level == 1 or len_mat <= 6 and ppt:
        if not ppt:
            # In some cases, the problem is *really* trivial.
            return is_positive_semidefinite(rho)

        # In this case, all they asked for is a 1-copy PPT symmetric extension
        # (i.e., they're asking if the state is PPT).
        return is_ppt(rho, 2, dim_arr) and is_positive_semidefinite(rho)

    # In the 2-qubit case, an analytic formula is known for whether or not a state has a
    # (2-copy, non-PPT) symmetric extension that is much faster to use than semidefinite
    # programming [CJKLZB14]_.
    if level == 2 and not ppt and dim_x == 2 and dim_y == 2:
        return np.trace(np.linalg.matrix_power(partial_trace(rho, [0], [dim_x, dim_y]), 2)) >= np.trace(
            np.linalg.matrix_power(rho, 2)
        ) - 4 * np.sqrt(np.linalg.det(rho))

    # Otherwise, use semidefinite programming to find a symmetric extension.
    # We solve a feasibility SDP: find sigma on X ⊗ Y^⊗level such that
    # tr_{Y_2,...,Y_level}(sigma) = rho, sigma >= 0, sigma is symmetric
    # under permutations of Y copies, and (optionally) PPT constraints hold.
    dim_list = np.int_([dim_x] + [dim_y] * level)
    sys_list = list(range(2, 2 + level - 1))
    sym = symmetric_projection(dim_y, level)
    dim_total = int(np.prod(dim_list))

    sigma = cvxpy.Variable((dim_total, dim_total), hermitian=True)

    constraints = [
        partial_trace(sigma, sys_list, dim_list) == rho,
        sigma >> 0,
        np.kron(np.identity(dim_x), sym) @ sigma @ np.kron(np.identity(dim_x), sym) == sigma,
    ]

    if ppt:
        constraints.append(partial_transpose(sigma, [0], dim_list) >> 0)
        for sys in range(level - 1):
            constraints.append(partial_transpose(sigma, [sys + 2], dim_list) >> 0)

    problem = cvxpy.Problem(cvxpy.Minimize(0), constraints)
    problem.solve()

    return problem.status == "optimal"