Faster binary mult

symmer/operators/base.py

@@ -4,12 +4,15 @@
 import networkx as nx
 import matplotlib.pyplot as plt
 from symmer import process
-from symmer.operators.utils import *
-from symmer.operators.utils import _cref_binary
+from symmer.operators.utils import (
+    matmul_GF2, random_symplectic_matrix, string_to_symplectic, QubitOperator_to_dict, SparsePauliOp_to_dict,
+    symplectic_to_string, cref_binary, check_independent, check_jordan_independent, symplectic_cleanup,
+    check_adjmat_noncontextual, symplectic_to_openfermion, binary_array_to_int, count1_in_int_bitstring
+)
 from tqdm.auto import tqdm
 from copy import deepcopy
 from functools import reduce
-from typing import List, Union, Optional
+from typing import List, Union, Optional, Dict, Tuple
 from numbers import Number
 from cached_property import cached_property
 from scipy.stats import unitary_group
@@ -458,6 +461,8 @@ def sort(self,
         """
         if by == 'magnitude':
             sort_order = np.argsort(-abs(self.coeff_vec))
+        elif by == 'lex':
+            sort_order = np.lexsort(self.symp_matrix.T)
         elif by == 'weight':
             sort_order = np.argsort(-np.sum(self.symp_matrix.astype(int), axis=1))
         elif by == 'support':
@@ -641,8 +646,8 @@ def __eq__(self,
         Returns:
             bool: True if the two PauliwordOp objects are equal, False otherwise.
         """
-        check_1 = self.cleanup()
-        check_2 = Pword.cleanup()
+        check_1 = self.cleanup().sort('lex')
+        check_2 = Pword.cleanup().sort('lex')
         if check_1.n_qubits != check_2.n_qubits:
             raise ValueError('Operators defined over differing numbers of qubits.')
         elif check_1.n_terms != check_2.n_terms:
@@ -751,43 +756,17 @@ def multiply_by_constant(self,
         """
         return PauliwordOp(self.symp_matrix, self.coeff_vec*const)
 
-    def _mul_symplectic(self, 
-            symp_vec: np.array, 
-            coeff: complex, 
-            Y_count_in: np.array
-        ) -> Tuple[np.array, np.array]:
-        """ 
-        Performs Pauli multiplication with phases at the level of the symplectic 
-        matrices to avoid superfluous PauliwordOp initializations. The phase compensation 
-        is implemented as per https://doi.org/10.1103/PhysRevA.68.042318.
-
-        Args:
-            symp_vec (np.array): The symplectic vector representing the Pauli operator to be multiplied with.
-            coeff (complex): The complex coefficient of the Pauli operator.
-            Y_count_in (np.array): The Y-counts of the input Pauli operator.
-
-        Returns:
-            Tuple[np.array, np.array]: The resulting symplectic vector and coefficient vector after multiplication.
-        """
-        # phaseless multiplication is binary addition in symplectic representation
-        phaseless_prod = np.bitwise_xor(self.symp_matrix, symp_vec)
-        # phase is determined by Y counts plus additional sign flip
-        Y_count_out = np.sum(np.bitwise_and(*np.hsplit(phaseless_prod,2)), axis=1)
-        sign_change = (-1) ** (
-            np.sum(np.bitwise_and(self.X_block, np.hsplit(symp_vec,2)[1]), axis=1) % 2
-        ) # mod 2 as only care about parity
-        # final phase modification
-        phase_mod = sign_change * (1j) ** ((3*Y_count_in + Y_count_out) % 4) # mod 4 as roots of unity
-        coeff_vec = phase_mod * self.coeff_vec * coeff
-        return phaseless_prod, coeff_vec
-
     def _multiply_by_operator(self, 
             PwordOp: Union["PauliwordOp", "QuantumState", complex],
             zero_threshold: float = 1e-15
         ) -> "PauliwordOp":
         """ 
         Right-multiplication of this PauliwordOp by another PauliwordOp or QuantumState ket.
 
+        Performs Pauli multiplication with phases at the level of the symplectic 
+        matrices to avoid superfluous PauliwordOp initializations. The phase compensation 
+        is implemented as per https://doi.org/10.1103/PhysRevA.68.042318.
+
         Args:
             PwordOp (Union["PauliwordOp", "QuantumState", complex]): The operator or ket to multiply by.
             zero_threshold (float): The threshold below which coefficients are considered negligible.
@@ -796,27 +775,18 @@ def _multiply_by_operator(self,
             PauliwordOp: The resulting PauliwordOp after multiplication.
         """
         assert (self.n_qubits == PwordOp.n_qubits), 'PauliwordOps defined for different number of qubits'
-
-        if PwordOp.n_terms == 1:
-            # no cleanup if multiplying by a single term (faster)
-            symp_stack, coeff_stack = self._mul_symplectic(
-                symp_vec=PwordOp.symp_matrix, 
-                coeff=PwordOp.coeff_vec, 
-                Y_count_in=self.Y_count+PwordOp.Y_count
-            )
-            pauli_mult_out = PauliwordOp(symp_stack, coeff_stack)
-        else:
-            # multiplication is performed at the symplectic level, before being stacked and cleaned
-            symp_stack, coeff_stack = zip(
-                *[self._mul_symplectic(symp_vec=symp_vec, coeff=coeff, Y_count_in=self.Y_count+Y_count) 
-                for symp_vec, coeff, Y_count in zip(PwordOp.symp_matrix, PwordOp.coeff_vec, PwordOp.Y_count)]
+        Q_symp_matrix = PwordOp.symp_matrix.reshape([PwordOp.n_terms, 1, 2*PwordOp.n_qubits])
+        termwise_phaseless_prod = self.symp_matrix ^ Q_symp_matrix
+        Y_count_in  = self.Y_count + PwordOp.Y_count.reshape(-1,1)
+        Y_count_out = np.sum(termwise_phaseless_prod[:,:,:PwordOp.n_qubits] & termwise_phaseless_prod[:,:,PwordOp.n_qubits:], axis=2)
+        sign_change = (-1) ** (np.sum(self.X_block & Q_symp_matrix[:,:,PwordOp.n_qubits:], axis=2) % 2)
+        phase_mod = sign_change * (1j) ** ((3*Y_count_in + Y_count_out) % 4)
+        return PauliwordOp(
+            *symplectic_cleanup(
+                termwise_phaseless_prod.reshape(-1,2*self.n_qubits), 
+                (phase_mod * np.outer(self.coeff_vec, PwordOp.coeff_vec).T).reshape(-1), zero_threshold=zero_threshold
             )
-            pauli_mult_out = PauliwordOp(
-                *symplectic_cleanup(
-                    np.vstack(symp_stack), np.hstack(coeff_stack), zero_threshold=zero_threshold
-                )
-            )
-        return pauli_mult_out
+        )
 
     def expval(self, psi: "QuantumState") -> complex:
         """ 
@@ -990,8 +960,10 @@ def commutes_termwise(self,
         # return np.logical_not(adjacency_matrix.toarray())
 
         ### dense code
-        Omega_PwordOp_symp = np.hstack((PwordOp.Z_block,  PwordOp.X_block)).astype(int)
-        return (self.symp_matrix @ Omega_PwordOp_symp.T) % 2 == 0
+        # Omega_PwordOp_symp = np.hstack((PwordOp.Z_block,  PwordOp.X_block)).astype(int)
+        # return (self.symp_matrix @ Omega_PwordOp_symp.T) % 2 == 0
+
+        return ~matmul_GF2(self.symp_matrix, np.hstack((PwordOp.Z_block,  PwordOp.X_block)).T)
 
     def anticommutes_termwise(self,
             PwordOp: "PauliwordOp"
@@ -1218,8 +1190,8 @@ def tensor(self, right_op: "PauliwordOp") -> "PauliwordOp":
         Returns:
             PauliwordOp: The resulting Pauli operator after the tensor product.
         """
-        identity_block_right = np.zeros([right_op.n_terms, self.n_qubits]).astype(int)
-        identity_block_left  = np.zeros([self.n_terms,  right_op.n_qubits]).astype(int)
+        identity_block_right = np.zeros([right_op.n_terms, self.n_qubits], dtype=bool)#.astype(int)
+        identity_block_left  = np.zeros([self.n_terms,  right_op.n_qubits], dtype=bool)#.astype(int)
         padded_left_symp = np.hstack([self.X_block, identity_block_left, self.Z_block, identity_block_left])
         padded_right_symp = np.hstack([identity_block_right, right_op.X_block, identity_block_right, right_op.Z_block])
         left_factor = PauliwordOp(padded_left_symp, self.coeff_vec)

symmer/operators/independent_op.py

@@ -137,7 +137,7 @@ def symmetry_generators(cls,
                 S_commuting = S.largest_clique(edge_relation='C')    
             else:
                 # greedy graph-colouring approach when symmetry basis is large
-                S_commuting = S.clique_cover(edge_relation='C')[0]
+                S_commuting = S.clique_cover(edge_relation='C', strategy='independent_set')[0]
                 warnings.warn('Greedy method may identify non-optimal commuting symmetry terms; might be able to taper again.')
 
             return cls(S_commuting.symp_matrix, np.ones(S_commuting.n_terms, dtype=complex))

symmer/operators/noncontextual_op.py

@@ -249,65 +249,70 @@ def _from_generators_noncontextual_op(cls,
             _, noncontextual_terms_mask = H.generator_reconstruction(generators, override_independence_check=True)
 
         return cls.from_PauliwordOp(H[noncontextual_terms_mask])
-
     @classmethod
     def random(cls,
             n_qubits: int,
             n_cliques:Optional[int]=3,
             complex_coeffs:Optional[bool]=False,
             n_commuting_terms:Optional[int]=None,
+            apply_clifford: Optional[bool] = True,
         ) -> "NoncontextualOp":
         """
         Generate a random Noncontextual operator with normally distributed coefficients.
         Note to maximise size choose number of n_cliques to be 3 (and for 2<= n_cliques <= 5 the operator
         will be larger than only using commuting generators).
+
         WARNING: this function can generates an exponentially large Hamiltonian unless n_terms set.
-        size when NOT set is: n_cliques * [ 2**(n_qubits -  int(np.ceil((n_cliques - 1) / 2))) ]
+        size when NOT set is: n_cliques * [ 2**(n_qubits -  int(np.ceil((n_cliques - 1) / 2))) ] + n_commuting_terms
 
-        Note: The number of terms in output will be: n_cliques*n_commuting_terms
+        Note: The number of terms in output will be: n_cliques*n_commuting_terms + n_commuting_terms
+              apart from when n_commuting_terms=0 in which case = n_cliques
 
         Args:
             n_qubits (int): Number of qubits noncontextual operator defined on
             n_cliques (int): Number of cliques representives in operator
             complex_coeffs (bool): Whether to generate complex coefficients (default: True).
             n_commuting_terms (int): Optional int for number of commuting terms. if not set then it will be: 2**(n_qubits -  int(np.ceil((n_cliques - 1) / 2))) (i.e. exponentially large)
-
+            apply_clifford (bool): Whether to apply clifford before output (to get rid of structure used to build Hnc)
         Returns:
             NoncontextualOp: A random NoncontextualOp object.
         """
-        assert n_cliques > 1, 'number of cliques must be set to 2 or more (cannot have one anticommuting term)'
+        assert n_cliques > 1 or n_cliques == 0, 'number of cliques must be zero or set to 2 or more (cannot have one anticommuting term)'
         n_clique_qubits = int(np.ceil((n_cliques - 1) / 2))
-        assert n_clique_qubits <= n_qubits, 'cannot have {n_cliques} anticommuting cliques on {n_qubits} qubits'
+        assert n_clique_qubits <= n_qubits, f'cannot have {n_cliques} anticommuting cliques on {n_qubits} qubits'
 
         remaining_qubits = n_qubits - n_clique_qubits
 
         if n_commuting_terms:
-            assert n_commuting_terms<= 2**remaining_qubits, f'cannot have {n_commuting_terms} commuting operators on {remaining_qubits} qubits'
+            assert n_commuting_terms <= 2 ** remaining_qubits, f'cannot have {n_commuting_terms} commuting operators on {remaining_qubits} qubits'
+        elif n_qubits == n_clique_qubits:
+            n_commuting_terms = 0
 
-        if remaining_qubits>=1:
-            if n_commuting_terms==None:
+        if remaining_qubits >= 1:
+            if n_commuting_terms == None:
                 n_commuting_terms = 2 ** (remaining_qubits)
                 XZ_block = (((np.arange(n_commuting_terms)[:, None] & (1 << np.arange(2 * remaining_qubits))[
-                                                                            ::-1])) > 0).astype(bool)
+                                                                      ::-1])) > 0).astype(bool)
+            elif n_commuting_terms == 0:
+                XZ_block = np.zeros(2 * remaining_qubits, dtype=bool).reshape([1, -1])
             else:
                 # randomly chooise Z bitstrings in symp matrix:
                 indices = np.unique(np.random.random_integers(0,
-                                                              high=2**remaining_qubits-1,
-                                                              size=10*n_commuting_terms))
+                                                              high=2 ** remaining_qubits - 1,
+                                                              size=10 * n_commuting_terms))
                 while len(indices) < n_commuting_terms:
                     indices = np.unique(np.append(indices,
                                                   np.unique(np.random.random_integers(0,
                                                                                       high=2 ** remaining_qubits - 1,
-                                                                                      size=10*n_commuting_terms)))
+                                                                                      size=10 * n_commuting_terms)))
                                         )
 
                 indices = indices[:n_commuting_terms]
                 XZ_block = (((indices[:, None] & (1 << np.arange(2 * remaining_qubits))[
-                                                                        ::-1])) > 0).astype(bool)
+                                                 ::-1])) > 0).astype(bool)
 
         if n_cliques == 0:
             H_nc = PauliwordOp(XZ_block, np.ones(XZ_block.shape[0]))
-
         else:
             AC = random_anitcomm_2n_1_PauliwordOp(n_clique_qubits,
                                                   apply_clifford=True)[:n_cliques]
@@ -318,22 +323,28 @@ def random(cls,
                 diag_H = PauliwordOp.from_list(['I' * remaining_qubits])
 
             AC_full = PauliwordOp.from_list(['I' * remaining_qubits]).tensor(AC)
+
             H_sym = diag_H.tensor(PauliwordOp.from_list(['I' * n_clique_qubits]))
-            H_nc = AC_full * H_sym
-            assert AC.n_terms * n_commuting_terms == H_nc.n_terms, 'operator not largest it can be'
+            if n_commuting_terms > 0:
+                H_nc = AC_full * H_sym + H_sym
+                assert n_commuting_terms*n_cliques + n_commuting_terms == H_nc.n_terms, 'operator not largest it can be'
+            else:
+                H_nc = AC_full * H_sym + H_sym
+                assert AC.n_terms + 1 == H_nc.n_terms, 'operator not largest it can be'
 
         coeff_vec = np.random.randn(H_nc.n_terms).astype(complex)
         if complex_coeffs:
             coeff_vec += 1j * np.random.randn(H_nc.n_terms)
 
-        # apply clifford rotations to get rid of some of generation structure
-        U_cliff_rotations = []
-        for _ in range(n_qubits * 5):
-            P_rand = PauliwordOp.random(H_nc.n_qubits, n_terms=1)
-            P_rand.coeff_vec = [1]
-            U_cliff_rotations.append((P_rand, np.random.choice([np.pi/2, -np.pi/2])))
+        if apply_clifford:
+            # apply clifford rotations to get rid of some of generation structure
+            U_cliff_rotations = []
+            for _ in range(n_qubits * 5):
+                P_rand = PauliwordOp.random(H_nc.n_qubits, n_terms=1)
+                P_rand.coeff_vec = [1]
+                U_cliff_rotations.append((P_rand, (np.pi/2)*np.random.choice([1,3])))
 
-        H_nc = H_nc.perform_rotations(U_cliff_rotations)
+            H_nc = H_nc.perform_rotations(U_cliff_rotations)
 
         return cls(H_nc.symp_matrix, coeff_vec)
 
@@ -575,7 +586,7 @@ def solve(self,
         if self.n_cliques > 0:
             self.update_clique_representative_operator()
 
-    def noncon_state(self, UP_method:Optional[str]= 'LCU') -> Tuple[QuantumState, np.array]:
+    def noncon_state(self, UP_method = 'LCU') -> Tuple[QuantumState, np.array]:
         """
         Method to generate noncontextual state for current symmetry generators assignments. Note by default
         UP_method is set to LCU as this avoids generating exponentially large states (which seq_rot can do!)
@@ -589,55 +600,43 @@ def noncon_state(self, UP_method:Optional[str]= 'LCU') -> Tuple[QuantumState, np
 
         """
         nu_assignment = self.symmetry_generators.coeff_vec.copy()
-
         ## update clique coeffs from nu assignment!
         _, si = self.get_symmetry_contributions(nu_assignment)
         self.clique_operator.coeff_vec = si
-
         assert UP_method in ['LCU', 'seq_rot']
-
         if UP_method == 'LCU':
             rotations_LCU = self.clique_operator.R_LCU
-            Ps, rotations_LCU, gamma_l, AC_normed = self.clique_operator.unitary_partitioning(s_index=0,
-                                                                                                  up_method='LCU')
+            Ps, rotations_LCU, gamma_l, AC_normed = self.clique_operator.unitary_partitioning(s_index=0,up_method='LCU')
         else:
-            Ps, rotations_SEQ, gamma_l, AC_normed = self.clique_operator.unitary_partitioning(s_index=0,
-                                                                                   up_method='seq_rot')
-
+            Ps, rotations_SEQ, gamma_l, AC_normed = self.clique_operator.unitary_partitioning(s_index=0,up_method='seq_rot')
         # choose negative value for clique operator (to minimize energy)
         Ps.coeff_vec[0] = -1
-
         ### to find ground state, need to map noncontextual stabilizers to single qubit Pauli Zs
         independent_stabilizers = self.symmetry_generators + Ps
-
         # rotate onto computational basis
         independent_stabilizers.target_sqp = 'Z'
-
         rotated_stabs = independent_stabilizers.rotate_onto_single_qubit_paulis()
         clifford_rots = independent_stabilizers.stabilizer_rotations
-
         ## get stabilizer state for the rotated stabilizers
-        Z_indices = np.sum(rotated_stabs.Z_block, axis=0)
-        Z_vals = np.sum(rotated_stabs.Z_block[:, Z_indices.astype(bool)] * rotated_stabs.coeff_vec, axis=1)
-        Z_indices[Z_indices.astype(bool)] = ((Z_vals - 1) * -0.5).astype(int)
-
-        state = QuantumState(Z_indices.reshape(1, -1))
-
+        nc_vec = np.zeros(self.n_qubits, dtype=int)
+        for val,row in zip(rotated_stabs.coeff_vec, rotated_stabs.Z_block):
+            assert np.count_nonzero(row) == 1
+            nc_vec[row] = (1-val)/2
+        state = QuantumState(nc_vec)
         ## undo clifford rotations
         from symmer.evolution.exponentiation import exponentiate_single_Pop
         for op, _ in clifford_rots[::-1]:
             rot = exponentiate_single_Pop(op.multiply_by_constant(1j * np.pi / 4))
             state = rot.dagger * state
-
         ## undo unitary partitioning step
         if UP_method == 'LCU':
             state = self.clique_operator.R_LCU.dagger * state
         else:
             for op, angle in rotations_SEQ[::-1]:
                 state = exponentiate_single_Pop(op.multiply_by_constant(1j * angle / 2)).dagger * state
-
         # TODO: could return clifford and UP rotations here too!
-        return state, nu_assignment    
+        return state, nu_assignment
+
 ###############################################################################
 ################### NONCONTEXTUAL SOLVERS BELOW ###############################
 ###############################################################################