Problems measuring logical T+ state of a surface QEC code

[sergei-mironov]

Let me report a possible problem in the near-Clifford simulation. The below code initializes a surface quantum error correction code in a T+ state and calculates 200 shots of its logical measurements out of raw 'physical' ones in X basis. The logical measurements are then saved into the logical_measurements.json file. If one puts the default 'lightning.qubit' device of PennyLane, the statistics would be something like 170 0s against 30 1s, which approaches the predicted value of sqrt(2)/2 ~= 0.7. However, with the shown QRack device I get inconsistent results which is far from expectations (e.g. 91/109 or 114/86).

I'm using pennylane 0.43.1, pennylane_catalyst 0.13.0, pennylane_qrack 0.25.0 (compiled from source). The simulation is done with isTensorNetwork flag set to False.

Can we improve the precision of this simulation mode? Does it have a maximum-precision mode which would match the statevector method's precision? Please let me know, I will be glad to debug or adjust parameters.

--

The runnable program:

import json
import time
from contextlib import contextmanager
from signal import ITIMER_REAL, setitimer
from functools import reduce

import pennylane as qml
from catalyst import qjit, measure as catalyst_measure

num_shots = 200
nattempts = 1
timeout = 3000.0
dry_run = False


@contextmanager
def with_alarm(timeout: float):
    timer_set = False
    try:
        if timeout is not None and 0 < timeout < float("inf"):
            setitimer(ITIMER_REAL, timeout)
            timer_set = True
        yield
    finally:
        if timer_set:
            setitimer(ITIMER_REAL, 0)


def dump_logical_observables(measurement_attempts) -> None:
    """ For each observable, calculate its value from measurements (by xor-ing them) and dump to the
    configured JSON file.  """
    logical_values_attempts = []
    logical_measurements = {'q1_X_0': ('c_(0, 0, 0)_0', 'c_(1, 0, 0)_0', 'c_(2, 0, 0)_0', 'c_(2, 1, 1)_0', 'c_(0, 1, 1)_0')}

    for attempt_shots in measurement_attempts:
        if not attempt_shots:
            continue
        nshots = len(attempt_shots)
        logical_values = {}

        for lol, meas_labels in logical_measurements.items():
            lo_shots = []
            for i in range(nshots):
                shot_dict = attempt_shots[i]
                samples = [int(shot_dict[l]) for l in meas_labels]
                lo_value = reduce(lambda a,b: a ^ b, samples)
                lo_shots.append(int(lo_value))
            logical_values[lol] = lo_shots

        logical_values_attempts.append(logical_values)

    with open("logical_measurements.json", "w") as f:
        json.dump(logical_values_attempts, f)

