From fa54ebbb45ef4efea5d02195dc8bfffccb713786 Mon Sep 17 00:00:00 2001 From: Weiguo Ma Date: Wed, 10 Sep 2025 15:06:38 +0800 Subject: [PATCH 1/7] add new tests for qudit(d>2) in quantum information and entanglement. --- tests/test_quditcircuit.py | 103 +++++++++++++++++++++++++++++++++++++ 1 file changed, 103 insertions(+) diff --git a/tests/test_quditcircuit.py b/tests/test_quditcircuit.py index 492071b5..5f929e1f 100644 --- a/tests/test_quditcircuit.py +++ b/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) From dffeeed6186cfce416eb4924d6b569c25aab1cba Mon Sep 17 00:00:00 2001 From: Weiguo Ma Date: Wed, 10 Sep 2025 16:38:04 +0800 Subject: [PATCH 2/7] format imports --- examples/vqe_qudit_example.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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: From 6789112b41c7d5666cf089b3057d1627f015fc7d Mon Sep 17 00:00:00 2001 From: Weiguo Ma Date: Wed, 10 Sep 2025 16:44:23 +0800 Subject: [PATCH 3/7] an example script for MPSCircuit using dim>2 and c.uitary for gates application to check that MPSCircuit is workingfor qudit case --- examples/qudit_mps.py | 143 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 143 insertions(+) create mode 100644 examples/qudit_mps.py diff --git a/examples/qudit_mps.py b/examples/qudit_mps.py new file mode 100644 index 00000000..dcd24d28 --- /dev/null +++ b/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() From 0c61a9dbce63969a87e02d6717d809a1652f2971 Mon Sep 17 00:00:00 2001 From: Weiguo Ma Date: Wed, 10 Sep 2025 18:03:53 +0800 Subject: [PATCH 4/7] remove cupy tests --- tests/test_quditcircuit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_quditcircuit.py b/tests/test_quditcircuit.py index 5f929e1f..da7e1007 100644 --- a/tests/test_quditcircuit.py +++ b/tests/test_quditcircuit.py @@ -447,7 +447,7 @@ def test_quditcircuit_amplitude_before_wrapper(): 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")]) +@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 @@ -479,7 +479,7 @@ def test_qudit_entanglement_measures_maximally_entangled(backend): 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")]) +@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})` From cee34e600d7950bf51bf0f9fe4f034778b057a81 Mon Sep 17 00:00:00 2001 From: Weiguo Ma Date: Wed, 10 Sep 2025 18:08:43 +0800 Subject: [PATCH 5/7] optimize code --- examples/qudit_mps.py | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/examples/qudit_mps.py b/examples/qudit_mps.py index dcd24d28..6dc63c68 100644 --- a/examples/qudit_mps.py +++ b/examples/qudit_mps.py @@ -119,23 +119,7 @@ def main(): 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}" + np.testing.assert_allclose(mps.wavefunction(), qc.wavefunction()) print("OK: MPSCircuit matches QuditCircuit for d=3 with unitary-applied gates.") From c920d69d666ae55d0c348ca4df69a9236bf3c920 Mon Sep 17 00:00:00 2001 From: Weiguo Ma Date: Wed, 10 Sep 2025 18:16:23 +0800 Subject: [PATCH 6/7] optimized tests in quditcircuit --- tests/test_quditcircuit.py | 74 ++++++++++++-------------------------- 1 file changed, 22 insertions(+), 52 deletions(-) diff --git a/tests/test_quditcircuit.py b/tests/test_quditcircuit.py index da7e1007..80fc34da 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")]) @@ -457,24 +458,27 @@ def test_qudit_entanglement_measures_maximally_entangled(backend): :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`; + - 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 # any d>2 is fine; use qutrit here + d = 3 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 + # Reduced density matrix of subsystem A (trace out qudit-1) + rho_A = qu.reduced_density_matrix(c.state(), cut=[1], dim=d) - evals = np.linalg.eigvalsh(rho_A) + # 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 = np.real_if_close(np.trace(rho_A @ rho_A)) + # 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) @@ -483,8 +487,7 @@ def test_qudit_entanglement_measures_maximally_entangled(backend): 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`): + (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. @@ -497,54 +500,21 @@ def test_qudit_mutual_information_product_vs_entangled(backend): """ d = 3 - # --- Case 1: explicit product state via local rotations (uses native ry) --- + # Case 1: Product state 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) + 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 |\phi_d> via H + CSUM --- + # Case 2: Maximally entangled state 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 - ) + 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) \ No newline at end of file From 4652d83898a5e151a64c184275a5232a0a0f3448 Mon Sep 17 00:00:00 2001 From: Weiguo Ma Date: Wed, 10 Sep 2025 18:20:25 +0800 Subject: [PATCH 7/7] black . --- tests/test_quditcircuit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_quditcircuit.py b/tests/test_quditcircuit.py index 80fc34da..e2118986 100644 --- a/tests/test_quditcircuit.py +++ b/tests/test_quditcircuit.py @@ -517,4 +517,4 @@ def test_qudit_mutual_information_product_vs_entangled(backend): # 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) \ No newline at end of file + np.testing.assert_allclose(SA, np.log(d), rtol=1e-6, atol=1e-7)