1
$\begingroup$

I can generate encoder/decoder circuits for arbitrary QECC in qiskit. I tested a number of codes : $[[5,1,3]], [[7,1,3]],[[8,3,3]],\cdots$ and everything seems to work fine. The problem is that the simulation time is very slow (~1 sec per codeword). I know stim has been interfaced to pymatching decoder and can run sims much faster; is there a way to develop a similar interface to my decoders? would there be a significant increase in speed or would the bottleneck still be the qiskit decoder part? As a side, I can also generate the encoder and syndrome part in stim; the decoder uses non-clifford gates so it can't be done in stim as far as I know.

Here's an example $[[5,1,3]]$ encoder/decoder circuit and the way I checked that it can correct single pauli errors on any of the 5 qubits.

import numpy as np
import itertools
import math
import qiskit

def Test():
 backend=qiskit.Aer.get_backend('qasm_simulator')
 n=5; k=1; m=4; t=1;
 qr=qiskit.QuantumRegister(n+m,name="qr")
 cr=qiskit.ClassicalRegister(n+m,name="cr")

 qcenc=qiskit.QuantumCircuit(qr,cr);
 qcenc.h(0); qcenc.cx(0,4);
 qcenc.h(1); qcenc.cz(1,0); qcenc.cz(1,4); qcenc.cx(1,4);
 qcenc.h(2); qcenc.cz(2,1); qcenc.cz(2,4); qcenc.cx(2,4);
 qcenc.h(3); qcenc.cz(3,0); qcenc.cz(3,1); qcenc.cx(3,4);
 qcenc.barrier()

 qcsyn=qiskit.QuantumCircuit(qr,cr);
 qcsyn.h(5); qcsyn.h(6); qcsyn.h(7); qcsyn.h(8);
 qcsyn.cx(5,0); qcsyn.cz(5,2); qcsyn.cz(5,3); qcsyn.cx(5,4);
 qcsyn.cz(6,0); qcsyn.cz(6,1); qcsyn.cx(6,1); qcsyn.cz(6,2); qcsyn.cz(6,4); qcsyn.cx(6,4);
 qcsyn.cz(7,1); qcsyn.cz(7,2); qcsyn.cx(7,2); qcsyn.cz(7,3); qcsyn.cz(7,4); qcsyn.cx(7,4);
 qcsyn.cz(8,0); qcsyn.cz(8,1); qcsyn.cx(8,3); qcsyn.cx(8,4);
 qcsyn.h(5); qcsyn.h(6); qcsyn.h(7); qcsyn.h(8);
 qcsyn.barrier()

 qccor=qiskit.QuantumCircuit(qr,cr);
 qccor.x(5); qccor.x(7);
 gate=qiskit.QuantumCircuit(5);  gate.x(0);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(5); qccor.x(7);
 qccor.x(6); qccor.x(7); qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.z(0);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(6); qccor.x(7); qccor.x(8);
 qccor.x(7);
 gate=qiskit.QuantumCircuit(5);  gate.x(0); gate.z(0);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(7);
 qccor.x(5);
 gate=qiskit.QuantumCircuit(5);  gate.x(1);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(5);
 qccor.x(5); qccor.x(7); qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.z(1);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(5); qccor.x(7); qccor.x(8);
 qccor.x(5); qccor.x(6);
 gate=qiskit.QuantumCircuit(5);  gate.x(1); gate.z(1);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(5); qccor.x(6);
 qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.x(2);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(8);
 qccor.x(5); qccor.x(6); qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.z(2);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(5); qccor.x(6); qccor.x(8);
 qccor.x(7); qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.x(2); gate.z(2);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(7); qccor.x(8);
 qccor.x(6); qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.x(3);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(6); qccor.x(8);
 qccor.x(5); qccor.x(6); qccor.x(7);
 gate=qiskit.QuantumCircuit(5);  gate.z(3);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(5); qccor.x(6); qccor.x(7);
 qccor.x(6);
 gate=qiskit.QuantumCircuit(5);  gate.x(3); gate.z(3);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(6);
 qccor.x(5); qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.x(4);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(5); qccor.x(8);
 gate=qiskit.QuantumCircuit(5);  gate.z(4);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(6); qccor.x(7);
 gate=qiskit.QuantumCircuit(5);  gate.x(4); gate.z(4);
 qccor.append(instruction=gate.control(num_ctrl_qubits=4),qargs=[ 5, 6, 7, 8, 0, 1, 2, 3, 4 ])
 qccor.x(6); qccor.x(7);
 qccor.barrier()

 qcenx=qcenc.inverse();

 WrdNum=0;
 WrdErr=0;
 locs=itertools.combinations(range(n),t)
 for i in range(math.comb(n,t)):
  loc=next(locs)
  errs=itertools.product([0,1,2], repeat=t)
  for j in range(3**t):
   err=next(errs);
   TxBits=itertools.product([0,1], repeat=k)
   for h in range(2**k):
    txbits=next(TxBits)
    qcini=qiskit.QuantumCircuit(qr,cr);
    for ii in range(m): qcini.z(i)
    for ii in range(k): 
     if(txbits[ii]==1):
      qcini.x(m+ii)
     else:
      qcini.z(m+ii)
    for ii in range(m): qcini.z(n+ii)
    chn=qiskit.QuantumCircuit(qr);
    for kk in range(t):
     if(err[kk]==0):chn.x(loc[kk])
     if(err[kk]==1):chn.z(loc[kk])
     if(err[kk]==2):chn.y(loc[kk])
    chn.barrier()
    qc=qcini+qcenc+chn+qcsyn+qccor+qcenx
    for i in range(n+m): qc.measure(qr[i],cr[i])
    job=qiskit.execute(qc,backend,shots=1)
    result=job.result()
    counts=result.get_counts(qc)
    Counts = [(kx[::-1],v) for kx,v in counts.items()]
    rxbits=[];
    for kx, v in Counts: rxbits.append( [ int(c) for c in kx ] )
    WrdNum=WrdNum+1;
    werr=0;
    for ii in range(k):
     if(rxbits[0][m+ii]!=txbits[ii]):werr=1;
    WrdErr=WrdErr+werr;
    print("measured=",counts, end =" ")
    print("tx=", end =" ")
    for ii in range(k): print(txbits[ii], end =" ")
    print("rx=", end =" ")
    for ii in range(k): print(rxbits[0][m+ii], end =" ")
    print("err loc=",loc,end =" ")
    print("err val=",err,end =" ")
    print("decoding errors=",WrdErr,"/",WrdNum,end =" ")
    print("");

