Machine learning non-Markovian quantum dynamics

TL;DR

This paper introduces a machine learning approach to extract environmental information from quantum measurement data, enabling better control of non-Markovian quantum systems without full process tomography.

Contribution

It develops a novel method to infer non-Markovian environment characteristics from measurement sequences using Markovian embedding and maximum likelihood estimation.

Findings

Successfully predicts non-Markovian dynamics

Enables environment characterization without process tomography

Improves quantum system control and manipulation

Abstract

Machine learning methods have proved to be useful for the recognition of patterns in statistical data. The measurement outcomes are intrinsically random in quantum physics, however, they do have a pattern when the measurements are performed successively on an open quantum system. This pattern is due to the system-environment interaction and contains information about the relaxation rates as well as non-Markovian memory effects. Here we develop a method to extract the information about the unknown environment from a series of projective single-shot measurements on the system (without resorting to the process tomography). The method is based on embedding the non-Markovian system dynamics into a Markovian dynamics of the system and the effective reservoir of finite dimension. The generator of Markovian embedding is learned by the maximum likelihood estimation. We verify the method by comparing its prediction with an exactly solvable non-Markovian dynamics. The developed algorithm to learn unknown quantum environments enables one to efficiently control and manipulate quantum systems.

Equations

$${\varrho }_S(t)$$

$$\frac{d{\varrho }_{S+ER}(t)}{dt}$$

$$p={\rm tr}\left\{E_{n}\ldots\Phi\big{[}E_{1}\Phi[\varrho_{S+ER}(0)]E_{1}\big{]}\ldots E_{n}\right\}.$$

$$p={\rm tr}\big{\{}\Phi[\widetilde{\varrho}_{S+ER}(t_{m-1})]E_{m}{\cal E}_{S+ER}(t_{m})E_{m}\big{\}},$$

$$\Phi [{\varrho }_{S+ER}]={\mathrm{tr}}_A\left[U(H){\varrho }_{S+ER}\otimes {\varrho }_A{U}^{†}(H)\right],$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!g_{\mu\nu}=\frac{1}{p}\sum_{m=1}^{n}{\rm tr}\bigg{\{}[E_{m}\,{\cal E}_{S+ER}(t_{m})\,E_{m}]\otimes I_{A}$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\times\bigg{[}\frac{\partial U(H)}{\partial H_{\mu\nu}}\,\,\widetilde{\varrho}_{S+ER}(t_{m-1})\otimes\varrho_{A}\,\,U^{{\dagger}}(H)+{\rm H.c.}\bigg{]}\bigg{\}},\quad$$

$$\frac{\partial U(H)}{\partial H_{\mu \nu }}=\sum _{k,l}\frac{e^{-i{\lambda }_k\tau }-e^{-i{\lambda }_l\tau }}{{\lambda }_k-{\lambda }_l}⟨{\psi }_k|\mu ⟩⟨\nu |{\psi }_l⟩|{\psi }_k⟩⟨{\psi }_l|.$$

$$\widetilde{\varrho}_{S+ER}(t_{i+1})=E_{i}\Phi[\widetilde{\varrho}_{S+ER}(t_{i})]E_{i}={\rm tr}_{A}\big{[}E_{i}\otimes I_{A}\,U(H)\,\widetilde{\varrho}_{S+ER}(t_{i})\otimes\varrho_{A}U^{{\dagger}}(H)\,E_{i}\otimes I_{A}\big{]}$$

$${\cal E}_{S+ER}(t_{i-1})=\Phi^{{\dagger}}[E_{i}{\cal E}_{S+ER}(t_{i})E_{i}]={\rm tr}_{A}\big{[}U^{{\dagger}}(H)\,(E_{i}{\cal E}_{S+ER}(t_{i})E_{i}\otimes I_{A})U(H)\,(I_{S+ER}\otimes\varrho_{A})\big{]}$$

$$p(\{E_i{\}}_{i=1}^n|H)$$

$$\frac{\partial p(\{E_{i}\}_{i=1}^{n}|H)}{\partial H_{\mu\nu}}=\sum_{m=1}^{n}{\rm tr}\bigg{\{}E_{m}{\cal E}_{S+ER}(t_{m})E_{m}\otimes I_{A}\bigg{[}\frac{\partial U(H)}{\partial H_{\mu\nu}}\,\widetilde{\varrho}_{S+ER}(t_{m-1})\otimes\varrho_{A}U^{{\dagger}}(H)+U(H)\widetilde{\varrho}_{S+ER}(t_{m-1})\otimes\varrho_{A}\frac{\partial U^{{\dagger}}(H)}{\partial H_{\mu\nu}}\bigg{]}\bigg{\}}.$$

$$\frac{\partial U(H)}{\partial H_{\mu \nu }}={\left(\frac{\partial U^{†}(H)}{\partial H_{\mu \nu }}\right)}^{†}=\sum _{k,l}\frac{e^{-i{\lambda }_k\tau }-e^{-i{\lambda }_l\tau }}{{\lambda }_k-{\lambda }_l}⟨{\psi }_k|\mu ⟩⟨\nu |{\psi }_l⟩|{\psi }_k⟩⟨{\psi }_l|.$$

$$\frac{\partial \log p(\{E_i{\}}_{i=1}^n|H)}{\partial H_{\mu \nu }}=\frac{1}{p(\{E_i{\}}_{i=1}^n|H)}\frac{\partial p(\{E_i{\}}_{i=1}^n|H)}{\partial H_{\mu \nu }}.$$

$${\varrho }_E=P_1\otimes P_2\otimes \dots \otimes P_n$$

$$W=\exp (-iH\tau )=\mathsf{SHIFT}\cdot \mathsf{SWAP},$$

$$\mathsf{SWAP}=\sum {i,j=1}^{d_S}{|i⟩}_S⟨j|\otimes |j⟩⟨i|\otimes I_{2\dots n},$$

$$\mathsf{SHIFT}=I_S\otimes \sum _{i_1,\dots ,i_n=1}^{d_S}|i_2⟩⟨i_1|\otimes |i_1⟩⟨i_2|\otimes |i_2⟩⟨i_3|\otimes \dots \otimes |i_1⟩⟨i_n|.$$

$${\varrho }_E(\tau )=P_2\otimes P_3\otimes \dots \otimes P_n\otimes {\varrho }_S(0).$$

$${\varrho }_E(2\tau )=P_3\otimes \dots \otimes P_n\otimes {\varrho }_S(0)\otimes P_1.$$

$$p(\{E_i{\}}_{i=1}^n)=\prod _{i=1}^n\mathrm{tr}[{\varrho }_S(i\tau )P_i]=1$$

$${\varrho }_E(n\tau )={\varrho }_S(0)\otimes P_1\dots \otimes P_{n-1}.$$

$$p(\{E_i{\}}_{i=n+1}^{2n})=\mathrm{tr}[{\varrho }_S(0)P_{n+1}]\prod _{i=1}^{n-1}\mathrm{tr}[P_i{P}_{n+i+1}].$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\log p(\{E_i{\}}_{i=n+1}^{2n})\le \log \left[\underset{n\to \infty }{\lim }\frac{1}{n}\left(\mathrm{tr}[{\varrho }_S(0)P_{n+1}]+\sum _{i=1}^{n-1}\mathrm{tr}[P_i{P}_{n+i+1}]\right)\right]=\log \frac{1}{d_S}.$$

