Merge 17bae17 into ff14aaa

src/openfermion/utils/__init__.py

@@ -32,7 +32,8 @@
                               is_hermitian, is_identity,
                               normal_ordered, prune_unused_indices,
                               reorder, up_then_down,
-                              load_operator, save_operator)
+                              load_operator, save_operator,
+                              group_into_tensor_product_basis_sets)
 
 from ._rdm_mapping_functions import (kronecker_delta,
                                      map_two_pdm_to_two_hole_dm,

src/openfermion/utils/_operator_utils.py

@@ -17,6 +17,7 @@
 import os
 import copy
 import marshal
+from numpy.random import RandomState
 
 import numpy
 from scipy.sparse import spmatrix
@@ -846,3 +847,86 @@ def normal_ordered(operator, hbar=1.):
         ordered_operator += order_fn(term, coefficient, **kwargs)
 
     return ordered_operator
+
+
+def _find_compatible_basis(term, bases):
+    for basis in bases:
+        basis_qubits = {op[0] for op in basis}
+        conflicts = ((i, P) for (i, P) in term
+                     if i in basis_qubits and (i, P) not in basis)
+        if any(conflicts):
+            continue
+        return basis
+    return None
+
+
+def group_into_tensor_product_basis_sets(operator, seed=None):
+    """
+    Split an operator (instance of QubitOperator) into `sub-operator`
+    QubitOperators, where each sub-operator has terms that are diagonal
+    in the same tensor product basis.
+
+    Each `sub-operator` can be measured using the same qubit post-rotations
+    in expectation estimation. Grouping into these tensor product basis
+    sets has been found to improve the efficiency of expectation estimation
+    significantly for some Hamiltonians in the context of
+    VQE (see section V(A) in the supplementary material of
+    https://arxiv.org/pdf/1704.05018v2.pdf). The more general problem
+    of grouping operators into commutitative groups is discussed in
+    section IV (B2) of https://arxiv.org/pdf/1509.04279v1.pdf. The
+    original input operator is the union of all output sub-operators,
+    and all sub-operators are disjoint (do not share any terms).
+
+    Args:
+        operator (QubitOperator): the operator that will be split into
+            sub-operators (tensor product basis sets).
+        seed (int): default None. Random seed used to initialize the
+            numpy.RandomState pseudo-random number generator.
+
+    Returns:
+        sub_operators (dict): a dictionary where each key defines a
+            tensor product basis, and each corresponding value is a
+            QubitOperator with terms that are all diagonal in
+            that basis.
+            **key** (tuple of tuples): Each key is a term, which defines
+                a tensor product basis. A term is a product of individual
+                factors; each factor is represented by a tuple of the form
+                (`index`, `action`), and these tuples are collected into a
+                larger tuple which represents the term as the product of
+                its factors. `action` is from the set {'X', 'Y', 'Z'} and
+                `index` is a non-negative integer corresponding to the
+                index of a qubit.
+            **value** (QubitOperator): A QubitOperator with terms that are
+                diagonal in the basis defined by the key it is stored in.
+
+    Raises:
+       TypeError: Operator of invalid type.
+    """
+
+    if not isinstance(operator, QubitOperator):
+        raise TypeError('Can only split QubitOperator into tensor product'
+                        ' basis sets. {} is not supported.'.format(
+                            type(operator).__name__))
+
+    sub_operators = {}
+    r = RandomState(seed)
+    for term, coefficient in operator.terms.items():
+        bases = list(sub_operators.keys())
+        r.shuffle(bases)
+
+        basis = _find_compatible_basis(term, bases)
+
+        if basis is None:
+            sub_operators[term] = QubitOperator(term, coefficient)
+        else:
+            additions = tuple(op for op in term if op not in basis)
+
+            if additions:
+                new_basis = tuple(sorted(basis + additions,
+                                         key=lambda factor: factor[0]))
+                sub_operators[new_basis] = sub_operators.pop(basis)
+                sub_operators[new_basis] += QubitOperator(term, coefficient)
+            else:
+                sub_operators[basis] += QubitOperator(term, coefficient)
+
+    return sub_operators

src/openfermion/utils/_operator_utils_test.py

@@ -843,3 +843,52 @@ def test_quad_triple(self):
     def test_exceptions(self):
         with self.assertRaises(TypeError):
             _ = normal_ordered(1)
+
+
+class GroupTensorProductBasisTest(unittest.TestCase):
+
+    def test_demo_qubit_operator(self):
+        for seed in [None, 0, 10000]:
+            op = QubitOperator('X0 Y1', 2.) + QubitOperator('X1 Y2', 3.j)
+            sub_operators = group_into_tensor_product_basis_sets(op, seed=seed)
+            expected = {((0, 'X'), (1, 'Y')): QubitOperator('X0 Y1', 2.),
+                        ((1, 'X'), (2, 'Y')): QubitOperator('X1 Y2', 3.j)}
+            self.assertEqual(sub_operators, expected)
+
+            op = QubitOperator('X0 Y1', 2.) + QubitOperator('Y1 Y2', 3.j)
+            sub_operators = group_into_tensor_product_basis_sets(op, seed=seed)
+            expected = {((0, 'X'), (1, 'Y'), (2, 'Y')): op}
+            self.assertEqual(sub_operators, expected)
+
+            op = QubitOperator('', 4.) + QubitOperator('X1', 2.j)
+            sub_operators = group_into_tensor_product_basis_sets(op, seed=seed)
+            expected = {((1, 'X'),): op}
+            self.assertEqual(sub_operators, expected)
+
+            op = (QubitOperator('X0 X1', 0.1) + QubitOperator('X1 X2', 2.j)
+                  + QubitOperator('Y2 Z3', 3.) + QubitOperator('X3 Z4', 5.))
+            sub_operators = group_into_tensor_product_basis_sets(op, seed=seed)
+            expected = {
+                ((0, 'X'), (1, 'X'), (2, 'X'),
+                 (3, 'X'), (4, 'Z')): (QubitOperator('X0 X1', 0.1)
+                                       + QubitOperator('X1 X2', 2.j)
+                                       + QubitOperator('X3 Z4', 5.)),
+                ((2, 'Y'), (3, 'Z')): QubitOperator('Y2 Z3', 3.)
+            }
+            self.assertEqual(sub_operators, expected)
+
+    def test_empty_qubit_operator(self):
+        sub_operators = group_into_tensor_product_basis_sets(QubitOperator())
+        self.assertTrue(sub_operators == {})
+
+    def test_fermion_operator_bad_type(self):
+        with self.assertRaises(TypeError):
+            _ = group_into_tensor_product_basis_sets(FermionOperator())
+
+    def test_boson_operator_bad_type(self):
+        with self.assertRaises(TypeError):
+            _ = group_into_tensor_product_basis_sets(BosonOperator())
+
+    def test_none_bad_type(self):
+        with self.assertRaises(TypeError):
+            _ = group_into_tensor_product_basis_sets(None)