Add "layout" to the SparsePauliOp (and related operators)

[mrossinek]

What should we add?

The Pitch (or TL;DR)

For certain applications and operations I keep myself wishing for a builtin "layout" inside the SparsePauliOp. Note, that I am using "layout" due to the lack for a better word (since qargs is already taken).

Example

I think my particular application is best described with an example:

Consider I have a 127-qubit QuantumCircuit and would like to define a single-qubit Z observable on qubit index 13. While there are many means to construct a SparsePauliOp to represent this, in the end, the resulting object will always contain 126 identities.

>>> SparsePauliOp.from_sparse_list([("Z", [13], 1.0)], 127)
SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIII'],
              coeffs=[1.+0.j])

Not only is it pretty hard to tell on which index the Z lies, but any subsequent mathematical operations with this object will also be "expensive".
For example. consider I want to compose this operator with an XX on qubits 12 and 13. This would "perform unneeded math" on all the identities surrounding the very local observables.

Granted, for many arithmetic operations, one can set the qargs to perform local arithmetics. However, keeping track of the qargs alongside the operators can become cumbersome and cluttering.

---

Simple proposal

Thus, I propose some sort of "layout" or "qargs" (but this is already used) attribute that is internal to the SparsePauliOp (and likely PauliList and maybe BasePauli by extension).

Again, I am using "layout" due to a lack for a better word. I am open to suggestions!

Basically, here is a proposed idea:

a = SparsePauliOp("Z")
a.layout = [13]
a.global_num_qubits = 127  # again, open for naming suggestions

This defines a SparsePauliOp with a builtin layout. This object has an internal PauliList of length 1 but understands that it acts on a qubit index 13 out of 127 (hence the "layout" name suggestion).

Later, one can also define the XX similarly and do math as follows:

b = SparsePauliOp("XX")
b.layout = [12, 13]
b.global_num_qubits = 127

ab = a @ b
print(ab)
print(ab.layout, ab.global_num_qubits)

which will result in:

SparsePauliOp(['YX'],
              coeffs=[-0.+1.j])
[12, 13] 127

The key takeaways are, that even though we performed arithmetic operations, the operator remained small and, thus, avoids unnecessary computations.
I would also expect this to work on an XX with a layout [7, 13] which should still output a 2-qubit overall sized object:

SparsePauliOp(['XY'],
              coeffs=[-0.+1.j])
[13, 7] 127

---

Naive PoC

SparsePauliOp.apply_layout already makes an implementation of this fairly straight forward. Here is a naive insertion that I added to SparsePauliOp.compose to make the example above actually functional:

        composed_layout = None
        global_num_qubits = None
        if qargs is None:
            do_layouting = self.layout is not None or other.layout is not None
            if do_layouting:
                composed_layout: set[int] | list[int] = set()
                if self.layout is not None:
                    composed_layout = set(self.layout)
                else:
                    composed_layout = set(range(self.num_qubits))
                if other.layout is not None:
                    composed_layout.update(set(other.layout))
                else:
                    composed_layout = set(range(other.num_qubits))
                composed_layout = list(composed_layout)

                global_num_qubits = max(
                    self.global_num_qubits or self.num_qubits,
                    other.global_num_qubits or other.num_qubits,
                )

                if self.layout is not None:
                    self = self.apply_layout(
                        [composed_layout.index(q) for q in self.layout],
                        len(composed_layout),
                    )
                if other.layout is not None:
                    other = other.apply_layout(
                        [composed_layout.index(q) for q in other.layout],
                        len(composed_layout),
                    )

This is, of course, by no means a production quality solution, but it shows the basic idea and it shows, that with the addition of apply_layout we have already taken a step into this direction.

---

Open Question

