Fix ryy gate definition

pyquil/paulis.py

@@ -36,6 +36,9 @@
  Union,
  cast,
 )
+from numpy.typing import NDArray
+from functools import reduce
+from scipy.linalg import expm
 
 from pyquil.quilatom import (
  QubitPlaceholder,
@@ -840,7 +843,6 @@ def commuting_sets(pauli_terms: PauliSum) -> List[List[PauliTerm]]:
  isAssigned_bool = False
  for p in range(m_s): # check if it commutes with each group
  if isAssigned_bool is False:
-
  if check_commutation(groups[p], pauli_terms.terms[j]):
  isAssigned_bool = True
  groups[p].append(pauli_terms.terms[j])
@@ -931,6 +933,64 @@ def combined_exp_wrap(param: Union[float, MemoryReference]) -> Program:
  return combined_exp_wrap
 
 
+def exponentiate_pauli_sum(
+ pauli_sum: Union[PauliSum, PauliTerm],
+) -> NDArray[np.complex_]:
+ r"""
+ Exponentiates a sequence of PauliTerms, which may or may not commute. The Pauliterms must
+ have fixed (non-parametric) coefficients. The coefficients are interpreted in cycles
+ rather than radians or degrees.
+
+ e^{-i\pi\sum_i \theta_i P_i}
+
+ Thus, 1.0 corresponds to one full rotation.
+
+ To produce a CZ gate:
+
+ >>> phi = pi
+ >>> coeff = phi/(-4*pi) # -0.25
+ >>> hamiltonian = PauliTerm("Z", 0) * PauliTerm("Z", 1) - 1*PauliTerm("Z", 0) - 1*PauliTerm("Z", 1)
+ >>> exponentiate_pauli_sum(coeff*hamiltonian)
+
+ To produce the Quil XY(theta) gate, you can use:
+
+ >>> coeff = theta/(-4*pi)
+ >>> hamiltonian = PauliTerm("X", 0) * PauliTerm("X", 1) + PauliTerm("Y", 0) * PauliTerm("Y", 1)
+ >>> exponentiate_pauli_sum(coeff*hamiltonian)
+
+ A global phase is applied to the unitary such that the [0,0] element is always real.
+
+ :param pauli_sum: PauliSum to exponentiate.
+ :returns: The matrix exponetial of the PauliSum
+ """
+ if isinstance(pauli_sum, PauliTerm):
+ pauli_sum = PauliSum([pauli_sum])
+
+ pauli_matrices = {
+ "X": np.array([[0.0, 1.0], [1.0, 0.0]]),
+ "Y": np.array([[0.0, 0.0 - 1.0j], [0.0 + 1.0j, 0.0]]),
+ "Z": np.array([[1.0, 0.0], [0.0, -1.0]]),
+ "I": np.array([[1.0, 0.0], [0.0, 1.0]]),
+ }
+
+ qubits = pauli_sum.get_qubits()
+ for q in qubits:
+ assert isinstance(q, int)
+ qubits.sort()
+
+ matrices = []
+ for term in pauli_sum.terms:
+ coeff = term.coefficient
+ assert isinstance(coeff, Number)
+ qubit_paulis = {qubit: pauli for qubit, pauli in term.operations_as_set()}
+ paulis = [qubit_paulis[q] if q in qubit_paulis else "I" for q in qubits]
+ matrix = float(np.real(coeff)) * reduce(np.kron, [pauli_matrices[p] for p in paulis]) # type: ignore
+ matrices.append(matrix)
+ generated_unitary = expm(-1j * np.pi * sum(matrices))
+ phase = np.exp(-1j * np.angle(generated_unitary[0, 0])) # type: ignore
+ return np.asarray(phase * generated_unitary, dtype=np.complex_)
+
+
 def _exponentiate_general_case(pauli_term: PauliTerm, param: float) -> Program:
  """
  Returns a Quil (Program()) object corresponding to the exponential of

pyquil/simulation/matrices.py

@@ -94,6 +94,16 @@
  PHASEDFSIM(:math:`\theta, \zeta, \chi, \gamma, \phi`) - XX+YY interaction with conditonal phase on \|11\>
  :math:`\begin{pmatrix} 1&0&0&0 \\ 0&\ e^{-i(\gamma+\zeta)}\cos(\frac{\theta}{2})&ie^{-i(\gamma-\chi)}\sin(\frac{\theta}{2})&0 \\ 0&ie^{-i(\gamma+\chi)}\sin(\frac{\theta}{2})&e^{-i(\gamma-\zeta)}\cos(\frac{\theta}{2})&0 \\ 0&0&0&e^{ i\phi - 2i\gamma} \end{pmatrix}`
 
+ RXX(:math:`\phi`) - XX-interaction
+ :math:`\begin{pmatrix} \cos(\phi/2)&0&0&-i\sin(\phi/2) \\ 0&\cos(\phi/2)&-i\sin(\phi/2)&0 \\ 0&-i\sin(\phi/2)&\cos(\phi/2)&0 \\ -i\sin(\phi/2)&0&0&cos(\phi/2) \end{pmatrix}`
+
+ RYY(:math:`\phi`) - YY-interaction
+ :math:`\begin{pmatrix} \cos(\phi/2)&0&0&i\sin(\phi/2) \\ 0&\cos(\phi/2)&-i\sin(\phi/2)&0 \\ 0&-i\sin(\phi/2)&\cos(\phi/2)&0 \\ i\sin(\phi/2)&0&0&cos(\phi/2) \end{pmatrix}`
+
+ RZZ(:math:`\phi`) - ZZ-interaction
+ :math:`\begin{pmatrix} 1&0&0&0 \\ 0&\cos(\phi/2)&-i\sin(\phi/2)&0 \\ 0&-i\sin(\phi/2)&\cos(\phi/2)&0 \\ 0&0&0&1 \end{pmatrix}`
+
+
 Specialized gates / internal utility gates:
  BARENCO(:math:`\alpha, \phi, \theta`) - Barenco gate
  :math:`\begin{pmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&e^{i\phi} \cos\theta & -i e^{i(\alpha-\phi)} \sin\theta \\ 0&0&-i e^{i(\alpha+\phi)} \sin\theta & e^{i\alpha} \cos\theta \end{pmatrix}`
@@ -278,10 +288,10 @@ def RXX(phi: float) -> np.ndarray:
 def RYY(phi: float) -> np.ndarray:
  return np.array(
  [
- [np.cos(phi / 2), 0, 0, -1j * np.sin(phi / 2)],
- [0, np.cos(phi / 2), +1j * np.sin(phi / 2), 0],
- [0, +1j * np.sin(phi / 2), np.cos(phi / 2), 0],
- [-1j * np.sin(phi / 2), 0, 0, np.cos(phi / 2)],
+ [np.cos(phi / 2), 0, 0, +1j * np.sin(phi / 2)],
+ [0, np.cos(phi / 2), -1j * np.sin(phi / 2), 0],
+ [0, -1j * np.sin(phi / 2), np.cos(phi / 2), 0],
+ [+1j * np.sin(phi / 2), 0, 0, np.cos(phi / 2)],
  ]
  )
 

pyquil/simulation/tools.py

@@ -445,3 +445,10 @@ def scale_out_phase(unitary1: np.ndarray, unitary2: np.ndarray) -> np.ndarray:
  rescale_value = unitary2[j, 0] / unitary1[j, 0]
 
  return rescale_value * unitary1
+
+
+def unitary_equal(A: np.ndarray, B: np.ndarray) -> bool:
+ """Check if two matrices are unitarily equal."""
+ assert A.shape == B.shape
+ dim = A.shape[0]
+ return np.allclose(np.abs(np.trace(A.T.conjugate() @ B) / dim), 1.0)

test/unit/test_paulis.py

@@ -30,6 +30,7 @@
  PauliSum,
  exponential_map,
  exponentiate_commuting_pauli_sum,
+ exponentiate_pauli_sum,
  ID,
  exponentiate,
  trotterize,
@@ -45,7 +46,8 @@
  is_identity,
 )
 from pyquil.quil import Program
-from pyquil.simulation.tools import program_unitary
+from pyquil.simulation.tools import program_unitary, unitary_equal
+from pyquil.simulation import matrices
 
 
 def isclose(a, b, rel_tol=1e-10, abs_tol=0.0):
@@ -397,6 +399,56 @@ def test_exponentiate_prog():
  assert prog == result_prog
 
 
+def test_exponentiate_pauli_sum_rxx():
+ """Test that exponentiate_exponentiate_pauli_sum generates the RXX gate"""
+ generators = PauliTerm("X", 0) * PauliTerm("X", 1)
+ for angle in np.linspace(-0.5, 0.5):
+ generated_unitary = exponentiate_pauli_sum(generators * angle)
+ assert unitary_equal(generated_unitary, matrices.RXX(2 * np.pi * angle))
+
+
+def test_exponentiate_pauli_sum_ryy():
+ """Test that exponentiate_exponentiate_pauli_sum generates the RYY gate"""
+ generators = PauliTerm("Y", 0) * PauliTerm("Y", 1)
+ for angle in np.linspace(-0.5, 0.5):
+ generated_unitary = exponentiate_pauli_sum(generators * angle)
+ assert unitary_equal(generated_unitary, matrices.RYY(2 * np.pi * angle))
+
+
+def test_exponentiate_pauli_sum_rzz():
+ """Test that exponentiate_exponentiate_pauli_sum generates the RZZ gate"""
+ generators = PauliTerm("Z", 0) * PauliTerm("Z", 1)
+ for angle in np.linspace(-0.5, 0.5):
+ generated_unitary = exponentiate_pauli_sum(generators * angle)
+ assert unitary_equal(generated_unitary, matrices.RZZ(2 * np.pi * angle))
+
+
+def test_exponentiate_pauli_sum_cphase():
+ """Test that exponentiate_exponentiate_pauli_sum generates the CZ gate"""
+ generators = PauliTerm("Z", 0) * PauliTerm("Z", 1) - 1 * PauliTerm("Z", 0) - 1 * PauliTerm("Z", 1)
+ for angle in np.linspace(-0.5, 0.5):
+ generated_unitary = exponentiate_pauli_sum(generators * angle)
+ assert unitary_equal(generated_unitary, matrices.CPHASE(2 * np.pi * (-2 * angle)))
+
+
+def test_exponentiate_pauli_sum_xy():
+ """Test that exponentiate_exponentiate_pauli_sum generates the XY gate"""
+ generators = PauliTerm("X", 0) * PauliTerm("X", 1) + PauliTerm("Y", 0) * PauliTerm("Y", 1)
+ for angle in np.linspace(-0.5, 0.5):
+ generated_unitary = exponentiate_pauli_sum(generators * angle)
+ assert unitary_equal(generated_unitary, matrices.XY(2 * np.pi * (-2 * angle)))
+
+
+def test_exponentiate_pauli_sum_fsim():
+ """Test that exponentiate_exponentiate_pauli_sum generates the FSIM gate"""
+ xy_generators = PauliTerm("X", 0) * PauliTerm("X", 1) + PauliTerm("Y", 0) * PauliTerm("Y", 1)
+ cphase_generators = PauliTerm("Z", 0) * PauliTerm("Z", 1) - 1 * PauliTerm("Z", 0) - 1 * PauliTerm("Z", 1)
+ for theta in np.linspace(-0.5, 0.5):
+ for phi in np.linspace(-0.5, 0.5):
+ generated_unitary = exponentiate_pauli_sum(xy_generators * theta + cphase_generators * phi)
+ assert unitary_equal(generated_unitary, matrices.FSIM(2 * np.pi * (-2 * theta), 2 * np.pi * (-2 * phi)))
+
+
 def test_exponentiate_identity():
  generator = PauliTerm("I", 1, 0.0)
  para_prog = exponential_map(generator)
@@ -566,7 +618,6 @@ def test_check_commutation_rigorous():
  commuting_pairs = []
  for x in range(len(pauli_ops_pq)):
  for y in range(x, len(pauli_ops_pq)):
-
  tmp_op = _commutator(pauli_ops_pq[x], pauli_ops_pq[y])
  assert len(tmp_op.terms) == 1
  if is_zero(tmp_op.terms[0]):
@@ -608,10 +659,10 @@ def test_term_powers():
  for qubit_id in range(2):
  pauli_terms = [sI(qubit_id), sX(qubit_id), sY(qubit_id), sZ(qubit_id)]
  for pauli_term in pauli_terms:
- assert pauli_term ** 0 == sI(qubit_id)
- assert pauli_term ** 1 == pauli_term
- assert pauli_term ** 2 == sI(qubit_id)
- assert pauli_term ** 3 == pauli_term
+ assert pauli_term**0 == sI(qubit_id)
+ assert pauli_term**1 == pauli_term
+ assert pauli_term**2 == sI(qubit_id)
+ assert pauli_term**3 == pauli_term
  with pytest.raises(ValueError):
  pauli_terms[0] ** -1
  # Test to make sure large powers can be computed
@@ -620,13 +671,13 @@ def test_term_powers():
 
 def test_sum_power():
  pauli_sum = (sY(0) - sX(0)) * (1.0 / np.sqrt(2))
- assert pauli_sum ** 2 == PauliSum([sI(0)])
+ assert pauli_sum**2 == PauliSum([sI(0)])
  with pytest.raises(ValueError):
- _ = pauli_sum ** -1
+ _ = pauli_sum**-1
  pauli_sum = sI(0) + sI(1)
- assert pauli_sum ** 0 == sI(0)
+ assert pauli_sum**0 == sI(0)
  # Test to make sure large powers can be computed
- pauli_sum ** 400
+ pauli_sum**400
 
 
 def test_term_equality():