Recently I'm reading the paper of Bivariate Bicycle(BB) QLDPC code proposed by IBM. The paper mentioned that:

The actual distance of each candidate code was computed using the integer linear programming method

and it refered to this paper for the mixed ILP method for computing the exact code distance. This referenced paper primarily discusses ILP for MLE decoding in color codes. I am struggling to understand how the ILP method for MLE decoding can be adapted to compute code distances. Is there any connection between these two problems?

  • $\begingroup$ For classical codes this paper has details : scholar.google.com/… I tried it on small codes with matlab and seems to work ok but I don't think it scales very well. $\endgroup$
    – unknown
    Mar 16 at 14:58
  • 2
    $\begingroup$ A very concrete example is this pull request being written for stim, which will convert a circuit into a WDIMACS description of the shortest error as a maxsat problem. $\endgroup$ Mar 16 at 17:42

2 Answers 2


Here is python code that generates the stabilizers and observables of the IBM code, then produces a stim circuit measuring the Z stabilizers and observables after applying noise, then turns the distance-finding problem into a maxSAT problem solved by pysat:

import time

import stim

def make_code() -> tuple[list[stim.PauliString], list[stim.PauliString], list[stim.PauliString]]:
    offsets = {1, -1, 1j, -1j, 3 + 6j, -6 + 3j}
    w = 24
    h = 12

    def wrap(c: complex) -> complex:
        return c.real % w + (c.imag % h)*1j

    def index_of(c: complex) -> int:
        c = wrap(c)
        return int(c.real + c.imag * w) // 2

    stabilizers: list[stim.PauliString] = []
    for x in range(w):
        for y in range(h):
            if x % 2 != y % 2:
                continue  # This is a data qubit.
            m = x + 1j*y
            basis = 'XZ'[x % 2]
            sign = -1 if basis == 'Z' else +1
            stabilizer = stim.PauliString(w * h // 2)
            for offset in offsets:
                stabilizer[index_of(m + offset * sign)] = basis

    completed_tableau = stim.Tableau.from_stabilizers(
    obs_indices = [
        for k in range(len(completed_tableau))
        if completed_tableau.z_output(k) not in stabilizers
    observable_xs: list[stim.PauliString] = [
        for k in obs_indices
    observable_zs: list[stim.PauliString] = [
        for k in obs_indices

    return stabilizers, observable_xs, observable_zs

def make_z_basis_circuit(stabilizers: list[stim.PauliString], obs_zs: list[stim.PauliString]):
    num_qubits = len(stabilizers[0])
    circuit = stim.Circuit()
    circuit.append("X_ERROR", range(num_qubits), 1e-3)
    for k, observable in enumerate(obs_zs):
        circuit.append("MPP", stim.target_combined_paulis(observable))
        circuit.append("OBSERVABLE_INCLUDE", stim.target_rec(-1), k)
    for stabilizer in stabilizers:
        if stabilizer.pauli_indices('Z'):
            circuit.append("MPP", stim.target_combined_paulis(stabilizer))
            circuit.append("DETECTOR", stim.target_rec(-1))
    return circuit

def main():
    t0 = time.monotonic()

    stabilizers, obs_xs, obs_zs = make_code()
    circuit = make_z_basis_circuit(stabilizers, obs_zs)
    wcnf_string = circuit.shortest_error_sat_problem(format='WDIMACS')
    t1 = time.monotonic()
    print(f"Problem created in {t1 - t0:0.3f}s")

    from pysat.examples.rc2 import RC2
    from pysat.formula import WCNF
    wcnf = WCNF(from_string=wcnf_string)
    with RC2(wcnf) as rc2:
        print(f"distance = {rc2.cost}")
    t2 = time.monotonic()
    print(f"Problem solved in {t2 - t1:0.3f}s")

if __name__ == '__main__':

This code will run for a minute or two then print:

Problem created in 0.012s
distance = 12
Problem solved in 112.547s

The author of the paper provided the implementation in the repository.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.