Simplify expectation_estimation_shadow code

[bdg221]

Fixes #1929

This is an attempt at simplifying the previous code in expectation_estimation_shadow. These changes have included removing unnecessary variables and finding alternative ways of getting the same results.

I do feel like there is a balance to be found between simplicity and readability. I tried to use the example in the #1929 as the guide for that balance.

I first only pulled the relevant indexes/qubits from the measurement_outcomes based on the pauli_str's qubits so that the data would not need to be trimmed later.
b_lists_shadow = np.array([list(u) for u in measurement_outcomes[0]])[:, target_locs]

Next, I looked at a different way to do the main part of the code which checked if the Pauli bases for the indexes in the split were equal between the measurement_outcomes and pauli_str and if the corresponding measurement_outcomes had an EVEN number of 1s (mean is positive) or an ODD number of 1s (mean is negative).

1. Pauli bases check from this split (idxes) results in True or False
   np.all(u_lists_shadow[idxes] == target_obs, axis=1)

2. Get only the True indexes:
   np.nonzero(np.all(u_lists_shadow[idxes] == target_obs, axis=1))

3. Use the True indexes to get the measurements (0 or 1)
   b_lists_shadow[np.nonzero(np.all(u_lists_shadow[idxes] == target_obs, axis=1))]

4. even number of 1s results in 1 and odd number of 1s results in -1
   4a) Since the values are either 0 or 1, summing the values gives the count of 1s. Since the values were originally strings, they are set to a type int as well
   np.sum(b_lists_shadow[np.nonzero(np.all(u_lists_shadow[idxes] == target_obs, axis=1))].astype(int),axis=1)

4b) for setting the result of -1 or 1, an exponent for -1 is set to either 1 or 0 which is accomplished with mod 2
product = (-1) ** np.mod(np.sum(b_lists_shadow[np.nonzero(np.all(u_lists_shadow[idxes] == target_obs, axis=1))].astype(int),axis=1), 2)

I would like to point out that I am not certain of the math or principles being used in the code. For example, when the issue was originally opened the formula used was: 3.0 * number of 1s in measurement_outcome. However in the new code it is: +/- 3 ** number of 1s from pauli_str with the +/- decided even or odd number of 1s from measurement_outcome.

I did check the timing difference between the implementations and found them to be negligible.

License

• I license this contribution under the terms of the GNU GPL, version 3 and grant Unitary Fund the right to provide additional permissions as described in section 7 of the GNU GPL, version 3.

Before opening the PR, please ensure you have completed the following where appropriate.

• I added unit tests for new code.
• I used type hints in function signatures.
• I used Google-style docstrings for functions.
• I updated the documentation where relevant.
• Added myself / the copyright holder to the AUTHORS file

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (d468dd2) 98.21% compared to head (a31a730) 98.21%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2113      +/-   ##
==========================================
- Coverage   98.21%   98.21%   -0.01%     
==========================================
  Files          87       87              
  Lines        4156     4149       -7     
==========================================
- Hits         4082     4075       -7     
  Misses         74       74              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[bdg221]

@andreamari I took a shot at simplifying the code. I was not certain what you expected, but I am available if you would like to have a call on Discord to go over what you were thinking.

Also, I believe there is a bug in the code in the pauli_twirling_calibration section. If f_val is None then means.append(0.0) is called. However, "continue" needs to be added otherwise another value with be appended below. If you agree, I can open a bug issue to keep it separate.

[andreamari]

The code looks good to me and much cleaner!

About the potential bug I agree on opening an independent issue to keep things separated.

About the underlying theory, I also feel it is very hard to check if formulas are correct. However I remember that there is a reference somewhere in the code about a theory paper from which all the mathematical expressions are taken. If you are interested or if you have any doubts, comparing the content of that paper with the code is a possible way of checking if formals are correct.

[andreamari]

Maybe the np.mod() function here is redundant ? If I am wrong ignore this comment.

[bdg221]

You are 100% correct.

[andreamari]

Thanks @bdg221 !

[Min-Li]

@bdg221 Please check if this causes the wrong output in the two notebooks in the operator expectation estimation part, which indeed happened... I suggest reverting back to this pull request. I think it's because it would give the wrong median calculation.

https://mitiq.readthedocs.io/en/stable/examples/shadows_tutorial.html
https://mitiq.readthedocs.io/en/stable/examples/rshadows_tutorial.html

@andreamari I could look into the code of what's wrong happened, but considering that it does not give the right result (plots of the last parts of expectation value estimation). I think we should go back to the previous version.

[andreamari]

Thanks @Min-Li for catching this problem!
We missed it since all unit tests passed but I agree that both tutorials show that something is not working.

@bdg221 do you have any idea how to fix the problem? Otherwise we'll need to revert.

FYI: Since mitiq 0.32 has been released, the links with correct tutorials are:
https://mitiq.readthedocs.io/en/v0.31.0/examples/shadows_tutorial.html
https://mitiq.readthedocs.io/en/v0.31.0/examples/rshadows_tutorial.html

Once things are fixed, we can do a small patch release in order to get the "stable" version of the documentation updated.

[andreamari]

What is the reason of [0] at the end of np.nonzero(np.all(u_lists_shadow[idxes] == target_obs, axis=1))[0] ?
@bdg221 could this be the origin of all the problems?

UPDATE: It's not. Tried removing [0] and still got bad plots in tutorials.

[bdg221]

@andreamari @Min-Li @natestemen The issue issue is with slice of b_lists_shadow on line 346. Currently, the code uses indexes relative to the split with the full b_list_shadow. Instead, the indexes relative to the split need to go again the split of b_list_shadow.
Therefore, the solution is to slice b_lists_shadow down to the current split (b_lists_shadow[idxes]) then use the indexes from u_lists_shadow where the values match target_obs.
The final code should look like:

                 b_lists_shadow[idxes][
                    np.nonzero(
                        np.all(u_lists_shadow[idxes] == target_obs, axis=1)
                    )
                ].astype(int),
                axis=1,
            )

[Min-Li]

@andreamari Ahh, I tried to change the original code to make it compatible with the derandomize protocol https://arxiv.org/pdf/2103.07510.pdf however, find something went wrong here. I'll create a pull request to make the code compatible with the randomized version.
@bdg221 Thanks for correcting the code!

[Min-Li]

@bdg221, @andreamari sorry, Andrea, I don't know who had modified the workflow, but there are some issues and I'll outline some of the problems here:

In the first workflow, specifically in step 5, one cannot obtain an estimation of the expectation value merely by averaging over the 'classical snapshots.' There is an additional step involved in solving for the expectation value. To be precise, it involves finding an exact match for this local version.

In the second workflow, at step 6, what is reconstructed is the shadow quantum channel. However, there's no way to mitigate the quantum channel (i.e., the z-basis measurement in the quantum processing, as mentioned in step 2 of the first workflow). What is calibrated is the classical post-processing. Instead of manually calculating the inverse of the ideal quantum channel (@bdg221 This is the source of the 3^n factor mentioned in your previous question), the calibrated version replaces the inverse channel with the reconstruction results.

[purva-thakre]

I don't know who had modified the workflow

That would be me :)

@Min-Li would you like to chat about this over a call sometime? I prefer conversations over going back and forth via messages. We can work on making the corrections then.

[bdg221]

@purva-thakre @Min-Li I'd love to be on that call too.

[Min-Li]

@bdg221 @purva-thakre Great! I could meet tomorrow afternoon after 3pm CT.

[purva-thakre]

Great! I could meet tomorrow afternoon after 3pm CT.

Sorry! I have got other plans tomorrow afternoon. I am available next Friday after next week.