$$d_{ER}(ϵ)=\underset{0<\alpha <1}{\min }\frac{\sqrt{1-\alpha }}{{ϵ}^{\alpha /2(1-\alpha )}}\exp \left[n\gamma T\frac{(\gamma \tau {)}^{\alpha -1}-\alpha }{1-\alpha }\right],$$

$$H_{S+S_1+R}={\sigma }_z\otimes I\otimes I+{\sigma }_x\otimes I\otimes I+I\otimes {\sigma }_z\otimes I+I\otimes {\sigma }_x\otimes I+{\sigma }_z\otimes {\sigma }_z\otimes I+0.3I\otimes {\sigma }_z\otimes {\sigma }_z+0.3I\otimes {\sigma }_y\otimes {\sigma }_y+0.3I\otimes {\sigma }_x\otimes {\sigma }_x.$$

$$p(H|\{E_i{\}}_{i=1}^n)=\frac{p(\{E_i{\}}_{i=1}^n|H)p(H)}{\int p(\{E_i{\}}_{i=1}^n|H)p(H)dH}=\frac{1}{Z}p(\{E_i{\}}_{i=1}^n|H)$$

$$Q_{\{{\kappa }_{\mu \nu },{\sigma }_{\mu \nu },{\varkappa }_{\mu \nu },{\varsigma }_{\mu \nu }\}}(H)=\prod _{\mu ,\nu }\left[\frac{1}{\sqrt{2\pi }{\sigma }_{\mu \nu }}\exp \left(-\frac{(\mathrm{Re}{H}_{\mu \nu }-{\kappa }_{\mu \nu }{)}^2}{2{\sigma }_{\mu \nu }^2}\right)\right]\prod _{\mu ,\nu }\left[\frac{1}{\sqrt{2\pi }{\varsigma }_{\mu \nu }}\exp \left(-\frac{(\mathrm{Im}{H}_{\mu \nu }-{\varkappa }_{\mu \nu }{)}^2}{2{\varsigma }_{\mu \nu }^2}\right)\right],$$

$$D_{\mathrm{KL}}\left(Q_{\{{\kappa }_{\mu \nu },{\sigma }_{\mu \nu },{\varkappa }_{\mu \nu },{\varsigma }_{\mu \nu }\}}(H)||p(H|\{E_i{\}}_{i=1}^n)\right)=D_{\mathrm{KL}}\left(Q_{\{{\kappa }_{\mu \nu },{\sigma }_{\mu \nu },{\varkappa }_{\mu \nu },{\varsigma }_{\mu \nu }\}}(H)||p(\{E_i{\}}_{i=1}^n|H)\right)+\log Z.$$

$$F(\kappa_{\mu\nu},\sigma_{\mu\nu},\varkappa_{\mu\nu},\varsigma_{\mu\nu})=-\sum_{\mu\nu}\log\sigma_{\mu\nu}-\sum_{\mu\nu}\log\varsigma_{\mu\nu}-\mathbb{E}_{\scriptsize\begin{array}[]{c}\xi_{\mu\nu}\sim{\cal N}(0,1)\\ \zeta_{\mu\nu}\sim{\cal N}(0,1)\\ \end{array}}\left[p\Big{(}\{E_{i}\}_{i=1}^{n}\Big{|}\{\kappa_{\mu\nu}+\xi_{\mu\nu}\sigma_{\mu\nu}+i(\varkappa_{\mu\nu}+\zeta_{\mu\nu}\varsigma_{\mu\nu})\}_{\mu,\nu}\Big{)}\right],$$

$$\displaystyle p\big{(}\varrho_{S}(t)\big{|}\{E_{i}\}_{i=1}^{n}\big{)}$$

$${\Phi }_S(t)[{\varrho }_S(0)]={\mathrm{tr}}_{ER}\left[\exp (t{\mathcal{L}}_{S+ER})[{\varrho }_S(0)\otimes {\varrho }_{ER}(0)]\right].$$

Full text

Introduction.—

Quantum systems are never perfectly isolated which makes the study of open quantum dynamics important for various disciplines including solid-state physics takahashi-2008 , quantum chemistry valkunas-2013 , quantum sensing degen-2017 , quantum information transmission wilde-2017 , and quantum computing nielsen-2000 . Open quantum dynamics is a result of interaction between the system of interest and its environment. It is usually assumed that the environment is an infinitely large reservoir in statistical equilibrium, which has a well-defined interaction with the system schoeller-2018 . However, the environments of many physical systems are rather complex and structured piilo-2011 ; cirac-2011 ; ma-2012 ; hoope-2012 ; yang-2013 ; hughes-2015 ; eisert-2015 ; cirac-2017 ; wittemer-2018 ; wang-2018 ; peng-2018 ; haase-2018 ; mascherpa-2019 . A model of the system-environment interaction is often heuristic and oversimplified (e.g., a harmonic environment), but even in this case the analysis is rater complicated and requires some elaborated analytical and numerical methods strathearn-2018 ; pollock-2019 ; altaisky-2017 . A theoretical model may also neglect some additional sources of decoherence and relaxation. The experimental analysis of the environmental degrees of freedom is difficult because of their inaccessibility in practice. In fact, one can only get some information about the actual environment by probing the system paris-2018 ; bennink-2019 . Therefore, one faces an important problem to learn the unknown environment and its interaction with the quantum system by probing and affecting the system only.

This problem can be partly solved within the assumption of fast bath relaxation, when the system density operator $\varrho_{S}$ experiences the semigroup dynamics $\varrho_{S}(t)=e^{{\cal L}_{S}t}\varrho_{S}(0)$ with the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) generator ${\cal L}_{S}$ gks-1976 ; lindblad-1976 . In this case, the generator is reconstructed by performing a process tomography of the channel $\Phi_{S}(t_{1})=e^{{\cal L}_{S}t_{1}}$ for a fixed time $t_{1}>0$ Tomography ; howard-2006 . The actual dynamics does not usually reduce to a semigroup though de-vega-2017 ; li-2018 ; fc-2018 . The problem of learning the environment is mostly attributed to memory effects accompanying the non-Markovian dynamics. In this case, one can still resort to the process tomography of channels $\Phi_{S}(t_{1})$, $\Phi_{S}(t_{2})$, $\ldots$, $\Phi_{S}(t_{K})$ by preparing various initial system states $\varrho_{S}(0)$ and performing different measurements on the system at time moments $t_{1}<t_{2}<\ldots<t_{K}$. This procedure is time consuming because one has to gather enough statistics for all time moments (the total number of required measurements is $Kd_{S}^{8}/\epsilon^{2}$ for a $d_{S}$-dimensional quantum system and the accuracy $\epsilon$ of statistical reconstruction bogdanov-2013 ; haah-2017 ). Moreover, the tomographic reconstruction of each channel $\Phi_{S}(t_{i})$ implies resetting the environment in the same initial state after each measurement, which is difficult to control in the experiment especially for a strong coupling between the system and environment.

Recently proposed methods exploit the transfer tensor techniques Cerrilo-2014 ; TTM ; TTMTomography to learn the Nakajima-Zwanzig equation nakajima-1958 ; zwanzig-1960 $\frac{d}{dt}\varrho_{S}(t)=\int_{0}^{t}{\cal K}(t-t^{\prime})\varrho_{S}(t^{\prime})dt^{\prime}$ and the recurrent neural networks RNN for defining Lindblad operators and learning the convolutionless master equation $\frac{d}{dt}\varrho_{S}(t)={\cal L}_{S}(t)\varrho_{S}(t)$. An implementation of the latter approach in practice encounters the same difficulties related with the necessity to perform state tomography at different time steps.

