Qudit information mps

examples/qudit_mps.py

@@ -0,0 +1,143 @@
+r"""
+This script performs a sanity check for qudit (d>2) circuits using the MPSCircuit class
+by explicitly applying unitary gates.
+
+It constructs a small qutrit (d=3) circuit and applies gates exclusively through the
+unitary method of MPSCircuit, using built-in qudit gate matrices from tensorcircuit.quditgates.
+
+The same circuit is built using QuditCircuit to provide a dense reference statevector.
+The resulting statevectors from both approaches are compared to verify correctness
+of MPS evolution for qudits.
+
+The tested 3-qutrit circuit consists of the following gates:
+1) Hadamard (Hd) on wire 0
+2) Generalized Pauli-X (Xd) on wire 1
+3) Controlled-sum (CSUM) on wires (0, 1)
+4) Controlled-phase (CPHASE) on wires (1, 2)
+5) RZ rotation on wire 2
+
+Numerical closeness between the MPS and dense states is asserted.
+"""
+
+# import os
+# import sys
+#
+# base_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
+# if base_dir not in sys.path:
+#     sys.path.insert(0, base_dir)
+
+import numpy as np
+import tensorcircuit as tc
+
+tc.set_backend("jax")
+tc.set_dtype("complex128")
+
+from tensorcircuit.quditgates import (
+    x_matrix_func,
+    h_matrix_func,
+    csum_matrix_func,
+    cphase_matrix_func,
+    rz_matrix_func,
+)
+
+
+def build_qutrit_circuit_mps(n=3, d=3, theta=np.pi / 7):
+    """
+    Build an MPSCircuit of size n with local dimension d=3 and apply qudit gates via the unitary method.
+
+    The circuit applies the following gates in order:
+    - Hadamard (Hd) on qudit 0
+    - Generalized Pauli-X (Xd) on qudit 1
+    - Controlled-sum (CSUM) on qudits (0, 1)
+    - Controlled-phase (CPHASE) on qudits (1, 2)
+    - RZ rotation with angle theta on qudit 2 (acting nontrivially on a chosen level pair)
+
+    :param int n: Number of qudits in the circuit.
+    :param int d: Local dimension of each qudit (default is 3 for qutrits).
+    :param float theta: Rotation angle for the RZ gate.
+    :return: The constructed MPSCircuit with applied gates.
+    :rtype: tc.MPSCircuit
+    """
+    omega = np.exp(2j * np.pi / d)
+
+    mps = tc.MPSCircuit(n, dim=d)
+
+    Hd = h_matrix_func(d, omega)
+    mps.unitary(0, unitary=Hd, name="H_d", dim=d)  # <-- pass dim=d
+
+    Xd = x_matrix_func(d)
+    mps.unitary(1, unitary=Xd, name="X_d", dim=d)  # <-- pass dim=d
+
+    CSUM = csum_matrix_func(d)  # (d*d, d*d)
+    mps.unitary(0, 1, unitary=CSUM, name="CSUM", dim=d)  # <-- pass dim=d
+
+    CPHASE = cphase_matrix_func(d, omega=omega)  # (d*d, d*d)
+    mps.unitary(1, 2, unitary=CPHASE, name="CPHASE", dim=d)  # <-- pass dim=d
+
+    RZ2 = rz_matrix_func(d, theta=theta, j=1)  # (d, d)
+    mps.unitary(2, unitary=RZ2, name="RZ_lvl1", dim=d)  # <-- pass dim=d
+
+    return mps
+
+
+def build_qutrit_circuit_dense(n=3, d=3, theta=np.pi / 7):
+    """
+    Build the same qudit circuit using QuditCircuit with high-level named qudit gates.
+
+    This circuit serves as a dense reference for cross-checking against the MPSCircuit.
+
+    :param int n: Number of qudits in the circuit.
+    :param int d: Local dimension of each qudit (default is 3 for qutrits).
+    :param float theta: Rotation angle for the RZ gate.
+    :return: The constructed QuditCircuit with applied gates.
+    :rtype: tc.QuditCircuit
+    """
+    qc = tc.QuditCircuit(n, dim=d)
+
+    qc.h(0)  # H_d
+    qc.x(1)  # X_d
+
+    qc.csum(0, 1)
+    qc.cphase(1, 2)
+
+    qc.rz(2, theta=tc.num_to_tensor(theta), j=1)
+
+    return qc
+
+
+def main():
+    """
+    Construct both MPSCircuit and QuditCircuit for a 3-qutrit circuit and verify their statevectors match.
+
+    This function checks that the maximum absolute difference between the MPS and dense statevectors
+    is below a small numerical tolerance, ensuring correctness of MPS evolution for qudits.
+    """
+    n, d = 3, 3
+    theta = np.pi / 7
+
+    # Build circuits
+    mps = build_qutrit_circuit_mps(n=n, d=d, theta=theta)
+    qc = build_qutrit_circuit_dense(n=n, d=d, theta=theta)
+
+    # Obtain statevectors (shape [-1])
+    psi_mps = tc.backend.numpy(mps.wavefunction(form="default"))
+    psi_dense = tc.backend.numpy(qc.wavefunction(form="default"))
+
+    # Normalize both just in case of tiny numerical drift
+    def _norm(v):
+        return v / np.linalg.norm(v)
+
+    psi_mps = _norm(psi_mps)
+    psi_dense = _norm(psi_dense)
+
+    # \ell_\infty comparison: :math:`\max_i |(\psi_{\mathrm{MPS}} - \psi_{\mathrm{dense}})_i|`
+    inf_err = np.max(np.abs(psi_mps - psi_dense))
+    print(r"Max $|\Delta|$ between MPS and dense:", inf_err)
+
+    tol = 1e-10
+    assert inf_err < tol, f"States differ too much: {inf_err} >= {tol}"
+    print("OK: MPSCircuit matches QuditCircuit for d=3 with unitary-applied gates.")
+
+
+if __name__ == "__main__":
+    main()

