From 21ceebd35e862f63ee56eb91ce915ba25bbb9e3e Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sun, 20 Oct 2024 22:38:44 +0800 Subject: [PATCH 01/11] a runable OSD --- src/LDPCDecoders.jl | 3 + src/decoders/belief_propagation_osd.jl | 115 +++++++++++++++++++++++++ 2 files changed, 118 insertions(+) create mode 100644 src/decoders/belief_propagation_osd.jl diff --git a/src/LDPCDecoders.jl b/src/LDPCDecoders.jl index 86b70df..a34265e 100644 --- a/src/LDPCDecoders.jl +++ b/src/LDPCDecoders.jl @@ -12,6 +12,7 @@ using RowEchelon export decode!, batchdecode!, BeliefPropagationDecoder, + BeliefPropagationOSDDecoder, BitFlipDecoder include("generator.jl") @@ -22,7 +23,9 @@ include("parity_generator.jl") include("decoders/abstract_decoder.jl") include("decoders/belief_propagation.jl") +include("decoders/belief_propagation_osd.jl") include("decoders/iterative_bitflip.jl") + include("syndrome_bp_decoder.jl") include("syndrome_simulator.jl") include("syndrome_it_decoder.jl") diff --git a/src/decoders/belief_propagation_osd.jl b/src/decoders/belief_propagation_osd.jl new file mode 100644 index 0000000..8e86445 --- /dev/null +++ b/src/decoders/belief_propagation_osd.jl @@ -0,0 +1,115 @@ +struct BeliefPropagationOSDDecoder <: AbstractDecoder + bp_decoder::BeliefPropagationDecoder + H::BitMatrix # use dense matrix + osd_order::Int +end + +function BeliefPropagationOSDDecoder(H, per::Float64, max_iters::Int; osd_order::Int=0) + bp_decoder = BeliefPropagationDecoder(H, per, max_iters) + return BeliefPropagationOSDDecoder(bp_decoder, H, osd_order) +end + +function rowswap!(H::BitMatrix, i, j) + @inbounds H[i, :], H[j, :] = H[j, :], H[i, :] # TODO This could be further optimized? +end + +function decode!(decoder::BeliefPropagationOSDDecoder, syndrome::AbstractVector) + # use BP to get hard and soft decisions + bp_err, converged = decode!(decoder.bp_decoder, syndrome) # hard decisions + bp_log_probabs = decoder.bp_decoder.scratch.log_probabs # soft decisions + bp_probabs = exp.(bp_log_probabs) + # sort colums by reliability, less reliable columns first + sort_by_reliability = sortperm(max.(bp_probabs, 1 .- bp_probabs), rev=true) + H_sorted = decoder.H[:, sort_by_reliability] + bp_err_sorted = bp_err[sort_by_reliability] + if decoder.osd_order == 0 + err = fastosd0(H_sorted, syndrome, bp_err_sorted) + else + err = osd(H_sorted, syndrome, bp_err_sorted, decoder.osd_order) + end + return err[invperm(sort_by_reliability)], converged # also return whether BP is converged +end + +function osd(H, syndrome, bp_err, osd_order) + m, n = size(H) + # diagnolize the submatrix corrsponding to idependent columns via Gaussian elimination + # first obtain the row echelon form + # and find least reliable indices, i.e., the first r pivot columns (assume H is rearranged by reliability) + least_reliable_rows = [] # row indices of pivot elements + least_reliable_cols = [] # column indices of pivot elements + r = 0 # compute rank of H + i, j = 1, 1 + s = copy(syndrome) # tranform syndrom along with H in Gaussian elimination + + while i <= m && j <= n + k = findfirst(H[i:end, j]) + if isnothing(k) # not an independent column + j += 1 + else + if k > 1 + ii = i + k - 1 + rowswap!(H, i, ii) + s[i], s[ii] = s[ii], s[i] + end + for ii in i+1:m + if H[ii, j] + H[ii, :] .⊻= H[i, :] + s[ii] ⊻= s[i] + end + end + push!(least_reliable_rows, i) + push!(least_reliable_cols, j) + i += 1 + j += 1 + r += 1 + end + end + + # then obtain a diagnol submatrix on the least reliable part + for (i, j) in zip(reverse(least_reliable_rows), reverse(least_reliable_cols)) + for ii in 1:i-1 + if H[ii, j] + H[ii, :] .⊻= H[i, :] + s[ii] ⊻= s[i] + end + end + end + + if osd_order > n - r + @warn "The order of OSD $osd_order is greater than the size of the information set $(n-r). We set osd_order = $(n-r)." + osd_order = n - r + end + + best_err = copy(bp_err) + err = Bool.(copy(bp_err)) # TODO why error is in Float in BP? + most_reliable_cols = setdiff(1:n, least_reliable_cols) + min_weight = n + 1 + + for x in 0:2^osd_order-1 + # first compute the `most_reliable_cols` part of errors + # try all possible errors on the first `osd_order` bits within `most_reliable_cols` + if x != 0 + trial_err = BitArray([x >> i & 1 for i in 0:osd_order-1]) + err[most_reliable_cols[1:osd_order]] = trial_err + end + # then based on the `most_reliable_cols` part of errors, compute the `least_reliable_cols` part of errors + for (i, j) in zip(least_reliable_rows, least_reliable_cols) + err[j] = s[j] + for k in most_reliable_cols + err[j] ⊻= H[i, k] * err[k] + end + end + weight = sum(err) # TODO Actually it should be a function depending on the noise model + if weight < min_weight + min_weight = weight + best_err = copy(err) + end + end + + return best_err +end + +function fastosd0(H, syndrome, bp_err) + # TODO a optmized special case of osd with osd_order = 0 to be implemented + osd(H, syndrome, bp_err, 0) +end From 8015972fff30e8fccba1a690296a9c04b743f04c Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Mon, 4 Nov 2024 12:35:46 +0800 Subject: [PATCH 02/11] comments --- src/decoders/belief_propagation_osd.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/decoders/belief_propagation_osd.jl b/src/decoders/belief_propagation_osd.jl index 8e86445..c7aa7e0 100644 --- a/src/decoders/belief_propagation_osd.jl +++ b/src/decoders/belief_propagation_osd.jl @@ -33,7 +33,7 @@ end function osd(H, syndrome, bp_err, osd_order) m, n = size(H) # diagnolize the submatrix corrsponding to idependent columns via Gaussian elimination - # first obtain the row echelon form + # first obtain the row canonical form # and find least reliable indices, i.e., the first r pivot columns (assume H is rearranged by reliability) least_reliable_rows = [] # row indices of pivot elements least_reliable_cols = [] # column indices of pivot elements @@ -47,8 +47,8 @@ function osd(H, syndrome, bp_err, osd_order) j += 1 else if k > 1 - ii = i + k - 1 - rowswap!(H, i, ii) + ii = i + k - 1 # the first row after `i` with 1 in column `j` + rowswap!(H, i, ii) # TODO Is this swap necessary? We may just track the row index s[i], s[ii] = s[ii], s[i] end for ii in i+1:m From 553d0e5f9798eb4e49944c3867ff1c2ebbde0917 Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sat, 9 Nov 2024 19:03:36 +0800 Subject: [PATCH 03/11] fix an index bug --- src/decoders/belief_propagation_osd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/decoders/belief_propagation_osd.jl b/src/decoders/belief_propagation_osd.jl index c7aa7e0..1cc54eb 100644 --- a/src/decoders/belief_propagation_osd.jl +++ b/src/decoders/belief_propagation_osd.jl @@ -94,7 +94,7 @@ function osd(H, syndrome, bp_err, osd_order) end # then based on the `most_reliable_cols` part of errors, compute the `least_reliable_cols` part of errors for (i, j) in zip(least_reliable_rows, least_reliable_cols) - err[j] = s[j] + err[j] = s[i] for k in most_reliable_cols err[j] ⊻= H[i, k] * err[k] end From 8fc81f5a9c8953ecc845930d84ce4276fc42a315 Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sat, 9 Nov 2024 14:52:14 +0100 Subject: [PATCH 04/11] add basic tests --- test/test_bposd_decoder.jl | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 test/test_bposd_decoder.jl diff --git a/test/test_bposd_decoder.jl b/test/test_bposd_decoder.jl new file mode 100644 index 0000000..c2b94f6 --- /dev/null +++ b/test/test_bposd_decoder.jl @@ -0,0 +1,34 @@ +using Test +using LDPCDecoders + +@testset "test_bposd_decoder.jl" begin + + """Test for BP-OSD decoder""" + function test_bposd_decoder() + H = LDPCDecoders.parity_check_matrix(1000, 10, 9) + per = 0.01 + err = rand(1000) .< per + syn = (H * err) .% 2 + + bposd = BeliefPropagationOSDDecoder(H, per, 100) + guess, success = decode!(bposd, syn) + + return guess == err + end + + """Test""" + function test_bposd_decoder_large_error_rate() + H = LDPCDecoders.parity_check_matrix(1000, 10, 9) + per = 0.2 + err = rand(1000) .< per + syn = (H * err) .% 2 + + bposd = BeliefPropagationOSDDecoder(H, per, 100) + guess, success = decode!(bposd, syn) + + return syn == (H * guess) .% 2 + end + + @test test_bp_decoder() + @test test_bposd_decoder_large_error_rate() +end From 543e98da6c948bae4385c70786dad5430e3433fb Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sat, 9 Nov 2024 22:03:55 +0800 Subject: [PATCH 05/11] update tests --- test/runtests.jl | 1 + test/test_bposd_decoder.jl | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 974db7c..969ce29 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,6 +30,7 @@ println("ENV[\"PYTHON\"] = \"$(get(ENV,"PYTHON",nothing))\"") @doset "oldtests" @doset "bp_decoder" +@doset "bposd_decoder" @doset "bf_decoder" VERSION >= v"1.10" && @doset "doctests" diff --git a/test/test_bposd_decoder.jl b/test/test_bposd_decoder.jl index c2b94f6..eb8ec42 100644 --- a/test/test_bposd_decoder.jl +++ b/test/test_bposd_decoder.jl @@ -16,7 +16,7 @@ using LDPCDecoders return guess == err end - """Test""" + """Test for BP-OSD decoder with large error rate. Even if the decoding is not accurate, OSD will still ensure consistency between guess and syndromes.""" function test_bposd_decoder_large_error_rate() H = LDPCDecoders.parity_check_matrix(1000, 10, 9) per = 0.2 @@ -29,6 +29,6 @@ using LDPCDecoders return syn == (H * guess) .% 2 end - @test test_bp_decoder() + @test test_bposd_decoder() @test test_bposd_decoder_large_error_rate() end From 911d18f80da4781e751e3ff3da47eac16bfc850b Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sat, 9 Nov 2024 22:13:09 +0800 Subject: [PATCH 06/11] update TODOs --- src/decoders/belief_propagation_osd.jl | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/src/decoders/belief_propagation_osd.jl b/src/decoders/belief_propagation_osd.jl index 1cc54eb..3668794 100644 --- a/src/decoders/belief_propagation_osd.jl +++ b/src/decoders/belief_propagation_osd.jl @@ -22,11 +22,8 @@ function decode!(decoder::BeliefPropagationOSDDecoder, syndrome::AbstractVector) sort_by_reliability = sortperm(max.(bp_probabs, 1 .- bp_probabs), rev=true) H_sorted = decoder.H[:, sort_by_reliability] bp_err_sorted = bp_err[sort_by_reliability] - if decoder.osd_order == 0 - err = fastosd0(H_sorted, syndrome, bp_err_sorted) - else - err = osd(H_sorted, syndrome, bp_err_sorted, decoder.osd_order) - end + # TODO an optmized version of OSD can be implemented when osd_order = 0, see Algorithm 2 in https://doi.org/10.22331/q-2021-11-22-585 + err = osd(H_sorted, syndrome, bp_err_sorted, decoder.osd_order) return err[invperm(sort_by_reliability)], converged # also return whether BP is converged end @@ -48,7 +45,7 @@ function osd(H, syndrome, bp_err, osd_order) else if k > 1 ii = i + k - 1 # the first row after `i` with 1 in column `j` - rowswap!(H, i, ii) # TODO Is this swap necessary? We may just track the row index + rowswap!(H, i, ii) # TODO For optimization: Is this swap necessary? We may just track the row index s[i], s[ii] = s[ii], s[i] end for ii in i+1:m @@ -99,7 +96,8 @@ function osd(H, syndrome, bp_err, osd_order) err[j] ⊻= H[i, k] * err[k] end end - weight = sum(err) # TODO Actually it should be a function depending on the noise model + weight = sum(err) # This weight is set for depolarizing noise + # TODO More generally, it should be a function depending on the noise model if weight < min_weight min_weight = weight best_err = copy(err) @@ -108,8 +106,3 @@ function osd(H, syndrome, bp_err, osd_order) return best_err end - -function fastosd0(H, syndrome, bp_err) - # TODO a optmized special case of osd with osd_order = 0 to be implemented - osd(H, syndrome, bp_err, 0) -end From cf861900f225924c7e77fa0b604d8b948be8e141 Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sat, 9 Nov 2024 22:16:45 +0800 Subject: [PATCH 07/11] update docs --- src/decoders/belief_propagation_osd.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/decoders/belief_propagation_osd.jl b/src/decoders/belief_propagation_osd.jl index 3668794..cc9f279 100644 --- a/src/decoders/belief_propagation_osd.jl +++ b/src/decoders/belief_propagation_osd.jl @@ -1,6 +1,9 @@ struct BeliefPropagationOSDDecoder <: AbstractDecoder + """A belief propagation decoder as a subroutine""" bp_decoder::BeliefPropagationDecoder - H::BitMatrix # use dense matrix + """Dense form of the parity check matrix""" + H::BitMatrix + """The order of OSD; defaulted to be 0 in the constructor""" osd_order::Int end From cf29b76afa95408f1ccfda7b4a140db851481d1d Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sat, 9 Nov 2024 22:34:52 +0800 Subject: [PATCH 08/11] add high osd_order tests --- test/test_bposd_decoder.jl | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/test/test_bposd_decoder.jl b/test/test_bposd_decoder.jl index eb8ec42..9bdf998 100644 --- a/test/test_bposd_decoder.jl +++ b/test/test_bposd_decoder.jl @@ -16,6 +16,24 @@ using LDPCDecoders return guess == err end + """Test high order OSD""" + function test_bposd_decoder_high_order() + H = LDPCDecoders.parity_check_matrix(1000, 10, 9) + per = 0.01 + err = rand(1000) .< per + syn = (H * err) .% 2 + + orders = 2:5 + succ = true + for osd_order in orders + bposd = BeliefPropagationOSDDecoder(H, per, 100; osd_order=osd_order) + guess, success = decode!(bposd, syn) + succ = succ & (guess == err) + end + + return succ + end + """Test for BP-OSD decoder with large error rate. Even if the decoding is not accurate, OSD will still ensure consistency between guess and syndromes.""" function test_bposd_decoder_large_error_rate() H = LDPCDecoders.parity_check_matrix(1000, 10, 9) @@ -30,5 +48,6 @@ using LDPCDecoders end @test test_bposd_decoder() + @test test_bposd_decoder_high_order() @test test_bposd_decoder_large_error_rate() end From 0afbc3dd1f85b8d14cf74a2a938b11592f038e91 Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Sat, 9 Nov 2024 22:37:37 +0800 Subject: [PATCH 09/11] fix typos --- src/decoders/belief_propagation_osd.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/decoders/belief_propagation_osd.jl b/src/decoders/belief_propagation_osd.jl index cc9f279..74052f5 100644 --- a/src/decoders/belief_propagation_osd.jl +++ b/src/decoders/belief_propagation_osd.jl @@ -21,25 +21,25 @@ function decode!(decoder::BeliefPropagationOSDDecoder, syndrome::AbstractVector) bp_err, converged = decode!(decoder.bp_decoder, syndrome) # hard decisions bp_log_probabs = decoder.bp_decoder.scratch.log_probabs # soft decisions bp_probabs = exp.(bp_log_probabs) - # sort colums by reliability, less reliable columns first + # sort columns by reliability, less reliable columns first sort_by_reliability = sortperm(max.(bp_probabs, 1 .- bp_probabs), rev=true) H_sorted = decoder.H[:, sort_by_reliability] bp_err_sorted = bp_err[sort_by_reliability] - # TODO an optmized version of OSD can be implemented when osd_order = 0, see Algorithm 2 in https://doi.org/10.22331/q-2021-11-22-585 + # TODO an optimized version of OSD can be implemented when osd_order = 0, see Algorithm 2 in https://doi.org/10.22331/q-2021-11-22-585 err = osd(H_sorted, syndrome, bp_err_sorted, decoder.osd_order) return err[invperm(sort_by_reliability)], converged # also return whether BP is converged end function osd(H, syndrome, bp_err, osd_order) m, n = size(H) - # diagnolize the submatrix corrsponding to idependent columns via Gaussian elimination + # diagnolize the submatrix corresponding to independent columns via Gaussian elimination # first obtain the row canonical form # and find least reliable indices, i.e., the first r pivot columns (assume H is rearranged by reliability) least_reliable_rows = [] # row indices of pivot elements least_reliable_cols = [] # column indices of pivot elements r = 0 # compute rank of H i, j = 1, 1 - s = copy(syndrome) # tranform syndrom along with H in Gaussian elimination + s = copy(syndrome) # transform syndrome along with H in Gaussian elimination while i <= m && j <= n k = findfirst(H[i:end, j]) @@ -65,7 +65,7 @@ function osd(H, syndrome, bp_err, osd_order) end end - # then obtain a diagnol submatrix on the least reliable part + # then obtain a diagonal submatrix on the least reliable part for (i, j) in zip(reverse(least_reliable_rows), reverse(least_reliable_cols)) for ii in 1:i-1 if H[ii, j] From 2fad88e194e256cc2b010abf4469fad2fe9c1723 Mon Sep 17 00:00:00 2001 From: Stefan Krastanov Date: Fri, 15 Nov 2024 15:13:57 -0500 Subject: [PATCH 10/11] changelog and version bump --- CHANGELOG.md | 7 +++++++ Project.toml | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) create mode 100644 CHANGELOG.md diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..099d941 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,7 @@ +# News + +## v0.3.2 - 2024-11-15 + +- Add a (still unoptimized) implementation of a BP OSD decoder. + +## Older - before 2021-10-28 unrecorded \ No newline at end of file diff --git a/Project.toml b/Project.toml index 80726a8..48437ec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LDPCDecoders" uuid = "3c486d74-64b9-4c60-8b1a-13a564e77efb" authors = ["Krishna Praneet Gudipaty", "Stefan Krastanov", "QuantumSavory contributors"] -version = "0.3.1" +version = "0.3.2" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" From 1570a2f9a5b11cab5b58204aa97bb25ba778ebcf Mon Sep 17 00:00:00 2001 From: Yuxuan Yan Date: Wed, 27 Nov 2024 09:19:21 +0800 Subject: [PATCH 11/11] fix JET warnings --- src/decoders/belief_propagation.jl | 8 ++++---- src/decoders/belief_propagation_osd.jl | 6 ++++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/decoders/belief_propagation.jl b/src/decoders/belief_propagation.jl index ef4e190..76c6351 100644 --- a/src/decoders/belief_propagation.jl +++ b/src/decoders/belief_propagation.jl @@ -44,7 +44,7 @@ struct BeliefPropagationDecoder <: AbstractDecoder scratch::BeliefPropagationScratchSpace end -function BeliefPropagationDecoder(H, per::Float64, max_iters::Int) +function BeliefPropagationDecoder(H::Union{SparseArrays.SparseMatrixCSC{Bool,Int}, BitMatrix}, per::Float64, max_iters::Int) s, n = size(H) sparse_H = sparse(H) sparse_HT = sparse(H') @@ -108,7 +108,7 @@ true function decode!(decoder::BeliefPropagationDecoder, syndrome::AbstractVector) # TODO check if casting to bitarrays helps with performance -- if it does, set up warnings to the user for cases where they have not done the casting reset!(decoder) rows::Vector{Int} = rowvals(decoder.sparse_H); - rowsT::Vector{Int} = rowvals(decoder.sparse_HT); + rowsT::Vector{Int} = rowvals(decoder.sparse_HT); setup = decoder.scratch for j in 1:decoder.n @@ -138,7 +138,7 @@ function decode!(decoder::BeliefPropagationDecoder, syndrome::AbstractVector) # for j in 1:decoder.n temp::Float64 = setup.channel_probs[j] / (1 - setup.channel_probs[j]) - + for k in nzrange(decoder.sparse_H, j) setup.bit_2_check[rows[k],j] = temp temp *= setup.check_2_bit[rows[k],j] @@ -166,7 +166,7 @@ function decode!(decoder::BeliefPropagationDecoder, syndrome::AbstractVector) # syndrome_decoded = (decoder.sparse_H * setup.err) .% 2 if all(syndrome_decoded .== syndrome) - converged = true + converged = true break # Break if converged end end diff --git a/src/decoders/belief_propagation_osd.jl b/src/decoders/belief_propagation_osd.jl index 74052f5..be6259a 100644 --- a/src/decoders/belief_propagation_osd.jl +++ b/src/decoders/belief_propagation_osd.jl @@ -7,7 +7,7 @@ struct BeliefPropagationOSDDecoder <: AbstractDecoder osd_order::Int end -function BeliefPropagationOSDDecoder(H, per::Float64, max_iters::Int; osd_order::Int=0) +function BeliefPropagationOSDDecoder(H::BitMatrix, per::Float64, max_iters::Int; osd_order::Int=0) bp_decoder = BeliefPropagationDecoder(H, per, max_iters) return BeliefPropagationOSDDecoder(bp_decoder, H, osd_order) end @@ -90,7 +90,9 @@ function osd(H, syndrome, bp_err, osd_order) # try all possible errors on the first `osd_order` bits within `most_reliable_cols` if x != 0 trial_err = BitArray([x >> i & 1 for i in 0:osd_order-1]) - err[most_reliable_cols[1:osd_order]] = trial_err + for j in 1:osd_order + err[most_reliable_cols[j]] = trial_err[j] + end end # then based on the `most_reliable_cols` part of errors, compute the `least_reliable_cols` part of errors for (i, j) in zip(least_reliable_rows, least_reliable_cols)