if not dry_run:
    @qjit
    @qml.qnode(qml.device("qrack.simulator",wires=17,isTensorNetwork=False,isStabilizerHybrid=True,isOpenCL=True,isPaged=False,isCpuGpuHybrid=False,), shots=num_shots)
    def pl_main():

        # PennyLane Program Generated from Eka
        # Initialize Wires from Eka quantum channels Ids
        measurements = {}
        qml.registers({'a643ceef-00bc-44cd-8514-b7b8cbe8bf95': 1, 'f0252ea6-ffa6-40f9-a8cf-aa886542e808': 1, '64a8a798-6799-4dde-8eb3-c4e1c3591cc8': 1, '18823f07-02b7-4965-9f2b-63737573cd92': 1, 'e299cd6f-5b93-4001-92bc-eff2ac056fff': 1, 'e6f6e48b-741f-471a-a0f6-20c0006971cd': 1, 'bf7ed026-3336-4705-a285-17c9b788d74d': 1, '3adeb324-cb13-4eaf-9c84-6d286ee08c4b': 1, '18788481-309e-4b4f-b8de-f81f5ebe0b49': 1, '123d3886-b07a-4655-a739-5f9ef84a1170': 1, '1fb8be97-58fa-48c3-812d-e946b2e6e58d': 1, 'c2e764ee-9cb3-4115-82e7-384bda96e261': 1, 'e6dd20b2-4b29-4abd-9bf5-25f109954f51': 1, '25b6c46e-948f-4ffd-a5a3-ca9c190d6598': 1, '7fb8aa79-b6a4-4a31-81fb-705992eee653': 1, 'd5c2ab4f-ac70-4df2-8413-60d7219cdc2c': 1, 'ca4f79db-08f6-4011-9d58-11f5a350911c': 1})
        
        # Start of circuit: final circuit
        # Start of circuit: inject t into block q1 and measure syndromes
        # Start of circuit: inject t into block q1
        # Start of circuit: resource state reset
        catalyst_measure(5, reset=True)
        qml.Hadamard([5])
        qml.T([5])
        #  End of circuit: resource state reset
        # Start of circuit: reset four quadrants
        catalyst_measure(0, reset=True)
        qml.Hadamard([0])
        catalyst_measure(1, reset=True)
        qml.Hadamard([1])
        catalyst_measure(12, reset=True)
        qml.Hadamard([12])
        catalyst_measure(14, reset=True)
        qml.Hadamard([14])
        catalyst_measure(4, reset=True)
        catalyst_measure(10, reset=True)
        catalyst_measure(3, reset=True)
        catalyst_measure(7, reset=True)
        #  End of circuit: reset four quadrants
        #  End of circuit: inject t into block q1
        # Start of circuit: measure q1 syndromes 1 time(s)
        # Start of circuit: measure q1 syndromes - cycle 0
        catalyst_measure(6, reset=True)
        catalyst_measure(15, reset=True)
        catalyst_measure(8, reset=True)
        catalyst_measure(13, reset=True)
        catalyst_measure(2, reset=True)
        catalyst_measure(16, reset=True)
        catalyst_measure(11, reset=True)
        catalyst_measure(9, reset=True)
        qml.Hadamard([6])
        qml.Hadamard([15])
        qml.Hadamard([8])
        qml.Hadamard([13])
        qml.Hadamard([2])
        qml.Hadamard([16])
        qml.Hadamard([11])
        qml.Hadamard([9])
        qml.CZ([6, 4])
        qml.CZ([15, 12])
        qml.CNOT([8, 5])
        qml.CNOT([13, 10])
        qml.CNOT([2, 0])
        qml.CZ([9, 7])
        qml.CZ([6, 0])
        qml.CZ([15, 5])
        qml.CNOT([8, 7])
        qml.CNOT([13, 12])
        qml.CNOT([2, 1])
        qml.CZ([9, 3])
        qml.CZ([6, 5])
        qml.CZ([15, 14])
        qml.CNOT([8, 1])
        qml.CNOT([13, 4])
        qml.CNOT([16, 12])
        qml.CZ([11, 10])
        qml.CZ([6, 1])
        qml.CZ([15, 7])
        qml.CNOT([8, 3])
        qml.CNOT([13, 5])
        qml.CNOT([16, 14])
        qml.CZ([11, 4])
        qml.Hadamard([6])
        qml.Hadamard([15])
        qml.Hadamard([8])
        qml.Hadamard([13])
        qml.Hadamard([2])
        qml.Hadamard([16])
        qml.Hadamard([11])
        qml.Hadamard([9])
        measurements["c_(1, 1, 1)_0"] = catalyst_measure(6)
        measurements["c_(2, 2, 1)_0"] = catalyst_measure(15)
        measurements["c_(1, 2, 1)_0"] = catalyst_measure(8)
        measurements["c_(2, 1, 1)_0"] = catalyst_measure(13)
        measurements["c_(0, 1, 1)_0"] = catalyst_measure(2)
        measurements["c_(3, 2, 1)_0"] = catalyst_measure(16)
        measurements["c_(2, 0, 1)_0"] = catalyst_measure(11)
        measurements["c_(1, 3, 1)_0"] = catalyst_measure(9)
        #  End of circuit: measure q1 syndromes - cycle 0
        #  End of circuit: measure q1 syndromes 1 time(s)
        #  End of circuit: inject t into block q1 and measure syndromes
        # Start of circuit: measure q1 syndromes 1 time(s)
        # Start of circuit: measure q1 syndromes - cycle 0
        catalyst_measure(6, reset=True)
        catalyst_measure(15, reset=True)
        catalyst_measure(8, reset=True)
        catalyst_measure(13, reset=True)
        catalyst_measure(2, reset=True)
        catalyst_measure(16, reset=True)
        catalyst_measure(11, reset=True)
        catalyst_measure(9, reset=True)
        qml.Hadamard([6])
        qml.Hadamard([15])
        qml.Hadamard([8])
        qml.Hadamard([13])
        qml.Hadamard([2])
        qml.Hadamard([16])
        qml.Hadamard([11])
        qml.Hadamard([9])
        qml.CZ([6, 4])
        qml.CZ([15, 12])
        qml.CNOT([8, 5])
        qml.CNOT([13, 10])
        qml.CNOT([2, 0])
        qml.CZ([9, 7])
        qml.CZ([6, 0])
        qml.CZ([15, 5])
        qml.CNOT([8, 7])
        qml.CNOT([13, 12])
        qml.CNOT([2, 1])
        qml.CZ([9, 3])
        qml.CZ([6, 5])
        qml.CZ([15, 14])
        qml.CNOT([8, 1])
        qml.CNOT([13, 4])
        qml.CNOT([16, 12])
        qml.CZ([11, 10])
        qml.CZ([6, 1])
        qml.CZ([15, 7])
        qml.CNOT([8, 3])
        qml.CNOT([13, 5])
        qml.CNOT([16, 14])
        qml.CZ([11, 4])
        qml.Hadamard([6])
        qml.Hadamard([15])
        qml.Hadamard([8])
        qml.Hadamard([13])
        qml.Hadamard([2])
        qml.Hadamard([16])
        qml.Hadamard([11])
        qml.Hadamard([9])
        measurements["c_(1, 1, 1)_1"] = catalyst_measure(6)
        measurements["c_(2, 2, 1)_1"] = catalyst_measure(15)
        measurements["c_(1, 2, 1)_1"] = catalyst_measure(8)
        measurements["c_(2, 1, 1)_1"] = catalyst_measure(13)
        measurements["c_(0, 1, 1)_1"] = catalyst_measure(2)
        measurements["c_(3, 2, 1)_1"] = catalyst_measure(16)
        measurements["c_(2, 0, 1)_1"] = catalyst_measure(11)
        measurements["c_(1, 3, 1)_1"] = catalyst_measure(9)
        #  End of circuit: measure q1 syndromes - cycle 0
        #  End of circuit: measure q1 syndromes 1 time(s)
        # Start of circuit: measure logical x of q1
        qml.Hadamard([0])
        qml.Hadamard([1])
        qml.Hadamard([3])
        qml.Hadamard([4])
        qml.Hadamard([5])
        qml.Hadamard([7])
        qml.Hadamard([10])
        qml.Hadamard([12])
        qml.Hadamard([14])
        measurements["c_(0, 0, 0)_0"] = catalyst_measure(0)
        measurements["c_(0, 1, 0)_0"] = catalyst_measure(1)
        measurements["c_(0, 2, 0)_0"] = catalyst_measure(3)
        measurements["c_(1, 0, 0)_0"] = catalyst_measure(4)
        measurements["c_(1, 1, 0)_0"] = catalyst_measure(5)
        measurements["c_(1, 2, 0)_0"] = catalyst_measure(7)
        measurements["c_(2, 0, 0)_0"] = catalyst_measure(10)
        measurements["c_(2, 1, 0)_0"] = catalyst_measure(12)
        measurements["c_(2, 2, 0)_0"] = catalyst_measure(14)
        #  End of circuit: measure logical x of q1
        #  End of circuit: final circuit
        return measurements


