Error Generator Propagation Module

[coreyostrove]

This PR introduces version 1.0 of a new pyGSTi module for the efficient approximate simulation of noisy Clifford circuits using a technique based on the analytic propagation of elementary error generators. This is a rather major update and will likely be challenging to review just due to the sheer amount of code that has been added or touched, so I will do my best to summarize notable changes and additions below. Where relevant I'll also discuss any notable architectural/design considerations that may be relevant for reviewers.

Overview of Changes

Error Generator Propagator Class

The central new class for the error generator propagation module is ErrorGeneratorPropagator. This class houses the main methods responsible for performing/orchestrating the propagation of elementary error generators to the end of a circuit, and managing any post-propagation processing related to producing effective end-or-circuit error generators (e.g. iterative application of the BCH approximation). This class wraps/contains a reference to a corresponding model, and responsibility for generating sparse representations of the per-circuit layer error generators (i.e. error generator coefficient dictionaries) needed to instantiate the propagation computation is delegated to this model. As such, all of the functionality described in this PR should be compatible with any error-generator parameterized model, or otherwise any model able to return per-circuit-layer error generator coefficient dictionaries).

Stim is a dependency now

The functionality of this module is inextricably built on top of stim as a backend for fast stabilizer tableau related computations, so with this module is now added to the requirements list for pyGSTi.

LocalStimErrorgenLabel Class