In this Letter, we develop a method to learn the effective Markovian embedding MarkovEmb ; xue-2015 ; xue-2017 ; campbell-2018 ; bennink-2019 ; luchnikov-2019 for non-Markovian processes instead of learning the master equation for the system ($S$). Within such an approach, the environment is effectively divided into two parts: the first one carries memory of the system and is responsible for non-Markovian dynamics [effective reservoir ($ER$)]; the second one is memoryless and causes Markovian decoherence and dissipation of $S+ER$. The system evolution reads

[TABLE]

where the generator ${\cal L}_{S+ER}$ governs dissipative and decoherence processes on the system and the effective reservoir.

A division of the environment into two parts is similar to the pseudomode method imamoglu-1994 ; garraway-1997 ; mazzola-2009 , the reaction coordinate model iles-smith-2014 ; iles-smith-2016 , and the non-Markovian core model tamascelli-2018 , where one derives a Markovian master equation in the GKSL form for the extended system comprising the system and a finite number of auxiliary modes. In spin-bosonic models, the Markovian embedding is justified if the bath correlation function has exponentially damped correlations mascherpa-2019 . However, for power-law bath correlation functions with long-range tails polyakov-2019 the number of auxiliary modes diverges, which limits applicability of the Markovian embedding at a long timescale.

The density operator $\varrho_{S}(t)$ is unaccessible in a single measurement though, so any relevant information about the system is only gained in a series of measurements. On the other hand, measurement interventions into the system evolution complicate the analysis due to the no-information-without-disturbance principle. Consider a series of projective measurements performed on the system at different times $t_{1}<t_{2}<\ldots<t_{n}$, with the measurement basis being chosen randomly, see Fig. 1. The measurement outcomes seem to be completely uninformative due to the intrinsic probabilistic nature of quantum mechanics and the wave function collapse at each measurement as an example in Fig. 1 suggests. However, such a series of measurement outcomes does contain some information because the outcomes at each time moment are not equiprobable but appear in accordance with the Born rule. In this Letter, we demonstrate that a sufficiently long series of measurement results has a pattern that can be recognized by a machine carleo-2019 . This is a sharp distinction from conventional tomographic approaches based on numerous repetitions of identical experiments to gather enough statistics.

Our algorithm maximizes the likelihood of observed measurement outcomes and provides the generator ${\cal L}_{S+ER}$ for any fixed dimension of the effective reservoir $d_{ER}$, which is a hyperparameter. Computationally, the optimal $d_{ER}$ corresponds to the maximal likelihood on the validation set, which prevents overfitting supplemental . Physically, the sufficient value of $d_{ER}$ can be estimated through a reduced set of parameters: the system-environment coupling strength, the reservoir correlation time, the cutoff frequency of the spectral function, and the system’s number of degrees of freedom interacting with the environment luchnikov-2019 ; supplemental . Alternatively, $d_{ER}$ can also be estimated via the ensemble learning method shrapnel-2018 .

If the system evolution is Markovian ($d_{ER}=1$), then the result of measurement at time $t_{k}$ depends on the measurement outcome at time $t_{k-1}$ only and does not depend on results of earlier measurements at times $t_{k-2},t_{k-3},\ldots$ lindblad-1979 . Instead, the non-Markovian dynamics is accompanied by correlations in the measurement outcomes lindblad-1979 ; pollock-2018 ; modi-2018 ; budini-2018 ; taranto-2019 , which can be analyzed via the process matrix costa-2016 and the process tensor milz-2017 . The process tensor is a particular form of a quantum network chiribella-2009 , which is defined through the generator ${\cal L}_{S+ER}$ in our model, see Fig. 2. The reconstruction of a general process tensor requires exponentially many projective measurements milz-2018 . However, the process tensor has a peculiar form in our model and depends on the generator ${\cal L}_{S+ER}$ only, so it can be reconstructed by maximizing the likelihood of getting the observed outcomes for a single series of measurements without resorting to the full quantum tomography.

Likelihood function and its gradient.—

Suppose the experimental setup allows for projective measurements of the system at times $t_{i}=i\tau$, $i=1,\ldots,n$, with the measurement basis $\{|\varphi_{k}^{(i)}\rangle\}_{k=1}^{d_{S}}$ being randomly chosen at each time moment $t_{i}$. Observation of the particular measurement outcome $k_{i}$ transforms the system state into $|\varphi_{k_{i}}^{(i)}\rangle\langle\varphi_{k_{i}}^{(i)}|$. Denote $E_{i}=|\varphi_{k_{i}}^{(i)}\rangle\langle\varphi_{k_{i}}^{(i)}|\otimes I_{ER}$ the projector acting on the system and effective reservoir. The collection of projectors $\{E_{i}\}_{i=1}^{n}$ is the training dataset that feeds the learning algorithm.

A superoperator $\Phi=\exp\left(\tau{\cal L}_{S+ER}\right)$ governs the system and the effective reservoir evolution in between two sequential measurements. The probability to get the particular sequence of measurement outcomes $\{k_{i}\}_{i=1}^{n}$ (the data $\{E_{i}\}_{i=1}^{n}$) equals luchnikov-2017

[TABLE]

The likelihood (3) admits alternative useful forms. Let $\Phi^{{\dagger}}$ be dual to $\Phi$ comment , then one can split Eq. (3) after $m$th measurement and get p={\rm tr}\big{[}\widetilde{\varrho}_{S+ER}(t_{m}){\cal E}_{S+ER}(t_{m})\big{]}, where the recurrence relation $\widetilde{\varrho}_{S+ER}(t_{i+1})=E_{i}\Phi[\widetilde{\varrho}_{S+ER}(t_{i})]E_{i}$ with $\widetilde{\varrho}_{S+ER}(0)=\varrho_{S+ER}(0)$ defines the forward propagation of the subnormalized density operator along the tensor network in Fig. 3(a) and the recurrence relation ${\cal E}_{S+ER}(t_{i-1})=\Phi^{{\dagger}}[E_{i}{\cal E}_{S+ER}(t_{i})E_{i}]$ with ${\cal E}_{S+ER}(t_{n})=I_{S+ER}$ defines the backward propagation molmer-2013 of effects in the Heisenberg picture along the tensor network in Fig. 3(b). This leads to a “sandwich” formula

[TABLE]

which is valid for all $m=1,\ldots,n$, see Fig. 3(c).

The likelihood function is to be maximized over parameters of the generator ${\cal L}_{S+ER}$ defining $\Phi=\exp(\tau{\cal L}_{S+ER})$. Such a maximization is the most common approach in supervised machine learning Bishop . The problem is that not every generator ${\cal L}_{S+ER}$ defines a legitimate (completely positive and trace preserving) map $\Phi$. To overcome this obstacle and simplify the implementation of the gradient ascent method CVX , we use the Stinespring dilation for the channel $\Phi$ [see, e.g., Holevo and Fig. 3(c)]:

[TABLE]

