Adaptive measurement filter: efficient strategy for optimal estimation of quantum Markov chains

Continuous-time measurements are instrumental for a multitude of tasks in quantum engineering and quantum control, including the estimation of dynamical parameters of open quantum systems monitored through the environment. However, such measurements do not extract the maximum amount of information available in the output state, so finding alternative optimal measurement strategies is a major open problem. In this paper we solve this problem in the setting of discrete-time input-output quantum Markov chains. We present an efficient algorithm for optimal estimation of one-dimensional dynamical parameters which consists of an iterative procedure for updating a `measurement filter' operator and determining successive measurement bases for the output units. A key ingredient of the scheme is the use of a coherent quantum absorber as a way to post-process the output after the interaction with the system. This is designed adaptively such that the joint system and absorber stationary state is pure at a reference parameter value. The scheme offers an exciting prospect for optimal continuous-time adaptive measurements, but more work is needed to find realistic practical implementations.


Introduction
The quantum input-output (I-O) formalism is an effective framework for describing the evolution, monitoring and control of Markovian quantum open systems [1][2][3].In this setting, the interaction with the environment is modelled by coupling the system of interest with a quantum transmission line (channel) represented by a Gaussian bosonic field.The output field carries information about the open system's dynamics which can be accessed by performing continuous-time measurements, and the corresponding conditional system evolution is described in terms of stochastic Schrödinger or filtering equations [4][5][6][7][8].
While these theories are key to quantum engineering applications, they rely on the precise knowledge of the system's dynamical parameters (e.g.Hamiltonian of field coupling), which are often uncertain, or completely unknown, and therefore need to be estimated from measurement data.The I-O formalism is ideally suited for this statistical inference task, and more generally for implementing online system identification methods [9].Unlike direct measurement techniques which require repeated system re-preparations and fast control operations [10][11][12][13][14][15][16], the parameters can be estimated continuously from the output measurement trajectory, even if the system is not directly accessible or it is involved in an information processing task.The first investigation in parameter estimation for continuously-observed quantum systems considered the estimation of the Rabi frequency of an atom in a cavity mode, while a photon counting measurement is performed on the cavity output [17].Subsequent works have addressed a variety of related problems including the dependence on measurement choice [18], adaptive estimation [19,20] filtering with uncertain parameters [21] particle filters for estimation [22,23] achieving Heisenberg scaling [24][25][26], sensing with error correction [27], Bayesian estimation [28][29][30][31][32], quantum smoothing [33][34][35][36], waveform estimation [37,38] estimation of linear systems [39][40][41][42], classical and quantum Fisher informations of the output channel [28,40,[43][44][45][46][47][48]].
An upshot of these studies is that standard measurements such as counting, homodyne or heterodyne generally do no achieve the ultimate limit given by the quantum Cramér-Rao bound [49][50][51], while the optimal measurement prescribed by the symmetric logarithmic derivative requires collective operations on the output state.In this work we make a first step towards addressing the key issue of devising realistic and statistically effective measurement strategies within the the framework of the I-O theory.By realistic we mean procedures which involve sequential continuoustime measurements (as opposed to general non-separable measurements on the output state), possibly combined with more advanced but theoretically well understood operations such as series connections and feedback [52].
For conceptual clarity we focus primarily on discrete-time dynamics, but we will indicate how the techniques may be extended to continuous-time.In the discrete-time setting, the I-O dynamics consists of a d-dimensional system of interest interacting sequentially with a chain of k-dimensional 'noise' input units, which are identically and independently prepared in a state |χ , cf. Figure 1.We assume that the interaction unitary U θ acting on C d ⊗ C k depends on a parameter θ ∈ R, C d < l a t e x i t s h a 1 _ b a s e 6 4 = " h 2 C i 2 1 P e g 6 b q L P O r / a d c i y 1 F W J g = " > A a k Y 4 n Z X 6 i a I x J h M 8 o j 1 D B Q 6 p c t N 5 6 h k 6 M 4 q P g k i a J z S a q 7 8 3 U h w q N Q 0 9 M 5 m F V M t e J v 7 n 9 R I d X L s p E 3 G i q S C L Q 0 H C k Y 5 Q V g H y m a R E 8 6 k h m E h m s i I y x h I T b Y o q m R K c 5 S + v k k 6 9 5 l z U L u / q l c Z N X k c R T u A U z s G B K 2 j A L b S g D Q Q k P M M r v F l P 1 o v 1 b n 0 s R g t W v n M M f 2 B 9 / g D e u Z I b < / l a t e x i t >

Input Output
System which we would like to estimate by measuring the output after n evolution steps.In principle, this can be done by applying the adaptive, separable measurement scheme developed in [53] to the joint pure system-output state; indeed this has been shown to attain the quantum Cramér-Rao bound.However, while theoretically applicable, the algorithm involves manipulating multi-partite operators, making it unsuitable for processing output states with a large number n of noise units.
In addition, it is not clear how the algorithm can be applied to continuous-time dynamics.Our main contribution is to eliminate these drawbacks by devising a scheme which exploits the intrinsic Markovian structure of the problem.Concretely, we propose an algorithm which finds optimal measurement bases for each of the output units by only performing computations on the space of a doubled-up system and a noise unit, i.e.C d ⊗ C d ⊗ C k .Our algorithm has a similar structure to that of the quantum state filter describing the system's conditional evolution, and can be run in real-time without having to specify the time length n in advance.
While our general algorithm requires measurements on both the output and system in order to achieve finite sample optimality, in Proposition 6.1 we prove that by measuring only the output we incur a loss of Fisher information which is bounded by a constant, independent of the time n.Since the quantum Fisher information scales linearly in time, this implies that output measurement is optimal in the leading contribution to the QFI.
We now describe our scheme in more detail.In the first stage of the protocol we use a small proportion of the output units (of sample size ñ ≈ n 1− with small ) in order to compute a preliminary 'rough estimator' θ 0 of the true parameter θ, by performing a standard sequential measurement.This step is necessary in any quantum estimation problem in which the optimal measurement depends on the unknown parameter [54].In particular, this means that strictly speaking one can only attain optimality in the limit of large sample sizes, as θ − θ 0 decays as n −(1− )/2 thanks to the preliminary estimation stage.In the second (main) stage of the protocol we use θ 0 to design a sequential measurement which achieves the output QFI at θ = θ 0 .Since the first stage insures that θ − θ 0 decays and the QFI is continuous with respect to θ, we find that the overall scheme is asymptotically optimal at any parameter value θ.
The second stage is illustrated in Figure 2: each output unit undergoes a physical transformation (which we call 'quantum post-processing') followed by an adaptive projective measurement whose basis is computed according to the 'measurement filter' algorithm described below.More specifically, after interacting with the system, the post-processing consists in applying a unitary V θ0 to the output noise unit together with an additional ancilla of the same size d as the system.The system and ancilla can be regarded as a single open system (denotes 's+a') of dimension D = d 2 C d < l a t e x i t s h a 1 _ b a s e 6 4 = " h 2 C i 2 1 P e g 6 b q L P O r / a d c i y 1 F W Figure 2: Adaptive measurement filter: the output units undergo post-processing with a coherent quantum absorber, followed by applying an adaptive measurement computed with the algorithm coupled to the noise units via the unitary W θ = V θ0 U θ .The unitary V θ0 is chosen such that s+a has a pure stationary state |ψ ∈ C D at θ = θ 0 , and the output state is identical to that of the input.This is a discrete-time analogue of the notion of coherent quantum absorber introduced in [55], and it insures that the 'reference' output state at θ = θ 0 is the same as the product input state (the 'vacuum'), while deviations from θ 0 produce 'excitations' in the output.After the interaction with the ancilla (absorber), the noise unit is measured in a basis determined by a simple iterative algorithm detailed in section 5.The iterative step consists of using the current value of a certain 'filter operator' on system+absorber to determine the next measurement basis, and then using the measurement outcome to update the filter operator.This simplification relies on the fact that the output is uncorrelated from system (and absorber), which is not the case in the original dynamical setup of Figure 1.
In section 7 we describe the results of two numerical investigations testing our theoretical results.The first investigation focuses on a simplified model where the system plus absorber are represented by a two-dimensional system with a pure stationary state.While this sidesteps the preliminary estimation stage of the protocol, it allows us to specifically test the key features of the adaptive measurement algorithm with a reasonably large trajectory length and a high number of repetitions.For this model, we can explicitly compute the system-output QFI (cf.Lemma 7.1), while the classical Fisher information of any output measurement strategy can be estimated by sampling techniques.The results confirm that the adaptive measurement attains the QFI when the system is measured at the end, while the output-only strategy is only worse by a constant independent of trajectory length.On the other hand, simple measurements (same fixed basis for each unit) perform strictly worse even when the measurement basis is optimised.While the improvement here is not dramatic, our preliminary investigations indicate that the gap increases significantly with the system dimension, depending on the chosen model.We further test the performance of the maximum likelihood estimator and find that its mean square error approaches the inverse of the classical Fisher information in the long-time limit, which agrees broadly with the Cramér-Rao bound.The second numerical investigation implements the full two-stage adaptive measurement algorithm including the use of the coherent absorber.
Finally, we note that our scheme can be extended to continuous-time dynamics by using standard time-discretisation techniques [56,57].Although we do not treat this in detail here, we comment on this extension at the end of the paper.
The paper is organised as follows.In section 2 we briefly review the adaptive algorithm for optimal separable measurements developed in [53].Section 3 introduces the Markov dynamics setting and reviews a key result on the asymptotic QFI of the output.Section 4 explains how the use of 'post-processing' by quantum absorber reduces the general estimation problem to one concerning a system with a pure stationary state.This is then used in section 5, which details the adaptive measurement procedure including the key 'measurement filter' algorithm.In section 6.2 we show that the proposed adaptive output measurement achieves the optimal QFI rate even if the system is not measured.We also devise a scheme to estimate the classical Fisher information of the measurement process by sampling over trajectories.Section 7 presents simulation results using an elaboration of an amplitude decay qubit model.

Optimal separable measurements
In this section we review a result by Zhou, Zou and Jiang [53] concerning optimal parameter estimation for multipartite pure states, using separable measurements (local measurements and classical communication).Their method will then be applied to the problem of estimating parameters of discrete time quantum input-output systems.By exploiting the Markovian nature of the dynamics, we will show that the algorithm can be recast in a simpler procedure akin to that of a quantum state filter.
Consider a one parameter quantum statistical model {ρ θ : θ ∈ R} where ρ θ is a state on a Hilbert space H which depends smoothly on the unknown parameter θ.To estimate θ we perform a measurement on the state ρ θ and compute an estimator θ based on the measurement outcome.According to the quantum Cramér-Rao bound (QCRB) [49][50][51], the variance of any unbiased estimator θ is lower bounded as where F θ is the quantum Fisher information (QFI) defined as F (θ) = Tr(ρ θ L 2 θ ), and L θ is the symmetric logarithmic derivative (SLD) satisfying d dθ ρ θ = L θ • ρ θ .In general, for any given parameter value θ 0 , the QCRB is saturated1 by measuring the SLD L θ0 and constructing a locally unbiased estimator θ = θ 0 + X/F (θ 0 ) where X is the measurement outcome.
While for full rank states the optimal measurement is essentially unique, for rank deficient states this is not the case and a necessary and sufficient condition for a measurement to saturate the QCRB has been derived in [51].This has practical relevance for multipartite systems where the measurement of the SLD may not be easy to implement.Motivated by this limitation, the saturability condition has been further investigated in [53] where it is shown that the QCRB for pure states of multipartite systems is achievable using separate measurements constructed in an adaptive fashion which we now proceed to describe.
Consider the pure state model ρ θ = |ψ θ ψ θ | with |ψ θ ∈ H.We denote | ψθ = d dθ |ψ θ and assume that ψ θ | ψθ = 0.This can generally be arranged by choosing the (unphysical) phase of |ψ θ to have an appropriate dependence on θ.In particular, in this case we have |ψ Further, we define the operator M which will play a key role in the analysis The authors of [53] note that if a projective rank-one measurement {E i = |e i e i |} satisfies the conditions e i |M |e i = 0, and then it fulfils the general criteria of [51] and therefore it achieves the QCRB.In fact, the second condition can be relaxed to p(i) = 0 for all i, but we will stick to the chosen expression for concreteness.The achievability can be understood as follows.The conditions e i |M |e i = 0 implies that e i |ψ θ ψθ |e i is real, so that the phase of the basis vectors |e i can be chosen such that both e i |ψ θ and e i | ψθ are real for all i.Together with the condition p θ (i) = 0, this means that in the first order of approximation, the quantum model is described by vectors with real coefficients with respect to the measurement basis.In this case the classical and quantum informations coincide We now assume that we deal with a multipartite system such that and follow [53] to show that a separable measurement satisfying the above conditions can be constructed by using the algorithm outlined below.For any set of indices A ⊂ {1, . . ., n} we denote its complement by A c and by ρ A = Tr A c (ρ θ ) the partial state of the sub-systems with indices in A. For m > 1 we denote by m the set {1, . . ., m}.Similarly, we denote M A = Tr A c (M ) and for single sub-systems (A = {i}) we use the notation ρ i and M i .
We measure the sub-systems sequentially, such that each individual measurement basis depends on the outcomes of the previous measurements, as follows.In the first step, the measurement basis {|e The existence of such a basis can be established by induction with respect to dimension, cf.proof of Lemma 1 in [53].The concrete construction in two dimensions is described in section 7.
After this, the following procedure is applied sequentially to determine the measurement basis for system j + 1 with j = 1, . . ., n − 1: given the outcomes i j := {i 1 , . . ., i j } of the first j measurements, we choose the basis |e , and e Note that the second condition means that for each j the outcome i j is independent of the others and has equal probabilities 1/k j .
After n steps we have defined in an adpative fashion a product measurement basis and one can verify that such a measurement satisfies the general condition (3).
3 Discrete quantum Markov chains and the output QFI In the input-output formalism the dynamics of a discrete-time quantum open system H s ∼ = C d is modeled by successive unitary interactions with independent 'noise units', identically prepared in a state |χ ∈ H u ∼ = C k .This can be pictured as a conveyor belt where the incoming 'noise units' constitute the input, while the outgoing 'noise units' make up the output of the process, cf. Figure 1.If |φ ∈ H s is the initial state of the system, and U is the unitary on H s ⊗ H u describing the interaction between system and a noise unit, then the state of the system and output after n time units is where U (i) is the unitary acting on the system and the i-th noise unit.From equation (4) we find that the reduced state of the system at time n is given by where the partial trace is taken over the output (noise units), and T : Fixing an orthonormal basis {|1 , . . ., |k } in H u , we can express the system-output state as a matrix product state where K i = i|U |χ are the Kraus operators of T , so that T (ρ) = i K i ρK * i .Now, let us assume that the dynamics depends on a parameter θ ∈ R which we would like to estimate, so that U = U θ and |Ψ(n) = |Ψ θ (n) .In the input-output formalism it is usually assumed that the experimenter can measure the output (noise units after the interaction) but may not have access to the system.In this case the relevant quantum statistical model is that of the (mixed) output state given by The problem of estimating θ in this formulation has been investigate in both the discrete time [43,47] and the continuous time [44,46,48] settings.For our purposes, we summarise here the relevant results of [47].We will assume that the Markov chain is primitive, i.e. the transition operator T has a unique full-rank steady state ρ ss (i.e.T (ρ ss ) = ρ ss ) , and is aperiodic (i.e. the only eigenvalue of T on the unit circle is 1).In particular, for any initial state ρ in , the system converges to the stationary state T n (ρ in ) → ρ ss in the large n limit.Therefore, for the asymptotic analysis we can assume that the dynamics is in the stationary regime and focus on the large time properties of the stationary output state.The following Theorem shows that the output QFI scales linearly with time and provides an explicit expression of the rate.Theorem 3.1.Consider a primitive discrete time Markov chain as described above, whose unitary depends smoothly on a one-dimensional parameter θ, so that U = U θ .The quantum Fisher information F θ (n) of the output state ρ out θ (n) scales linearly with n and its rate is equal to where R is the Moore-Penrose inverse of Id − T θ .
Following standard quantum Cramér-Rao theory [49][50][51], the theorem implies that the variance of any (unbiased) output-based estimator is bounded from below by n −1 /f (θ) for large n.A more in depth analysis [47] shows that the output model satisfies the property of local asymptotic normality which pertains to a certain quantum Gaussian approximation of the output state and implies that there exists an estimator θn which achieves the CR bound asymptotically and has normally distributed errors: where the convergence is in distribution to a normal variable with variance f (θ) −1 .Below, we will make use of an extension of Theorem 3.1 which shows that the same result holds for rank-deficient stationary states of ergodic chains, and in particular for pure states [59].
Having identified the output QFI rate, we would like to investigate measurement schemes which can provide good accuracy for estimating the parameter θ.As noted before, the QCRB can be achieved by measuring the SLD of the statistical model.However, the SLD of the output state is generally a complicated operator whose measurement requires collective operations on the noise units.On the other hand, one can consider separate measurements of the same observable on the different noise units and system.The average statistic may provide an efficient estimator and its (asymptotic) Fisher information can be computed explicitly [44].However, such measurements are in general not optimal.Here would would like to ask the more fundamental question: is it possible to achieve the QCRB using simpler 'local' manipulation of the output units which involve operations on single, rather than multiple units.

Output post-processing using quantum coherent absorber
In this section we introduce a key tool which will allow us to recast the estimation problem for a general primitive Markov chain ( which typically has a mixed stationary state) into one concerning a Markov chain with a 'doubled-up' system having a pure stationary state.The construction is a discrete-time adaptation on the concept of coherent quantum absorber introduced in [55] for continuous-time dynamics.In section 5 we then show how the absorber can be used to compute an adaptive, separable measurement in a simple recursive algorithm.
Consider the input-output system of section 3 characterised by a unitary U on H s ⊗ H u .We now modify the setup as illustrated in Figure 2 by inserting an additional physical d-dimensional system H a called absorber which interacts with each of the noise units via a fixed unitary V on H a ⊗ H u , applied immediately after U .This can be seen as a type of quantum post-processing of the output prior to the measurement.The original system and the absorber can be considered a single open system with space H s ⊗ H a which interacts with the same conveyor belt of noise units via the unitary W θ = V • U where V, U are now understood as the ampliations of the unitaries to the tensor product H s ⊗ H a ⊗ H u .The following lemma shows that for certain choices of V the auxiliary system forms a pure stationary state together with the original one, and the noise units pass unperturbed from input to output.This explains the 'absorber' terminology, which was originally introduced in the context of continuous-time input-output dynamics [55].Lemma 4.1.Any given primitive Markov with unitary U can be extended to a quantum Markov chain including an absorber with unitary V , such that the doubled-up system has a pure stationary state | ψ ∈ H s ⊗ H a and In particular, if the initial state of the doubled-up system is | ψ , then the n-steps output state of the doubled-up system is identical to the input state |χ ⊗n .
Proof.Let ρ ss = i λ i |f i f i | be the spectral decomposition of the stationary state of the original system with unitary U .We construct the purification which will play the role of stationary state of the extended system.Let |φ := U | ψ ⊗ χ ∈ H s ⊗ H a ⊗ H u be the state after applying U .We therefore look for unitary V on H a ⊗ H u (ampliated by identity on H s ) such that V reverts the action of U Since Tr a (| ψ ψ|) = ρ ss this means that the reduced state of the system after applying U is still the stationary state, so that where |g i are mutually orthogonal unit vectors in H a ⊗ H u .We now choose a unitary V such that V |g i = |f i ⊗ χ , for all i, which is always possible due to orthogonality.

Adaptive measurement algorithm
In this section we describe our adaptive output measurement protocol for estimating an unknown one-dimensional dynamical parameter θ of a discrete time quantum Markov chain with unitary U θ , as described in section 3. The protocol has two stages.In the first stage we use a small proportion (e.g.ñ = n 1− , with 0 < 1) of the output units in order to compute a preliminary 'rough estimator' θ 0 of the true parameter θ by performing a standard sequential measurement.This step is necessary in any quantum estimation problem in which the optimal measurement depends on the unknown parameter [54,58], and will inform the second stage of the protocol.In the second stage we use θ 0 to design a optimal sequential measurement for θ = θ 0 , i.e. one that achieves the output QFI at θ = θ 0 .Since δθ = θ − θ 0 = O(n −1/2+ ) [44], this implies that the procedure is asymptotically optimal for any parameter value θ, in that the classical Fisher information has the same linear scaling as the QFI F θ (n).This stage has two key ingredients (cf. Figure 2): a quantum 'post-processing' operation on output units immediately after interacting with the system, followed by an adaptive projective measurement whose basis is computed according to a 'measurement filter' algorithm inspired by [53].We now describe the two steps in detail.
Quantum post-processing.After the interaction U θ with the system, the output units interact sequentially with a d-dimensional ancillary system, cf. Figure 2. The interaction unitary V θ0 is chosen such that the ancillary system is a coherent quantum absorber for θ = θ 0 , see section 4 and Lemma 4.1 for construction.This means that the system plus absorber (s+a) can be regarded as a single D = d 2 dimensional open system with associated unitary W θ = V θ0 U θ , which has a pure stationary state at θ = θ 0 denoted |ψ ∈ C D , and whose output state is identical to the input.The general estimation problem for U θ has been reduced to a special one for a doubled-up system with unitary W θ which features a pure stationary state at θ = θ 0 .Remark 1.Since the absorber transformation V θ0 does not depend on θ and is applied after U θ , the overall effect over an n steps interval is to rotate the absorber plus output state by a fixed unitary θ0 .This means that the total QFI does not change by introducing the absorber.Adaptive measurement algorithm.We will assume for simplicity that the initial state of s + a is the stationary state |ψ such that the full (s+a)-output state is |Ψ θ (n) = W θ (n)|ψ⊗χ ⊗n as defined in (4).One could obtain similar results for different initial states by simply waiting long enought for the system and absorber to converge to the stationary state |ψ .In principle we could now apply the algorithm described in section 2 to construct and adaptive measurement whose classical Fisher information is equal to the system-output QFI.In fact we could have done this without using the absorber.However, this procedure has some drawbacks.Indeed, in order to compute the measurement bases one needs to work with large dimensional spaces which becomes unfeasible in an asymptotic setting.Secondly, it is not clear a priori whether the output units can be measured immediately after the interaction with the system, and whether the measurements depend on the length of the output (sample size).In addition, the procedure requires a final measurement on the system, which may be impractical in the context of input-output dynamics.We will show that all these issues can be addressed by taking into account the Markovian structure of our model, and exploiting the pure stationary state property.Let us denote W = W θ0 , Ẇ = dW dθ θ0 and and In Appendix 9 we show that the adaptive measurement of [53] reduces to the following iterative algorithm which is conceptually similar to quantum state filtering [4,5], and involves individual measurements on the output units immediately after the interaction with the system, and computations with operators on C D ⊗ C k at each step.

Initialisation step (j=1). The first measurement basis |e
[1] i in C k is chosen such that the following conditions are fulfilled: The first noise unit is measured in this basis and the outcome X 1 = i 1 is obtained.The filter at time j = 1 is defined as the (trace zero) s+a operator i1 A 1 e [1] i1 .
Iterative step.The following step is iterated for j = 2, . . ., n.Given the filter operator Π j−1 of the previous step, we define The j-th measurement basis |e [j] i is chosen to fulfill the conditions e [j] i B j e [j] i = 0, and e [j] We measure the j-th noise unit in the basis |e [j] i and obtain the result X j = i j .The filter at time j is updated to Final s+a measurement.This is an optional step which involves a final joint measurement on system and absorber.The basis |e is determined by the following conditions The system and absorber is measured in this basis and the outcome X = i 0 is obtained.The output measurement record {i 1 , . . ., i n , i 0 } is collected and used for estimating the parameter θ.The likelihood function is given by For later use, we denote by p θ (i 1 , . . .i n ) the marginal distribution of the output measurement record only.

Fisher informations considerations
In this section we investigate the relationship between the classical Fisher information (CFI) of the output adaptive measurement process and the system-output QFI.We prove that both scale with the same rate and the latter may be larger than the former by at most a constant, independent of time.
We also provide an expression of the CFI of sequential (adaptive or standard) output measurements, which is amenable to estimation by sampling.This tool will be used to confirm the optimality of our adaptive algorithm in numerical simulations.

Achievability of the QFI with adaptive output measurements
The adaptive measurement scheme described in section 5 insures the CFI of the full measurement (output and s+a) is equal to the QFI of the full pure state model where the last equality follows from the fact that the absorber acts as an additional rotation which does not change the QFI, cf.Remark 1.However, in certain physical implementations the system may not be accessible for measurements, so the more interesting scenario is that in which only the output state is measured.In this case the CFI will generally be strictly smaller that the QFI, and the question is whether by measuring only the output we incur a significant loss of information.In proposition 6.1 we show that this is not the case: the difference between the QFI and the output CFI is bounded by a constant, so for large times the loss of information is negligible compared to both QFI and output CFI, which scale linearly with time.In section 7 we will illustrate the result on a specific model.Proposition 6.1.Consider the setup described in section 5, and let F (n) be the system-absoberoutput QFI, and I (o) (n) be the output CFI for the optimal adaptive measurement, at θ = θ 0 .Then F (n) − I (o) (n) < c for all n where c is a constant depending only on the model U θ .Consequently, where f = f θ0 is the QFI rate (6).

Computing the classical Fisher information of the output
The classical Fisher information of the output measurement process at θ 0 is where the sum runs over indices such that p θ (i 1 , . . ., i n ) > 0. In general where the successive upper bounds are output and system-output QFIs respectively.
In our simulation study we will be interested to study to what extent these bounds are saturated in the adaptive and non-adaptive scenarios, and in particular, to verify the prediction of Proposition 6.1.Since the classical Fisher information is difficult to compute for long trajectories, we will recast it as an expectation which can be estimated by sampling measurement trajectories.In Lemma 6.2 below, we will use the fact that at θ = θ 0 the vector |ψ is the stationary state, and therefore for any Kraus decomposition K j |W |χ (for simplicity we use the same notation for s+a Kraus operators as in section 3).The proof of Lemma 6.2 can be found in Appendix 11.Lemma 6.2.Consider the setup described in section 5.The output CFI at θ = θ 0 is given by where f is the function [n] in and the constants c We now consider the case where a (projective) measurement {P (s+a) i } is performed on s+a, after obtaining the output measurement record (i 1 , . . ., i n ).We denote the additional outcome by i 0 .The CFI of the full process is is the likelihood function of a trajectory augmented by the system measurement outcome i 0 .The relevant upper bound in this case is A similar computation to that of Lemma 6.2 gives the system-output classical Fisher information where f is the function For fixed (non-adaptive) measurements, the bound (12) is generally not saturated except for special models (e.g. if the state coefficients in the measurement basis are real for all θ).In contrast, the system-absorber-output classical Fisher information for adaptive measurements is equal to the QFI thanks to the optimality of the adaptive measurement procedure 9.This is confirmed by our simulation study which also investigates the performance of the fixed measurement scenario.

Numerical simulations
We now test the key properties of the adaptive measurement scheme developed in section 5, in two separate numerical investigations.
The first investigation described in subsections 7.1 and 7.2 employs a simplified Markov model which bypasses stage-one of the scheme (computing a rough estimator θ 0 ) and simulates data at θ = θ 0 .This allows us to directly study the performance of the algorithm itself (stage two), rather than that of the combination of the two stages.The second simplification of this investigation is that we choose a system which has a pure stationary state at θ 0 ; this means that no absorber is required, so the system can be seen as a surrogate for the system+absorber in the general scheme.The reason for this is mainly practical, as it allows us to use a two-dimensional system while system+absorber would have dimension at least four.
The second numerical investigation consists of a full simulation study including the use of the coherent absorber and the two stage estimation procedure, and its results are presented in subsection 7.3.Here we can see the overall performance of the estimation method, but it is harder to estimate the Fisher information of the measurement process and to separate the contribution of the two stages in the overall estimation error.

Simplified Markov model for the first numerical investigation
We consider a dynamical model consisting of a two-dimensional system coupled via a unitary U θ to two dimensional noise units in state |χ = |0 .The input state and the unitary are designed such that the stationary state at θ 0 = 0 is |ψ = |0 .Since the input is prepared in a fixed state, we only need to define the action of U θ on the basis vectors |0 ⊗ |0 and |1 ⊗ |0 .The following choice has unknown parameter θ and two known parameters λ and φ A non-zero value of the phase parameter φ ensures that the system-output state does not have real coefficients in the standard basis, in which case the standard basis measurement would be optimal.Note that at θ 0 the dynamics provides a simple model for 'photon-decay' with decay parameter λ.We compare two measurements scenarios.In the first, non-adaptive scenario, the noise units are measured in a fixed orthonormal basis 2} while in the second scenario they are measured adaptively following the algorithm described in section 5.In both cases, given the measurement record {i 1 , . . ., i n } ∈ {0, 1} n , the conditional state of the system is where the Kraus operators are i |U |χ in the second one.The likelihood of the output measurement trajectory is Recall that the optimal measurement basis needs to satisfy the conditions (7) in section 5.For two dimensional noise units the computation reduces to the following scheme.We express the traceless, anti-Hermitian matrix B j defined in section 5 as B j = i r [j] • σ where r [j] = (r x , r y , 0) is its Bloch vector, and similarly, we let ± s [j] be the Bloch vectors of the basis vectors {|e satisfying the conditions (7).Then one finds that the conditions (7) are satisfied if s [j] is taken to be s x , 0).We study both the scenario where system of interest is measured after obtaining the output trajectory, as well as the one where only the output measurement is considered.In the former case, the algorithm guarantees that the classical Fisher information of the full measurement record is equal to the QFI of the system-output state.This claim is verified numerically by comparing the estimated classical Fisher information computed using the method described in section 6.2 with the quantum Fisher information of the system-output state.In the latter scenario, Proposition 6.1 insures that the loss of information compared to a 'full measurement' is bounded by a constant which does not depend on time.
As shown in Theorem 3.1, the system-output QFI scales linearly with time, i.e.
where the QFI rate is given by the equation ( 6) and the o(n) term depends on the specific system parameters.Applying this to our model we obtain However, it turns out that for this specific model and parameter value, the QFI can be computed explicitly for any fixed n.The proofs of (16) and of the following lemma can be found in Appendix 12.
Lemma 7.1.The system-output QFI at θ = θ 0 is given by the formula where a = √ 1 − λ and b = √ λ.In particular, the leading term in n is given by (16) while the remaining terms are bounded.

Simulation studies for the simplified model
We now present the results of our first numerical investigation consisting of 3 simulation studies using the simplified Markov model described above.
The first simulation study focuses on the comparison between the different notions of Fisher information: the system-output QFI, the CFI of the output trajectory in the non-adapted and adapted scenarios, and the CFI of the system-output measurement process in the adapted measurement scenario.The QFI is computed using the formula in Lemma 7.1 while the CFIs are estimated by sampling using the expression in Lemma 6.2.
The results are illustrated in Figure 3 where the different informations are plotted as a function of time (trajectory length) n, for λ = 0.8, φ = π/4.The simulation confirms the fact that the adaptive algorithm achieves the QFI when the system is measured together with the output, while the CFI of the output trajectory (without measuring the system) provides a close approximation which only differs by a constant factor.In contrast, the CFI of the standard measurement has a smaller rate of increase; additional numerical work shows that the CFI rate can be improved by optimising the basis of the standard measurement but it does not achieve the QFI.
The second simulation study focuses on the 'measurement trajectory', i.e. the sequence of measurement settings produced in the adaptive measurement scenario.Since all measurement bases consist of vectors in the equatorial plane on the Bloch sphere, one can parametrise each basis by the polar coordinate ϕ of the basis vector (|0 ± e iϕ |1 )/ √ 2 for which ϕ belongs to a specified interval of length π.This is illustrated in the left panel of Figure 4.This parametrisation has the disadvantage that it does not reproduce the topology of the space of measurements which is that of a circle, leading to some jumps in measurement angles appearing to be larger than the actual 'distance' between measurement bases.To remedy this, in the right panel of Figure 4 we plot 2ϕ on circles of radius increasing linearly with time.We note that the initial steps of the trajectory  show large variations which then 'stabilise' around a certain range of values.This has to do with the fact that the initial angle can be chosen arbitrarily in this model as B 1 = 0, cf.section 5. Understanding the nature of this stochastic process remains an interesting topic of future research.
The third simulation study concerns the performance of the maximum likelihood estimator (MLE) in an adapted and non-adapted output measurement scenarios.The MLE is defined by where the likelihood p τ (i 1 , . . .i n ) is computed as in equation (15).In numerics, we maximise the log-likelihood function which can be computed as the sum where ψ(i 1 , . . .i j−1 ) is the system's state conditional on the output trajectory (filter), cf.equation (14).The MLE accuracy is quantified by the mean square error (MSE) E(n) = E θ ( θn − θ) 2 , which is estimated empirically by averaging over a number of simulations runs to obtain Ê(n).
To verify that the MSE scales as n −1 we plot the inverse empirical error Ê−1 (n) as a function of time.Figure 5 shows the inverse error for the adaptive and non-adaptive measurements together with the QFI (which is equal to the CFI of the adaptive measurement) and the CFI of the nonadaptive measurement.We note that for small values of n the inverse error is significantly lower  that the corresponding Fisher information, but it approaches the latter for larger values of n.This suggests that the MLE achieves the Cramér-Rao bound asymptotically, which is not surprising since the MLE is known to be asymptotically optimal for independent samples as well as for certain classes of hidden Markov chains [60,61].However, proving its optimality for the adaptive measurement process remains an open problem.In addition to the mean square error, we looked at the distribution of the MLE. Figure 6 shows a histogram of the MLE based on N = 10000 simulations with n = 200, and it indicates that the MLE is approximately normally distributed.Based on this, it is reasonable to conjecture that the MLE is an efficient estimator [60,61] (i.e. has asymptotically normal distribution with variance equal to the inverse of the Fisher information).

The second numerical investigation using the full adaptive protocol
In this section we present the results of the second numerical investigation which implements the full estimation scheme proposed in this paper, including the preliminary estimation stage and the use of the coherent absorber.We use the same input-output model as described by the unitary in equation (13), but the data will be simulated at a true value of θ = 0.2 instead of θ = 0.While at θ = 0 the system has a pure stationary state and the coherent absorber is not needed, away from this value the stationary state is mixed and we will apply the full protocol described in section 4.
In the first stage of the adaptive estimation procedure we use a fixed proportion q of the total sample size n to obtain a preliminary estimator θ 0 of θ by measuring each output unit in the standard basis.The parameter is estimated using the maximum likelihood method.In the second stage we apply the adaptive scheme with an absorber 'tuned' to the parameter value θ 0 (cf.section 4).The system and absorber are prepared in the pure stationary state at θ 0 , the dynamics is run for time (1 − q)n and the output is measured according to the adaptive measurement filter algorithm.The maximum likelihood estimator for the stage one and two data is then computed.
For comparison, we also run a non-adaptive scheme where the output is measured in the standard basis for the whole duration n, without using an absorber.The maximum likelihood estimator is again computed from the measurement data.In both experiments, the mean square error of the MLE is estimated by averaging over 1000 repetitions.
Unlike the setup of the previous numerical study, the estimation of the classical Fisher information of the adaptive measurement process was too costly and is not included in the study.As a proxy for the classical Fisher information we plot the average value of the observed Fisher information for each of the two simulations.For a given measurement run (i 1 , . . .i n ), the observed Fisher information is defined as the second derivative of the log-likelihood function evaluated at the maximum likelihood estimator θn : While a full theoretical justification of I obs goes beyond the scope of this paper, we note that the use of the observed Fisher information for independent identically distributed data is well grounded in statistical methodology and is closely relate to the asymptotic normality property [62].
Figure 7 shows the results of the numerical experiments, for two values q = 0.15 and q = 0.25 of the proportion of samples used in the preliminary estimation stage.As before, we plot the inverse of the estimated mean square error for the simple measurement (green line) and the adaptive measurement (blue line).In addition we plot the leading contribution to the quantum Fisher information nf (θ) (purple line) computed using the results in Theorem 3.1, which provides the asymptotic slope of the actual output quantum Fisher information.The average observed Fisher information is plotted for both the simple (red line) and adaptive (orange line) measurement setups.We note a good agreement between the inverse mean square error of the MLEs and the observed Fisher informations.For the adaptive measurements, the observed Fisher information also shows the same slope as the asymptotic Fisher information, as expected.We also note that in the adaptive measurement, the mean square error does not quite achieve the (observed) quantum Fisher information for the range of times we considered, although it gets closer to in the case q = 0.25 (right panel).We speculate that this may be related to a number of factors such as the choice of q for different sample sizes, the details of the actual measurement procedure implemented in practice (in contrast to the theoretical prescription) and even the implementation of the MLE.For instance, in practice it may be beneficial to implement several adaptive estimation stages where the preliminary estimator is gradually improved and used in tuning the absorber in the next stage.All these remain interesting questions which are worth investigating in more detail and on a case by case basis.However, these are somewhat separate issues from that of designing adaptive measurements that achieve the QFI, which was the main focus of this work.

Conclusions and Outlook
In this paper we developed an efficient iterative algorithm for optimal estimation of dynamical parameters of a discrete-time quantum Markov chain, using adaptive sequential measurements on the output.The algorithm builds on the general measurement scheme of [53] which achieves the quantum Fisher information for pure state models of multipartite systems with one dimensional unknown parameters.However, unlike the scheme of [53] which requires manipulations involving the full multipartite state, the proposed algorithm only involves computations on d 2 •k-dimensional systems where d and k are the dimensions of the open system and noise unit respectively.Therefore, the method can be readily applied to Markov parameter estimation for large output sizes.The algorithm exploits the Markovian structure of the dynamics to sequentially compute optimal measurement bases in terms of a single time-dependent 'measurement filter' operator, which is updated in a way that is reminiscent of a state filter.One of the key ingredients of the proposed scheme is the use of a coherent quantum absorber [55] which reduces the estimation problem to one concerning a system with a pure stationary state.We considered both output and outputsystem measurements scenarios and we showed that while the former achieve the full quantum Fisher information, the latter is short of this by just a fixed constant, and in particular both have the same scaling with time.Our theoretical results are confirmed by numerical simulations using a simplified model related to the amplitude decay channel.We also presented results from 'full simulation' study involving the use of the coherent absorber.
Our discrete-time procedure raises interesting questions about the the possibility to design realistic optimal sequential measurements in continuous-time dynamics.In principle the scheme can be applied to continuous time by using time-discretisation techniques [56,57].Although we did not treat this in detail here, we can readily make several observations in the special case of a single Bosonic I-O channel.For small enough time intervals δt, the field can be approximated by a noise unit C 2 with basis vectors representing the vacuum and a one-photon state respectively.The proposed measurements consist of projections whose corresponding Bloch vectors are in the equatorial plane; loosely speaking, this corresponds to an adaptive homodyne measurement with a time-dependent angle.Our preliminary investigations indicate that the behavior of the timedependent angle ranges from deterministic evolution to a noisy stochastic process which does not appear to be of diffusive type.This raises the question whether the remaining freedom in choosing the measurement bases can be used to improve the adaptive algorithm and produce a more regular measurement processes.Indeed, the second defining condition for the measurement vectors can be relaxed, allowing for more general classes of optimal measurements, which may be more suitable for continuous-time direction.This will be the topic of a future investigation.
Another important open question concerns the extension to chains with mixed input states, or more multiple inputs, of which only some are observed.We speculate that for small departures from the current scheme, the algorithm will be quasi-optimal for some time interval but will be sub-optimal in the long time limit.In this case, restarting the evolution and measurement filter at regular intervals may be more efficient.
From a theoretical viewpoint, it is important to understand the mathematical properties of the stochastic processes introduced here, the adaptive measurement and the measurement trajectory.Finally, it is intriguing to consider to what extent the proposed method can adapted to multi-parameter estimation and to general (time-dependent) matrix product states, as opposed to stationary output states of Markov processes.

Appendix: Derivation of the measurement filter
We start by applying the algorithm [53] to our estimation problem.The natural setup is to measure the noise units in the order in which they emerge in the output, followed by a final measurement on the system+absorber.For simplicity we refer to the latter as the system.We start by defining where the index n keeps track of the output length.In this section we will use the label 0 for the system, and 1, . . ., n for the noise units of the output.The algorithm prescribes measurement bases |e for each of the output units l = 1, . . ., n in an adaptive, sequential fashion.The first basis satisfies the equations e where M (n) 1 = Tr 0,2,...,n M (n) .The second equality follows from the fact that at θ 0 we have |Ψ(n) = |ψ ⊗ χ ⊗n .The next basis depends on the outcome i 1 of the first measurement and satisfies the constraints .
The measurement basis e The last step consists of measuring the system using the same procedure as for the output units.A priori, the procedure depends on the size n of the output.The following lemma shows that the optimal bases obtained for different output sizes coincide.Lemma 9.1.Let j, n be two output lengths with j < n and consider applying the above procedure to the corresponding states |Ψ (j) and respectively |Ψ (n) .The optimal measurement bases {|e [l] i } for the units l = 1, . . ., j satisfy the same constraints and can be chosen to be the same.In addition we have M (n) j (i j−1 ) = M (j) j (i j−1 ).Proof.Recall that |Ψ(n) = W (n) . . .W (1) |ψ ⊗ χ ⊗n , and let us denote W (n) . . .Ẇ (i) . . .W (1) |ψ ⊗ χ ⊗n = n i=1 W (n) . . .Ẇ (i) |ψ ⊗ χ ⊗n where we used the fact that W leaves |ψ ⊗ χ invariant.
We first show that the matrix M (n) 1 does not depend on n, and therefore the first measurement basis does not depend on the length of the output.

Let e
[1] i be the measurement basis determined from 1 .We now show that, conditional on the outcome i 1 of this measurement, the second basis e [2] i does not depend on the length of the output, which can be take to be equal to 2.     where the first equality is due to measurement optimality, while the second is the inequality between classical and quantum information.This implies that Let | Ψ(1) := Ẇ |ψ ⊗ χ and let τ := Tr 1 ( Ψ(1) Ψ(1)|).Then (P ⊥ ψ ⊗ I)W (n) . . .W (j+1) Ẇ (j) ψ ⊗ χ ⊗n 2 = Tr 0,1,...n−j (P ⊥ ψ ⊗ I)W (n−j) . . .W (1) τ W (1) * . . .W (n−j) * = Tr 0 (P ⊥ ψ T n−j (τ )) where in the last equality we have uses the definition of the transition operator T of the sys-tem+absorber.Assuming that T is ergodic we have T n (τ ) → P ψ exponentially fast with n so that Tr 0 (P ⊥ ψ T n−j (τ )) ≤ a 2(n−j) for some a < 1. Therefore Note that if the spectral gap of T becomes small then the convergence to stationarity is slower and the upper bound increases.
11 Appendix: Proof of Lemma 6.2 The CFI of any output measurement can be computed explicitly by writing where f is the function where K i are the (fixed) Kraus operators with respect to the standard basis, and |psi = |0 .Our specific model ha the feature that both K i and Ki map the basis vectors into each other: This allows to compute the QFI explicitly by noting the terms in the sum (25) with more that two indices equal to 1 have zero contribution.The remaining terms can be computed as follows.The term with only zero indices is where a = √ 1 − λ.Consider now a sequence (0, . . ., , 0, 1, 0, . . .0) with a single one on position l.Since K 0 . . .K0 . . .K 1 . . .K 0 |0 = 0 the only contributing terms will be those with derivative on the first (l − 1) K 0 s or on K 1 .This gives where b = √ λ.Finally, consider the sequences of the type (0, . . .0, 1, 0, . . .0, 1, 0 . . .0) with 1s on positions 1 ≤ i < k ≤ n.In this case the nonzero contributions come from terms where the derivative is on positions j = i.The Fisher contribution is Adding together the contributions (26), ( 27) and (28) we obtain the total QFI where the leading term is consistent with the QFI rate formula (16).
C k< l a t e x i t s h a 1 _ b a s e 6 4 = " / k l e C e G 0t d R z A r U G / R q r A e z H J d w = " > A A A B 9 X i c b V D L T g I x F L 2 D L 8 Q X 6 t J N I 5 i 4 I j M E o 0 s i G 5 e Y y C O B g X R K B x o 6 n U n b 0 Z A J / + H G h c a 4 9 V / c + T d 2 Y B Y K n q T J y T n 3 5 p 4 e L + J M a d v + t n I b m 1 v b O / n d w t 7 + w e F R 8 f i k r c J Y E t o i I Q 9 l 1 8 O K c i Z o S z P N a T e S F A c e p x 1 v 2 k j 9 z i O V i o X i Q c 8 i 6 g Z 4 L J j P C N Z G G pT 7 A d Y T z 0 s a 8 8 G 0 P C y W 7 I q 9 H M 1 Z 2 c 4 p / I H 1 + Q P p X J I i < / l a t e x i t >C k < l a t e x i t s h a 1 _ b a s e 6 4 = " / k l e C e G 0 t d R z A r U G / R q r A e z H J d w = " > A A A B 9 X i c b V D L T g I x F L 2 D L 8 Q X 6 t J N I 5 i 4 I j M E o 0 s i G 5 e Y y C O B g X R K B x o 6 n U n b 0 Z A J / + H G h c a 4 9 V / c + T d 2 Y B Y K n q T J y T n 3 5 p 4 e L + J M a d v + t n I b m 1 v b O / n d w t 7 + w e F R 8 f i k r c J Y E t o i I Q 9 l 1 8 O K c i Z o S z P N a T e S F A c e p x 1 v 2 k j 9 z i O V i o X i Q c 8 i6 g Z 4 L J j P C N Z G G p T 7 A d Y T z 0 s a 8 8 G 0 P C y W 7 I q 9 H M 1 Z 2 c 4 p / I H 1 + Q P p X J I i < / l a t e x i t >C k < l a t e x i t s h a 1 _ b a s e 6 4 = " / k l e C e G 0 t d R z A r U G / R q r A e z H J d w = " > A A A B 9 X i c b V D L T g I x F L 2 D L 8 Q X 6 t J N I 5 i 4 I j M E o 0 s i G 5 e Y y C O B g X R K B x o 6 n U n b 0 Z A J / + H G h c a 4 9 V / c + T d 2 Y B Y K n q T J y T n 3 5 p 4 e L + J M a d v + t n I b m 1 v b O / n d w t 7 + w e F R 8 f i k r c J Y E t o i I Q 9 l 1 8 O K c i Z o S z P N a T e S F A c e p x 1 v 2 k j 9 z i O V i o X i Q c 8 i6 g Z 4 L J j P C N Z G G p T 7 A d Y T z 0 s a 8 8 G 0 P C y W 7 I q 9

Figure 1 :
Figure 1: Quantum input-output discrete-time dynamics with θ dependent unitary interaction U θ 2 y S w 6 I S 4 5 J k 5 y T C 9 I i j N y T R / J M X q w H 6 8 l 6 t d 4 m r T P W d G a b / I D 1 / g U W R 5 s V < / l a t e x i t > |ψ θ ψ θ |) | ψθ = | ψθ .Under this assumption the QFI is given by ij are defined by equation 10.In particular, I (o) (n) can be estimated by computing the empirical average of f 2 over sampled trajectories.

Figure 3 :
Figure 3: Fisher informations as function of output lenght n: quantum Fisher information (blue), classical Fisher information for the adaptive measurement with/without system measurement (orange/green), classical Fisher information for a regular (non-adaptive) measurement.

Figure 4 :
Figure 4: Adaptive measurement trajectory with basis angle ϕ on left and ϕ plotted on the circle on the right.

Figure 5 :
Figure 5: Comparison of different Fisher informations of the output state as function of trajectory length n, obtained by averaging over N = 10000 trajectories, at λ = 0.8, ϕ = π/4.The quantum Fisher information (QFI) (blue) and the classical Fisher information (CFI) of adapted measurement (green) completely overlap due to optimality.The inverse MLE error of the adapted measurement (orange) approaches the QFI for large n.Similarly, for the standard measurement, the inverse MLE error (red) approaches the CFI (purple) for large n.

Figure 7 :
Figure 7: Comparison of different Fisher informations of the output state as function of trajectory length n, in the fully adaptive numerical investigation, for two proportions of samples used in the preliminary stage: q = 0.15 (left panel) and q = 0.25 (right panel).The results are obtained by averaging over N = 1000 trajectories, at θ = 0.2 and λ = 0.8, ϕ = π/4.We plot: the asymptotic QFI (purple line), observed Fisher information for adaptive measurement (orange line) and simple measurement (red line), inverse MSE of the MLE for adaptive measurement (blue line) and simple measurement (green line).