QubitOperator term simplification

src/openfermion/hamiltonians/_hubbard_test.py

@@ -461,6 +461,7 @@ def test_bose_hubbard_2x3():
 2.0 [0^ 0 0^ 0] +
 0.3 [0^ 0 1^ 1] +
 0.3 [0^ 0 2^ 2] +
+0.3 [0^ 0 4^ 4] +
 -1.0 [0^ 1] +
 -1.0 [0^ 2] +
 -1.0 [0^ 4] +
@@ -469,6 +470,7 @@ def test_bose_hubbard_2x3():
 -2.5 [1^ 1] +
 2.0 [1^ 1 1^ 1] +
 0.3 [1^ 1 3^ 3] +
+0.3 [1^ 1 5^ 5] +
 -1.0 [1^ 3] +
 -1.0 [1^ 5] +
 -1.0 [2 3^] +
@@ -486,12 +488,10 @@ def test_bose_hubbard_2x3():
 -1.0 [3^ 5] +
 -1.0 [4 5^] +
 -2.5 [4^ 4] +
-0.3 [4^ 4 0^ 0] +
 2.0 [4^ 4 4^ 4] +
 0.3 [4^ 4 5^ 5] +
 -1.0 [4^ 5] +
 -2.5 [5^ 5] +
-0.3 [5^ 5 1^ 1] +
 2.0 [5^ 5 5^ 5]
 """.strip()
 
@@ -508,6 +508,7 @@ def test_bose_hubbard_3x2():
 -2.5 [0^ 0] +
 2.0 [0^ 0 0^ 0] +
 0.3 [0^ 0 1^ 1] +
+0.3 [0^ 0 2^ 2] +
 0.3 [0^ 0 3^ 3] +
 -1.0 [0^ 1] +
 -1.0 [0^ 2] +
@@ -522,7 +523,6 @@ def test_bose_hubbard_3x2():
 -1.0 [1^ 4] +
 -1.0 [2 5^] +
 -2.5 [2^ 2] +
-0.3 [2^ 2 0^ 0] +
 2.0 [2^ 2 2^ 2] +
 0.3 [2^ 2 5^ 5] +
 -1.0 [2^ 5] +
@@ -531,6 +531,7 @@ def test_bose_hubbard_3x2():
 -2.5 [3^ 3] +
 2.0 [3^ 3 3^ 3] +
 0.3 [3^ 3 4^ 4] +
+0.3 [3^ 3 5^ 5] +
 -1.0 [3^ 4] +
 -1.0 [3^ 5] +
 -1.0 [4 5^] +
@@ -539,7 +540,6 @@ def test_bose_hubbard_3x2():
 0.3 [4^ 4 5^ 5] +
 -1.0 [4^ 5] +
 -2.5 [5^ 5] +
-0.3 [5^ 5 3^ 3] +
 2.0 [5^ 5 5^ 5]
 """.strip()
 

src/openfermion/ops/_qubit_operator.py

@@ -108,82 +108,6 @@ def different_indices_commute(self):
         """Whether factors acting on different indices commute."""
         return True
 
