Replace calls with .flatten() to .ravel() when appropriate

[rileyjmurray]

This PR aims to resolve #431. It also resolves #424 because "why not."

More specific changes:

• removed a giant commented-out code block in algorithms/randomcircuit.py.
• replaced some calls to np.dot with @.
• simplified some code branching in a handful of cases.
• reduced use of multiple statements per line. Example: rewrite x = []; y = [] as x, y = [], [].
• I've left a note asking if it's cool to remove commented-out code in a function function in matrixtools.py.
• removed matrixtools.py::residuals(a, b) and inlined equivalent code where needed.
• replaced instances of np.vdot(a.flatten(), b.flatten()) with np.vdot(a, b).
• changed code like list(array.flatten()) with list(array.flat) instead of list(array.ravel()).
  • Apparently ravel is often faster than flat, so I'm actually sticking with ravel.

[robinbk]

Quick comment re: Corey's observation about Hamiltonian error rates. The part about axis and angle of rotations is only true for d=2, ie a single qubit. A Hamilton defines a unitary rotation, represented as a superoperator that is an orthogonal rotation in SO(d^2-1). For d=2, this is SO(3) which is basically SU(2), and so every rotation corresponds to an axis and an angle.

For d>2, we're in SO(>>3), and the rules are different. A rotation can be around multiple axes (the entire algebra of operators that commute with the rotation), and has multiple angles. So caveat emptor!!

I don't know if this has any effect on Corey's actual conjecture -- it may still be true -- but if the intuition about axes and angles is critical, there's a risk of GIGO.

-R

On Tue, Jun 11, 2024, 5:54 PM coreyostrove @.***> wrote:

@.**** requested changes on this pull request.

Great work, @rileyjmurray https://github.com/rileyjmurray! I have left some comments answering the questions that you had asked and pointed out a few (minor) inefficiencies. @sserita https://github.com/sserita, none of my feedback is blocking, so feel free to approve this whenever you and @rileyjmurray https://github.com/rileyjmurray think it is ready.

In pygsti/modelmembers/operations/eigpdenseop.py https://github.com/sandialabs/pyGSTi/pull/451#discussion_r1635639512:

@@ -432,13 +432,13 @@ def deriv_wrt_params(self, wrt_filter=None): dMx = _np.zeros((self.dim, self.dim), 'complex') for prefactor, (i, j) in pdesc: dMx[i, j] = prefactor

• tmp = _np.dot(self.B, _np.dot(dMx, self.Bi))
• tmp = self.B @ (dMx @ self.Bi)

Note to future selves (no immediate action item): I am pretty sure that dMx should either be diagonal for non-degenerate matrices, in which case we can do the first product in O(n^2) time, or else block-diagonal (degenerate case) in which case we should still be able to get a speedup by leveraging that structure.

In pygsti/report/reportables.py https://github.com/sandialabs/pyGSTi/pull/451#discussion_r1635642837:

@@ -2021,7 +2020,7 @@ def error_generator_jacobian(opstr): for i, gl in enumerate(opLabels): for k, errOnGate in enumerate(error_superops): noise = first_order_noise(opstr, errOnGate, gl)

• jac[:, i * nSuperOps + k] = [_np.vdot(errOut.flatten(), noise.flatten()) for errOut in error_superops]
• jac[:, i * nSuperOps + k] = [_np.vdot(errOut.ravel(), noise.ravel()) for errOut in error_superops]

I think these ravels are superfluous as vdot already flattens its inputs.

In pygsti/report/reportables.py https://github.com/sandialabs/pyGSTi/pull/451#discussion_r1635655648:

@@ -2161,10 +2160,10 @@ def general_decomposition(model_a, model_b): if gl == gl_other or abs(rotnAngle) < 1e-4 or abs(rotnAngle_other) < 1e-4: decomp[str(gl) + "," + str(gl_other) + " axis angle"] = 10000.0 # sentinel for irrelevant angle

• real_dot = _np.clip(
• _np.real(_np.dot(decomp[str(gl) + ' axis'].flatten(),
• decomp[str(gl_other) + ' axis'].flatten())),
• -1.0, 1.0)
• real_dot = _np.real(_np.dot(
• decomp[str(gl) + ' axis'].ravel(), decomp[str(gl_other) + ' axis'].ravel()
• )) # Riley question: should this be vdot instead of dot?
• real_dot = _np.clip(real_dot, -1.0, 1.0)

The short answer is no.

The longer answer: I traced through this function, and what is being calculated here is the inner product between the vectors of Hamiltonian rates for two different operators. This is of interest because the hamiltonian error generator rates for a psuedovector whose direction defines an axis of rotation, and whose length gives the angle. Assuming I understand everything here correctly the following should be true.

1. logG is always real-valued (this is enforced on line 2155 and within tools.matrixtools.approximate_matrix_log.
2. The hamiltonian projections are hardcoded to be with respect to the pauli-product basis. I believe that combining these two facts it should be the case that the hamiltonian projection vectors are always going to be real-valued. Ergo, there is no need to take a conjugate transpose. (as a corollary to that I also think the np.real call is unnecessary).

Fwiw I also think the ravel itself is uneeded. These vectors are references to a LindbladCoefficientBlock's block_data attribute, and for a Hamiltonian coefficient block this should already by a 1D array.

All that said, I might be missing something here. If I were you I'd stick in a couple assert statements and run the unit test suites (which should include report generation which will hit this codepath) and see if the stated assumptions are violated in some context.

In pygsti/report/workspaceplots.py https://github.com/sandialabs/pyGSTi/pull/451#discussion_r1635667859:

@@ -3033,13 +3034,13 @@ def _create(self, dataset, target, maxlen, fixed_lists, scale):

     xs = list(range(svals.size))
     trace1 = go.Bar(

• x=xs, y=list(svals.flatten()),
• x=xs, y=list(svals.flat),

the flat property for ndarrays returns an instance of a flatiter object, so I had a suspicion this might be a bit slower. TLDR: For small arrays (4x4) flat is about 40% faster. For larger arrays (16x16 and 64x64) flat is about 15% slower. May well not matter here (depends on the typical size of svals, which I don't know off hand), but worth keeping in mind.

In pygsti/tools/matrixtools.py https://github.com/sandialabs/pyGSTi/pull/451#discussion_r1635671089:

@@ -682,11 +682,16 @@ def approximate_matrix_log(m, target_logm, target_weight=10.0, tol=1e-6): assert(_np.linalg.norm(m.imag) < 1e-8), "Argument m must be a real matrix!" mx_shape = m.shape

• #

• Riley note: I'd like to remove all commented-out code in this function.

• @Corey or @Stefan -- you okay with that?

• #

If you're referring to that block on the bottom beginning with 'OLD DEBUG' then go for it.

— Reply to this email directly, view it on GitHub https://github.com/sandialabs/pyGSTi/pull/451#pullrequestreview-2111687978, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNUZ5P42HJUDBUVMHCEMADZG6S3XAVCNFSM6AAAAABI5M42SCVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCMJRGY4DOOJXHA . You are receiving this because you are subscribed to this thread.Message ID: @.***>