fix: schmidt_rank for unequal subsystem dimensions

toqito/state_props/is_separable.py

@@ -493,7 +493,11 @@ def is_separable(
     # if they are determined to be "weakly irreducible", as per Cariello [@cariello2013separability]
     # and QETLAB's implementation. This is distinct from this pure state check.)
     if state_rank == 1:
-        s_rank = schmidt_rank(current_state, dims_list)
+        # Recover the state vector and compute its Schmidt rank. Passing the density matrix itself
+        # to `schmidt_rank` would return the operator Schmidt rank, which for a pure state is the
+        # square of the Schmidt rank (e.g. 4 instead of 2 for a Bell state).
+        _, eigenvectors = np.linalg.eigh(current_state)
+        s_rank = schmidt_rank(eigenvectors[:, -1], dims_list)
         if s_rank == 1:
             return True, "pure state with Schmidt rank 1"
         return False, f"pure state with Schmidt rank {int(s_rank)} > 1"

toqito/state_props/schmidt_rank.py

@@ -96,7 +96,7 @@ def schmidt_rank(rho: np.ndarray, dim: int | list[int] | np.ndarray | None = Non
         dim = np.array([dim, len(rho) / dim], dtype=int)
         dim[1] = np.round(dim[1])
 
-    return np.linalg.matrix_rank(np.reshape(rho, dim[::-1]))
+    return np.linalg.matrix_rank(np.reshape(rho, dim))
 
 
 def _operator_schmidt_rank(rho: np.ndarray, dim: int | list[int] | np.ndarray | None = None) -> int | float:

toqito/state_props/tests/test_is_separable.py

@@ -2147,3 +2147,38 @@ def test_sorted_real_eigs_desc_returns_none_when_both_solvers_fail():
         mock.patch("numpy.linalg.eigvals", side_effect=np.linalg.LinAlgError("no convergence")),
     ):
         assert _sorted_real_eigs_desc(np.diag([0.5, 0.3, 0.2])) is None
+
+
+@pytest.mark.parametrize(
+    "vec, dims, expected_sep, expected_rank",
+    [
+        # Product state |0>|1>.
+        (np.kron(basis(2, 0), basis(2, 1)), [2, 2], True, 1),
+        # Bell state on 2x2.
+        (bell(0), [2, 2], False, 2),
+        # Maximally entangled state on 3x3.
+        (max_entangled(3, False, False), [3, 3], False, 3),
+        # Product state on 2x3 (unequal dimensions), previously reported as entangled.
+        (np.kron(np.array([[1], [2]]) / np.sqrt(5), np.array([[1], [2], [3]]) / np.sqrt(14)), [2, 3], True, 1),
+        # Partially entangled state on 2x3 with two non-zero Schmidt coefficients.
+        (
+            (np.sqrt(0.7) * np.kron(basis(2, 0), basis(3, 0)) + np.sqrt(0.3) * np.kron(basis(2, 1), basis(3, 2))),
+            [2, 3],
+            False,
+            2,
+        ),
+    ],
+)
+def test_pure_state_reports_schmidt_rank_not_operator_schmidt_rank(vec, dims, expected_sep, expected_rank):
+    """The pure-state branch reports the Schmidt rank of the vector, not the operator Schmidt rank.
+
+    The operator Schmidt rank of a pure state is the square of its Schmidt rank, so a Bell state was
+    previously reported as having "Schmidt rank 4".
+    """
+    rho = vec @ vec.conj().T
+    separable, reason = is_separable(rho, dims)
+    assert separable is expected_sep
+    if expected_sep:
+        assert reason == "pure state with Schmidt rank 1"
+    else:
+        assert reason == f"pure state with Schmidt rank {expected_rank} > 1"

toqito/state_props/tests/test_schmidt_rank.py

@@ -68,3 +68,37 @@ def test_schmidt_rank_one_dimensional_vector():
     """A 1D (non-column) vector input is handled (covers the vector branch)."""
     bell_1d = np.array([1, 0, 0, 1]) / np.sqrt(2)
     assert schmidt_rank(bell_1d) == 2
+
+
+@pytest.mark.parametrize(
+    "vec, dim, expected_result",
+    [
+        # Product states with unequal subsystem dimensions have Schmidt rank 1.
+        (np.kron(np.array([[1], [2]]) / np.sqrt(5), np.array([[1], [0], [0]])), [2, 3], 1),
+        (np.kron(np.array([[1], [0]]), np.array([[1], [2], [3]]) / np.sqrt(14)), [2, 3], 1),
+        (np.kron(np.array([[1], [2], [3]]) / np.sqrt(14), np.array([[1], [1j]]) / np.sqrt(2)), [3, 2], 1),
+        # (|00> + |11>) / sqrt(2) embedded in 2x3 and 3x2 has Schmidt rank 2.
+        ((np.kron([[1], [0]], [[1], [0], [0]]) + np.kron([[0], [1]], [[0], [1], [0]])) / np.sqrt(2), [2, 3], 2),
+        ((np.kron([[1], [0], [0]], [[1], [0]]) + np.kron([[0], [1], [0]], [[0], [1]])) / np.sqrt(2), [3, 2], 2),
+        # Integer `dim` gives the dimension of the first subsystem.
+        ((np.kron([[1], [0]], [[1], [0], [0]]) + np.kron([[0], [1]], [[0], [1], [0]])) / np.sqrt(2), 2, 2),
+    ],
+)
+def test_schmidt_rank_unequal_dimensions(vec, dim, expected_result):
+    """Schmidt rank is correct when the two subsystems have different dimensions."""
+    assert schmidt_rank(vec, dim) == expected_result
+
+
+@pytest.mark.parametrize(
+    "vec, dims, expected_result",
+    [
+        # Product state on 2x3: operator Schmidt rank 1.
+        (np.kron(np.array([[1], [2]]) / np.sqrt(5), np.array([[1], [2], [3]]) / np.sqrt(14)), [2, 3], 1),
+        # Pure state with Schmidt rank 2 on 2x3: operator Schmidt rank 2**2 = 4.
+        ((np.kron([[1], [0]], [[1], [0], [0]]) + np.kron([[0], [1]], [[0], [1], [0]])) / np.sqrt(2), [2, 3], 4),
+    ],
+)
+def test_operator_schmidt_rank_unequal_dimensions(vec, dims, expected_result):
+    """Operator Schmidt rank is correct when the two subsystems have different dimensions."""
+    rho = vec @ vec.conj().T
+    assert schmidt_rank(rho, dims) == expected_result