Add method for tapering gate operations

pennylane/qchem/__init__.py

@@ -66,4 +66,5 @@
     optimal_sector,
     taper,
     taper_hf,
+    taper_excitations,
 )

pennylane/qchem/tapering.py

@@ -538,3 +538,77 @@ def taper_hf(generators, paulixops, paulix_sector, num_electrons, num_wires):
         tapered_hartree_fock.append(0)
 
     return np.array(tapered_hartree_fock).astype(int)
+
+
+def _is_commuting_obs(ham_a, ham_b, wire_map=None):
+    r"""Check for commutivity between two Pauli observables.
+
+    Args:
+        ham_a (Hamiltonian): first observable
+        ham_b (Hamiltonian): second observable
+
+    Returns:
+        Bool: representing whether ham_a and ham_b commutes or not.
+    """
+    for op1 in ham_a.ops:
+        for op2 in ham_b.ops:
+            if not qml.grouping.is_commuting(op1, op2, wire_map):
+                return False
+    return True
+
+
+def taper_excitations(generators, paulixops, paulix_sector, singles, doubles):
+    r"""Transform excitations with a Clifford operator and taper qubits.
+
+    The qubit operators for single and double excitations are first generated using the generators of
+    :func:`~.SingleExcitation` and :func:`~.DoubleExcitation` operations. Each of these operators that commutes
+    with all :math:`\mathbb{Z}_2` symmetries of the molecular Hamiltonian are then tranformed using the
+    Clifford operators :math:`U` and then tapered, while rest of the other non-commuting operators are discarded.
+    These new tapered excitation operators can be exponentiated using :func:`~.PauliRot` for building a
+    tapered UCCSD-like circuit ansatze.
+
+    Args:
+        generators (list[Hamiltonian]): list of generators of symmetries, taus, for the Hamiltonian
+        paulixops (list[Operation]):  list of single-qubit Pauli-X operators
+        paulix_sector (list[int]): list of eigenvalues of Pauli-X operators
+        singles (list(list(int))): list with the indices `r`, `p` of the two qubits representing the single excitation :math:`\vert r, p \rangle = \hat{c}_p^\dagger \hat{c}_r \vert \mathrm{HF}\rangle`
+        doubles (list(list(int))): list with the indices `s`, `r`, `q`, `p` of the four qubits representing the double excitation :math:`\vert s, r, q, p \rangle = \hat{c}_p^\dagger \hat{c}_q^\dagger \hat{c}_r \hat{c}_s \vert \mathrm{HF}\rangle`
+
+    Returns:
+        tuple(list, list): tapered single and double excitation operators
+
+    **Example**
+
+    >>> symbols = ['He', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4588684632]])
+    >>> mol = qml.qchem.Molecule(symbols, geometry, charge=1)
+    >>> H, n_qubits = qml.qchem.molecular_hamiltonian(symbols, geometry)
+    >>> n_elec = mol.n_electrons
+    >>> generators = qml.qchem.symmetry_generators(H)
+    >>> paulixops = qml.qchem.paulix_ops(generators, 4)
+    >>> paulix_sector = qml.qchem.optimal_sector(H, generators, n_elec)
+    >>> singles, doubles = qml.qchem.excitations(n_elec, n_qubits)
+    >>> singles_tap, doubles_tap = taper_excitations(generators, paulixops,
+                                            paulix_sector, singles, doubles)
+    >>> print(singles_tap[0], doubles_tap[0])
+    ((0.5+0j)) [Y0]
+    ((-0.25+0j)) [X0 Y1] + ((-0.25+0j)) [Y0 X1]
+    """
+
+    singles_tapered, doubles_tapered = [], []
+
+    for excitation in singles:
+        hamil_gen = qml.SingleExcitation(1, wires=excitation).generator()
+        if np.all([_is_commuting_obs(generator, hamil_gen) for generator in generators]):
+            excitation_tapered_op = qml.taper(hamil_gen, generators, paulixops, paulix_sector)
+            qml.simplify(excitation_tapered_op)
+            singles_tapered.append(excitation_tapered_op)
+
+    for excitation in doubles:
+        hamil_gen = qml.DoubleExcitation(1, wires=excitation).generator()
+        if np.all([_is_commuting_obs(generator, hamil_gen) for generator in generators]):
+            excitation_tapered_op = qml.taper(hamil_gen, generators, paulixops, paulix_sector)
+            qml.simplify(excitation_tapered_op)
+            doubles_tapered.append(excitation_tapered_op)
+
+    return singles_tapered, doubles_tapered

