Some group theory tools for Quantum error correction (#293)

Co-authored-by: Kenneth Goodenough <kdgoodenough@gmail.com>
Co-authored-by: Stefan Krastanov <github.acc@krastanov.org>
Co-authored-by: Stefan Krastanov <stefan@krastanov.org>

CHANGELOG.md

@@ -5,8 +5,16 @@
 
 # News
 
-## v0.9.8 - 2024-07-24
-
+## v0.9.8 - 2024-08-03
+
+- New group-theoretical tools:
+ - `groupify` to generate full stabilizer group from generating set
+ - `minimal_generating_set` function to find the minimal generating set of a set of operators
+ - `pauligroup` to generate the full Pauli group of a certain number of qubits
+ - `normalizer` to generate all Paulis that commute with a set of Paulis
+ - `centralizer` to find a subset of a set of Paulis such that each element in the subset commutes with each element in the set
+ - `contractor` to find a subset of Paulis in a tableau that have an identity operator on a certain qubit
+ - `delete_columns` to remove the operators corresponding to a certain qubit from all Paulis in a Stabilizer
 - `PauliError` can now encode biased noise during Pauli frame simulation, i.e. one can simulate only X errors, or only Y errors, or only Z errors, or some weighted combination of these.
 
 ## v0.9.7 - 2024-07-23

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org> and QuantumSavory community members"]
-version = "0.9.7"
+version = "0.9.8"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/QuantumClifford.jl

@@ -74,6 +74,8 @@ export
  single_z, single_x, single_y,
  # Graphs
  graphstate, graphstate!, graph_gatesequence, graph_gate,
+ # Group theory tools
+ groupify, minimal_generating_set, pauligroup, normalizer, centralizer, contractor, delete_columns,
  # Clipped Gauge
  canonicalize_clip!, bigram, entanglement_entropy,
  # mctrajectories
@@ -1202,6 +1204,7 @@ include("sumtypes.jl")
 include("precompiles.jl")
 include("ecc/ECC.jl")
 include("nonclifford.jl")
+include("grouptableaux.jl")
 include("plotting_extensions.jl")
 #
 include("gpu_adapters.jl")

src/grouptableaux.jl

@@ -0,0 +1,211 @@
+using Graphs
+using LinearAlgebra
+
+"""
+Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)).
+
+The returned object is exponentially long.
+
+```jldoctest
+julia> groupify(S"XZ ZX")
++ __
++ XZ
++ ZX
++ YY
+```
+"""
+function groupify(s::Stabilizer)
+ # Create a `Tableau` of 2ⁿ n-qubit identity Pauli operators(where n is the size of 
+ # `Stabilizer` s), then multiply each one by a different subset of the elements in s to 
+ # create all 2ⁿ unique elements in the group generated by s, then return the `Tableau`.
+ n = length(s)::Int
+ group = zero(Tableau, 2^n, nqubits(s))
+ for i in 0:2^n-1
+ for (digit_order, j) in enumerate(digits(i, base=2, pad=n))
+ if j == 1
+ group[i+1] *= s[digit_order]
+ end
+ end
+ end
+ return group 
+end
+
+
+"""
+For a not-necessarily-minimal generating set,
+return the minimal generating set.
+
+The input has to have only real phases.
+
+```jldoctest
+julia> minimal_generating_set(S"__ XZ ZX YY")
++ XZ
++ ZX
+```
+"""
+function minimal_generating_set(s::Stabilizer)
+ # Canonicalize `Stabilizer` s, then return a `Stabilizer` with all non-identity Pauli operators 
+ # in the result. If s consists of only identity operators, return the negative 
+ # identity operator if one is contained in s, and the positive identity operator otherwise.
+ s, _, r = canonicalize!(copy(s), ranks=true)
+ if r == 0
+ gs = zero(Stabilizer, 1, nqubits(s))
+ if 0x02 in phases(s)
+ gs[1] = -1 * gs[1]
+ end
+ return gs
+ else
+ return s[1:r, :]
+ end
+end
+
+"""
+Return the full Pauli group of a given length. Phases are ignored by default, 
+but can be included by setting `phases=true`.
+
+```jldoctest
+julia> pauligroup(1)
++ _
++ X
++ Z
++ Y
+
+julia> pauligroup(1, phases=true)
++ _
++ X
++ Z
++ Y
+- _
+- X
+- Z
+- Y
++i_
++iX
++iZ
++iY
+-i_
+-iX
+-iZ
+-iY
+```
+"""
+function pauligroup(n::Int; phases=false)
+ if phases
+ s = zero(Tableau, 4^(n + 1), n)
+ paulis = ((false, false), (true, false), (false, true), (true, true))
+ for (i, P) in enumerate(Iterators.product(Iterators.repeated(paulis, n)...))
+ for (j, p) in enumerate(P)
+ s[i, j] = p
+ end
+ end
+ for i in 1:4^n
+ s[i+4^n] = -1 * s[i]
+ end
+ for i in 4^n+1:2*4^n
+ s[i+4^n] = -1im * s[i]
+ end
+ for i in 2*4^n+1:3*4^n
+ s[i+4^n] = -1 * s[i]
+ end
+ end
+ if !phases
+ s = zero(Tableau, 4^n, n)
+ paulis = ((false, false), (true, false), (false, true), (true, true))
+ for (i, P) in enumerate(Iterators.product(Iterators.repeated(paulis, n)...))
+ for (j, p) in enumerate(P)
+ s[i, j] = p
+ end
+ end
+ end
+ return s
+end
+
+"""
+Return all Pauli operators with the same number of qubits as the given `Tableau` `t`
+that commute with all operators in `t`.
+
+```jldoctest
+julia> normalizer(T"X")
++ _
++ X
+```
+"""
+function normalizer(t::Tableau)
+ # For each `PauliOperator` p in the with same number of qubits as the `Stabilizer` s, iterate through s and check each 
+ # operator's commutivity with p. If they all commute, add p a vector of `PauliOperators`. Return the vector 
+ # converted to `Tableau`.
+ n = nqubits(t)
+ pgroup = pauligroup(n, phases=false)
+ ptype = typeof(t[1])
+ normalizer = ptype[]
+ for p in pgroup
+ commutes = true
+ for q in t
+ if comm(p, q) == 0x01
+ commutes = false
+ end
+ end
+ if commutes
+ push!(normalizer, p)
+ end
+ end
+
+ return Tableau(normalizer)
+end
+
+"""
+For a given set of Paulis (in the form of a `Tableau`), return the subset of Paulis that commute with all Paulis in set.
+
+```jldoctest
+julia> centralizer(T"XX ZZ _Z")
++ ZZ
+```
+"""
+function centralizer(t::Tableau) 
+ center = typeof(t[1])[]
+ for P in t
+ commutes = 0
+ for Q in t
+ if comm(P, Q) == 0x01
+ commutes = 1
+ break
+ end
+ end
+ if commutes == 0
+ push!(center, P)
+ end
+ end
+ if length(center) == 0 
+ return Tableau(zeros(Bool, 1,1))
+ end
+ c = Tableau(center)
+ return c
+end
+
+"""
+Return the subset of Paulis in a Stabilizer that have identity operators on all qubits corresponding to 
+the given subset, without the entries corresponding to subset.
+
+```jldoctest
+julia> contractor(S"_X X_", [1])
++ X
+```
+"""
+function contractor(s::Stabilizer, subset)
+ result = typeof(s[1])[]
+ for p in s
+ contractable = true
+ for i in subset
+ if p[i] != (false, false) 
+ contractable = false 
+ break 
+ end
+ end
+ if contractable push!(result, p[setdiff(1:length(p), subset)]) end
+ end
+ if length(result) > 0
+ return Tableau(result)
+ else
+ return Tableau(zeros(Bool, 1,1))
+ end
+end

src/pauli_operator.jl

@@ -122,6 +122,11 @@ Base.hash(p::PauliOperator, h::UInt) = hash(p.phase,hash(p.nqubits,hash(p.xz, h)
 
 Base.copy(p::PauliOperator) = PauliOperator(copy(p.phase),p.nqubits,copy(p.xz))
 
+function Base.deleteat!(p::PauliOperator, subset) 
+ p =p[setdiff(1:length(p), subset)]
+ return p
+end
+
 _nchunks(i::Int,T::Type{<:Unsigned}) = 2*( (i-1) ÷ (8*sizeof(T)) + 1 )
 Base.zero(::Type{PauliOperator{Tₚ, Tᵥ}}, q) where {Tₚ,T<:Unsigned,Tᵥ<:AbstractVector{T}} = PauliOperator(zeros(UInt8), q, zeros(T, _nchunks(q,T)))
 Base.zero(::Type{PauliOperator}, q) = zero(PauliOperator{Array{UInt8, 0}, Vector{UInt}}, q)

src/project_trace_reset.jl

@@ -627,6 +627,8 @@ end
 $TYPEDSIGNATURES
 
 Trace out a qubit.
+
+See also: [`delete_columns`](@ref)
 """ # TODO all of these should raise an error if length(qubits)>rank
 function traceout!(s::Stabilizer, qubits; phases=true, rank=false)
  _,i = canonicalize_rref!(s,qubits;phases=phases)
@@ -866,3 +868,22 @@ function traceoutremove!(s::MixedDestabilizer, qubit)
  traceout!(s,[qubit]) # TODO this can be optimized thanks to the information already known from projfunc
  s = _remove_rowcol!(s, nqubits(s), qubit)
 end
+
+
+"""
+Return the given stabilizer without all the qubits in the given iterable.
+
+The resulting tableaux is not guaranteed to be valid (to retain its commutation relationships).
+
+```jldoctest
+julia> delete_columns(S"XYZ YZX ZXY", [1,3])
++ Y
++ Z
++ X
+```
+
+See also: [`traceout!`](@ref)
+"""
+function delete_columns(𝒮::Stabilizer, subset)
+ return 𝒮[:, setdiff(1:nqubits(𝒮), subset)]
+end