Merge pull request #243 from QuantumSavory/decodersspeedup

Decoder speedup improvements and breaking changes to how we set noise parameters to avoid confusion

CHANGELOG.md

@@ -2,10 +2,18 @@
 
 - `permute` will be a wrapper around to `QuantumInterface.permutesubsystems`. Documentation for `permute!` would be similarly updated
 - reworking the rest of `NoisyCircuits` and moving it out of `Experimental`
-- `random_pauli(...; realphase=true)` should be the default
 
 # News
 
+## v0.9.0 - 2024-03-19
+
+- **(breaking)** The defaults in `random_pauli` are now `realphase=true` and `nophase=true`.
+- **(breaking)** The convention for for flip probability in `random_pauli`.
+- **(breaking)** The convention for noise probability in `UnbiasedUncorrelatedNoise` changed. The input number is the total probability for an error to occur.
+- Implement an inplace `random_pauli!`, a non-allocating alternative to `random_pauli`.
+- Significant improvement in the performance of the ECC decoder pipeline (but many low-hanging fruits still remain).
+
+
 ## v0.8.21 - 2024-03-17
 
 - Implemented the Gottesman code family, also known as [[2^j, 2^j - j - 2, 3]] quantum Hamming codes.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <[email protected]>"]
-version = "0.8.21"
+version = "0.9.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

docs/src/commonstates.md

@@ -44,10 +44,10 @@ Random Pauli operators are implemented as well (with or without a random phase).
 julia> using StableRNGs; rng = StableRNG(42);
 
 julia> random_pauli(rng, 4)
-+i_ZZ_
++ ZYY_
 
-julia> random_pauli(rng, 4; nophase=true)
-+ ZXZY
+julia> random_pauli(rng, 4; nophase=false)
+- YZ_X
 ```
 
 ## Stabilizer States

ext/QuantumCliffordGPUExt/apply_noise.jl

@@ -3,7 +3,7 @@ using QuantumClifford: _div, _mod
 #according to https://github.com/JuliaGPU/CUDA.jl/blob/ac1bc29a118e7be56d9edb084a4dea4224c1d707/test/core/device/random.jl#L33
 #CUDA.jl supports calling rand() inside kernel
 function applynoise!(frame::PauliFrameGPU{T},noise::UnbiasedUncorrelatedNoise,i::Int) where {T <: Unsigned}
-    p = noise.errprobthird
+    p = noise.p
     lowbit = T(1)
     ibig = _div(T,i-1)+1
     ismall = _mod(T,i-1)
@@ -22,14 +22,15 @@ function applynoise_kernel(xzs::DeviceMatrix{Tme},
     p::Real,
     ibig::Int,
     ismallm::Tme,
-    rows::Int) where {Tme <: Unsigned} 
+    rows::Int) where {Tme <: Unsigned}
 
     f = (blockIdx().x - 1) * blockDim().x + threadIdx().x;
     if f > rows
         return nothing
     end
 
     r = rand()
+    p = p/3
     if  r < p # X error
         xzs[ibig,f] ⊻= ismallm
     elseif r < 2p # Z error

ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl

@@ -4,7 +4,7 @@ using PyQDecoders: np, sps, ldpc, pm, PythonCall
 using SparseArrays
 using QuantumClifford
 using QuantumClifford.ECC
-import QuantumClifford.ECC: AbstractSyndromeDecoder, decode, parity_checks
+import QuantumClifford.ECC: AbstractSyndromeDecoder, decode, batchdecode, parity_checks
 
 abstract type PyBP <: AbstractSyndromeDecoder end
 
@@ -102,4 +102,12 @@ function decode(d::PyMatchingDecoder, syndrome_sample)
     vcat(guess_x_errors, guess_z_errors)
 end
 
+function batchdecode(d::PyMatchingDecoder, syndrome_samples)
+    row_x = syndrome_samples[:,1:d.nx] # TODO This copy is costly!
+    row_z = syndrome_samples[:,d.nx+1:end]
+    guess_z_errors = PythonCall.PyArray(d.pyx.decode_batch(row_x))
+    guess_x_errors = PythonCall.PyArray(d.pyz.decode_batch(row_z))
+    hcat(guess_x_errors, guess_z_errors)
+end
+
 end

src/QuantumClifford.jl

@@ -22,7 +22,7 @@ export
     @T_str, xbit, zbit, xview, zview,
     @S_str, Stabilizer,
     Destabilizer, MixedStabilizer, MixedDestabilizer,
-    prodphase, comm,
+    prodphase, comm, comm!,
     nqubits,
     stabilizerview, destabilizerview, logicalxview, logicalzview, phases,
     fastcolumn, fastrow,
@@ -61,7 +61,8 @@ export
     enumerate_single_qubit_gates, random_clifford1,
     enumerate_cliffords, symplecticGS, clifford_cardinality, enumerate_phases,
     random_invertible_gf2,
-    random_pauli, random_stabilizer, random_destabilizer, random_clifford,
+    random_pauli, random_pauli!,
+    random_stabilizer, random_destabilizer, random_clifford,
     # Noise
     applynoise!, UnbiasedUncorrelatedNoise, NoiseOp, NoiseOpAll, NoisyGate,
     PauliNoise, PauliError,
@@ -701,7 +702,11 @@ julia> comm(P"ZZ", P"XX")
 julia> comm(P"IZ", P"XX")
 0x01
 ```
