New QEM benchmarking method: "rotated" RB circuits

[Misty-W]

Description

Fixes #2005.

Add a function generate_rotated_rb_circuits which generates a list of randomized benchmarking circuits, each with a random $R_z$ rotation inserted in the middle.

This benchmarking method enables testing QEM techniques in more general scenarios, closer to real-world applications in which expectation values can take arbitrary values.

Rotated randomized bencmarking circuits are randomized benchmarking circuits in which an $R_z(\theta)$ rotation is inserted in the middle, such that:

$$ C(\theta) = G_n \dots G_{n/2 +1} R_z(\theta)G_{n/2} \dots G_2 G_1 $$

where $G_j$ are Clifford elements or Clifford gates.

The circuits generate expectation values which are sinusoidal functions of $\theta$, which in the ideal (noiseless) case vary in a continuous interval of $E_{\rm ideal} \in [-1, 1]$.

Since (up to factors of 2) we have $R_z(\theta) =cos(\theta) I + i \ sin(\theta) Z$, the rotated Clifford circuit $C(\theta)$ can be written as a linear combination of just two Clifford circuits and, therefore, it is still easy to classically simulate.

---

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 lines are covered by tests ✅

Comparison is base (b6cb5c1) 98.28% compared to head (b6639e1) 98.32%.
Report is 1 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2028      +/-   ##
==========================================
+ Coverage   98.28%   98.32%   +0.03%     
==========================================
  Files          86       87       +1     
  Lines        4084     4118      +34     
==========================================
+ Hits         4014     4049      +35     
+ Misses         70       69       -1     

Files	Coverage Δ
mitiq/benchmarks/__init__.py	100.00% <100.00%> (ø)
mitiq/benchmarks/randomized_benchmarking.py	97.29% <100.00%> (+2.29%)	⬆️
...itiq/benchmarks/rotated_randomized_benchmarking.py	100.00% <100.00%> (ø)

... and 1 file with indirect coverage changes

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

[natestemen]

Nice, quick, work! Looking forward to using this to do some more testing! I've got a few higher level comments before diving into the details.

Math

The math from the docstring isn't being displayed properly. You can either use double backslashes in place of the existings \s, or you can use the following syntax.

:math:`equation`

for inline math, and

..math::
    equation

for display math. An example of these in use can be found here:

mitiq/mitiq/pec/channels.py