-    def __imul__(self, multiplier):
-        """
-        Override in-place multiply of SymbolicOperator
-
-        Args:
-          multiplier(complex float, or QubitOperator): multiplier
-        """
-        # Handle scalars.
-        if isinstance(multiplier, (int, float, complex)):
-            for term in self.terms:
-                self.terms[term] *= multiplier
-            return self
-
-        # Handle QubitOperator.
-        if not isinstance(multiplier, QubitOperator):
-            raise TypeError('Cannot in-place multiply term of invalid type to '
-                            'QubitTerm.')
-
-        result_terms = dict()
-        for left_term in self.terms:
-            for right_term in multiplier.terms:
-                new_coefficient = (self.terms[left_term] *
-                                   multiplier.terms[right_term])
-
-                # Loop through local operators and create new sorted list
-                # of representing the product local operator:
-                product_operators = []
-                left_operator_index = 0
-                right_operator_index = 0
-                n_operators_left = len(left_term)
-                n_operators_right = len(right_term)
-                while (left_operator_index < n_operators_left and
-                       right_operator_index < n_operators_right):
-                    (left_qubit, left_loc_op) = (
-                        left_term[left_operator_index])
-                    (right_qubit, right_loc_op) = (
-                        right_term[right_operator_index])
-
-                    # Multiply local operators acting on the same qubit
-                    if left_qubit == right_qubit:
-                        left_operator_index += 1
-                        right_operator_index += 1
-                        (scalar, loc_op) = _PAULI_OPERATOR_PRODUCTS[
-                            (left_loc_op, right_loc_op)]
-
-                        # Add new term.
-                        if loc_op != 'I':
-                            product_operators += [(left_qubit, loc_op)]
-                            new_coefficient *= scalar
-                        # Note if loc_op == 'I', then scalar == 1.0
-
-                    # If left_qubit > right_qubit, add right_loc_op; else,
-                    # add left_loc_op.
-                    elif left_qubit > right_qubit:
-                        product_operators += [(right_qubit, right_loc_op)]
-                        right_operator_index += 1
-                    else:
-                        product_operators += [(left_qubit, left_loc_op)]
-                        left_operator_index += 1
-
-                # Finish the remaining operators:
-                if left_operator_index == n_operators_left:
-                    product_operators += right_term[
-                        right_operator_index::]
-                elif right_operator_index == n_operators_right:
-                    product_operators += left_term[left_operator_index::]
-
-                # Add to result dict
-                tmp_key = tuple(product_operators)
-                if tmp_key in result_terms:
-                    result_terms[tmp_key] += new_coefficient
-                else:
-                    result_terms[tmp_key] = new_coefficient
-        self.terms = result_terms
-        return self
-
     def renormalize(self):
         """Fix the trace norm of an operator to 1"""
         norm = self.induced_norm(2)
@@ -228,3 +152,35 @@ def get_operator_groups(self, num_groups):
         for i in range(num_groups):
             yield QubitOperator.accumulate(itertools.islice(
                 operators, len(range(i, len(self.terms), num_groups))))
+
+    def _simplify(self, term, coefficient=1.0):
+        """Simplify a term using commutator and anti-commutator relations."""
+        if not term:
+            return coefficient, term
+
+        term = sorted(term, key=lambda factor: factor[0])
+
+        new_term = []
+        left_factor = term[0]
+        for right_factor in term[1:]:
+            left_index, left_action = left_factor
+            right_index, right_action = right_factor
+
+            # Still on the same qubit, keep simplifying.
+            if left_index == right_index:
+                new_coefficient, new_action = _PAULI_OPERATOR_PRODUCTS[
+                        left_action, right_action]
+                left_factor = (left_index, new_action)
+                coefficient *= new_coefficient
+
+            # Reached different qubit, save result and re-initialize.
+            else:
+                if left_action != 'I':
+                    new_term.append(left_factor)
+                left_factor = right_factor
+
+        # Save result of final iteration.
+        if left_factor[1] != 'I':
+            new_term.append(left_factor)
+
+        return coefficient, tuple(new_term)

src/openfermion/ops/_qubit_operator_test.py

@@ -39,6 +39,26 @@ def test_pauli_operator_product():
     assert _PAULI_OPERATOR_PRODUCTS == correct
 
 