If the above is in file "code5.py" then the results would like this

>>> import code5
>>> code5.Test()
measured= {'101000000': 1} tx= 0 rx= 0 err loc= (0,) err val= (0,) decoding errors= 0 / 1
measured= {'101010000': 1} tx= 1 rx= 1 err loc= (0,) err val= (0,) decoding errors= 0 / 2
measured= {'000100000': 1} tx= 0 rx= 0 err loc= (0,) err val= (1,) decoding errors= 0 / 3
measured= {'000110000': 1} tx= 1 rx= 1 err loc= (0,) err val= (1,) decoding errors= 0 / 4
measured= {'101100000': 1} tx= 0 rx= 0 err loc= (0,) err val= (2,) decoding errors= 0 / 5
measured= {'101110000': 1} tx= 1 rx= 1 err loc= (0,) err val= (2,) decoding errors= 0 / 6
measured= {'111000000': 1} tx= 0 rx= 0 err loc= (1,) err val= (0,) decoding errors= 0 / 7
measured= {'111010000': 1} tx= 1 rx= 1 err loc= (1,) err val= (0,) decoding errors= 0 / 8
measured= {'001000000': 1} tx= 0 rx= 0 err loc= (1,) err val= (1,) decoding errors= 0 / 9
measured= {'001010000': 1} tx= 1 rx= 1 err loc= (1,) err val= (1,) decoding errors= 0 / 10
measured= {'110000000': 1} tx= 0 rx= 0 err loc= (1,) err val= (2,) decoding errors= 0 / 11
measured= {'110010000': 1} tx= 1 rx= 1 err loc= (1,) err val= (2,) decoding errors= 0 / 12
measured= {'011100000': 1} tx= 0 rx= 0 err loc= (2,) err val= (0,) decoding errors= 0 / 13
measured= {'011110000': 1} tx= 1 rx= 1 err loc= (2,) err val= (0,) decoding errors= 0 / 14
measured= {'010000000': 1} tx= 0 rx= 0 err loc= (2,) err val= (1,) decoding errors= 0 / 15
measured= {'010010000': 1} tx= 1 rx= 1 err loc= (2,) err val= (1,) decoding errors= 0 / 16
measured= {'001100000': 1} tx= 0 rx= 0 err loc= (2,) err val= (2,) decoding errors= 0 / 17
measured= {'001110000': 1} tx= 1 rx= 1 err loc= (2,) err val= (2,) decoding errors= 0 / 18
measured= {'010100000': 1} tx= 0 rx= 0 err loc= (3,) err val= (0,) decoding errors= 0 / 19
measured= {'010110000': 1} tx= 1 rx= 1 err loc= (3,) err val= (0,) decoding errors= 0 / 20
measured= {'100000000': 1} tx= 0 rx= 0 err loc= (3,) err val= (1,) decoding errors= 0 / 21
measured= {'100010000': 1} tx= 1 rx= 1 err loc= (3,) err val= (1,) decoding errors= 0 / 22
measured= {'110100000': 1} tx= 0 rx= 0 err loc= (3,) err val= (2,) decoding errors= 0 / 23
measured= {'110110000': 1} tx= 1 rx= 1 err loc= (3,) err val= (2,) decoding errors= 0 / 24
measured= {'011000000': 1} tx= 0 rx= 0 err loc= (4,) err val= (0,) decoding errors= 0 / 25
measured= {'011010000': 1} tx= 1 rx= 1 err loc= (4,) err val= (0,) decoding errors= 0 / 26
measured= {'111100000': 1} tx= 0 rx= 0 err loc= (4,) err val= (1,) decoding errors= 0 / 27
measured= {'111110000': 1} tx= 1 rx= 1 err loc= (4,) err val= (1,) decoding errors= 0 / 28
measured= {'100100000': 1} tx= 0 rx= 0 err loc= (4,) err val= (2,) decoding errors= 0 / 29
measured= {'100110000': 1} tx= 1 rx= 1 err loc= (4,) err val= (2,) decoding errors= 0 / 30
>>>
Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