where $\varrho_{A}$ is a fixed pure state of the $d_{A}$-dimensional ancilla ($A$), $d_{A}=(d_{S}d_{ER})^{2}$, $U(H)=\exp(-iH\tau)$ is a unitary evolution operator, and $H$ is the effective Hamiltonian of $S+ER+A$. Eq. (5) guarantees $\Phi$ is completely positive and trace preserving provided $H$ is Hermitian. The ancillary operator $\varrho_{A}$ plays the role of a renewable subenvironment in quantum collision models rau-1963 ; scarani-2002 ; filippov-2017 and memoryless (Markovian) part of the environment shrapnel-2018 .

Because of the Stinespring dilation, the likelihood function is now to be maximized over parameters of the effective Hamiltonian, i.e., matrix elements $H_{\mu\nu}=\langle\mu|H|\nu\rangle$ of $H$ in some computational basis $\{|\mu\rangle\}_{\mu=1}^{d_{S}d_{ER}d_{A}}$. This means that parameters $H_{\mu\nu}$ are iteratively adjusted in the direction of the gradient of the logarithmic likelihood $g_{\mu\nu}=\frac{\partial\log p}{\partial H_{\mu\nu}}$. Since the likelihood function is the $n$-degree monomial with respect to both operators $U(H)$ and $U^{{\dagger}}(H)$, we readily get supplemental

[TABLE]

where the derivative $\frac{\partial U(H)}{\partial H_{\mu\nu}}$ is expressed through the spectral decomposition $H=\sum_{k}\lambda_{k}|\psi_{k}\rangle\langle\psi_{k}|$ as supplemental

[TABLE]

Keeping in a computer memory the operators $\widetilde{\varrho}_{S+ER}(t_{i})$ and ${\cal E}_{S+ER}(t_{i})$ for forward and backward propagations, respectively, we efficiently calculate the gradient in $O(n)$ steps. Since $\log p$ is a highly nonlinear and nonconvex function with respect to parameters $H_{\mu\nu}$, its optimization is accompanied with overcoming the convergence to local extremums and the slow convergence rate. In what follows, we use techniques that were shown to perform well in such nonconvex optimization problems as neural network learning jain-2017 .

Learning algorithm.—

The learning algorithm, which estimates the generator ${\cal L}_{S+ER}$ based on the training dataset $\{E_{i}\}_{i=1}^{n}$, is as follows grigoriev_github :

Fix the hyperparameter $d_{ER}$. Initialize the model by randomly choosing the factorized state $\varrho_{S+ER}(0)=\varrho_{S}(0)\otimes\varrho_{ER}(0)$ and the factorized Hamiltonian $H=H_{S+ER}\otimes I_{A}$. 2. 2.

Calculate the forward-propagation operators $\{\widetilde{\varrho}_{S+ER}(t_{i})\}_{i=1}^{n}$ and the backward-propagation operators $\{{\cal E}_{S+ER}(t_{i})\}_{i=0}^{n-1}$. 3. 3.

Calculate the likelihood (4). 4. 4.

Find the spectral decomposition of the $d_{S}d_{ER}d_{A}$-dimensional operator $H$ and calculate $\frac{\partial U(H)}{\partial H_{\mu\nu}}$ via (7). 5. 5.

Estimate the gradient (Likelihood function and its gradient.—) via a batch of summands and results of items 2, 3, 4. 6. 6.

Feed the estimated gradient to an advanced optimization method [e.g., the adaptive moment estimation (Adam) algorithm Adam ] and get the increment $\Delta H$. 7. 7.

Update the Hamiltonian $H\rightarrow H+\Delta H$ and repeat items 2–6 until the likelihood converges. 8. 8.

Make use of the final update of $H$ to find the channel $\Phi$ and the generator ${\cal L}_{S+ER}=\frac{1}{\tau}\ln\Phi$.

Synthetic data generation.—

We apply the learning algorithm above to the in silico training set $\{E_{i}\}_{i=1}^{n}$ generated in a non-Markovian composite bipartite collision model lorenzo-2017 . We consider a bipartite system $S+S_{1}$ composed of the very open qubit system under study $S$ and one auxiliary qubit system $S_{1}$. The bipartite system successively interacts with identical subenvironments during some collision time supplemental . Such a model is quite rich and describes, e.g., a qubit subject to random telegraph noise. The benefit of this model is that the measurement interventions into the system evolution are explicitly taken into account supplemental .

Validation.—

We run the learning algorithm for various values of the hyperparameter $d_{ER}=1,2,4,6$ on the generated training set $\{E_{i}\}_{i=1}^{n}$, $n=10^{5}$ supplemental . The value $d_{ER}=1$ corresponds to the best Markovian approximation for the dynamics that is most compatible with the observed measurement outcomes. However, the likelihood for $d_{ER}=1$ is less than that for non-Markovian models with $d_{ER}\geq 2$, see Fig. 4(a). The greater $d_{ER}$, the wider the complexity class of possible dynamics luchnikov-2019 . If $d_{ER}=d_{S}^{n}$, then any series of projectors $\{E_{i}\}_{i=1}^{n}$ can be perfectly reconstructed with the likelihood $p(\{E_{i}\}_{i=1}^{n})=1$, which is an ultimate case of overfitting supplemental , see Fig. 4(b). The maximally achieved values of the logarithmic likelihood $\log p(\{E_{i}\}_{i=1}^{n})$ on the training set monotonically increase with the increase of $d_{ER}$. To avoid overfitting, we calculate the likelihood (3) on a separate validation set of projectors $\{E_{i}\}_{i=n+1}^{2n}$. Fig. 4(a) shows that, for the data analyzed, the logarithmic likelihood $\log p(\{E_{i}\}_{i=n+1}^{2n})$ on the validation set increases up to $d_{ER}=2$ and then diminishes. The Markovian embedding with $d_{ER}=2$ is the simplest model that is the most compatible with the observed series of measurement outcomes. This is an expected result because we used the synthetic data generated within a collision model with qubits, $d_{S_{1}}=2$. For real experimental data, the hyperparameter $d_{ER}$ is tuned in such a way that the likelihood on the validation set achieves its maximum. Tuning is reasonable to perform in the vicinity of the physical estimate for $d_{ER}$ derived in Ref. luchnikov-2019 .

Results.—

With the estimated generator ${\cal L}_{S+ER}$ at hand, we predict the open system dynamics $\varrho_{S}(t)$ by Eqs. (1) and (2) and compare it with the exact theoretical model (with no measurement interventions). The missing initial state of the effective reservoir is chosen to be the equilibrium state ${\rm tr}_{S}[\varrho_{S+ER}^{\infty}]$ such that ${\cal L}_{S+ER}[\varrho_{S+ER}^{\infty}]=0$. The results are depicted in Fig. 5. Good agreement between the estimated dynamics and the exact one demonstrates that the presented learning algorithm actually extracts useful information from the correlation pattern in a sequence of measurements on the open quantum system.

The quality of the estimated dynamics is assessed in two ways. (i) If the exact dynamics $\Phi_{S}(t)$ is known, we calculate the distinguishability between the estimated dynamics and the exact one, then average over time moments within the interval $[0,T]$. The result is $\varepsilon=0.03$ for $T=50$ supplemental . (ii) If the exact solution is not known, the quality of the estimated dynamics is assessed within the variational Bayesian inference approach. This approach yields $\varepsilon=0.05$ for $T=50$ and the standard deviation 0.025 for matrix elements of the estimated density operator $\varrho_{S}(T)$ supplemental .