Lines 94 to 104 in b6cb5c1

	r"""Maps a set of Kraus operators into a single superoperator
	matrix acting by matrix multiplication on vectorized
	density matrices.
	The returned matrix :math:`S` is obtained with the formula:
	.. math::
	S = \sum_j K_j \otimes K_j^*,
	where :math:`\{K_j\}` are the Kraus operators.
	The mapping is based on the following isomorphism:

API

The other thing that I could imagine we might want to use this function for is testing EM methods as a function of this $\theta$ parameter. Instead of using a seed to set the angle, what if we let the user pass a value for \theta, falling back on random value when it's not passed. That would allow us to write code such as

for theta in np.linspace(0, 2*np.pi, 50):
    circuit = generate_rotated_rb_circuits(...) # passing theta in here
    # calculate circuits noisy, ideal, and mitigated expval
    # plot the error as a function of theta

Does this seem like something that would be interesting?

Combining this with generate_rb_circuits

Did you consider merging this function with the existing function? I could imagine we could allow that function to take a parameter theta that would rotate the circuit midway through if it was passed in. Otherwise, returning the "standard" RB circuit. We probably want to avoid overloading functions with too many options, but this may be okay? Not suggesting we have to go this route, just curious if you considered it.

---

Benchmarking utilities consistency

This PR made me realize that our utilities to generate benchmarking circuits all return a single QPROGRAM except for RB circuits (and now rotated RB circuits) where we return List[QPROGRAM]. @andreamari is there a reason for this? I feel like we should try to be consistent here. This is not necessarily something we should address in this PR, but maybe it's something we want to open an issue for if we want to address it?

[Misty-W]

Thanks @natestemen! All good points to consider.

I've fixed the math rendering. That's on me for looking too quickly at the math italics and centering in the API doc.

As for combining the functions, I thought about it, but decided against it because when we insert the rotation, it's no longer really an RB circuit anymore, but rather a linear combination of Clifford groups whose ideal expectation values can vary in the interval $[-1, 1]$.

As for setting $\theta$, by default it should be randomly sampled for benchmarking purposes and generally not specified by the user. The seed is for troubleshooting purposes or other cases where we want to generate the same circuit repeatedly (e.g. in a tutorial where we want consistent results).

Changed my mind on the last point above because $\theta$ should be consistent across the trials. Updating code accordingly (good catch!)

[Misty-W]

As for List[QPROGRAM] vs. QPROGRAM, it's because a GHZ, W state, or mirror circuit for a given number of qubits is a unique circuit, but an RB circuit for a given number of qubits and Clifford depth is a random combination of Clifford groups. Typically multiple RB circuits are generated and executed, and then the results are averaged.

[andreamari]

Thanks @Misty-W for the clean and simple solution!
I have only one doubt: are you sure about adding a layer many rotations (one for each qubit)?

Initially I was thinking of a single rotation gate applied to just one qubit, possibly a qubit approximately in the middle of the quantum register. If I am not wrong, this should also work and should produce (with very high probability) an expectation value which is exactly given by $E(\theta)=cos(\theta)$.

With a layer of many rotations, if I am not wrong, the expectation value will contain a product of many sinusoidal functions. So $E(\theta)$ will be more complicated.

Mine are just untested guesses. So, a practical strategy is probably to just plot $E(\theta)$ and check which of the two solutions looks better.

[nathanshammah]

Great work @Misty-W! @Misty-W @andreamari what do you think of adding this functionality to the main RB function, e.g., as a theta optional parameter? Is it too different of a functionality from the standard RB tests?

[Misty-W]

thanks for the suggestions @andreamari !

I've changed the function to only insert an $R_z$ gate on the first qubit, in (or near) the middle of the circuit.

Here are some plots of the output- seems to be what we want in terms of the rotation, but any idea why they are so "spiky"?

thetas = np.linspace(0, 2 * np.pi, 100)
zero_probs = []
for theta in thetas:
    circuit = generate_rotated_rb_circuits(
            n_qubits=2, num_cliffords=5, theta=theta,
          )[0]
    wvf = circuit.final_state_vector()
    zero_probs.append(wvf[0] ** 2)

plt.plot(thetas, zero_probs)

or alternatively,

thetas = np.linspace(0, 2 * np.pi, 100)
results = []
for theta in thetas:
    circuit = generate_rotated_rb_circuits(n_qubits=2, num_cliffords=5, theta=theta)[0]
    results.append(cirq.DensityMatrixSimulator().simulate(circuit).final_density_matrix[0, 0].real)

plt.plot(thetas, results)

[andreamari]

Here are some plots of the output- seems to be what we want in terms of the rotation, but any idea why they are so "spiky"?

About different plots with different simulation methods:
In the first plot you computed the square of the wave function amplitude, instead of the square of the absolute value. This is probably the reason why the two plots look very different.

About the spikes:
I think they are the consequence of two facts:

• For each theta a completely new random RB circuit is generated. This explains the non-continuous jumps. Adding the possibility of setting the same seed should help.

• Since the circuit is small (2 qubits), there is a non-negligible probability that the rotation has no effect on the expectation value. This only happens in the particular cases in which
  $$C_1 C_2... C_{d/2} Z C_{d/2 + 1} ... C_d = P$$ where P is a Pauli operator diagonal in the Z basis, in those cases the expectation value is always 1 for any theta. This is very unlikely when the number of qubits is large, but it can happen quite often for 2 qubits.

At the moment it's not easy to increase the number of qubits of RB circuits (this could be a new interesting issue), so my suggestions to avoid strange effects are:

• Add the possibility of setting a seed for the RB circuit. This is always a good thing to have.
• [Alternative A] If theta!=0, discard a pathological RB circuit if the expectation value is equal to 1 and regenerate a new random RB circuit.
• [Alternative B] In addition to the circuit also return the ideal expectation value. In this way, the user can easily detect and discard pathological cases.

You may have better solutions. I think anything is fine as long as the user is aware of the expected behavior and doesn't get surprised by unexpected results (as we did! 😃 ).

[Misty-W]

@andreamari, since the Cirq random Clifford functions don't accept a random seed, I went with [Alternative A]:

If theta!=0, discard a pathological RB circuit if the expectation value is equal to 1 and regenerate a new random RB circuit.

[Misty-W]

RTD problem is unrelated to this PR- seems GitHub automatically shows a failure when RTD takes longer than 1 hour to finish building, even though it built successfully.

[andreamari]

Hi Misty,

could you test if the new function produces a smooth Cos(theta) plot?
As you tried in this message.

I see some problems with the current code:

• Main problem: it seems that we replace a pathological circuit, i.e. the case with E(theta)=1, with a new random RB circuit but without the intermediate RZ gate. We should instead add the RZ gate otherwise E(theta) would still be = to 1, so still pathological.
• Second problem: the current code doesn't check if the second attempt is pathological as well.

Maybe a good pseudo-workflow is the following:

1. Generate a RB circuit (rb_circuit) and save a copy of it (rb_copy).
2. Insert a Z gate (equivalent to RZ(\pi) ) in the middle of the circuit.
3. If the expectation value is close to 1, we have a pathological circuit, and we restart from step 1. Else, the case is not pathological, and we can go to next step.
4. Insert the R_Z(theta) rotation in rb_copy and return it.

Whatever solution you decide to take, the final goal is having a smooth cos(theta) plot for the expectation value.

[Misty-W]

hi @andreamari,

could you test if the new function produces a smooth Cos(theta) plot? As you tried in this message.

Yes, I'll regenerate the plots with the new code.

I see some problems with the current code:

• Main problem: it seems that we replace a pathological circuit, i.e. the case with E(theta)=1, with a new random RB circuit but without the intermediate RZ gate. We should instead add the RZ gate otherwise E(theta) would still be = to 1, so still pathological.
• Second problem: the current code doesn't check if the second attempt is pathological as well.

I'm using generate_rotated_rb_circuits recursively, so the intermediate RZ gate is inserted and the second attempt is also checked. If I had used generate_rb_circuits the above problems would have materialized.

[andreamari]

I'm using generate_rotated_rb_circuits recursively, so the intermediate RZ gate is inserted and the second attempt is also checked. If I had used generate_rb_circuits the above problems would have materialized.

Ah, I missed the "rotated" word. Now I see!
However, the second recursive call to generate_rotated_rb_circuits is missing the argument theta.

Note also that you may enter an infinite recursion when theta is a multiple of 2pi.

[Misty-W]

In the end, I decided to use a seed. It lets the user generate a set of rotated RB circuits with different values of theta but otherwise the same circuit, which is useful for testing and troubleshooting. Once the user's workflow is set up, then they can randomly choose or sweep over theta with or without the seed, and e.g. for comparison plot the error rate instead of the raw expectation value.

To do this efficiently I used some code from the Cirq functions _random_single_qubit_clifford and _random_two_qubit_clifford. I have credited the developers in a similar manner to what is done in the shadows module.

[Misty-W]

Plotting the expectation values as a function of theta at three seed values:

import cirq
import numpy as np
import matplotlib.pyplot as plt

from mitiq.benchmarks.rotated_randomized_benchmarking import (
    generate_rotated_rb_circuits,
)

thetas = np.linspace(0, 2 * np.pi, 100)
seeds = (1, 2, 3)
results = []
for seed in seeds:
  zero_probs = []
  for theta in thetas:
      circuit = generate_rotated_rb_circuits(
              n_qubits=2, num_cliffords=5, theta=theta, seed=seed,
            )[0]
      zero_probs.append(cirq.DensityMatrixSimulator()
          .simulate(circuit)
          .final_density_matrix[0, 0]
          .real
      )
  results.append(zero_probs)
plt.plot(thetas, results[0], label="seed = 1")
plt.plot(thetas, results[1], label="seed = 2")
plt.plot(thetas, results[2], label="seed = 3")
plt.legend()
plt.xlabel(r"$\theta$")
plt.ylabel("E")
plt.show()

[andreamari]

Cool! Very nice to have the seed option.
This LGTM, once the comments below are addressed.

[Misty-W]

thanks for the helpful comments @andreamari, @natestemen, and @nathanshammah!

I think I've addressed everything, so once the tests pass here I'll merge this.