tests/qchem/test_tapering.py

@@ -29,6 +29,8 @@
     clifford,
     optimal_sector,
     taper_hf,
+    _is_commuting_obs,
+    taper_excitations,
 )
 
 
@@ -629,3 +631,167 @@ def test_taper_obs(symbols, geometry, charge):
             @ scipy.sparse.coo_matrix(state_tapered).T
         ).toarray()
         assert np.isclose(obs_val, obs_val_tapered)
+
+
+@pytest.mark.parametrize(
+    ("obs_1", "obs_2", "wire_map", "commute_status"),
+    [
+        (
+            qml.Hamiltonian([1.0], [qml.Identity(0)]),
+            qml.Hamiltonian([1.0], [qml.PauliZ(0)]),
+            {0: 0},
+            True,
+        ),
+        (
+            qml.Hamiltonian([1.0, 1.0], [qml.PauliY(0), qml.PauliX(1)]),
+            qml.Hamiltonian([1.0, 1.0], [qml.PauliZ(0), qml.PauliX(1)]),
+            {0: 0, 1: 1},
+            False,
+        ),
+        (
+            qml.Hamiltonian(
+                [1.0, 1.0], [qml.PauliY(0) @ qml.PauliX(2), qml.PauliX(0) @ qml.PauliY(2)]
+            ),
+            qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+            {0: 0, 1: 1, 2: 2},
+            False,
+        ),
+        (
+            qml.Hamiltonian(
+                [1.0, 1.0], [qml.PauliY(0) @ qml.PauliX(2), qml.PauliX(0) @ qml.PauliY(2)]
+            ),
+            qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+            {0: 0, 1: 1, 2: 2},
+            True,
+        ),
+        (
+            qml.Hamiltonian(
+                [1.0, 1.0], [qml.PauliY("b") @ qml.PauliX("d"), qml.PauliX("b") @ qml.PauliY("d")]
+            ),
+            qml.Hamiltonian([1.0], [qml.PauliZ("b") @ qml.PauliZ("c")]),
+            {"a": 0, "b": 1, "c": 2, "d": 3},
+            False,
+        ),
+        (
+            qml.Hamiltonian(
+                [1.0, 1.0], [qml.PauliY("b") @ qml.PauliX("d"), qml.PauliX("b") @ qml.PauliY("d")]
+            ),
+            qml.Hamiltonian([1.0], [qml.PauliZ("b") @ qml.PauliZ("d")]),
+            {"a": 0, "b": 1, "c": 2, "d": 3},
+            True,
+        ),
+    ],
+)
+def test_is_commuting_obs(obs_1, obs_2, wire_map, commute_status):
+    """Test that (non)-commuting Pauli observables are correctly identified."""
+    do_they_commute = _is_commuting_obs(obs_1, obs_2, wire_map=wire_map)
+    assert do_they_commute == commute_status
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "charge", "num_commuting"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.40104295]], requires_grad=True),
+            0,
+            (0, 1),
+        ),
+        (
+            ["He", "H"],
+            np.array(
+                [[0.0, 0.0, 0.0], [0.0, 0.0, 1.4588684632]],
+                requires_grad=True,
+            ),
+            1,
+            (2, 1),
+        ),
+        (
+            ["H", "H", "H"],
+            np.array(
+                [[-0.84586466, 0.0, 0.0], [0.84586466, 0.0, 0.0], [0.0, 1.46508057, 0.0]],
+                requires_grad=True,
+            ),
+            1,
+            (4, 4),
+        ),
+        (
+            ["H", "H", "H", "H"],
+            np.array(
+                [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0], [0.0, 0.0, 3.0]],
+                requires_grad=True,
+            ),
+            0,
+            (4, 10),
+        ),
+    ],
+)
+def test_taper_excitations(symbols, geometry, charge, num_commuting):
+    r"""Test that the tapered excitation operators are consistent with the
+    tapered Hartree-Fock state."""
+    mol = qml.qchem.Molecule(symbols, geometry, charge)
+    hamiltonian = qml.qchem.diff_hamiltonian(mol)(geometry)
+    hf_state = np.where(np.arange(len(hamiltonian.wires)) < mol.n_electrons, 1, 0)
+    generators = qml.symmetry_generators(hamiltonian)
+    paulixops = qml.paulix_ops(generators, len(hamiltonian.wires))
+    paulix_sector = optimal_sector(hamiltonian, generators, mol.n_electrons)
+    hf_state_tapered = taper_hf(
+        generators, paulixops, paulix_sector, mol.n_electrons, len(hamiltonian.wires)
+    )
+    particle_num = qml.qchem.particle_number(len(hamiltonian.wires))
+    particle_num_tapered = qml.taper(particle_num, generators, paulixops, paulix_sector)
+
+    o = np.array([1, 0])
+    l = np.array([0, 1])
+    state = functools.reduce(lambda i, j: np.kron(i, j), [l if s else o for s in hf_state])
+    state_tapered = functools.reduce(
+        lambda i, j: np.kron(i, j), [l if s else o for s in hf_state_tapered]
+    )
+    singles, doubles = qml.qchem.excitations(mol.n_electrons, 2 * mol.n_orbitals)
+    singles_obs = [qml.SingleExcitation(1, wires=single).generator() for single in singles]
+    doubles_obs = [qml.DoubleExcitation(1, wires=double).generator() for double in doubles]
+
+    singles_tap, doubles_tap = taper_excitations(
+        generators, paulixops, paulix_sector, singles, doubles
+    )
+    assert len(singles_tap) == num_commuting[0] and len(doubles_tap) == num_commuting[1]
+
+    exc_iter = iter(singles_tap + doubles_tap)
+    for observable in singles_obs + doubles_obs:
+        if np.all([_is_commuting_obs(generator, observable) for generator in generators]):
+            excitation_obs = next(exc_iter)
+            # sanity check for consistent expectation values
+            obs_val = (
+                scipy.sparse.coo_matrix(state)
+                @ qml.utils.sparse_hamiltonian(observable, wires=range(len(hf_state)))
+                @ scipy.sparse.coo_matrix(state).T
+            ).toarray()
+            obs_val_tapered = (
+                scipy.sparse.coo_matrix(state_tapered)
+                @ qml.utils.sparse_hamiltonian(excitation_obs, wires=range(len(hf_state_tapered)))
+                @ scipy.sparse.coo_matrix(state_tapered).T
+            ).toarray()
+            assert np.isclose(obs_val, obs_val_tapered)
+
+            # check if tapered excitation are particle number conserving
+            excited_state = (
+                scipy.linalg.expm(1j * qml.matrix(observable, wire_order=range(len(hf_state))))
+                @ state
+            )
+            excited_state_tapered = (
+                scipy.linalg.expm(
+                    1j * qml.matrix(excitation_obs, wire_order=range(len(hf_state_tapered)))
+                )
+                @ state_tapered
+            )
+            pnum_val = (
+                scipy.sparse.coo_matrix(excited_state)
+                @ qml.utils.sparse_hamiltonian(particle_num)
+                @ scipy.sparse.coo_matrix(excited_state).T
+            ).toarray()
+            pnum_val_tapered = (
+                scipy.sparse.coo_matrix(excited_state_tapered)
+                @ qml.utils.sparse_hamiltonian(particle_num_tapered)
+                @ scipy.sparse.coo_matrix(excited_state_tapered).T
+            ).toarray()
+            assert np.isclose(pnum_val, pnum_val_tapered)