else:
    pl_main = None

if not dry_run:
    with with_alarm(timeout or float("inf")):
        pl_main.compile()

times = []
measurements = []

print("Running", end="")
for _ in range(nattempts):
    begin = time.time()
    shot_results = []
    if not dry_run:
        with with_alarm(timeout or float("inf")):
            raw = pl_main()
            per_key_lists = {
                k: v.astype(int).tolist() for k, v in raw.items()
            }
            nshots = (
                len(next(iter(per_key_lists.values())))
                if per_key_lists else 0
            )
            for i in range(nshots):
                shot_results.append({
                    k: vals[i] for k, vals in per_key_lists.items()
                })
    end = time.time()
    times.append(end - begin)
    measurements.append(shot_results)

    with open("running_times.json", "w") as f:
        json.dump(times, f)
    with open("measurements.json", "w") as f:
        json.dump(measurements, f)

    print(f", t={times[-1]}", end="", flush=True)
print(" (dry-run)" if dry_run else "")

if not dry_run:
    dump_logical_observables(measurements)

[WrathfulSpatula]

Good impulse to try to be careful about selecting just the optimization options you think are appropriate, but these options likely don't work the way you think, and you're missing one that's important. If you have an accelerator, isCpuGpuHybrid and isPaged "stack" in a "transparent" way with isOpenCL. In fact, almost all constructor options "stack," which is why the interface presents them as boolean rather than an overall enum. (That's not a mistake of the design.) These options you're turning off are not mutually exclusive and should not be treated as such in any way. In fact, when isSchmidtDecompose=True (which is the option you missed that's causing the biggest problem, here), all the underlying state vector starts out as single-qubit substates, which isOpenCL on its own actually can't handle gracefully due to work-group sizing. That's why you hit this edge case.

Try this:
@qml.qnode(qml.device("qrack.simulator",wires=17,isSchmidtDecompose=False), shots=num_shots)

I just returned 175/25 split on 0/1 on the first try, with that. You might also want environment variables QRACK_MAX_PAGING_QB=-1 and QRACK_MAX_CPU_QB=-1 to force near-Clifford simulation, so long as you know your case will only have about 30 or fewer non-Clifford phase gates, typically.

That's the simplest way to get at the case you actually want. Writing it out entirely explicitly, argument-by-argument, that actually means this:

