Merge pull request #176 from UCL-CCS/v1_upgrade

V1 upgrade

pyproject.toml

@@ -17,7 +17,6 @@ flake8 = "^4.0.1"
 py3Dmol = "^1.8.0"
 ncon = "^1.0.0"
 opt-einsum = "^3.3.0"
-qubovert = "^1.2.5"
 pytest = "^7.2.2"
 numba = "0.56"
 quimb = "^1.5.1"

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,31 @@ 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
+
+ @classmethod
+ def random(cls, n_qubits: int, n_terms: Union[None, int]=None, apply_clifford=True) -> "AntiCommutingOp":
+ """
+ generate a random real coefficient anticommuting op
+
+ """
+ from symmer.utils import random_anitcomm_2n_1_PauliwordOp
+ if n_terms is None:
+ n_terms = 2*n_qubits+1
+
+ assert n_terms<= 2*n_qubits+1, f'cannot have {n_terms} Pops on {n_qubits} qubits'
+ return cls.from_PauliwordOp( random_anitcomm_2n_1_PauliwordOp(n_qubits, apply_clifford=apply_clifford)[:n_terms])
 
  def generate_LCU_operator(self, AC_op) -> PauliwordOp:
  """
@@ -235,74 +255,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
+ ## 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]
 
- if (before_cleanup>1 and post_cleanup==1):
- if AC_op.coeff_vec[0]<0:
- # need to fix neg sign (use Pauli multiplication)
+ # ∑ β_k 𝑃_k ... note this doesn't contain 𝛽_s 𝑃_s
+ no_βsPs = AC_op - (Ps_LCU.multiply_by_constant(βs))
 
- # 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))
+ # Ω_𝑙 ∑ 𝛿_k 𝑃_k ... renormalized!
+ omega_l = np.linalg.norm(no_βsPs.coeff_vec)
+ no_βsPs.coeff_vec = no_βsPs.coeff_vec / omega_l
 
- 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]]
+ 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)})
 
-  sign_correction = PauliwordOp(AC_op_cpy.symp_matrix[1],[1])
+ sin_term = -np.sin(alpha / 2)
 
- 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
-
- 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)})
-
- 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 +317,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

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]
@@ -1098,7 +1098,8 @@ def adjacency_matrix_qwc(self) -> np.array:
  def is_noncontextual(self) -> bool:
  """ 
  Returns True if the operator is noncontextual, False if contextual
- Scales as O(N^2), compared with the O(N^3) algorithm of https://doi.org/10.1103/PhysRevLett.123.200501
+ Scales as O(M^2), compared with the O(M^3) algorithm of https://doi.org/10.1103/PhysRevLett.123.200501
+ where M is the number of terms in the operator.
  Constructing the adjacency matrix is by far the most expensive part - very fast once that has been built.
 
  Returns:
@@ -1111,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:
@@ -1127,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)
@@ -1150,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)))
@@ -1182,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":
  """