updated to_sparse_matrix

UCL-CCS · Oct 8, 2024 · dc8cf6f · dc8cf6f
1 parent 7fe767d
commit dc8cf6f
Show file tree

Hide file tree

Showing 3 changed files with 39 additions and 22 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/symmer/.DS_Store b/symmer/.DS_Store
diff --git a/symmer/operators/base.py b/symmer/operators/base.py
@@ -21,6 +21,11 @@
 from qiskit.quantum_info import SparsePauliOp
 warnings.simplefilter('always', UserWarning)
 
+from qiskit._accelerate.sparse_pauli_op import (
+    ZXPaulis,
+    to_matrix_sparse,
+)
+
 from numba.core.errors import NumbaDeprecationWarning, NumbaPendingDeprecationWarning
 warnings.simplefilter('ignore', category=NumbaDeprecationWarning)
 warnings.simplefilter('ignore', category=NumbaPendingDeprecationWarning)
@@ -1464,33 +1469,45 @@ def to_sparse_matrix(self) -> csr_matrix:
         if self.n_qubits == 0:
             return csr_matrix(self.coeff_vec)
 
-        if self.n_qubits>15:
-            from symmer.utils import get_sparse_matrix_large_pauliwordop
-            sparse_matrix = get_sparse_matrix_large_pauliwordop(self)
-            return sparse_matrix
-        else:
-            x_int = binary_array_to_int(self.X_block).reshape(-1, 1)
-            z_int = binary_array_to_int(self.Z_block).reshape(-1, 1)
+        # if self.n_qubits>15:
+        #     from symmer.utils import get_sparse_matrix_large_pauliwordop
+        #     sparse_matrix = get_sparse_matrix_large_pauliwordop(self)
+        #     return sparse_matrix
+        # else:
+        #     x_int = binary_array_to_int(self.X_block).reshape(-1, 1)
+        #     z_int = binary_array_to_int(self.Z_block).reshape(-1, 1)
 
-            Y_number = np.sum(np.bitwise_and(self.X_block, self.Z_block), axis=1)
-            global_phase = (-1j) ** Y_number
+        #     Y_number = np.sum(np.bitwise_and(self.X_block, self.Z_block), axis=1)
+        #     global_phase = (-1j) ** Y_number
 
-            dimension = 2 ** self.n_qubits
-            row_ind = np.repeat(np.arange(dimension).reshape(1, -1), self.X_block.shape[0], axis=0)
-            col_ind = np.bitwise_xor(row_ind, x_int)
+        #     dimension = 2 ** self.n_qubits
+        #     row_ind = np.repeat(np.arange(dimension).reshape(1, -1), self.X_block.shape[0], axis=0)
+        #     col_ind = np.bitwise_xor(row_ind, x_int)
 
-            row_inds_and_Zint = np.bitwise_and(row_ind, z_int)
-            vals = global_phase.reshape(-1, 1) * (-1) ** (
-                        count1_in_int_bitstring(row_inds_and_Zint) % 2)  # .astype(complex))
+        #     row_inds_and_Zint = np.bitwise_and(row_ind, z_int)
+        #     vals = global_phase.reshape(-1, 1) * (-1) ** (
+        #                 count1_in_int_bitstring(row_inds_and_Zint) % 2)  # .astype(complex))
 
-            values_and_coeff = np.einsum('ij,i->ij', vals, self.coeff_vec)
+        #     values_and_coeff = np.einsum('ij,i->ij', vals, self.coeff_vec)
 
-            sparse_matrix = csr_matrix(
-                (values_and_coeff.flatten(), (row_ind.flatten(), col_ind.flatten())),
-                shape=(dimension, dimension),
-                dtype=complex
-            )
-            return sparse_matrix
+        #     sparse_matrix = csr_matrix(
+        #         (values_and_coeff.flatten(), (row_ind.flatten(), col_ind.flatten())),
+        #         shape=(dimension, dimension),
+        #         dtype=complex
+        #     )
+        #     return sparse_matrix
+
+        phase = np.zeros(self.n_terms, dtype=np.uint8)
+        zx = ZXPaulis(
+                    self.X_block[:,::-1],
+                    self.Z_block[:,::-1],
+                    phase,
+                    self.coeff_vec,
+                )
+
+        data, indices, indptr = to_matrix_sparse(zx, force_serial=False)
+        side = 1 << self.n_qubits
+        return csr_matrix((data, indices, indptr), shape=(side, side))
 
     def conjugate_op(self, R: 'PauliwordOp') -> 'PauliwordOp':
         """