The average error in estimating the discretized process $\{\Phi_{S}(t_{i})\}_{i=1}^{K}$ scales as $1/\sqrt{n}$ and is essentially independent of $K$ in the proposed algorithm supplemental because all the channels $\{\exp({\cal L}_{E+ER}\Delta t_{i})\}_{i=1}^{K}$ in the process tensor in Fig. 1 are defined by a single generator ${\cal L}_{S+ER}$ independent of time moments $\{t_{i}\}_{i=1}^{K}$ (parameter sharing). On the other hand, the full process tomography yields the error scaling as $\sqrt{K/n}$ with the same total number of measurements $n$ supplemental . Therefore, the proposed method is $\sqrt{K}$ times more efficient as compared to the full process tomography for large $K$.

Importantly, the formalism of Markovian embedding is compatible with a control operation on system $S$, say, a quick unitary transformation $\varrho_{S}(t^{\prime})\rightarrow V\varrho_{S}(t^{\prime})V^{{\dagger}}$ at time moment $t^{\prime}$. After the operation, \varrho_{S}(t)={\rm tr}_{ER}\big{\{}\exp[(t-t^{\prime}){\cal L}_{S+ER}]\varrho_{S+ER}(t^{\prime})\big{\}}. The result is in good agreement with the exact dynamics (Fig. 6), thus opening an avenue toward efficient control and manipulation of non-Markovian quantum systems. In contrast, the conventional process tomography cannot take such a control operation into account: its prediction $\Phi_{S}(t)\Phi_{S}(t^{\prime})^{-1}[V\varrho_{S}(t^{\prime})V^{{\dagger}}]$ differs from $\varrho_{S}(t)$ because of the system-environment correlations supplemental ; gessner-2011 ; rivas-2014 ; milz-2019 , see Fig. 6.

Conclusions.—

We proposed a method to learn the Markovian embedding for non-Markovian quantum evolution. The primary information needed is the outcomes of successive projective measurements on the system. Correlations in the measurements at different times indicate non-Markovianity and allow for the reconstruction of memory effects. The decay of correlations between spaced-in-time measurements enables the reconstruction of relaxation effects. Both memory and relaxation phenomena are taken into account by the generator ${\cal L}_{S+ER}$ acting on the system and the effective reservoir of finite dimension. Our algorithm estimates ${\cal L}_{S+ER}$ and does not exploit the full tomography of either states or processes. Learnability of the algorithm is tested on a dataset for the non-Markovian qubit dynamics. The presented approach enables to take control on the system into consideration, which is impossible with conventional tomographic techniques.

Acknowledgements.

The authors thank Peter Staňo for useful comments. Conceptualization, implementation, and validation of the learning algorithm is supported by the Russian Foundation for Basic Research under Project Nos. 18-37-00282 and 18-37-20073 and is performed in Moscow Institute of Physics and Technology, Skolkovo Institute of Science and Technology, and Valiev Institute of Physics and Technology of Russian Academy of Sciences, where S.N.F. was partially supported by Program No. 0066-2019-0005 of the Russian Ministry of Science and Higher Education. Synthetic data generation is performed in Moscow Institute of Physics and Technology, where I.A.L. and S.N.F. were partially supported by the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS” under Grant No. 19-1-2-66-1. The study of the likelihood function and quantum control is supported by the Russian Science Foundation under Project No. 17-11-01388 and is performed in Steklov Mathematical Institute of Russian Academy of Sciences.

.1 Derivation of the likelihood function gradient

The likelihood function $p({E_{i}}_{i=1}^{n})$ can be rewritten in many alternative ways with the help of the forward propagation operators $\widetilde{\varrho}_{S+ER}(t_{i})$ and the backward propagation operators ${\cal E}_{S+ER}(t_{i})$. In fact, $\widetilde{\varrho}_{S+ER}(t_{i})$ is the subnormalized state of $S+ER$ at time moment $t_{i}$ such that ${\rm tr}[\widetilde{\varrho}_{S+ER}(t_{m})]=p(\{E_{i}\}_{i=1}^{m})$. The total likelihood for $n$ outcomes equals $p(\{E_{i}\}_{i=1}^{n})={\rm tr}[\widetilde{\varrho}_{S+ER}(t_{n})]={\rm tr}[\widetilde{\varrho}_{S+ER}(t_{n})I_{S+ER}]={\rm tr}[\widetilde{\varrho}_{S+ER}(t_{n}){\cal E}_{S+ER}(t_{n})]$, where we have introduced the “initial” condition for backward propagation ${\cal E}_{S+ER}(t_{n})=I_{S+ER}$.

Forward propagation is given by the recurrence relation

[TABLE]

with $\widetilde{\varrho}_{S+ER}(0)=\varrho_{S+ER}(0)$. The last equality in (8) is due to the Stinespring dilation for the channel $\Phi$ exploiting an ancilla $A$ of dimension $d_{A}=(d_{S}d_{ER})^{2}$ (the maximum Kraus rank for the channel $\Phi$ Holevo ). Tensor diagram for Eq. (8) is depicted in Fig. 7. Ancillary systems $A$ play the role of particles colliding with $S+ER$ (see, e.g., fsp-2020 ).

Similarly, the operator ${\cal E}_{S+ER}(t_{m})$ propagates backward, i.e., in the Heisenberg picture, and is given by the recurrence relation

[TABLE]

with ${\cal E}_{S+ER}(t_{n})=I_{S+ER}$. The last equality in (9) explicitly uses the Stinesping dilation for the dual channel $\Phi^{{\dagger}}$. Tensor diagram for Eq. (9) is depicted in Fig. 8.

Merging the forward and backward propagations at a fixed time $t_{m}$, we express the likelihood function in many various though equivalent forms, namely,

[TABLE]

The latter expression (10) is a “sandwich” composed of the forward propogation till time $t_{m-1}$ [the state $\widetilde{\varrho}_{S+ER}(t_{m-1})$], the backward propagation till time $t_{m}$ [the operator ${\cal E}_{S+ER}(t_{m})$], and the unitary transformation $U(H)\,\cdot\,U^{{\dagger}}(H)$ followed by the $m$th measurement in between. Since the likelihood function is the $n$-degree monomial with respect to both operators $U(H)$ and $U^{{\dagger}}(H)$, we readily express its gradient operator as follows:

[TABLE]

The operator $\frac{\partial p(\{E_{i}\}_{i=1}^{n}|H)}{\partial H}$ is Hermitian provided the Hamiltonian $H$ is Hermitian. The derivative $\frac{\partial U(H)}{\partial H}=\left.\frac{\partial U(H+V)}{\partial V}\right|_{V=0}$. We use the perturbation expansion $\exp[-i(H+V)\tau]=\exp(-iH\tau)T_{\leftarrow}\exp\left[-i\int_{0}^{\tau}\exp(iHt^{\prime})V\exp(-iHt^{\prime})dt^{\prime}\right]=\exp(-iH\tau)-i\exp(-iH\tau)\int_{0}^{\tau}\exp(iHt^{\prime})V\exp(-iHt^{\prime})dt^{\prime}+o(V)$ and the spectral decomposition $H=\sum_{k}\lambda_{k}|\psi_{k}\rangle\langle\psi_{k}|$ to get

[TABLE]

