Qudit information mps

examples/qudit_mps.py

@@ -0,0 +1,127 @@
+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)
+
+    np.testing.assert_allclose(mps.wavefunction(), qc.wavefunction())
+    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

@@ -14,6 +14,7 @@
 
 sys.path.insert(0, modulepath)
 import tensorcircuit as tc
+import tensorcircuit.quantum as qu
 
 
 @pytest.mark.parametrize("backend", [lf("npb"), lf("cpb")])
@@ -445,3 +446,75 @@ 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")])
+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` (via Rényi entropy);
+    - linear entropy :math:`S_L = 1 - \operatorname{Tr}(\rho_A^2) = 1 - 1/d`.
+    """
+    d = 3
+    c = tc.QuditCircuit(2, d)
+    c.h(0)
+    c.csum(0, 1)
+
+    # Reduced density matrix of subsystem A (trace out qudit-1)
+    rho_A = qu.reduced_density_matrix(c.state(), cut=[1], dim=d)
+
+    # Spectrum check
+    evals = tc.backend.eigh(rho_A)[0]
+    np.testing.assert_allclose(evals, np.ones(d) / d, rtol=1e-6, atol=1e-7)
+
+    # Purity from Rényi entropy of order 2
+    S2 = qu.renyi_entropy(rho_A, k=2)
+    purity = float(np.exp(-S2))
+    np.testing.assert_allclose(purity, 1.0 / d, rtol=1e-6, atol=1e-7)
+
+    # Linear entropy check
+    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")])
+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 from built-in functions) 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: Product state
+    c1 = tc.QuditCircuit(2, d)
+    c1.ry(0, theta=0.37, j=0, k=1)
+    c1.ry(1, theta=-0.59, j=0, k=1)
+    I_AB_1 = qu.mutual_information(c1.state(), cut=[1], dim=d)
+    np.testing.assert_allclose(I_AB_1, 0.0, atol=1e-7)
+
+    # Case 2: Maximally entangled state
+    c2 = tc.QuditCircuit(2, d)
+    c2.h(0)
+    c2.csum(0, 1)
+    I_AB_2 = qu.mutual_information(c2.state(), cut=[1], dim=d)
+    np.testing.assert_allclose(I_AB_2, 2.0 * np.log(d), rtol=1e-6, atol=1e-7)
+
+    # Optional: confirm single-subsystem entropy equals log(d)
+    rho_A = qu.reduced_density_matrix(c2.state(), cut=[1], dim=d)
+    SA = qu.entropy(rho_A)
+    np.testing.assert_allclose(SA, np.log(d), rtol=1e-6, atol=1e-7)