sandialabs / pyGSTi

A python implementation of Gate Set Tomography
http://www.pygsti.info
Apache License 2.0
132 stars 55 forks source link

Remove wildcard optimization code that called cvxopt directly #444

Closed rileyjmurray closed 3 weeks ago

rileyjmurray commented 1 month ago

This is a follow-up on https://github.com/ionq/pyGSTi/issues/1. There are currently three functions that call cvxopt's nonlinear optimization interface directly,

As far as user exposure goes, these functions could be called if a GSTBadFitOptions constructor was called with certain values for wildcard_methods ("cvxopt", "cvxopt_smoothed", and "cvxopt_small"). I looked throughout pygsti's codebase -- including tests and notebooks -- and there was not a single place were we hit this codepath.

My original plan for this PR was to rewrite these functions to avoid relying on cvxopt explicitly, and instead rely on cvxpy. But once I started messing with the code it was clear (see details below) that this would take an inordinate amount of work. Given that these functions seem to have low user exposure I suggest we just remove them entirely in pyGSTi 0.9.13.

Details

I spent time trying to figure out how optimize_wildcard_budget_cvxopt could be modified to use cvxpy instead. Right now it's relying on cvxopt's general nonlinear (convex) optimization interface. This got me nervous because it's possible to define extremely complicated problems using that interface.

Luckily, optimize_wildcard_budget_cvxopt only uses cvxopt's advanced interface for a single nonlinear constraint! The objective function and all other constraints are linear. So the question is, can I model the nonlinear constraint using cvxpy's primitives? Here's the desired constraint at a mathematical level, where probs is a function of the optimization variable and everything else in the expression is constant:

2 * sum(total_counts * freqs * log(freqs / probs)) <= two_dlogl_threshold.

This was promising at first. CVXPY is perfectly happy to work with these constraints as long as probs is an affine function of the optimization variable.

Unfortunately, the code for evaluating probs as a function of the decision variable is very complicated. Here's a taste of the call sequence, in terms of the optimization variable x:

  1. budget.from_vector(_np.array(x)), where budget is an instance of PrimitiveOpsWildcardBudgetBase. This basically just sets budget.wildcard_vector = _np.array(x) and budget.per_op_wildcard_vector = _np.array(x).
  2. The value of probs that we care about gets written to the second argument of budget.update_probs(...). So the question question is how budget.update_probs(...) writes to its second argument as a function of self.wildcard_vector or self.per_op_wildcard_vector.
  3. budget.update_probs(...) only accesses self. twice.
    • The first time is in a function call circuit_budgets = self.circuit_budgets(...). There are two codepaths here depending on whether or not something called "precomp" has been created. The codepath when precomp is provided is very simple (a matrix-vector product with self.wildcard_vector). The codepath with precomp is None is less simple, but at the end of the day it's an affine function of self.wildcard_error_per_op (which is a property that basically exposes self.per_op_wildcard_vector as a dictionary keyed by Label objects).
    • The second time is to call a function self.precompute_for_same_probs_freqs that seems to not depend on self.wildcard_vector or self.per_op_wildcard_vector.
  4. There's a for-loop over enumerate(zip(circuits, circuit_budgets, probs_freqs_precomp)) that's about 130 lines long. We care about how the variable W is used in that for-loop. It's mentioned 19 times between the code and comments. The for loop makes calls to several other functions defined in the file -- both instance methods and free functions (that is, functions that aren't attached to a class in any way).

I decided to give up once I saw the complexity in step 4 of that call sequence.

coreyostrove commented 1 month ago

I am down with all of this. Excellent work sorting through all of this, @rileyjmurray.

Edit: I think this is fine to approve, btw, but will hold off temporarily to allow other to chime in with opinions on this.

rileyjmurray commented 1 month ago

Ping @enielse

coreyostrove commented 1 month ago

Quick follow up following today's meeting. It sounds like we agreed to add a note in the place of the original code for these functions pointing to the latest SHA and or pyGSTi version for which these functions were present. Once you've added those notes I think we're good to merge this.

enielse commented 3 weeks ago

We can just remove the wildcard optimization code using cvxopt. These cvxopt-specific routines were just used as a part of research when the QPL was exploring various options for wildcard optimization. It's possible we'd want to reference the optimization code at some later time, so we could either comment it out or replace the routines with a comment referencing the last commit where the code exists to leave this door open.