Here's the stim equivalent. A few notes:

  • It finishes checking all 64 cases in 70 milliseconds. The reason it's so slow is because I'm driving the tableau simulator step by step using python, instead of properly turning the problem into a stim.Circuit that can be bulk sampled.
  • I had to correct your bad usage of $CX \cdot CZ$ instead of CY in order to make the circuit work. $CX \cdot CZ$ creates a phase kickback of $i$ onto the control which randomizes stabilizer measurements. You were only getting away with this because your mistakes happened to be cancelling each other out. Always always use $Y$ instead of $X \cdot Z$.
  • I corrected the fact that the code was only checking the Z logical observable, not both X and Z logical.
  • I replaced your block of hard coded quantum operations implementing the corrections with the conversion of the failing stabilizer bits into an integer and a lookup into a list of corrections. The "decoder" is basically just the contents of that list.
  • There are open issues to add the pauli string iteration and pauli string measurement as built-in methods.
from typing import Iterator

import numpy as np
import stim


def iter_pauli_strings(n: int) -> Iterator[stim.PauliString]:
    assert n <= 8
    for x in range(1 << n):
        for z in range(1 << n):
            yield stim.PauliString.from_numpy(
                xs=np.array([x], dtype=np.uint8),
                zs=np.array([z], dtype=np.uint8),
                num_qubits=n,
            )


