This tutorial will demonstrate how perform tomography on models which, in addition to normal gates, contain quantum instruments. Quantum instruments are maps that act on a qubit state (density matrix) and produce a qubit state along with a classical outcome. That is, instruments are maps from $\mathcal{B}(\mathcal{H})$, the space of density matrices, to $\mathcal{B}(\mathcal{H}) \otimes K(n)$, where $K(n)$ is a classical space of $n$ elements.
In pyGSTi, instruments are represented as collections of gates, one for each classical "outcome" of the instrument. This tutorial will demonstrate how to add instruments to Model
objects, compute probabilities using such Model
s, and ultimately perform tomography on them. We'll start with a few familiar imports:
import pygsti
from pygsti.construction import std1Q_XYI as std
import numpy as np
Next, we'll add an instrument to our "standard" model - a 1-qubit model containing $I$, $X(\pi/2)$, and $Y(\pi/2)$ gates. The ideal instrument will be named "Iz"
(all instrument names must begin with "I"
), and consist of perfect projectors onto the 0 and 1 states. Instead of labelling the associated outcomes "0" and "1", which might me most logical, we'll name them "p0" and "p1" so it's easier to distinguish them from the final POVM outcomes which are labelled "0" and "1".
#Make a copy so we don't modify the original
target_model = std.target_model()
#Create and add the ideal instrument
E0 = target_model.effects['0']
E1 = target_model.effects['1']
# Alternate indexing that uses POVM label explicitly
# E0 = target_model['Mdefault']['0'] # 'Mdefault' = POVM label, '0' = effect label
# E1 = target_model['Mdefault']['1']
Gmz_plus = np.dot(E0,E0.T) #note effect vectors are stored as column vectors
Gmz_minus = np.dot(E1,E1.T)
target_model['Iz'] = pygsti.obj.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})
#For later use, record the identity POVM vector
povm_ident = target_model.effects['0'] + target_model.effects['1']
In order to generate some simulated data later on, we'll now create a noisy version of target_model
by depolarizing the state preparation, gates, and POVM, and also rotating the basis that is measured by the instrument and POVM.
mdl_noisy = target_model.depolarize(op_noise=0.01, spam_noise=0.01)
mdl_noisy.effects.depolarize(0.01) #because above call only depolarizes the state prep, not the POVMs
# add a rotation error to the POVM
Uerr = pygsti.rotation_gate_mx([0,0.02,0])
mdl_noisy.effects['0'] = np.dot(mdl_noisy.effects['0'].T,Uerr).T
mdl_noisy.effects['1'] = povm_ident - mdl_noisy.effects['0']
#Could also do this:
#E0 = np.dot(mdl_noisy['Mdefault']['0'].T,Uerr).T
#E1 = povm_ident - E0
#mdl_noisy['Mdefault'] = pygsti.obj.UnconstrainedPOVM({'0': E0, '1': E1})
# Use the same rotated effect vectors to "rotate" the instrument Iz too
E0 = mdl_noisy.effects['0']
E1 = mdl_noisy.effects['1']
Gmz_plus = np.dot(E0,E0.T)
Gmz_minus = np.dot(E1,E1.T)
mdl_noisy['Iz'] = pygsti.obj.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})
#print(mdl_noisy) #print the model
Instrument labels (e.g. "Iz"
) may be included within Circuit
objects, and Model
objects are able to compute probabilities for them just like normal (non-instrument) operation sequences. The difference is that probabilities are labeled by tuples of instrument and POVM outcomes - referred to as "outcome tuples" - one for each instrument and one for the final POVM:
dict(target_model.probs( pygsti.obj.Circuit(('Gx','Iz')) ))
{('p0', '0'): 0.5000000000000003, ('p0', '1'): 0.0, ('p1', '0'): 0.0, ('p1', '1'): 0.4999999999999999}
dict(target_model.probs( pygsti.obj.Circuit(('Iz','Gx','Iz')) ))
{('p0', 'p0', '0'): 0.5000000000000003, ('p0', 'p0', '1'): 0.0, ('p0', 'p1', '0'): 0.0, ('p0', 'p1', '1'): 0.5, ('p1', 'p0', '0'): 0.0, ('p1', 'p0', '1'): 0.0, ('p1', 'p1', '0'): 0.0, ('p1', 'p1', '1'): 0.0}
In fact, pyGSTi always labels probabilties using outcome tuples, it's just that in the non-instrument case they're always 1-tuples and by OutcomeLabelDict
magic can be treated as if they were just strings:
probs = target_model.probs( pygsti.obj.Circuit(('Gx',)) )
print("probs = ",dict(probs))
print("probs['0'] = ", probs['0']) #This works...
print("probs[('0',)] = ", probs[('0',)]) # and so does this.
probs = {('0',): 0.5000000000000002, ('1',): 0.4999999999999998} probs['0'] = 0.5000000000000002 probs[('0',)] = 0.5000000000000002
germs = std.germs
fiducials = std.fiducials
max_lengths = [1] # keep it simple & fast
lsgst_list = pygsti.construction.make_lsgst_experiment_list(
mdl_noisy,fiducials,fiducials,germs,max_lengths)
#print("LinearOperator sequences:")
#print(lsgst_list) #note that this contains LGST strings with "Iz"
#Create the DataSet
ds = pygsti.construction.generate_fake_data(mdl_noisy,lsgst_list,1000,'multinomial',seed=2018)
#Write it to a text file to demonstrate the format:
pygsti.io.write_dataset("../../tutorial_files/intermediate_meas_dataset.txt",ds)
Notice the format of intermediate_meas_dataset.txt, which includes a column for each distinct outcome tuple. Since not all experiments contain data for all outcome tuples, the "--"
is used as a placeholder. Now that the data is generated, we run LGST or LSGST just like we would for any other model:
#Run LGST
mdl_lgst = pygsti.do_lgst(ds, fiducials,fiducials, target_model)
#print(mdl_lgst)
#Gauge optimize the result to the true data-generating model (mdl_noisy),
# and compare. Mismatch is due to finite sample noise.
mdl_lgst_opt = pygsti.gaugeopt_to_target(mdl_lgst,mdl_noisy)
print(mdl_noisy.strdiff(mdl_lgst_opt))
print("Frobdiff after GOpt = ",mdl_noisy.frobeniusdist(mdl_lgst_opt))
Model Difference: Preps: rho0 = 0.0265047 POVMs: Mdefault: 0 = 0.0235047 1 = 0.0235564 Gates: Gi = 0.0595052 Gx = 0.038965 Gy = 0.0281307 Instruments: Iz: p0 = 0.0359309 p1 = 0.0377959 Frobdiff after GOpt = 0.010621914779793969
Instruments just add parameters to a Model
like gates, state preparations, and POVMs do. The total number of parameters in our model is
$4$ (prep) + $2\times 4$ (2 effects) + $5\times 16$ (3 gates and 2 instrument members) $ = 92$.
target_model.num_params()
92
#Run long sequence GST
results = pygsti.do_long_sequence_gst(ds,target_model,fiducials,fiducials,germs,max_lengths)
--- Circuit Creation --- 92 sequences created Dataset has 128 entries: 92 utilized, 0 requested sequences were missing --- LGST --- Singular values of I_tilde (truncating to first 4 of 6) = 4.242909389225678 1.3657661576181168 1.3394292109983001 1.329422077166059 0.0521487661472659 0.013767831452193883 Singular values of target I_tilde (truncating to first 4 of 6) = 4.242640687119286 1.414213562373096 1.414213562373096 1.4142135623730954 2.484037189058858e-16 1.506337939585075e-16 --- Iterative MLGST: Iter 1 of 1 92 operation sequences ---: --- Minimum Chi^2 GST --- Sum of Chi^2 = 57.6583 (92 data params - 76 model params = expected mean of 16; p-value = 1.29168e-06) Completed in 0.4s 2*Delta(log(L)) = 57.3792 Iteration 1 took 0.4s Switching to ML objective (last iteration) --- MLGST --- Maximum log(L) = 28.6444 below upper bound of -138546 2*Delta(log(L)) = 57.2887 (92 data params - 76 model params = expected mean of 16; p-value = 1.48819e-06) Completed in 0.1s 2*Delta(log(L)) = 57.2887 Final MLGST took 0.2s Iterative MLGST Total Time: 0.6s -- Adding Gauge Optimized (go0) --
/Users/enielse/research/pyGSTi/packages/pygsti/objects/estimate.py:525: UserWarning: Max-model params (92) <= model params (92)! Using k == 1.
#Compare estimated model (after gauge opt) to data-generating one
mdl_est = results.estimates['default'].models['go0']
mdl_est_opt = pygsti.gaugeopt_to_target(mdl_est,mdl_noisy)
print("Frobdiff after GOpt = ", mdl_noisy.frobeniusdist(mdl_est_opt))
Frobdiff after GOpt = 0.008577217079914283
The same analysis can be done for a trace-preserving model, whose instruments are constrained to add to a perfectly trace-preserving map. The number of parameters in the model are now:
$3$ (prep) + $1\times 4$ (effect and complement) + $3\times 12$ (3 gates) + $(2\times 16 - 3)$ (TP instrument) $ = 71$
mdl_targetTP = target_model.copy()
mdl_targetTP.set_all_parameterizations("TP")
print("POVM type = ",type(mdl_targetTP["Mdefault"])," Np=",mdl_targetTP["Mdefault"].num_params())
print("Instrument type = ",type(mdl_targetTP["Iz"])," Np=",mdl_targetTP["Iz"].num_params())
print("Number of model parameters = ", mdl_targetTP.num_params())
POVM type = <class 'pygsti.objects.povm.TPPOVM'> Np= 4 Instrument type = <class 'pygsti.objects.instrument.TPInstrument'> Np= 28 Number of model parameters = 71
resultsTP = pygsti.do_long_sequence_gst(ds,mdl_targetTP,fiducials,fiducials,germs,max_lengths)
--- Circuit Creation --- 92 sequences created Dataset has 128 entries: 92 utilized, 0 requested sequences were missing --- LGST --- Singular values of I_tilde (truncating to first 4 of 6) = 4.242909389225678 1.3657661576181168 1.3394292109983001 1.329422077166059 0.0521487661472659 0.013767831452193883 Singular values of target I_tilde (truncating to first 4 of 6) = 4.242640687119286 1.414213562373096 1.414213562373096 1.4142135623730954 2.484037189058858e-16 1.506337939585075e-16 --- Iterative MLGST: Iter 1 of 1 92 operation sequences ---: --- Minimum Chi^2 GST --- Sum of Chi^2 = 59.1464 (92 data params - 63 model params = expected mean of 29; p-value = 0.000787803) Completed in 0.3s 2*Delta(log(L)) = 58.8821 Iteration 1 took 0.4s Switching to ML objective (last iteration) --- MLGST --- Maximum log(L) = 29.4003 below upper bound of -138546 2*Delta(log(L)) = 58.8006 (92 data params - 63 model params = expected mean of 29; p-value = 0.000868812) Completed in 0.1s 2*Delta(log(L)) = 58.8006 Final MLGST took 0.1s Iterative MLGST Total Time: 0.5s -- Adding Gauge Optimized (go0) --
#Again compare estimated model (after gauge opt) to data-generating one
mdl_est = resultsTP.estimates['default'].models['go0']
mdl_est_opt = pygsti.gaugeopt_to_target(mdl_est,mdl_noisy)
print("Frobdiff after GOpt = ", mdl_noisy.frobeniusdist(mdl_est_opt))
Frobdiff after GOpt = 0.0084068958685607
Thats it! You've done tomography on a model with intermediate measurments (instruments).