QInfer: Statistical inference software for quantum applications

QInfer: Statistical inference soft w a re fo r quantum applications Christopher Granade 1 2 , Christopher F errie 3 , Ian Hincks 4 5 , Steven Casagrande , Thomas Alexander 5 6 , Jonathan Gross 7 , Michal K ononenk o 5 6 , and Y uval Sanders 5 6 8 1 School of Physics, University of Sydney , Sydney , NSW, Australia 2 Centre for Engineered Quantum Systems, Universit y of Sydney , Sydney , NSW, Australia 3 University of T echnology Sydney , Centre for Quantum Soft w are and Info rmation, Ultimo NSW 2007, Australia 4 Department of Applied Mathematics, Universit y of Waterloo, Waterloo, ON, Canada 5 Institute for Quantum Computing, Universit y of Waterloo, Waterloo, ON, Canada 6 Department of Physics and Astronomy , University of W aterloo, Waterloo, ON, Canada 7 Center for Quantum Info rmation and Control, University of New Mexico, Albuquerque, NM 87131-0001, USA 8 Department of Physics and Astronomy , Macquarie Universit y , Sydney , NSW, A ustralia April 14, 2017 Characterizing quan tum systems through ex- p erimen tal data is critical to applications as div erse as metrology and quan tum computing. Analyzing this exp erimen tal data in a robust and repro ducible manner is made challenging, ho w ev er, b y the lack of readily-av ailable soft- w are for p erforming principled statistical anal- ysis. W e impro v e the robustness and repro- ducibilit y of c haracterization b y in tro ducing an op en-source library , QInfer, to address this need. Our library mak es it easy to analyze data from tomography , randomized b enc hmark- ing, and Hamiltonian learning exp erimen ts ei- ther in p ost-processing, or online as data is ac- quired. QInfer also pro vides functionalit y for predicting the p erformance of prop osed exp eri- men tal proto cols from simulated runs. By deliv- ering easy-to-use c haracterization tools based on principled statistical analysis, QInfer helps ad- dress many outstanding c hallenges facing quan- tum tec hnology . Contents 1 Introduction 1 2 Inference and Particle Filtering 2 Christopher Granade: cgranade@cgranade.com , www.cgranade.com Complete source co de available at DOI 10.5281/zeno do.157007 . 3 Applications in Quantum Information 4 3.1 Phase and Frequency Learning . . . . . 4 3.2 State and Process T omography . . . . . 5 3.3 Randomized Benchmarking . . . . . . . 7 4 Additional Functionality 7 4.1 Region Estimation and Error Bars . . . . 7 4.2 Derived Models . . . . . . . . . . . . . . 8 4.3 T ime-Dependent Models . . . . . . . . . 9 4.4 Performance and Robustness T esting . . 10 4.5 Parallelization . . . . . . . . . . . . . . . 11 4.6 Other Features . . . . . . . . . . . . . . . 13 5 Conclusions 14 Acknowledgments 14 References 14 A Custom Model Example 17 B Resampling 17 1 Intro duction Statistical modeling and parameter estimation play a critical role in many quantum applications. In quan- tum information in particular , the pursuit of large- scale quantum information processing devices has mo- tivated a range of different characterization protocols, and in turn, new statistical models. For example, quan- tum state and process tomography are widely used to characterize quantum systems, and are in essence 1 matrix-valued parameter estimation pr oblems [ 1 , 2 ]. Similarly , randomized benchmarking is now a main- stay in assessing quantum devices, motivating the use of rigorous statistical analysis [ 3 ] and algorithms [ 4 ]. Quantum metrology , meanwhile, is intimately con- cerned with what parameters can be extracted from measurements of a physical system, immediately ne- cessitating a statistical view [ 5 , 6 ]. The prevalence of statistical modeling in quantum applications should not be surprising: quantum me- chanics is an inherently statistical theory , thus infer- ence is an integral part of both experimental and theo- retical practice. In the former , experimentalists need to model their systems and infer the value of parameters for the purpose of improving control as well as validat- ing performance. In the latter , numerical experiments utilizing simulated data are now commonplace in the- oretical studies, such that the same infer ence pr oblems are encountered usually as a necessary step to answer questions about optimal data processing protocols or experiment design. In both cases, we lack tools to rapidly prototype and access inference strategies; QIn- fer addresses this need by providing a modular inter- face to a Monte Carlo algorithm for performing statis- tical inference. Critically , in doing so, QInfer also supports and en- ables open and repr oducible resear ch practices. Par- allel to the challenges faced in many other disciplines [ 7 ], physics resear ch cannot long survive its own cur- rent practices. Open access, open source, and open data provide an indispensable means for r esearch to be repr oducible, ensuring that resear ch work is useful to the communities invested in that resear ch [ 8 ]. In the particular context of quantum information research, open methods are especially critical given the impact of statistical errors that can undermine the claims of published resear ch [ 9 , 10 ]. Ensuring the repr oducibil- ity of resear ch is critical for evaluating the extent to which statistical and methodological errors undermine the credibility of published r esearch [ 11 ]. QInfer also constitutes an important step towards a more general framework for quantum verification and validation (QCVV). As quantum information proces- sor prototypes become more complex, the challenge of ensuring that noise processes affecting these de- vices conform to some agr eed-upon standard becomes more dif ficult. This challenge can be managed, at least in principle, by developing confidence in the truth of certain simplifying assumptions and approximations. The value of randomized benchmarking, for example, depends strongly upon the extent to which noise is approximately Pauli [ 12 ]. QInfer provides a valuable framework for the design of automated and efficient noise assessment methods that will enable the compar - ison of actual device performance to the specifications demanded by theory . T o the end of enabling reproducible and accessible resear ch, and hence providing a r eliable pr ocess for in- terpreting advances in quantum information process- ing, we base QInfer using openly-available tools such as the Python programming language, the IPython in- terpreter , and Jupyter [ 13 , 14 ]. Jupyter in particular has already proven to be an invaluable tool for repr o- ducible resear ch, in that it provides a powerful frame- work for describing and explaining resear ch software [ 15 ]. W e provide our library under an open-source li- cense along with examples [ 16 ] of how to use QInfer to support r eproducible resear ch practices. In this way , our library builds on and supports recent efforts to de- velop repr oducible methods for physics resear ch [ 17 ]. QInfer is a mature open-source softwar e library writ- ten in the Python programming language which has now been extensively tested in a wide range of infer- ential problems by various r esearch groups. Recogniz- ing its maturity through its continuing development, we now formally release version 1.0. This maturity has given its developers the opportunity to step back and focus on the accessibility of QInfer such that other re- searchers can benefit from its utility . This short paper is the culmination of that effort. A full Users’ Guide is available in the ancillary files. W e proceed as following. In Section 2 , we give a brief introduction to Bayesian inference and particle fil- tering, the numerical algorithm we use to implement Bayesian updates. In Section 3 , we describe applica- tions of QInfer to common tasks in quantum informa- tion processing. Next, we describe in Section 4 addi- tional features of QInfer before concluding in Section 5 . 2 Inference a nd Pa rticle Filtering QInfer is primarily intended to serve as a toolset for implementing Bayesian approaches to statistical infer- ence. In this section, we provide a brief review of the Bayesian formalism for statistical inference. This sec- tion is not intended to be comprehensive; our aim is rather to establish the language needed to describe the QInfer codebase. In the Bayesian paradigm, statistical inference is the process of evaluating data obtained by sampling an unknown member of a family of r elated pr obability distributions, then using these samples to assign a rel- ative plausibility to each distribution. Colloquially , we think of this family of distributions as a model param- eterized by a vector x of model parameters. W e then express the probability that a dataset D was obtained from the model parameters x as Pr ( D | x ) and read it as 2 “the probability of D given that the model specified by x is the correct model.” The function Pr ( ·| x ) is called the likelihood function, and computing it is equivalent to simulating an experiment 1 . For example, the Born rule is a likelihood function, in that it maps a known or hypothetical quantum density matrix x ≡ ρ to a distri- bution over measurement outcomes of a measurement D ∈ { E , 1 − E } via Pr ( D = E | x ) = T r ( E ρ ) . (1) The problem of estimating model parameters is as follows. Suppose an agent is provided with a dataset D and is tasked with judging the probability that the model specified by a given vector x is in fact the correct one. According to Bayes’ rule, Pr ( x | D ) = Pr ( D | x ) Pr ( x ) Pr ( D ) , (2) where Pr ( x ) is a probability distribution called the prior distribution and Pr ( x | D ) is called the posterior distribution. If the agent is provided with a prior distribution, then they can estimate parameters using Bayes’ rule. Note that Pr ( D ) can be computed through marginalization, which is to say that the value can in principle be calculated via the equation Pr ( D ) = Z x Pr ( D | x ) Pr ( x ) d x . (3) For the inference algorithm used by QInfer, Pr ( D ) is an easily computed normalization constant and there is no need to compute a possibly complicated integral. Importantly , we will demand that the agent’s data processing approach works in an iterative manner . Consider the example in which the data D is in fact a set D = { d 1 , . . . , d N } of individual observations. In most if not all classical applications, each individual datum is distributed independently of the r est of the data set, conditioned on the true state. Formally , we write that for all j and k such that j 6 = k , d j ⊥ d k | x . This may not hold in quantum models where measure- ment back-action can alter the state. In such cases, we can simply redefine what the parameters x label, such that this independence property can be taken as a con- vention, instead of as an assumption. Then, we have that Pr ( x | d 1 , . . . , d N ) = Pr ( d N | x ) Pr ( x | d 1 , . . . , d N − 1 ) Pr ( d N ) . (4) 1 Here, we use the word “simulation” in the sense of what V an den Nest [ 18 ] terms “strong simulation,” as opp osed to drawing data consistent with a given mo del (“weak simulation”). In other words, the agent can process the data sequen- tially where the prior for each successive datum is the posterior from the last. This Bayes update can be solved analytically in some important special cases, such as frequency estimation [ 19 , 20 ], but is more generally intractable. Thus, to de- velop a robust and generically useful framework for parameter estimation, the agent relies on numerical al- gorithms. In particular , QInfer is largely concerned with the particle filtering algorithm [ 21 ], also known as the sequential Monte Carlo (SMC) algorithm. In the context of quantum information, SMC was first proposed for learning from continuous measurement recor ds [ 22 ], and has since been used to learn from state tomography [ 23 ], Hamiltonian learning [ 24 ], and randomized benchmarking [ 4 ], as well as other appli- cations. The aim of particle filtering is to replace a continuous probability distribution Pr ( x ) with a discrete approxi- mation X k w k δ ( x − x k ) , (5) where w = ( w k ) is a vector of probabilities. The entry w k is called the weight of the particle , labeled k , and x k is the location of particle k . Of course, the particle filter P k w k δ ( x − x k ) does not directly approximate Pr ( x ) as a distribution; the particle filter , if considered as a distribution, is sup- ported only a discrete set of points. Instead, the par- ticle filter is used to approximate expectation values: if f is a function whose domain is the set of model vec- tors x , we want the particle filter to satisfy Z f ( x ) Pr ( x ) d x ≈ X k w k f ( x k ) . (6) The posterior distribution can also be approximated using a particle filter . In fact, a posterior particle filter can be computed directly from a particle filter for the prior distribution as follows. Let { ( w k , x k ) } be the set of weights and locations for a particle filter for some prior distribution Pr ( x ) . W e then compute a particle filter { ( w 0 k , x 0 k ) } for the posterior distribution by set- ting x 0 k = x k and w 0 k = w k Pr ( D | x k ) P j w j Pr ( D | x j ) , (7) where D is the data set used in the Bayesian update. In practice, updating the weights in this fashion causes the particle filter to become unstable as data is col- lected; by default, QInfer will periodically apply the Liu–W est algorithm to restor e stability [ 25 ]. See Ap- pendix B for details. 3 At any point during the processing of data, the ex- pectation of any function with respect to the posterior is approximated as E [ f ( x ) | D ] ≈ X k w k ( D ) f ( x k ) . (8) In particular , the expected error in x is given by the posterior covariance, Cov ( x | D ) : = E [ xx T | D ] − E [ x | D ] E T [ x | D ] . This can be used, for instance, to adaptively choose experiments which minimize the posterior variance [ 24 ]. This approach has been used to exponentially improve the number of samples re- quired in frequency estimation problems [ 19 , 20 ], and in phase estimation [ 26 , 27 ]. Alternatively , other cost functions can be considered, such as the information gain [ 23 , 28 ]. QInfer allows for quickly computing ei- ther the expected posterior variance or the informa- tion gain for proposed experiments, making it straight- forward to develop adaptive experiment design proto- cols. The functionality exposed by QInfer follows a sim- ple object model, in which the experiment is described in terms of a model , and background information is de- scribed in terms of a prior distribution . Each of these classes is abstract , meaning that they define what be- havior a QInfer user must specify in order to fully specify an inference procedure. For convenience, QIn- fer provides several pre-defined implementations of each, as we will see in the following examples. Con- crete implementations of a model and a prior distribu- tion are then used with SMC to update the prior based on data. In summary , the iterative approach described above is formalized in terms of the following Python object model: class qinfer.Distribution: abstr act sample( n ): Returns n samples from the repr esented distribution. class qinfer.Mo del: abstr act lik eliho o d( d , x , e ): Returns an evalu- ation of the likelihood function Pr ( d | x ; e ) for a single datum d , a vector of model parame- ters x and an experiment e . abstr act are mo dels v alid( x ): Evaluates whether x is a valid assignment of model parameters. class qinfer.SMCUp dater: up date( d , e ): Computes the Bayes update (7) for a single datum (that is, D = { d } ). est mean(): Returns the current estimate ˆ x = E [ x ] . A complete description of the QInfer object model can be found in the Users’ Guide. Notably qinfer . SMCUpdater r elies only on the behavior specified by each of the abstract classes in this object model. Thus, it is straightforwar d for the user to specify their own prior and likelihood function by either im- plementing these classes (as in the example of Ap- pendix A ), or by using one of the many concrete im- plementations provided with QInfer . The concrete implementations provided with QIn- fer are useful in a range of common applications, as de- scribed in the next Section. W e will demonstrate how these classes ar e used in practice with examples drawn from quantum information applications. W e will also consider the qinfer . Heuristic class, which is useful in con- texts such as online adaptive experiments and simu- lated experiments. 3 Applications in Quantum Information In this Section, we describe various possible applica- tions of QInfer to existing experimental protocols. In doing so, we highlight both functionality built-in to QInfer and how this functionality can be readily ex- tended with custom models and distributions. W e be- gin with the problems of phase and frequency learning, then describe the use of QInfer for state and process to- mography , and conclude with applications to random- ized benchmarking. 3.1 Phase and F requency Learning One of the primary applications for particle filtering is for learning the Hamiltonian H under which a quan- tum system evolves [ 24 ]. For instance, consider the single-qubit Hamiltonian H = ω σ z / 2 for an unknown parameter ω . An experiment on this qubit may then consist of preparing a state | + i = ( | 0 i + | 1 i ) / √ 2 , evolving for a time t and then measuring in the σ x ba- sis. This model commonly arises from Ramsey inter- ferometry , and gives a likelihood function Pr ( 0 | ω ; t ) =   h + | e − i ω tσ z / 2 | + i   2 = cos 2 ( ω t / 2 ) . (9) Note that this is also the same model for Rabi inter- ferometry as well, with the interpretation of H as drive term rather than the internal Hamiltonian for a system. Similarly , this model forms the basis of Bayesian and maximum likelihood approaches to phase estimation. In any case, QInfer implements (9) as the SimplePrecessionModel class, making it easy to quickly perform Bayesian inference for Ramsey or Rabi esti- mation problems. W e demonstrate this in Listing 1 , 4 0 20 40 60 80 100 # of Experiments Performed 0.4 0.5 0.6 0.7 0.8 0.9 ω Est. True Figure 1: Frequency estimate obtained using Listing 1 as a function of the numb er of exp eriments p erformed. using ExpSparseHeuristic to select the k th measurement time t k = ( 9 / 8 ) k , as suggested by analytic arguments [ 19 ]. Listing 1: Frequency estimation example using SimplePrecessionMo del . > > > fr o m q i n f e r i m p o r t * > > > m o d e l = S i m p l e P r e c e s s i o n M o d e l ( ) > > > p r i o r = U n i f o r m D i s t r i b u t i o n ( [ 0 , 1 ] ) > > > n p a r t i c l e s = 2 0 0 0 5 > > > n e x p e r i m e n t s = 1 0 0 > > > u p d a t e r = SM CUp dat er ( m o d el , n p a r t i c l e s , p r i o r ) > > > h e u r i s t i c = E x p S p a r s e H e u r i s t i c ( u p d a t e r ) > > > t r u e p a r a m s = p r i o r . s a m p l e ( ) > > > f o r i d x e x p e r i m e n t i n r a n g e ( n e x p e r i m e n t s ) : 10 . . . e x p e r i m e n t = h e u r i s t i c ( ) . . . dat um = mo d e l . s i m u l a t e e x p e r i m e n t ( t r u e p a r a m s , e x p e r i m e n t ) . . . u p d a t e r . u p d a t e ( da tum , e x p e r i m e n t ) > > > p r i n t ( u p d a t e r . e s t m e a n ( ) ) More complicated models for learning Hamiltonians with particle filtering have also been considered [ 29 – 32 ]; these can be readily implemented in QInfer as cus- tom models by deriving from the Model class, as de- scribed in Appendix A . 3.2 State and Pro cess T omography Though originally conceived of as a algebraic inverse problem, quantum tomography is also a pr oblem of parameter estimation. Many have also considered the problem in a Bayesian framework [ 33 , 34 ] and the se- quential Monte Carlo algorithm has been used in both theoretical and experimental studies [ 23 , 28 , 35 – 37 ]. T o define the model, we start with a basis for trace- less Hermitian operators { B j } d 2 − 1 j = 1 . In the case of a qubit, this could be the basis of Pauli matrices, for ex- ample. Then, any state ρ can be written ρ = 1 d + d 2 − 1 X j = 1 θ j B j , (10) for some vector of parameters θ . These parameters must be constrained such that ρ ≥ 0 . In the simplest case, we can consider two-outcome measurements represented by the pair { E , 1 − E } . The Born rule defines the likelihood function Pr ( E | ρ ) = T r ( ρE ) . (11) For multiple measurements, we simply iterate. For many trials of the same measurement, we can use a derived model as discussed below . QInfer’s T omographyModel abstracts many of the im- plementation details of this problem, exposing tomo- graphic models and estimates in terms of QuT iP’s Qobj 5 1 0 1 T r ( σ x ρ ) 1 0 1 T r ( σ z ρ ) C r e d i b l e R e g i o n ( α = 0 . 9 5 ) Posterior True Posterior Mean Figure 2: Posterio r over rebit states after 100 random Pauli measurements, each rep eated five times, as implemented by Listing 2 . class [ 38 ]. This allows for readily integrating QIn- fer functionality with that of QuT iP, such as fidelity metrics, diamond norm calculation, and other such manipulations. T omography support in QInfer requir es one of the bases mentioned above in order to parameterize the state. Many common choices of basis are included as T omographyBasis objects, such as the Pauli or Gell-Mann bases. Many of the most commonly used priors are al- ready implemented as a QInfer Distribution . Listing 2: Rebit state tomography example using T omographyModel . > > > fr o m q i n f e r i m p o r t * > > > fr o m q i n f e r . t o m o g r a p h y i m p o r t * > > > b a s i s = p a u l i b a s i s ( 1 ) # S i n g l e − q u b i t P a u l i b a s i s . > > > m o d e l = T om o gr a ph y Mo d el ( b a s i s ) 5 > > > p r i o r = G i n i b r e R e d i t D i s t r i b u t i o n ( b a s i s ) > > > u p d a t e r = SM CUp dat er ( m o d el , 8 0 0 0 , p r i o r ) > > > h e u r i s t i c = R a n d o m P a u l i H e u r i s t i c ( u p d a t e r ) > > > t r u e s t a t e = p r i o r . s a m p l e ( ) > > > 10 > > > f o r i d x e x p e r i m e n t i n r a n g e ( 5 0 0 ) : > > > e x p e r i m e n t = h e u r i s t i c ( ) > > > dat um = mo d e l . s i m u l a t e e x p e r i m e n t ( t r u e s t a t e , e x p e r i m e n t ) > > > u p d a t e r . u p d a t e ( da tum , e x p e r i m e n t ) For simulations, common randomized measure- ment choices are already implemented. For exam- ple, RandomPauliHeuristic chooses random Pauli mea- surements for qubit tomography . In Listing 2 , we demonstrate QInfer ’s tomography support for a rebit . By analogy to the Bloch sphere, a rebit may be represented by a point in the unit disk, making rebit tomography useful for plotting exam- ples. More generally , with different choices of basis, QInfer can be used for qubits or higher-dimensional states. For example, recent work has demonstrated the use of QInfer for tomography procedures on seven- dimensional states [ 39 ]. Critically , QInfer provides a region estimate for this example, describing a region that has a 95% probability of containing the true state. W e will explore r egion estimation further in Section 4.1 . Finally , we note that pr ocess tomography is a spe- cial case of state tomography [ 36 ], such that the same functionality described above can also be used to ana- lyze process tomography experiments. In particular , the qinfer . ProcessT omographyHeuristic class repr esents the experiment design constraints imposed by process to- mography , while BCSZChoiDistribution uses the distribu- tion over completely positive trace-preserving maps proposed by Bruzda et al. [ 40 ] to repr esent a prior dis- tribution over the Choi states of random channels. 6 3.3 Randomized Benchmarking In recent years, randomized benchmarking (RB) has reached a critical role in evaluating candidate quantum information processing systems. By using random se- quences of gates drawn from the Clifford group, RB provides a likelihood function that depends on the fi- delity with which each Cliffor d group element is im- plemented, allowing for estimates of that fidelity to be drawn from experimental data [ 4 ]. In particular , suppose that each gate is implemented with fidelity F , and consider a fixed initial state and measurement. Then, the survival probability over se- quences of length m is given by [ 41 ] Pr ( survival | p , A , B ; m ) = Ap m + B , (12) where p : = ( dF − 1 ) / ( d − 1 ) , d is the dimension of the system under consideration, and where A and B de- scribe the state preparation and measurement (SP AM) errors. Learning the model x = ( p , A , B ) thus provides an estimate of the fidelity of interest F . The likelihood function for randomized benchmark- ing is extremely simple, and requires only scalar arithmetic to compute, making it especially useful for avoiding the computational overhead typically re- quired to characterize large quantum systems with classical r esources. Multiple generalizations of RB have been recently developed which extend these ben- efits to estimating cr osstalk [ 42 ], coher ence [ 43 ], and to estimating fidelities of non-Cliffor d gates [ 44 , 45 ]. RB has also been extended to provide tomographic infor- mation as well [ 46 ]. The estimates provided by ran- domized benchmarking have also been applied to de- sign improved contr ol sequences [ 47 , 48 ]. QInfer supports RB experiments through the qinfer . RandomizedBenchmarkingModel class. For common priors, QInfer also provides a simplified interface, qinfer . simple est rb , that reports the mean and covariance over an RB model given experimental data. W e pro- vide an example in Listing 3 . Listing 3: Randomized b enchmarking example using simple est rb . > > > fr o m q i n f e r i m p o r t * > > > i m p o r t nu mp y a s n p > > > p , A , B = 0 . 9 5 , 0 . 5 , 0 . 5 > > > m s = np . l i n s p a c e ( 1 , 8 0 0 , 2 0 1 ) . a s t y p e ( i n t ) 5 > > > s i g n a l = A * p * * m s + B > > > n s h o t s = 2 5 > > > c o u n t s = n p . rand om . b i n o m i a l ( p= s i g n a l , n= n s h o t s ) > > > d a t a = np . c o l u m n s t a c k ( [ c o u n t s , ms , n s h o t s * np . o n e s l i k e ( c o u n t s ) ] ) > > > mean , c o v = s i m p l e e s t r b ( d a t a , n p a r t i c l e s = 1 2 0 0 0 , p mi n = 0 . 8 ) 10 > > > p r i n t ( mea n , np . s q r t ( n p . d i a g ( c o v ) ) ) 4 A dditional Functionalit y Having introduced common applications for QInfer, in this Section we describe additional functionality which can be used with each of these applications, or with custom models. 4.1 Region Estimation and Error Ba rs As an alternative to specifying the entire posterior dis- tribution approximated by qinfer . SMCUpdater , we pro- vide methods for reporting credible regions over the posterior , based on covariance ellipsoids, convex hulls, and minimum volume enclosing ellipsoids [ 35 ]. These region estimators provide a rigor ous way of summa- rizing one’s uncertainty following an experiment (col- loquially referred to as “error bars”), and owing to the Bayesian approach, do so in a manner consistent with experimental experience. Posterior credible regions can be found by using the SMCUdater.est credible region method. This method returns a set of particles such that the sum of their weights corresponding weights is at least a specified ratio of the total weight. For example, a 95% credible regions is repr esented as a collection of particles whose weight sums to at least 0.95. This does not necessarily admit a very compact de- scription since many of the particles would be interior to the regions. In such cases, it is useful to find region estimators containing all of the particles describing a credible region. The SMCUpdater.region est hull method does this by finding a convex hull of the credible parti- cles. Such a hull is depicted in Figure 2 . The convex hull of an otherwise random set of points is also not necessarily easy to describe or intuit. In such cases, SMCUdpater.region est ellipsoid finds the minimum- volume enclosing ellipse (MVEE) of the convex hull region estimator . As the name suggests, this is the smallest ellipsoid containing the credible particles. It is strictly larger than the hull and thus maintains cred- ibility . Ellipsoids are specified by their center and covariance matrix. V isualizing the covariance ma- trix can also usually provide important diagnostic in- formation, as in Figure 3 . In that example, we can quickly see that the p and A parameters estimated from a randomized benchmarking experiment are anti- 7 0.86 0.88 0.90 0.92 0.94 0.96 0.98 1.00 p 0 10 20 30 40 50 Posterior True p A B p A B Figure 3: (Left) Posterio r over the randomized b enchmarking decay rate parameter p after measuring 25 sequences at each of 201 sequence lengths, as describ ed in Listing 3 . (Right) The posterior covariance matrix over all three randomized benchmarking pa rameters x = ( p , A , B ) , rep resented as a Hinton diagram. White squa res indicate p ositive elements, while black squares indicate negative elements, and the relative sizes indicate magnitude of each element. correlated, such that we can explain more preparation and measurement errors by assuming better gates and vice versa. 4.2 Derived Mo dels QInfer allows for the notion of a model chain , wher e the likelihood of a given model in the chain is a function of the likelihoods of models below it, and possibly new model or experiment parameters. This abstraction is useful for a couple of reasons. It encourages more ro- bust programs, since models in the chain will often be debugged independently . It also often makes writing new models easier since part of the chain may be in- cluded by default in the QInfer library , or may overlap with other similar models the user is implementing. Fi- nally , in quantum systems, it is common to have a like- lihood function which is most naturally expr essed as a hierarchical probability distribution, with base mod- els describing quantum physics, and overlying models describing measurement pr ocesses. Model chains ar e typically implemented through the use of the abstract class DerivedModel . Since this class itself inherits fr om the Model class, subclass instances must provide standard model properties and meth- ods such as likelihood , n outcomes , and modelparam names . Additionally , DerivedModel accepts an ar gument model , referring to the underlying model directly below it in the model chain. Class properties exist for refer enc- ing models at arbitrary depths in the chain, all the way down to the base model . As an example, consider a base model which is the precession model discussed in Section 3.1 . This is a two-outcome model whose outcomes correspond to measuring the state | + i or the orthogonal state |−i , which can be viewed as flipping a biased coin. Per- haps an actual experiment of this system consists of flipping the coin N times with identical settings, wher e the individual results are not recorded, only the total number n + of | + i results. In this case, we can con- catenate this base model with the built-in DerivedModel called BinomialModel . This model adds an additional ex- periment parameter n meas specifying how many times the underlying model’s coin is flipped in a single ex- periment. Listing 4: Frequency estimation with the derived BinomialMo del and a linear time-sampling heuristic. > > > fr o m q i n f e r i m p o r t * > > > i m p o r t nu mp y a s n p > > > m o d e l = B i n o m i a l M o d e l ( S i m p l e P r e c e s s i o n M o d e l ( ) ) 8 0 5 10 15 20 # of Times Sampled (25 measurements/ea) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 ω Est. True Figure 4: F requency estimate after 25 measurements at each of 20 linearly-spaced times, using qinfer . BinomialMo del as in Listing 4 . > > > n m e a s = 2 5 5 > > > p r i o r = U n i f o r m D i s t r i b u t i o n ( [ 0 , 1 ] ) > > > u p d a t e r = SM CUp dat er ( m o d el , 2 0 0 0 , p r i o r ) > > > t r u e p a r a m s = p r i o r . s a m p l e ( ) > > > f o r t i n np . l i n s p a c e ( 0 . 1 , 2 0 , 2 0 ) : . . . e x p e r i m e n t = n p . a r r a y ( [ ( t , n m e a s ) ] , d t y p e = m o d e l . e x p p a r a m s d t y p e ) 10 . . . dat um = m o de l . s i m u l a t e e x p e r i m e n t ( t r u e p a r a m s , e x p e r i m e n t ) . . . u p d a t e r . u p d a t e ( da tum , e x p e r i m e n t ) > > > p r i n t ( u p d a t e r . e s t m e a n ( ) ) Note that parallelization, discussed in Section 4.5 , is implemented as a DerivedModel whose likelihood batches the underlying model’s likelihood function across pr ocessors. 4.3 Time-Dep endent Mo dels So far , we have only consider ed time-independent (pa- rameter estimation) models, but particle filtering is useful more generally for estimating time-dependent (state-space) models. Following the work of Isar d and Blake [ 49 ], when performing a Bayes update, we may also incorporate state-space dynamics by adding a time-step update. For example, to follow a W iener process, we move each particle x i ( t k ) at time t k to its new position x i ( t k + 1 ) = x t ( t k ) + ( t k + 1 − t k ) η , (13) with η ∼ N ( 0, Σ ) for a covariance matrix Σ . Importantly , we need not assume that time- dependence in x follows specifically a W iener process. For instance, one may consider timestep increments describing stochastic evolution of a system undergo- ing weak measurement [ 22 ], such as an atomic en- semble undergoing probing by an optical interferom- eter [ 50 ]. In each case, QInfer uses the timestep incre- ment implemented by the Model.update timestep method, which specifies the time step that SMCUpdater should perform after each datum. This design allows for the specification of more complicated time step updates than the representative example of (13) . For instance, the qinfer . RandomW alkModel class adds diffusive steps to existing models and can be used to quickly learn time-dependent properties, such as shown in Listing 5 . Moreover , QInfer provides the DiffusiveT omographyModel for including time-dependence in tomography by truncating time step updates to lie within the space of valid states [ 36 ]. A video example of time-dependent tomography can be found on Y ouT ube [ 51 ]. In this way , by following the Isard and Blake [ 49 ] al- gorithm, we obtain a very general solution for time- dependence. Importantly , other approaches exist that 9 0 10 20 30 40 50 60 70 t ( µ s ) 0.2 0.3 0.4 0.5 0.6 0.7 0.8 ω ( G H z ) True Estimated Figure 5: Time-dep endent frequency estimation, using qinfer . RandomWalkModel as in Listing 5 . may be better suited for individual problems, includ- ing modifying resampling procedures to incorporate additional noise [ 52 , 53 ], or adding hyperparameters to describe deterministic time-dependence [ 36 ]. Listing 5: Frequency estimation with a time-dep endent mo del. > > > fr o m q i n f e r i m p o r t * > > > i m p o r t nu mp y a s n p > > > p r i o r = U n i f o r m D i s t r i b u t i o n ( [ 0 , 1 ] ) > > > t r u e p a r a m s = np . a r r a y ( [ [ 0 . 5 ] ] ) 5 > > > n p a r t i c l e s = 2 0 0 0 > > > m o d e l = RandomWalkModel ( . . . B i n o m i a l M o d e l ( S i m p l e P r e c e s s i o n M o d e l ( ) ) , N o r m a l D i s t r i b u t i o n ( 0 , 0 . 0 1 * * 2 ) ) > > > u p d a t e r = SM CUp dat er ( m o d el , n p a r t i c l e s , p r i o r ) > > > t = np . p i / 2 10 > > > n m e a s = 4 0 > > > e x p p a r a m s = np . a r r a y ( [ ( t , n m e a s ) ] , d t y p e = m o d e l . e x p p a r a m s d t y p e ) > > > f o r i d x i n r a n g e ( 1 0 0 0 ) : . . . dat um = mo d e l . s i m u l a t e e x p e r i m e n t ( t r u e p a r a m s , e x p p a r a m s ) . . . t r u e p a r a m s = n p . c l i p ( m o d e l . u p d a t e t i m e s t e p ( t r u e p a r a m s , e x p p a r a m s ) [ : , : , 0 ] , 0 , 1 ) 15 . . . u p d a t e r . u p d a t e ( da tum , e x p p a r a m s ) 4.4 P erfo rmance and Robustness T esting One important application of QInfer is pr edicting how well a particular parameter estimation experiment will work in practice. This can be formalized by consider- ing the risk R ( x ) : = E D [( ˆ x ( D ) − x ) T ( ˆ x ( D ) − x ) ] in- curred by the estimate ˆ x ( D ) as a function of some true model x . The risk can be estimated by drawing many differ ent data sets D , computing the estimates for each, and reporting the average error . Similarly , one can es- timate the Bayes risk r ( π ) : = E x ∼ π [ R ( x ) ] by drawing a new “true” model x from a prior π along with each data set. In both cases, QInfer automates the process of per- forming many independent estimation trials thr ough the perf test multiple function. This function will run an updater loop for a given model, prior , and exper - iment design heuristic, returning the errors incurred after each measurement in each trial. T aking an expec- tation value with numpy.mean returns the risk or Bayes risk, depending if the true model keyword argument is set. For example, Listing 6 finds the Bayes risk for a fre- quency estimation experiment ( Section 3.1 ) as a func- tion of the number of measurements performed. Performance evaluation can also easily be paral- 10 0 10 20 30 40 50 # of Experiments Performed 1 0 - 6 1 0 - 5 1 0 - 4 1 0 - 3 1 0 - 2 1 0 - 1 Bayes Risk Figure 6: Bay es risk of a frequency estimation mo del with exp onentially sparse sampling as a function of the numb er of exp eriments p erfo rmed, and as calculated by Listing 6 . lelized over trials, as discussed in Section 4.5 , allow- ing for efficient use of computational resour ces. This is especially important when comparing performance for a range of differ ent parameters. For instance, one might want to consider how the risk and Bayes risk of an estimation procedur e scale with errors in a faulty simulator; QInfer supports this usecase with the qinfer . PoisonedModel derived model, which adds errors to an underlying “valid” model. In this way , QInfer en- ables quickly reasoning about how much approxima- tion error can be tolerated by an estimation pr ocedure. Listing 6: Bay es risk of frequency estimation as a function of the numb er of measurements, calculated using p erf test multiple . > > > p e r f o r m a n c e = p e r f t e s t m u l t i p l e ( . . . # U se 1 0 0 t r i a l s t o e s t i m a t e e x p e c t a t i o n o v e r d a t a . . . . 1 0 0 , . . . # U se a s i m p l e p r e c e s s i o n m o d e l b o t h t o g e n e r a t e , 5 . . . # d a t a , a n d t o p e r f o r m e s t i m a t i o n . . . . S i m p l e P r e c e s s i o n M o d e l ( ) , . . . # U se 2 , 0 0 0 p a r t i c l e s a n d a u n i f o r m p r i o r . . . . 2 0 0 0 , U n i f o r m D i s t r i b u t i o n ( [ 0 , 1 ] ) , . . . # T a k e 5 0 m e a s u r e m e n t s w i t h t k = ab k . 10 . . . 5 0 , E x p S p a r s e H e u r i s t i c . . . ) > > > # T he r e t u r n e d p e r f o r m a n c e d a t a h a s a n i n d e x f o r t h e t r i a l , an d a n i n d e x f o r t h e m e a s u r e m e n t num be r . > > > p r i n t ( p e r f o r m a n c e . s h a p e ) ( 1 0 0 , 5 0 ) 15 > > > # C a l c u l a t e t h e B a y e s r i s k b y t a k i n g a me an o v e r t h e t r i a l i n d e x . > > > r i s k = np . mean ( p e r f o r m a n c e [ ’ l o s s ’ ] , a x i s = 0 ) 4.5 P a rallelization At each step of the SMC algorithm, the likelihood Pr ( d n | x ) of an experimental datum d n is computed for every particle x k in the distribution. T ypically , the to- tal running time of the algorithm is overwhelmingly spent calculating these likelihoods. However , individ- ual likelihood computations are independent of each other and therefore may be performed in parallel. On a single computational node with multiple cores, lim- ited parallelization is performed automatically by rely- ing on NumPy’s vectorization primitives [ 54 ]. More generally , however , if the running time of Pr ( d n | x ) is largely independent of x , we may divide 11 1 2 4 8 12 16 24 # Engines 2 - 4 2 - 3 2 - 2 2 - 1 2 0 Normalized Computation Time Figure 7: Pa rallelization of the likelihoo d function being tested on a single computer with 12 physical Intel Xeon co res. 5000 pa rticles are sha red over a varying numb er of ipypa rallel engines. The linear unit slop e indicates that overhead is negligible in this example. This holds until the number of physical cores is reached, past which hyper-threading continues to give diminishing returns. The single-engine running time w as about 37 seconds, including ten different experiment values, and 5 p ossible outcomes. our particles into L disjoint groups, { x ( 1 ) 1 , ..., x ( 1 ) k 1 } t · · · t { x ( L ) 1 , ..., x ( L ) k L } , (14) and send each group along with d n to a separate pro- cessor to be computed in parallel. In QInfer , this is handled by the derived model ( Sec- tion 4.2 ) qinfer . DirectV iewParallelizedModel which uses the Python library ipyparallel [ 55 ]. This library supports everything from simple parallelization over the cores of a single pr ocessor , to make-shift clusters set up over SSH, to professional clusters using standard job sched- ulers. Passing the model of interest as well as an ipyparallel . DirectV iew of the processing engines is all that is necessary to parallelize a model. Listing 7: Example of parallelizing likelihoo d calls with DirectViewP a rallelizedMo del . > > > fr o m q i n f e r i m p o r t * > > > fr o m i p y p a r a l l e l i m p o r t C l i e n t > > > r c = C l i e n t ( p r o f i l e = ” m y c o r e s ” ) > > > m o d e l = D i r e c t V i e w P a r a l l e l i z e d M o d e l ( S i m p l e P r e c e s s i o n M o d e l ( ) , r c [ : ] ) In Figure 7 , a roughly 12 × speed-up is demonstrated by parallelizing a model over the 12 cores of a sin- gle computer . This model was contrived to demon- strate the parallelization potential of a generic Hamil- tonian learning problem which uses dense operators and states. A single likelihood call generates a random 16 × 16 anti-hermitian matrix (repr esenting the gener- ator of a four qubit system), exponentiates it, and re- turns overlap with the | 0000 i state. Implementation details can be found in the QInfer examples r epository [ 16 ], or in the ancillary files. So far , we have discussed parallelization from the perspective of traditional processors (CPUs), which typically have a small number of processing cores on each chip. By contrast, moderately-priced desktop graphical processing units (GPUs) will often contain thousands of cores, while GPU hosts tailored for sci- entific use can have tens of thousands. This massive parallelization makes GPUs attractive for particle fil- tering [ 56 ]. Using libraries such as PyOpenCL and Py- CUDA [ 57 ] or Numba [ 58 ], custom models can be writ- ten which take advantage of GPUs within QInfer [ 52 ]. For example, qinfer . AcceleratedPrecessionModel offloads its computation of cos 2 to GPUs using PyOpenCL. 12 4.6 Other Features In addition to the functionality described above, QIn- fer has a wide range of other features that we describe more briefly here. A complete description can be found in the provided Users’ Guide (see ancillary files or do cs.qinfer.o rg ). Plotting and Visualization QInfer provides plot- ting and visualization support based on matplotlib [ 59 ] and mpltools [ 60 ]. In particular , qinfer . SMCUpdater pro- vides methods for plotting posterior distributions and covariance matrices. These methods make it straight- forward to visually diagnose the operation of and re- sults obtained from particle filtering. Similarly , the qinfer . tomography module provides sev- eral functions for producing plots of states and distri- butions over rebits (qubits restricted to real numbers). Rebit visualization is in particular useful for demon- strating the conceptual operation of particle filter– based tomography in a clear and attractive manner . Fisher Information Calculation In evaluating es- timation protocols, it is important to establish a base- line of how accurately one can estimate a model even in principle. Similarly , such a baseline can be used to compare between protocols by informing as to how much information can be extracted from a proposed experiment. The Cram ´ er–Rao bound and its Bayesian analog, the van T rees inequality ( a.k.a. the Bayesian Cram ´ er–Rao bound), formalize this notion in terms of the Fisher information matrix [ 61 , 62 ]. For any model which specifies its derivative in terms of a score , QIn- fer will calculate each of these bounds, providing use- ful information about proposed experimental and es- timation protocols. The qinfer . Scor eMixin class builds on this by calculating the score of an arbitrary model us- ing numerical differ entiation. Mo del Selection and A veraging Statistical infer - ence does not requir e asserting a priori the correctness of a particular model (that is, likelihood function), but allows a model to be taken as a hypothesis and com- pared to other models. This is made formal by model selection . From a Bayesian perspective, the ratio of the posterior normalizations for two differ ent models gives a natural and principled model selection crite- rion, known as the Bayes factor [ 63 ]. The Bayes fac- tor provides a model selection rule that is significantly more r obust to outlying data than conventional hy- pothesis testing approaches [ 64 ]. For example, in quan- tum applications, the Bayes factor is particularly use- ful in tomography , and can be used to decide the rank or dimension of a state [ 35 ]. QInfer implements this criterion as the SMCUpdater.normalization recor d property , allowing for model selection and averaging to be per- formed in a straightforward manner . Appro ximate Maximum-Lik eliho o d Estimation As opposed to the Bayesian approach, one may also consider maximum likelihood estimation (MLE), in which a model is estimated as ˆ x MLE : = arg max x Pr ( D | x ) . MLE can be appr oximated as the mean of an artificially tight posterior distribution ob- tained by performing Bayesian inference with a likeli- hood function Pr 0 ( D | x ) related to the true likelihood by Pr 0 ( D | x ) = ( Pr ( D | x ) ) γ (15) for a quality parameter γ > 1 [ 65 ]. Similarly , taking γ < 1 with appropriate resampling parameters allows the user to anneal updates [ 66 ], avoiding the dangers posed by strongly multimodal likelihood functions. In this way , taking γ < 1 is roughly analogous to the use of “reset rule” techniques employed in other filtering algorithms [ 53 ]. In QInfer, both cases are implemented by the class qinfer . MLEModel , which decorates another model in the manner of Section 4.2 . Lik eliho o d-F ree Estimation For some models, ex- plicitly calculating the likelihood function Pr ( D | x ) is intractable, but good approaches may exist for draw- ing new data sets consistent with a hypothesis. This is the case, for instance, if a quantum simulator is used in place of a classical algorithm, as recently proposed for learning in lar ge quantum systems [ 29 ]. In the absence of an explicit likelihood function, Bayesian inference must be implemented in a likelihood-free manner , us- ing hypothetical data sets consistent to form a approx- imate likelihood instead [ 67 ]. This introduces an esti- mation err or which can be modeled in QInfer by using the qinfer . PoisonedModel class discussed in Section 4.4 . Simplified Estimation For the frequency estima- tion and randomized benchmarking examples de- scribed in Section 3 , QInfer provides functions to perform estimation using a “standard” updater loop, making it easy to load data from NumPy-, MA TLAB- or CSV -formatted files. Jup yter In tegration Several QInfer classes, in- cluding qinfer . Model and qinfer . SMCUpdater , integrate with Jupyter Notebook to provide additional in- formation formatted using HTML. Moreover , the qinfer . IPythonProgressBar class provides a progress bar as a Jupyter Notebook widget with a QuT iP-compatible interface, making it easy to r eport on performance test- ing progr ess. 13 MA TLAB/Julia Interoperability Finally , QIn- fer functionality is also compatible with MA TLAB 2016a and later , and with Julia (using the PyCall. jl pack- age [ 68 ]), enabling integration both with legacy code and with new developments in scientific computing. 5 Conclusions In this work, we have presented QInfer , our open- source library for statistical inference in quantum in- formation processing. QInfer is useful for a range of differ ent applications, and can be readily used for cus- tom problems due to its modular and extensible de- sign, addressing a pressing need in both quantum in- formation theory and in experimental practice. Impor- tantly , our library is also accessible, in part due to the extensive documentation that we provide (see ancil- lary files or do cs.qinfer.org ). In this way , QInfer sup- ports the goal of reproducible r esearch by providing open-source tools for data analysis in a clear and un- derstandable manner . A ckno wledgments CG and CF acknowledge funding from the Army Research Office grant numbers W911NF-14-1-0098 and W911NF-14-1-0103, from the Australian Research Council Centre of Excellence for Engineer ed Quantum Systems. CG, CF , IH, and T A acknowledge funding from Canadian Excellence Research Chairs (CERC), Natural Sciences and Engineering Research Council of Canada (NSERC), the Province of Ontario, and In- dustry Canada. JG acknowledges funding from ONR Grant No. N00014-15-1-2167. YS acknowledges fund- ing from ARC Discovery Project DP160102426. CG greatly appreciates help in testing and feedback from Nathan W iebe, Joshua Combes, Alan Robertson, and Sarah Kaiser . References [1] G. M. D’Ariano, M. D. Laurentis, M. G. A. Paris, A. Porzio, and S. Solimeno, “Quantum tomog- raphy as a tool for the characterization of optical devices,” Journal of Optics B: Quantum and Semi- classical Optics 4 , S127 (2002) . [2] J. B. Altepeter , D. Branning, E. Jeffr ey , T . C. W ei, P . G. Kwiat, R. T . Thew , J. L. O’Brien, M. A. Nielsen, and A. G. White, “Ancilla-assisted quan- tum process tomography ,” Phys. Rev . Lett. 90 , 193601 (2003) . [3] J. J. W allman and S. T . Flammia, “Randomized benchmarking with confidence,” New Journal of Physics 16 , 103032 (2014) . [4] C. Granade, C. Ferrie, and D. G. Cory , “Acceler- ated randomized benchmarking,” New Journal of Physics 17 , 013042 (2015) . [5] A. S. Holevo, Statistical Structure of Quantum Theory , edited by R. Beig, J. Ehlers, U. Frisch, K. Hepp, W . Hillebrandt, D. Imboden, R. L. Jaffe, R. Kippenhahn, R. Lipowsky , H. v . L ¨ ohneysen, I. Ojima, H. A. W eidenm ¨ uller , J. W ess, J. Zit- tartz, and W . Beiglb ¨ ock, Lecture Notes in Physics Monographs, V ol. 67 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2001). [6] C. W . Helstr om, Quantum Detection and Estimation Theory (Academic Press, 1976). [7] B. Goldacre, “Scientists ar e hoarding data and it’s ruining medical r esearch,” BuzzFeed (2015) . [8] V . Stodden and S. Miguez, “Best practices for computational science: Software infrastructur e and environments for reproducible and extensi- ble resear ch,” Journal of Open Research Software 2 , e21 (2014) . [9] J. P . A. Ioannidis, “Why Most Published Research Findings Are False,” PLOS Med 2 , e124 (2005) . [10] R. Hoekstra, R. D. Morey , J. N. Rouder , and E.-J. W agenmakers, “Robust misinterpretation of con- fidence intervals,” Psychonomic Bulletin & Re- view , 1 (2014) . [11] J. P . A. Ioannidis, “How to Make More Published Research T rue,” PLOS Med 11 , e1001747 (2014) . [12] Y . R. Sanders, J. J. W allman, and B. C. Sanders, “Bounding quantum gate error rate based on re- ported average fidelity ,” New Journal of Physics 18 , 012002 (2016) . [13] F . P ´ erez and B. E. Granger , “IPython: A System for Interactive Scientific Computing,” Computing in Science and Engineering 9 , 21 (2007) . [14] Jupyter Development T eam, “Jupyter ,” (2016). [15] S. R. Piccolo and M. B. Frampton, “T ools and tech- niques for computational reproducibility ,” Giga- Science 5 , 30 (2016) ; A. de V ries, “Using R with Jupyter Notebooks,” (2015); D. Donoho and V . Stodden, “Reproducible resear ch in the math- ematical sciences,” in The Princeton Companion to Applied Mathematics , edited by N. J. Higham (2015). [16] C. Granade, C. Ferrie, I. Hincks, S. Casagrande, T . Alexander , J. Gross, M. Kononenko, and Y . Sanders, “QInfer Examples,” http://go o.gl/ 4sXY1t . [17] M. Dolfi, J. Gukelberger , A. Hehn, J. Imri ˇ ska, K. Pakrouski, T . F . Rønnow , M. T royer , I. Zintchenko, F . Chirigati, J. Freir e, and 14 D. Shasha, “A model project for repr oducible papers: Critical temperature for the Ising model on a squar e lattice,” (2014) , [cond-mat, physics:physics] . [18] M. V an den Nest, “Simulating quantum com- puters with probabilistic methods,” Quant. Inf. Comp. 11 , 784 (2011), arXiv:0911.1624 . [19] C. Ferrie, C. E. Granade, and D. G. Cory , “How to best sample a periodic probability distribution, or on the accuracy of Hamiltonian finding strate- gies,” Quantum Information Processing 12 , 611 (2013) . [20] A. Sergeevich, A. Chandran, J. Combes, S. D. Bartlett, and H. M. W iseman, “Characterization of a qubit Hamiltonian using adaptive measure- ments in a fixed basis,” Physical Review A 84 , 052315 (2011) . [21] A. Doucet and A. M. Johansen, A T utorial on Par- ticle Filtering and Smoothing: Fifteen Y ears Later (2011). [22] B. A. Chase and J. M. Geremia, “Single-shot pa- rameter estimation via continuous quantum mea- surement,” Physical Review A 79 , 022314 (2009) . [23] F . Husz ´ ar and N. M. T . Houlsby , “Adaptive Bayesian quantum tomography ,” Physical Re- view A 85 , 052120 (2012) . [24] C. E. Granade, C. Ferrie, N. W iebe, and D. G. Cory , “Robust online Hamiltonian learn- ing,” New Journal of Physics 14 , 103013 (2012) . [25] J. Liu and M. W est, “Combined parameter and state estimation in simulation-based filtering,” in Sequential Monte Carlo Methods in Practice , edited by De Freitas and N. Gordon (Springer-V erlag, New Y ork, 2001). [26] D. W . Berry and H. M. W iseman, “Optimal States and Almost Optimal Adaptive Measurements for Quantum Interferometry ,” Physical Review Let- ters 85 , 5098 (2000) . [27] B. L. Higgins, D. W . Berry , S. D. Bartlett, H. M. W iseman, and G. J. Pryde, “Entanglement- free Heisenberg-limited phase estimation,” Na- ture 450 , 393 (2007) . [28] G. I. Struchalin, I. A. Pogorelov , S. S. Straupe, K. S. Kravtsov , I. V . Radchenko, and S. P . Kulik, “Experimental adaptive quantum tomography of two-qubit states,” Physical Review A 93 , 012103 (2016) . [29] N. W iebe, C. Granade, C. Ferrie, and D. G. Cory , “Hamiltonian learning and certification us- ing quantum resour ces,” Physical Review Letters 112 , 190501 (2014) . [30] N. W iebe, C. Granade, C. Ferrie, and D. Cory , “Quantum Hamiltonian learning using imper- fect quantum resour ces,” Physical Review A 89 , 042314 (2014) . [31] N. W iebe, C. Granade, and D. G. Cory , “Quantum bootstrapping via compressed quantum Hamilto- nian learning,” New Journal of Physics 17 , 022005 (2015) . [32] M. P . V . Stenberg, Y . R. Sanders, and F . K. W ilhelm, “Efficient Estimation of Resonant Cou- pling between Quantum Systems,” Physical Re- view Letters 113 , 210404 (2014) . [33] K. R. W . Jones, “Principles of quantum infer ence,” Annals of Physics 207 , 140 (1991) . [34] R. Blume-Kohout, “Optimal, reliable estimation of quantum states,” New Journal of Physics 12 , 043034 (2010) . [35] C. Ferrie, “Quantum model averaging,” New Journal of Physics 16 , 093035 (2014) . [36] C. Granade, J. Combes, and D. G. Cory , “Practical Bayesian tomography ,” New Journal of Physics 18 , 033024 (2016) . [37] K. S. Kravtsov , S. S. Straupe, I. V . Radchenko, N. M. T . Houlsby , F . Husz ´ ar , and S. P . Kulik, “Experimental adaptive Bayesian tomography ,” Physical Review A 87 , 062122 (2013) . [38] A. Pitchford, C. Granade, P . D. Nation, and R. J. Johansson, “QuT iP 4.0.0,” (2015–). [39] C. Granade, C. Ferrie, and S. T . Flammia, “Prac- tical adaptive quantum tomography ,” (2016), arXiv:1605.05039 [quant-ph] . [40] W . Bruzda, V . Cappellini, H.-J. Sommers, and K. ˙ Zyczkowski, “Random quantum operations,” Physics Letters A 373 , 320 (2009) . [41] E. Magesan, J. M. Gambetta, and J. Emerson, “Scalable and robust randomized benchmarking of quantum processes,” Physical Review Letters 106 , 180504 (2011) . [42] J. M. Gambetta, A. D. C ´ orcoles, S. T . Merkel, B. R. Johnson, J. A. Smolin, J. M. Chow , C. A. R yan, C. Rigetti, S. Poletto, T . A. Ohki, M. B. Ketchen, and M. Steffen, “Characterization of Ad- dressability by Simultaneous Randomized Bench- marking,” Physical Review Letters 109 , 240504 (2012) . [43] J. W allman, C. Granade, R. Harper , and S. T . Flammia, “Estimating the coherence of noise,” New Journal of Physics 17 , 113020 (2015) . [44] A. W . Cross, E. Magesan, L. S. Bishop, J. A. Smolin, and J. M. Gambetta, “Scalable ran- domized benchmarking of non-Clif ford gates,” npj Quantum Information 2 , 16012 (2016) , arXiv:1510.02720 . [45] R. Harper and S. T . Flammia, “Estimating the fidelity of T gates using standard interleaved randomized benchmarking,” Quantum Science 15 and T echnology 2 , 015008 (2017) , [quant-ph] . [46] S. Kimmel, M. P . da Silva, C. A. Ryan, B. R. John- son, and T . Ohki, “Robust Extraction of T omo- graphic Information via Randomized Benchmark- ing,” Physical Review X 4 , 011050 (2014) . [47] C. Ferrie and O. Moussa, “Robust and efficient in situ quantum control,” Physical Review A 91 , 052306 (2015) . [48] D. J. Egger and F . K. W ilhelm, “Adaptive Hybrid Optimal Quantum Control for Imprecisely Char- acterized Systems,” Physical Review Letters 112 , 240503 (2014) . [49] M. Isard and A. Blake, “CONDENSA- TION—Conditional Density Propagation for V isual T racking,” International Journal of Com- puter V ision 29 , 5 (1998) . [50] B. A. Chase, B. Q. Baragiola, H. L. Partner , B. D. Black, and J. M. Geremia, “Magnetometry via a double-pass continuous quantum measurement of atomic spin,” Physical Review A 79 , 062107 (2009) . [51] C. Granade, J. Combes, and D. G. Cory , “Prac- tical Bayesian tomography supplementary video: State-space state tomography ,” https://go o.gl/ mkibti (2015). [52] C. E. Granade, Characterization, V erification and Control for Large Quantum Systems , Ph.D. thesis (2015). [53] N. W iebe and C. Granade, “Efficient Bayesian phase estimation,” Physical Review Letters 117 , 010503 (2016) . [54] S. van der W alt, S. C. Colbert, and G. V aroquaux, “The NumPy Array: A Structure for Ef ficient Nu- merical Computation,” Computing in Science & Engineering 13 , 22 (2011) . [55] IPython Development T eam, “Ipyparallel,” (2016). [56] A. Lee, C. Y au, M. B. Giles, A. Doucet, and C. C. Holmes, “On the Utility of Graphics Cards to Per- form Massively Parallel Simulation of Advanced Monte Carlo Methods,” Journal of Computational and Graphical Statistics 19 , 769 (2010) . [57] A. Kl ¨ ockner , N. Pinto, Y . Lee, B. Catanzaro, P . Ivanov , and A. Fasih, “PyCUDA and Py- OpenCL: A Scripting-Based Approach to GPU Run-T ime Code Generation,” Parallel Computing 38 , 157 (2012) . [58] S. K. Lam, A. Pitrou, and S. Seibert, “Numba: A LL VM-based Python JIT Compiler ,” in Proceedings of the Second Workshop on the LL VM Compiler Infras- tructure in HPC , LL VM ’15 (ACM, New Y ork, NY , USA, 2015) pp. 7:1–7:6. [59] J. D. Hunter , “Matplotlib: A 2D graphics environ- ment,” Computing In Science & Engineering 9 , 90 (2007) . [60] T . S. Y u, “Mpltools,” (2015). [61] T . M. Cover and J. A. Thomas, Elements of Infor- mation Theory (W iley-Interscience, Hoboken, N.J., 2006). [62] R. D. Gill and B. Y . Levit, “Applications of the van T rees inequality: A Bayesian Cram ´ er-Rao bound,” Bernoulli 1 , 59 (1995) , mathematical Re- views number (MathSciNet): MR1354456. [63] H. Jeffr eys, The Theory of Probability (Oxford Uni- versity Press, Oxfor d, 1998). [64] W . Edwards, H. Lindman, and L. J. Savage, “Bayesian statistical inference for psychological resear ch,” Psychological Review 70 , 193 (1963) . [65] A. M. Johansen, A. Doucet, and M. Davy , “Par- ticle methods for maximum likelihood estimation in latent variable models,” Statistics and Comput- ing 18 , 47 (2008) . [66] J. Deutscher , A. Blake, and I. Reid, “Articulated body motion capture by annealed particle filter- ing,” in Proceedings IEEE Conference on Computer V ision and Pattern Recognition. CVPR 2000 (Cat. No.PR00662) , V ol. 2 (2000) pp. 126–133 vol.2. [67] C. Ferrie and C. E. Granade, “Likelihood-free methods for quantum parameter estimation,” Physical Review Letters 112 , 130402 (2014) . [68] S. G. Johnson, “PyCall.jl,” (2016). [69] A. Beskos, D. Crisan, and A. Jasra, “On the sta- bility of sequential Monte Carlo methods in high dimensions,” The Annals of Applied Probability 24 , 1396 (2014) . [70] M. W est, “Approximating posterior distributions by mixture,” Journal of the Royal Statistical Soci- ety . Series B (Methodological) 55 , 409 (1993) . [71] Y . Sanders, Characterizing Errors in Quantum Infor- mation Processors , Ph.D. thesis , University of W a- terloo (2016). [72] T . Minka, A Family of Algorithms for Approximate Bayesian Inference , Ph.D. thesis , Massachusetts In- stitute of T echnology (2001). [73] C. Granade, “Robust online Hamiltonian learn- ing: Multi- cos 2 model resampling,” https://go o. gl/O2KmEQ (2015). 16 A Custom Mo del Example In Listing 8 , below , we provide an example of a custom subclass of qinfer . FiniteOutcomeModel that implements the likelihood function Pr ( 0 | ω 1 , ω 2 ; t 1 , t 2 ) = cos 2 ( ω 1 t 1 / 2 ) cos 2 ( ω 2 t 2 / 2 ) (16) for model parameters x = ( ω 1 , ω 2 ) and experiment parameters e = ( t 1 , t 2 ) . A more efficient implementation of this model using NumPy vectorization is presented in mor e detail in the Users’ Guide. Listing 8: Example of a custom FiniteOutcomeMo del sub class implementing the multi- cos lik eliho o d (16) . f r o m q i n f e r i m p o r t F i n i t e O u t c o m e M o d e l i m p o r t nu m p y a s np c l a s s M u l t i C o s M o d e l ( F i n i t e O u t c o m e M o d e l ) : 5 @ p r o p e r t y d e f n m o d e l p a r a m s ( s e l f ) : r e t u r n 2 10 @ p r o p e r t y d e f i s n o u t c o m e s c o n s t a n t ( s e l f ) : r e t u r n T r u e d e f n o u t c o m e s ( s e l f , e x p p a r a m s ) : 15 r e t u r n 2 d e f a r e m o d e l s v a l i d ( s e l f , m o d e l p a r a m s ) : r e t u r n n p . a l l ( np . l o g i c a l a n d ( m o d e l p a r a m s > 0 , m o d e l p a r a m s < = 1 ) , a x i s = 1 ) 20 @ p r o p e r t y d e f e x p p a r a m s d t y p e ( s e l f ) : r e t u r n [ ( ’ t s ’ , ’ f l o a t ’ , 2 ) ] d e f l i k e l i h o o d ( s e l f , o u t c o m e s , m o d e l p a r a m s , e x p p a r a m s ) : 25 s u p e r ( M u l t i C o s M o d e l , s e l f ) . l i k e l i h o o d ( o u t c o m e s , m o d e l p a r a m s , e x p p a r a m s ) p r 0 = n p . em pt y ( ( mo d e l p a r a m s . s h a p e [ 0 ] , ex p p a r a m s . s h a p e [ 0 ] ) ) w 1 , w 2 = m o d e l p a r a m s . T t 1 , t 2 = e x p p a r a m s [ ’ t s ’ ] . T 30 f o r i d x m o d e l i n r a n g e ( m o d e l p a r a m s . s h a p e [ 0 ] ) : f o r i d x e x p e r i m e n t i n r a n g e ( e x p p a r a m s . s h a p e [ 0 ] ) : p r 0 [ i d x m o d e l , i d x e x p e r i m e n t ] = ( np . c o s ( w1 [ i d x m o d e l ] * t 1 [ i d x e x p e r i m e n t ] / 2 ) * 35 np . c o s ( w2 [ i d x m o d e l ] * t 2 [ i d x e x p e r i m e n t ] / 2 ) ) * * 2 r e t u r n F i n i t e O u t c o m e M o d e l . p r 0 t o l i k e l i h o o d a r r a y ( o u t c o m e s , p r 0 ) B Resampling The purpose of this appendix is to of fer a brief discussion of resampling for particle filters. In QInfer, the standard resampler is the one proposed by Liu and W est [ 25 ]. W e begin by motivating the development of resamplers by explaining the problem of impoverishment in a particle filter . W e then describe resamplers by explaining that their effect is to produce a particle filter that approximates a smoothed version of the underlying probability distribution. Finally , we explain the Liu–W est r esampling algorithm. Particle filters are intended to allow for the approximation of the expectation values of functions, but admit an ambiguity between using the locations and the weights to do so. Assuming that the particle locations are primarily responsible for representing the particle filtering approximation, and assuming that we want to approximate the expectation value of a function that is not pathological in some way , the number of particles then serves as a rea- sonable pr oxy for the quality of the particle filter . This follows fr om exactly the same ar gument as for Monte Carlo 17 integration, as in this case, the particle locations can be seen to dir ectly correspond to samples of an integrand. On the other hand, if either of these assumptions are violated, one cannot trust the numerical answers obtained using the particle filter method without an additional argument. Since the weights will in general become less even as Bayesian inference pr oceeds, we will rely on r esampling to provide us with precisely such an ar gument. More precisely , the purpose of resampling is to mitigate against the loss in numerical stability caused by having a large number of low-weight particles. If a particle has a small weight, we could neglect it from the computation of an expectation value without introducing a large difference in the result, such that it no longer contributes to the approximation quality of (5) . That is, the particle’s ef fectiveness at contributing to the stability of the algorithm decreases as its weight decr eases. This observation then motivates using the effective sample size n ess : = 1 P k w 2 k (17) as a criterion to ensure the numerical stability of particle filtering [ 69 ]. In the case of roughly equal weights, the effective sample size n ess is roughly equal to the actual number of particles n . If the weights are distributed rather unevenly , however , n ess  n . In the latter case, we say that the particle filter is impoverished . Notably , if a particle filter becomes extremely impoverished, it may not be feasible to effectively recover numerical stability , such that the minimum observed value of n ess serves as a diagnostic criteria for when the numerical approximations used by particle filtering have failed [ 39 ]. For this reason, QInfer will by default call the resampler when n ess falls below half of its initial value, and will warn if n ess ≤ 10 is ever observed. Impoverishment is the result of choosing particle locations according to prior information alone; with near certainty , the posterior distribution will be tightly centered away from all initial particle locations such that some re-discr etization will be needed to represent the final posterior . The goal of resampling is therefor e to modify the choice of particle locations not based on the interpretation of new data, but rather to ensure that the particles are concentrated so as to accurately repr esent the probability distribution of interest. A resampling algorithm is then any algorithm designed to replace an input particle filter that may be impoverished with a new particle filter that is not impoverished but approximates the same probability distribution. Since such a procedur e cannot exist in full generality , each resampling procedure works by enforcing that a particular set of invariants remains true before and after the resampling step. Returning to the particular case of the Liu–W est resampler [ 25 ] employed by default in QInfer, we note that the Liu–W est algorithm is particularly simple to understand from the perspective of kernel density estimation [ 70 ]. W e will not provide a complete algorithm here; such explicit algorithms can be found many places in the literature, including Algorithm 2 of Granade [ 52 ] and Algorithm 2.5 of Sanders [ 71 ]. Rather , we explain that the Liu–W est algorithm acts by preserving the first two moments (the expectation and the variance) of its input distribution. The Liu–W est algorithm starts by using that, as in the celebrated method of kernel density estimation, each sample from a continuous distribution can be used to infer properties about the neighborhood around that sample, given moderate assumptions on the smoothness of a distribution. Thus, the particle filter might be used to define a direct approximation to the true probability distribution by first defining some function K , called a kernel , such that the estimated distribution X k w k K ( x − x k ) (18) is a good approximation of the true probability distribution. This then leaves open the question of what kernel function K should be chosen. Owing to its generality , choosing K to be a normal distribution works well under only mild assumptions, leaving us to choose the variance of the kernel. Specializing to the single-parameter case for simplicity of notation, we denote Gaussian kernels as K ( x ; σ 2 ) , wher e σ 2 is the variance to be chosen. The multi-parameter case follows by replacing this variance by a covariance matrix. The key insight at the core of the Liu–W est algorithm is the observation that the posterior distribution will narrow over time, so that we can choose the variance of the kernel to narrow over time in proportion to the variance of the distribution being approximated. In particular , let a ∈ [ 0, 1 ] be a parameter and let h such that a 2 + h 2 = 1 . Then choose the kernel to be K ( x ; h 2 V [ x ] ) . Since this alone would increase the variance of the distribution under approximation to 1 + h 2 , the Liu–W est resampler instead draws new particle locations x 0 from the distribution x 0 ∼ X k w k K ( x − ( ax k + ( 1 − a ) E [ x ] ; h 2 V [ x ] ) . (19) 18 This distribution contracts each original particle towards the mean by 1 − a , such that the mean and variance of the post-resampling distribution ar e identical to the distribution being approximated. Importantly , the Liu–W est generalizes existing r esampling algorithms, such that the bootstrap [ 21 ] and assumed density filtering resamplers [ 53 , 72 ] are given by a = 1 and a = 0 , respectively . W e also note that violating the invariant that V [ x ] is preserved can allow for some robustness to multimodal distributions and time-dependence [ 52 ]. Finally , a video of the Liu–W est r esampler applied to Bayesian infer ence on the model of Listing 8 is available online [ 73 ]. 19