Repository: entropicalabs/openqaoa (language: python)

Issue:
Implement `as_matrix` method for the `Hamiltonian` class
### Prerequisites
Before raising this feature request, I made sure:

- [x] It does not lead to duplication.
- [x] It has not already been requested.

### Is your feature request related to a problem? Please describe.
The class `Hamiltonian` in OpenQAOA defined in `src/openqaoa-src/openqaoa/qaoa_components/ansatz_constructor/operators.py` is responsible for constructing the `Cost` and `Mixer` Hamiltonians in QAOA. It fundamentally describes the Hamiltonian object via list of `PauliOp`s defined on a set of qubits and a corresponding list of their coefficients. Implicitly, the definition of the Hamiltonian is then assumed as follows: $$H = \sum_i \alpha_{i} P_i + \sum_{ij} J_{ij} P_{ij},$$ where $\ P_i, P_{ij}$ are one-qubit and two-qubit `PauliOp`s respectively. 

Implement a method in this class that returns the matrix representation of the `Hamiltonian` object. This will be a $2^n x 2^n$ dimensional sum over all Pauli entries as shown above. 

### Describe the solution you'd like
- Implement the method
- Use numpy sparse matrix where relevant
- Implement checks preventing users from creating matrices beyond 20 qubits to prevent memory overflow
- Bonus task: If the Hamiltonian can be expressed as a sparse matrix, shift the qubit restriction to a higher number depending on the sparsity of the matrix

```{python}
class Hamiltonian:
    def as_matrix(self):
        "Implements a numpy array (sparse where relevant) representation
         of the Hamiltonian object."
        #TODO
        return mat
```


Discussion / hints:
Hello @vishal-ph! I am new to OpenQAOA and woud like to implement this functionality. Could you please assign it to me and if possible show any relavant literature to show matrix rep?

@manulpatel, thanks for taking the initiative! 
In OpenQAOA, the Hamiltonian object is a sum over the terms of the Pauli and Identity operators set tensored over n-qubits. The goal is to represent the sum of these terms in a matrix form. To get started, you may want to start looking at the matrix representation of the Pauli operators here https://en.wikipedia.org/wiki/Pauli_matrices and think about their tensor products.

Fix patch:
diff --git a/src/openqaoa-core/qaoa_components/ansatz_constructor/operators.py b/src/openqaoa-core/qaoa_components/ansatz_constructor/operators.py
index 113efb0a1..08871517b 100644
--- a/src/openqaoa-core/qaoa_components/ansatz_constructor/operators.py
+++ b/src/openqaoa-core/qaoa_components/ansatz_constructor/operators.py
@@ -5,6 +5,9 @@
 from typing import List, Union, Tuple
 from sympy import Symbol
 import numpy as np
+from scipy.sparse import kron as sparse_kron
+from functools import reduce
+from scipy.sparse import csr_matrix
 
 Identity = np.array(([1, 0], [0, 1]), dtype=complex)
 PauliX = np.array(([0, 1], [1, 0]), dtype=complex)
@@ -686,3 +689,45 @@ def classical_hamiltonian(
                 raise ValueError("Hamiltonian only supports Linear and Quadratic terms")
 
         return cls(pauli_ops, pauli_coeffs, constant)
+
+    @property
+    def as_matrix(self):
+        """
+        Build sparse matrix of Hamiltonian. Note this includes self.constant value.
+
+        Returns
+        -------
+        H_mat:
+            sparse matrix of Hamiltonian.
+        """
+
+        if self.n_qubits>20:
+            raise ValueError('number of qubits exceeds maximum of 20')
+
+        p_mat_dict = {'X': csr_matrix([[0, 1],
+                                       [1, 0]]),
+                      'Y': csr_matrix([[0., -1j],
+                                       [1j, 0.]]),
+                      'Z': csr_matrix([[1, 0],
+                                       [0, -1]]),
+                      'I': csr_matrix([[1, 0], [0, 1]])}
+
+        converter = lambda s: p_mat_dict[s]
+
+        # vectorized function to map single qubit pauli strings to sparse matrix
+        vfunc = np.vectorize(converter)
+
+        base = np.array(['I' for _ in range(self.n_qubits)])
+
+        H_mat = reduce(sparse_kron, vfunc(base)) * self.constant
+        for P, coeff in zip(self.terms, self.coeffs):
+            temp = base.copy()
+
+            # replace single qubit pauli I in base term
+            # with single qubit pauli operators according to qubit_indices and pauli_strs in P
+            temp[np.array(P.qubit_indices)] = np.array(list(P.pauli_str))
+
+            # take kronecker product of sparse list of matrices
+            H_mat += reduce(sparse_kron, vfunc(temp)) * coeff
+
+        return H_mat
\ No newline at end of file


Test patch:
diff --git a/tests/test_operators.py b/tests/test_operators.py
index a3d9501a4..28f37cc86 100644
--- a/tests/test_operators.py
+++ b/tests/test_operators.py
@@ -1314,6 +1314,74 @@ def test_hamiltonian_dict(self):
             hamiltonian_dict_pauliop == correct_dict_pauliop
         ), "Hamiltonian dictionary did not yield correct dictionary"
 
