Adaptive Bayesian and frequentist data processing for quantum tomography Giacomo M. D’Arianoa,b, Daniele F. Magnania,c, Paolo Perinottia,b 9 0 a QUIT Group, Dipartimento di Fisica ”A. Volta”, via Bassi 6, 27100 Pavia, Italy 0 b Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Pavia, via Bassi 6, 27100 2 Pavia, Italy n c Consorzio Nazionale Interuniversitario per le Scienze Fisiche della Materia (CNISM), a J unita` di Pavia, via Bassi 6, 27100 Pavia, Italy 4 1 ] h Abstract p - t Theoutcomestatistics ofaninformationally completequantummeasurement n a forasystem inagivenstatecanbeusedtoevaluatetheensemble expectation u of any linear operator in the same state, by averaging a function of the q [ outcomes that depends on the specific operator. Here we introduce two novel data-processing strategies, non-linear in the frequencies, which lead to 2 v faster convergence to theoretical expectations. 8 5 Key words: quantum estimation, quantum tomography 0 PACS: 03.65.Wj, 03.65.-w, 03.67.-a 5 . 7 0 8 1. Introduction 0 : v In Quantum Mechanics measuring a single observable provides only par- i tial information about the state of the measured system. According to the X Born interpretation, the quantum state is a rule for evaluating the outcome r a probabilities in all conceivable measurements, and a complete information about the quantum state requires a thorough outcome statistics for a quo- rum of observables, or for a suitable informationally complete measurement (shortly info-complete)[1, 2], in conjunction with a suitable data-processing, as it is done in quantum tomography (for a review see Ref. [3]). There are two main classes of approaches in quantum tomography: a) averaging ”patterns functions” a method initiated in Ref. [4]; b) Maximum likelihood techniques [5] Preprint submitted to Elsevier January 14, 2009 Method a) has the advantage of providing any expectation value, e.g. a single density matrix element, without the need of estimating the entire den- sity operator. However, the estimated full matrix is not necessarily positive, which is not a serious drawback, since the non positivity falls within a small fluctuation for large numbers of data. Method b) has the advantage of providing a positive density operator, with smaller fluctuations, however, it has the more serious drawback of need- ing to estimate the full density matrix, while is exponentially large versus the number of systems, and, in the infinite dimensional case needs a dimension- ality cutoff which introduce a bias that is under control only if there is some prior knowledge of the state. In a recent paper [6] the optimal data-processing for evaluating ensemble averages from experimental outcomes was derived for a completely general setting within a Bayesian scheme that assumes a prior probability distri- bution of states. Using as optimality criterion the rate of estimated-to- theoretical convergence of averages, the optimal data-processing itself de- pends on the prior distribution of states. The purpose of the present paper is to exploit the dependence of the optimal data-processing on the prior distribution of states, in order to im- prove the convergence rate using an adaptive data-processing scheme. We will consider info-complete measurements—more generally than a quorum of observables—whose statistics allows to reconstruct all possible ensemble av- erages. Estimation of the quantum state itself is equivalent to the estimation of all possible ensemble averages. We will adopt the natural figure of merit used in Ref. [6], which, in the present context, represents the estimated- to-theoretical convergence rate (in Hilbert-Schmidt distance) of the state. As we will see, exploiting the dependence of the optimal data-processing on the prior state leads to two different data processing strategies, which both improve the convergence rate compared to the standard tomographic procedures, and are easily implementable and computationally efficient: Method 1 (Bayesian iterative procedure): Bayesian update of the prior distribution after the first state reconstruction, with subsequent itera- tion of the optimization. Method 2 (Frequentist approach): replace the theoretical probability distribution of the info-complete in the optimal data-processing with the experimental frequencies. 2 We will see that numerical simulations carried out with both methods show relevant improvement of convergence compared to the plain non adaptive processing of Ref. [6]. The paper is organized as follows. In Section 2 we re-derive the optimal data-processing for given prior distribution of Ref. [6] within an improved theoretical framework. In Sections 3 and 4 we introduce Methods 1 and 2, respectively. Finally, in Section 5, we present numerical simulations for testing both methods in comparison with the original plain non adaptive data-processing, and in Section 6 we end the paper with concluding remarks. 2. Optimization of the data processing InthemodernformulationofQuantumMechanics, thestateofaquantum system associated to a d-dimensional Hilbert space H ∼ Cd is represented by a density matrix, namely a positive operator ρ ≥ 0 with Tr[ρ] = 1. The Born formula provides the probabilities of outcomes in a quantum measurement in terms of the state ρ as follows p(i|ρ) := Tr[ρP ], (1) i where the POVM (P ) (Positive Operator Valued Measure) is a set of (gen- i erally non orthogonal) positive operators P ≥ 0 resolving the identity as i N P = I, thus guaranteeing positivity and normalization of probabilities. i=1 i P ThelinearspanofthePOVMelements Pi, definedasS := Span{Pi}16i6n, is a linear subspace of the space L(H) of linear operators on H, and we will take as a canonical basis in L(H) the operators |mihn|, where |ni is an orthonormal basis thus representing operators X by the vectors of their matrix elements X = hm|x|ni. A POVM is info-complete if S ≡ L(H), m,n namely all operators X ∈ L(H) can be expanded on the POVM elements, and it is possible to determine all ensemble averages hXi , as in Quantum ρ Tomography. For each complex operator X ∈ S the following decomposition holds N X = f [X]P , (2) i i Xi=1 where f [X] is not unique if the set {P } is over-complete. i i With the above expressions we can write the ensemble average of X as follows: N hXi := Tr[Xρ] = f [X]p(i|ρ), (3) ρ i Xi=1 3 with the following statistical error N δX2 := |f [X]|2p(i|ρ)−|hXi |2 (4) ρ i ρ (cid:0) (cid:1) Xi=1 InaBayesianschemeonehasana priori ensembleE := {ρ ,p }ofpossible i i states ρ of the quantum system occurring with probability p . We want i i to minimize the average statistical error on E in the determination of the expectation value of X, namely the variance N δX2 := |f [X]|2p(i|ρ )−|hXi|2 , (5) ε i ε ε (cid:0) (cid:1) Xi=1 where ρ = p ρ and |hXi|2 = p |Tr[ρ X]|2 is the squared modulus ε i i i ε i i i of the expectPation of X averaged ovPer the states in the ensemble (since this term depends only on the ensemble it will be neglected from now on). Using Eq.(1) the first term in Eq.(5) can be rewritten as N Σ (X) := |f [X]|2Tr[P ρ ]. (6) f i i ε Xi=1 Given a POVM (P ), it is possible to define a linear map Λ from an i abstract N-dimensional space K of coefficient vectors c ∈ K to L(H), with range S: N Λc = c P , (7) i i Xi=1 so that using thecanonical basis inK, Λ has matrix elements Λ = (P ) . mn,i i mn A generalized inverse (shortly g-inverse) of Λ is any matrix Γ representing linear operators from L(H) to K such that the following identity holds ΛΓΛ = Λ (8) Notice that the matrix elements (Γ ) of Γ define a set of operators D with i,mn i matrix elements (D ) := Γ∗ . The role of g-inverse Γ is assessed by the i mn i,mn two following important theorems Theorem 1. The following statements are equivalent 4 1. Γ is a g-inverse of Λ 2. For all y ∈ Rng(Λ), x = Γy is a solution of the equation Λx = y. Proof. See Ref. [7]. Theorem 2. For all g-inverse Γ of Λ all solutions of Λx = y are of the form x = Γy+(I −ΓΛ)z, (9) with arbitrary z. Proof. See Ref. [7]. We now define a norm in K as follows N 2 2 ||c||π := |ci| πii, (10) Xi=1 where π = δ π is a positive matrix which is diagonal in the canonical basis ij ij ii in K. In terms of π we define the minimum norm g-inverses Γ that satisfy [8] πΓΛ = Λ†Γ†π. (11) Notice that the present definition of minimum norm g-inverse requires that the norm is induced by a scalar product (in our case ~a·~b := N a∗π b ). i=1 i ii i We will now prove the following crucial theorem P Theorem 3. The following assertions are equivalent 1. Γ is a minimum norm g-inverse of Λ. 2. For all y ∈ Rng(Λ), x = Γy is a solution of the equation Λx = y with minimum norm. Proof. We first prove that 1 ⇒ 2. For Γ g-inverse of Λ, one has due to Theorem 2 ||Γy+(I −ΓΛ)z||2 π =[y†Γ† +z†(I −Λ†Γ†)]π[Γy+(I −ΓΛ)z] (12) =||Γy||2 +||(I −ΓΛ)z||2 π π +z†(I −Λ†Γ†)πΓy+y†Γ†π(I −ΓΛ)z. 5 Since by hypothesis y ∈ Rng(Λ), then y = Λu for some u in K. For a minimum norm g-inverse Γ as in the hypothesis, due to Eq. (11) one has z†(I −Λ†Γ†)πΓΛu+u†Λ†Γ†π(I −ΓΛ)z = (13) z†(I −Λ†Γ†)Λ†Γ†πu+u†πΓΛ(I −ΓΛ)z = 0. where the last equality is due to Eq. (8). Finally, this proves that ||Γy+(I −ΓΛ)z||2 = ||Γy||2 +||(I −ΓΛ)z||2 ≥ ||Γy||2. (14) π π π π namely the solution x = Γy is minimum-norm. Now we prove 2 ⇒ 1. If x = Γy is a solution of Λx = y for all y ∈ Rng(Λ), byTheorem 1 Γ isag-inverse ofΛ, namely ΛΓΛ = Λ. Then ifΓy isminimum norm solution of |Λx = y| then due to Theorem 2 ||Γy||2 ≤ ||Γy+(I −ΓΛ)z||2 (15) π π for all y ∈ Rng(Λ) and for all z one has 0 ≤ ||(I −ΓΛ)z||2 +z†(I −Λ†Γ†)πΓy+y†Γ†π(I −ΓΛ)z. (16) π Since an arbitrary y ∈ Rng(Λ) is Λu for arbitrary u, the second term in Eq. (16) becomes z†(I −Λ†Γ†)πΓΛu+u†Λ†Γ†π(I −ΓΛ)z = (17) 2ℜ z†(I −Λ†Γ†)πΓΛu . (cid:0) (cid:1) Let us keep z fixed and multiply u by an arbitrary α. If the expression in Eq. (17) is not vanishing then taking |α| sufficiently large, for suitable phase one can contradict the bound in Eq. (16), hence ℜ z†(I −Λ†Γ†)πΓΛu = 0 for all u and z and by the same reasoning ℑ z†(I −(cid:0)Λ†Γ†)πΓΛu = 0 f(cid:1)or all u and z. We can then conclude that (I −Λ†(cid:0)Γ†)πΓΛ = Λ†Γ†π(I(cid:1)−ΓΛ) = 0, and consequently πΓΛ = Λ†Γ†π (cid:4) Using Eq. (11), and considering that Σ (X) is the norm of the vector f of coefficients f[X] with π = Tr[P ρ ], it has been proved in [6] that the ii i ε minimum noise is achieved by Γ corresponding to the set of operators D i given by N Dopt := ∆ − {[(I −M)π(I −M)]‡πM} ∆ (18) i i ij j Xj=1 6 where∆ isthesetofoperatorscorresponding totheMoore-Penroseg-inverse i Γ of Λ, satisfying the properties mp Γ Λ = Λ†Γ† , Γ ΛΓ = Γ , Γ† Λ† = ΛΓ , (19) mp mp mp mp mp mp mp and M := Γ Λ = M† = M2. The symbol X‡ denotes the Moore-Penrose mp g-inverse of X. It is indeed easy to verify that Γ := Γ −[(I −M)π(I −M)]‡πMΓ (20) opt mp mp satisfies Eq. (11). Notice that being Γ minimum norm independently of opt X, the statistical error is minimized by the same choice Dopt for all operators i X. When a N-outcomes POVM on a d-dimensional Hilbert space H ∼ Cd is info-complete the state ρ can be written as N ρ = D p(i|ρ), (21) i Xi=1 where D corresponds to any g-inverse Γ. It is then possible to reconstruct i any state ρ using the statistics from measurements: N N ρ = p(i|ρ)D ∼= ν Dopt, (22) i i i Xi=1 Xi=1 where ν = ni is the experimental frequency of the i-th outcome, n being i ntot i the number of occurrence of the i-th outcome, and n = n . By the law tot i i of large numbers we have that lim νi = p(i|ρ). HoweverP, the convergence ntot→∞ rateofρ˜toρdependsonthechoiceofD . Itturnsout[9]thatthechoiceDopt, i i corresponding to Γ , istheonewiththefastest convergence (inaverage over opt all possible experimental outcomes) in the Hilbert-Schmidt distance, defined as follows ||ρ˜−ρ ||2 := Tr[(ρ˜−ρ)2]. (23) 2 2 This can be easily proved considering that the Hilbert-Schmidt distance can be written as the sum of the variances δ(|mihn|)2, and all of the summands are minimized by the choice of minimum-norm Γ = Γ . opt 7 3. The Bayesian iterative procedure In this Section we describe the iterative estimation procedure based on the update of the prior information by means of the state reconstruction provided by experimental data. Here we provide an algorithmic description of the procedure, that yields a self-consistent solution: 1. The protocol starts with the choice of a priori ensemble E := {ρ ,p } i i (where ρ are states and p are their prior probabilities), with the cor- i i responding density matrix ρ(0) := ρ(0) = p ρ , e. g. the one of the E i i i uniform ensemble of all pure states ρ(0) = PI/d. 2. Using ρ(0) it is possible to calculate the diagonal matrix with the prob- ability of the different outcomes: π := δ Tr[P ρ(0)] (24) ij ij i 3. Usingπ inEq.(18)wecanfindtheoptimalg-inverseΓ corresponding ij opt to Dopt associated with ρ(0). i 4. Now the initial a priori density matrix ρ(0) ≡ ρE will be updated as follows: N ρ(1) = ν Dopt (25) i i Xi=1 5. If ρ(1) ∼= ρ(0) within a given tolerable error ε then the average input state is ρ˜:= ρ(1) and the procedure stops. 6. Otherwise after setting ρ(0) := ρ(1) the procedure will go back to the step 2. It is important to remark that at each step the matrices ρ(1) and Dopt i are automatically self-adjoint and normalized: Tr[ρ(1)] = 1 since for all i: Tr[Dopt] = 1 [6], however, they are not necessarily positive. i Thisprotocolinprincipleprovidesreliablestatereconstructions, however, its iterative character makes it less efficient than the one introduced in next Section, since at any iterative step one has to calculate the Moore-Penrose g- inverse in Eq. (18), which is typically a time-consuming operation, especially for POVM’s with a large number N of outcomes. 8 4. The frequentist approach In this Section we introduce the second processing strategy, based on the substitution of prior probabilities by experimental frequencies in Eq. (11). While the previous protocol is essentially a Bayesian update, in this case the the processing relies on the law of large numbers, namely on the fact that lim ν = p(i|ρ),wherethelimithastobeunderstoodinprobability. We ntot→∞ i name this approach frequentist because it fits the frequentist interpretation of probabilities as approximations of experimental frequencies, avoiding prior probabilities, which are the signature of the Bayesian approach. If we substitute the metric matrix π in the Eq. (10) with the diagonal matrix of the frequencies ν , we get: i νΓΛ = Λ†Γ†ν (26) andfollowingthesameproofasforEq.(18)weobtainthefollowingexpression of the optimal g-inverse Γ satisfying condition Eq. (26), in terms of the ν (ν) corresponding operators D i N D(ν) := ∆ − {[(I −M)ν(I −M)]‡νM} ∆ (27) i i ij j Xj=1 that is non linear in the outcomes frequencies due to the Moore-Penrose g-inverse of (I −M)ν(I −M). This protocol has the advantage that it requires only one evaluation of Moore-Penrose g-inverse, and it is then much faster—in terms of computa- tional resources—than the iterative one introduced in the previous Section. (ν) However, here generally Tr[D ] 6= 1, whence in addition to positivity of the i estimated state ρ˜, also the normalization constraint is lost (but not hermitic- ity). 5. Numerical simulations In order to test these two methods and to compare their performances with the plain un-updated procedure some Monte Carlo simulation have ben performed. As an example, we considered the info-complete POVM com- posed by the following six elements 1 P = (I ±σ ), (28) ±i i 6 9 σ = I and ~σ = (σ ,σ ,σ ) denoting the usual Pauli matrices. The theoreti- 0 x y z cal state is 4 1 + i 1 2 2 3 ρ = 5 7 3 = I + σ − σ + σ . (29) (cid:18)1 − i 1 (cid:19) 2 (cid:18) 7 x 3 y 5 z(cid:19) 7 3 5 The simulation consists in 1000 experiments, each consisting in 1000 single-shot measurements, simulated by POVM events extraction according to the theoretical probabilities p(±i|ρ) := Tr[P ρ]. The number of iterations ±i in the Bayesian processing is 10. Figure 1: Histograms representing the number of experiments versus the Hilbert-Schmidt distance of the resulting state from the theoretical one. Upper plot: the light gray bars correspond to the Bayesian processing, the dark grey correspond to the plain processing withoutupdating,thewhitepartistheoverlap. Lowerplot: thedarkgreybarscorrespond tothefrequentistprocessingmethod. Bothplotsshowawellvisibleshiftofthehistograms corresponding to the new adaptive methods towards small errors compared to the plain processing without update. [For other data concerning plots see text.] In Fig. 1 we show the histograms representing the number of experiments 10