def measure_pauli_product(simulator: stim.TableauSimulator, paulis: stim.PauliString) -> bool:
    targets = []
    for k, p in enumerate(paulis):
        if p == 1:
            t = stim.target_x
        elif p == 2:
            t = stim.target_y
        elif p == 3:
            t = stim.target_z
        else:
            continue
        targets.append(t(k))
        targets.append(stim.target_combiner())
    targets.pop()
    simulator.do(stim.CircuitInstruction(name='MPP', targets=targets))
    return simulator.current_measurement_record()[-1]


encoder = stim.Circuit("""
    H 0
    CX 0 4
    H 1
    CZ 1 0
    CY 1 4
    H 2
    CZ 2 1
    CY 2 4
    H 3
    CZ 3 0
    CZ 3 1
    CX 3 4
""")
stabilizers = [
    stim.PauliString("ZXZ_Y"),
    stim.PauliString("_ZXZY"),
    stim.PauliString("ZZ_XX"),
    stim.PauliString("X_ZZX"),
]
corrections = [
    stim.PauliString("+_____"),
    stim.PauliString("+Z____"),
    stim.PauliString("+___Z_"),
    stim.PauliString("+____Y"),
    stim.PauliString("+__Z__"),
    stim.PauliString("+___X_"),
    stim.PauliString("+_X___"),
    stim.PauliString("+___Y_"),
    stim.PauliString("+_Z___"),
    stim.PauliString("+__X__"),
    stim.PauliString("+X____"),
    stim.PauliString("+Y____"),
    stim.PauliString("+____X"),
    stim.PauliString("+__Y__"),
    stim.PauliString("+_Y___"),
    stim.PauliString("+____Z"),
]

max_error_weight = 1
num_qubits = 5
checked_qubits = [0, 1, 2, 3, 4]

any_failures = False
for state in '01+-':
    for err in iter_pauli_strings(num_qubits):
        if sum(p != 0 for p in err) > max_error_weight:
            continue
        sim = stim.TableauSimulator()

        if state == '0':
            init = stim.Circuit()
        elif state == '1':
            init = stim.Circuit("""
                X 4
            """)
        elif state == '+':
            init = stim.Circuit("""
                H 4
            """)
        elif state == '-':
            init = stim.Circuit("""
                H 4
                Z 4
            """)
        else:
            raise NotImplementedError(f'{state=}')

        sim.do(init)
        sim.do(encoder)

        # Apply chosen noise.
        sim.do(err)

        # Pick correction based on flipped stabilizers.
        fault_index = 0
        for stabilizer in stabilizers:
            fault_index *= 2
            fault_index += measure_pauli_product(sim, stabilizer)
        sim.do(corrections[fault_index])

        # caution: assumes each operation in the circuit is self-inverse
        sim.do(encoder[::-1])
        sim.do(init[::-1])

        if any(sim.measure_many(*checked_qubits)):
            print("Failed to correct", err, 'on logical state', state)
            any_failures = True
if not any_failures:
    print("All corrected")
Improve this answer
$\endgroup$
3
  • $\begingroup$ That's pretty fast! The CXCZ vs CY is a recurring issue which I'll resolve later. I should be able to automate the generation of the correction table for different codes. One thing to note is that my circuit performs the correction without measurement; the measurement is done to verify things; your correction is done based on the classical bits you get from measurements. For checking the code performance either approach will work. $\endgroup$
    – unknown
    Oct 28, 2022 at 20:15
  • 1
    $\begingroup$ @unknown I generated the corrections list by printing out the fault index and the error, sorted by the fault index. $\endgroup$
    – Craig Gidney
    Oct 28, 2022 at 20:18
  • $\begingroup$ All the circuits (including the hard coded decoder) are automatically generated from a high level description of the code. Generating the corresponding correction table should be straight forward...I'm more concerned about a clean way to handle the XZ vs Y convention... $\endgroup$
    – unknown
    Oct 28, 2022 at 20:25

Your Answer

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

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