performance improvements

ext/QuantumCliffordJuMPExt/QuantumCliffordJuMPExt.jl

@@ -10,9 +10,10 @@ import GLPK: Optimizer
 
 import QuantumClifford
 import QuantumClifford: stab_to_gf2, logicalxview, logicalzview, canonicalize!,
-    MixedDestabilizer, Stabilizer
+    MixedDestabilizer, Stabilizer, gf2_row_echelon_with_pivots!
 import QuantumClifford.ECC
-import QuantumClifford.ECC: distance, AbstractECC, code_n, code_k, parity_checks
+import QuantumClifford.ECC: distance, AbstractECC, code_n, code_k, parity_checks,
+    parity_checks_x, parity_checks_z
 
 import SparseArrays
 import SparseArrays: SparseMatrixCSC, sparse, spzeros, findnz, sparsevec

ext/QuantumCliffordJuMPExt/min_distance_mixed_integer_programming.jl

@@ -147,17 +147,17 @@ to a mixed integer linear program and using the GNU Linear Programming Kit.
 the code distance was calculated using the mixed integer programming approach.
 
 """
-function distance(c::Stabilizer; upper_bound=false, logical_qubit=code_k(c), all_logical_qubits=false, logical_operator_type=:X)
+function distance(c::AbstractECC; upper_bound=false, logical_qubit=code_k(c), all_logical_qubits=false, logical_operator_type=:X)
     1 <= logical_qubit <= code_k(c) || throw(ArgumentError("The number of logical qubit must be between 1 and $(code_k(c)) inclusive"))
     logical_operator_type == :X || logical_operator_type == :Z || throw(ArgumentError("Invalid type of logical operator: Use :X or :Z"))
-    l = get_lx_lz(c)[1]
-    H = SparseMatrixCSC{Int, Int}(stab_to_gf2(parity_checks(c)))
-    h = get_stab(H, :X)
+    hx = Matrix{Int}(parity_checks_x(c))
+    hz = Matrix{Int}(parity_checks_z(c))
+    l = get_logicals(hz, hx)
     if logical_operator_type == :Z
-        l = get_lx_lz(c)[2]
-        H = SparseMatrixCSC{Int, Int}(stab_to_gf2(parity_checks(c)))
-        h = get_stab(H, :Z)
+        l = get_logicals(hx, hz)
+        h = SparseMatrixCSC{Int, Int}(hz)
     end
+    h = SparseMatrixCSC{Int, Int}(hx)
     weights = Dict{Int, Int}()
     for i in 1:logical_qubit
         w = _minimum_distance(h, l[i, :])

ext/QuantumCliffordJuMPExt/util.jl

@@ -1,20 +1,24 @@
-function get_stab(matrix::SparseMatrixCSC{Int, Int}, logical_operator_type::Symbol)
-    rows, cols = size(matrix)
-    if logical_operator_type == :Z
-        col_range = 1:div(cols, 2)
-    else logical_operator_type == :X
-        col_range = div(cols, 2) + 1:cols
-    end
-    submatrix = matrix[:, col_range]
-    non_zero_rows = unique(submatrix.rowval)
-    zero_rows = setdiff(1:rows, non_zero_rows)
-    return matrix[zero_rows, :]
-end
+"""
+A logical X-type operator \$L_x\$ satisfies the
+following two conditions:
+1. \$L_x \\in \\ker(H_z)\$
+2. \$L_x \\notin \\text{Im}(H_x^T)\$
 
-function get_lx_lz(c::Stabilizer)
-    lx = stab_to_gf2(logicalxview(canonicalize!(MixedDestabilizer(c))))
-    lz = stab_to_gf2(logicalzview(canonicalize!(MixedDestabilizer(c))))
-    lx = SparseMatrixCSC{Int, Int}(lx)
-    lz = SparseMatrixCSC{Int, Int}(lz)
-    return lx, lz
+A logical Z-type operator \$L_z\$ satisfies the
+following two conditions:
+1. \$L_z \\in \\ker(H_x)\$
+2. \$L_z \\notin \\text{Im}(H_z^T)\$
+"""
+function get_logicals(hx, hz)
+    m = size(copy(hx)',1)
+    ker_hx = (gf2_row_echelon_with_pivots!(copy(hx)')[3])[gf2_row_echelon_with_pivots!(copy(hx)')[2]+1:m, :] # nullspace
+    im_hzT = copy(hz)[gf2_row_echelon_with_pivots!(copy(hz)')[4],:] # row-basis
+    log_stack = vcat(im_hzT, ker_hx)
+    log_stackT = log_stack'
+    pivots = gf2_row_echelon_with_pivots!(copy(log_stackT))[4]
+    r1 = size(copy(im_hzT),1)
+    r2 = size(copy(log_stack),1)
+    log_op_indices = [i for i in r1+1:r2 if i in copy(pivots)]
+    log_ops = copy(log_stack)[log_op_indices, :]
+    SparseMatrixCSC{Int, Int}(log_ops)
 end

src/QuantumClifford.jl

@@ -31,6 +31,7 @@ export
     affectedqubits, #TODO move to QuantumInterface?
     # GF2
     stab_to_gf2, gf2_gausselim!, gf2_isinvertible, gf2_invert, gf2_H_to_G,
+    gf2_row_echelon_with_pivots!,
     # Canonicalization
     canonicalize!, canonicalize_rref!, canonicalize_gott!,
     # Linear Algebra
@@ -1168,6 +1169,50 @@ function gf2_H_to_G(H)
     G[:,invperm(sindx)]
 end
 
+"""Performs in-place Gaussian elimination on a binary matrix and returns
+its *row echelon form*,*rank*, the *transformation matrix*, and the *pivot
+columns*. The transformation matrix that converts the original matrix into
+the row echelon form. The `full` parameter controls the extent of elimination:
+if `true`, only rows below the pivot are affected; if `false`, both above and
+below the pivot are eliminated."""
+function gf2_row_echelon_with_pivots!(M::AbstractMatrix{Int}; full=false)
+    r, c = size(M)
+    N = Matrix{Int}(LinearAlgebra.I, r, r)
+    p = 1
+    pivots = Int[]
+    for col in 1:c
+        @inbounds for row in p:r
+            if M[row, col] == 1
+                if row != p
+                    M[[row, p], :] .= M[[p, row], :]
+                    N[[row, p], :] .= N[[p, row], :]
+                end
+                break
+            end
+        end
+        if M[p, col] == 1
+            if !full
+                elim_range = p+1:r
+            else
+                elim_range = 1:r
+            end
+            @simd for j in elim_range
+                @inbounds if j != p && M[j, col] == 1
+                    M[j, :] .= (M[j, :] .+ M[p, :]) .% 2
+                    N[j, :] .= (N[j, :] .+ N[p, :]) .% 2
+                end
+            end
+            p += 1
+            push!(pivots, col)
+        end
+        if p > r
+            break
+        end
+    end
+    rank = p - 1
+    return M, rank, N, pivots
+end
+
 ##############################
 # Error classes
 ##############################

src/ecc/ECC.jl

@@ -128,8 +128,6 @@ The minimum distance of a qLDPC code requires:
 """
 function distance end
 
-distance(c::AbstractECC; kwargs...) = distance(parity_checks(c); kwargs...)
-
 """Parity matrix of a code, given as a stabilizer tableau."""
 function parity_matrix(c::AbstractECC)
     paritym = stab_to_gf2(parity_checks(c::AbstractECC))