Merge 74f4729 into a626465

src/openfermion/ops/_majorana_operator.py

@@ -15,6 +15,8 @@
 from __future__ import division
 from future.utils import viewkeys
 
+import itertools
+
 import numpy
 
 
@@ -88,6 +90,47 @@ def commutes_with(self, other):
                                            list(other.terms.keys())[0])
         return self*other == other*self
 
+    def with_basis_rotated_by(self, transformation_matrix):
+        r"""Change to a basis of new Majorana operators.
+
+        The input to this method is a real orthogonal matrix :math:`O`.
+        It returns a new MajoranaOperator which is equivalent to the old one
+        but rewritten in terms of a new basis of Majorana operators.
+        Let the original Majorana operators be labeled by 
+        :math:`\gamma_i` and the new operators be labeled by
+        :math:`\tilde{\gamma_i}`. Then they are related by the equation
+
+        .. math::
+
+            O
+            \begin{pmatrix}
+                 \gamma_1 \\
+                 \vdots \\
+                 \gamma_N \\
+            \end{pmatrix}
+            =
+            \begin{pmatrix}
+                 \tilde{\gamma_1} \\
+                 \vdots \\
+                 \tilde{\gamma_N} \\
+            \end{pmatrix}.
+
+        Args:
+            transformation_matrix: A real orthogonal matrix representing
+                the basis transformation.
+        Returns:
+            The rotated operator.
+        """
+        if not _is_real_orthogonal(transformation_matrix):
+            raise ValueError("Transformation matrix is not real orthogonal.")
+
+        rotated_op = MajoranaOperator()
+        for term, coeff in self.terms.items():
+            rotated_term = _rotate_basis(term, transformation_matrix)
+            rotated_term *= coeff
+            rotated_op += rotated_term
+        return rotated_op
+
     def __iadd__(self, other):
         if not isinstance(other, type(self)):
             return NotImplemented
@@ -326,3 +369,21 @@ def _majorana_terms_commute(term_a, term_b):
             j += 1
     parity = (len(term_a)*len(term_b) - intersection) % 2
     return not parity
+
+
+def _rotate_basis(term, transformation_matrix):
+    n = transformation_matrix.shape[0]
+    rotated_op = MajoranaOperator()
+    for tup in itertools.product(range(n), repeat=len(term)):
+        coeff = 1.0
+        for i, j in zip(term, tup):
+            coeff *= transformation_matrix[j, i]
+        rotated_op += MajoranaOperator(tup, coeff)
+    return rotated_op
+
+
+def _is_real_orthogonal(matrix):
+    n, m = matrix.shape
+    return (n == m
+            and numpy.allclose(numpy.imag(matrix), 0.0)
+            and numpy.allclose(numpy.dot(matrix.T, matrix), numpy.eye(n)))

src/openfermion/ops/_majorana_operator_test.py

@@ -11,6 +11,7 @@
 #   limitations under the License.
 
 import pytest
+import numpy
 
 from openfermion import MajoranaOperator
 
@@ -52,6 +53,33 @@ def test_majorana_operator_commutes_with():
         _ = e.commutes_with(0)
 
 
+def test_majorana_operator_with_basis_rotated_by():
+    H = numpy.array([[1, 1], [1, -1]]) / numpy.sqrt(2)
+
+    a = MajoranaOperator((0, 1), 2.0)
+    op = a.with_basis_rotated_by(H)
+    assert op == MajoranaOperator.from_dict({(0, 1): -2.0})
+
+    b = MajoranaOperator((0,), 2.0)
+    op = b.with_basis_rotated_by(H)
+    assert op == MajoranaOperator.from_dict({(0,): numpy.sqrt(2),
+                                             (1,): numpy.sqrt(2)})
+
+    c = MajoranaOperator((1,), 2.0)
+    op = c.with_basis_rotated_by(H)
+    assert op == MajoranaOperator.from_dict({(0,): numpy.sqrt(2),
+                                             (1,): -numpy.sqrt(2)})
+
+    P = numpy.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
+
+    d = MajoranaOperator((0, 1, 2)) + MajoranaOperator((1, 2))
+    op = d.with_basis_rotated_by(P)
+    assert op == MajoranaOperator.from_dict({(0, 1, 2): 1.0,
+                                             (0, 1): 1.0})
+
+    with pytest.raises(ValueError):
+        _ = a.with_basis_rotated_by(2 * H)
+
 
 def test_majorana_operator_add_subtract():
     a = MajoranaOperator((0, 2, 3), -1.25)