How to Pauli decompose a sparse matrix

PennyLane Help
8 
import pennylane as qml
import numpy as np
import scipy as sp
np.set_printoptions(linewidth=400)
arr = np.array([[0, 0, 0, 1, 0, 0, 0, 0],
       [0, 0, 1, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 0],
       [1, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 1, 0, 0],
       [0, 0, 0, 0, 1, 0, 0, 0]])
qml.pauli_decompose(arr)

row = np.array([0, 1, 2, 3, 4, 5, 6, 7])
col = np.array([3, 2, 1, 0, 7, 6, 5, 4])
data = np.array([1, 1, 1, 1, 1, 1, 1, 1])
sparse = sp.sparse.csr_array((data, (row, col)), shape=(8, 8))

qml.pauli_decompose(sparse)

I borrowed the example from @mchau, say I create some sparse (hermitian) matrix using scipy. I then want to turn this matrix into a sum of Pauli strings, as you can do with a numpy array. However, with the sparse matrix you cannot directly apply qml.pauli_decompose, it raises

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[6], line 20
     17 data = np.array([1, 1, 1, 1, 1, 1, 1, 1])
     18 sparse = sp.sparse.csr_array((data, (row, col)), shape=(8, 8))
---> 20 qml.pauli_decompose(sparse)

File /opt/conda/lib/python3.10/site-packages/pennylane/pauli/conversion.py:316, in pauli_decompose(H, hide_identity, wire_order, pauli, check_hermitian)
    313     if not is_abstract(H) and not qml.math.allclose(H, qml.math.conj(qml.math.transpose(H))):
    314         raise ValueError("The matrix is not Hermitian")
--> 316 coeffs, obs = _generalized_pauli_decompose(
    317     H, hide_identity=hide_identity, wire_order=wire_order, pauli=pauli, padding=True
    318 )
    320 if check_hermitian:
    321     coeffs = qml.math.real(coeffs)

File /opt/conda/lib/python3.10/site-packages/pennylane/pauli/conversion.py:129, in _generalized_pauli_decompose(matrix, hide_identity, wire_order, pauli, padding)
     34 r"""Decomposes any matrix into a linear combination of Pauli operators.
     35 
     36 This method converts any matrix to a weighted sum of Pauli words acting on :math:`n` qubits
   (...)
    126 
    127 """
    128 # Ensuring original matrix is not manipulated and we support builtin types.
--> 129 matrix = qml.math.convert_like(matrix, next(iter([*matrix[0]]), []))
    131 # Pad with zeros to make the matrix shape equal and a power of two.
    132 if padding:

File /opt/conda/lib/python3.10/site-packages/scipy/sparse/_base.py:262, in _spbase.__iter__(self)
    260 def __iter__(self):
    261     for r in range(self.shape[0]):
--> 262         yield self[r]

TypeError: 'coo_array' object is not subscriptable

So I was wondering if there is maybe some alternate way to make a linear combination of Pauli operators out of a sparse matrix.
Some extra things:

  • I am specifically interested the pauli sum itself, not some the expectation value, so from what I gather SparseHamiltonian does not work in this case.
  • I have considered turning the sparse matrix into a numpy array, but I want to avoid doing this, as the array becomes too large (I was looking at a 16 qubit system), so I would like to avoid this if at all possible.

Can you help me with this?