T1 fitter is bounded in amplitude, offset, and t1

[oliverdial]

Informations

• Qiskit Experiments version:
• commit 112f4ce
• Python version:
• 3
• Operating system:
• osx

What is the current behavior?

The t1 fit is constrained by default to give a value between [0,1]. This poses a few problems; if the actual data is not so constrained (because for example readout is poorly calibrated), the resulting fit is poor. If the bound becomes active, it violates the assumptions in linear error propagation and the error bars become invalid.

In extreme cases, if the inital guess is outside [0,1] for amplitude or offset, or outside some other range for t1, the fitter raises an exception, does not generate a plot, and ultimately will not put a fit result (even if a poor one) in the results db. (this is hard to reproduce w/ the T1 fake backend since it uses discretized experimental data).

These experiments (all) need to be able to give some sort of result even for non-ideal data, and should not require manual intervention to get data to fit if it is at all possible to fit it. (This was even included as an example in the t1 fitter).

Steps to reproduce the problem

What is the expected behavior?

For an admittedly extreme example, try

from qiskit_experiments.characterization import T1
from qiskit_experiments.composite import ParallelExperiment

# A T1 simulator
from qiskit_experiments.test.t1_backend import T1Backend

# Simulate T1 of 25 microseconds
t1 = 1
backend = T1Backend(t1=[t1*1e-6], initial_prob1=[-2])

# Time intervals to wait before measurement
import numpy as np
delays = np.linspace(2,4,51)
    
# Create an experiment for qubit 0,
# setting the unit to microseconds,
# with the specified time intervals
exp = T1(qubit=0, 
         delays=delays,
         unit="us")

# Run the experiment circuits with 1000 shots each,
# and analyze the result
exp_data = exp.run(backend=backend,
                   shots=1000)

# Print the result
res = exp_data.analysis_result(0)

Suggested solutions

The default parameters on experiments should be "correct" for common use cases.

never use a bounded optimizer on experimental data. They're simply too good at hiding defects in the data (ie, hiding that tomography data is intrinsically non-physical)

Instead use an unbounded one and set the result to bad if the MLE is outside physical bounds by more than a few error bars.

[eggerdj]

No sure if related but I'm getting some strange behavior when running T1 on actual backends. So far I tried Sydney. Here is the code (taken from the tutorial)

from qiskit_experiments.characterization import T1

# Create an experiment for qubit 0,
# setting the unit to microseconds,
# with the specified time intervals
exp = T1(qubit=0, delays=list(range(1, 500, 5)), unit="us")

# Run the experiment circuits with 1000 shots each,
# and analyze the result
exp_data = exp.run(backend=backend, shots=1024)

# Print the result
res = exp_data.analysis_result(0)
res

and the result

From the backend I would expect approx. 39us.

[coruscating]

@eggerdj that seems related to Qiskit/qiskit#6678

[eggerdj]

Indeed, I get a reasonable T1 with multiple of 16 samples.

[oliverdial]

Perhaps a more convincing view of the problem: using the nice fake T1 backend again, look for bias in the T1 estimate

from qiskit_experiments.characterization import T1
from qiskit_experiments.composite import ParallelExperiment

# A T1 simulator
from qiskit_experiments.test.t1_backend import T1Backend

# Simulate T1 of 25 microseconds
t1 = 25
backend = T1Backend(t1=[t1*1e-6], initial_prob1=[0])

# Time intervals to wait before measurement
import numpy as np
delays = np.linspace(2,50,51)
    
# Create an experiment for qubit 0,
# setting the unit to microseconds,
# with the specified time intervals
exp = T1(qubit=0, 
         delays=delays,
         unit="us")

# Run the experiment circuits with 1000 shots each,
# and analyze the result
results=[]
for i in range(100):
    exp_data = exp.run(backend=backend,
                       shots=1000)

    # Print the result
    res = exp_data.analysis_result(0)
    print(res)
    results.append(res)

import scipy.stats as stats
import matplotlib.pyplot as plt

fig = plt.figure()

stats.probplot([r['value']-t1*1e-6 for r in results], plot=plt.gca())
print(np.mean([r['value']-t1*1e-6 for r in results]))

As written, there's a 500 ns bias in the T1 estimate, and the distribution is non-normal

If I disable the bounds by removing it from curve_fit line 174:

    # Run curve fit                                                                                                                                                                                         
    try:
        # pylint: disable = unbalanced-tuple-unpacking                                                                                                                                                      
        popt, pcov = opt.curve_fit(
            fit_func, xdata, ydata, sigma=sigma, p0=param_p0, **kwargs
        )
    except Exception as ex:
        raise AnalysisError(
            "scipy.optimize.curve_fit failed with error: {}".format(str(ex))
        ) from ex

the bias is gone.

In general putting these hard wallls on the fit adds bias if the true value is within the measurement noise of the constraint.

[nkanazawa1989]

I agree this bound assumes discriminator is working correctly. Since this experiment uses the probability processor with level2 measurement, the experiment data value is always bound to [0, 1] even if discriminator is broken (we divide a result count by total shots, count never becomes negative and exceed total shots), and assuming a fit bound based on this is not quite bad assumption. However, of course user can override experiment setting with level1 (and with IQ data processor), and this hard-coded bounds can cause an issue.

Perhaps we can automatically generate fit bounds, or we can remove it. If data processor is broken, amp parameter can be negative (if 0/1 is flipped), for example. I think this is very rare event, so users are still able to retrieve the experimental result from the db and rerun analysis by overriding bounds. I don't have concrete idea what is the best approach here.

[oliverdial]

In part this is why I added the second example; the bounded fit is also biased even if the data is discriminated and in fact perfect in every way just because of the presence of shot noise. (However, I assume that a to-do item is to make these classes generic to if users want to use averaged IQ data or discriminated data, so the first need doesn't go away).

The bounds complicate the code, and make the fits biased and less reliable. What's the argument for keeping them?

[coruscating]

Addressed by #168.