Finally, we find the gradient of the logarithmic likelihood $\log p({E_{i}}_{i=1}^{n}|H)$ with respect to the unknown parameters $H_{\mu\nu}$:

[TABLE]

Note that (13) is insensitive to the normalization of $p(\{E_{i}\}_{i=1}^{n}|H)$, which significantly simplifies the calculation.

.2 Details on the learning algorithm

At item 1 of the algorithm, we initialize the model by randomly choosing the factorized state $\varrho_{S+ER}(0)=\varrho_{S}(0)\otimes\varrho_{ER}(0)$ and the factorized Hamiltonian $H=H_{S+ER}\otimes I_{A}$. Starting with a Hamiltonian, which is factorized with respect to $S+ER$ and $A$, fastens the learning process of memory effects. Otherwise, the correlations between $S+ER$ and $A$ induce irreducible decoherence and dissipation on $S+ER$ that smear out the memory effects.

Typical learning curves in Fig. 9 show how the logarithmic likelihood for the training set increases during the learning process and approaches the theoretical prediction, whereas the likelihood for the validation set starts to decrease after some point in the case of overfitting ($d_{ER}=6$). We test multiple variations of the batch size and the Adam optimizer parameters Adam (used at step 6 of the learning algorithm) to determine the fastest algorithm convergence. The tuned parameters are $\beta_{1}=0.9$, $\beta_{2}=0.95$, $\epsilon=10^{-4}$, the learning rate is $10^{-3}$, the batch size is $10^{3}$.

.3 Details on overfitting

The greater the dimension of the effective reservoir $d_{ER}$, the greater the likelihood for the training set [see Fig. 4(a) in the main text]. However, this leads to overfitting because the likelihood for the validation set starts decreasing with the growth of $d_{ER}$ above the optimal value. A discrepancy between the likelihood for the training set and that for the validation set is a direct indication of the bias-variance tradeoff in machine learning Bishop . Fig. 10 demonstrates the effect of overfitting on the quality of the estimated open dynamics for the example studied. Non-optimal hyperparameter $d_{ER}=6$ leads to redundant oscillations in the predicted dynamics (as opposed to the optimal value $d_{ER}=2$).

Fig. 4(b) in the main text illustrates the perfect overfitting of the learning algorithm if the dimension of the effective reservoir $d_{ER}=d_{S}^{n}$, where $n$ is the number of projectors in the data set $\{E_{i}\}_{i=1}^{n}$. Recall that $E_{i}=P_{i}\otimes I_{ER}$, with $P_{i}=\ket{\varphi_{k_{i}}^{(i)}}\bra{\varphi_{k_{i}}^{(i)}}$ being a pure state of $S$. Consider the $n$-partite environment of the form

[TABLE]

and the time-independent unitary transformation for $S+E$

[TABLE]

Note that $\varrho_{S}(\tau)={\rm tr}_{E}(W\varrho_{S}(0)\otimes\varrho_{E}W^{{\dagger}})=P_{1}$ so the first projective measurement on the system in the basis $\ket{\varphi_{k_{1}}^{(1)}}$ would produce the specific outcome ($P_{1}$) with certainty, i.e., with probability 1. The state of environment after the first measurement on the system is

[TABLE]

Therefore, $\varrho_{S}(2\tau)={\rm tr}_{E}(W\varrho_{S}(\tau)\otimes\varrho_{E}(\tau)W^{{\dagger}})=P_{2}$, i.e., the outcome for effect $P_{2}$ will be observed with certainty. The state of environment after the first measurement on the system is

[TABLE]

One can continue the same line of reasoning until all $n$ measurements are performed. As a result we get

[TABLE]

and the final state of the environment is

[TABLE]

This scenario corresponds to the perfect overfitting and yields the logarithmic likelihood $\log p(\{E_{i}\}_{i=1}^{n})=0$ [depicted in the top right corner of Fig. 4(a) in the main text].

However, if more than $n$ measurements in random bases are performed, then the proposed effective reservoir of dimension $d_{S}^{n}$ fails in reproducing the results perfectly. In fact, continue the process above with an extended series of measurements $\{E_{i}\}_{i=n+1}^{2n}$, then

[TABLE]

If the measurement bases are chosen randomly (Haar measure), then the average $\braket{{\rm tr}[P_{i}P_{n+i+1}]}=\frac{1}{d_{S}}$. Concavity of the logarithm implies

[TABLE]

In other words, the regularized logarithmic likelihood on the validation set $\{E_{i}\}_{i=n+1}^{2n}$ tends to a value not exceeding $\log\frac{1}{d_{S}}$, which is approximately $-0.69$ for $d_{S}=2$ [see the bottom right corner of Fig. 4(a) in the main text].

.4 Details on the hyperparameter $d_{ER}$

In general, the hyperparameter $d_{ER}$ is tuned in such a way that the likelihood on the validation set achieves its maximum. Tuning is reasonable to perform in the vicinity of the physical estimate luchnikov-2019

[TABLE]

where $\epsilon$ is the desired accuracy of open dynamics, $n$ is the effective subsystem’s number of degrees of freedom interacting with the reservoir ($n\leq d_{S}^{2}$), $\gamma$ is the coupling strength (between the open system and the reservoir), $T$ is the reservoir correlation time, and $\tau$ is the minimal timescale for the open system dynamics (the inverse of the cutoff frequency for the reservoir spectral function).

For some physical systems $d_{ER}$ is defined by the very reservoir structure. For instance, the nuclear spin $I=1$ in nitrogen ${}^{14}N$ is as an effective reservoir for the electronic spin qubit in a nitrogen-vacancy center in diamond haase-2018 , which implies $d_{ER}=2I+1=3$.

.5 Details on the data generation

The synthetic training set $\{E_{i}\}_{i=1}^{n}$ is generated in a non-Markovian composite bipartite collision model lorenzo-2017 . We consider a bipartite system that successively interacts with qubit subenvironments $R$ during collision time $\Delta t$, with the bipartite system being composed of the very open qubit system under study $S$ and one auxiliary qubit system $S_{1}$, see Fig. 11. The composite bipartite collision model lorenzo-2017 allows to find the system evolution $\varrho_{S}(t)={\rm tr}_{S_{1}}[\varrho_{S+S_{1}}(t)]$ intervened by measurements on the system.

We fix the (dimensionless) interaction Hamiltonian between $S$, $S_{1}$, and $R$ in the form

[TABLE]

The coefficients in the interaction Hamiltonian correspond to the case when strong memory effects are present in the evolution whereas the relaxation time is much longer than the recurrence time of memory effects — the hardest open dynamics to reconstruct. Each collision results in the transformation \varrho_{S+S_{1}}(t+\Delta t)={\rm tr}_{R}\big{[}\exp(-iH_{S+S_{1}+R}\Delta t)\varrho_{S+S_{1}}(t)\otimes\varrho_{R}\exp(iH_{S+S_{1}+R}\Delta t)\big{]}.

