Add functions to find correct sector for qubit tapering

doc/releases/changelog-dev.md

@@ -39,7 +39,8 @@
 
 * Functions for tapering qubits based on molecular symmetries is added.
   [(#1966)](https://github.com/PennyLaneAI/pennylane/pull/1966)
-* [(#1974)](https://github.com/PennyLaneAI/pennylane/pull/1974)
+  [(#1974)](https://github.com/PennyLaneAI/pennylane/pull/1974)
+  [(#2041)](https://github.com/PennyLaneAI/pennylane/pull/2041)
 
   With this functionality, a molecular Hamiltonian can be transformed to a new Hamiltonian that acts
   on a reduced number of qubits.

pennylane/hf/__init__.py

@@ -50,4 +50,5 @@
     generate_symmetries,
     get_generators,
     transform_hamiltonian,
+    find_optimal_sector,
 )

pennylane/hf/tapering.py

@@ -15,8 +15,11 @@
 This module contains the functions needed for tapering qubits using symmetries.
 """
 # pylint: disable=unnecessary-lambda
-import functools
 
+import itertools as it
+import functools
+import scipy as sp
+import numpy
 import autograd.numpy as anp
 import pennylane as qml
 from pennylane import numpy as np
@@ -464,3 +467,106 @@ def transform_hamiltonian(h, generators, paulix_ops, paulix_sector):
     c = qml.math.stack(c)
 
     return _simplify(qml.Hamiltonian(c, o))
+
+
+def _sector_energy(sector):
+    """Computes energy of the given tapered Hamiltonian.
+
+    Args:
+        sector (pennylane.Hamiltonian): tapered Hamiltonian whose energy has to be calculated
+
+    Returns:
+        energy (float): energy of the tapered hamiltonian
+    """
+
+    if len(sector.wires) < 2:
+        energy = min(sp.linalg.eigvals(qml.utils.sparse_hamiltonian(sector).toarray()))
+    else:
+        energy = sp.sparse.linalg.eigs(qml.utils.sparse_hamiltonian(sector), k=2)[0].min()
+
+    return energy
+
+
+def find_optimal_sector(qubit_op, generators, paulix_ops, brute_force=True, num_electrons=0):
+    r"""Get the optimal sector which contains the ground state.
+
+    To obtain the optimal sector one can brute force through all the permutation or by utilize the following
+    relation between the Pauli-Z qubit operator and the occupation number under Jordan-Wigner transform.
+
+    .. math::
+
+        \sigma_z^{i} = a_{i}^{\dagger}a_{i} - 1
+
+    This relation allows us to figure out whether the orbital is occupied or unoccupied in a given symmetry sector,
+    to build the correct eigensector.
+
+    Args:
+        qubit_op (pennylane.Hamiltonian): Hamiltonian for which symmetries are being generated for performing tapering
+        generators (list[pennylane.Hamiltonian]): list of generators of symmetries, taus, for the Hamiltonian
+        paulix_ops (list[pennylane.operation.Observable]):  list of single-qubit Pauli X operators
+        brute_force (bool): determines whether to use brute-force strategy to pick the correct sector or not
+        num_electrons (int): If `brute_force = True`, user must provide the number of active electrons in
+                              the system for generating the Hartree-Fock bitstring
+
+    Returns:
+        tuple (list[int], float):
+
+            * list[int]: eigenvalues corresponding to the optimal sector which contains the ground state
+            * float: energy of the computed optimal sector
+
+    .. code-block:: python
+
+        >>> symbols = ['H', 'H']
+        >>> coordinates = np.array([0., 0., -0.66140414, 0., 0., 0.66140414]))
+        >>> mol = qml.hf.Molecule(symbols, coordinates)
+        >>> H, qubits = qml.hf.generate_hamiltonian(mol)(), 4
+        >>> generators, paulix_ops = generate_symmetries(H, qubits)
+        >>> find_optimal_sector(H, generators, paulix_ops, False, 2)
+          ([1, -1, -1], -1.1372701746609024)
+
+    """
+
+    if not brute_force:
+
+        if num_electrons < 1:
+            raise ValueError(
+                f"Brute force search is disabled;"
+                f"the number of electrons must be provided and should be greater than zero;"
+                f"got 'electrons'={num_electrons}"
+            )
+
+        num_orbitals = len(qubit_op.wires)
+
+        if num_electrons > num_orbitals:
+            raise ValueError(
+                f"Number of active orbitals cannot be smaller than number of active electrons;"
+                f" got 'orbitals'={num_orbitals} < 'electrons'={num_electrons}."
+            )
+
+        hf_str = np.where(np.arange(num_orbitals) < num_electrons, 1, 0)
+        perm = []
+
+        for op in generators:
+            symmstr = np.zeros(len(qubit_op.wires), dtype=int)
+            for wire in op.wires:
+                symmstr[wire] = 1
+            coeff = -1 if numpy.logical_xor.reduce(numpy.logical_and(symmstr, hf_str)) else 1
+            perm.append(coeff)
+
+        sector = transform_hamiltonian(qubit_op, generators, paulix_ops, perm)
+        energy = _sector_energy(sector)
+
+        return perm, energy
+
+    sectors = []
+    energies = []
+    perms = list(it.product([1, -1], repeat=len(generators)))
+
+    for perm in perms:
+        sector = transform_hamiltonian(qubit_op, generators, paulix_ops, perm)
+        energy = _sector_energy(sector)
+        sectors.append(sector)
+        energies.append(float(energy.real))
+
+    index = energies.index(min(energies))
+    return list(perms[index]), energies[index]

tests/hf/test_tapering.py

@@ -14,8 +14,8 @@
 """
 Unit tests for functions needed for qubit tapering.
 """
-import pennylane as qml
 import pytest
+import pennylane as qml
 from pennylane import numpy as np
 from pennylane.hf.tapering import (
     _binary_matrix,
@@ -28,6 +28,8 @@
     generate_symmetries,
     get_generators,
     transform_hamiltonian,
+    _sector_energy,
+    find_optimal_sector,
 )
 
 
@@ -286,59 +288,11 @@ def test_generate_paulis(generators, num_qubits, result):
 
 
 @pytest.mark.parametrize(
-    ("hamiltonian", "num_qubits", "res_generators", "res_pauli_ops"),
+    ("symbols", "geometry", "num_qubits", "res_generators", "res_pauli_ops"),
     [
         (
-            qml.Hamiltonian(
-                np.array(
-                    [
-                        -0.09886397,
-                        0.17119775,
-                        0.17119775,
-                        -0.22278593,
-                        -0.22278593,
-                        0.16862219,
-                        0.0453222,
-                        -0.0453222,
-                        -0.0453222,
-                        0.0453222,
-                        0.12054482,
-                        0.16586702,
-                        0.16586702,
-                        0.12054482,
-                        0.17434844,
-                    ]
-                ),
-                [
-                    qml.Identity(wires=[0]),
-                    qml.PauliZ(wires=[0]),
-                    qml.PauliZ(wires=[1]),
-                    qml.PauliZ(wires=[2]),
-                    qml.PauliZ(wires=[3]),
-                    qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[1]),
-                    qml.PauliY(wires=[0])
-                    @ qml.PauliX(wires=[1])
-                    @ qml.PauliX(wires=[2])
-                    @ qml.PauliY(wires=[3]),
-                    qml.PauliY(wires=[0])
-                    @ qml.PauliY(wires=[1])
-                    @ qml.PauliX(wires=[2])
-                    @ qml.PauliX(wires=[3]),
-                    qml.PauliX(wires=[0])
-                    @ qml.PauliX(wires=[1])
-                    @ qml.PauliY(wires=[2])
-                    @ qml.PauliY(wires=[3]),
-                    qml.PauliX(wires=[0])
-                    @ qml.PauliY(wires=[1])
-                    @ qml.PauliY(wires=[2])
-                    @ qml.PauliX(wires=[3]),
-                    qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[2]),
-                    qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[3]),
-                    qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[2]),
-                    qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[3]),
-                    qml.PauliZ(wires=[2]) @ qml.PauliZ(wires=[3]),
-                ],
-            ),
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.40104295]], requires_grad=False),
             4,
             [
                 qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
@@ -349,10 +303,13 @@ def test_generate_paulis(generators, num_qubits, result):
         ),
     ],
 )
-def test_generate_symmetries(hamiltonian, num_qubits, res_generators, res_pauli_ops):
+def test_generate_symmetries(symbols, geometry, num_qubits, res_generators, res_pauli_ops):
     r"""Test that generate_symmetries returns the correct result."""
 
+    mol = qml.hf.Molecule(symbols, geometry)
+    hamiltonian = qml.hf.generate_hamiltonian(mol)()
     generators, pauli_ops = generate_symmetries(hamiltonian, num_qubits)
+
     for g1, g2 in zip(generators, res_generators):
         assert g1.compare(g2)
 
@@ -499,3 +456,89 @@ def test_transform_hamiltonian(symbols, geometry, generator, paulix_ops, paulix_
     for i, term in enumerate(sorted_terms):
         assert np.allclose(term[0], ham_ref.terms[0][i])
         assert term[1].compare(ham_ref.terms[1][i])
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "result"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.40104295]], requires_grad=False),
+            -1.137270174,
+        ),
+    ],
+)
+def test_sector_energy(symbols, geometry, result):
+    r"""Test that _sector_energy returns the correct result."""
+    mol = qml.hf.Molecule(symbols, geometry)
+    hamiltonian = qml.hf.generate_hamiltonian(mol)()
+    energy = _sector_energy(hamiltonian)
+    assert np.isclose(energy, result)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "generators", "paulix_ops", "result"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.40104295]], requires_grad=False),
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(3)]),
+            ],
+            [qml.PauliX(wires=[1]), qml.PauliX(wires=[2]), qml.PauliX(wires=[3])],
+            ([1, -1, -1], -1.137270174),
+        ),
+    ],
+)
+@pytest.mark.parametrize(("brute_force", "num_electrons"), [(True, 0), (False, 2)])
+def test_find_optimal_sector(
+    symbols, geometry, generators, paulix_ops, brute_force, num_electrons, result
+):
+    r"""Test that find_optimal_sector returns the correct result."""
+    mol = qml.hf.Molecule(symbols, geometry)
+    hamiltonian = qml.hf.generate_hamiltonian(mol)()
+
+    perm, energy = find_optimal_sector(
+        hamiltonian, generators, paulix_ops, brute_force, num_electrons
+    )
+
+    for sec in zip(perm, result[0]):
+        assert sec[0] == sec[1]
+    assert np.isclose(energy, result[1])
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "generators", "paulix_ops"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.40104295]], requires_grad=False),
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(3)]),
+            ],
+            [qml.PauliX(wires=[1]), qml.PauliX(wires=[2]), qml.PauliX(wires=[3])],
+        ),
+    ],
+)
+@pytest.mark.parametrize(
+    ("brute_force", "num_electrons", "msg_match"),
+    [
+        (False, 0, f"Brute force search is disabled"),
+        (False, 5, f"Number of active orbitals cannot be smaller than number of active electrons"),
+    ],
+)
+def test_exceptions_find_optimal_sector(
+    symbols, geometry, generators, paulix_ops, brute_force, num_electrons, msg_match
+):
+    r"""Test that find_optimal_sector returns the correct result."""
+    mol = qml.hf.Molecule(symbols, geometry)
+    hamiltonian = qml.hf.generate_hamiltonian(mol)()
+
+    with pytest.raises(ValueError, match=msg_match):
+        perm, energy = find_optimal_sector(
+            hamiltonian, generators, paulix_ops, brute_force, num_electrons
+        )