• isTensorNetwork=True ("Stackable," optimizes the circuit diagram and subsets it by light-cone on local operations)
• isSchmidtDecompose=False (Stackable, but you usually don't want it for near-Clifford, and the real problem might be that, with your other settings, OpenCL kernels were in an edge case where group size was larger than your smallest subsystems.)
• isStabilizerHybrid=True (Stackable, and this is where near-Clifford simulation happens, with conservative fallback to state vector mostly controlled by those environment variables given above)
• isCpuGpuHybrid=True (Stackable: there's no drawback to using this with OpenCL or CUDA. It'll also prevent that group-size-based edge case when isSchmidtDecompose=True)
• isPaged=True (Stackable: again, there's no drawback, here. It's for simulation coordinated across more than 1 of the 4 maximum memory allocation segments in your GPU at a time, if you have an NVIDIA GPU, and it also coordinates multi-GPU. But it just won't engage to cause any kind of overhead, if you're below single-segment GPU memory allocation limits anyway.)
• isOpenCL=True (Stackable: yes, use the default accelerator as appropriate, too, similar to the last two options)
• isBinaryDecisionTree=False (This switches out state vector, at "layer stack" base, to use quantum binary decision diagrams, QBDD, which are CPU-only in Qrack, so it will "stack" only with everything except isCpuGpuHybrid, isPaged, isOpenCL, and isSparse, which end up ignored.)
• isSparse=False. (This is a different CPU-only method, similar to QBDD in that regard. It "kicks-in" when isOpenCL=False and isBinaryDecisionTree=False.)

But everything from isStabilizerHybrid "upward" in this list "stacks" with everything else. Usually, you want to disable the Schmidt-decomposition layer specifically only for near-Clifford.

[sergei-mironov]

Hihi! In this post I put the results of the same T|+> experiments for QEC distance 3 and 5 where
QRack is set like in this snippet:

@qjit
@qml.qnode(qml.device("qrack.simulator",wires=...,isSchmidtDecompose=False,), shots=1)
def main():
  ...

The below archive contains several folders which correspond to the runs of the experiment. Each folder contains runnable pennylane_runner.py files similar to the one in the initial post. Also the combined stdout and stderr outputs are in the output.log file. Other files contain the remaining details that might be unnecessary. Also, the report1.md contains a short description and links which should work from a comprehensive Markdown viewer.

report1.tar.gz

In short, the errors are:

• ValueError: QEngineOCL::CopyStateVec argument size differs from this! for the Laptop d=3 case,
• RuntimeError: PennyLane runner failed with return code -14 for other cases (with different error codes)

The first one makes me think that something might be wrong with the compatibilities, could you please check?

[WrathfulSpatula]

I'm bemused to see that you're on Intel integrated graphics:

Let me explain: we devote significant energy to making sure that Qrack supports integrated graphics, and we're one of the only major simulators who offer any such support, to my knowledge. However, the drivers are often just not robust. There are some steps you can take to help improve their stability, though.

As the README states, you probably want to set environment variable (as on Linux or Mac) export PYQRACK_HOST_POINTER_DEFAULT_ON=1. With most or all discrete GPUs (like NVIDIA devices), they have their own internal RAM, so we allocate for GPU simulation on that RAM on the device; the Intel UHD, however, actually uses something like "direct memory access" (DMA) on general system RAM (heap) and does not have its own dedicated on-board RAM banks. So, this environment variable I just gave above respects the fact that device RAM is host RAM, which can improve performance and prevent expensive duplication of the state vector. This is a good first step to try.

One of the tricky aspects of working with integrated graphics is that, since they're using general heap, they have access to more memory (even more than a recent high-end NVIDIA consumer-line GPU) than practical throughput to utilize it. Because Qrack tries to automatically detect your maximum RAM available for state vector simulation as a fallback for near-Clifford, it can end up expecting to be able allocate far more RAM than an integrated graphics can reasonably handle without practical failure to allocate. So, what can happen with near-Clifford circuits is that Qrack mistakenly detects that it should be able to allocate something like 31-to-33 qubits worth of state vector, hence it caches more T gates, but, when it encounters even more T gates than that later in the circuit, it tries to flush the T-gadget cache at around the ~33 qubit mark, for auxiliary state vector, and the Intel UHD just can't handle it. You can check the maximum OpenCL allocation with the clinfo command line tool, but I believe that recent Intel UHDs allow a maximum single allocation at once of 2 GB at a time. On a NVIDIA device, you can fit up to 4 such maximum allocation segments, but the UHD often struggles with any more than 1.

So, the upshot of all that is, try export QRACK_MAX_PAGING_QB=28 or 27. This controls multi-segment allocation limits, but the slightly different QRACK_MAX_PAGE_QB controls single maximum allocation, and maybe you can set that to 27 or 28 and raise the global limit for QRACK_MAX_PAGING_QB 1 or 2 qubits. (This requires clinfo and knowledge about your system to solve exactly, but you can likely figure it out by trial-and-error somewhere about that 27-to-28 qubit mark.) Hope this helps!

[WrathfulSpatula]

By the way, another option is just turn off OpenCL or install POCL to use your CPU as an accelerator for kernels. We have a lot of love for the Intel HD, but, often, you're better off running on the CPU itself.

[WrathfulSpatula]

One more point to disambiguate: even if you tell Qrack to use 17 logical qubits, the ~32-qubit mark I'm talking about is for auxiliary qubits, where, in near-Clifford simulation, Qrack will try to use as many of these as it detects it can fit, or whatever you tell it your desired limit to be, with those environment variables I explained above.

[sergei-mironov]

Hi @WrathfulSpatula . First of all, thank you for your very detailed response. I must say sorry for leading you to the conclusion that the Intel UHD is our main simulation device. Let me explain the situation. My ultimate goal is to learn how I can run near-Clifford experiments like the above T|+> using as high code distance as possible. So I am synthesizing the code like I sent initially for d = 3, 5, 7, etc. (the number of qubits goes like 17, 49, 99, ...). I expect to learn a regime where I would be able to run high code distances with the degraded precision and then estimate how bad this degradation is.

Apparently, I have two machines at my disposal: my own laptop running the Intel UHD card and a hired Runpod.IO container offering a large set of GPUs (and slow conventional RAM btw), where I often pick RTX 3090 for price reasons.

Thanks to your answer I see that the real-GPU setup would be very different from the one for Intel cards (I do understand the point about hardware pointers and DMA channels). But since we consider the Nvidia one as a more scalable solution, let's focus on the Runpod-style machine equipped with an Nvidia RTX 3090 GPU further in this thread!

---

Now let me provide an update on the issue. I ran a grid-search among the 128 combinations of QRack parameters applied to a distance 5 code (49 qubits) of a similar nature as above. The experiment took place on the Runpod machine I also described above. The main loop of the experiment looks like this

num_attempt = 5 # Later changed to 1 (see below)
num_shots = 1 # I put 1 to check for basic runnability, 100 for a statistical test.
all_times = []
@qjit
@qml.qnode(
    qml.device("qrack.simulator",wires=49,
        .... # varying parameters
    ),
    shots=num_shots
)
def pennylane_qrack_kernel():
    ....

for i in range(num_attempt):
   t0 = time()
   results = pennylane_qrack_kernel()
   t1 = time()
   all_times.append(t1 - t0)
   save_to_file('running_times.json', all_times)

Initially I saw no success (each and every run returned an error), but when I took a closer look I noticed that in some cases I had both errors and measurement results. I discovered that sometimes QRack runs the first attempt successfully and crashes on a subsequent one. In total, I counted four such semi-successful cases. For these I conducted a statistical test (with num_shots equal 100) and found that in two cases I had good results and nice timings. In another two cases I got statistically improbable results, but one needs to check them for coincidence.

Apparently, the fastest case goes like:

os.environ["QRACK_MAX_PAGING_QB"] = "-1"
os.environ["QRACK_MAX_CPU_QB"] = "-1"
@qjit
@qml.qnode(
    qml.device("qrack.simulator",
        wires=49,
        isTensorNetwork=False,
        isSchmidtDecompose=False,
        isStabilizerHybrid=False,
        isOpenCL=True,
        isPaged=False,
        isCpuGpuHybrid=False,
    ),
    shots=100
)
def pl_main():
    ...

... which might exactly match your earlier suggestion depending on the default values for False options!

I am attaching an archive with details on all four cases, with comments on what worked and what did not. Please let me know what you think about the second attempt failures in the Basic correctness tests section! Could we fix it? Also, let me know what do you think about the low precision cases. Both don't look as totally incorrect, maybe we can do something? What setup do you think should scale better for d>=7 ?

Details:

report2.tar.gz

[WrathfulSpatula]

49 qubits! Depending on the exact case, that might be a tall order. So, actually, if you do this in a script:

os.environ["QRACK_MAX_PAGING_QB"] = "-1"
os.environ["QRACK_MAX_CPU_QB"] = "-1"

Python will not pick up the environment variables for Qrack, or any C++-based shared library wrapper. But I assume you're doing this either in the CLI shell or an environment like Jupyter, rather. However, with these settings:

@qjit
@qml.qnode(
    qml.device("qrack.simulator",
        wires=49,
        isTensorNetwork=False,
        isSchmidtDecompose=False,
        isStabilizerHybrid=False,
        isOpenCL=True,
        isPaged=False,
        isCpuGpuHybrid=False,
    ),
    shots=100
)

those environment variables won't do anything at all, because you've disabled isPaged and isCpuGpuHybrid, and you're not using CPU-based simulation. You've selected just state-vector-based OpenCL simulation on a single GPU memory segment (out of 4). Whatever results you might have produced in this case, realize that, with these simulation method options, you would need 4 PB of GPU memory (in a single maximum allocation segment, so, rather, the situation is already unrealistic, but you'd actually need a GPU with something more like 16 PB of total RAM with a 4 PB segment size, but anything follows from a contradiction in reasoning).

You might be seeing some semi-plausible results, with default settings rather like this:

@qjit
@qml.qnode(qml.device("qrack.simulator", wires=49), shots=100)

because of Qrack's automatic circuit elision (ACE). When isSchmidtDecompose=True is engaged (as by default), if the first-principles fidelity estimate stays high enough, or if you disable the "fidelity guard," Qrack will actually try to automatically segment your simulator state into subsystems no larger than your QRACK_MAX_PAGING_QB setting and act semi-classical gate shadows across the separable subsystem boundaries, to approximate the effect of "inter-patch" gates without exceeding system resources. This will decrement a first-principles, worst-case fidelity estimate that unfortunately can't be queried from the PennyLane Catalyst API (though it's sim.get_unitary_fidelity() in pyqrack), but, honestly, it tends to be quite pessimistic, anyway. We've found that, for 4 separable quadrant patches, nearest-neighbor random circuit sampling (RCS) on a rectangular 2D chip topology validated on cross-entropy benchmark (XEB) fidelity up to at least 32 qubits (i.e., at least as high as we can directly validate at desktop scale) returns ~50% XEB independent of width and depth, interestingly, except for specifically only perfect-square qubit dimensions, like 2-by-2, 3-by-3, 4-by-4, etc., which return close to 0%. Fully-connected RCS does not share this XEB performance characteristic, but it keeps some XEB through relatively low circuit layer depths. So, if results look "plausible," that might be automatic ACE approximation.

There are some other approaches you can try. Honestly, they're probably not going to work. (Like I said already, state-vector simulation would require ~4 PB.)

For alternative ideal simulation methods, your options are near-Clifford or quantum binary decision diagrams (QBDD).

For near-Clifford, in CLI,

export QRACK_MAX_PAGING_QB=-1
export QRACK_MAX_CPU_QB=-1

then, in the same shell,

@qml.qnode(qml.device("qrack.simulator", wires=49, isSchmidtDecompose=False), shots=100)

This is "exact," but it's very slow. I tried this before, with one of your smaller circuits, and we hit your time-out condition. In theory, if it weren't too slow, with an RTX 3090's memory limits, you could handle up to about 30 qubits (give-or-take) of non-Clifford phase gates, like your T gates. Beyond that many non-Clifford gates, you would run out of memory, even in principle.

For QBDD (which is best for low-entanglement scenarios, sparsity, or low-depth, so probably not your situation):

@qml.qnode(qml.device("qrack.simulator", wires=49, isBinaryDecisionTree=True), shots=100)

Then, there are approximate approaches. I explained "ACE," already. To make sure the "fidelity guard" doesn't kick in, you use this in the shell session:

export QRACK_DISABLE_QUNIT_FIDELITY_GUARD=1

Then, just use default device settings. Or, you could try sparse simulation, which, working in concert with our Schmidt-decomposition layer, behaves somewhat similarly to MPS with a slight bias towards measurement basis (likely). You can set a maximum memory limit before truncation of low-probability amplitudes starts to occur:

export QRACK_SPARSE_MAX_ALLOC_MB=[n]`

[n] is whatever you're willing to wait for, for MB of memory limit, because this is CPU-based and slow. I've had surprisingly good results on fully-connected RCS results around roughly the 20 qb range with a setting as low as even 1 MB. And you can also set a probability floor to "proactively" truncate small amplitudes based on probability itself, instead of just "reactive" probability-ranked priority to respect the memory limit:

export QRACK_SPARSE_TRUNCATION_THRESHOLD=0.001

(just using 1 part in 1000 probability as an example, here). Then, disable OpenCL and choose sparse simulation on CPU:

@qml.qnode(qml.device("qrack.simulator", wires=49, isOpenCL=False, isSparse=True), shots=100)

Beyond that, I can suggest trying an alternative simulator, like Quimb for tensor network trees, or the alternative near-Clifford approach in IBM Qiskit. Like I said earlier, we don't look to replace the work and approaches of others, just to supplement the rest of the open-source quantum software ecosystem, to work better together. But, the simple fact is, you can approximate 49 qubits to any necessary depth, to some indeterminate quality level, but naive expectation for ideal simulation requires something like 4 PB of memory. That's just how it is.

[sergei-mironov]

os.environ["QRACK_MAX_PAGING_QB"] = "-1"
os.environ["QRACK_MAX_CPU_QB"] = "-1"

Python will not pick up the environment variables for Qrack, or any C++-based shared library wrapper.

It was actually an issue in my runs, the variables did not worked. Somehow these instructions managed to get into the main program from the external runner where they were supposed to be. Now it is fixed, thank you! I have re-run the grid-search (also I set number of attempts to 1 by default to workaround the multiple-attempt crashes). But for the distance 5 version (49 qubits), the outcome was the same set of four runnable configs as I posted. FYI I am attaching the full results for 128 configs report3.tar.gz (the .md file classifies them using tables, links should work locally.)

with these settings:

@qjit
@qml.qnode(
    qml.device("qrack.simulator", wires=49, ... isOpenCL=True ...  ), shots=100
)

those environment variables won't do anything at all, because you've disabled isPaged and isCpuGpuHybrid, and you're not using CPU-based simulation. You've selected just state-vector-based OpenCL simulation on a single GPU memory segment (out of 4).

In my new run I see the same config works again. Yep, I got it that it effectively runs a state-vector, this is not what I want for 49 qubits. I hoped to end up with an approximate method instead. Now I see that the near-Clifford might also work for me.

Did I understand you correctly that the near-Clifford simulation method is an exact method in QRack? Is it still slow even in the case of exactly one T gate, as in all my programs? But it is a GPU-based, right? Note that while the code distance (d=3/d=5) parameter changes the number of Clifford gates, there is only a single T-gate in my programs, it is applied to the middle qubit of the qec code square.

I will continue experimenting with methods you have described.

I wonder might Qrack treat the absence of a parameter differently from True or False? If this is the case, I might need to include the third "not set" parameter option for my grid search runs.

[WrathfulSpatula]

Did I understand you correctly that the near-Clifford simulation method is an exact method in QRack? Is it still slow even in the case of exactly one T gate, as in all my programs? But it is a GPU-based, right? Note that while the code distance (d=3/d=5) parameter changes the number of Clifford gates, there is only a single T-gate in my programs, it is applied to the middle qubit of the qec code square.

Yes, our near-Clifford method is exact, if other approximation methods like "ACE" aren't mixed in. For each T gate that can't be simplified away through local buffering and gate fusion (like your single T gates), we act a "reverse phase-injection gadget" as described in Appendix A in https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.3.020361. The general idea is that, with such a gadget constructed at the time the T gate is demanded, but with measurement deferred, can you continue along with more Clifford gates in an exact manner so long as you can guarantee post-selection success on the gadget ancilla, before terminal or mid-circuit measurement, which can always be guaranteed in simulo, in this case. (Hence, it's exact.) My biggest contribution to this general idea, in the Qrack implementation, was to notice that we could calculate single amplitudes at a time in the large near-Clifford circuits by finding the overall post-selection |0> state amplitude with an auxiliary state vector simulation based on just the ancillary qubit states here in QStabilizerHybrid. This is reasonable for single amplitudes, but sampling still requires calculating half of all amplitudes, average case, so you're still left waiting on ~2^48 repetitions of the single-amplitude calculation. (That's painful.) This method is distinct from IBM's approximate approach based on "stabilizer rank," though, and I have absolutely no qualms about suggesting our users try that other method for sampling in cases Qrack can't handle.

I wonder might Qrack treat the absence of a parameter differently from True or False? If this is the case, I might need to include the third "not set" parameter option for my grid search runs.

No, it does not. Each constructor parameter defaults to either True or False, and "the middle is excluded."

You say you only use a single T gate: this might sound heretical, but have you considered an approach like trying a kind of branch cut on pure stabilizer simulation by applying either S gate or no-operation in its place, then somehow stitching back together the two pure-Clifford branches? The effect of the single T gate likely can't propagate chaotically downstream to ruin the entire behavior. We demonstrate a similar idea in example scripts like this, and the others in the folder that contains that, in the examples repository.

[WrathfulSpatula]

In fact, hear me out: the "reverse phase-injection gadgets" method would already work great, for that kind of "branch cut," with only a tiny bit of modification, to linearly superpose the branch cuts by relative phase angle of T or S gate, by the end of the circuit. I'll hack up an alternate sampling mode, and I'll validate on small circuits against Qiskit Aer as gold standard, in that same examples script folder. If it validates, happy day: everyone can use it, at least as an alternative approximation method. But the idea is simply that independent stabilizer states are linearly superposable in arbitrary ratios just like any other proper ket. (This is the basic idea behind IBM's "stabilizer rank" method, but this approach I propose is still distinct.) I'll give it a try, and I should have something later today, or by tomorrow.

[WrathfulSpatula]

Update on this: those changes are published. I just ran your 49-qubit circuit (approximately) with "ACE," in about a minute, and I'm currently waiting on approximate near-Clifford, but it takes quite a bit longer. For approximate near-Clifford, you'll want to experiment with the following options and settings:

@qjit
@qml.qnode(qml.device("qrack.simulator", wires=49, isTensorNetwork=False, isStabilizerHybrid=True, isSchmidtDecompose=False), shots=num_shots)

Before calling any scripts, set the following environment variables in your shell:

export QRACK_MAX_PAGING_QB=-1
export QRACK_MAX_CPU_QB=-1
export QRACK_NONCLIFFORD_ROUNDING_THRESHOLD=1.0
export QRACK_USE_APPROX_NEAR_CLIFFORD=1

Then, run your scripts in the same shell. This is potentially orders of magnitude faster than exact near-Clifford simulation, but I can't say that it will be "fast," in your specific case. The idea, here, is that we stochastically decide whether to replace your T gate with either S (or adjoint S, as appropriate) or no-operation, and the pure-Clifford behavior otherwise always "lines up" on the closest Clifford circuit to within a gadget angle of the same magnitude of a T gate or less, whatever the closest phase angle quadrant is by S gates. The probability of S-vs.no-op' stochastic correction is weighted by the fraction of S-gate angle for each buffered non-Clifford angle in each ancilla (so it's 0.5/0.5 for any T gate).

I think that's the best I can offer, for now. Your QEC case is even particularly expensive, beyond simple near-Clifford case considerations, because of your mid-circuit measurement (which really makes this "brutal"). This truly might be a good case for the Qiskit Aer near-Clifford method, otherwise. (But thanks for motivating me to think about whether we could come up with any further improvements on approximate near-Clifford cases, because I realized we could!)

[WrathfulSpatula]

Here's a quick theoretical note: what's definitely lost in this stochastic approximation is interference effects between multiple non-Clifford phase gates. What isn't clear to me is whether anything else is lost for accuracy or precision, i.e., any effect that can't be modeled as interference between just the non-Clifford gates. So, I'll speculate, for just one non-Clifford gate at a time, it might be possible to fully recover sampling statistics with this weighted stochastic correction approach. (For multiple non-Clifford gates, we effectively model the interaction with each individual gadget ancilla as "separable," which it typically won't be, which is why I say we lose "interference" between them.) There are some caveats that come to mind, though:

• Terminal sampling without mid-circuit measurement is the optimal case. The single T gate doesn't get injected until the point you actually request sampling, so the gadget remains consistent with any (non-Clifford) phase injected until the last step, right at sampling.
• Mid-circuit measurement ruins this, because, in any case, you would need to inject all buffered non-Clifford gates early. (I'm not even sure I've covered this in the library: I actually can't recall what happens to ancillary gadgets in cases of mid-circuit measurement, but I'll look into this.)
• You could potentially correct the statistics for mid-circuit measurement by using the equivalent of PennyLane's "one-shot" mode: this re-simulates the entire circuit top-to-bottom for every single measurement shot (just like a hardware quantum computer must).
• However, if you use one-shot mode ("weak sampling condition"), then you can't calculate statistics for particular intermediate collapse results without post-selection (but, again, neither can a genuine hardware quantum computer).

So, you see, the real "monkey wrench" that's left in the machinery is mid-circuit measurement. I'm still thinking about ways around this, but no one said every arbitrary quantum circuit was "easy."

[sergei-mironov]

I think that's the best I can offer, for now. Your QEC case is even particularly expensive, beyond simple near-Clifford case considerations, because of your mid-circuit measurement (which really makes this "brutal"). This truly might be a good case for the Qiskit Aer near-Clifford method, otherwise

Thank you, I'll check the new version soon and come back!

Let me share thoughts about Qiskit: Indeed, we already found out that Aer works best on the circuits from this issue: it runs our d=5 circuit in ~2 sec per shot while QRack did it in ~4 sec in our experiments (plus some workarounds were required and it went into a wrong regime as we found). Still, we probably cant go with Qiskit because of the nature of the Qec workflow which requires calling a classical "error syndrome decoder" routine having all these measurements as input. Importantly, we also need decoder results before the program is finished so we can apply corrections! This is not the part of this issue, but in a full experiment PennyLane+Catalyst (and thus QRack) would allow us to wire this kind of callback into the main simulation loop (we hope so!), while in Qiskit it would be hard to do so without patching the package.

Some theory I learned says that I probably should want a Monte-Carlo near-Clifford implementation so I can explode T-gates into a sums like $T = p_1 Pauli_1 + p_2 Pauli_2 + ...$ where $p_i$ are probabilities so I can sample from branches each shot rather than simulate this explosion and sum results together. I am not sure however how this combines with the measurements. Yet I was not able to find a suitable implementation beside few proof-of-concepts.

[sergei-mironov]

You could potentially correct the statistics for mid-circuit measurement by using the equivalent of PennyLane's "one-shot" mode: this re-simulates the entire circuit top-to-bottom for every single measurement shot (just like a hardware quantum computer must).
However, if you use one-shot mode ("weak sampling condition"), then you can't calculate statistics for particular intermediate collapse results without post-selection (but, again, neither can a genuine hardware quantum computer).

I just grasped these last points. This looks exactly what I need I think. I'd say we indeed need top-to-bottom sampling, something which would emulate quantum computer's operational behavior, in order to run the error correction routine. But did it work differently up to now despite giving me all the mid-circuit measurements as bits rather than as statistics? I think I need to learn the difference between weak and strong sampling conditions.

[WrathfulSpatula]

But did it work differently up to now despite giving me all the mid-circuit measurements as bits rather than as statistics?

Yes, and this is the difference between "weak" vs. "strong" sampling conditions. In any case with probabilistic mid-circuit measurement, you can think of actual quantum hardware as effectively weak due to wave-function collapse. What "strong" means is that you have something like a Schrödinger-picture (classical and fully-specified) representation of the in-flight unitary state, and you can simulate measurement on it without collapsing the state. (Hence, sampling becomes inexpensive to simulate, since you don't need to recompute the unitary portion of the circuit.)

It's not strictly impossible, at all, to produce statistics for final results with mid-circuit bits. However, the only way you can do it on actual quantum hardware is to subset by mid-circuit bits on single-shot end results, which is effectively post-selection. If you're trying to cover a specific bit sequence, you just have to repeat until random success, which is effectively exponential costly, like post-selection in general on a hardware quantum computer.

Note also that some cases of classical simulation must use top-to-bottom one-shot to be statistically correct. Notably, noisy simulation with Qrack in PennyLane Catalyst is based on "weak Pauli noise." In that case, correct sampling statistics are guaranteed only for one-shot. The trade-off is that, in this case, with weak sampling, we can still simulate as state vector instead of density matrix, hence we can simulate about twice as many qubits in circuit width compared to "strong" density matrix.