Implement LinearOperator protocol for SparsePauliOp

[kevinsung]

What should we add?

Doing linear algebra with a SparsePauliOp typically involves converting it to a sparse matrix, which takes a very long time. However, many linear algebra tasks don't require explicit storage of all of the matrix elements, but merely require knowledge of how the operator transforms a vector. SciPy has a protocol to capture this concept called LinearOperator, which is a level of abstraction above sparse matrix. We should implement this for SparsePauliOp. There are at least two possible ways to do this:

1. Have SparsePauliOp subclass LinearOperator, which would involve implementing at least the method _matvec that defines how the operator acts on a vector.
2. Make a separate function that converts a SparsePauliOp to a LinearOperator.

[jakelishman]

Copied from my Slack messages, for further context:

LinearOperator is applied to dense vectors (as Numpy arrays), so it’s going to have the same 2^n complexity implication (for storage of the statevector) as the CSR form of to_matrix. If that’s still useful, we can look into it.

I think we might not be able to have SparsePauliOp directly subclass LinearOperator because of incompatible subclass typing in the methods (though I could be wrong - I haven’t looked in detail, but SparsePauliOp definitely already provides some methods of those names, and the behaviours might not be compatible). If it’s not possible to do directly, then one thing we could do would be to add a to_linear_operator method to SparsePauliOp that returns an opaque LinearOperator that fulfils the interface.

That said, it’s still not going to let SparsePauliOp (or SparseObservable) scale to 30q for such linear algebra operations, because the statevector representation will still be bounded.

[T]he matrix/vector product will take a similar amount of time as the conversion to CSR, because it scales the same way.

Actually, tbf, there’s better ways of writing that algorithm so it doesn’t scale as poorly, especially if all the contained Paulis are low weight,
though that relies on the matvec method allowing you to mutate the vector in place, rather than returning a new one. If we have to return a new one, there’s no way to escape the 2^n scaling of the method.

[willkirby1]

FWIW I have an existing workaround for this (private repo) that does not use LinearOperators but just directly does matrix-free multiplication of SparsePauliOps on Scipy dok_arrays. I don't know if that is the optimal sparse vector format for performance, but what I have is scalable. In terms of practicality, it does allow multiplying e.g. a linear combination of 100000 computational basis vectors by a 40ish qubit SparsePauliOp. I haven't carefully checked how far it can be pushed. LMK if this is of interest, or if we are really wedded to using LinearOperators.

[jakelishman]

This is one more reason to use a to_linear_operator approach, so we don't cloud the main class with less efficient methods - dense arrays of course don't scale, even if the operator is matrix free.

Fwiw, dok_array will be totally fine if for computational basis vectors if you're manipulating them from Python and they stay very sparse. It's not so good for linear-algebra manipulation from Rust, though, because it stores all its data in Python formats, which makes it very slow for us to access. DOK is the cleanest for scalar indexing operations, since all of those are single hash lookups. For iterative update access (like matvec), the most efficient will be something that stores the indices in one array and the corresponding coefficients in another - a slight generalisation of the compressed-sparse formats, though iirc, scipy's implementation of those is limited to 2D. A csc format statevector with explicit shape (1, 2**n) is the closest to the ideal format, most likely (I didn't put a huge amount of thought into it), but indexing operations on that go as $\mathcal O\bigl(\lg(\text{nnz})\bigr)$ rather than constant-time (assuming the statevector is in canonical order - if it's unsorted, it's $\mathcal O\bigl(\text{nnz}\bigr)$, short of using Grover's lol).

[jakelishman]

When I say "very slow", I mean relative to what we should be able to achieve in Rust. Writing the same thing from Rust would still be approximately as fast as a pure Python-space implementation, maybe even a hair faster because of reduced loop overhead, but we'd still be bottlenecked on accesses to Python objects.