diff --git a/symmer/operators/base.py b/symmer/operators/base.py
index f39ff269..c2f16aeb 100644
--- a/symmer/operators/base.py
+++ b/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,36 +756,6 @@ 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
@@ -788,6 +763,10 @@ def _multiply_by_operator(self,
         """ 
         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)
diff --git a/symmer/operators/independent_op.py b/symmer/operators/independent_op.py
index 1e339a23..afaa2a1f 100644
--- a/symmer/operators/independent_op.py
+++ b/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))
diff --git a/symmer/operators/noncontextual_op.py b/symmer/operators/noncontextual_op.py
index c24749eb..9b1fac36 100644
--- a/symmer/operators/noncontextual_op.py
+++ b/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 ###############################
 ###############################################################################
diff --git a/symmer/operators/utils.py b/symmer/operators/utils.py
index 31174c65..a9eb7f9a 100644
--- a/symmer/operators/utils.py
+++ b/symmer/operators/utils.py
@@ -3,6 +3,43 @@
 from typing import Tuple, Dict
 from openfermion import QubitOperator
 from qiskit.quantum_info import SparsePauliOp
+from qiskit._accelerate.sparse_pauli_op import unordered_unique
+import numba as nb
+
+@nb.njit(fastmath=True, parallel=True, cache=True)
+def matmul_GF2(A: np.array, B: np.array) -> np.array:
+    """
+    custom GF(2) multiplication using numba. i.e. C = (A@B) mod 2.
+    
+    Note function uses fact that multiplication over GF2 doesn't require sums for each matrix element in C
+    instead it uses and AND operation (same as elementwise multiplicaiton of row,col pairs) 
+    followed by XOR operation (same as taking sum of row,col multiplied pairs)
+    
+    Args:
+        A (np.array): numpy boolean array
+        B (np.array): numpy boolean array
+    Returns:
+        C (np.array): numpy boolean array of (A@B) mod 2
+    """
+    C = np.empty((A.shape[0], B.shape[1]), dtype=np.bool_)
+    for i in nb.prange(C.shape[0]):
+        for j in range(C.shape[1]):
+
+            ## logical version 1 (slower) 
+            # C[i, j] = bool(np.sum(np.logical_and(A[i, :], B[:, j]))%2)
+
+            # # logical version 2 (slower) 
+            # parity = False
+            # for bit in np.logical_and(A[i, :], B[:, j]):
+            #     parity^=bit
+            # C[i, j] = parity
+
+            ## faster loop
+            acc = False
+            for k in range(B.shape[0]):
+                acc ^= A[i, k] & B[k, j]
+            C[i, j] = acc
+    return C
 
 def symplectic_to_string(symp_vec) -> str:
     """
@@ -171,28 +208,38 @@ def symplectic_cleanup(
     Returns: 
         Reduced symplectic matrix and reduced coefficient vector.
     """
-    # order lexicographically using a fast void view implementation...
-    # this scales to large numbers of qubits more favourably than np.lexsort
-    symp_matrix_view = np.ascontiguousarray(symp_matrix).view(
-        np.dtype((np.void, symp_matrix.dtype.itemsize * symp_matrix.shape[1]))
-    )
-    re_order_indices = np.argsort(symp_matrix_view.ravel())
-    # sort the symplectic matrix and vector of coefficients accordingly
-    sorted_terms = symp_matrix[re_order_indices]
-    sorted_coeff = coeff_vec[re_order_indices]
-    # unique terms are those with non-zero entries in the adjacent row difference array
-    diff_adjacent = np.diff(sorted_terms, axis=0)
-    mask_unique_terms = np.append(True, np.any(diff_adjacent, axis=1))
-    reduced_symp_matrix = sorted_terms[mask_unique_terms]
-    # mask the term indices such that those which are skipped are summed under np.reduceat
-    summing_indices = np.arange(symp_matrix.shape[0])[mask_unique_terms]
-    reduced_coeff_vec = np.add.reduceat(sorted_coeff, summing_indices, axis=0)
-    # if a zero threshold is specified terms with sufficiently small coefficient will be dropped
+    # # order lexicographically using a fast void view implementation...
+    # # this scales to large numbers of qubits more favourably than np.lexsort
+    # symp_matrix_view = np.ascontiguousarray(symp_matrix).view(
+    #     np.dtype((np.void, symp_matrix.dtype.itemsize * symp_matrix.shape[1]))
+    # )
+    # re_order_indices = np.argsort(symp_matrix_view.ravel())
+    # # sort the symplectic matrix and vector of coefficients accordingly
+    # sorted_terms = symp_matrix[re_order_indices]
+    # sorted_coeff = coeff_vec[re_order_indices]
+    # # unique terms are those with non-zero entries in the adjacent row difference array
+    # diff_adjacent = np.diff(sorted_terms, axis=0)
+    # mask_unique_terms = np.append(True, np.any(diff_adjacent, axis=1))
+    # reduced_symp_matrix = sorted_terms[mask_unique_terms]
+    # # mask the term indices such that those which are skipped are summed under np.reduceat
+    # summing_indices = np.arange(symp_matrix.shape[0])[mask_unique_terms]
+    # reduced_coeff_vec = np.add.reduceat(sorted_coeff, summing_indices, axis=0)
+    # # if a zero threshold is specified terms with sufficiently small coefficient will be dropped
+    # if zero_threshold is not None:
+    #     mask_nonzero = abs(reduced_coeff_vec)>zero_threshold
+    #     reduced_symp_matrix = reduced_symp_matrix[mask_nonzero]
+    #     reduced_coeff_vec = reduced_coeff_vec[mask_nonzero]
+
+    # return reduced_symp_matrix, reduced_coeff_vec
+
+    unique_locations, inverse_map = unordered_unique(symp_matrix.astype('uint16'))
+    reduced_symp_matrix = symp_matrix[unique_locations]
+    reduced_coeff_vec = np.zeros(unique_locations.shape[0], dtype=complex)
+    np.add.at(reduced_coeff_vec, inverse_map, coeff_vec)
     if zero_threshold is not None:
         mask_nonzero = abs(reduced_coeff_vec)>zero_threshold
         reduced_symp_matrix = reduced_symp_matrix[mask_nonzero]
         reduced_coeff_vec = reduced_coeff_vec[mask_nonzero]
-
     return reduced_symp_matrix, reduced_coeff_vec
 
 def random_symplectic_matrix(n_qubits,n_terms, diagonal=False, density=0.3):
diff --git a/tests/test_operators/test_noncontextual_op.py b/tests/test_operators/test_noncontextual_op.py
index d93ff2bf..1888b3ca 100644
--- a/tests/test_operators/test_noncontextual_op.py
+++ b/tests/test_operators/test_noncontextual_op.py
@@ -240,8 +240,7 @@ def test_noncon_noncommuting_Z2():
 
 def test_noncontextual_objective_function():
     H_noncon = NoncontextualOp.from_dictionary(noncon_problem['H_dict'])
-    H_noncon.symmetry_generators = H_noncon.symmetry_generators.sort()
-    nu = [1, -1,  -1]
+    nu = [-1, -1,  +1]
     e_noncon = H_noncon.get_energy(nu)
     assert np.isclose(e_noncon, noncon_problem['E'])
 
@@ -301,7 +300,7 @@ def test_noncon_state():
     """
     Check that state generated by NoncontexualOp matches the noncontextual energy
     """
-    for _ in range(10):
+    for _ in range(50):
         Hnc = NoncontextualOp.random(6)
         Hnc.solve()
         state, nu = Hnc.noncon_state(UP_method='LCU')
diff --git a/tests/test_projection/test_contextual_subspace.py b/tests/test_projection/test_contextual_subspace.py
index 24bd448e..33d60482 100644
--- a/tests/test_projection/test_contextual_subspace.py
+++ b/tests/test_projection/test_contextual_subspace.py
@@ -45,7 +45,7 @@ def test_manual_stabilizers():
     H_cs = CS.project_onto_subspace()
     assert CS.n_qubits_in_subspace == 3
     assert H_cs.n_qubits == 3
-    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0004
+    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0005
     
 def test_update_stabilizers_aux_preserving():
     CS = ContextualSubspace(H_taper, noncontextual_strategy='SingleSweep_magnitude')
@@ -53,7 +53,7 @@ def test_update_stabilizers_aux_preserving():
     H_cs = CS.project_onto_subspace()
     assert CS.n_qubits_in_subspace == 3
     assert H_cs.n_qubits == 3
-    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0004
+    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0005
 
 def test_update_stabilizers_unrecognised_strategy():
     CS = ContextualSubspace(H_taper, noncontextual_strategy='SingleSweep_magnitude')
@@ -89,7 +89,7 @@ def test_StabilizeFirst_strategy_correct_usage():
     CS.update_stabilizers(3, aux_operator=CC_taper, strategy='aux_preserving')
     H_cs = CS.project_onto_subspace()
     assert H_cs.n_qubits == 3
-    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0004
+    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0005
 
 @pytest.mark.parametrize("ref_state", [QT.tapered_ref_state, QT.tapered_ref_state.state_matrix[0]])
 def test_reference_state(ref_state):
@@ -100,14 +100,14 @@ def test_reference_state(ref_state):
     CS.update_stabilizers(3, aux_operator=CC_taper, strategy='aux_preserving')
     H_cs = CS.project_onto_subspace()
     assert H_cs.n_qubits == 3
-    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0004
+    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0005
 
 def test_StabilizeFirst_strategy_correct_usage():
     CS = ContextualSubspace(H_taper, noncontextual_strategy='StabilizeFirst')
     CS.update_stabilizers(3, aux_operator=CC_taper, strategy='aux_preserving')
     H_cs = CS.project_onto_subspace()
     assert H_cs.n_qubits == 3
-    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0004
+    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0005
 
 def test_project_auxiliary_operator():
     CS = ContextualSubspace(H_taper, noncontextual_strategy='SingleSweep_magnitude')
@@ -116,7 +116,7 @@ def test_project_auxiliary_operator():
     H_cs = CS.project_onto_subspace()
     CC_cs = CS.project_onto_subspace(operator_to_project=CC_taper)
     assert CC_cs.n_qubits == 3
-    assert abs(H_cs.expval(trotter(CC_cs*1j, trotnum=10) * QuantumState([0,0,0])) - fci_energy) < 0.0004
+    assert abs(H_cs.expval(trotter(CC_cs*1j, trotnum=10) * QuantumState([0,0,0])) - fci_energy) < 0.0005
 
 def test_no_aux_operator_provided():
     CS = ContextualSubspace(H_taper, noncontextual_strategy='SingleSweep_magnitude')
@@ -139,7 +139,7 @@ def test_unitary_partitioning_method(up_method):
     CS.update_stabilizers(3, aux_operator=CC_taper, strategy='aux_preserving')
     H_cs = CS.project_onto_subspace()
     assert H_cs.n_qubits == 3
-    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0004
+    assert abs(exact_gs_energy(H_cs.to_sparse_matrix)[0] - fci_energy) < 0.0005
 
 @pytest.mark.parametrize("up_method", ['LCU', 'seq_rot'])
 def test_project_state_onto_subspace(up_method):