Skip to content

Commit

Permalink
QuantumCliffordJuMPExt: compute the minimum distance of QLDPC using M…
Browse files Browse the repository at this point in the history
…ixed Integer Programming (MIP)
  • Loading branch information
Fe-r-oz committed Nov 30, 2024
1 parent 5fcf353 commit 71b5f9e
Show file tree
Hide file tree
Showing 9 changed files with 234 additions and 11 deletions.
5 changes: 5 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,9 @@ SumTypes = "8e1ec7a9-0e02-4297-b0fe-6433085c89f2"
[weakdeps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
LDPCDecoders = "3c486d74-64b9-4c60-8b1a-13a564e77efb"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PyQDecoders = "17f5de1a-9b79-4409-a58d-4d45812840f7"
Expand All @@ -35,6 +37,7 @@ QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
[extensions]
QuantumCliffordGPUExt = "CUDA"
QuantumCliffordHeckeExt = "Hecke"
QuantumCliffordJuMPExt = ["GLPK", "JuMP"]
QuantumCliffordLDPCDecodersExt = "LDPCDecoders"
QuantumCliffordMakieExt = "Makie"
QuantumCliffordPlotsExt = "Plots"
Expand All @@ -48,10 +51,12 @@ Combinatorics = "1.0"
DataStructures = "0.18"
DocStringExtensions = "0.9"
Graphs = "1.9"
GLPK = "1.2"
Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34"
HostCPUFeatures = "0.1.6"
ILog2 = "0.2.3, 1, 2"
InteractiveUtils = "1.9"
JuMP = "1.23"
LDPCDecoders = "0.3.1"
LinearAlgebra = "1.9"
MacroTools = "0.5.9"
Expand Down
2 changes: 2 additions & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,12 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21"
LDPCDecoders = "3c486d74-64b9-4c60-8b1a-13a564e77efb"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PyQDecoders = "17f5de1a-9b79-4409-a58d-4d45812840f7"
Expand Down
7 changes: 7 additions & 0 deletions docs/src/references.bib
Original file line number Diff line number Diff line change
Expand Up @@ -561,3 +561,10 @@ @article{haah2011local
pages={042330},
year={2011},
}

@article{landahl2011fault,
title={Fault-tolerant quantum computing with color codes},
author={Landahl, Andrew J and Anderson, Jonas T and Rice, Patrick R},
journal={arXiv preprint arXiv:1108.5738},
year={2011}
}
2 changes: 1 addition & 1 deletion ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ import Nemo: characteristic, matrix_repr, GF, ZZ, lift
import QuantumClifford.ECC: AbstractECC, CSS, ClassicalCode,
hgp, code_k, code_n, code_s, iscss, parity_checks, parity_checks_x, parity_checks_z, parity_checks_xz,
two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes, check_repr_commutation_relation,
haah_cubic_codes
haah_cubic_codes, minimum_distance

include("util.jl")
include("types.jl")
Expand Down
32 changes: 23 additions & 9 deletions ext/QuantumCliffordHeckeExt/lifted_product.jl
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,9 @@ Here is an example of a [[56, 28, 2]] 2BGA code from Table 2 of [lin2024quantum]
with direct product of `C₄ x C₂`.
```jldoctest
julia> import Hecke: group_algebra, GF, abelian_group, gens
julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes; # hide
julia> GA = group_algebra(GF(2), abelian_group([14,2]));
Expand All @@ -174,8 +176,8 @@ julia> B = 1 + x^7 + s + x^8 + s*x^7 + x
julia> c = two_block_group_algebra_codes(A,B);
julia> code_n(c), code_k(c)
(56, 28)
julia> code_n(c), code_k(c), minimum_distance(c)
(56, 28, 2)
```
### Bivariate Bicycle codes
Expand All @@ -189,7 +191,7 @@ The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/qcga
A [[756, 16, ≤ 34]] code from Table 3 of [bravyi2024high](@cite):
```jldoctest
julia> import Hecke: group_algebra, GF, abelian_group, gens
julia> import Hecke: group_algebra, GF, abelian_group, gens;
julia> l=21; m=18;
Expand All @@ -215,10 +217,16 @@ where `𝐺ᵣ = ℤ/l₁ × ℤ/l₂ × ... × ℤ/lᵣ`.
A [[48, 4, 6]] Weight-6 TB-QLDPC code from Appendix A Table 2 of [voss2024multivariatebicyclecodes](@cite).
```jldoctest
julia> import Hecke: group_algebra, GF, abelian_group, gens; # hide
julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes; # hide
julia> l=4; m=6;
julia> GA = group_algebra(GF(2), abelian_group([l, m]));
julia> x, y = gens(GA);
julia> z = x*y;
julia> A = x^3 + y^5;
Expand All @@ -227,8 +235,8 @@ julia> B = x + z^5 + y^5 + y^2;
julia> c = two_block_group_algebra_codes(A, B);
julia> code_n(c), code_k(c)
(48, 4)
julia> code_n(c), code_k(c), minimum_distance(c)
(48, 4, 6)
```
### Coprime Bivariate Bicycle code
Expand All @@ -241,7 +249,9 @@ based on abelian group `ℤₗ x ℤₘ` where `ℤⱼ` cyclic group of order `j
[108, 12, 6]] coprime-bivariate bicycle (BB) code from Table 2 of [wang2024coprime](@cite).
```jldoctest
julia> import Hecke: group_algebra, GF, abelian_group, gens;
julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes; # hide
julia> l=2; m=27;
Expand All @@ -252,7 +262,11 @@ julia> 𝜋 = gens(GA)[1];
julia> A = 𝜋^2 + 𝜋^5 + 𝜋^44;
julia> B = 𝜋^8 + 𝜋^14 + 𝜋^47;
(108, 12)
julia> c = two_block_group_algebra_codes(A, B);
julia> code_n(c), code_k(c), minimum_distance(c)
(108, 12, 6)
```
See also: [`LPCode`](@ref), [`generalized_bicycle_codes`](@ref), [`bicycle_codes`](@ref), [`haah_cubic_codes`](@ref).
Expand Down
24 changes: 24 additions & 0 deletions ext/QuantumCliffordJuMPExt/QuantumCliffordJuMPExt.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
module QuantumCliffordJuMPExt

