Apply a PauliMeasurement to a PauliFrame

[J-C-Q]

Describe the bug 🐞

It is not possible to apply a PauliMeasurement to a PauliFrame.

Expected behavior

Using the pftrajectories() function, I thought I would be able to simulate a circuit that includes PauliMeasurements, like PauliMeasurement(P"ZZ", 1)

Minimal Reproducible Example 👇

circuit = [sHadamard(1),sHadamard(2),PauliMeasurement(P"ZZ", 1)]
pftrajectories(circuit; trajectories=1000)

Error & Stacktrace or other complete output produced by the MRE ⚠️

┌ Warning: Could not compactify the circuit, falling back to a slower version of the simulation. Consider reporting this issue to the package maintainers to improve performance. The offending gate was `PauliMeasurement{Array{UInt8, 0}, Vector{UInt64}}(+ ZZ, 1)`.
└ @ QuantumClifford ~/.julia/packages/QuantumClifford/haClW/src/pauli_frames.jl:203
ERROR: MethodError: no method matching apply!(::PauliFrame{Stabilizer{…}, Matrix{…}}, ::PauliMeasurement{Array{…}, Vector{…}})
The function `apply!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  apply!(::Register, ::PauliMeasurement{A, B}) where {A, B}
   @ QuantumClifford ~/.julia/packages/QuantumClifford/haClW/src/classical_register.jl:86
  apply!(::Register, ::Any, Any...; kwargs...)
   @ QuantumClifford ~/.julia/packages/QuantumClifford/haClW/src/classical_register.jl:36
  apply!(::QuantumClifford.AbstractQCState, ::NoisyGate)
   @ QuantumClifford ~/.julia/packages/QuantumClifford/haClW/src/noise.jl:125
  ...

Stacktrace:
 [1] pftrajectories(state::PauliFrame{Stabilizer{QuantumClifford.Tableau{Vector{UInt8}, LinearAlgebra.Adjoint{UInt64, Matrix{…}}}}, Matrix{Bool}}, circuit::Vector{QuantumClifford.AbstractOperation})
   @ QuantumClifford ~/.julia/packages/QuantumClifford/haClW/src/pauli_frames.jl:229
 [2] _pftrajectories(circuit::Vector{QuantumClifford.AbstractOperation}; trajectories::Int64, threads::Bool)
   @ QuantumClifford ~/.julia/packages/QuantumClifford/haClW/src/pauli_frames.jl:217
 [3] _pftrajectories
   @ ~/.julia/packages/QuantumClifford/haClW/src/pauli_frames.jl:196 [inlined]
 [4] #pftrajectories#198
   @ ~/.julia/packages/QuantumClifford/haClW/src/pauli_frames.jl:187 [inlined]
 [5] top-level scope
   @ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

[0525e862] QuantumClifford v0.9.12

• Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)

[c3fe647b] AbstractAlgebra v0.43.9
  [ec485272] ArnoldiMethod v0.4.0
  [62783981] BitTwiddlingConvenienceFunctions v0.1.6
  [861a8166] Combinatorics v1.0.2
  [f70d9fcc] CommonWorldInvalidations v1.0.0
  [34da2185] Compat v4.16.0
  [864edb3b] DataStructures v0.18.20
  [ffbed154] DocStringExtensions v0.9.3
  [86223c79] Graphs v1.12.0
  [3e5b6fbb] HostCPUFeatures v0.1.17
⌅ [2cd5bd5f] ILog2 v0.2.4
  [615f187c] IfElse v0.1.1
  [d25df0c9] Inflate v0.1.5
  [692b3bcd] JLLWrappers v1.6.1
  [1914dd2f] MacroTools v0.5.13
  [2edaba10] Nemo v0.47.3
  [bac558e1] OrderedCollections v1.6.3
  [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [0525e862] QuantumClifford v0.9.13
  [5717a53b] QuantumInterface v0.3.6
  [fb686558] RandomExtensions v0.4.4
  [fdea26ae] SIMD v3.6.0
  [699a6c99] SimpleTraits v0.9.4
  [aedffcd0] Static v1.1.1
  [90137ffa] StaticArrays v1.9.8
  [1e83bf80] StaticArraysCore v1.4.3
  [10745b16] Statistics v1.11.1
  [8e1ec7a9] SumTypes v0.5.8
  [e134572f] FLINT_jll v300.100.300+0
  [656ef2d0] OpenBLAS32_jll v0.3.28+3
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.11.0
  [de0858da] Printf v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [1a1011a3] SharedArrays v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test v1.11.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [781609d7] GMP_jll v6.3.0+0
  [deac9b47] LibCURL_jll v8.6.0+0
  [e37daf67] LibGit2_jll v1.7.2+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [3a97d323] MPFR_jll v4.2.1+0
  [c8ffd9c3] MbedTLS_jll v2.28.6+0
  [14a3606d] MozillaCACerts_jll v2023.12.12
  [4536629a] OpenBLAS_jll v0.3.27+1
  [bea87d4a] SuiteSparse_jll v7.7.0+0
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.59.0+0
  [3f19e933] p7zip_jll v17.4.0+2

• Output of versioninfo()

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 24 × 13th Gen Intel(R) Core(TM) i7-13700
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 1 default, 0 interactive, 1 GC (on 24 virtual cores)

Additional context

I'm guessing there just isn't a function for this type yet? How difficult do you think it would be to implement this? Or am I not getting something here why this is not possible?

Thanks!

[Fe-r-oz]

Thank you for highlighting this issue.

Based on the stack trace, the closest match is apply!(::Register, ::PauliMeasurement{A, B}) where {A, B}. You might try the following approach in the meantime:

julia> bell_state = ghz(2);

julia> m = PauliMeasurement(P"ZZ", 1);

julia> g1, g2 = sHadamard(1),sHadamard(2);

julia> v = VerifyOp(bell_state,[1,2]);

julia> reg = Register(bell_state, zeros(Bool, 1));

julia> circuit =  [g1,g2,m,v];

julia> petrajectories(reg, circuit)
Dict{CircuitStatus, Float64} with 3 entries:
  false_success:CircuitStatus(2) => 0.0
  true_success:CircuitStatus(1)  => 1.0
  failure:CircuitStatus(3)       => 0.0

julia> mctrajectories(reg, circuit, trajectories=1000)
Dict{CircuitStatus, Float64} with 4 entries:
  false_success:CircuitStatus(2) => 0.0
  true_success:CircuitStatus(1)  => 1000.0
  failure:CircuitStatus(3)       => 0.0
  continue:CircuitStatus(0)      => 0.0

[Krastanov]

Hi, @J-C-Q, Quinten. You are right, this is currently not supported for a few non-technical reasons. Below I describe our reasoning and the way we would currently perform such measurements in simulation. Please let me know if they address your needs -- we can also discus potential additions to the library that more closely match what you want to set up.

Option 1: Alternative simulators (slower but providing additional functionality)

As @Fe-r-oz mentioned, we have a Monte Carlo simulator that supports this feature (but it is naturally slower as it does not use Pauli frames) and an interesting "Symbolic Perturbative Expansion" simulator that can give you symbolic results instead of just numerical estimates (also much slower though).

Option 2: Using an ancilla qubit

The Pauli frame simulator supports single qubit measurements with the sMRX/sMRY/sMRZ operations. Usually any multi-qubit Pauli measurement would be decomposed into a single-qubit measurement of an ancillary qubit which was entangled with the qubits being measured. Usually we would want to be able to simulate the noise in the circuit implementing this multi-qubit measurement. Because of these requirements we have a tool to automatically convert multi-qubit measurements to the necessary single-qubit-ancilla measurement. E.g. when we run an error-correcting code, we do not want to assume that all multi-qubit measurements are perfect or even equally-good. Below is an example how one would do that:

julia> using QuantumClifford, QuantumClifford.ECC, Quantikz

julia> parity_checks = S"ZZ" # add more rows if you want more checks
+ ZZ

julia> circuit, nbanc, resultbits = naive_syndrome_circuit(parity_checks) # you can change the indices of the aux qubits and bit-storage locations
(QuantumClifford.AbstractOperation[sCNOT(1,3), sCNOT(2,3), sMRZ(3, 1)], 1, 1:1)

julia> full_circuit = [sHadamard(1), sHadamard(2), circuit...]

julia> traj_result = pftrajectories(full_circuit; trajectories=1000);

julia> pfmeasurements(traj_result) |> mean # from Statistics
0.535

Option 3: Implementing a Pauli frame runner for PauliMeasurement

We can certainly do that, it just has not been a priority of ours because of other low-hanging fruit. If you are interested in contributing this, we would be happy to help you navigate the pull request process.

Option 4: Just a commutativity check

If you do not care about simulating the circuit-level noise (which requires the decomposition shown in option 2), then maybe just checking the commutativity between the current Pauli frame and the Pauli operator you want to measure is all you need. That is also much faster:

julia> pre_measurement_circuit = [sHadamard(1), sHadamard(2)]

julia> traj_result = pftrajectories(pre_measurement_circuit; trajectories=1000);

julia> commutativities = comm(traj_result.frame, P"ZZ"); # TODO we should have a getter for traj_result.frame

julia> mean(commutativities)
0.491

---

@Fe-r-oz , Quinten is asking about Pauli frame simulations, which are much faster than other methods of simulation. Your approach (using the other two types of simulators we have, the "Monte Carlo over Stabilizer states" and the "Symbolic Perturbative Expansion" simulators) happen to support multi-qubit Pauli measurements, but they would be many orders of magnitude slower and does not necessarily address Quinten's request.

[J-C-Q]

Amazing!
Thank you very much for this super detailed answer.

Option 2 is what I thought I could avoid, but is indeed what I want. Since I want mid-circuit measurements, option 4 sadly doesn't work for my use case.

I would be thrilled to give option 3 a try. However, a low-level implementation seems a bit over my head right now. What do you think about a layer of abstraction to make it work, where we basically insert option 2 here under the hood? Similar to how apply!(frame::PauliFrame, op::sMX) is implemented at the moment. Although # TODO implement a faster direct version makes me think this is not desired. But maybe I could start with that so that the interface works out for now.

[Krastanov]

I like your suggestion. It is indeed possible to do "option 2 under the hood", but the problem would be that ancillary qubits will need to be made available in the Pauli frame object. That is doable, but would be annoying to track.

I would suggest the following:

1. In the PauliFrame constructor, always make sure that you have one extra qubit in the internal tableau
2. Implement your "option 2 under-the-hood" using that extra qubit
3. Submit a PR without worrying about performance or tests to discuss
4. Potentially we might want to discuss some changes to how the extra qubit is counted, to make sure the rest of the library does not crash and burn
5. Potentially we might discuss performance improvements and additional tests
6. Merging :)