+def test_init_simplify():
+    assert QubitOperator("X0 X0") == QubitOperator.identity()
+    assert QubitOperator("X1 X1") == QubitOperator.identity()
+    assert QubitOperator("Y0 Y0") == QubitOperator.identity()
+    assert QubitOperator("Z0 Z0") == QubitOperator.identity()
+    assert QubitOperator("X0 Y0") == QubitOperator("Z0", coefficient=1j)
+    assert QubitOperator("Y0 X0") == QubitOperator("Z0", coefficient=-1j)
+    assert QubitOperator("X0 Y0 Z0") == 1j * QubitOperator.identity()
+    assert QubitOperator("Y0 Z0 X0") == 1j * QubitOperator.identity()
+    assert QubitOperator("Z0 Y0 X0") == -1j * QubitOperator.identity()
+    assert QubitOperator("X1 Y0 X1") == QubitOperator("Y0")
+    assert QubitOperator("Y1 Y1 Y1") == QubitOperator("Y1")
+    assert QubitOperator("Y2 Y2 Y2 Y2") == QubitOperator.identity()
+    assert QubitOperator("Y3 Y3 Y3 Y3 Y3") == QubitOperator("Y3")
+    assert QubitOperator("Y4 Y4 Y4 Y4 Y4 Y4") == QubitOperator.identity()
+    assert QubitOperator("X0 Y1 Y0 X1") == QubitOperator("Z0 Z1")
+    assert QubitOperator("X0 Y1 Z3 X2 Z3 Y0") == QubitOperator("Z0 Y1 X2",
+            coefficient=1j)
+
+
 def test_imul_inplace():
     qubit_op = QubitOperator("X1")
     prev_id = id(qubit_op)

src/openfermion/ops/_symbolic_operator.py

@@ -129,6 +129,9 @@ def __init__(self, term=None, coefficient=1.):
         else:
             raise ValueError('term specified incorrectly.')
 
+        # Simplify the term
+        coefficient, term = self._simplify(term, coefficient=coefficient)
+
         # Add the term to the dictionary
         self.terms[term] = coefficient
 
@@ -165,8 +168,9 @@ def _long_string_init(self, long_string, coefficient):
                             'Invalid coefficient {}.'.format(coef_string))
             coef *= coefficient
 
-            # Parse the term and add it to the dict
+            # Parse the term, simpify it and add to the dict
             term = self._parse_string(match[1])
+            coef, term = self._simplify(term, coefficient=coef)
             if term not in self.terms:
                 self.terms[term] = coef
             else:
@@ -189,6 +193,12 @@ def _validate_factor(self, factor):
                              'The index should be a non-negative '
                              'integer.'.format(factor))
 
+    def _simplify(self, term, coefficient=1.0):
+        """Simplifies a term."""
+        if self.different_indices_commute:
+            term = sorted(term, key=lambda factor: factor[0])
+        return coefficient, tuple(term)
+
     def _parse_sequence(self, term):
         """Parse a term given as a sequence type (i.e., list, tuple, etc.).
 
@@ -207,11 +217,6 @@ def _parse_sequence(self, term):
             for factor in term:
                 self._validate_factor(factor)
 
-            # If factors with different indices commute, sort the factors
-            # by index
-            if self.different_indices_commute:
-                term = sorted(term, key=lambda factor: factor[0])
-
             # Return a tuple
             return tuple(term)
 
@@ -265,12 +270,6 @@ def _parse_string(self, term):
             # Add the factor to the list as a tuple
             processed_term.append((index, action))
 
-        # If factors with different indices commute, sort the factors
-        # by index
-        if self.different_indices_commute:
-            processed_term = sorted(processed_term,
-                                    key=lambda factor: factor[0])
-
         # Return a tuple
         return tuple(processed_term)
 
@@ -343,16 +342,20 @@ def __imul__(self, multiplier):
             result_terms = dict()
             for left_term in self.terms:
                 for right_term in multiplier.terms:
-                    new_coefficient = (self.terms[left_term] *
-                                       multiplier.terms[right_term])
-                    product_operators = left_term + right_term
+                    left_coefficient = self.terms[left_term]
+                    right_coefficient = multiplier.terms[right_term]
+
+                    new_coefficient = left_coefficient * right_coefficient
+                    new_term = left_term + right_term
+
+                    new_coefficient, new_term = self._simplify(
+                            new_term, coefficient=new_coefficient)
 
                     # Update result dict.
-                    product_operators = tuple(product_operators)
-                    if product_operators in result_terms:
-                        result_terms[product_operators] += new_coefficient
+                    if new_term in result_terms:
+                        result_terms[new_term] += new_coefficient
                     else:
-                        result_terms[product_operators] = new_coefficient
+                        result_terms[new_term] = new_coefficient
             self.terms = result_terms
             return self