Add support for sparse measurements

[J-C-Q]

Motivation

Measurements of PauliOperators on an n qubit state are currently either performed by checking the commutation relation between each row and the whole n qubit PauliOperator, or for single qubit measurements using the specialized functions projectX!, ... This specialized algorithm is great. However, since stabilizer measurements usually contain more than one qubit, a version optimized for more generalized sparse PauliOperators could be useful.

While the code for single qubit measurement can (more or less) easily be generalized to two qubits, the memory layout of the tableau makes this sparse checking less efficient.

Implementation

This PR tries to find a balance between efficiently using the memory layout of the tableau and sparsely computing commutation.
This is done by only computing the xor (in the comm function) for UInt64s in the tableau row that contain the sparsely selected qubits. This results in sparse scenarios where the qubits are close together (meaning packed into the same UInt64) to be very fast. For more uniformly distributed qubit selections, the speed-up is not as significant, or even non-existent. Nevertheless, for very large system sizes, measurements with a few qubits (e.g. 3 or 4) qubits are still faster this way. Especially, two-qubit parity checks.

However, right now, this requires computing the indices of the UInt64s in which the qubits are packed, which does allocate. Maybe there is a better way to do this.

TODO

There is still some cleanup to be done. As well as documentation and testing. However, I thought I'd share this idea already to get feedback if you are even interested in integrating this. I'm very open to any criticism and ideas to improve my code, as well as happy to work on this more.

I could also try to implement a two qubit parity measurement specialized function (analogous to the singe qubit functions), that would probably beat this one, if that is something you are interested in.

Benchmark

I did some rudimentary benchmarking. Attached you can find a plot of the results for a system of 2048 qubits. The time is in seconds.

The base case represents the n qubit PauliOperator, measured using the existing project! function.

n = 2048
k = 2
m = MixedDestabilizer(one(Stabilizer, n))
p = embed(n,rand(1:n,k),random_pauli(k))
@belapsed  QuantumClifford.project!($m, $p; keep_result=false,phases=false)

The best case represents the results for densely packed spares measurements.

n = 2048
k = 2
m = MixedDestabilizer(one(Stabilizer, n))
qubits = collect(1:k)
local_pauli = random_pauli(k)
@belapsed  QuantumClifford.project_sparse!($m, $qubits, $local_pauli;keep_result=false,phases=false)

The worst case for uniformly (maybe not the best way of doing this) distributed measurements.

n = 2048
k = 2
m = MixedDestabilizer(one(Stabilizer, n))
qubits = ceil.(Int,collect(range(1,n,k)))
local_pauli = random_pauli(k)
@belapsed  QuantumClifford.project_sparse!($m, $qubits, $local_pauli;keep_result=false,phases=false)

I'm not sure what causes the jumps in the worst case.

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 81 lines in your changes missing coverage. Please review.

Project coverage is 72.39%. Comparing base (2d199a8) to head (6c2525b).

Files with missing lines	Patch %	Lines
src/project_trace_reset.jl	0.00%	81 Missing ⚠️

❗ There is a different number of reports uploaded between BASE (2d199a8) and HEAD (6c2525b). Click for more details.

HEAD has 3 uploads less than BASE

Flag	BASE (2d199a8)	HEAD (6c2525b)
	11	8

Additional details and impacted files

@@             Coverage Diff             @@
##           master     #544       +/-   ##
===========================================
- Coverage   83.24%   72.39%   -10.86%     
===========================================
  Files          79       79               
  Lines        5235     5311       +76     
===========================================
- Hits         4358     3845      -513     
- Misses        877     1466      +589     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Krastanov]

This is interesting. Do you have a use case you are developing this for, or just an interesting optimization?

I would suggest introducing a new comm method that takes an extra argument qubits, and then the indexing calculation happens on the fly in the loop.

Overall though, I think any SIMD speedup gets removed here, so in many situations I would expect the naive dense approach to be faster.

[J-C-Q]

Well, my use case are measurement only circuits with lots of parity measurements of a few (2-6) qubits. But I'm guessing something like this could also be useful for QEC stabilizer measurements?!

Indeed, I think this is only an optimization for very large systems.