+
+See also: [`comm!`](@ref)
 """
+function comm end
+
 @inline function comm(l::AbstractVector{T}, r::AbstractVector{T})::UInt8 where T<:Unsigned
     res = T(0)
     len = length(l)÷2
@@ -739,6 +744,24 @@ comm(l::Tableau, r::PauliOperator) = comm(r, l)
 @inline comm(l::Stabilizer, r::PauliOperator) = comm(tab(l), r)
 @inline comm(s::Stabilizer, l::Int, r::Int) = comm(tab(s), l, r)
 
+"""An in-place version of [`comm`](@ref), storing its output in the given buffer."""
+function comm! end
+function comm!(v, l::PauliOperator, r::Tableau)
+    length(v) == length(r) || throw(DimensionMismatch(lazy"The dimensions of the output buffer and the input tableau have to match in `comm!`"))
+    nqubits(l) == nqubits(r) || throw(DimensionMismatch(lazy"The number of qubits of the input Pauli operator and the input tableau have to match in `comm!`"))
+    for i in 1:length(r)
+        v[i] = comm(l,r,i)
+    end
+    v
+end
+comm!(v, l::Tableau, r::PauliOperator) = comm!(v, r, l)
+@inline comm!(v, l::PauliOperator, r::Stabilizer, i::Int) = comm!(v, l, tab(r), i)
+@inline comm!(v, l::Stabilizer, r::PauliOperator, i::Int) = comm!(v, tab(l), r, i)
+@inline comm!(v, l::PauliOperator, r::Stabilizer) = comm!(v, l, tab(r))
+@inline comm!(v, l::Stabilizer, r::PauliOperator) = comm!(v, tab(l), r)
+@inline comm!(v, s::Stabilizer, l::Int, r::Int) = comm!(v, tab(s), l, r)
+
+
 Base.:(*)(l::PauliOperator, r::PauliOperator) = mul_left!(copy(r),l)
 
 (⊗)(l::PauliOperator, r::PauliOperator) = PauliOperator((l.phase[]+r.phase[])&0x3, vcat(xbit(l),xbit(r)), vcat(zbit(l),zbit(r)))

src/ecc/ECC.jl

@@ -3,7 +3,7 @@ module ECC
 using LinearAlgebra
 using QuantumClifford
 using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer
-import QuantumClifford: Stabilizer, MixedDestabilizer
+import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits
 using DocStringExtensions
 using Combinatorics: combinations
 using SparseArrays: sparse
@@ -56,6 +56,8 @@ MixedDestabilizer(c::AbstractECC; kwarg...) = MixedDestabilizer(Stabilizer(c); k
 """The number of physical qubits in a code."""
 function code_n end
 
+nqubits(c::AbstractECC) = code_n(c::AbstractECC)
+
 code_n(c::AbstractECC) = code_n(parity_checks(c))
 
 code_n(s::Stabilizer) = nqubits(s)

src/ecc/decoder_pipeline.jl

@@ -6,6 +6,24 @@ All `AbstractSyndromeDecoder` types are expected to:
 - have an `evaluate_decoder` method which runs a full simulation but it supports only a small number of ECC protocols"""
 abstract type AbstractSyndromeDecoder end
 