To help manage metadata and other functionality which is helpful in the context of error generator propagation we've introduced a new ElementaryErrorgenLabel subclass called LocalStimErrorgenLabel. In many ways this class is similar to the existing LocalElementaryErrorgenLabel class, both store their basis element subscript labels in a 'local' format for example. I.e. as full n-qubit pauli strings, but for LocalStimErrorgenLabel these are now instances of stim.PauliString objects instead of python strings (heavily reducing the overhead of conversion back and forth throughout the propagation code). This class is also responsible for one of the core subroutines in propagation, being delegated the responsibility of mapping its basis element labels to their new value under conjugation by a Clifford operation. To simplify interoperability between this new class and other existing parts of pyGSTi which don't presently natively support these new labels (perhaps at some point we'll change that) this class has built-in methods for casting itself to instances of LocalElementaryErrorgenLabel or GlobalElementaryErrogenLabel, as well as for casting from instances of these classes.

New errgenproptools module

To support the functionality of the ErrorGeneratorPropagator class, as well as extend those capabilities in other directions, there is a new module in pygsti.tools with a large number of related utilities. Here is an overview:

• Functions for analytically computing the commutators and compositions of elementary error generators. These underpin the efficient analytic computation of higher-order BCH approximations. By volume these functions are by far the biggest ones in this module (they are essentially massive complicated table lookups and due to the complexity my main goal for v1.0 was correctness and not conciseness), so apologies to anyone attempting to actually read that code end-to-end. Additionally there are functions for iteratively performing multiple (nested) analytic error generator compositions.
• A function for computing the BCH approximation for the composition of two exponentiation error generators. This function takes as input two error generators written in their sparse (coefficient dictionary) representation and manages the logic for applying and combining all of the needed commutators (including iterative commutators) for combining these error generators at the specified order. At the moment (and for the foreseeable future) this supports up to the 5th order BCH approximation. Ain't no one got time to implement Dynkin's formula.
• Functions for efficiently and analytically computing corrections to the output probabilities of computational basis state readout, and expectation values of pauli observables. This is based on performing an analytic taylor series approximation to the matrix exponential of a sparse error generator and then looking at the linear sensitivities (the code for computing these sensitivities is likely to be of utility in other contexts as well, looking at you idle tomography) of the probability/expectation value in question. This approximation does (in principle) support the computation of arbitrary order corrections to the probabilities/expectation values.
• Most of the core functions in this module have mirrored implementations which perform the various computations purely numerically. These versions of the function are not scalable, of course, but are used heavily in the testing infrastructure for this module.

Lots of supporting backend changes

In no particular order (I'm just going to go top to bottom through the diff summary and may skip any de minimis or obvious ones). This section will probably be the longest of this write up, mostly because backend changes touch a lot of pyGSTi code and are the hardest to understand without context.

• Basis performance improvements. Previously certain frequently accessed basis properties corresponding to properties that should be static were recomputed on the fly (e.g. the dimension, or the component bases). This slowed down certain subroutines (particularly hashing). We now precompute these quantities at initialization time (the actual computations are cheap, the slow part was redoing them thousands of times) and cache their values internally. There is also some rudimentary caching added for the basis transformation matrix construction (just a fairly small LRU cache) which significantly reduces the overhead of converting multiple matrices from one matrix to another. (It looks like I must have also found an edge case in the original hashing function but I don't recall the details on that anymore.)
• CompleteElementaryErrorgenBasis and ExplcitElementaryErrorgenBasis have been notably extended. This includes:
  • significant performance improvements by adding in caching for error generator basis elements, labels, and other expensive to compute attributes (previously these were being reconstructed each time a user asked for them).
  • Reorganization/splitting of certain methods. E.g. previously there was only a single function for jointly returning the error generator matrices and the supports together, but we rarely need the support part, so that is now available separately (the old functions are still available for compatibility).
  • More flexibility for CompleteElementaryErrorgenBasis. Users can now specify separate maximum weights for all 4 sectors (previously you could only separately control the weights of the H block and the SCA block as a whole). There is also now an option for controlling the format of the elementary error generators used by the basis (either local or global style).
  • Miscellaneous bugfixes and (nearly) complete documentation coverage.
• LocalElementaryErrorgenLabel and GlobalElementaryErrorgenLabel:
  • Hash caching for better dictionary lookup performance (like other pyGSTi label objects these are treated as static).
  • Minor bugfix for casting.
  • Lots of additional documentation.
• StateSpace is now effectively static. Certain attributes that were previous exposed as part of the public API and potentially mutable are now protected as properties. This change was made for two reasons: 1. We were pretty much already treating state spaces as immutable in practice, I didn't find any instance in pyGSTi where we do inplace state space manipulation and all of the state space methods that involve mutation already returned new objects and avoided inplace manipulation. 2. Enforcing this means we can start caching hashes to improve performance, which in-turn improves the performance of certain operations involving the Basis class.
• The Circuit class now has two new methods for conversion from a pyGSTi circuit (clifford only) to either a list of stim.Tableau
  objects (one per circuit layer), or a single stim.Tableau corresponding to the combined clifford operation for the overall circuit. There have also been requests for conversion in the other direction (stim to pyGSTi), so that is on the roadmap for a future version.
• First-class support for LocalElementaryErrorgenLabel. Previously most functions which involved working with elementary error generator labels assumed that these labels were instances of GlobalElementaryErrorgenLabel and these were the labels used internally in many/most error generator related model members, for example. The exception to this is the code which actually constructs elementary error generator matrices which requires basis element labels which are in the local-style and so casted the error generator coefficients to LocalElementaryErrorgenLabel. Anyhow, rather than go on and on, all of this conversion back and forth in different contexts and doing lots of casting (sometimes multiple casts for a single operation) was a real performance bottleneck when instantiating error generator dictionaries for a circuit layer when instantiating propagation. To address this I've rewritten these code paths to allow working interchangeably with both label types, and all of the relevant methods now have the option of specifying which type of label you'd like returned. I've also added some caching to reduce the overhead of switching back and forth between the requested label types in different contexts.
  -LindbladErrorgen and LindbladCoefficientBlock:
  • For LindbladCoefficientBlock the performance of constructing error generator coefficient dictionaries was improved, and caching was added for these so that the dictionaries are only reconstructed when block data has been changed. Previously it would reconstruct this from scratch each time it was requested which was a big bottleneck..
  • For LindbladErrorgen aside from changes related to adding first-class support for local labels as described in the previous point, we now finally have nearly completely coverage and up-to-date documentation for this module! (GPT helped me out a bit).
• lindbladtools: Faster dense elementary error generator construction. This is achieved two ways.
  • First, these functions always return these matrices in the standard/matrix-unit basis. This fact can be leveraged to replace explicit matrix multiplications with lower-complexity vector outer products, or in some cases just a combination of appropriate array slicing operations and additions.
  • Second, the original implementation was written for the most general case (i.e. it doesn't assume we're working in the pauli basis), and as such computed additional terms that aren't required when working in the pauli basis (and we can assume hermiticity and idempotency). As such I've added a new specialized construction function for the pauli basis case which is used whenever pyGSTi detects we're building error generators using the pauli basis. Side note: I'm pretty sure that we should be able to get rid of even any explicit outer products in the pauli basis case by using the fact that paulis matrices are generalized permutation matrices to set the non-zero elements of the elementary error generators directly using just indexing operations, but I haven't implemented this.

Miscellaneous

• There is a new function in lindbladtools called for random_error_generator_rates for generating random (CP-constrained) error generators. This has been a frequently requested utility and includes a number of requested features from end-users:
  • Control over which sectors are included.
  • Per-sector control over the maximum error generator rate.
  • Control over random sampling parameters for the H and SCA blocks.
  • The ability to specify a target error metric value for the random error generator. (Currently supported options are generator infidelity and total generator error). When using this option users can also set the relative contributions of H and S errors to the error metric value in the final randomly sampled error generator.
  • You can set specific error generator rates to be fixed values in addition to the randomly sampled ones.
  • You can specify a particular subset of state space labels upon which the error generators must be supported.

Unit tests and Tutorial

Finally, the new error generator propagation module is backed up by a relatively complete suite of unit tests. I've also written brand new test modules for the new/existing error generator basis code and the error generator label code which we previously didn't have any coverage for.

In addition there is a new tutorial covering the main features and usage of this module located in the tutorials sub-folder on algorithms.

Final comment

When you're reviewing this you'll notice some commented out code in a module related to low-frequency Hamiltonian models. When I originally started writing this I thought that these two modules would be inextricably linked and worked on their development in parallel. But that turned out not to be true and rewriting the git history to extricate these changes from the rest would be very hard as far as I can tell. Expect a separate PR related to this LFH code in the next month or two (I'm willing to commit to cleaning this up in advance of the next major release), but for now I'd ask you to pretend it isn't there.

Special thanks go out to @ashenmill for getting the ball rolling and for consultation throughout the development of this module. Thanks also go out to @tjproct, @jordanh6 and @kevincyoung for their help working out theory, and to @sserita and @enielse for their help with debugging and design advice.

[sserita]

Thanks as always for the very thorough work @coreyostrove! I have a few questions/comments to discuss, but nothing super critical. Having this with docstrings/unit tests/Jupyter notebook working for me locally + your PR explanation goes a long way in making this go smoothly :)

[sserita]

Just noting another unguarded STIM.

[coreyostrove]

Quick note. Took a look to see why the unit tests failed and this looks like an innocuous (but interestingly platform dependent) false-postive due to some comparison logic:

AssertionError:
E Not equal to tolerance rtol=1e-05, atol=0
E
E Mismatched elements: 1 / 1 (100%)
E Max absolute difference: 2.01233315e-17
E Max relative difference: 1.
E x: array(0)
E y: array(2.012333e-17)

TLDR: The test didn't like that 2e-17 wasn't exactly zero. That level of difference isn't actually meaningful though so I'll rework this test to use different comparison logic that doesn't trigger for something like this.

[coreyostrove]

Quick note. Took a look to see why the unit tests failed and this looks like an innocuous (but interestingly platform dependent) false-postive due to some comparison logic:

AssertionError:
E Not equal to tolerance rtol=1e-05, atol=0
E
E Mismatched elements: 1 / 1 (100%)
E Max absolute difference: 2.01233315e-17
E Max relative difference: 1.
E x: array(0)
E y: array(2.012333e-17)

TLDR: The test didn't like that 2e-17 wasn't exactly zero. That level of difference isn't actually meaningful though so I'll rework this test to use different comparison logic that doesn't trigger for something like this.

Addendum: I was kind of confused by the fact that this triggered at all. The numpy allclose function uses the formula:

absolute(a - b) <= (atol + rtol * absolute(b))

Where the default value of atol is 1e-8. But it turns out that the function np.testing.assert_allclose which uses the same exact formula for the condition instead uses a default atol of 0. So assert np.allclose(something) != np.testing.assert_allclose(something). That is a pretty silly design choice...

[sserita]

After discussion in dev meeting, we will keep this in its own submodule for now and will revisit that as it integrates in/stuff gets built on top of it.