Leakage-aware GST

[rileyjmurray]

This PR significantly expands pyGSTi's leakage modeling capabilities.

Note. This PR description was written after Corey's approval and a few follow-up commits from me.

Important stuff in new files

jupyter_notebooks/Examples/Leakage_automagic.ipynb. This illustrates how data from a 1-qubit GST experiment can be used to fit a 1-qutrit model, and how to generate reports that include leakage analysis of that 1-qutrit model.

pygsti/tools/leakage.py. I've divided this file into four regions.

• Report generation. Right now there's one function here,construct_leakage_report, which handles gauge optimization and report generation. It returns a Report object and an updated ModelEstimateResults object.
• Metrics. We have several new gate metrics of use when we don't care about how gates differ on "leakage space." These functions are called with n_leak=1 when a user calls construct_leakage_report. Users are also free to call these functions themselves.
• Model construction. Right now this consists of one function, leaky_qubit_model_from_pspec, which creates a three-level ExplicitOpModel by lifting a qubit processor spec.
• Gauge optimization. Of the three functions here, only add_lago_models is particularly user-facing. It adds additional "leakage-aware gauge optimized" or "LAGO" models to the provided ModelEstimateResults object. It's called automatically in construct_leakage_report.

pygsti/tools/sdptools.py. This has pyGSTi's semidefinite programming modeling code. Some of its contents were pulled from pygsti/tools/optools.py. Much of its contents are new. The diamond_distance_projection_model function lets us compute a superoperator's diamond-distance projection onto certain convex sets of superoperators (namely, any intersection of the sets of CPTP operators, leakage-free operators, and seepage-free operators); this same function lets us compute the projection with respect to the standard diamond norm and a "subspace" diamond norm.

Overview of API changes in existing files

Things with higher user-exposure

1. Calling b = BuiltinBasis('l2p1', 9) gets you a special leakage basis for nine-dimensional Hilbert-Schmidt space.
   The basis reflects a canonical separation of its underlying Hilbert space into a direct sum of a 2-dimensional space and a 1-dimensional space. In the future we extend this to l{x}p{y} for integers x >= 2 and y >= 1.

2. pygsti/models/explicitmodel.py has a new free function:

    transform_composed_model(mdl : ExplicitOpModel, s: GaugeGroupElement) -> ExplicitOpModel

   which uses ComposedOp, ComposedState, and ComposedPOVM to build a gauge-transformed copy of mdl that has the same parameterization as mdl.

3. The jamiolkowski_iso and inv_jamiolkowski_iso functions in pygsti/tools/jamiolkowski.py can now operate on CVXPY Expression objects. Separately, all Jamiolkowski functions in this file new accept an optional normalized argument that defaults to True, indicating if we're working with the normalized (inverse) Choi matrix.

4. The LazyBasis class (which is a parent of BuiltinBasis) has a new member called elindlookup. This is a lazily-constructed dict where

   element = basis.elements[basis.elindlookup[element_label]]

   is the basis element corresponding to element_label.

5. Our LinearOperator class (a subclass of ModelMember) now has a .shape parameter, in order to provide API similarity with common Python array types.

Things with low user-exposure that could be higher

I added several new functions to pygsti/report/reportables.py. Nearly all of these functions also has a class representation generated by the opfn_factory function from pygsti/report/modelfunction.py. The descriptions of these functions can be found in the docstring and implementation of info_of_opfn_by_name. Functions that appear in pygsti/tools/leakage.py have their own documentation. I'm leaving it to future work to include all leakage-related metrics in pygsti/tools/leakage.py.

Things with low exposure that should stay that way

The ConfidenceRegionFactoryView class has a new instance method called compute_grad_f. This is for computing the gradient of a ModelFunction object at a given point. Evaluating compute_grad_f constitutes most of the work in a call to ConfidenceRegionFactoryView.compute_confidence_interval.

pygsti/tools/matrixtools.py has a new class called IdentityOperator. Objects of this class implement @, .T, and .conj() as though they were numpy ndarrays. However, these operations are "no-ops," and so cost basically nothing. This class is used in new gauge optimization code. It comes with a free function called is_operatorlike for checking if something can plausibly be used as a linear operator with a numpy-like API.

The gaugeopt_to_target function in pygsti/algorithms/gaugeopt.py accepts an n_leak argument that's a nonnegative integer and that defaults to 0. Choices of n_leak > 0 cannot be used with the default optimizer ("ls"); one should request "L-BFGS-B" instead.

The instance method ExplicitOpModelCalc.frobeniusdist now has a different behavior for the transform_mx argument. You can pass in None or an ndarray (as in previous code) or you can pass in a pair of objects where matrixtools.is_operatorlike(obj) returns True.

Other notes

We now test for correctness of gauge optimization code in test/unit/algorithms/test_gaugeopt_correctness.py. Old tests (which were just smoke tests) were moved to test/unit/algorithms/test_gaugeopt_noerrors.py.

The code in the "Leakage-automagic" example notebook is run as a smoke test in test/unit/tools/test_leakage_gst_pipeline.py.

Some tests were added to test/unit/tools/test_optools.py.

Unresolved comments:

• https://github.com/sandialabs/pyGSTi/pull/410/files#r2255468471
• https://github.com/sandialabs/pyGSTi/pull/410/files#r2255477604 (spelling)
• https://github.com/sandialabs/pyGSTi/pull/410/files#r2255482265 (spelling)
  *https://github.com/sandialabs/pyGSTi/pull/410/files#r2255546335 (request for test to trigger a new error handling codepath).

Future work:

• systems with > 1 qubit.
• systems with > 1 leakage level.
• CPTP kite models

[rileyjmurray]

EDIT: we now have an open issue about handling this sort of thing #633.

---

Text from an email Piper sent (not showing plots since they contained real data):

I was looking at the model violation tab, trying to work out why the smaller L models weren’t fitting the data as well as might be expected as quantified by Nsigma. I think this maybe be a function of how we calculate Nsigma = (2delta logl -k)/sqrt(2k). The number of non-gauge model parameters N_p is reported as 351 in the leakage-enhanced GST report (I don’t have any particular sense of whether this is correct). The number of data points N_s is 193 for L=1. This means that k (the degrees of freedom) is negative if you just calculate it as N_s-N_p. In cases like this, pyGSTi chooses k=1. k remains small until L is fairly large, which means two delta logl must be quite small for Nsigma to be small. This explains why we don’t see any red boxes in the per sequence detail for the L=1 iteration (and Nsigma is large enough that I would expect to see some or at least more dark gray squares than can be seen).

If you directly compare two delta logl between the two fits, you can see that the two delta logl is in fact smaller in the leakage-enhanced report than in the standard report.

I think the solution here is just to make sure that the first iteration of leakage-aware GST uses L=2 instead of L=1. (Or, more generally, we use the smallest L where N_s > N_p.)

[rileyjmurray]

Leaving file-level comments so @coreyostrove and I can figure out what to split out from this PR.

[coreyostrove]

Text from an email Piper sent (not showing plots since they contained real data):

I was looking at the model violation tab, trying to work out why the smaller L models weren’t fitting the data as well as might be expected as quantified by Nsigma. I think this maybe be a function of how we calculate Nsigma = (2delta logl -k)/sqrt(2k). The number of non-gauge model parameters N_p is reported as 351 in the leakage-enhanced GST report (I don’t have any particular sense of whether this is correct). The number of data points N_s is 193 for L=1. This means that k (the degrees of freedom) is negative if you just calculate it as N_s-N_p. In cases like this, pyGSTi chooses k=1. k remains small until L is fairly large, which means two delta logl must be quite small for Nsigma to be small. This explains why we don’t see any red boxes in the per sequence detail for the L=1 iteration (and Nsigma is large enough that I would expect to see some or at least more dark gray squares than can be seen).
If you directly compare two delta logl between the two fits, you can see that the two delta logl is in fact smaller in the leakage-enhanced report than in the standard report.

I think the solution here is just to make sure that the first iteration of leakage-aware GST uses L=2 instead of L=1. (Or, more generally, we use the smallest L where N_s > N_p.)

This isn’t a uniquely leakage-related shortcoming, I’ve encountered this in other settings too (with really aggressive FPR schemes when not forcing inclusion of LGST circuits, for example) and I suggest we open up a separate github issue for this and give some additional thought to how we should be approaching the statistics in this case.

[coreyostrove]

Just finished my first pass through. This is a ton of effort and these contributions are sincerely appreciated. Thank you as always for the excellent work, @rileyjmurray!

Of the comments I've left, the only one I am deeply concerned about is the potential bug that I've pointed out in the cost function for gauge optimization as it pertains to instruments. I have good reason to suspect (I mean, the git history is pretty clear and this present on develop so it isn't really a suspicion) that in the event there is a bug here it is very likely you had nothing to do with introducing it, but since you're already knee deep in this code I'd ask that you take a look and if necessary clean this up on our behalf.

Rest of my comments are either minor clarifying questions or otherwise nice-to-have but non-blocking requests. I've also added my concurrence with a couple comments/suggestions from @nkoskelo.

[coreyostrove]

Please add an error message for the AssertionError indicating that these metrics are not supported when doing leakage-aware gauge optimization (i.e. n_leak>0).

[coreyostrove]

These look functionally equivalent to me. I’ll take your word for it, but for my own edification what about the original version was incorrect?

[coreyostrove]

This is a good suggestion and I think it would be an independently useful primitive for circuits to return their unique gate labels, probably with some flags to control matching behavior (e.g. just label names versus full labels with sslbls).

[coreyostrove]

It looks like most/all of these reportable functions hard code in that n_leak is one. I don't think there is anything to flag a user specifying something different in construct_standard_report, though. That is fine for now, but maybe add an assertion for this to avoid unexpected behavior? (i.e. a user being confused by why this was quietly changed under the hood).

[rileyjmurray]

Good point. I'll update the construct_standard_report docstring and raise an exception if the user passes in n_leak other than 0 or 1.

[coreyostrove]

Similar comment to the one I left on the leakage model regarding docstrings and API-level.

[rileyjmurray]

I'm going to say that everything here is in the private API for now.

[rileyjmurray]

@coreyostrove I think I've addressed your comments. Just pushed what should be the last round of changes.

[coreyostrove]

@coreyostrove I think I've addressed your comments. Just pushed what should be the last round of changes.

I think the only major outstanding concern here is the bug I alluded to with gauge optimization for instruments. Since this is a pre-existing bug I am alright with tabling this for the purposes of this PR, but we should open up a new github issue to track it. If you have time today can you open this issue? If not I can do so myself later this evening.

[coreyostrove]

Fantastic work, @rileyjmurray! Thank you for the significant effort this work represents and your contributions advancing our leakage modeling capabilities. Merge at will.

One last request as you merge things in, can you flesh out the summary in the PR? That’d be helpful as I prepare release notes (and for posterity).

[rileyjmurray]

The single failing test is caused by implementation details of /test/unit/mpi/test_mpi.py::MPITester::test_all. Merging as-is.