To simulate projective measurements in random bases we proceed as follows. Suppose the qubit system is in the state $\varrho_{S}(t_{i})$ at time $t_{i}=i\tau$. We randomly choose a direction ${\bf r}^{(i)}\in\mathbb{R}^{3}$, $|{\bf r}^{(i)}|=1$, on a Bloch ball and calculate eigenvectors $|\varphi_{+}^{(i)}\rangle$ and $|\varphi_{-}^{(i)}\rangle$ of the polarization operator $r_{x}^{(i)}\sigma_{x}+r_{y}^{(i)}\sigma_{y}+r_{z}^{(i)}\sigma_{z}$, where $(\sigma_{x},\sigma_{y},\sigma_{z})$ is the conventional set of Pauli operators. The transformation $\{\pm\}\rightarrow|\varphi_{\pm}^{(i)}\rangle\langle\varphi_{\pm}^{(i)}|$ is an observable at time $t_{i}=i\tau$. One of the two measurement outcomes $\{\pm\}$ is accepted, with the probability to accept the result $+$ being $\langle\varphi_{+}^{(i)}|\varrho_{S}(t)|\varphi_{+}^{(i)}\rangle$. As a result, one of the operators $|\varphi_{\pm}^{(i)}\rangle\langle\varphi_{\pm}^{(i)}|\otimes I_{ER}$ is accepted as $E_{i}$. Observation of the outcome $\pm$ in the $i$th measurement of the system results in the transformation $\varrho_{S+S_{1}}\rightarrow|\varphi_{\pm}^{(i)}\rangle\langle\varphi_{\pm}^{(i)}|\otimes\varrho_{S_{1}}^{\pm}/{\rm tr}[\varrho_{S_{1}}^{\pm}]$, where $\varrho_{S_{1}}^{\pm}=(\langle\varphi_{\pm}^{(i)}|\otimes I_{S_{1}})\varrho_{S+S_{1}}(t)(|\varphi_{\pm}^{(i)}\rangle\otimes I_{S_{1}})$. The measurement is followed by another collision described above, which in turn is followed by a measurement, and so on until the set $\{E_{i}\}_{i=1}^{n}$ is completed.

.6 Variational Bayesian inference approach

In the presented learning algorithm, we maximize the likelihood function and find parameters $H_{\mu\nu}$ encoding the desired generator ${\cal L}_{S+ER}$ for the Markovian embedding. However, the algorithm itself does not provide the error (variance) of parameters $H_{\mu\nu}$. This error can, however, be estimated via the variational Bayesian method as follows. Let $p(H)$ be an a priori distribution for Hamiltonian $H$, say, a uniform distribution on ${\rm Re}H_{\mu\nu}$ and ${\rm Im}H_{\mu\nu}$ within a wide range. For an observed sequence of operators $\{E_{i}\}_{i=1}^{n}$ we get the a posteriori distribution

[TABLE]

with some constant $Z$. Although $p(H|\{E_{i}\}_{i=1}^{n})$ is not known precisely, we expect that for sufficiently big data set this distribution can be approximated by a factorized Gaussian distribution

[TABLE]

where the parameters $\kappa_{\mu\nu}+i\varkappa_{\mu\nu}$ define the optimal values $H_{\mu\nu}$ maximizing the likelihood, and the standard deviations $\sigma_{\mu\nu}$ and $\varsigma_{\mu\nu}$ define the accuracy of parameter estimation for the real and imaginary part of $H_{\mu\nu}$, respectively. If the number of measurements $n\rightarrow\infty$, then $\sigma_{\mu\nu}\rightarrow 0$ and $\varsigma_{\mu\nu}\rightarrow 0$. Our goal is to find $\sigma_{\mu\nu}$ and $\varsigma_{\mu\nu}$ for a finite $n$. To do so we minimize the Kullback–Leibler divergence $D_{\rm KL}$ of $p(H|\{E_{i}\}_{i=1}^{n})$ from $Q_{\{\kappa_{\mu\nu},\sigma_{\mu\nu},\varkappa_{\mu\nu},\varsigma_{\mu\nu}\}}(H)$, where $D_{\rm KL}(q(x)\,||\,p(x))=\int dx\,q(x)[\log q(x)-\log p(x)]$. One can readily see that

[TABLE]

Since $Z$ is independent of $H$ or either of $\kappa_{\mu\nu},\sigma_{\mu\nu},\varkappa_{\mu\nu},\varsigma_{\mu\nu}$, the minimization of (28) reduces to the minimization of $D_{\rm KL}\left(Q_{\{\kappa_{\mu\nu},\sigma_{\mu\nu},\varkappa_{\mu\nu},\varsigma_{\mu\nu}\}}(H)\,||\,p(\{E_{i}\}_{i=1}^{n}|H)\right)$. Using the explicit form of the Gaussian distribution, we come the so-called reparameterization trick molchanov-2017 ; kingma-2013 : the minimization of $D_{\rm KL}\left(Q_{\{\kappa_{\mu\nu},\sigma_{\mu\nu},\varkappa_{\mu\nu},\varsigma_{\mu\nu}\}}(H)\,||\,p(\{E_{i}\}_{i=1}^{n}|H)\right)$ is equivalent to the minimization of the functional

[TABLE]

where ${\cal N}(0,1)$ is the standard normal distribution. The expectation value in the right hand side of Eq. (29) is readily estimated by sampling from the standard normal distribution, i.e., $\mathbb{E}_{\xi\sim{\cal N}(0,1)}f(\xi)\approx\frac{1}{M}\sum_{j=1}^{M}f(\xi^{(j)})$, where $\xi^{(j)}$ is a sample from ${\cal N}(0,1)$. With such an estimation at hand, the minimization of (29) is readily performed numerically, which yields the optimal approximation $Q_{\{\kappa_{\mu\nu}^{\rm opt},\sigma_{\mu\nu}^{\rm opt},\varkappa_{\mu\nu}^{\rm opt},\varsigma_{\mu\nu}^{\rm opt}\}}(H)$.

A posteriori distribution of the system density operator $\varrho_{S}(t)$ at time $t$ is

[TABLE]

Sampling from the distribution (30) for a given time $t$, we get the standard deviation for matrix elements of the system density operator $\varrho_{S}(t)$. The results are depicted in Fig. 12 for short and long timescales. The maximum standard deviation for matrix elements of the estimated density operator equals 0.025 at time moments $t=50$.

.7 Error estimation

The presented learning algorithm uses $n$ projective measurements on the system to estimate the Markovian embedding for the non-Markovian system dynamics. This reconstruction results in the process tensor depicted in Fig. 2 in the main text. As a result, we can infer a desired number of channels $\Phi_{S}(t_{i})$, $i=1,\ldots,K$, for the system dynamics from time [math] to time $t_{i}$ by formula

[TABLE]

By the Choi–Jamiołkowski isomorphism, all the information about the channel $\Phi_{S}(t)$ is contained in the matrix $\Omega_{\Phi_{S}(t)}$ defined through (see, e.g., Holevo ; filippov-jms-2019 )

[TABLE]

where $\ket{\psi_{+}}=\frac{1}{\sqrt{d_{S}}}\sum_{i=1}^{d_{S}}\ket{i}\otimes\ket{i}$ is the maximally entangled state.

Provided the exact dynamical map $\Phi_{S}^{\rm exact}(t)$ is known, the error in estimating the channel $\Phi_{S}(t)$ can therefore be expressed as

[TABLE]

where $\|A\|_{1}={\rm tr}\sqrt{A^{{\dagger}}A}$. The average error in estimating a set of channels $\{\Phi_{S}(t_{i})\}_{i=1}^{K}$ equals

[TABLE]

Fig. 13 shows the error $\varepsilon(\{\Phi_{S}(t_{i})\}_{i=1}^{K})$ scales with $n$ as $\frac{1}{\sqrt{n}}$ and is almost independent of $K$.

