-
Notifications
You must be signed in to change notification settings - Fork 56
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Replace calls with .flatten() to .ravel() when appropriate #451
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks good to me for my file! One thing I'd keep in mind that you will have to do this file (workspacesplots.ipynb) again when Corey merges in his substantial reporting changes since I doubt the auto-merge will catch this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Great work, @rileyjmurray! I have left some comments answering the questions that you had asked and pointed out a few (minor) inefficiencies. @sserita, none of my feedback is blocking, so feel free to approve this whenever you and @rileyjmurray think it is ready.
@@ -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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
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
<#451 (comment)>:
> @@ -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
<#451 (comment)>:
> @@ -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
<#451 (comment)>:
> @@ -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
<#451 (comment)>:
> @@ -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
<#451 (comment)>:
> @@ -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
<#451 (review)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ACNUZ5P42HJUDBUVMHCEMADZG6S3XAVCNFSM6AAAAABI5M42SCVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCMJRGY4DOOJXHA>
.
You are receiving this because you are subscribed to this thread.Message
ID: ***@***.***>
|
…ere originally .flatten(), which is strictly worse than ravel()), and simplify calculations in pygsti/report/reportables.py::general_decomposition
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Everything looks good to me, thanks @rileyjmurray!
This PR aims to resolve #431. It also resolves #424 because "why not."
More specific changes:
np.dot
with@
.x = []; y = []
asx, y = [], []
.matrixtools.py
.matrixtools.py::residuals(a, b)
and inlined equivalent code where needed.np.vdot(a.flatten(), b.flatten())
withnp.vdot(a, b)
.changed code like.list(array.flatten())
withlist(array.flat)
instead oflist(array.ravel())