Fix convention of PolynomialTensor basis change (#805)

In general, we rotate a matrix M to M' is to apply a rotation matrix R:
M' = R @ M @ R.T
The corresponding Einstein notation is
M'^{p_1p_2} = R^{p_1}_{a_1} R^{p_2}_{b_1} M^{a_1a_2}

---------

Co-authored-by: Matthew Harrigan <mpharrigan@google.com>
Co-authored-by: Michael Hucka <mhucka@google.com>
Co-authored-by: arettig <arettig@google.com>

Author: 4 people

src/openfermion/ops/representations/polynomial_tensor.py

@@ -32,8 +32,7 @@ def general_basis_change(general_tensor, rotation_matrix, key):
     r"""Change the basis of a general interaction tensor.
 
     M'^{p_1p_2...p_n} = R^{p_1}_{a_1} R^{p_2}_{a_2} ...
-                        R^{p_n}_{a_n} M^{a_1a_2...a_n} R^{p_n}_{a_n}^T ...
-                        R^{p_2}_{a_2}^T R_{p_1}_{a_1}^T
+                        R^{p_n}_{a_n} M^{a_1a_2...a_n}
 
     where R is the rotation matrix, M is the general tensor, M' is the
     transformed general tensor, and a_k and p_k are indices. The formula uses
@@ -68,7 +67,7 @@ def general_basis_change(general_tensor, rotation_matrix, key):
 
     # Do the basis change through a single call of numpy.einsum. For example,
     # for the (1, 1, 0, 0) tensor, the call is:
-    #     numpy.einsum('abcd,aA,bB,cC,dD',
+    #     numpy.einsum('abcd,Aa,Bb,Cc,Dd',
     #                  general_tensor,
     #                  rotation_matrix.conj(),
     #                  rotation_matrix.conj(),
@@ -79,7 +78,7 @@ def general_basis_change(general_tensor, rotation_matrix, key):
     subscripts_first = ''.join(chr(ord('a') + i) for i in range(order))
 
     # The 'Aa,Bb,Cc,Dd' part of the subscripts
-    subscripts_rest = ','.join(chr(ord('a') + i) + chr(ord('A') + i) for i in range(order))
+    subscripts_rest = ','.join(chr(ord('A') + i) + chr(ord('a') + i) for i in range(order))
 
     subscripts = subscripts_first + ',' + subscripts_rest
 

src/openfermion/ops/representations/polynomial_tensor_test.py

@@ -474,6 +474,34 @@ def test_rotate_basis_reverse(self):
         polynomial_tensor.rotate_basis(rotation_matrix_reverse)
         self.assertEqual(polynomial_tensor, want_polynomial_tensor)
 
+    def test_rotate_basis_cyclic(self):
+        n_qubits = 3
+        # Cyclic permutation matrix: 0->1, 1->2, 2->0
+        rotation_matrix = numpy.zeros((n_qubits, n_qubits))
+        rotation_matrix[0, 1] = 1
+        rotation_matrix[1, 2] = 1
+        rotation_matrix[2, 0] = 1
+
+        one_body = numpy.arange(n_qubits**2, dtype=float).reshape((n_qubits, n_qubits))
+        two_body = numpy.arange(n_qubits**4, dtype=float).reshape(
+            (n_qubits, n_qubits, n_qubits, n_qubits)
+        )
+        polynomial_tensor = PolynomialTensor(
+            {(): self.constant, (1, 0): one_body, (1, 1, 0, 0): two_body}
+        )
+
+        # map from rotated index to old index rotation_map[new_ind] = old_ind
+        # R maps 1->0, 2->1, 0->2. So new index 0 comes from old 1, etc.
+        rotation_map = [1, 2, 0]
+        ref_one_body = one_body[numpy.ix_(rotation_map, rotation_map)]
+        ref_two_body = two_body[numpy.ix_(rotation_map, rotation_map, rotation_map, rotation_map)]
+        ref_polynomial_tensor = PolynomialTensor(
+            {(): self.constant, (1, 0): ref_one_body, (1, 1, 0, 0): ref_two_body}
+        )
+
+        polynomial_tensor.rotate_basis(rotation_matrix)
+        self.assertEqual(polynomial_tensor, ref_polynomial_tensor)
+
     def test_rotate_basis_quadratic_hamiltonian_real(self):
         self.do_rotate_basis_quadratic_hamiltonian(True)
 
@@ -492,7 +520,7 @@ def do_rotate_basis_quadratic_hamiltonian(self, real):
 
         # Rotate a basis where the Hamiltonian is diagonal
         _, diagonalizing_unitary, _ = quad_ham.diagonalizing_bogoliubov_transform()
-        quad_ham.rotate_basis(diagonalizing_unitary.T)
+        quad_ham.rotate_basis(diagonalizing_unitary)
 
         # Check that the rotated Hamiltonian is diagonal with the correct
         # orbital energies
@@ -515,6 +543,8 @@ def test_rotate_basis_max_order(self):
             self.assertEqual(tensor, want_tensor)
         # I originally wanted to test 25 and 26, but it turns out that
         # numpy.einsum complains "too many subscripts in einsum" before 26.
+        # Currently, the sum of the number of output labels and combined labels
+        # can't not exceed NPY_MAXDIMS(=32).
 
         for order in [27, 28]:
             with self.assertRaises(ValueError):