Add (state) Schmidt decomposition to quantum_info.linalg_operations

doc/source/api-reference/qibo.rst

@@ -2005,6 +2005,12 @@ Singular value decomposition
 .. autofunction:: qibo.quantum_info.singular_value_decomposition
 
 
+Schmidt decomposition
+"""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.schmidt_decomposition
+
+
 Quantum Networks
 ^^^^^^^^^^^^^^^^
 

src/qibo/quantum_info/linalg_operations.py

@@ -325,3 +325,60 @@ def singular_value_decomposition(matrix, backend=None):
     backend = _check_backend(backend)
 
     return backend.calculate_singular_value_decomposition(matrix)
+
+
+def schmidt_decomposition(
+    state, partition: Union[List[int], Tuple[int, ...]], backend=None
+):
+    """Return the Schmidt decomposition of a :math:`n`-qubit bipartite pure quantum ``state``.
+
+    Given a bipartite pure state :math:`\\ket{\\psi}\\in\\mathcal{H}_{A}\\otimes\\mathcal{H}_{B}`,
+    its Schmidt decomposition is given by
+
+    .. math::
+        \\ket{\\psi} = \\sum_{k = 1}^{\\min\\{a, \\, b\\}} \\, c_{k} \\,
+            \\ket{\\phi_{k}} \\otimes \\ket{\\nu_{k}} \\, ,
+
+    with :math:`a` and :math:`b` being the respective cardinalities of :math:`\\mathcal{H}_{A}`
+    and :math:`\\mathcal{H}_{B}`, and :math:`\\{\\phi_{k}\\}_{k\\in[\\min\\{a, \\, b\\}]}
+    \\subset \\mathcal{H}_{A}` and :math:`\\{\\nu_{k}\\}_{k\\in[\\min\\{a, \\, b\\}]}
+    \\subset \\mathcal{H}_{B}` being orthonormal sets. The coefficients
+    :math:`\\{c_{k}\\}_{k\\in[\\min\\{a, \\, b\\}]}` are real, non-negative, and unique
+    up to re-ordering.
+
+    The decomposition is calculated using :func:`qibo.quantum_info.singular_value_decomposition`,
+    resulting in
+
+    .. math::
+        \\ketbra{\\psi}{\\psi} = U \\, S \\, V^{\\dagger} \\, ,
+
+    where :math:`U` is an :math:`a \\times a` unitary matrix, :math:`V` is an :math:`b \\times b`
+    unitary matrix, and :math:`S` is an :math:`a \\times b` positive semidefinite diagonal matrix
+    that contains the singular values of :math:`\\ketbra{\\psi}{\\psi}`.
+
+    Args:
+        state (ndarray): stevector or density matrix.
+        partition (Union[List[int], Tuple[int, ...]]): indices of qubits in one of the two
+            partitions. The other partition is inferred as the remaining qubits.
+        backend (:class:`qibo.backends.abstract.Backend`, optional): backend
+            to be used in the execution. If ``None``, it uses
+            :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
+
+    Returns:
+        ndarray, ndarray, ndarray: Respectively, the matrices :math:`U`, :math:`S`,
+        and :math:`V^{\\dagger}`.
+    """
+    backend = _check_backend(backend)
+
+    nqubits = math.log2(state.shape[-1])
+    if not nqubits.is_integer():
+        raise_error(ValueError, f"dimensions of ``state`` must be a power of 2.")
+
+    nqubits = int(nqubits)
+    partition_2 = partition.__class__(set(list(range(nqubits))) ^ set(partition))
+
+    tensor = backend.np.reshape(state, [2] * nqubits)
+    tensor = backend.np.transpose(tensor, partition + partition_2)
+    tensor = backend.np.reshape(tensor, (2 ** len(partition), -1))
+
+    return singular_value_decomposition(tensor, backend=backend)

src/qibo/transpiler/unitary_decompositions.py

@@ -2,6 +2,7 @@
 
 from qibo import gates, matrices
 from qibo.config import raise_error
+from qibo.quantum_info.linalg_operations import schmidt_decomposition
 
 magic_basis = np.array(
     [[1, -1j, 0, 0], [0, 0, 1, -1j], [0, 0, -1, -1j], [1, 1j, 0, 0]]
@@ -83,30 +84,7 @@ def calculate_psi(unitary, backend, magic_basis=magic_basis):
     return psi, eigvals
 
 
-def schmidt_decompose(state, backend):
-    """Decomposes a two-qubit product state to its single-qubit parts.
-
-    Args:
-        state (ndarray): product state to be decomposed.
-
-    Returns:
-        (ndarray, ndarray): decomposition
-
-    """
-    # tf.linalg.svd has a different behaviour
-    if backend.__class__.__name__ == "TensorflowBackend":
-        u, d, v = np.linalg.svd(backend.np.reshape(state, (2, 2)))
-    else:
-        u, d, v = backend.np.linalg.svd(backend.np.reshape(state, (2, 2)))
-    if not np.allclose(backend.to_numpy(d), [1, 0]):  # pragma: no cover
-        raise_error(
-            ValueError,
-            f"Unexpected singular values: {d}\nCan only decompose product states.",
-        )
-    return u[:, 0], v[0]
-
-
-def calculate_single_qubit_unitaries(psi, backend):
+def calculate_single_qubit_unitaries(psi, backend=None):
     """Calculates local unitaries that maps a maximally entangled basis to the magic basis.
 
     See Lemma 1 of Appendix A in arXiv:quant-ph/0011050.
