Repository: PennyLaneAI/pennylane
Pull request #2041 (2021-12-17 12:54:11)

Add functions to find correct sector for qubit tapering
**Context:**
This PR adds functions to find the correct sector for which the tapered Hamiltonian will correspond to the ground state of the system.

**Description of the Change:**
1. `_sector_energy` function returns the energy of the Hamiltonian in a given sector.
2. `find_optimal_sector` function returns the correct sector containing the ground state. It offers two strategies to do so - (i) by performing a brute force through all the possible permutations (exponential), and (ii) by building the sector using a relation between Pauli-Z operator and occupation number under JW transform (linear).

**Benefits:**
**Possible Drawbacks:**
**Related GitHub Issues:**

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#2041](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (158ae78) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/ebebb5e975374a0e98bf6b0d26fea5a8b1f7a29d?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (ebebb5e) will **decrease** coverage by `2.08%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #2041      +/-   ##
==========================================
- Coverage   99.17%   97.09%   -2.09%     
==========================================
  Files         225      225              
  Lines       17204    17240      +36     
==========================================
- Hits        17062    16739     -323     
- Misses        142      501     +359     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/hf/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2hmL19faW5pdF9fLnB5) | `100.00% <ø> (ø)` | |
| [pennylane/hf/tapering.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2hmL3RhcGVyaW5nLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/interfaces/batch/jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvYmF0Y2gvamF4LnB5) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/interfaces/jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvamF4LnB5) | `0.00% <0.00%> (-97.44%)` | :arrow_down: |
| [pennylane/devices/default\_qubit\_jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvZGVmYXVsdF9xdWJpdF9qYXgucHk=) | `0.00% <0.00%> (-96.21%)` | :arrow_down: |
| [pennylane/qnn/cost.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3Fubi9jb3N0LnB5) | `42.85% <0.00%> (-57.15%)` | :arrow_down: |
| [pennylane/transforms/classical\_jacobian.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvY2xhc3NpY2FsX2phY29iaWFuLnB5) | `80.00% <0.00%> (-20.00%)` | :arrow_down: |
| [pennylane/math/utils.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL21hdGgvdXRpbHMucHk=) | `81.81% <0.00%> (-18.19%)` | :arrow_down: |
| [pennylane/math/multi\_dispatch.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL21hdGgvbXVsdGlfZGlzcGF0Y2gucHk=) | `86.84% <0.00%> (-13.16%)` | :arrow_down: |
| [pennylane/math/is\_independent.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL21hdGgvaXNfaW5kZXBlbmRlbnQucHk=) | `88.50% <0.00%> (-11.50%)` | :arrow_down: |
| ... and [14 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [ebebb5e...158ae78](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2041?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

Jaybsoni: Looks good to me. Sorry for the delay! 

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index b01b28dbc..ab4d6fc8d 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -70,7 +70,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.
@@ -81,7 +82,7 @@
   mol = qml.hf.Molecule(symbols, geometry)
   H = qml.hf.generate_hamiltonian(mol)(geometry)
   generators, paulix_ops = qml.hf.generate_symmetries(H, len(H.wires))
-  paulix_sector = [1, -1, -1]
+  paulix_sector = qml.hf.optimal_sector(H, generators, mol.n_electrons)
   H_tapered = qml.hf.transform_hamiltonian(H, generators, paulix_ops, paulix_sector)
   ```
 
diff --git a/pennylane/hf/__init__.py b/pennylane/hf/__init__.py
index 0fdda9ec6..4d6410ce2 100644
--- a/pennylane/hf/__init__.py
+++ b/pennylane/hf/__init__.py
@@ -50,4 +50,5 @@ from .tapering import (
     generate_symmetries,
     get_generators,
     transform_hamiltonian,
+    optimal_sector,
 )
diff --git a/pennylane/hf/tapering.py b/pennylane/hf/tapering.py
index 11af063ff..12ab39da6 100644
--- a/pennylane/hf/tapering.py
+++ b/pennylane/hf/tapering.py
@@ -15,8 +15,9 @@
 This module contains the functions needed for tapering qubits using symmetries.
 """
 # pylint: disable=unnecessary-lambda
-import functools
 
+import functools
+import numpy
 import autograd.numpy as anp
 import pennylane as qml
 from pennylane import numpy as np
@@ -246,7 +247,7 @@ def generate_symmetries(qubit_op, num_qubits):
         num_qubits (int): number of wires required to define the Hamiltonian
 
     Returns:
-        tuple (list[Hamiltonian], list[Observable]):
+        tuple (list[Hamiltonian], list[Operation]):
 
             * list[Hamiltonian]: list of generators of symmetries, :math:`\mathbf{\tau}`,
               for the Hamiltonian.
@@ -464,3 +465,64 @@ def transform_hamiltonian(h, generators, paulix_ops, paulix_sector):
     c = qml.math.stack(c)
 
     return _simplify(qml.Hamiltonian(c, o))
+
+
+def optimal_sector(qubit_op, generators, active_electrons):
+    r"""Get the optimal sector which contains the ground state.
+
+    To obtain the optimal sector, we need to choose the right eigenvalues for the symmetry generators :math:`\bm{\tau}`.
+    We can do so by using the following relation between the Pauli-Z qubit operator and the occupation number under a
+    Jordan-Wigner transform.
+
+    .. math::
+
+        \sigma_{i}^{z} = I - 2a_{i}^{\dagger}a_{i}
+
+    According to this relation, the occupied and unoccupied fermionic modes correspond to the -1 and +1 eigenvalues of
+    the Pauli-Z operator, respectively. Since all of the generators :math:`\bm{\tau}` consist only of :math:`I` and
+    Pauli-Z operators, the correct eigenvalue for each :math:`\tau` operator can be simply obtained by applying it on
+    the reference Hartree-Fock (HF) state, and looking at the overlap between the wires on which the Pauli-Z operators
+    act and the wires that correspond to occupied orbitals in the HF state.
+
+    Args:
+        qubit_op (Hamiltonian): Hamiltonian for which symmetries are being generated to perform tapering
+        generators (list[Hamiltonian]): list of symmetry generators for the Hamiltonian
+        active_electrons (int): The number of active electrons in the system for generating the Hartree-Fock bitstring
+
+    Returns:
+        list[int]: eigenvalues corresponding to the optimal sector which contains the ground state
+
+    **Example**
+
+    >>> symbols = ['H', 'H']
+    >>> geometry = np.array([[0., 0., -0.66140414], [0., 0., 0.66140414]])
+    >>> mol = qml.hf.Molecule(symbols, geometry)
+    >>> H = qml.hf.generate_hamiltonian(mol)(geometry)
+    >>> generators, paulix_ops = qml.hf.generate_symmetries(H, len(H.wires))
+    >>> qml.hf.optimal_sector(H, generators, 2)
+        [1, -1, -1]
+    """
+
+    if active_electrons < 1:
+        raise ValueError(
+            f"The number of active electrons must be greater than zero;"
+            f"got 'electrons'={active_electrons}"
+        )
+
+    num_orbitals = len(qubit_op.wires)
+
+    if active_electrons > num_orbitals:
+        raise ValueError(
+            f"Number of active orbitals cannot be smaller than number of active electrons;"
+            f" got 'orbitals'={num_orbitals} < 'electrons'={active_electrons}."
+        )
+
+    hf_str = np.where(np.arange(num_orbitals) < active_electrons, 1, 0)
+
+    perm = []
+    for tau in generators:
+        symmstr = np.array([1 if wire in tau.ops[0].wires else 0 for wire in qubit_op.wires])
+        coeff = -1 if numpy.logical_xor.reduce(numpy.logical_and(symmstr, hf_str)) else 1
+        perm.append(coeff)
+
+    return perm


Test patch:
diff --git a/tests/hf/test_tapering.py b/tests/hf/test_tapering.py
index 6511e91a3..aa88fbb87 100644
--- a/tests/hf/test_tapering.py
+++ b/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,7 @@ from pennylane.hf.tapering import (
     generate_symmetries,
     get_generators,
     transform_hamiltonian,
+    optimal_sector,
 )
 
 
@@ -286,59 +287,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 +302,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 +455,85 @@ 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", "charge", "generators", "num_electrons", "result"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.40104295]], requires_grad=False),
+            0,
+            [
+                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)]),
+            ],
+            2,
+            [1, -1, -1],
+        ),
+        (
+            ["He", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4588684632]], requires_grad=False),
+            1,
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(1) @ qml.PauliZ(3)]),
+            ],
+            1,
+            [-1, 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=False,
+            ),
+            1,
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2) @ qml.PauliZ(4)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(1) @ qml.PauliZ(3) @ qml.PauliZ(5)]),
+            ],
+            2,
+            [-1, -1],
+        ),
+    ],
+)
+def test_optimal_sector(symbols, geometry, charge, generators, num_electrons, result):
+    r"""Test that find_optimal_sector returns the correct result."""
+    mol = qml.hf.Molecule(symbols, geometry, charge)
+    hamiltonian = qml.hf.generate_hamiltonian(mol)()
+
+    perm = optimal_sector(hamiltonian, generators, num_electrons)
+
+    assert perm == result
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "generators"),
+    [
+        (
+            ["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)]),
+            ],
+        ),
+    ],
+)
+@pytest.mark.parametrize(
+    ("num_electrons", "msg_match"),
+    [
+        (0, f"The number of active electrons must be greater than zero"),
+        (5, f"Number of active orbitals cannot be smaller than number of active electrons"),
+    ],
+)
+def test_exceptions_optimal_sector(symbols, geometry, generators, 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):
+        optimal_sector(hamiltonian, generators, num_electrons)