implemented commutator function + tests

Author: positr0nium

src/qrisp/operators/qubit/qubit_hamiltonian.py

@@ -425,6 +425,23 @@ def subs(self, subs_dict):
     #
     # Miscellaneous
     #
+    
+    def commutator(self, other):
+        
+        res = 0
+        
+        for term_self, coeff_self in self.terms_dict.items():
+            for term_other, coeff_other in other.terms_dict.items():
+                res += coeff_self*coeff_other*term_self.commutator(term_other)
+                
+        min_coeff_self = min([abs(coeff) for coeff in self.terms_dict.values()])
+        min_coeff_other = min([abs(coeff) for coeff in other.terms_dict.values()])
+        
+        res.apply_threshold(min_coeff_self*min_coeff_other/2)
+        
+        return res
+        
+    
 
     def apply_threshold(self,threshold):
         """
@@ -501,7 +518,7 @@ def recursive_TP(keys,term_dict):
                 M = sp.csr_matrix((1,1))
                 M[0,0] = res
                 return M
-        elif factor_amount < max(participating_indices):
+        elif participating_indices and factor_amount < max(participating_indices):
             raise Exception("Tried to compute Hermitian matrix with factor_amount variable lower than the largest factor index")
 
         keys = list(range(factor_amount))
@@ -512,10 +529,9 @@ def recursive_TP(keys,term_dict):
         M = sp.csr_matrix((2**dim, 2**dim))
         for k in range(m):
             M += complex(coeffs[k])*recursive_TP(keys.copy(),term_dicts[k])
-
-        res = ((M + M.transpose().conjugate())/2)
-        res.sum_duplicates()
-        return res
+        # res = ((M + M.transpose().conjugate())/2)
+        # res.sum_duplicates()
+        return M
 
     def ground_state_energy(self):
         """

src/qrisp/operators/qubit/qubit_term.py

@@ -55,25 +55,6 @@ def is_identity(self):
     #
     # Simulation
     #
-    
-    # Assume that the operator is diagonal after change of basis
-    # Implements exp(i*coeff*\prod_j Z_j) where the product goes over all indices j in self.factor_dict
-    # @lifted
-    def simulate(self, coeff, qarg):
-
-        def parity(qarg, indices):
-            n = len(indices)
-            for i in range(n-1):
-                cx(qarg[indices[i]],qarg[indices[i+1]])
-
-        if not self.is_identity():
-            indices = list(self.factor_dict.keys())
-            with conjugate(parity)(qarg, indices):
-                rz(-2*coeff,qarg[indices[-1]])
-        else:
-            gphase(coeff,qarg[0])
-
-    # @lifted
     def simulate(self, coeff, qv, do_change_of_basis = True):
 
         from qrisp import h, cx, rz, mcp, conjugate, control, QuantumBool, mcx, x, p, s, QuantumEnvironment, gphase
@@ -393,6 +374,8 @@ def __mul__(self, other):
         keys = set(a.keys()) | set(b.keys())
         for key in keys:
             factor, coeff = PAULI_TABLE[a.get(key,"I"),b.get(key,"I")]
+            if coeff == 0:
+                return QubitTerm({}), 0
             if factor != "I":
                 result_factor_dict[key]=factor
                 result_coeff *= coeff
@@ -425,8 +408,18 @@ def subs(self, subs_dict):
         return QubitTerm(result_factor_dict), result_coeff
 
     #
-    # Commutativity checks
+    # Commutativity
     #
+    def commutator(self, other):
+        from qrisp.operators.qubit import QubitHamiltonian
+        
+        term_0, coeff_0 = self*other
+        term_1, coeff_1 = other*self
+        
+        temp =  QubitHamiltonian({term_0 : coeff_0}) - QubitHamiltonian({term_1 : coeff_1})
+        temp.apply_threshold(0.5)
+        
+        return temp
 
     def commute(self, other):
         """
@@ -461,7 +454,7 @@ def commute_qw(self, other):
         keys.update(set(b.keys()))
 
         for key in keys:
-            if a.get(key,"I")!="I" and b.get(key,"I")!="I" and a.get(key,"I")!=b.get(key,"I"):
+            if not PAULI_TABLE[a.get(key,"I"), b.get(key,"I")] == PAULI_TABLE[b.get(key,"I"), a.get(key,"I")]:
                 return False
         return True
 

tests/hamiltonian_tests/test_qubit_hamiltonian_commutator.py

@@ -0,0 +1,51 @@
+"""
+\********************************************************************************
+* Copyright (c) 2024 the Qrisp authors
+*
+* This program and the accompanying materials are made available under the
+* terms of the Eclipse Public License 2.0 which is available at
+* http://www.eclipse.org/legal/epl-2.0.
+*
+* This Source Code may also be made available under the following Secondary
+* Licenses when the conditions for such availability set forth in the Eclipse
+* Public License, v. 2.0 are satisfied: GNU General Public License, version 2
+* with the GNU Classpath Exception which is
+* available at https://www.gnu.org/software/classpath/license.html.
+*
+* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0
+********************************************************************************/
+"""
+
+from qrisp import QuantumVariable, x, QPE
+from qrisp.operators import X, Y, Z, A, C, P0, P1
+import numpy as np
+from numpy.linalg import norm
+
+def test_qubit_hamiltonian_commutator():
+
+    def verify_commutator(H1, H2):
+        commutator_hamiltonian = H1.commutator(H2)
+        
+        commutator_matrix = commutator_hamiltonian.to_sparse_matrix(factor_amount = 3).todense()
+        
+        H1_matrix = H1.to_sparse_matrix(factor_amount = 3).todense()
+        H2_matrix = H2.to_sparse_matrix(factor_amount = 3).todense()
+        
+        assert norm(commutator_matrix - (H1_matrix @ H2_matrix - H2_matrix @ H1_matrix)) < 1E-5
+
+    operator_list = [lambda x : 1, X, Y, Z, A, C, P0, P1]
+
+    counter = 0
+    operator_list = [lambda x : 1, X, Y, Z, A, C, P0, P1]
+    for O0 in operator_list: 
+        for O1 in operator_list:
+            for O2 in operator_list:
+                for O3 in operator_list:
+                    H1 = O0(0)*O1(1)
+                    H2 = O2(0)*O3(1)
+                    
+                    if 1 in [H1, H2]:
+                        continue
+                    verify_commutator(H1, H2)
+                    counter += 1
+

tests/hamiltonian_tests/test_qubit_hamiltonian_simulation.py

@@ -52,7 +52,7 @@ def verify_trotterization(H):
 
         # Compute hermitean matrix
         H_matrix = H.to_sparse_matrix().todense()
-        
+        H_matrix = (H_matrix + H_matrix.transpose().conjugate())/2
         # Compute unitary matrix
         U_matrix = expm(1j*H_matrix)