If the exact quantum dynamical map $\Phi_{S}(t)$ is not known (as it takes place for experimental data), the error is estimated via the variational Bayesian inference approach. In full analogy with the previous section, we get the posterior distribution of channels $\Phi_{S}(t)$ and the corresponding Choi operators, namely,

[TABLE]

By sampling from the latter distribution, we calculate the average error for the proposed learning algorithm

[TABLE]

where $\mathbb{E}\Phi_{S}(t_{i})$ is the mean Bayesian inference (obtained via averaging over samples). Fig. 14 shows that, in this case, the error $\varepsilon(\{\Phi_{S}(t_{i})\}_{i=1}^{K})$ scales with $n$ as $\frac{1}{\sqrt{n}}$ and is almost independent of $K$.

.8 Comparison with the full process tomography

The standard quantum-process tomography exploits an ensemble of identically prepared quantum systems corresponding to a given experimental setting (see, e.g., the review knee-2018 ). Suppose the total number of available projective measurements is $n$. As we are interested in reconstructing $K$ channels $\{\Phi_{S}(t_{i})\}_{i=1}^{K}$ with the minimal possible average error, the number of projective measurements per each channel equals $\frac{n}{K}$. The theory of process tomography bogdanov-2013 ; haah-2017 predicts the reconstruction error $\epsilon\sim\frac{d_{S}^{4}}{\sqrt{n/K}}$ in this case. In what follows, we confirm this prediction numerically.

Ref. knee-2018 proposes an algorithm that maximizes the likelihood for the observed measurement outcomes and provides a legitimate (trace preserving and completely positive) estimate for $\Phi_{S}^{\rm tom}(t_{i})$. Dealing with a qubit dynamical map ($d_{S}=2$), one needs to prepare the system in one of four pure initial states $\{\varrho_{S}^{(j)}(0)\}_{j=1}^{4}$, where $\varrho_{S}^{(1)}(0)=\ket{0}\bra{0}$, $\varrho_{S}^{(2)}(0)=\ket{1}\bra{1}$, $\varrho_{S}^{(3)}(0)=\frac{1}{2}(\ket{0}+\ket{1})(\bra{0}+\bra{1})$, $\varrho_{S}^{(4)}(0)=\frac{1}{2}(\ket{0}+i\ket{1})(\bra{0}-i\bra{1})$. Then the randomly chosen state $\varrho_{S}^{(j)}(0)$ is evolved through the channel $\Phi_{S}^{\rm tom}(t_{i})$ and is measured with the help of an 8-outcome positive operator-valued measure (POVM) with effects $\{F_{k}\}_{k=1}^{8}$, where $F_{2m-1}=\frac{1}{4}\varrho_{S}^{(m)}(0)$ and $F_{2m}=\frac{1}{4}[I-\varrho_{S}^{(m)}(0)]$, $m=1,\ldots,4$. After the measurement outcome is read out, the environment should be reset to the initial (thermal equilibrium) state and the system should be again prepared in one of the states $\{\varrho_{S}^{(j)}(0)\}_{j=1}^{4}$. This is a challenge in real experimental setup (especially in the case of strong coupling between the system and environment) and a disadvantage as compared to our proposed scheme of sequential measurements with no environment resets. Suppose, however, that the experiment is repeated $n/K$ times. This results in integers $\{n_{jk}\}_{j=1,\ldots,4,\ k=1,\ldots,8}$, which quantify how many times the outcome $k$ is observed provided the system is prepared in the state $\varrho_{S}^{(j)}(0)$. Clearly, the relative frequencies $\frac{4Kn_{jk}}{n}$ tend to probabilities p_{jk}={\rm tr}\Big{[}\Phi_{S}^{\rm tom}(t_{i})[\varrho_{S}^{(j)}(0)]F_{k}\Big{]} if $n\rightarrow\infty$. The authors of Ref. knee-2018 maximize the likelihood $\prod_{jk}p_{jk}^{n_{jk}}$ with respect to $\Phi_{S}(t_{i})$ and find the best estimate $\Phi_{S}^{\rm tom}(t_{i})$ for the quantum channel. We use the solver in Ref. knee-2018 to find $\Phi_{S}^{\rm tom}(t_{i})$ for a given number of measurements $n/K$.

Suppose the exact dynamical map $\Phi_{S}(t)$ is known, then the average reconstruction error equals

[TABLE]

Fig. 15 shows the error $\varepsilon^{\rm tom}(\{\Phi_{S}(t_{i})\}_{i=1}^{K})$ scales as $\sqrt{K/n}$, which is in strong contrast to Fig. 13.

.9 Details on the coherent control

Within the process tensor formalism for Markovian embedding, we are able to describe the action of coherent control pulses on the system experiencing a non-Markovian dynamics. Suppose the system is subjected to a quick unitary transformation $\varrho_{S}(t^{\prime})\rightarrow V\varrho_{S}(t^{\prime})V^{{\dagger}}$ at time moment $t^{\prime}$. The corresponding tensor diagram is depicted in Fig. 16.

The estimated system dynamics after the operation is given by the equation \varrho_{S}(t)={\rm tr}_{ER}\big{\{}\exp[(t-t^{\prime}){\cal L}_{S+ER}]\varrho_{S+ER}(t^{\prime})\big{\}}. Fig. 17 illustrates the qubit evolution for $V=\sigma_{x}$ and $t^{\prime}=20$. The exact dynamics and the estimated dynamics are in good agreement with each other.

Suppose the full process tomography is performed for time moments $t_{1},t_{2},\ldots,t_{K}$ and the maps $\Phi_{S}(t_{1}),\Phi_{S}(t_{2}),\ldots,\Phi_{S}(t_{K})$ are reconstructed precisely. The dynamics within the time interval $[t_{i},t_{i+1}]$ is given by the intermediate map $\Lambda(t_{i+1},t_{i})=\Phi_{S}(t_{i+1})\Phi_{S}^{-1}(t_{i})$. Denoting $t_{0}=0$, we note that $\Phi_{S}(t_{j})=\Lambda(t_{j},t_{j-1})\cdots\Lambda(t_{2},t_{1})\Lambda(t_{1},t_{0})=:\bigcirc_{i=0}^{j-1}\Lambda(t_{i+1},t_{i})$, i.e., the dynamics is described by concatenation of intermediate maps. If a coherent control gate $V$ is applied at time moment $t_{m}$, then the concatenation approach yields

[TABLE]

This approach results in the dynamics depicted in Fig. 17 by dots. Clearly, Eq. (38) is not able to reproduce the system dynamics after the control gate is applied because $\Phi_{S}(t)\Phi_{S}(t^{\prime})^{-1}[V\varrho_{S}(t^{\prime})V^{{\dagger}}]\neq\varrho_{S}(t)$ for $t>t^{\prime}$ due to the system-environment correlations gessner-2011 ; rivas-2014 ; milz-2019 .

Non-monotonicity of the trace distance $\frac{1}{2}\|\varrho_{S}^{\prime}(t)-\varrho_{S}^{\prime\prime}(t)\|_{1}$ for some initial states $\varrho_{S}^{\prime}(0)$ and $\varrho_{S}^{\prime\prime}(0)$ is a clear indication of non-Markovianity BLP , and the learned Markovian embedding reproduces such a non-monotonic behavior quite well lvgf-2019 .