diff --git a/pyquil/paulis.py b/pyquil/paulis.py index 5325c2e8c..b43580ba7 100644 --- a/pyquil/paulis.py +++ b/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 diff --git a/pyquil/simulation/matrices.py b/pyquil/simulation/matrices.py index 1ae425e36..9b8a61180 100644 --- a/pyquil/simulation/matrices.py +++ b/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)], ] ) diff --git a/pyquil/simulation/tools.py b/pyquil/simulation/tools.py index a6b7d7294..50404c29d 100644 --- a/pyquil/simulation/tools.py +++ b/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) diff --git a/test/unit/test_paulis.py b/test/unit/test_paulis.py index 5807b1bc2..520b58806 100644 --- a/test/unit/test_paulis.py +++ b/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():