examples/vqe_qudit_example.py

@@ -16,7 +16,8 @@
 The code defaults to a 2-qutrit (d=3) problem but can be changed via CLI flags.
 """
 
-# import os, sys
+# import os
+# import sys
 #
 # base_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
 # if base_dir not in sys.path:

tests/test_quditcircuit.py

@@ -445,3 +445,106 @@ def test_quditcircuit_amplitude_before_wrapper():
     nodes = c.amplitude_before("00")
     assert isinstance(nodes, list)
     assert len(nodes) == 5  # one gate (X on qudit 0) -> single node in the traced path
+
+
+@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb"), lf("cpb")])
+def test_qudit_entanglement_measures_maximally_entangled(backend):
+    r"""
+    Prepare the two-qudit maximally entangled state
+    :math:`\lvert \Phi_d \rangle = \frac{1}{\sqrt{d}} \sum_{j=0}^{d-1} \lvert j \rangle \otimes \lvert j \rangle`
+    for :math:`d>2` using a generalized Hadamard :math:`H` on qudit-0 followed by
+    :math:`\mathrm{CSUM}(0\!\to\!1)`. For this state,
+    :math:`\rho_A = \operatorname{Tr}_B(\lvert \Phi_d \rangle\langle \Phi_d \rvert) = \mathbb{I}_d/d`.
+    We check:
+    - eigenvalues of :math:`\rho_A` are all :math:`1/d`;
+    - purity :math:`\operatorname{Tr}(\rho_A^2) = 1/d`;
+    - linear entropy :math:`S_L = 1 - \operatorname{Tr}(\rho_A^2) = 1 - 1/d`.
+    """
+    d = 3  # any d>2 is fine; use qutrit here
+    c = tc.QuditCircuit(2, d)
+    c.h(0)
+    c.csum(0, 1)
+
+    # |\psi> as a d*d amplitude matrix: \psi[j_A, j_B]
+    psi = tc.backend.numpy(c.state()).reshape(d, d)
+    rho_A = psi @ psi.conj().T  # partial trace over B
+
+    evals = np.linalg.eigvalsh(rho_A)
+    np.testing.assert_allclose(evals, np.ones(d) / d, rtol=1e-6, atol=1e-7)
+
+    purity = np.real_if_close(np.trace(rho_A @ rho_A))
+    np.testing.assert_allclose(purity, 1.0 / d, rtol=1e-6, atol=1e-7)
+
+    linear_entropy = 1.0 - purity
+    np.testing.assert_allclose(linear_entropy, 1.0 - 1.0 / d, rtol=1e-6, atol=1e-7)
+
+
+@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb"), lf("cpb")])
+def test_qudit_mutual_information_product_vs_entangled(backend):
+    r"""
+    Compare quantum mutual information :math:`I(A\!:\!B) = S(\rho_A)+S(\rho_B)-S(\rho_{AB})`
+    (with von Neumann entropy :math:`S(\rho)=-\operatorname{Tr}[\rho \ln \rho]`, natural log)
+    for two two-qudit states (local dimension :math:`d>2`):
+
+    1) **Product state** prepared *with gates* by applying single-qudit rotations
+       :math:`\mathrm{RY}(\theta)` in the two-level subspace :math:`(j,k)=(0,1)` on **each** qudit.
+       This yields a separable pure state, hence :math:`I(A\!:\!B)=0`.
+
+    2) **Maximally entangled state** :math:`\lvert \Phi_d \rangle` prepared by :math:`H` on qudit-0
+       then :math:`\mathrm{CSUM}(0\!\to\!1)`. For this pure bipartite state,
+       :math:`\rho_A=\rho_B=\mathbb{I}_d/d`, so :math:`S(\rho_A)=S(\rho_B)=\ln d`,
+       :math:`S(\rho_{AB})=0`, and :math:`I(A\!:\!B)=2\ln d`.
+    """
+    d = 3
+
+    # --- Case 1: explicit product state via local rotations (uses native ry) ---
+    c1 = tc.QuditCircuit(2, d)
+    # rotate each qudit in subspace (0,1) by different angles to ensure it's not trivial
+    c1.ry(0, theta=0.37, j=0, k=1)
+    c1.ry(1, theta=-0.59, j=0, k=1)
+    psi1 = tc.backend.numpy(c1.state()).reshape(d, d)
+
+    # --- Case 2: maximally entangled |\phi_d> via H + CSUM ---
+    c2 = tc.QuditCircuit(2, d)
+    c2.h(0)
+    c2.csum(0, 1)
+    psi2 = tc.backend.numpy(c2.state()).reshape(d, d)
+
+    def reduced_density(psi, trace_out="B"):
+        # \rho_A = \psi \psi^\dagger ; \rho_B = \psi^T \psi^*
+        return psi @ psi.conj().T if trace_out == "B" else psi.T @ psi.conj()
+
+    def von_neumann_entropy(rho):
+        # Hermitize for stability, drop ~0 eigenvalues before log.
+        rh = (rho + rho.conj().T) / 2
+        vals = np.linalg.eigvalsh(rh)
+        vals = np.clip(np.real_if_close(vals), 0.0, 1.0)
+        nz = vals > 1e-12
+        return float(-np.sum(vals[nz] * np.log(vals[nz])))
+
+    # product state mutual information
+    rhoA1, rhoB1 = reduced_density(psi1, "B"), reduced_density(psi1, "A")
+    rhoAB1 = np.outer(psi1.reshape(-1), psi1.conj().reshape(-1)).reshape(d * d, d * d)
+    I_AB_1 = (
+        von_neumann_entropy(rhoA1)
+        + von_neumann_entropy(rhoB1)
+        - von_neumann_entropy(rhoAB1)
+    )
+    np.testing.assert_allclose(I_AB_1, 0.0, atol=1e-7)
+
+    # maximally entangled mutual information
+    rhoA2, rhoB2 = reduced_density(psi2, "B"), reduced_density(psi2, "A")
+    rhoAB2 = np.outer(psi2.reshape(-1), psi2.conj().reshape(-1)).reshape(d * d, d * d)
+    I_AB_2 = (
+        von_neumann_entropy(rhoA2)
+        + von_neumann_entropy(rhoB2)
+        - von_neumann_entropy(rhoAB2)
+    )
+    np.testing.assert_allclose(von_neumann_entropy(rhoAB2), 0.0, atol=1e-7)
+    np.testing.assert_allclose(
+        von_neumann_entropy(rhoA2), np.log(d), rtol=1e-6, atol=1e-7
+    )
+    np.testing.assert_allclose(
+        von_neumann_entropy(rhoB2), np.log(d), rtol=1e-6, atol=1e-7
+    )
+    np.testing.assert_allclose(I_AB_2, 2.0 * np.log(d), rtol=1e-6, atol=1e-7)