3
$\begingroup$

I am currently trying to compute logical error rates for the surface code using Stim's detector error models and PyMatching for different distances and noise strengths.

tl;dr : What is the best strategy for parallelising this process over different cores?

In doing so, I am trying to optimise the number of cores at my disposal by parallelising the different computations across multiple cores. My current workflow leverages joblib for this task

from joblib import Parallel, delayed, parallel_backend

distances = [...] # A set of distances
noise_strengths = [...] # A set of noise strengths
num_shots = ... # Some number of shots to sample
num_cycles = ... # Number of cycles over which to perform syndrome extraction

threads = ... # Maximum number of threads allowed
jobs = ... # Maximum number of jobs allowed

with parallel_backend("loky", inner_max_num_threads=threads):
        joblib_results = Parallel(n_jobs=jobs)(delayed(compute_logical_error_rate)(distance,noise,n_shots) for distance in distances for noise in noise_strengths)

where the logical error rate is computed using stim and pymatching in the following way

def compute_logical_error_rate(distance,noise,num_shots,num_cycles):

       sc_circuit = surface_code(distance,noise,num_cycles) # Define circuit
       sc_model = sc_circuit.detector_error_model(decompose_errors=True) # Detector error model
       
       sampler = sc_circuit.compile_detector_sampler() # Define sampler
       syndrome, measured_obs = sampler.sample(shots=num_shots, separate_observables=True) # Extract samples

       matching = Matching.from_detector_error_model(sc_model) # Matching graph
       expected_obs = matching.decode_batch(syndrome) # Perform decoding

       num_errors = np.sum(np.any(expected_obs != measured_obs, axis=1)) # Compute errors
       logical_error_rate = num_errors/num_shots # Compute logical error rate

       return logical_error_rate

This workflow allows me to assign a different simulation to each core, i.e. each core will perform the computation of the logical error rate for a given distance and noise strength. However, upon increasing the number of samples that I take, defined as num_shots, I start running into memory problems. I have identified the reason to be the syndrome variable, which is a numpy array whose size scales linearly with number of samples considered syndrome.shape = (num_shots,...). The way I am parallelising naturally results in storing this variable multiple times (once for each job/core), thus leading to a memory overload.

I am interested in understanding how I can best leverage multiple cores to make this whole set of simulations more efficient. Would it make more sense to tackle one job at a time (fixed distance and noise strength) and parallelise the sampling across different cores? (If so, how can that be done?) - There seems to be a comment about this approach on the original Stim paper, but I do not see how to replicate this behaviour with my code.

$\endgroup$

1 Answer 1

4
$\begingroup$

(Sounds like you want sinter. In particular you should check out sinter.collect. It's for exactly this task, but supports additional features such as the ability to save progress to a file so it can resume work later if your machine lost power.)

First, you should divide up the work into batches. Stim is optimized to work best when you keep asking for batches of shots of the same size, where that size is a multiple of 256.

Second, you should ask for bit packed data. This makes the data 8x denser, so there's a lot less copying.

Third, you should return how many samples were taken and how many errors were seen, instead of their ratio. This allows you to run the same circuit on multiple threads and easily combine the result afterwards by adding the returned numbers of shots and errors.

def count_logical_errors(distance,noise,num_shots,num_cycles):
   sc_circuit = surface_code(distance,noise,num_cycles) # Define circuit
   sc_model = sc_circuit.detector_error_model(decompose_errors=True) # Detector error model       
   sampler = sc_circuit.compile_detector_sampler() # Define sampler
   matching = Matching.from_detector_error_model(sc_model) # Matching graph

   shots_left = num_shots
   num_errors = 0
   while shots_left > 0:
       batch_size = min(1024, shots_left)
       shots_left -= batch_size
       syndrome, measured_obs = sampler.sample(
           shots=batch_size,
           separate_observables=True,
           bit_packed=True,
       ) # Extract samples

       expected_obs = matching.decode_batch(syndrome) # Perform decoding

       num_errors += np.sum(np.any(expected_obs != measured_obs, axis=1)) # Compute errors

   return num_errors, num_shots

Incidentally, your method will be much more flexible if you pass in the circuit instead of the parameters used to make the circuit.

$\endgroup$
4
  • $\begingroup$ Thank you very much for the answer! Indeed Sinter did the job - I do have two follow up questions on it: (1) How does the choice of max_error affect the final logical error rate obtained? (2) I noticed that before leveraging all the num_workers (e.g. 16), sinter does some pre-processing over a notably smaller amount of cores (e.g. 2), for a large system size the pre-processing naturally becomes longer, should I be enabling something else to facilitate a speed up of this pre-processing? $\endgroup$ Feb 27 at 3:26
  • $\begingroup$ @EzequielRodriguezChiacchio max_error shouldn't change the logical error rate you obtain in expectation; it changes how much variance there is in the estimate. Because the variance of a binomial distribution drops as you get close to p=0, the number of errors you observe basically corresponds to how wide the uncertainty region looks on a log plot regardless of whether it's at 10^-8 or 10^-3. For sinter preprocessing there's not much to be done except maybe include the detector error models in the task from the start. $\endgroup$ Feb 27 at 5:17
  • $\begingroup$ That makes sense, thank you. My last question would be, is it possible to recover the samples from the Detector Error Model when using sinter.collect? - I am thinking about the list of detectors that did and did not fire up, as when calling circuit.compile_detector_sampler().sample() in Stim. $\endgroup$ Feb 29 at 4:47
  • $\begingroup$ @EzequielRodriguezChiacchio No, it doesn't support that. You'll have to do the sampling yourself for that. $\endgroup$ Feb 29 at 6:11

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.