using DocStringExtensions

import JuMP
import JuMP: @variable, @objective, @constraint, optimize!, Model, set_silent,
sum, value
import GLPK
import GLPK: Optimizer

import QuantumClifford
import QuantumClifford: stab_to_gf2, logicalxview, logicalzview, canonicalize!,
MixedDestabilizer, Stabilizer
import QuantumClifford.ECC
import QuantumClifford.ECC: minimum_distance, AbstractECC, code_n, code_k,
parity_checks

import SparseArrays
import SparseArrays: SparseMatrixCSC, sparse

include("util.jl")
include("min_distance_mixed_integer_programming.jl")

end
147 changes: 147 additions & 0 deletions ext/QuantumCliffordJuMPExt/min_distance_mixed_integer_programming.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,147 @@
"""
TYPEDSIGNATURES
Computes the minimum Hamming weight of a binary vector `x` by solving an **mixed
integer program (MIP)** that satisfies the following constraints:
- \$\\text{stab} \\cdot x \\equiv 0 \\pmod{2}\$: The binary vector x must have an
even overlap with each X-check of the stabilizer binary representation `stab`.
- \$\\text{logicOp} \\cdot x \\equiv 1 \\pmod{2}\$: The binary vector x must have
an odd overlap with logical-X operator `logicOp` on the i-th logical qubit.
It computes the minimum Hamming weight \$d_{Z}\$ for the Z-type logical
operator. The distance for X-type logical operators is the same.
### Mixed Integer Programming
The MIP minimizes the Hamming weight `w(x)`, defined as the number of nonzero
elements in `x`, while satisfying the constraints:
\$\$
\\begin{aligned}
\\text{Minimize} \\quad & w(x) = \\sum_{i=1}^n x_i, \\\\
\\text{subject to} \\quad & \\text{stab} \\cdot x \\equiv 0 \\pmod{2}, \\\\
& \\text{logicOp} \\cdot x \\equiv 1 \\pmod{2}, \\\\
& x_i \\in \\{0, 1\\} \\quad \\text{for all } i.
\\end{aligned}
\$\$
Here:
- \$\\text{stab}\$ is the binary matrix representing the stabilizer group.
- \$\\text{logicOp}\$ is the binary vector representing the logical-X operator.
- `x` is the binary vector (decision variable) being optimized.
The optimal solution \$w_{i}\$ for each logical-X operator corresponds to the
minimum weight of a Pauli Z-type operator satisfying the above conditions.
The Z-type distance is given by:
\$\$
\\begin{aligned}
d_Z = \\min(w_1, w_2, \\dots, w_k),
\\end{aligned}
\$\$
where k is the number of logical qubits.
# Example
A [[40, 8, 5]] 2BGA code with the minimum distance of 5.
```jldoctest
julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes;
julia> l = 10; m = 2;
julia> GA = group_algebra(GF(2), abelian_group([l,m]));
julia> x, s = gens(GA);
julia> A = 1 + x^6;
julia> B = 1 + x^5 + s + x^6 + x + s*x^2;
julia> c = two_block_group_algebra_codes(A,B);
julia> minimum_distance(c)
5
```
### Applications
- The first usecase of the MIP approach was the code capacity Most Likely
Error (MLE) decoder for color codes introduced in [landahl2011fault](@cite).
- For all quantum LDPC codes presented in [panteleev2021degenerate](@cite),
the lower and upper bounds on the minimum distance was obtained by reduction
to a mixed integer linear program and using the GNU Linear Programming Kit.
- For all the Bivariate Bicycle (BB) codes presented in [bravyi2024high](@cite),
the code distance was calculated using the mixed integer programming approach.
"""
function minimum_distance(c::Stabilizer; num_logical_qubits=code_k(c))
lx, _ = get_lx_lz(c)
H = stab_to_gf2(parity_checks(c))
H = sparse(H)
hx = Matrix{Int}(get_stab_hx(H))
num_logical_qubits == 1 && return _minimum_distance(hx, lx[num_logical_qubits, :])
if 2 <= num_logical_qubits && num_logical_qubits <= code_k(c)
w_values = []
for i in 1:num_logical_qubits
w = _minimum_distance(hx, lx[num_logical_qubits, :])
push!(w_values, w)
end
return round(Int, sum(w_values)/num_logical_qubits)
else
throw(ArgumentError("The number of logical qubits must be between 1 to $(code_k(c)) inclusive"))
end
end