+    def test_hamiltonian_as_matrix(self):
+        """Test if the method produces
+        the correct matrix corresponding to given Hamiltonian.
+
+        """
+        P_list = ['XIIZ', 'YIII', 'IZXI', 'IXIY', 'IIIZ', 'IIXI', 'ZZII']
+        coeff_list = list(range(1, 8))
+
+        Pauli_list = []
+        for P in P_list:
+            P = np.array(list(P))
+            non_I = np.where(P != 'I')[0]
+            Pauli_list.append(PauliOp(P[non_I], tuple(non_I)))
+
+        H = Hamiltonian(Pauli_list, coeff_list, 0)
+
+        solution = np.array([[12. + 0.j, 0. + 0.j, 9. + 0.j, 0. + 0.j, 0. + 0.j, 0. - 4.j,
+                              0. + 0.j, 0. + 0.j, 1. - 2.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, 2. + 0.j, 0. + 0.j, 9. + 0.j, 0. + 4.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, -1. - 2.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j],
+                             [9. + 0.j, 0. + 0.j, 12. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. - 4.j, 0. + 0.j, 0. + 0.j, 1. - 2.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j],
+                             [0. + 0.j, 9. + 0.j, 0. + 0.j, 2. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 4.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, -1. - 2.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j],
+                             [0. + 0.j, 0. - 4.j, 0. + 0.j, 0. + 0.j, -2. + 0.j, 0. + 0.j,
+                              3. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              1. - 2.j, 0. + 0.j, 0. + 0.j, 0. + 0.j],
+                             [0. + 4.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, -12. + 0.j,
+                              0. + 0.j, 3. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, -1. - 2.j, 0. + 0.j, 0. + 0.j],
+                             [0. + 0.j, 0. + 0.j, 0. + 0.j, 0. - 4.j, 3. + 0.j, 0. + 0.j,
+                              -2. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 1. - 2.j, 0. + 0.j],
+                             [0. + 0.j, 0. + 0.j, 0. + 4.j, 0. + 0.j, 0. + 0.j, 3. + 0.j,
+                              0. + 0.j, -12. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, -1. - 2.j],
+                             [1. + 2.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, -2. + 0.j, 0. + 0.j, 9. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. - 4.j, 0. + 0.j, 0. + 0.j],
+                             [0. + 0.j, -1. + 2.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, -12. + 0.j, 0. + 0.j, 9. + 0.j,
+                              0. + 4.j, 0. + 0.j, 0. + 0.j, 0. + 0.j],
+                             [0. + 0.j, 0. + 0.j, 1. + 2.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 9. + 0.j, 0. + 0.j, -2. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, 0. - 4.j],
+                             [0. + 0.j, 0. + 0.j, 0. + 0.j, -1. + 2.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, 9. + 0.j, 0. + 0.j, -12. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 4.j, 0. + 0.j],
+                             [0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 1. + 2.j, 0. + 0.j,
+                              0. + 0.j, 0. + 0.j, 0. + 0.j, 0. - 4.j, 0. + 0.j, 0. + 0.j,
+                              12. + 0.j, 0. + 0.j, 3. + 0.j, 0. + 0.j],
+                             [0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, -1. + 2.j,
+                              0. + 0.j, 0. + 0.j, 0. + 4.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              0. + 0.j, 2. + 0.j, 0. + 0.j, 3. + 0.j],
+                             [0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j,
+                              1. + 2.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j, 0. - 4.j,
+                              3. + 0.j, 0. + 0.j, 12. + 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. + 2.j, 0. + 0.j, 0. + 0.j, 0. + 4.j, 0. + 0.j,
+                              0. + 0.j, 3. + 0.j, 0. + 0.j, 2. + 0.j]])
+
+        assert np.allclose(solution,
+                    H.as_matrix.toarray()), 'as_matrix method not correct'
+
 
 if __name__ == "__main__":
     unittest.main()