diff --git a/symmer/operators/anticommuting_op.py b/symmer/operators/anticommuting_op.py index fae497ee..e401e784 100644 --- a/symmer/operators/anticommuting_op.py +++ b/symmer/operators/anticommuting_op.py @@ -167,14 +167,31 @@ def unitary_partitioning(self, s_index: int=None, up_method: Optional[str]='seq_ AC_op (AntiCommutingOp): normalized clique - i.e. self == gamma_l * AC_op """ assert up_method in ['LCU', 'seq_rot'], f'unknown unitary partitioning method: {up_method}' - AC_op = self.copy() + + if s_index is None: + s_index = self.get_least_dense_term_index() + + if np.isclose(self.coeff_vec[s_index], 0): + # need to correct for s_index having zero coeff... + s_index = np.argmax(abs(self.coeff_vec)) + warnings.warn(f's indexed term has zero coeff, s_index set to {s_index} so that nonzero operator is rotated onto') + + s_index = int(s_index) + BsPs = self[s_index] + + # NOTE: term to reduce to is at the top of sym matrix i.e. s_index of ZERO now! + no_BsPs = (self - BsPs).cleanup() + if (len(no_BsPs.coeff_vec)==1 and no_BsPs.coeff_vec[0]==0): + AC_op = BsPs + else: + AC_op = BsPs.append(no_BsPs) if AC_op.n_terms == 1: - rotations = None + rotations = [] gamma_l = np.linalg.norm(AC_op.coeff_vec) AC_op.coeff_vec = AC_op.coeff_vec / gamma_l - Ps = PauliwordOp(AC_op.symp_matrix, [1]) - return Ps, rotations, gamma_l, AC_op + Ps = AC_op + return Ps, rotations, gamma_l, self.multiply_by_constant(1/gamma_l) else: @@ -183,23 +200,6 @@ def unitary_partitioning(self, s_index: int=None, up_method: Optional[str]='seq_ gamma_l = np.linalg.norm(AC_op.coeff_vec) AC_op.coeff_vec = AC_op.coeff_vec / gamma_l - if s_index is None: - s_index = self.get_least_dense_term_index() - - if s_index!=0: - # re-order so s term is ALWAYS at top of symplectic matrix and thus is index as 0! - ### assert s_index <= AC_op.n_terms-1, 's_index out of range' - AC_op.coeff_vec[[0, s_index]] = AC_op.coeff_vec[[s_index, 0]] - AC_op.symp_matrix[[0, s_index]] = AC_op.symp_matrix[[s_index, 0]] - AC_op = AntiCommutingOp(AC_op.symp_matrix, AC_op.coeff_vec) # need to reinit otherwise Z and X blocks wrong - - # assert not np.isclose(AC_op.coeff_vec[0], 0), f's_index cannot have zero coefficent: {AC_op.coeff_vec[0]}' - if np.isclose(AC_op[0].coeff_vec, 0): - # need to correct for s_index having zero coeff... then need to swap to nonzero index - non_zero_index = np.argmax(abs(AC_op.coeff_vec)) - AC_op.coeff_vec[[0, non_zero_index]] = AC_op.coeff_vec[[non_zero_index, 0]] - AC_op.symp_matrix[[0, non_zero_index]] = AC_op.symp_matrix[[non_zero_index, 0]] - if up_method=='seq_rot': if len(self.X_sk_rotations)!=0: self.X_sk_rotations = [] @@ -210,11 +210,18 @@ def unitary_partitioning(self, s_index: int=None, up_method: Optional[str]='seq_ self.R_LCU = None Ps = self.generate_LCU_operator(AC_op) - rotations = self.R_LCU + rotations = LCU_as_seq_rot(self.R_LCU) else: raise ValueError(f'unknown unitary partitioning method: {up_method}!') - return Ps, rotations, gamma_l, AC_op + return Ps, rotations, gamma_l, self.multiply_by_constant(1/gamma_l) + + def multiply_by_constant(self, constant: float) -> "AntiCommutingOp": + """ Return AntiCommutingOp under constant multiplication + """ + AC_op_copy = self.copy() + AC_op_copy.coeff_vec *= constant + return AC_op_copy def generate_LCU_operator(self, AC_op) -> PauliwordOp: """ @@ -235,74 +242,45 @@ def generate_LCU_operator(self, AC_op) -> PauliwordOp: R_LCU (PauliwordOp): PauliwordOp that is a linear combination of unitaries P_s (PauliwordOp): single PauliwordOp that has been reduced too. """ - # need to remove zero coeff terms - AC_op_cpy = AC_op.copy() - before_cleanup = AC_op_cpy.n_terms - AC_op = AC_op_cpy[np.where(abs(AC_op.coeff_vec)>1e-15)[0]] - post_cleanup = AC_op.n_terms - # AC_op = AC_op.cleanup(zero_threshold=1e-15) ## cleanup re-orders which is BAD for s_index - - - if (before_cleanup>1 and post_cleanup==1): - if AC_op.coeff_vec[0]<0: - # need to fix neg sign (use Pauli multiplication) - - # as s index defaults to 0, take the next term (in CS-VQE this will commute with symmetries)! - if np.isclose(AC_op_cpy[0].coeff_vec, 0): - # need to correct for s_index having zero coeff... then need to swap to nonzero index - non_zero_index = np.argmax(abs(AC_op_cpy.coeff_vec)) - - AC_op_cpy.coeff_vec[[0, non_zero_index]] = AC_op_cpy.coeff_vec[[non_zero_index, 0]] - AC_op_cpy.symp_matrix[[0, non_zero_index]] = AC_op_cpy.symp_matrix[[non_zero_index, 0]] + ## s_index is ensured to be in zero position in unitary_partitioning method! + ## if using function without this method need to ensure term to rotate onto is the zeroth index of AC_op! + s_index=0 + # note gamma_l norm applied on init! + Ps_LCU = PauliwordOp(AC_op.symp_matrix[s_index], [1]) + βs = AC_op.coeff_vec[s_index] - sign_correction = PauliwordOp(AC_op_cpy.symp_matrix[1],[1]) + # ∑ β_k 𝑃_k ... note this doesn't contain 𝛽_s 𝑃_s + no_βsPs = AC_op - (Ps_LCU.multiply_by_constant(βs)) - self.R_LCU = sign_correction - Ps_LCU = PauliwordOp(AC_op.symp_matrix, [1]) - else: - self.R_LCU = PauliwordOp.from_list(['I'*AC_op.n_qubits]) - Ps_LCU = PauliwordOp(AC_op.symp_matrix, AC_op.coeff_vec) - else: - s_index=0 - - # note gamma_l norm applied on init! - Ps_LCU = PauliwordOp(AC_op.symp_matrix[s_index], [1]) - βs = AC_op.coeff_vec[s_index] - - # ∑ β_k 𝑃_k ... note this doesn't contain 𝛽_s 𝑃_s - no_βsPs = AC_op - (Ps_LCU.multiply_by_constant(βs)) - - # Ω_𝑙 ∑ 𝛿_k 𝑃_k ... renormalized! - omega_l = np.linalg.norm(no_βsPs.coeff_vec) - no_βsPs.coeff_vec = no_βsPs.coeff_vec / omega_l + # Ω_𝑙 ∑ 𝛿_k 𝑃_k ... renormalized! + omega_l = np.linalg.norm(no_βsPs.coeff_vec) + no_βsPs.coeff_vec = no_βsPs.coeff_vec / omega_l - phi_n_1 = np.arccos(βs) - # require sin(𝜙_{𝑛−1}) to be positive... - if (phi_n_1 > np.pi): - phi_n_1 = 2 * np.pi - phi_n_1 + phi_n_1 = np.arccos(βs) + # require sin(𝜙_{𝑛−1}) to be positive... + if (phi_n_1 > np.pi): + phi_n_1 = 2 * np.pi - phi_n_1 - alpha = phi_n_1 - I_term = 'I' * Ps_LCU.n_qubits - self.R_LCU = PauliwordOp.from_dictionary({I_term: np.cos(alpha / 2)}) + alpha = phi_n_1 + I_term = 'I' * Ps_LCU.n_qubits + self.R_LCU = PauliwordOp.from_dictionary({I_term: np.cos(alpha / 2)}) - sin_term = -np.sin(alpha / 2) + sin_term = -np.sin(alpha / 2) - for dkPk in no_βsPs: - dk_PkPs = dkPk * Ps_LCU - self.R_LCU += dk_PkPs.multiply_by_constant(sin_term) + for dkPk in no_βsPs: + dk_PkPs = dkPk * Ps_LCU + self.R_LCU += dk_PkPs.multiply_by_constant(sin_term) return Ps_LCU -def LCU_as_seq_rot(AC_op: PauliwordOp, include_global_phase_correction=False): +def LCU_as_seq_rot(R_LCU: PauliwordOp): """ Convert a unitary composed of a See equations 18 and 19 of https://arxiv.org/pdf/1907.09040.pdf - Note global phase fix is not necessary Args: - AC_op (PauliwordOp): unitary composed as a linear combination of anticommuting Pauli operators (excluding identity) - include_global_phase_correction (optional): whether to fix global phase to matrix input operator matrix exactly + R_LCU (PauliwordOp): unitary composed as a normalized linear combination of imaginary anticommuting Pauli operators (excluding identity) Returns: expon_p_terms (list): list of rotations generated by Pauli operators to implement AC_op unitary @@ -326,40 +304,34 @@ def LCU_as_seq_rot(AC_op: PauliwordOp, include_global_phase_correction=False): check = reduce(lambda a,b: a*b, [exponentiate_single_Pop(x.multiply_by_constant(1j*y/2)) for x, y in exp_terms]) print(check == rotations_LCU) """ - - assert AC_op.n_terms > 1, 'AC_op must have more than 1 term' - assert np.isclose(np.linalg.norm(AC_op.coeff_vec), 1), 'AC_op must be l2 normalized' + if isinstance(R_LCU, list) and len(R_LCU)==0: + # case where there are no rotations + return list() + + assert R_LCU.n_terms > 1, 'AC_op must have more than 1 term' + assert np.isclose(np.linalg.norm(R_LCU.coeff_vec), 1), 'AC_op must be l2 normalized' expon_p_terms = [] # IF imaginary components the this makes real (but need phase correction later!) - coeff_vec = AC_op.coeff_vec.real + AC_op.coeff_vec.imag + coeff_vec = R_LCU.coeff_vec.real + R_LCU.coeff_vec.imag for k, c_k in enumerate(coeff_vec): - P_k = AC_op[k] + P_k = R_LCU[k] theta_k = np.arcsin(c_k / np.linalg.norm(coeff_vec[:(k + 1)])) P_k.coeff_vec[0] = 1 expon_p_terms.append(tuple((P_k, theta_k))) expon_p_terms = [*expon_p_terms, *expon_p_terms[::-1]] - ### check - # from symmer.evolution.exponentiation import exponentiate_single_Pop - # terms = [exponentiate_single_Pop(op.multiply_by_constant(1j*angle/2)) for op,angle in expon_p_terms] - # final_op = reduce(lambda x,y: x*y, terms) * PauliwordOp.from_dictionary({'I'*AC_op.n_qubits: -1j}) - # assert AC_op == final_op - - # in circuit this would be done with Z * Y * X gate series: - # global phase correction (not necessary) - ## phase_correction = PauliwordOp.from_dictionary({'I'*AC_op.n_qubits: -1j}) - - if include_global_phase_correction: - ## multiply by -1j Identity term! - phase_rot = (PauliwordOp.from_dictionary({'I' * AC_op.n_qubits: 1}), -np.pi) - expon_p_terms.append(phase_rot) - - # check1 = reduce(lambda a,b: a*b, [exponentiate_single_Pop(x.multiply_by_constant(1j*y/2)) for x, y in expon_p_terms]) - # assert check1 == AC_op - + ##### manual phase correction with a rotation! + # if include_global_phase_correction: + # ## multiply by -1j Identity term! + # phase_rot = (PauliwordOp.from_dictionary({'I' * R_LCU.n_qubits: 1}), -np.pi) + # expon_p_terms.append(phase_rot) + + ## phase correction - change angle by -pi in first rotation! + expon_p_terms[0] = (expon_p_terms[0][0], expon_p_terms[0][1]-np.pi) + return expon_p_terms # from symmer.operators.utils import mul_symplectic diff --git a/symmer/operators/base.py b/symmer/operators/base.py index 7c9e7f2a..a5c643c1 100644 --- a/symmer/operators/base.py +++ b/symmer/operators/base.py @@ -945,7 +945,7 @@ def __getitem__(self, elif isinstance(key, (list, np.ndarray)): mask = np.asarray(key) else: - raise ValueError('Unrecognised input, must be an integer, slice, list or np.array') + raise ValueError(f'Unrecognised input {type(key)}, must be an integer, slice, list or np.array') symp_items = self.symp_matrix[mask] coeff_items = self.coeff_vec[mask] @@ -1112,7 +1112,8 @@ def is_noncontextual(self) -> bool: def _rotate_by_single_Pword(self, Pword: "PauliwordOp", - angle: float = None + angle: float = None, + threshold: float = 1e-18 ) -> "PauliwordOp": """ Let R(t) = e^{i t/2 Q} = cos(t/2)*I + i*sin(t/2)*Q, then one of the following can occur: @@ -1128,10 +1129,17 @@ def _rotate_by_single_Pword(self, Args: Pword (PauliwordOp): The Pauliword to rotate by. angle (float): The rotation angle in radians. If None, a Clifford rotation (angle=pi/2) is assumed. - + threshold (float): Angle threshold for Clifford rotation (precision to which the angle is a multiple of pi/2) Returns: PauliwordOp: The rotated operator. """ + if angle is None: # for legacy compatibility + angle = np.pi/2 + + if angle.imag != 0: + warnings.warn('Complex component in angle: this will be ignored.') + angle = angle.real + assert(Pword.n_terms==1), 'Only rotation by single Pauliword allowed here' if Pword.coeff_vec[0] != 1: # non-1 coefficients will affect the sign and angle in the exponent of R(t) @@ -1151,15 +1159,25 @@ def _rotate_by_single_Pword(self, commute_self = PauliwordOp(self.symp_matrix[commute_vec], self.coeff_vec[commute_vec]) anticom_self = PauliwordOp(self.symp_matrix[~commute_vec], self.coeff_vec[~commute_vec]) - if angle is None: - # assumes pi/2 rotation so Clifford - anticom_part = (anticom_self*Pword_copy).multiply_by_constant(-1j) + multiple = angle * 2 / np.pi + int_part = round(multiple) + if abs(int_part - multiple)<=threshold: + if int_part % 2 == 0: + # no rotation for angle congruent to 0 or pi disregarding sign, fixed in next line + anticom_part = anticom_self + else: + # rotation of -pi/2 disregarding sign, fixed in next line + anticom_part = (anticom_self*Pword_copy).multiply_by_constant(-1j) + if int_part in [2,3]: + anticom_part = anticom_part.multiply_by_constant(-1) # if rotation is Clifford cannot produce duplicate terms so cleanup not necessary return PauliwordOp( np.vstack([anticom_part.symp_matrix, commute_self.symp_matrix]), np.hstack([anticom_part.coeff_vec, commute_self.coeff_vec]) ) else: + if abs(angle)>1e6: + warnings.warn('Large angle can lead to precision errors: recommend using high-precision math library such as mpmath or redefine angle in range [-pi, pi]') # if angle is specified, performs non-Clifford rotation anticom_part = (anticom_self.multiply_by_constant(np.cos(angle)) + (anticom_self*Pword_copy).multiply_by_constant(-1j*np.sin(angle))) @@ -1183,9 +1201,12 @@ def perform_rotations(self, PauliwordOp: The operator after performing the rotations. """ op_copy = self.copy() - for pauli_rotation, angle in rotations: - op_copy = op_copy._rotate_by_single_Pword(pauli_rotation, angle).cleanup() - return op_copy + if rotations == []: + return op_copy.cleanup() + else: + for pauli_rotation, angle in rotations: + op_copy = op_copy._rotate_by_single_Pword(pauli_rotation, angle).cleanup() + return op_copy def tensor(self, right_op: "PauliwordOp") -> "PauliwordOp": """ diff --git a/symmer/operators/noncontextual_op.py b/symmer/operators/noncontextual_op.py index 20eb3870..fbf23b22 100644 --- a/symmer/operators/noncontextual_op.py +++ b/symmer/operators/noncontextual_op.py @@ -597,6 +597,7 @@ def noncon_state(self, UP_method:Optional[str]= 'LCU') -> Tuple[QuantumState, np 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') else: @@ -624,13 +625,13 @@ def noncon_state(self, UP_method:Optional[str]= 'LCU') -> Tuple[QuantumState, np ## undo clifford rotations from symmer.evolution.exponentiation import exponentiate_single_Pop - for op, _ in clifford_rots: + 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 = rotations_LCU.dagger * state + 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 diff --git a/symmer/projection/contextual_subspace.py b/symmer/projection/contextual_subspace.py index d9f8fdd5..30ea3a06 100644 --- a/symmer/projection/contextual_subspace.py +++ b/symmer/projection/contextual_subspace.py @@ -319,16 +319,7 @@ def project_onto_subspace(self, operator_to_project:PauliwordOp=None) -> Pauliwo self.S3_initialized = True # perform unitary partitioning if self.perform_unitary_partitioning: - # the rotation is implemented differently depending on the choice of LCU or seq_rot - if self.noncontextual_operator.up_method=='LCU': - # linear-combination-of-unitaries approach - rotated_op = (self.noncontextual_operator.unitary_partitioning_rotations * operator_to_project - * self.noncontextual_operator.unitary_partitioning_rotations.dagger).cleanup() - elif self.noncontextual_operator.up_method=='seq_rot': - # sequence-of-rotations approach - rotated_op = operator_to_project.perform_rotations(self.noncontextual_operator.unitary_partitioning_rotations) - else: - raise ValueError('Unrecognised unitary partitioning rotation method, must be one of LCU or seq_rot.') + rotated_op = operator_to_project.perform_rotations(self.noncontextual_operator.unitary_partitioning_rotations) else: rotated_op = operator_to_project # finally, project the operator before returning @@ -365,10 +356,9 @@ def project_state_onto_subspace(self, state_to_project = self.ref_state if self.perform_unitary_partitioning: - # behaviour is different whether using the LCU or seq_rot UP methods - if self.noncontextual_operator.up_method == 'LCU': - rotation = self.noncontextual_operator.unitary_partitioning_rotations - elif self.noncontextual_operator.up_method == 'seq_rot': + if self.noncontextual_operator.unitary_partitioning_rotations == []: + rotation = PauliwordOp.from_list(['I'*self.operator.n_qubits]) + else: rotation_generator = sum([R*angle*.5*1j for R,angle in self.noncontextual_operator.unitary_partitioning_rotations]) rotation = trotter(rotation_generator) return self.project_state(rotation * state_to_project) diff --git a/tests/test_operators/test_anticommuting_op.py b/tests/test_operators/test_anticommuting_op.py index 09ecee7f..15e018c5 100644 --- a/tests/test_operators/test_anticommuting_op.py +++ b/tests/test_operators/test_anticommuting_op.py @@ -62,7 +62,7 @@ def test_unitary_partitioning_no_s_index_seq_rot(): Ps, rotations, gamma_l, AC_op = AcOp_real.unitary_partitioning(s_index=None, up_method='seq_rot') assert Ps.n_terms==1, 'can only rotate onto single Pauli operator' - assert gamma_l == np.linalg.norm(list(anti_commuting_real.values())), 'normalization wrong' + assert np.isclose(gamma_l, np.linalg.norm(list(anti_commuting_real.values()))), 'normalization wrong' assert len(rotations) == AcOp_real.n_terms-1, 'seq of rotation number incorrect' P_red = AcOp_real.perform_rotations(rotations) @@ -91,15 +91,10 @@ def test_unitary_partitioning_no_s_index_LCU(): Ps, rotations, gamma_l, AC_op = AcOp_real.unitary_partitioning(s_index=None, up_method='LCU') assert Ps.n_terms==1, 'can only rotate onto single Pauli operator' - assert gamma_l == np.linalg.norm(list(anti_commuting_real.values())), 'normalization wrong' + assert np.isclose(gamma_l, np.linalg.norm(list(anti_commuting_real.values()))), 'normalization wrong' - assert isinstance(rotations, PauliwordOp) - assert rotations.n_terms == AcOp_real.n_terms, 'LCU op contains too many pauli operators' - - R_AC_op_Rdag = rotations * AcOp_real * rotations.dagger + R_AC_op_Rdag = AcOp_real.perform_rotations(rotations) assert R_AC_op_Rdag == Ps.multiply_by_constant(gamma_l) - assert rotations*rotations.dagger == PauliwordOp.from_dictionary({'I'*AcOp_real.n_qubits : 1}), 'R not unitary' - def test_unitary_partitioning_s_index_seq_rot(): AcOp_real = AntiCommutingOp.from_dictionary(anti_commuting_real) @@ -109,7 +104,7 @@ def test_unitary_partitioning_s_index_seq_rot(): up_method='seq_rot') assert Ps.n_terms==1, 'can only rotate onto single Pauli operator' - assert gamma_l == np.linalg.norm(list(anti_commuting_real.values())), 'normalization wrong' + assert np.isclose(gamma_l, np.linalg.norm(list(anti_commuting_real.values()))), 'normalization wrong' assert len(rotations) == AcOp_real.n_terms-1, 'seq of rotation number incorrect' P_red = AcOp_real.perform_rotations(rotations) @@ -141,15 +136,10 @@ def test_unitary_partitioning_s_index_LCU(): up_method='LCU') assert Ps.n_terms==1, 'can only rotate onto single Pauli operator' - assert gamma_l == np.linalg.norm(list(anti_commuting_real.values())), 'normalization wrong' - - assert isinstance(rotations, PauliwordOp) - assert rotations.n_terms == AcOp_real.n_terms, 'LCU op contains too many pauli operators' + assert np.isclose(gamma_l, np.linalg.norm(list(anti_commuting_real.values()))), 'normalization wrong' - R_AC_op_Rdag = rotations * AcOp_real * rotations.dagger + R_AC_op_Rdag = AcOp_real.perform_rotations(rotations) assert R_AC_op_Rdag == Ps.multiply_by_constant(gamma_l) - assert rotations*rotations.dagger == PauliwordOp.from_dictionary({'I'*AcOp_real.n_qubits : 1}), 'R not unitary' - def test_unitary_partitioning_seq_rot_complex(): AcOp_comp = AntiCommutingOp.from_dictionary(anti_commuting_complex) @@ -206,7 +196,7 @@ def test_ac_set_with_zero_ceoffs(): Ps_LCU, rotations_LCU, gamma_l, AC_normed = AcOp_real.unitary_partitioning(s_index=1, up_method='LCU') - LCU_output = (rotations_LCU * AC_normed * rotations_LCU.dagger).cleanup() + LCU_output = AC_normed.perform_rotations(rotations_LCU).cleanup() assert LCU_output.n_terms==1 assert np.isclose(LCU_output.coeff_vec[0], 1) assert Ps_LCU == LCU_output @@ -218,13 +208,13 @@ def test_ac_set_with_negative_and_zero_ceoffs(): up_method='seq_rot') seq_rot_output = AC_normed.perform_rotations(rotations_seq_rot) assert seq_rot_output.n_terms==1 - assert np.isclose(seq_rot_output.coeff_vec[0], 1) + assert np.isclose(seq_rot_output.coeff_vec[0], -1) assert Ps_seq_rot == seq_rot_output Ps_LCU, rotations_LCU, gamma_l, AC_normed = AcOp_real.unitary_partitioning(s_index=0, up_method='LCU') - LCU_output = (rotations_LCU * AC_normed * rotations_LCU.dagger).cleanup() + LCU_output = AC_normed.perform_rotations(rotations_LCU) assert LCU_output.n_terms==1 - assert np.isclose(LCU_output.coeff_vec[0], 1) + assert np.isclose(LCU_output.coeff_vec[0], -1) assert Ps_LCU == LCU_output \ No newline at end of file