diff --git a/examples/qudit_mps.py b/examples/qudit_mps.py new file mode 100644 index 00000000..6dc63c68 --- /dev/null +++ b/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() diff --git a/examples/vqe_qudit_example.py b/examples/vqe_qudit_example.py index 46c416a4..ab680e29 100644 --- a/examples/vqe_qudit_example.py +++ b/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: diff --git a/tests/test_quditcircuit.py b/tests/test_quditcircuit.py index 492071b5..e2118986 100644 --- a/tests/test_quditcircuit.py +++ b/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)