# Computing minimum distance for quantum LDPC codes is a NP-Hard problem,
# but we can solve mixed integer program (MIP) for small instance sizes.
# inspired from [bravyi2024high](@cite) and [landahl2011fault](@cite).
function _minimum_distance(hx, lx)
n = size(hx, 2) # number of qubits
m = size(hx, 1) # number of stabilizers
# Maximum stabilizer weight
whx = maximum(sum(hx[i, :] for i in 1:m)) # Maximum row sum of hx
# Weight of the logical operator
wlx = count(!iszero, lx) # Count non-zero elements in lx
# Number of slack variables needed to express orthogonality constraints modulo two
num_anc_hx = ceil(Int, log2(whx))
num_anc_logical = ceil(Int, log2(wlx))
num_var = n + m * num_anc_hx + num_anc_logical
model = Model(GLPK.Optimizer) # model
set_silent(model)
@variable(model, x[1:num_var], Bin) # binary variables
@objective(model, Min, sum(x[i] for i in 1:n)) # objective function
# Orthogonality to rows of hx constraints
for row in 1:m
weight = zeros(Int, num_var)
supp = findall(hx[row, :] .!= 0) # Non-zero indices in hx[row, :]
for q in supp
weight[q] = 1
end
cnt = 1
for q in 1:num_anc_hx
weight[n + (row - 1) * num_anc_hx + q] = -(1 << cnt)
cnt += 1
end
@constraint(model, sum(weight[i] * x[i] for i in 1:num_var) == 0)
end
# Odd overlap with lx constraint
supp = findall(lx .!= 0) # Non-zero indices in lx
weight = zeros(Int, num_var)
for q in supp
weight[q] = 1
end
cnt = 1
for q in 1:num_anc_logical
weight[n + m * num_anc_hx + q] = -(1 << cnt)
cnt += 1
end
@constraint(model, sum(weight[i] * x[i] for i in 1:num_var) == 1)
optimize!(model)
opt_val = sum(value(x[i]) for i in 1:n)
return Int(opt_val)
end
14 changes: 14 additions & 0 deletions ext/QuantumCliffordJuMPExt/util.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
function get_stab_hx(matrix::SparseMatrixCSC{T, Int}) where T
rows, cols = size(matrix)
rhs_start = div(cols, 2) + 1
rhs_cols = matrix[:, rhs_start:cols]
non_zero_rows_rhs = unique(rhs_cols.rowval)
zero_rows_rhs = setdiff(1:rows, non_zero_rows_rhs)
return matrix[zero_rows_rhs, :]
end

function get_lx_lz(c::Stabilizer)
lx = stab_to_gf2(logicalxview(canonicalize!(MixedDestabilizer(c))))
lz = stab_to_gf2(logicalzview(canonicalize!(MixedDestabilizer(c))))
return lx, lz
end
12 changes: 11 additions & 1 deletion src/ecc/ECC.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqP
abstract type AbstractECC end

export parity_checks, parity_checks_x, parity_checks_z, iscss,
code_n, code_s, code_k, rate, distance,
code_n, code_s, code_k, rate, distance, minimum_distance,
isdegenerate, faults_matrix,
naive_syndrome_circuit, shor_syndrome_circuit, naive_encoding_circuit,
RepCode, LiftedCode,
Expand Down Expand Up @@ -123,6 +123,16 @@ end
"""The distance of a code."""
function distance end

"""The minimum distance of a qLDPC code.
Requires a JuMP.jl backend to be loaded, e.g. `using JuMP`.
Requires a GLPK.jl backend to be loaded, e.g. `using GLPK`.
"""
function minimum_distance end

minimum_distance(c::AbstractECC; kwargs...) = minimum_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))
Expand Down

0 comments on commit 71b5f9e

Please sign in to comment.