+"""Decode a syndrome using a given decoding algorithm."""
+function decode end
+
+"""Decode a batch of syndromes using a given decoding algorithm."""
+function batchdecode(d::AbstractSyndromeDecoder, syndrome_samples)
+    H = parity_checks(d)
+    s, n = size(H)
+    samples, _s = size(syndrome_samples)
+    s == _s || throw(ArgumentError(lazy"The syndromes given to `batchdecode` have the wrong dimensions. The syndrome length is $(_s) while it should be $(s)"))
+    results = falses(samples, 2n)
+    for (i,syndrome_sample) in enumerate(eachrow(syndrome_samples))
+        guess = decode(d, syndrome_sample)# TODO use `decode!`
+        isnothing(guess) || (results[i,:] = guess)
+    end
+    return results
+end
+
+
 """An abstract type mostly used by [`evaluate_decoder`](@ref) to specify in what context to evaluate an ECC."""
 abstract type AbstractECCSetup end
 
@@ -119,37 +137,37 @@ function evaluate_decoder(d::AbstractSyndromeDecoder, nsamples, circuit, syndrom
 
     syndromes = @view pfmeasurements(frames)[:, syndrome_bits]
     measured_faults = @view pfmeasurements(frames)[:, logical_bits]
+    guesses = batchdecode(d, syndromes)
+    evaluate_guesses(measured_faults, guesses, faults_submatrix)
+end
+
+function evaluate_guesses(measured_faults, guesses, faults_matrix)
+    nsamples = size(guesses, 1)
+    guess_faults = (faults_matrix * guesses') .% 2 # TODO this can be faster and non-allocating by turning it into a loop
     decoded = 0
-    for i in 1:nsamples
-        guess = decode(d, @view syndromes[i,:])
-        isnothing(guess) && continue
-        guess_faults = faults_submatrix * guess .% 2
-        if guess_faults == @view measured_faults[i,:]
-            decoded += 1
-        end
+    for i in 1:nsamples # TODO this can be faster and non-allocating by having the loop and the matrix multiplication on the line above work together and not store anything
+        (@view guess_faults[:,i]) == (@view measured_faults[i,:]) && (decoded += 1)
     end
-
     return (nsamples - decoded) / nsamples
 end
 
 function evaluate_decoder(d::AbstractSyndromeDecoder, setup::CommutationCheckECCSetup, nsamples::Int)
     H = parity_checks(d)
     fm = faults_matrix(H)
+    fmtab = QuantumClifford.Tableau(fm[:,end÷2+1:end],fm[:,1:end÷2]) # TODO there should be a special method for this
+    k, s = size(fm)
     n = nqubits(H)
-    decoded = 0
-    for i in 1:nsamples # TODO fix all this casting and allocation
-        err = random_pauli(n, setup.xz_noise, nophase=true)
-        syndrome = Bool.(comm(H,err))
-        guess = decode(d, syndrome)
-        isnothing(guess) && continue
-        guess_faults = fm * guess .% 2
-        measured_faults = fm * stab_to_gf2(err) .% 2
-        if guess_faults == measured_faults
-            decoded += 1
-        end
+    err = zero(PauliOperator, n)
+    syndromes = zeros(Bool, nsamples, length(H)) # TODO will this be better and faster if we use bitmatrices
+    measured_faults = zeros(UInt8, nsamples, k)
+    for i in 1:nsamples
+        err = random_pauli!(err, setup.xz_noise, nophase=true)
+        comm!((@view syndromes[i,:]), H,err)
+        comm!((@view measured_faults[i,:]), fmtab,err)
     end
-
-    return (nsamples - decoded) / nsamples
+    measured_faults .%= 2
+    guesses = batchdecode(d, syndromes)
+    evaluate_guesses(measured_faults, guesses, fm)
 end
 
 """A simple look-up table decoder for error correcting codes.
@@ -164,13 +182,15 @@ struct TableDecoder <: AbstractSyndromeDecoder
     """Faults matrix corresponding to the code"""
     faults_matrix
     """The number of qubits in the code"""
-    n
+    n::Int
     """The depth of the code"""
-    s
+    s::Int
     """The number of encoded qubits"""
-    k
+    k::Int
     """The lookup table corresponding to the code, slow to create"""
-    lookup_table
+    lookup_table::Dict{Vector{Bool},Vector{Bool}}
+    lookup_buffer::Vector{Bool}
+    TableDecoder(H, faults_matrix, n, s, k, lookup_table) = new(H, faults_matrix, n, s, k, lookup_table, fill(false, s))
 end
 
 function TableDecoder(c)
@@ -180,13 +200,13 @@ function TableDecoder(c)
     k = n - r
     lookup_table = create_lookup_table(H)
     fm = faults_matrix(H)
-    return TableDecoder(H, n, s, k, fm, lookup_table)
+    return TableDecoder(H, fm, n, s, k, lookup_table)
 end
 
 parity_checks(d::TableDecoder) = d.H
 
 function create_lookup_table(code::Stabilizer)
-    lookup_table = Dict()
+    lookup_table = Dict{Vector{Bool},Vector{Bool}}()
     constraints, qubits = size(code)
     # In the case of no errors
     lookup_table[ zeros(UInt8, constraints) ] = stab_to_gf2(zero(PauliOperator, qubits))
@@ -206,7 +226,8 @@ function create_lookup_table(code::Stabilizer)
 end;
 
 function decode(d::TableDecoder, syndrome_sample)
-    return get(d.lookup_table, syndrome_sample, nothing)
+    d.lookup_buffer .= syndrome_sample # TODO have this work without data copying, by supporting the correct types, especially in the batch decode case
+    return get(d.lookup_table, d.lookup_buffer, nothing)
 end
 
 # From extensions:

src/noise.jl

@@ -22,23 +22,24 @@ function applynoise!(r::Register, n, indices::Base.AbstractVecOrTuple)
     return r
 end
 
-"""Depolarization noise model with total probability of error `3*errprobthird`."""
+"""Depolarization noise model with total probability of error `p`."""
 struct UnbiasedUncorrelatedNoise{T} <: AbstractNoise
-    errprobthird::T
+    p::T
 end
+UnbiasedUncorrelatedNoise(p::Integer) = UnbiasedUncorrelatedNoise(float(p))
 
 """A convenient constructor for various types of Pauli noise models.
 Returns more specific types when necessary."""
 function PauliNoise end
 
 """Constructs an unbiased Pauli noise model with total probability of error `p`."""
 function PauliNoise(p)
-    UnbiasedUncorrelatedNoise(p/3)
+    UnbiasedUncorrelatedNoise(p)
 end
 
 function applynoise!(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,i::Int)
     n = nqubits(s)
-    infid = noise.errprobthird
+    infid = noise.p/3
     r = rand()
     if r<infid
         apply_single_x!(s,i)
@@ -121,11 +122,11 @@ end
 function applynoise_branches(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,indices; max_order=1)
     n = nqubits(s)
     l = length(indices)
-    infid = noise.errprobthird
+    infid = noise.p/3
     if l==0
         return [s,one(infid)]
     end
-    error1 = 3*infid
+    error1 = noise.p
     no_error1 = 1-error1
     no_error = no_error1^l
     results = [(copy(s),no_error,0)] # state, prob, order

src/pauli_frames.jl

@@ -120,13 +120,14 @@ function apply!(frame::PauliFrame, op::sMRZ) # TODO sMRY, and faster sMRX
 end
 
 function applynoise!(frame::PauliFrame,noise::UnbiasedUncorrelatedNoise,i::Int)
-    p = noise.errprobthird
+    p = noise.p
     T = eltype(frame.frame.tab.xzs)
 
     lowbit = T(1)
     ibig = _div(T,i-1)+1
     ismall = _mod(T,i-1)
     ismallm = lowbit<<(ismall)
+    p = p/3
 
     @inbounds @simd for f in eachindex(frame)
         r = rand()