• how would this need to interact with OpShape?
• what is the "best" level to implement this on? Should the layout and global_num_qubits better be placed on PauliList and leave SparsePauliOp mostly untouched? I think this should be possible to a large degree. The same could probably be said for the apply_layout method (see also Add apply_layout to other operators #11824)
• should we wrap layout and global_num_qubits into a single object of some sort?

---

Use cases

• easier "local" arithmetics
• step towards "sparse SparsePauliOp" (caveat being that all terms inside a SparsePauliOp are constrained to the same sparse qubit indices; so this is not quite fully sparse yet)
• routing of observables (possibly of interest to the runtime team)

---

I am happy to work on a PR to implement this. @caleb-johnson had also expressed interest in helping out with this.

[garrison]

a = SparsePauliOp("Z")
a.layout = [13]
a.global_num_qubits = 127  # again, open for naming suggestions

I think it is important for any solution to this problem to be elegantly composable. One should be able to create an operator in some succinct way with a Pauli Z on the 13th qubit (as you have done), stored sparsely. But one should also be able to compose that with a layout (e.g. a TranspiledLayout) to get an even bigger operator in the same form. For instance, I should be able to apply a layout on a 133 qubit device to the operator a above and work with it as a new, well defined SparsePauliOp. In this case, should layout be simply a list, or some object that represents this history? I suppose a list-like object is probably best because it would make working with it the most efficient.

[jakelishman]

Sorry this didn't have much of a response. To link to other threads of work that we have going on, we're going to be implementing #12282 as part of our six-month performance push, which heavily touches on the same kinds of concerns as this, as I understand it.

Let me preface this with: having sparse-per-qubit operators, as mooted by this issue, is absolutely something I agree is a good idea and super useful. Beyond that, though, I'm not certain that putting it onto SparsePauliOp would be super successful/useful:

1. just from an engineering perspective, there's a lot of code out there that uses SparsePauliOp and the symplectic structure directly, and adding an extra attribute onto it that changes the meaning of the operator is going to be hard to do in-place without causing problems for existing code.
2. the tendency to do arithmetic on these operators is likely to stop them from being qubit-sparse very quickly, which means we then end up with an operator that's not optimised to store the qubit-sparse components, and because of the partial support, is no longer well-optimised to operate in the general case fast either.

Point one I think is relatively self-explanatory - there's ways to mitigate it, but it'll always be tricky given then extant body of code that already relies on SparsePauliOp, including as an input format (primitives).

To explain my worries around point two more: this makes it such that now most operations on SparsePauliOp can no longer be direct bitwise operations on the two sets of symplectic matrices, but first need some sort of temporary expansion (say if there's a layout [0, 1] interacting with [0, 2]) or qubit re-ordering ([0, 1] interacting with [1, 0]), and even operations that don't need to do that would still have to run through the layouts to compare them before doing anything. There'd also be the risk that things like parity operators like $Z_0 + Z_1 + Z_2 + \dotsb$ would erase the sparsity structure entirely (since all summands need to be the same width in the SparsePauliOp symplectic format).

---

From this, then, the way I'm looking at the problem of "we need an operator that's efficient to manipulate where the terms are sparse over qubits" is more in favour of a new operator that builds in the sparsity structure to individual summands, and that that's something that would be cleaner and more performant to build into a solution to #12282 than to put on the top of SparsePauliOp and risk performance and affecting existing interfaces.

If we wanted to move ahead with the general principles laid out in this issue, I think something to consider is whether it would be better to make a wrapper around SparsePauliOp that adds the concept of layout, similar to how a Qiskit Gate is only defined as the natural width of the gate, but QuantumCircuit stores a (gate, operands) context pair, and interacts gates that way. A similar thing could be done with a wrapper object around SparsePauliOp (though more likely, (Pauli, coeff) pairs) . That said, that's kind of just coming down to how I'm thinking about #12282 as well.

[jakelishman]

Fwiw, at the ~100-qubit scale, I think the memory savings of a fully qubit-sparse operator vs a bitpacked full symplectic representation of the Pauli terms (i.e. just SparsePauliOp, but bitpacked in Rust space rather than byte-aligned Booleans like here) are rather less than one might think once things like the extra sparse overhead and internal data alignment have been fully considered. But at the several-thousand qubit scale it's obviously much more of an improvement, and the real kicker with #12282 is the extended alphabet to avoid exponential explosion of the terms for things like magnetisation operators.