Repository: PennyLaneAI/pennylane (python)

Issue:
Add matrix support for fermionic objects
### Feature details

⚠️ This issue is part of an internal assignment and not meant for external contributors.

The [fermi](https://docs.pennylane.ai/en/stable/code/qml_fermi.html) module has several functions and classes for creating and manipulating fermionic operators. It will be useful to have tools for obtaining the matrix representation of these operators. The [FermiWord](https://docs.pennylane.ai/en/stable/code/api/pennylane.fermi.FermiWord.html) and [FermiSentence](https://docs.pennylane.ai/en/stable/code/api/pennylane.fermi.FermiSentence.html) objects can be upgraded to have a `to_mat()` method that allows computing their matrix representation, similar to the `PauliWord` and `PauliSentence` classes. The users can then compute the matrix as

```
>>> import pennylane as qml
>>> op = qml.fermi.FermiWord({(0, 0): '+', (1, 1): '-'})
>>> op.to_mat()
array([[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]])
```

### Implementation
The `to_mat()` method should be added as an instance method to `FermiWord` and `FermiSentence`. The method should have an optional keyword argument, `n`,  for specifying the total number of orbitals which determines the matrix size. 

### Requirements
1. The computed matrix should properly describe the behaviour of fermionic operators.
2. The computed matrix should satisfy the commutation relationship between fermionic operators.

Pull request description:
Closes #5895 Add matrix support for fermionic objects
**Context:**
The [fermi](https://docs.pennylane.ai/en/stable/code/qml_fermi.html) module includes functions and classes for creating and manipulating fermionic operators: `FermiWord` and `FermiSentence` class. It would be beneficial to add support for obtaining the matrix representation of these operators, similar to the `to_mat` instance method in the `PauliWord` and `PauliSentence` classes.

**Description of the Change:**

- `to_mat` instance method is added to `FermiWord` and `FermiSentence` classes; the implementation for two methods is nearly identical -- using Jordan-Wigner transform to convert an instance of `FermiWord` or `FermiSentence` to a `PauliWord` or `PauliSentence` instance, then using the existing `to_mat` method in the two classes to compute the matrix representation. 
- `to_mat` takes `n_orbitals` as an optional input, which determines the matrix size $2^{n_\mathrm{orbitals}} \times 2^{n_\mathrm{orbitals}}$

**Benefits:**
The new method `to_mat` in `FermiWord` and `FermiSentence` class allows computing the matrix representation of a given Fermi operator.

**Possible Drawbacks:**
- Currently, `to_mat` in `FermiWord` and `FermiSentence` class only supports dense matrices. The exponential growth of matrix dimensions (a.k.a Fock space) with respect to `n_orbitals` (or the largest orbital index in the input fermionic operator) may limit its applicability.

**Related GitHub Issues:**
[#5895](https://github.com/PennyLaneAI/pennylane/issues/5895)

### Why Jordan-Wigner transform?
The choice of Jordan-Wigner mapping was primarily motivated by its direct and straightforward representation, where **each Slater determinant has a one-to-one mapping with the computational basis states of qubits**. This feature simplifies the interpretation of matrix elements significantly, as **each matrix element directly corresponds to the expectation values of the fermionic operator between these determinants**. Such a direct correspondence facilitates easier understanding and analysis of quantum states, making it intuitive and useful in the context of quantum chemical applications, as well as in educational and debugging scenarios.

In contrast, while Parity and Bravyi-Kitaev mappings offer computational efficiencies for larger systems and specific quantum hardware considerations, they encode fermionic information into qubit states in ways that are more complex and less intuitive. Parity mapping encodes the parity of the number of occupied orbitals up to each position, and Bravyi-Kitaev uses a hierarchical structure. 
- _Parity_: Matrix elements in the parity-mapped Hamiltonian correspond to interactions that conserve or flip these parity properties. Physically, each element can be thought of as involving transitions that maintain or alter the parity characteristics of the fermionic system, rather than direct particle occupations.
- _Bravyi-Kitaev_: Matrix elements in a Bravyi-Kitaev transformed Hamiltonian reflect operations on a hierarchical structure of fermionic occupations. Each qubit state relates to a specific combination of fermion occupation and non-occupation across different levels of this hierarchy, with matrix elements representing more complex interactions across these levels.

The physical interpretation of each matrix element in the two mappings requires a deeper understanding of the encoded structures, making them less transparent for straightforward interpretations (e.g. A transformation matrix may be necessary to establish the correspondence between qubit computational basis states and Slater determinants).

Including the option to select different mappings could potentially offer greater flexibility and optimization for users dealing with various types of quantum hardware or larger systems. However, it would also introduce additional complexity into our codebase and could complicate the user experience for those who are not as familiar with the subtleties of each mapping method. Therefore, I chose to focus exclusively on the Jordan-Wigner mapping for its clarity and ease of interpretation.


#### Notes on the design (OUTDATED, `_to_mat` function has been removed)
The primary motivation for centralizing the main logic in the `_to_mat` function, with both `FermiWord.to_mat()` and `FermiSentence.to_mat()` methods calling this function, is to reduce the code duplication and enhance maintainability. The implementation for these methods is nearly identical; hence, modularizing the code into a separate, reusable function that accepts either a `FermiWord` or `FermiSentence` instance seemed like the most efficient approach. Here are several advantages of the current implementation:

1. **Modularity and Reduced Code Duplication**: As mentioned, this approach helps in avoiding repetition of similar code across multiple classes, adhering to the DRY (Don't Repeat Yourself) principle.

2. **Flexibility in Usage**: By potentially making this function public (renaming it to `to_mat` and importing it in the `__init__.py`), users gain flexibility. They can either use the method associated with an instance:
   ```python
   word = FermiWord({(0, 0): "+", (1, 1): "-"})
   word.to_mat()
   ```
   Or directly call the function:
   ```python
   from pennylane.fermi.fermionic import to_mat
   to_mat(word)
   ```
Extending this idea, we could further develop it into a versatile `to_mat` function that accepts both Fermi and Pauli operators. This function would have the signature:
```python
def to_mat(operator: Union[FermiWord, FermiSentence, PauliWord, PauliSentence], **kwargs) -> ndarray
```
Such a design aligns with Python best practices of maintaining lightweight classes, as functions are generally easier to test than instance methods because they are less dependent on instance state. This approach not only simplifies testing but also enhances code usability and readability.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 6cd04adf4..8559de76b 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -21,6 +21,10 @@
 
 * `SProd.terms` now flattens out the terms if the base is a multi-term observable.
   [(#5885)](https://github.com/PennyLaneAI/pennylane/pull/5885)
+  
+* A new method `to_mat` has been added to the `FermiWord` and `FermiSentence` classes, which allows
+  computing the matrix representation of these Fermi operators.
+  [(#5920)](https://github.com/PennyLaneAI/pennylane/pull/5920)
 
 <h3>Improvements 🛠</h3>
 
@@ -146,6 +150,7 @@ Ahmed Darwish,
 Lillian M. A. Frederiksen,
 Pietropaolo Frisoni,
 Emiliano Godinez,
+Renke Huang,
 Christina Lee,
 Austin Huang,
 Christina Lee,
diff --git a/pennylane/fermi/fermionic.py b/pennylane/fermi/fermionic.py
index 01b5049e4..73d8935ce 100644
--- a/pennylane/fermi/fermionic.py
+++ b/pennylane/fermi/fermionic.py
@@ -273,6 +273,36 @@ class FermiWord(dict):
 
         return operator
 
+    def to_mat(self, n_orbitals=None):
+        r"""Return the matrix representation.
+
+        Args:
+            n_orbitals (int or None): Number of orbitals. If not provided, it will be inferred from
+                the largest orbital index in the Fermi operator.
+
+        Returns:
+            NumpyArray: Matrix representation of the :class:`~.FermiWord`.
+
+        **Example**
+
+        >>> w = FermiWord({(0, 0): '+', (1, 1): '-'})
+        >>> w.to_mat()
+        array([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
+              [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
+              [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
+              [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])
+        """
+        largest_orb_id = max(key[1] for key in self.keys()) + 1
+        if n_orbitals and n_orbitals < largest_orb_id:
+            raise ValueError(
+                f"n_orbitals cannot be smaller than {largest_orb_id}, got: {n_orbitals}."
+            )
+
+        largest_order = n_orbitals or largest_orb_id
+        mat = qml.jordan_wigner(self, ps=True).to_mat(wire_order=list(range(largest_order)))
+
+        return mat
+
 
 # pylint: disable=useless-super-delegation
 class FermiSentence(dict):
@@ -463,6 +493,36 @@ class FermiSentence(dict):
             if abs(coeff) <= tol:
                 del self[fw]
 
+    def to_mat(self, n_orbitals=None):
+        r"""Return the matrix representation.
+
+        Args:
+            n_orbitals (int or None): Number of orbitals. If not provided, it will be inferred from
+                the largest orbital index in the Fermi operator
+
+        Returns:
+            NumpyArray: Matrix representation of the :class:`~.FermiSentence`.
+
+        **Example**
+
+        >>> fs = FermiSentence({FermiWord({(0, 0): "+", (1, 1): "-"}): 1.2, FermiWord({(0, 0): "+", (1, 0): "-"}): 3.1})
+        >>> fs.to_mat()
+        array([0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+              [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+              [0.0 + 0.0j, 1.2 + 0.0j, 3.1 + 0.0j, 0.0 + 0.0j],
+              [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 3.1 + 0.0j])
+        """
+        largest_orb_id = max(key[1] for fermi_word in self.keys() for key in fermi_word.keys()) + 1
+        if n_orbitals and n_orbitals < largest_orb_id:
+            raise ValueError(
+                f"n_orbitals cannot be smaller than {largest_orb_id}, got: {n_orbitals}."
+            )
+
+        largest_order = n_orbitals or largest_orb_id
+        mat = qml.jordan_wigner(self, ps=True).to_mat(wire_order=list(range(largest_order)))
+
+        return mat
+
 
 def from_string(fermi_string):
     r"""Return a fermionic operator object from its string representation.


Test patch:
diff --git a/tests/fermi/test_fermionic.py b/tests/fermi/test_fermionic.py
index b5a9b237f..3e0c324e4 100644
--- a/tests/fermi/test_fermionic.py
+++ b/tests/fermi/test_fermionic.py
@@ -125,6 +125,22 @@ class TestFermiWord:
         with pytest.raises(ValueError, match="The operator indices must belong to the set"):
             FermiWord(operator)
 
+    def test_to_mat(self):
+        """Test that the matrix representation of FermiWord is correct."""
+
+        expected_mat = np.zeros((4, 4), dtype=complex)
+        expected_mat[2, 1] = 1.0
+
+        mat = fw1.to_mat()
+        assert np.allclose(mat, expected_mat)
+
+    def test_to_mat_error(self):
+        """Test that an error is raised if the requested matrix dimension is smaller than the
+        dimension inferred from the largest orbital index.
+        """
+        with pytest.raises(ValueError, match="n_orbitals cannot be smaller than 2"):
+            fw1.to_mat(n_orbitals=1)
+
 
 class TestFermiWordArithmetic:
     WORDS_MUL = (
@@ -442,6 +458,14 @@ fs2_hamiltonian = FermiSentence({fw1: -1.23, fw2: -4, fw3: 0.5})
 fs3 = FermiSentence({fw3: -0.5, fw4: 1})
 fs4 = FermiSentence({fw4: 1})
 fs5 = FermiSentence({})
+fs6 = FermiSentence({fw1: 1.2, fw2: 3.1})
+fs7 = FermiSentence(
+    {
+        FermiWord({(0, 0): "+", (1, 1): "-"}): 1.23,  # a+(0) a(1)
+        FermiWord({(0, 0): "+", (1, 0): "-"}): 4.0j,  # a+(0) a(0) = n(0) (number operator)
+        FermiWord({(0, 0): "+", (1, 2): "-", (2, 1): "+"}): -0.5,  # a+(0) a(2) a+(1)
+    }
+)
 
 fs1_x_fs2 = FermiSentence(  # fs1 * fs1, computed by hand
     {
@@ -600,6 +624,25 @@ class TestFermiSentence:
         new_fs = pickle.loads(serialization)
         assert fs == new_fs
 
+    def test_to_mat(self):
+        """Test that the matrix representation of FermiSentence is correct."""
+        expected_mat = np.zeros((8, 8), dtype=complex)
+        expected_mat[4, 2] = 1.23 + 0j
+        expected_mat[5, 3] = 1.23 + 0j
+        for i in [4, 5, 6, 7]:
+            expected_mat[i, i] = 4.0j
+        expected_mat[6, 1] = 0.5 + 0j
+
+        mat = fs7.to_mat()
+        assert np.allclose(mat, expected_mat)
+
+    def test_to_mat_error(self):
+        """Test that an error is raised if the requested matrix dimension is smaller than the
+        dimension inferred from the largest orbital index.
+        """
+        with pytest.raises(ValueError, match="n_orbitals cannot be smaller than 3"):
+            fs7.to_mat(n_orbitals=2)
+
 
 class TestFermiSentenceArithmetic:
     tup_fs_mult = (  # computed by hand