@@ -130,8 +108,10 @@ def calculate_single_qubit_unitaries(psi, backend):
     # find e and f by inverting (A3), (A4)
     ef = (psi_bar[0] + 1j * psi_bar[1]) / np.sqrt(2)
     e_f_ = (psi_bar[0] - 1j * psi_bar[1]) / np.sqrt(2)
-    e, f = schmidt_decompose(ef, backend=backend)
-    e_, f_ = schmidt_decompose(e_f_, backend=backend)
+    e, _, f = schmidt_decomposition(ef, [0], backend=backend)
+    e, f = e[:, 0], f[0]
+    e_, _, f_ = schmidt_decomposition(e_f_, [0], backend=backend)
+    e_, f_ = e_[:, 0], f_[0]
     # find exp(1j * delta) using (A5a)
     ef_ = backend.np.kron(e, f_)
     phase = (

tests/test_quantum_info_operations.py

@@ -9,6 +9,7 @@
     matrix_power,
     partial_trace,
     partial_transpose,
+    schmidt_decomposition,
     singular_value_decomposition,
 )
 from qibo.quantum_info.metrics import purity
@@ -263,3 +264,31 @@ def test_singular_value_decomposition(backend):
         S_sorted, coeffs_sorted = S_sorted[0], coeffs_sorted[0]
 
     backend.assert_allclose(S_sorted, coeffs_sorted)
+
+
+def test_schmidt_decomposition(backend):
+    with pytest.raises(ValueError):
+        test = random_statevector(3, backend=backend)
+        test = schmidt_decomposition(test, [0], backend=backend)
+
+    state_A = random_statevector(4, seed=10, backend=backend)
+    state_B = random_statevector(4, seed=11, backend=backend)
+    state = backend.np.kron(state_A, state_B)
+
+    U, S, Vh = schmidt_decomposition(state, [0, 1], backend=backend)
+
+    # recovering original state
+    recovered = np.zeros_like(state.shape, dtype=complex)
+    recovered = backend.cast(recovered, dtype=recovered.dtype)
+    for coeff, u, vh in zip(S, U.T, Vh):
+        if abs(coeff) > 1e-10:
+            recovered = recovered + coeff * backend.np.kron(u, vh)
+
+    backend.assert_allclose(recovered, state)
+
+    # entropy test
+    coeffs = backend.np.abs(S) ** 2
+    entropy = backend.np.where(backend.np.abs(S) < 1e-10, 0.0, backend.np.log(coeffs))
+    entropy = -backend.np.sum(coeffs * entropy)
+
+    backend.assert_allclose(entropy, 0.0, atol=1e-14)