Learning Quantum Processes with Memory - Quantum Recurrent Neural Networks
Dmytro Bondarenko
Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstraße 2, 30167 Hannover, Germany
Stewart Blusson Quantum Matter Institute, University of British Columbia, 2355 East Mall, V6T 1Z4 Vancouver, BC, Canada<EMAIL_ADDRESS>Robert Salzmann
Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstraße 2, 30167 Hannover, Germany
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, United Kingdom<EMAIL_ADDRESS>Viktoria-S. Schmiesing
Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstraße 2, 30167 Hannover, Germany<EMAIL_ADDRESS>
Recurrent neural networks play an important role in both research and industry. With the advent of quantum machine learning, the quantisation of recurrent neural networks has become recently relevant. We propose fully quantum recurrent neural networks, based on dissipative quantum neural networks, capable of learning general causal quantum automata. A quantum training algorithm is proposed and classical simulations for the case of product outputs with the fidelity as cost function are carried out. We thereby demonstrate the potential of these algorithms to learn complex quantum processes with memory in terms of the exemplary delay channel, the time evolution of quantum states governed by a time-dependent Hamiltonian, and high- and low-frequency noise mitigation. Numerical simulations indicate that our quantum recurrent neural networks exhibit a striking ability to generalise from small training sets.

Machine learning (ML) has made remarkable progress during the last decades [1, 2, 3]. Central among the many key successes of ML are recurrent neural networks (RNNs): These networks are especially suited to the learning of sequential data generated in a wide variety of natural settings, with use cases including time evolution [4], speech recognition [5], the prediction of electric power demand [6], and machine translation [7]. Various RNN architectures have been proposed, ranging from fundamental RNNs to the more complicated, such as LSTMs [8], and gated recurrent units (GRUs) [9, 10]. However, challenges remain as ML algorithms are typically computationally expensive to implement.

Quantum computing promises new techniques to bolster the faltering Moore's law [11]. This has led, in the context of ML, to the advent of
quantum machine learning (QML) [12]. A rough taxonomy of QML algorithms may be provided by categorising them according to whether they exploit either quantum processing (CQ), quantum data (QC), or both (QQ) [13, 14]. CQ ML is built on quantum algorithms to learn classical data, mainly aimed at speeding up classical learning algorithms [15, 16, 17, 18]. QC ML utilises classical machine learning for quantum problems, such as quantum metrology [19], simulating quantum many-body systems [20], or adaptive quantum computation [21]. QQ ML directly targets the development of quantum learning algorithms to analyse quantum data [22, 23, 24, 25, 26, 27, 28, 29, 30, 31].

Many QML architectures would benefit from the ability to repeatedly query a (quantum) memory storage via a recurrent network structure. We say the system has memory if the output at time $t$ depends on the inputs at times $\tau \leq t$.

(0,0.75) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U000);
(0,0) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U001);
(1,1.125) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U010);
(1,0.375) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U011);
(1,-0.375) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U012);
(2,0.75) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U020);
(2,0) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U021);
(U000.north east)–(U010.west);
(U001.north east)–(U010.south west);
(U000.east)–(U011.north west);
(U001.east)–(U011.south west);
(U000.south east)–(U012.north west);
(U001.south east)–(U012.west);
(U010.east)–(U020.north west);
(U011.north east)–(U020.west);
(U012.north east)–(U020.south west);
(U010.south east)–(U021.north west);
(U011.south east)–(U021.west);
(U012.east)–(U021.south west);
[->,orange,thick] (-0.75,0)–(U001.west);
[->,blue,thick] (-0.75,0.75)–(U000.west);
[->,purple,thick] (U021.east)–(10.75,0);
(3,1.5) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U100);
(3,0.75) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U101);
(4,1.875) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U110);
(4,1.125) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U111);
(4,0.375) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U112);
(5,1.5) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U120);
(5,0.75) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (U121);
(U100.north east)–(U110.west);
(U101.north east)–(U110.south west);
(U100.east)–(U111.north west);
(U101.east)–(U111.south west);
(U100.south east)–(U112.north west);
(U101.south east)–(U112.west);
(U110.east)–(U120.north west);
(U111.north east)–(U120.west);
(U112.north east)–(U120.south west);
(U110.south east)–(U121.north west);
(U111.south east)–(U121.west);
(U112.east)–(U121.south west);
[->,orange,thick] (U020.east)–(U101.west);
[->,blue,thick] (-0.75,1.5)–(U100.west);
[->,purple,thick] (U121.east)–(10.75,0.75);
[->,orange,thick] (U120.east)–(5.75,1.5);
(8,3.75) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (Un00);
(8,3) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (Un01);
(9,4.125) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (Un10);
(9,3.375) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (Un11);
(9,2.625) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (Un12);
(10,3.75) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (Un20);
(10,3) node[circle,minimum height=0.3cm,minimum width=.3cm,draw] (Un21);
(Un00.north east)–(Un10.west);
(Un01.north east)–(Un10.south west);
(Un00.east)–(Un11.north west);
(Un01.east)–(Un11.south west);
(Un00.south east)–(Un12.north west);
(Un01.south east)–(Un12.west);
(Un10.east)–(Un20.north west);
(Un11.north east)–(Un20.west);
(Un12.north east)–(Un20.south west);
(Un10.south east)–(Un21.north west);
(Un11.south east)–(Un21.west);
(Un12.east)–(Un21.south west);
[->,orange,thick] (7.25,3)–(Un01.west);
[->,blue,thick] (-0.75,3.75)–(Un00.west);
[->,purple,thick] (Un21.east)–(10.75,3);
[->,orange,thick] (Un20.east)–(10.75,3.75);
[dotted, thick, purple] (10.5,1)–(10.5,2.75);
[dotted, thick, blue] (-0.25,1.75)–(-0.25,3.5);
[dotted, thick, orange] (6,1.75)–(7,2.75);
(-1.125,1.875) node[]()\(\rho^\text{IN}\);
[decorate, decoration=brace, yshift=2ex]  (11,3.75) – node[right=0.4ex] \(\rho^\text{OUT}\)  (11,-0.75);

A general Quantum Recurrent Neural Network. A QRNN is an iteration over the memory system (orange) of a feed-forward QNN with \(L\) hidden layers (\(L=1\) in Figure) where the total input is split into an input (blue) and the memory and the total output is split into an output (purple) and the memory. For a concrete number of iterations we can input any state \(\rho^\text{IN}\) on \(\h^\text{m}\otimes\h^{\text{in}}\otimes \dots\otimes\h^{\text{in}}\) and get out a state \(\rho^\text{OUT}\) on \(\h^{\text{out}}\otimes \dots\otimes\h^{\text{out}}\otimes \h^\text{m}\).

In the CQ context, there are attempts to use quantum computers, quantum channels [32] and states [33, 34, 35] for tasks where classical RNNs are often used, such as natural language processing [36, 37, 38, 39]. In this context, there are already several quantisations of RNNs [40, 41, 42, 43, 44], usually exploiting a substitution of the feed-forward NNs in RNNs, LSTMs or GRUs with some form of QNN. Typically such QNNs encode classical information into a quantum state, apply a variational circuit [45, 46, 47, 48, 37, 49] or a tunable evolution [50, 51, 52, 53, 54, 55, 44], and then carry out a measurement, followed by further classical processing. A similar approach can be pursued with non-RNN variational classes used for machine learning with memory, such as matrix-product operators [56, 57].
Non-universal CQ RNNs that exhibit contextuality were proven capable of better scaling of memory resources than their entirely classical counterparts [44]. Amazingly, these networks were also able to outperform widely-used classical architectures of similar size in a natural language translation task [44]. This result suggests a wide industrial appeal of QRNNs when fault-tolerant quantum computers become available.

Unfortunately, with state-of-the-art classical natural language processing ML models training hundreds of billions of parameters and employing highly connected sophisticated architectures [58, 59], it is unclear if CQ RNNs can be competitive in the current noisy intermediate-scale quantum (NISQ) era.

Unlike for CQ learning, there is no classical alternative to QQ ML. If the task calls for direct handling of quantum data, QQ ML can immediately improve the state-of-the-art even if applied to modestly-sized systems.
Many tasks require quantum data spanning simulation, metrology, error correction and denoising, and communication [56]. It is known that entangled states are essential to represent many-body correlations in materials, able to enhance variance scaling of an $n$-shot measurement from $\sim 1/n$ to $\sim 1/n^2$, allow for unprecedentedly secure key distribution, and can significantly speed up specific algorithms if properly shielded from noise. For complex systems in hard-to-characterise environments, data is often easier to procure than exact solutions. Thus, QQ ML is a natural approach to searching for good entangled states for many applications.

In the QQ scenario, some exciting tasks also require quantum memory. We call the memory quantum if the time dependence is facilitated by a subsystem that is able to store entangled states. One of the simplest examples where quantum memory enlarges the capabilities of quantum protocols is when an extra state does not change during a particular process yet facilitates it. This is precisely what a catalyst does — a widely applied concept in chemistry [60] and quantum information processing [61]. There are also many practically relevant tasks where it is desirable to add (quantum or classical) memory that can change during the process. These QQ tasks include modelling the time-dependent evolution of a quantum system in order to steer it in a desired direction. Further application areas arise in chemistry and molecular dynamics [62], quantum control [63, 64, 65], scheduling for quantum annealers [66], and the design of classical and quantum adaptive filters [67, 68, 69, 56].

However, the QQ scenario is targeted by comparatively fewer architectures. One example of such a fully quantum RNN was presented in [70] for systems of continuous variables. In [71], a graph quantum RNN was used to represent a Trotterised Hamiltonian evolution. For all but the first and last steps, the memory in [71] is treated as the entire input and output. Instead, we want to focus on qudits and having in- and outputs in every time step.

To build a fully quantum recurrent neural network (QRNN), we must pick a suitable neuron. Just like classical RNNs, neurons should not modify subsystems not connected to it, and big enough networks should be universal [72, 73, 74, 75]. We say that a QRNN architecture is universal if it can implement any causal quantum automaton, i.e. any quantum process where the output at time \(t\) only depends on the input at times \(t^\prime \leq t\) [76]. Implementing a dissipative quantum channel might be desirable even if the data is pure, as the ability to forget is quintessential for learning long-term dependencies [77] and robustness to data corruption [78]. These requirements are conveniently satisfied by dissipative QNNs from [79], which can represent any quantum channel. These QNNs are then connected as depicted in Fig. <ref> to build the QRNN. We apply our QRNN to learn sequential quantum data generated in various contexts, including the delay channel, the time evolution of a state governed by a time-dependent Hamiltonian, and high- and low-frequency noise mitigation. Our classical simulations demonstrate a striking capability of the QRNN architecture to generalise from small training sets. If desired, our QRNNs can also be restricted to classical processing or classical data.

The network architecture.
A feed-forward (dissipative) quantum neural network (QNN) for qudits as proposed in [79] is built from connected quantum perceptrons arranged in \(L+2\) layers, i.e.\(\, L\) hidden layers, one input and one output layer, with \(m_l\) perceptrons in layer \(l\). To the perceptron in the \(l^\text{th}\) layer at place \(j\) we assign a qubit and a unitary \(U_j^l\). The unitary acts on the corresponding perceptron qubit and the qubits in the \((l-1)^{\text{th}}\) layer connected to it. The unitaries of layer \(l\) can be pooled together to create a layer unitary \(U^l=U^l_{m_l}\dots U^l_1\), and these can again be pooled together to yield the network unitary \(\mathcal{U}=U^{L+1}\dots U^1\). For a total input state \(\rho^0\), the total output is

\begin{align}
\label{eq:a}
\rho^{L+1}=& \mathcal{N}_{\U}\left( \rho^{0} \right)\\ =&\mathrm{tr}_{0,...,L}(\mathcal{U}(\rho^{0}\otimes \ket{0\dots 0}_{1,...,L+1}\bra{0\dots 0})\mathcal{U}^\dagger),
\end{align}

which also can be written as \(\mathcal{E}^{L+1}\left(\dots \mathcal{E}^2 \left(\mathcal{E}^1\left(\rho^{0}\right)\right)\right)\) with the layer-to-layer channels \(\mathcal{E}^l\) being defined by \(\mathcal{E}^l(X^{l-1})=\mathrm{tr}_{l-1}(U^l(X^{l-1}\otimes \ket{0\dots 0}_{l}\bra{0\dots 0})U^{l^\dagger})\). This definition of a perceptron is naturally connected to the open-system representation of a quantum channel  [80, 56, 81].

The defining feature of a recurrent network is that some perceptron outputs are used again as part of the input in the same or previous layers. Without loss of generality, we exploit a subsystem of the total output as part of the input, which is often called the memory. Hence, the total input Hilbert space \(\h^\text{intotal}\) is split into an input Hilbert space \(\h^\text{in}\) and a memory Hilbert space \(\h^\text{m}\), and accordingly the total output Hilbert space \(\h^\text{outtotal}\) is split into an output Hilbert space \(\h^\text{out}\) and \(\h^\text{m}\). The memory part of the total output of a QNN is then input into the memory part of the total input of the QNN in the next iteration step.

We index the Hilbert spaces, states, channels and unitaries with the number of the iteration.
In general, we can allow an arbitrary input state \(\rho^\text{IN}\) on \(\h^\text{IN}_{N}=\h_0^\text{m}\otimes \h_1^\text{in}\otimes \dots \otimes \h_{N}^\text{in}\) giving rise to the output state

\begin{align}
\rho^\text{OUT} &:=\left(\left(\id^\text{out}_{1,...,N-1}\otimes \N_{\U}\right) \circ \dots \circ  \left(\N_{\U}\otimes \id^\text{in}_{2,...,N}\right)\right)\left( \rho^\text{IN} \right)\nonumber\\
&:=\left( \N_{\U_{N}}\circ \dots \circ  \N_{\U_1}\right)\left( \rho^\text{IN} \right)
\end{align}

on \(\h^\text{OUT}_{N}=\h_1^\text{out}\otimes \dots \otimes \h_{N}^\text{out} \otimes \h_{N}^\text{m} \). See Fig. <ref> for an illustration.
It is a direct consequence of the quantum-circuit structure of our QRNNs that they can be used to describe all causal quantum automata with finite in- and output systems. To argue this firstly note that the strategy above directly gives rise to concatenated quantum channels with memory, as described in, e.g. [76, 82, 83], only that the QNN with channel \(\N_{\U}\) in the above description is replaced with an arbitrary channel. Secondly, quantum channels with memory can be used to describe general causal quantum automata [76] where the memory system is finite in the case of finite in- and output systems [84]. Finally, a QNN, as described above, can, if it is large enough, represent any quantum channel on a finite input or output system [79]. Together, this means that our QRNNs can represent any causal quantum automaton if the underlying QNN is large enough.

The training process.
Now that we have an architecture for our QRNN, we can specify the learning task. To quantise the classical scenario, one might first replace each classical input and output with a quantum state. However, in the quantum mechanical setting, the different in- and outputs could also be entangled. In general, the training set is given by \(S=\left\{ S_\alpha \right\}_{\alpha=1}^{M}=  \left\{ \left( \rho^\text{IN}_\alpha, \sigma^\text{OUT}_\alpha \right) \right\}_{{\alpha}=1}^{M} \in \bigtimes_{{{\alpha}}=1}^{M} \left(\D (\h^\text{IN}_{N_\alpha})\times \D (\h^\text{OUT}_{N_\alpha} \right)) \) for \(M\) runs of a QRNN where we iterate over the underlying QNN \(N_{\alpha}\) times, and \(\D (\h)\) denotes the set of density operators on a Hilbert space \(\h\).  Note that the training set can contain states of varying dimensionality. Nevertheless, for the quantum algorithm to work, we need to be able to request multiple copies of every in- and output state.

We focus on the case of product in- and output states here. That is, the input state is given by \(\rho^\text{IN}_\alpha=\rho_{0 \, \alpha}^\text{m}\otimes \rho_{1\, \alpha}^\text{in}\otimes \dots \otimes \rho_{N_{\alpha} \, \alpha}^\text{in}\) and the output state by \(\sigma_\alpha^\text{OUT}=\sigma_{1 \, \alpha}^\text{out}\otimes \dots \sigma_{N_\alpha \, \alpha}^\text{out}\otimes \sigma_{N_{{\alpha}} \, \alpha}^\text{m}\). In this case, we often write \(S_\alpha=\left(\rho_{0 \, \alpha}^\text{m},\left(\rho_{1 \, \alpha}^\text{in},\rho_{1 \, \alpha}^\text{out}\right),\dots,  \left(\rho_{N_{{\alpha}} \, \alpha}^\text{in},\rho_{N_{{\alpha}} \, \alpha}^\text{out}\right) \right)\). The case of entangled inputs is left to a future publication.
To evaluate the performance of our QRNN for learning training data, i.e., to evaluate how close the network output $\rho_{\alpha}^{\text{OUT}}$ for the input $\rho_\alpha^\text{IN}$ is to the correct output $\sigma_\alpha^\text{OUT}$, we need a cost function. If $\sigma_\alpha^\text{OUT}$ is a product state, we can either compare $\rho_{\alpha}^{\text{OUT}}$ with $\sigma_\alpha^\text{OUT}$ locally or globally.
The local case corresponds to trying to locally learn the different states \(\sigma_{x \, \alpha}^\text{out}\). We compare the states \(\sigma_{x \, \alpha}^\text{out}\) with the reduced density operators on the \(x^\text{th}\) output in the \(\alpha ^\text{th}\) run, i.e.
\(\rho_{x \, \alpha}^\text{out}=\tilde{\mathrm{tr}}_{\text{out}}^{x\, {{\alpha}}}  \left(\rho_\alpha^\text{OUT}\right) =\tilde{\mathrm{tr}}_{\text{out}}^{x\, {{\alpha}}} \left( \M_{\U}\left(\rho_\alpha^\text{IN}\right)\right) , \)
where \(\tilde{\mathrm{tr}}\) traces out everything but the indicated layer. Operationally, the fidelity is the essentially unique measure of closeness for pure quantum states. Therefore, we use it to define our cost function for pure output states. Hence, as a direct generalisation of the risk function considered in training classical RNNs, for \(\sigma_{x \, \alpha}^\text{out}=\ket{\phi_{x \, \alpha}^\text{out}} \bra{\phi_{x \, \alpha}^\text{out}}\) the local cost is given by

\begin{equation}
C_\text{local,pure}=\frac{1}{M}\sum_{\alpha=1}^{M}{\frac{1}{N_{\alpha}}}\sum_{x=1}^{N_{{\alpha}}}\bra{\phi_{x \, \alpha}^\text{out}} \rho_{x \, \alpha}^\text{out} \ket{\phi_{x \, \alpha}^\text{out}}.
\end{equation}

We ignore the last memory state as it is unimportant to us in which state the memory is at the end. For mixed states the fidelity is given by \(F(\rho,\sigma)=\left(\mathrm{tr}\sqrt{\rho^{1/2}\sigma \rho^{1/2}}\right)^2\). A polynomial-time algorithm exists to estimate mixed state fidelity [85]. Nevertheless, the Hilbert-Schmidt distance is much easier to estimate and minimise in classical simulation or quantum implementation than the mixed state fidelity [85, 86]. Hence, we use the square of the Hilbert-Schmidt distance $\mathrm{tr}(\rho -\sigma)^2$ instead if \(\sigma_{x\, \alpha}^\text{out}\) is not pure in general.

The global case corresponds to trying to learn the overall product state. We compare the global output of the training set with the global output of the QRNN. For pure outputs, we use the fidelity between the global output state of the QRNN and the correct one. Again, the memory state at the end is not important to us. As a state is a product state, if the partial trace over a part of the system is a pure state, we can also neglect the state of the memory at the end and choose

\begin{align}
C_\text{global, pure}=& \frac{1}{M}\sum_{\alpha=1}^M \left( \bra{\phi_{1 \, \alpha}^\text{out}}\otimes \dots\otimes \bra{\phi_{N_{{\alpha}} \, \alpha}^\text{out}} \right) \nonumber\\ &\mathrm{tr}_\text{m}^N \left(\rho_\alpha^\text{OUT}\right)
\left( \ket{\phi_{1 \, \alpha}^\text{out}}\otimes \dots\otimes \ket{\phi_{N_{{\alpha}} \, \alpha}^\text{out}} \right)
\end{align}

as the cost function. We call this the global cost.

at (-1,3.6) a;
xlabel= number of training pairs,
ylabel= Cost,
legend columns=2,
legend style=at=(0.5,-0.3),anchor=north,
+[only marks,color=red,mark=triangle,mark size=2.9pt]  table[x=N,y=Cwomround,col sep=comma] Generalisationbehaviour.csv;
[scale=0.1, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.15) ;
+[only marks,color=blue,mark=triangle,mark size=2.9pt]   table[x=N,y=C1mround,col sep=comma] Generalisationbehaviour.csv;
[xscale=0.1,yscale=0.175, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5 pt] at (-1,0.3) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(5) – (0.3,0.3);
[out=0,in=0] (0.3,0.3) to (0.3,0.5);
(0.3,0.5) to (-1.3,0.5);
[out=180,in=180] (-1.3,0.5) to (-1.3,0.3);
(-1.3,0.3) – (2);
+[only marks,color=red,mark=square,mark size=2pt]   table[x=N,y=Cwomtestround,col sep=comma] Generalisationbehaviour.csv;
[scale=0.1, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.15) ;
+[only marks,color=blue,mark=square,mark size=2pt]   table[x=N,y=C1mtestround,col sep=comma] Generalisationbehaviour.csv;
[xscale=0.1,yscale=0.175, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5 pt] at (-1,0.3) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(5) – (0.3,0.3);
[out=0,in=0] (0.3,0.3) to (0.3,0.5);
(0.3,0.5) to (-1.3,0.5);
[out=180,in=180] (-1.3,0.5) to (-1.3,0.3);
(-1.3,0.3) – (2);
[domain=0:11, samples=10, color=gray,]0.5;
[domain=0:11, samples=10, color=gray,]1;

at (-1,6.6) b;
[] (3,3) circle (2.675);
[->] (3,3)–(3,6) node[above]\(tr(\rho \sigma^z)\);
[->] (3,3)–(6,3) node[right]\(tr(\rho \sigma^y)\);
[white] (-1,-1)–(6,-1);
[blue](4) [circle,draw,inner sep=0pt,minimum size=12pt] at (3.094,0.328) ;
[blue](5) [circle,draw,inner sep=0pt,minimum size=12pt] at (0.325,3.005) ;
[blue](6) [circle,draw,inner sep=0pt,minimum size=12pt] at (3.013,5.675) ;
[thick](1) [star,draw,inner sep=0pt,minimum size=12pt] at (3,5.675) \(1\);
[thick](2) [star,draw,inner sep=0pt,minimum size=12pt] at (0.325,3) \(2\);
[thick](3) [star,draw,inner sep=0pt,minimum size=12pt] at (3,0.325) \(3\);
[red](4) [circle,draw,inner sep=0pt,minimum size=12pt] at (0.526,1.994) \(1\);
[red](5) [circle,draw,inner sep=0pt,minimum size=12pt] at (0.389,2.609) \(3\);
[red](6) [circle,draw,inner sep=0pt,minimum size=12pt] at (0.469,2.35) \(2\);
[->] (2.8,5.77) arc (95:174:2.85);
[->] (0.23,2.8) arc (185:264:2.85);
[->] (0.12,3.25) arc (175:93.5:2.87);
[->] (2.75,0.12) arc (265:183.5:2.87);

[thick](7) [star,draw,inner sep=0pt,minimum size=12pt] at (2,-0.75);
(8) at (4,-0.75)training data;

[blue](10) [circle,draw,inner sep=0pt,minimum size=12pt] at (1.1,-1.7) ;
(11) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1.7,-1.85) ;
(12) [circle,draw,inner sep=0pt,minimum size=3.5 pt] at (1.7,-1.55) ;
(14) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (2.45,-1.85) ;
(15) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (2.45,-1.55) ;
(16) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (3.2,-1.85) ;
(17) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (3.2,-1.55) ;
(17) – (3.3,-1.55);
[out=0,in=0] (3.3,-1.55) to (3.3,-1.35);
(3.3,-1.35) to (1.6,-1.35);
[out=180,in=180] (1.6,-1.35) to (1.6,-1.55);
(1.6,-1.55) – (12);

[red](20) [circle,draw,inner sep=0pt,minimum size=12pt] at (4.1,-1.7) ;
(21) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (4.7,-1.7) ;
(24) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (5.45,-1.85) ;
(25) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (5.45,-1.55) ;
(26) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (6.2,-1.7) ;

at (-1,3.5) c;
xlabel= Sequence length,
ylabel= Cost,
legend columns=3,
legend style=at=(0.4,-0.3),anchor=north,
[thick,color=red,mark=None]  table[x=x,y=i1m0round,col sep=comma] fidelity_DriftAndSmooth.csv;
[scale=0.1, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.15) ;

[thick,color=magenta,mark=None]  table[x=x,y=i1m0h2round,col sep=comma] fidelity_DriftAndSmooth.csv;
[xscale=0.1,yscale=0.175, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(6) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0.15) ;

[thick,color=orange,mark=None]  table[x=x,y=i1m0h2h2round,col sep=comma] fidelity_DriftAndSmooth.csv;
[xscale=0.1,yscale=0.175, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(3) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0.3) ;
(6) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (2,0.15) ;

[thick,color=blue,mark=None]  table[x=x,y=i1m1round,col sep=comma] fidelity_DriftAndSmooth.csv;
[xscale=0.1,yscale=0.175, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5 pt] at (-1,0.3) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(5) – (0.3,0.3);
[out=0,in=0] (0.3,0.3) to (0.3,0.5);
(0.3,0.5) to (-1.3,0.5);
[out=180,in=180] (-1.3,0.5) to (-1.3,0.3);
(-1.3,0.3) – (2);

[ thick,color=cyan,mark=None]  table[x=x,y=i1m1h2round,col sep=comma] fidelity_DriftAndSmooth.csv;
[xscale=0.1,yscale=0.175, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.3) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(6) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0) ;
(7) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0.3) ;
(7) – (1.3,0.3);
[out=0,in=0] (1.3,0.3) to (1.3,0.5);
(1.3,0.5) to (-1.3,0.5);
[out=180,in=180] (-1.3,0.5) to (-1.3,0.3);
(-1.3,0.3) – (2);

[thick,color=ForestGreen,mark=None]  table[x=x,y=i1m2round,col sep=comma] fidelity_DriftAndSmooth.csv;
[xscale=0.1,yscale=0.175, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,-0.15) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(3) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.45) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,-0.15) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.15) ;
(6) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.45) ;
(6) – (0.3,0.45);
[out=0,in=0] (0.3,0.45) to (0.3,0.65);
(0.3,0.65) to (-1.3,0.65);
[out=180,in=180] (-1.3,0.65) to (-1.3,0.45);
(-1.3,0.45) – (3);
(5) – (0.35,0.15);
[out=0,in=0] (0.35,0.15) to (0.35,0.75);
(0.35,0.75) to (-1.35,0.75);
[out=180,in=180] (-1.35,0.75) to (-1.35,0.15);
(-1.35,0.15) – (2);

[domain=9:101, samples=10, color=gray,]0.5;

[]Comparison of QRNN and QNN.  Panel (a) shows the generalisation behaviour of a
[xscale=0.4,yscale=0.6, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5 pt] at (-1,0.3) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(5) – (0.3,0.3);
[out=0,in=0] (0.3,0.3) to (0.3,0.5);
(0.3,0.5) to (-1.3,0.5);
[out=180,in=180] (-1.3,0.5) to (-1.3,0.3);
(-1.3,0.3) – (2);
QRNN with one memory qubit (blue) and a [scale=0.5, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.15) ;
feed-forward QNN (red) for the delay by one channel. We trained the QNNs for $1000$ rounds with a learning rate of $\eta\epsilon=0.06$ and \(M=1\) sequence with $N = 1, 2, …, 10$ training pairs and evaluated the cost function on the training set (triangles) and for a set of $30$ test pairs (squares) afterwards,  we then averaged this over $50$ rounds. Panel (b) shows the $y,z$-plane of the Bloch sphere. A [xscale=0.4,yscale=0.6, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0) ;
(2) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.3) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(6) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0) ;
(7) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0.3) ;
(7) – (1.3,0.3);
[out=0,in=0] (1.3,0.3) to (1.3,0.5);
(1.3,0.5) to (-1.3,0.5);
[out=180,in=180] (-1.3,0.5) to (-1.3,0.3);
(-1.3,0.3) – (2);
QRNN with one memory qubit (blue) and a [xscale=0.4,yscale=0.6, baseline]
(1) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (-1,0.15) ;
(4) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0) ;
(5) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (0,0.3) ;
(6) [circle,draw,inner sep=0pt,minimum size=3.5pt] at (1,0.15) ;
feed-forward QNN (red) are trained over $1000$ training rounds with a learning rate of $\eta\epsilon=0.05$ to learn the time evolution, where the state goes from $\ket{0}$ (black point 1) over $\ket{\phi^-}$ (black star 2) to \(\ket{1}\) (black star 3), over $\ket{\phi^-}$ (black star 2) back to $\ket{0}$ (black star 3) and so on. To train the networks, we used an \(M=1\) sequence with $N=10$ training pairs and predicted the time evolution for $N_\text{test}=20$ points. The blue circles show the prediction for the QRNN and the red of the feed-forward QNN. In the latter case, the evolution approximately goes from 1 to 2 to 1 to 3 to 1 to 2 to 1 to 3 to 1, etc. Panel (c) shows performance of various networks for the low- and high-frequency noise mitigation task, see Eq. <ref>. The QRNNs are trained using $M=25$ sequences of $N_{\alpha } =25$ pairs each and are tested on sequences with length from $10$ to $100$ in the increments of $2$. We used $150$ training rounds with the Nadam [87] optimiser; hyperparameters are learning rate—$555.(5)$, momentum—$0.9$, RMSprop—$0.999$.

Thus the network can be trained with several cost functions. The cost function for pure outputs varies between $0$ (worst) and $1$ (best), and the one for mixed outputs between \(4\) (worst) and \(0\) (best). We can use the cost functions in different cases. Both have advantages and drawbacks: While the global cost is more accurate when it comes to high costs (this is because \(C_\text{global, pure}\leq C_\text{local, pure}\), as shown in Appendix <ref>) the local cost function is easier to train with, and fewer qubits and gates are needed to calculate it on a quantum computer. Hence, although using the global cost will also result in learning the correlations between the individual input and output systems, using the local cost is easier to realise on present and next-generation quantum information processing devices.

As in the classical case, we obtain the training algorithm by maximising/minimising the derivative of the cost function. The algorithms for classical simulations of the learning process can be found in Appendices <ref>, <ref>, <ref> and <ref> for the different costs. Mainly the one for the local cost with pure output states (Appendix <ref>) is used to obtain numerical results. The updated unitaries in this algorithm can be calculated in a layer-by-layer manner. Hence, only unitaries on \(m_l+m_{l+1}\) and not all qubits have to be applied in one go, i.e., the size of the matrices used to simulate the learning process only scales with the width and not the depth of the network. This reduces the requirements of the simulation and allows us to classically simulate the training of modest-sized QRNNs.

The quantum training algorithm may be found in Appendix <ref>.

Simulation of learning tasks.
For simplicity, we focus on pure qubit outputs and training with the local cost in the following. Extensive numerical results may be found in Appendix <ref>. The exponential growth of Hilbert space dimension with the width of the network restricts us to QRNNs with small widths in our classical simulations. We classically simulated the learning process for different tasks while we used total input and output spaces up to \(4\) qubits. We illustrate the results in terms of three exemplary tasks.

The first learning task is focused on the delay, SWAP or shift channel, which often serves as the first example in papers on quantum channels with memory, see, e.g. [82, 83, 76]. This channel is defined as follows: The \(x^\text{th}\) output equals the \((x-1)^\text{th}\) input, i.e. \(\rho_x^\text{in}=\ket{\phi_x^\text{in}} \bra{\phi_x^\text{in}}\) and  \(\ket{\phi_x^\text{out}} = \ket{\phi_{x-1}^\text{in}} \). We thus studied the ability of a QRNN to learn a training set of form \(S=\left(\ket{\phi_0^\text{m}},\left( \ket{\phi_1^\text{in}} ,\ket{\phi_0^\text{m}}\right),\left( \ket{\phi_2^\text{in}}, \ket{\phi_1^\text{in}} \right),\dots,  \left(\ket{\phi_N^\text{in}},\ket{\phi_{N-1}^\text{in}}\right) \right)\) with varying size \(N\) and to generalise to (bigger) test sets of the same form with different input states. We chose the input states at random with respect to the Haar measure. The generalisation behaviour is shown in Figure <ref> (a). The QRNN learns the SWAP channel perfectly and generalises well, even for only eight training pairs. A feed-forward QNN, by contrast, learns worse for every training pair added and does not generalise at all. This is what we expected as a feed-forward QNN can only guess the \(x^\text{th}\) output since it only has access to the \(x^\text{th}\) input and not the \((x-1)^\text{th}\).

The second task we consider is the learning of the time evolution in a time interval \(\left( T_2, T_3\right]\) from a time evolution of a state \(\ket{\psi (t) }\) in a time interval \(t\in\left[T_1, T_2\right]\), where \(T_1< T_2 < T_3\). Classically this is a prototypical example where RNNs are often more efficient than feed-forward NNs [88, 89, 4]. To carry out the quantum task, we exploited the training set \( S = ( \ket{\psi(T_1)} , ( \ket{\psi(T_1+\tau)},  \ket{\psi(T_1+2\tau)}  ), \dots , ( \ket{\psi(T_2-\tau)} , \allowbreak \ket{\psi(T_2)} ) )\) to train our network. We then use \( \ket{\psi(T_2)} \) as input of the network after \( \ket{\psi(T_2-\tau)} \). The goal is for the output state to be close to  \( \ket{\psi(T_2+\tau)} \). This state is then used as input again and so forth. Then the time predictions are compared to the correct time evolution via the fidelity. As an example we started in state \(\ket{0}\), then went over \(\ket{\phi^-}=\frac{1}{\sqrt{2}}\ket{0}-\frac{i}{\sqrt{2}}\ket{1}\) to \(\ket{1}\), back to \(\ket{0}\) over \(\ket{\phi^-}\) and so on. This is shown in Figure <ref> (b). As one can see, the training and prediction process works better with QRNNs than with feed-forward QNNs, as the blue circles for QRNNs are nearly coincident with the black stars (target states) and the red ones for feed-forward QNNs are completely different.

Learning how to filter low- and high-frequency noise in quantum data is the third task we considered. Consider an interferometric experiment that encodes information of interest into a relative phase $\phi$ in the state $| \phi \rangle \equiv \frac{| 0 \rangle + e^{i \phi} | 1 \rangle}{\sqrt{2}}$. Suppose that one desires a cheaper, less complex, or more transportable version of said experiment. Naturally, the relative phase in the output states will become noisier. We can circumvent this by designing an appropriate post-processing device trained with the supervision of an expensive state-of-the-art machine. Moreover, we might want to use denoised states for further quantum computation. In this case, QNNs are a capable platform for learning denoisers [90, 91, 56]. Shot-to-shot denoisers can have a feed-forward structure. However, operating in the frequency domain, e.g. to combat drifts, requires memory [56].

We train QRNNs on the data set $S = \{ S_\alpha \} = \left\{|0\rangle, \left(|\phi^\alpha_t \rangle, |\psi^\alpha_t\rangle \right)_{t=1}^{T_\alpha} \right\}_{\alpha=1}^{N}$. We construct a simple example of $\phi^\alpha_t$ that require bandwidth filtering. Consider noise consisting of two components. The low-frequency noise can be cancelled by subtracting a drift with speed $v$. The high-frequency noise is cancelled by exponential smoothing, i.e. taking a convex combination with previous correct outcomes. That is for all \(\alpha\)

\begin{align} \label{eq:filter}
\psi^\alpha_0 &= \phi^\alpha_0 \nonumber\\
\psi^\alpha_t &= w \cdot (\phi^\alpha_t - v \cdot t) + (1-w) \cdot \psi^\alpha_{t-1}, \ t>1.
\end{align}

This task can not be solved with fidelity one due to the no-cloning theorem. Nevertheless, QRNNs with sufficient memory can significantly improve the quality of the outputs, see Figure <ref>(c).

For an in-depth discussion and more numerical results, see Appendix <ref>.

In this paper, we have proposed a natural quantisation of recurrent neural networks. These are universal, i.e. they can be used to simulate general causal quantum automata. Our approach is scalable – only a few qubits must be stored to train the network. The number of qubits needed only scales with the width and not the network's depth. One can view our QRNNs as trainable matrix-product quantum channels. The experimental requirements for implementing QRNNs are comparable with those demonstrated in producing matrix-product states and density matrices [92, 93, 94, 95, 96]. Indeed, in [92], the authors have experimentally implemented a quantum channel with memory for 6 time steps on an ion-trap quantum computer. QRNNs are also suitable for some existing measurement-based photonic devices. These devices were used to demonstrate up to 100 sequential operations on a single logical qubit [97] and generate more than 1 million mode 1d cluster states [98]. As such, while suited even for current devices, QRNNs can leverage the rapid progress of quantum hardware [99].

We demonstrated that QRNNs can outperform feed-forward QNNs with respect to several exemplary learning tasks. For this, we have designed some benchmarking tasks and classically simulated the training of QRNNs. We started with simple yet paradigmatic examples, such as learning a delay channel or cancelling a linear drift. We then transitioned to relevant tasks for applications, such as learning time-dependent dynamics and bandwidth filters.

Many intriguing questions remain in the study of QRNNs: We have concentrated on the case of product in- and outputs; are there interesting tasks involving data entangled at different points in time (or space), and how well can QRNNs cope with such data? Theoretical results show the benefit of quantum memory for learning algorithms [100], so how large is this set of tasks? When is quantum memory necessary, and when can it be substituted by a classical memory [44]? Are QRNNs particularly suited for matrix-product states and density matrices, and is there a deeper connection to CC ML tensor network methods [101, 102]? If so, can QRNNs be used for state preparation and to extend the notion of the parent Hamiltonian [103, 104, 105, 106, 107] and Lindbladian [108, 109, 110]? Can one use QRNNs for the discovery of catalysts? Does recurrency influence robustness to noise of QNNs?

Helpful correspondence and discussions with Tobias J. Osborne, Kerstin Beer, Polina Feldmann, Thomas Cope, Gabriel Müller, Taras Kucherenko, Terry Farrelly, Yihui Quek, Ramona Wolf, Luis Mantilla Calderon, Robert Raußendorf are gratefully acknowledged. This work was funded by the Deutsche Forschungsgemeinschaft (German Research Foundation) under Germany's Excellence Strategy (EXC-2123, QuantumFrontiers Grant No. 390837967), through CRC 1227 (DQ-mat) within project No. A06, by the Quantum Valley Lower Saxony
(QVLS) through the Volkswagen Foundation, Canada First Research Excellence Fund, Quantum Materials and Future Technologies Program, and
the Cambridge Commonwealth, European and International Trust.

Data availability: We obtained all results using Python and Matlab. The code is available at <https://github.com/qigitphannover/DeepQuantumNeuralNetworks/tree/master/QRNN>.

§ COMPARISON OF THE DIFFERENT COSTS

Generally, the global cost will be smaller or equal to the average cost, which is the subject of the following lemma:

For all global input states \(\rho^\text{IN}\) with corresponding output \(\rho^\text{OUT}\) being a pure product state, i.e. \(\rho^{OUT}=  \kb{\Phi^\text{OUT}}=\bigotimes_{x=1}^N \kb{\phi_x^\text{out}}\in\h^{\text{out}^{\otimes n}}\), and network unitaries \(\U\) we have

\begin{equation}
\label{eq:costcomparison}
C_\text{global,pure}\leq C_\text{local,pure} .
\end{equation}

The fidelity of two quantum states, \(F(\rho,\sigma)=\left(\mathrm{tr}\sqrt{\rho^{1/2}\sigma \rho^{1/2}}\right)^2\), satisfies the data-processing inequality (see e.g. <cit.>), i.e. for all quantum channels \(\E\) we have

\begin{equation}
\label{eq:datafidelity}
F\left( \E(\rho),\E(\sigma)\right)\geq F\left(\rho,\sigma\right).
\end{equation}

Equation (<ref>) then simply follows by using (<ref>) for the quantum channel being the partial trace:

\begin{align*}
C_\text{local,pure}&=\frac{1}{N} \sum_{x=1}^{N}\bra{\phi_x^\text{out}} \rho_x^\text{out} \ket{\phi_x^\text{out}}= \frac{1}{N} \sum_{x=1}^{N} F\left(  \kb{\phi_x^\text{out}}, \tilde{\mathrm{tr}}_{\text{out}}^x  \left(\rho^\text{OUT}\right) \right)\\
&=\frac{1}{N} \sum_{x=1}^{N} F\left( \tilde{\mathrm{tr}}_{\text{out}}^x  \left(\kb{\Phi^\text{Out}}\right), \tilde{\mathrm{tr}}_{\text{out}}^x   \left(\mathrm{tr}_\text{m}^N \left(\rho^\text{OUT}\right)\right) \right) \\
&\geq \frac{1}{N} \sum_{x=1}^{N} F\left( \kb{\Phi^\text{Out}}, \mathrm{tr}_\text{m}^N \left(\rho^\text{OUT}\right)\right) = F\left( \kb{\Phi^\text{Out}}, \mathrm{tr}_\text{m}^N \left(\rho^\text{OUT}\right)\right) =C_\text{global,pure}.
\end{align*}

§ OPTIMISING THE LOCAL COST WITH PURE OUTPUT STATES CLASSICALLY

We now want to derive the algorithm to maximise the local cost with pure outputs. In a first step, we only consider \(M=1\) run of the QRNN in Appendices <ref> to <ref>, and only consider many runs with \(M > 1\) in Appendix <ref>. Hence, we can ignore \(\alpha\) for now. It does not matter for the classical simulation if we trace out everything but the \(x^\text{th}\) output at the end or if we calculate the reduced density operators \(\rho_x^\text{out}\) and \(\rho_x^\text{m}\) after each iteration step and use \(\rho_x^\text{m}\otimes \rho_{x+1}^\text{in}\) as total input in the next iteration step. The state \(\rho_x^\text{m}\) is then given by

\begin{align*}
\rho^\text{m}_x=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^x(\mathrm{tr}_{\text{m}}^{x-1}(\mathcal{U}_x(\rho_x^{\text{in}}\otimes \rho_{x-1}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^x\bra{0\dots 0})\mathcal{U}_x^\dagger))\\
=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^x(\mathrm{tr}_{\text{m}}^{x-1}(\mathcal{U}_x(\rho_x^{\text{in}}\otimes \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x-1}(\mathrm{tr}_{\text{m}}^{x-2}(\mathcal{U}_{x-1}(\rho_{x-1}^{\text{in}}\otimes \rho_{x-2}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{x-1}\bra{0\dots 0})\mathcal{U}_{x-1}^\dagger))\\
&\otimes \ket{0\dots 0}_{ {1:L+1}}^x\bra{0\dots 0})\mathcal{U}_x^\dagger)).
\end{align*}

For \(A\in \B(\h_A)\) and \(B\in \B(\h_A\otimes \h_B)\) for some Hilbert spaces \(\h_A\) and \(\h_B\) it is

\begin{align*}
\mathrm{tr}_B((A\otimes \id_B)B)&%=\sum_k (\id_A\otimes\bra{k})(A\otimes \id_B)B(\id_A\otimes\ket{k}) =\sum_k A(\id_A\otimes\bra{k})B(\id_A\otimes\ket{k})\\
%&= A\sum_k(\id_A\otimes\bra{k})B(\id_A\otimes\ket{k})
\end{align*}

Hence, as \(\U_{x}=\U_{x}\otimes \id^\text{m}_{x-2}\otimes \id^\text{in, hid, out}_{x-1}\) it follows that

\begin{align*}
\rho^\text{m}_x=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x,x-1}(\mathrm{tr}_{\text{m}}^{x-1,x-2}(\mathcal{U}_x(\rho_x^{\text{in}}\otimes (\mathcal{U}_{x-1}(\rho_{x-1}^{\text{in}}\otimes \rho_{x-2}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{x-1}\bra{0\dots 0})\mathcal{U}_{x-1}^\dagger))\\
&\otimes \ket{0\dots 0}_{ {1:L+1}}^x\bra{0\dots 0})\mathcal{U}_x^\dagger)).
\end{align*}

Also, \(\U_{x-1}=\U_{x-1}\otimes  \id^\text{in, hid, out,m}_{x}\) which leads to

\begin{align}
\rho^\text{m}_x=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x,x-1}(\mathrm{tr}_{\text{m}}^{x-1,x-2}(\mathcal{U}_x\mathcal{U}_{x-1}(\rho_x^{\text{in}}\otimes \rho_{x-1}^{\text{in}}\otimes \rho_{x-2}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{x-1,x}\bra{0\dots 0})\mathcal{U}_{x-1}^\dagger\mathcal{U}_x^\dagger))\nonumber\\
=& \dots\nonumber\\
=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:1}(\mathrm{tr}_{\text{m}}^{x-1:0}(\mathcal{U}_x\dots\mathcal{U}_{1}(\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{1:x}\bra{0\dots 0})\mathcal{U}_{1}^\dagger\dots\mathcal{U}_x^\dagger)), \label{eq:d}
\end{align}

where \(i:j\) is a short notation for \(i,\dots,j\).

The loss function is given by

\begin{align}
C=&\frac{1}{N}\sum_{x=1}^N\bra{\phi_x^\text{out}}\rho_x^\text{out}\ket{\phi_x^\text{out}} \nonumber\\
=&\frac{1}{N}\sum_{x=1}^N\bra{\phi_x^\text{out}}\mathrm{tr}_{\text{in,}1:L\text{,m}}^x(\mathrm{tr}_{\text{m}}^{x-1}(\mathcal{U}_x(\rho_x^{\text{in}}\otimes \rho_{x-1}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^x\bra{0\dots 0})\mathcal{U}_x^\dagger))\ket{\phi_x^\text{out}} \nonumber\\	=&\frac{1}{N}\sum_{x=1}^N \mathrm{tr}\left((\id_x^{\text{in,}1:L\text{,m}}\otimes \id_{x-1}^\text{m}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x(\rho_x^{\text{in}}\otimes \rho_{x-1}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^x\bra{0\dots 0})\mathcal{U}_x^\dagger\right) \label{eq:c}
\end{align}

and we write \(C=C_\text{local,pure}\) here for simplicity. In order to train the QRNN, we have to find the perceptron unitaries \(U_j^l\) that maximize the cost function. To do so, we use the same ansatz as in [79] and the unitaries are updated in each step according to

\[U_j^l\mapsto U_j^{l \prime}= e^{i\epsilon K_j^l}U_j^l\]

where the matrices \(K_j^l\) are hermitian and depend on the matrices \(U_j^l\). We now have to determine the matrices \(K_j^l\) that make the biggest positive change in the cost function. To do so, we have to maximise the derivative of the cost function at \(U_j^l\) when changing the matrices as presented above

\[\delta C= \lim_{\epsilon\to 0}\frac{\Delta C}{\epsilon}= \lim_{\epsilon\to 0}\frac{C^\prime- C}{\epsilon}. \]

In the expression for the cost function in equation \(\eqref{eq:c}\) the unitaries \(U_j^l\), that are the object of change, appear directly in the equation and indirectly in the state \(\rho_{x-1}^{\text{m}}\) that also depends on the unitaries (for \(x>1\)). So we first look at the change in \(\rho_{x-1}^{\text{m}}\) after a step of size \(\epsilon\). As described in equation \(\eqref{eq:d}\) it is

\begin{align*}
\rho^{\text{m}}_x=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:1}(\mathrm{tr}_{\text{m}}^{x-1:0}(U_{m_{L+1}\, x}^{L+1}\dots U_{1\,x}^1\dots U_{j\, z}^l\dots U_{m_{L+1}\, 1}^{L+1} \dots U^1_{1\,1}\\
&(\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ 1:L+1}^{1:x}\bra{0\dots 0})U_{1\,1}^{1^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger})),
\end{align*}

\begin{align*}
\rho^{\text{m}\prime}_x=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:1}(\mathrm{tr}_{\text{m}}^{x-1:0}(e^{i\epsilon K_{m_{L+1}}^{L+1}}U_{m_{L+1}\, x}^{L+1}\dots e^{i\epsilon K_{1}^{1}} U^1_{1\,1}(\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \\
&\otimes \ket{0\dots 0}_{ 1:L+1}^{1:x}\bra{0\dots 0}) U_{1\,1}^{1^\dagger}e^{-i\epsilon K_{1}^{1}} \dots U_{m_{L+1}\, x}^{L+1^\dagger}e^{-i\epsilon K_{m_{L+1}}^{L+1}})).
\end{align*}

Using the Taylor expansion of the exponential function it follows that

\begin{align*}
\rho^{\text{m}\prime}_x=&\mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:1}(\mathrm{tr}_{\text{m}}^{x-1:0}(U_{m_{L+1}\, x}^{L+1}\dots  U^1_{1\,1}(\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ 1:L+1}^{1:x}\bra{0\dots 0}) U_{1\,1}^{1^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger}))\\
&+i\epsilon \sum_{z=1}^{x}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{1:L+1}^{x:1}\left(\mathrm{tr}_{\text{m}}^{x-1:0} \left( \mathcal{U}_x\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l} K_{j\, z}^l U_{j\, z}^{l} \dots U_{1\,z}^1\mathcal{U}_{z-1}\dots \mathcal{U}_1 \right.\right.\\
&\cdot  (\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ 1:L+1}^{1:x}\bra{0\dots 0}) \U_1^\dagger \dots \U_{z-1}^\dagger U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger\\
&-\mathcal{U}_x\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}  U_{j\, z}^{l} \dots U_{1\,z}^1\mathcal{U}_{z-1}\dots \mathcal{U}_1 (\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{1:x}\bra{0\dots 0}) \\
&\cdot  \left.\left.\U_1^\dagger \dots \U_{z-1}^\dagger U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }K_{j\, z}^l U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger \right)\right)+\mathcal{O}(\epsilon^2).
\end{align*}

\begin{align}
\Delta \rho_x^m=&\rho^{\text{m}\prime}_x-\rho^\text{m}_x\nonumber \\
=&i\epsilon \sum_{z=1}^{x}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:1}\left(\mathrm{tr}_{\text{m}}^{x-1:0} \left( \mathcal{U}_x\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}  \right.\right.\nonumber\\
&\cdot \left[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1\mathcal{U}_{z-1}\dots \mathcal{U}_1 (\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{1:x}\bra{0\dots 0}) \U_1^\dagger \dots \U_{z-1}^\dagger U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\right]\nonumber\\
&\cdot\left.\left.U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger\right)\right)+\mathcal{O}(\epsilon^2). \nonumber
\end{align}

As \(\U_x...U^1_{1\, z}\) acts as the identity on  \(\h_{z-1:1}^{\text{in,}1:L\text{,out}}\) and \(\h_{z-2:0}^\text{m}\), we can write

\begin{align}
\Delta \rho_x^m=&i\epsilon \sum_{z=1}^{x}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:z}\Big(\mathrm{tr}_{\text{m}}^{x-1:z-1} \Big( \mathcal{U}_x\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Big[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1  \nonumber\\
&\cdot \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{z-1:1}\left(\mathrm{tr}_{\text{m}}^{z-2:0} \left( \mathcal{U}_{z-1}\dots \mathcal{U}_1 \left(\bigotimes_{y=x}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{1:x}\bra{0\dots 0}\right) \U_1^\dagger \dots \U_{z-1}^\dagger\right)\right)\nonumber\\
&\cdot U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Big]U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger\Big)\Big)+\mathcal{O}(\epsilon^2). \nonumber\\
=&i\epsilon \sum_{z=1}^{x}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:z}\Bigg(\mathrm{tr}_{\text{m}}^{x-1:z-1} \Bigg( \mathcal{U}_x\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x}^z\rho_y^{\text{in}} \right. \nonumber\\
&\otimes \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{z-1:1}  \left(\mathrm{tr}_{\text{m}}^{z-2:0} \left( \mathcal{U}_{z-1}\dots \mathcal{U}_1 \left(\bigotimes_{y=z-1}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{1:z-1}\bra{0\dots 0}\right) \U_1^\dagger \dots \U_{z-1}^\dagger\right)\right)  \nonumber\\
&\otimes \ket{0\dots 0}_{ {1:L+1}}^{z:x}\bra{0\dots 0}\Bigg) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger\Bigg)\Bigg)+\mathcal{O}(\epsilon^2). \nonumber
\end{align}

Actually, it is

\[\rho_{z-1}^\text{m}= \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{z-1:1}  \left(\mathrm{tr}_{\text{m}}^{z-2:0} \left( \mathcal{U}_{z-1}\dots \mathcal{U}_1 \left(\bigotimes_{y=z-1}^1\rho_y^{\text{in}}\otimes \rho_{0}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^{1:z-1}\bra{0\dots 0}\right) \U_1^\dagger \dots \U_{z-1}^\dagger\right)\right), \]

\begin{align}
\Delta \rho_x^m=&i\epsilon \sum_{z=1}^{x}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:z}\Bigg(\mathrm{tr}_{\text{m}}^{x-1:z-1} \Bigg( \mathcal{U}_x\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x}^z\rho_y^{\text{in}} \right. \nonumber\\
&\otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ {1:L+1}}^{z:x}\bra{0\dots 0}\Bigg) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger\Bigg)\Bigg)+\mathcal{O}(\epsilon^2). \label{eq:e}
\end{align}

By keeping in mind, that the the unitaries \(U_j^l\) appear directly in the cost function (equation \(\eqref{eq:c}\)) and indirectly in the state \(\rho_{x-1}^{\text{m}}\), we analogously get

\begin{align*}
\Delta C=& C^\prime-C\\
=& \frac{1}{N}\sum_{x=1}^{N}\Bigg[ i\epsilon\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr} \left( \left( \id_x^{{\text{in,}1:L\text{,m}}} \otimes \id_{x-1}^\text{m} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\left[ K_{j\,x}^l,\right.\right.\\
&\left. \left. U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \right)\\
&+\mathrm{tr}\left((\id_x^{\text{in,}1:L\text{,m}}\otimes \id_{x-1}^\text{m}\otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x(\rho_x^{\text{in}}\otimes \Delta \rho_{x-1}^{\text{m}} \otimes \ket{0\dots 0}_{ 1:L+1}^x\bra{0\dots 0})\mathcal{U}_x^\dagger\right) \Bigg]+\mathcal{O}(\epsilon^2).
\end{align*}

With equation \(\eqref{eq:e}\) it follows that

\begin{align*}
\Delta C=& \frac{i\epsilon}{N}\sum_{x=1}^{N}\Bigg[ \sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr} \left( \left( \id_x^{{\text{in,}1:L\text{,m}}} \otimes \id_{x-1}^\text{m} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\left[ K_{j\,x}^l,\right.\right.\\
&\left. \left. U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \right)\\
&+\mathrm{tr}\left((\id_x^{\text{in,}1:L\text{,m}}\otimes \id_{x-1}^\text{m}\otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x\Bigg(\rho_x^{\text{in}}\otimes \Bigg(\sum_{z=1}^{x-1}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x-1:z}\Bigg(\mathrm{tr}_{\text{m}}^{x-2:z-1} \Bigg( \mathcal{U}_{x-1}\dots \mathcal{U}_{z+1}\right. \nonumber\\
&\cdot U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x-1}^z\rho_y^{\text{in}} \otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ 1:L+1}^{z:x-1}\bra{0\dots 0}\right) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]\\
&\cdot U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x-1}^\dagger\Bigg)\Bigg)\Bigg) \otimes \ket{0\dots 0}_{ 1:L+1}^x\bra{0\dots 0}\Bigg)\mathcal{U}_x^\dagger \Bigg) \Bigg]+\mathcal{O}(\epsilon^2).
\end{align*}

Hence, we get up to first order

\begin{align*}
\delta C=& \lim_{\epsilon\to 0}\frac{\Delta C}{\epsilon}\\
=& \frac{i}{N}\sum_{x=1}^{N}\Bigg[ \sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr} \left( \left( \id_x^{{\text{in,}1:L\text{,m}}} \otimes \id_{x-1}^\text{m} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\left[ K_{j\,x}^l,\right.\right.\\
&\left. \left. U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \right)\\
&+\sum_{z=1}^{x-1}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\mathrm{tr}\left((\id_{x:z}^{\text{in,}1:L\text{,m}}\otimes \id_{z-1}^\text{m}\otimes \id_{x-1:z}^\text{out}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x\Bigg(\rho_x^{\text{in}}\otimes \Bigg( \mathcal{U}_{x-1}\dots \mathcal{U}_{z+1}\right. \nonumber\\
&\cdot U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x-1}^z\rho_y^{\text{in}} \otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ 1:L+1}^{z:x-1}\bra{0\dots 0}\right) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]\\
&\cdot U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x-1}^\dagger\Bigg)\Bigg)\Bigg) \otimes \ket{0\dots 0}_{ 1:L+1}^x\bra{0\dots 0}\Bigg)\mathcal{U}_x^\dagger \Bigg) \Bigg]\\
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left( \left( \id_x^{{\text{in,}1:L\text{,m}}} \otimes \id_{x-1}^\text{m} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\left[ K_{j\,x}^l,\right.\right.\\
&\left. \left. U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \right)\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\Bigg((\id_{x:z}^{\text{in,}1:L\text{,m}}\otimes \id_{z-1}^\text{m}\otimes \id_{x-1:z}^\text{out}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x\dots \mathcal{U}_{z+1}U_{m_{L+1}\, z}^{L+1}\dots U_{j+1\, z}^{l} \nonumber\\
&\cdot  \Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x}^z\rho_y^{\text{in}} \otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ 1:L+1}^{z:x}\bra{0\dots 0}\right) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]\\
&\cdot U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \mathcal{U}_x^\dagger \Bigg) \Bigg].
\end{align*}

By using the cyclic rule of trace we get

\begin{align*}
\delta C
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left( U_{x}^{l+1^\dagger}\dots U_{x}^{L+1^\dagger} \left( \id_x^{{\text{in,}1:L\text{,m}}} \otimes \id_{x-1}^\text{m} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right) U_{x}^{L+1}\dots U_{x}^{l+1}   \right.\\
&\left. U_{m_{l}\,x}^{l}\dots U_{j+1\,x}^{l} \left[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{l}\, x}^{l^\dagger} \right)\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\Bigg(U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_x^\dagger(\id_{x:z}^{\text{in,}1:L\text{,m}}\otimes \id_{z-1}^\text{m}\otimes \id_{x-1:z}^\text{out}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1} \nonumber\\
&\cdot U_{m_{l}\, z}^{l}\dots U_{j+1\, z}^{l} \Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x}^z\rho_y^{\text{in}} \otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ 1:L+1}^{z:x}\bra{0\dots 0}\right) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]U_{j+1\, z}^{l^\dagger }\dots U_{m_{l}\, z}^{l^\dagger }  \Bigg) \Bigg].
\end{align*}

\begin{align}
A=&U_{x}^{l+1^\dagger}\dots U_{x}^{L+1^\dagger} \left( \id_x^{l:L} \otimes \id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right) U_{x}^{L+1}\dots U_{x}^{l+1}\label{eq:f} \\
B_x=&U_{m_{l}\,x}^{l}\dots U_{j+1\,x}^{l} \left[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:l}^x \bra{0\dots 0}\right)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{l}\, x}^{l^\dagger}\label{eq:g}\\
C=&U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_x^\dagger(\id_{x:z+1}^{\text{in,}1:L\text{,m}}\otimes \id_{z}^\text{l:L+1}\otimes \id_{x-1:z+1}^\text{out}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\label{eq:h}
\end{align}

it is

\begin{align}
\delta C
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left( \left(\id_x^{0:l-1}\otimes A\right)  \left(B_x\otimes \ket{0\dots 0}_{l+1:L+1}^x\bra{0\dots 0}\right) \right)\nonumber\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\left( \left(\id_z^{0:l-1}\otimes C\right)  \left( \bigotimes_{y=x}^{z+1}\rho_y^\text{in}\otimes  B_z\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\otimes \ket{0\dots 0}_{1:L+1}^{z+1:x}\bra{0\dots 0} \right) \right)\Bigg]\nonumber\\
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left( A  \left(\mathrm{tr}_{0:l-1}^x(B_x)\otimes \ket{0\dots 0}_{l+1:L+1}^x\bra{0\dots 0}\right) \right)\nonumber\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\left( C  \left( \bigotimes_{y=x}^{z+1}\rho_y^\text{in}\otimes \mathrm{tr}_{0:l-1}^z( B_z)\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\otimes \ket{0\dots 0}_{1:L+1}^{z+1:x}\bra{0\dots 0} \right) \right)\Bigg]\nonumber\\
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left( A  \left(\mathrm{tr}_{0:l-1}^x(B_x)\otimes \id^{l+1:L+1}_x\right) \left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L+1}^x\bra{0\dots 0}\right)\right) \nonumber\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\Bigg( C \left(\id^\text{in}_{x:z+1}\otimes  \mathrm{tr}_{0:l-1}^z( B_z)\otimes \id_z^{l+1:L+1}\otimes \id_{z+1:x}^{1:L+1}\right)  \nonumber \\
&\cdot\left( \bigotimes_{y=x}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\otimes \ket{0\dots 0}_{1:L+1}^{z+1:x}\bra{0\dots 0} \right) \Bigg)\Bigg]\nonumber\\
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L+1}^x\bra{0\dots 0}\right) A  \left(\mathrm{tr}_{0:l-1}^x(B_x)\otimes \id^{l+1:L+1}_x\right) \right) \nonumber\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\Bigg( \left( \bigotimes_{y=x}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\otimes \ket{0\dots 0}_{1:L+1}^{z+1:x}\bra{0\dots 0} \right) C   \nonumber \\
&\cdot \left(\id^\text{in}_{x:z+1}\otimes  \mathrm{tr}_{0:l-1}^z( B_z)\otimes \id_z^{l+1:L+1}\otimes \id_{z+1:x}^{1:L+1}\right)\Bigg)\Bigg]\nonumber
\end{align}

which leads to

\begin{align}
\delta C
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left(\mathrm{tr}_{l+1:L+1}^x\left(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L+1}^x\bra{0\dots 0}\right) A \right) \mathrm{tr}_{0:l-1}^x(B_x) \right) \nonumber\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\Bigg(\mathrm{tr}_\text{in}^{x:z+1} \Bigg(\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{z+1:x}_{1:L+1}\Bigg(\Bigg( \bigotimes_{y=x}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \\
&\otimes \ket{0\dots 0}_{1:L+1}^{z+1:x}\bra{0\dots 0} \Bigg) C  \Bigg)\Bigg)\Bigg)  \mathrm{tr}_{0:l-1}^z( B_z)\Bigg)\Bigg].
\end{align}

These partial traces now have to be expressed in terms of layer-to-layer channels. With (<ref>) it follows that

\begin{align*}
\mathrm{tr}_{0:l-1}^x&(B_x) \\
=&\mathrm{tr}_{0:l-1}^x \Big(U_{m_{l}\,x}^{l}\dots U_{j+1\,x}^{l} \Big[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^1 \Big( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:l}^x \bra{0\dots 0}\Big)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\Big] U_{j+1\, x}^{l^\dagger} \dots U_{m_{l}\, x}^{l^\dagger}\Big)\\
=&\mathrm{tr}_{l-1}^x \Big(U_{m_{l}\,x}^{l}\dots U_{j+1\,x}^{l} \Big[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^l \Big( \mathrm{tr}_{0:l-2}^x \Big(U_{x}^{l-1}\dots  U_{x}^1 \Big( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:l-1}^x \bra{0\dots 0}\Big)U_{x}^{1^\dagger} \dots U_{ x}^{l-1^\dagger}\Big)\\
&\otimes \ket{0\dots 0}_l^x \bra{0\dots 0} \Big) U_{1\, x}^{l^\dagger} \dots U_{j\, x}^{l^\dagger}\Big] U_{j+1\, x}^{l^\dagger} \dots U_{m_{l}\, x}^{l^\dagger}\Big)
\end{align*}

With the layer to layer channel

\begin{align}
\mathcal{E}^l:B(\h^{l-1})&\to  \B(\h^{l})\\
X^{l-1}&\mapsto\mathrm{tr}_{l-1}\left(U^l\left(X^{l-1}\otimes \ket{0\dots 0}_l\bra{0\dots 0}\right)U^{l^\dagger}\right)\nonumber
\end{align}

it is

\begin{align}
\rho_x^{l-1}&=\mathcal{E}^{l-1}\left(\dots \mathcal{E}^1\left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \right)\dots\right)
\end{align}

as shown in [79],  hence

\begin{align*}
\mathrm{tr}&_{0:l-2}^x \Big(U_{x}^{l-1}\dots  U_{x}^1 \Big( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:l-1}^x \bra{0\dots 0}\Big)U_{x}^{1^\dagger} \dots U_{ x}^{l-1^\dagger}\Big)\\
&= \mathrm{tr}_{1:l-2}^x \Big(U_{x}^{l-1}\dots U_{x}^1 \Big(\mathrm{tr}_{0}^x \Big( U_{x}^1 \Big( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1}^x \bra{0\dots 0}\Big)U_{x}^{1^\dagger}\Big) \otimes \ket{0\dots 0}_\text{2:l-1}^x \bra{0\dots 0}\Big) U_{x}^{2^\dagger }\dots U_{ x}^{l-1^\dagger}\Big)\\
&= \mathrm{tr}_{1:l-2}^x \Big(U_{x}^{l-1}\dots U_{x}^1 \Big(\mathcal{E}^1\Big(\rho_x^\text{in}\otimes \rho_{x-1}^\text{m}\Big) \otimes \ket{0\dots 0}_\text{2:l-1}^x \bra{0\dots 0}\Big) U_{x}^{2^\dagger }\dots U_{ x}^{l-1^\dagger}\Big)\\
&=\mathcal{E}^{l-1}\left(\dots \mathcal{E}^1\left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \right)\dots\right)\\
\end{align*}

Thus we get

\begin{align}
\mathrm{tr}_{0:l-1}^x&(B_x) \nonumber\\
=&\mathrm{tr}_{l-1}^x \Big(U_{m_{l}\,x}^{l}\dots U_{j+1\,x}^{l} \Big[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^l \Big( \rho_x^{l-1}\otimes \ket{0\dots 0}_l^x \bra{0\dots 0} \Big) U_{1\, x}^{l^\dagger} \dots U_{j\, x}^{l^\dagger}\Big] U_{j+1\, x}^{l^\dagger} \dots U_{m_{l}\, x}^{l^\dagger}\Big).
\end{align}

Now let

\[\sigma_x^l:=\mathrm{tr}_{l+1:L+1}^x\left(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L+1}^x\bra{0\dots 0}\right) A \right).\]

With equation \(\eqref{eq:f}\) it follows that

\begin{align*}
\sigma_x^l=&\mathrm{tr}_{l+1:L+1}^x\left(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L+1}^x\bra{0\dots 0}\right) U_{x}^{l+1^\dagger}\dots U_{x}^{L+1^\dagger} \left( \id_x^{l:L} \otimes \id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right) U_{x}^{L+1}\dots U_{x}^{l+1} \right)\\
=&\mathrm{tr}_{l+1:L+1}^x\Big(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L}^x\bra{0\dots 0} \otimes \id_x^{L+1}\right)\left(\id_{x}^{l:L}\otimes\ket{0\dots 0}_{L+1}^x\bra{0\dots 0}\right) \left(U_{x}^{l+1^\dagger}\dots U_{x}^{L^\dagger}\otimes \id_x^{L+1}\right)\\
&\cdot \left(\id_x^{l:L-1}\otimes  U_{x}^{L+1^\dagger}\right) \left( \id_x^{l:L} \otimes \id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right) \left(\id_x^{l:L-1}\otimes  U_{x}^{L+1}\right) \left( U_{x}^{L}\dots U_{x}^{l+1}\otimes \id_x^{L+1}\right)\Big)
\\
=&\mathrm{tr}_{l+1:L+1}^x\Big(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L}^x\bra{0\dots 0} \otimes \id_x^{L+1}\right)\left(U_{x}^{l+1^\dagger}\dots U_{x}^{L^\dagger}\otimes \id_x^{L+1}\right)\left(\id_{x}^{l:L}\otimes\ket{0\dots 0}_{L+1}^x\bra{0\dots 0}\right) \\
&\cdot \left(\id_x^{l:L-1}\otimes  U_{x}^{L+1^\dagger}\right) \left( \id_x^{l:L} \otimes \id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right) \left(\id_x^{l:L-1}\otimes  U_{x}^{L+1}\right) \left( U_{x}^{L}\dots U_{x}^{l+1}\otimes \id_x^{L+1}\right)\Big)\\
=&\mathrm{tr}_{l+1:L}^x\Big(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L}^x\bra{0\dots 0}\right)U_{x}^{l+1^\dagger}\dots U_{x}^{L^\dagger}\Big(\id_x^{l:L-1}\otimes \mathrm{tr}_{L+1}^x\Big(\left(\id_{x}^{L}\otimes\ket{0\dots 0}_{L+1}^x\bra{0\dots 0}\right) \\
&\cdot   U_{x}^{L+1^\dagger} \left( \id_x^{L} \otimes \id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)  U_{x}^{L+1}\Big)\Big)  U_{x}^{L}\dots U_{x}^{l+1}\Big).
\end{align*}

By defining the channel

\[\mathcal{F}^l(X^l)=\mathrm{tr}_l\left( \left(\id^{l-1}\otimes \ket{0\dots 0}_l\bra{0\dots 0}\right)U^{l^\dagger}\left(\id^{l-1}\otimes X^l\right)U^l\right),\]

that is the adjoint channel of \(\mathcal{E}^l\) as shown in [79], actually, we get

\begin{align}
\sigma_x^l=&\mathrm{tr}_{l+1:L}^x\Big(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1:L}^x\bra{0\dots 0}\right)U_{x}^{l+1^\dagger}\dots U_{x}^{L^\dagger}\Big(\id_x^{l:L-1}\otimes \mathcal{F}^{L+1}\left(\id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)\Big) U_{x}^{L}\dots U_{x}^{l+1}\Big)\nonumber \\
=&\mathrm{tr}_{l+1}^x\Big(\left(\id_{x}^l\otimes\ket{0\dots 0}_{l+1}^x\bra{0\dots 0}\right)U_{x}^{l+1^\dagger}\Big(\id_x^{l}\otimes \mathcal{F}^{l+2}\left(\dots\mathcal{F}^{L+1}\left(\id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)\dots\right)\Big)  U_{x}^{l+1}\Big)\nonumber \\
=& \mathcal{F}^{l+1}\left(\dots\mathcal{F}^{L+1}\left(\id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)\dots\right).
\end{align}

Now, the only partial trace left is

\begin{align*}
\omega_{zx}^l:=&\mathrm{tr}_\text{in}^{x:z+1} \Bigg(\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x:z+1}_{1:L+1}\Bigg(\Bigg( \bigotimes_{y=x}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x}\bra{0\dots 0} \Bigg) C  \Bigg)\Bigg)\Bigg)
\end{align*}

With equation \(\eqref{eq:h}\) it follows that

\begin{align*}
\omega_{zx}^l=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg( \bigotimes_{y=x}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x}\bra{0\dots 0} \Bigg)\\
&\cdot U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_x^\dagger(\id_{x:z+1}^{\text{in,}1:L\text{,m}}\otimes \id_{z}^\text{l:L+1}\otimes \id_{x-1:z+1}^\text{out}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}})\mathcal{U}_x\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}
\Bigg)\Bigg)\\
=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg(\id_x^\text{in}\otimes  \bigotimes_{y=x-1}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-1}\bra{0\dots 0}\otimes \id_x^{1:L+1} \Bigg)\\
&\cdot \Big(\rho_x^\text{in}\otimes  \id_{z+1:x-1}^\text{in}\otimes \id^{l:L+1}_z \otimes \id^{1:L+1}_{z+1:x-1}\otimes \ket{0\dots 0}_{1:L+1}^{x}\bra{0\dots 0} \Big)\Big(U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_{x-1}^\dagger\otimes \id_x^{1:L+1}\Big)\\
&\cdot \Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_x^\dagger\Big)\Big(\id_{x:z+1}^{\text{in,}1:L\text{,m}}\otimes \id_{z}^\text{l:L+1}\otimes \id_{x-1:z+1}^\text{out}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\Big)\\
&\cdot \Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_x\Big)\Big(\mathcal{U}_{x-1}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\otimes \id_x^{1:L+1} \Big)
\Bigg)\Bigg)\\
=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg(\id_x^\text{in}\otimes  \bigotimes_{y=x-1}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-1}\bra{0\dots 0}\otimes \id_x^{1:L+1} \Bigg)\\
&\cdot \Big(U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_{x-1}^\dagger\otimes \id_x^{1:L+1}\Big)\Big(\rho_x^\text{in}\otimes  \id_{z+1:x-1}^\text{in}\otimes \id^{l:L+1}_z \otimes \id^{1:L+1}_{z+1:x-1}\otimes \ket{0\dots 0}_{1:L+1}^{x}\bra{0\dots 0} \Big)\\
&\cdot \Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_x^\dagger\Big)\Big(\id_{x:z+1}^{\text{in,}1:L\text{,m}}\otimes \id_{z}^\text{l:L+1}\otimes \id_{x-1:z+1}^\text{out}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\Big)\\
&\cdot \Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_x\Big)\Big(\mathcal{U}_{x-1}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\otimes \id_x^{1:L+1} \Big)
\Bigg)\Bigg)\\
=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x-1:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg( \bigotimes_{y=x-1}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-1}\bra{0\dots 0} \Bigg)\\
&\cdot U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_{x-1}^\dagger\Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes\mathrm{tr}^{x}_{\text{in,}1:L+1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}^{x}\bra{0\dots 0} \Big)\\
&\cdot \mathcal{U}_x^\dagger\Big(\id_{x-1}^\text{m}\otimes \id_{x}^{\text{in,}1:L\text{,m}}\otimes\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\Big)\mathcal{U}_x\Big)\Big)\mathcal{U}_{x-1}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\Bigg)\Bigg).
\end{align*}

\begin{align}
\tilde{\mathcal{S}}_x:\B(\h_x^\text{out})&\to \B(\h_{x-1}^\text{m})\\
X_x^\text{out}&\mapsto\mathrm{tr}^{x}_{\text{in,}1:L+1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}^{x}\bra{0\dots 0} \Big) \mathcal{U}_x^\dagger\Big(\id_{x-1}^\text{m}\otimes \id_{x}^{\text{in,}1:L\text{,m}}\otimes X_x^\text{out}\Big)\mathcal{U}_x\Big)\nonumber
\end{align}

it is

\begin{align*}
\omega_{zx}^l=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x-1:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg( \bigotimes_{y=x-1}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-1}\bra{0\dots 0} \Bigg) U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\\
&\cdot\U_{z+1}^\dagger\dots \mathcal{U}_{x-1}^\dagger\Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\Big)\mathcal{U}_{x-1}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\Bigg)\Bigg)\\
=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x-1:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg( \bigotimes_{y=x-1}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-1}\bra{0\dots 0} \Bigg) U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\\
&\cdot\U_{z+1}^\dagger\dots \mathcal{U}_{x-1}^\dagger\Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\Big)\mathcal{U}_{x-1}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\Bigg)\Bigg)\\
=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x-1:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg(\id_{x-1}^\text{in}\otimes  \bigotimes_{y=x-2}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-2}\bra{0\dots 0}\otimes \id_{x-1}^{1:L+1} \Bigg)\\
&\cdot \Big(\rho_{x-1}^\text{in}\otimes  \id_{z+1:x-2}^\text{in}\otimes\id_z^{l:L+1} \otimes \id^{1:L+1}_{z+1:x-2}\otimes \ket{0\dots 0}_{1:L+1}^{x-1}\bra{0\dots 0} \Big)\Big( U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_{x-2}^\dagger\otimes \id_{x-1}^{\text{in},1:L+1}\Big)\\
&\cdot \Big(\id_z^{l:L+1}\otimes \id_{z+1:x-3}^{\text{in},1:L+1}\otimes \id_{x-2}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_{x-1}^\dagger\Big)\Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\Big)\\
&\cdot\Big(\id_z^{l:L+1}\otimes \id_{z+1:x-3}^{\text{in},1:L+1}\otimes \id_{x-2}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_{x-1}\Big)\Big(\mathcal{U}_{x-2}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\otimes \id_{x-1}^{\text{in},1:L+1}\Big)\Bigg)\Bigg)\\
=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x-1:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg(\id_{x-1}^\text{in}\otimes  \bigotimes_{y=x-2}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-2}\bra{0\dots 0}\otimes \id_{x-1}^{1:L+1} \Bigg)\\
&\cdot \Big( U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\U_{z+1}^\dagger\dots \mathcal{U}_{x-2}^\dagger\otimes \id_{x-1}^{\text{in},1:L+1}\Big)\Big(\rho_{x-1}^\text{in}\otimes  \id_{z+1:x-2}^\text{in}\otimes\id_z^{l:L+1} \otimes \id^{1:L+1}_{z+1:x-2}\otimes \ket{0\dots 0}_{1:L+1}^{x-1}\bra{0\dots 0} \Big)\\
&\cdot \Big(\id_z^{l:L+1}\otimes \id_{z+1:x-3}^{\text{in},1:L+1}\otimes \id_{x-2}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_{x-1}^\dagger\Big)\Big(\id_z^{l:L+1}\otimes\id_{z+1:x-2}^{\text{in},1:L+1}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\Big)\\
&\cdot\Big(\id_z^{l:L+1}\otimes \id_{z+1:x-3}^{\text{in},1:L+1}\otimes \id_{x-2}^{\text{in},1:L,\text{out}}\otimes \mathcal{U}_{x-1}\Big)\Big(\mathcal{U}_{x-2}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1}\otimes \id_{x-1}^{\text{in},1:L+1}\Big)\Bigg)\Bigg)\\
=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x-2:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg( \bigotimes_{y=x-2}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-2}\bra{0\dots 0} \Bigg)U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\\
&\cdot  \U_{z+1}^\dagger\dots \mathcal{U}_{x-2}^\dagger\Big(\id_z^{l:L+1}\otimes \id_{z+1:x-3}^{\text{in},1:L+1}\otimes \id_{x-2}^{\text{in},1:L,\text{out}}\otimes\mathrm{tr}^{x-1}_{\text{in,}1:L+1}\Big(\Big(\rho_{x-1}^\text{in} \otimes \id^{\text{m}}_{x-2}\otimes \ket{0\dots 0}_{1:L+1}^{x-1}\bra{0\dots 0} \Big)\\
&\cdot \mathcal{U}_{x-1}^\dagger\Big(\id_{x-2}^{\text{m}}\otimes \id_{x-1}^{\text{in},1:L,\text{out}}\otimes\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\Big)\mathcal{U}_{x-1}\Big)\Big)\mathcal{U}_{x-2}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1} \Bigg)\Bigg)
\end{align*}

By defining

\begin{align}
\mathcal{S}_x:\B(\h_x^\text{m})&\to \B(\h_{x-1}^\text{m})\\
X_x^\text{m}&\mapsto\mathrm{tr}^{x}_{\text{in,}1:L+1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}^{x}\bra{0\dots 0} \Big) \mathcal{U}_x^\dagger\Big(\id_{x-1}^\text{m}\otimes \id_{x}^{\text{in,}1:L\text{,out}}\otimes X_x^\text{m}\Big)\mathcal{U}_x\Big)\nonumber
\end{align}

we get

\begin{align}
\omega_{zx}^l=&\mathrm{tr}^z_{l+1:L+1}\Bigg(\mathrm{tr}^{x-2:z+1}_{\text{in,}1:L+1}\Bigg(\Bigg( \bigotimes_{y=x-2}^{z+1}\rho_y^\text{in}\otimes\id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\nonumber \otimes \ket{0\dots 0}_{1:L+1}^{z+1:x-2}\bra{0\dots 0} \Bigg)U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger}\\
&\cdot  \U_{z+1}^\dagger\dots \mathcal{U}_{x-2}^\dagger\Big(\id_z^{l:L+1}\otimes \id_{z+1:x-3}^{\text{in},1:L+1}\otimes \id_{x-2}^{\text{in},1:L,\text{out}}\otimes \mathcal{S}_{x-1}\big(\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\big)\Big)\mathcal{U}_{x-2}\dots \mathcal{U}_{z+1} U_{z}^{L+1}\dots U_{z}^{l+1} \Bigg)\Bigg)\nonumber\\
=&\dots \nonumber\\
=&\mathrm{tr}^z_{l+1:L+1}\Big(\Big( \id_z^l\otimes \ket{0\dots 0}_{l+1:L+1}^z\bra{0\dots 0}\Big) U_{z}^{l+1^\dagger}\dots U_{z}^{L+1^\dagger} \Big(\id_z^{l:L,\text{out}}\otimes \mathcal{S}_{z+1}\big(\dots \mathcal{S}_{x-1}\big(\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\big)\dots\big)\Big)\nonumber \\
&\cdot U_{z}^{L+1}\dots U_{z}^{l+1} \Big)\nonumber
\end{align}

Analogously to the calculation for \(\sigma_x^l\) this is just

\begin{align}
\omega_{zx}^l=\mathcal{F}_z^{l+1}\big(\dots \mathcal{F}_z^{L+1}\big(\mathcal{S}_{z+1}\big(\dots \mathcal{S}_{x-1}\big(\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\big)\dots\big)\big)\dots\big).
\end{align}

The channels \(\mathcal{S}\) and \(\tilde{\mathcal{S}}\) also have to be written in terms of layer-to-layer channels. It is

\begin{align}
\mathcal{S}_x(X_x^\text{m})=&\mathrm{tr}^{x}_{\text{in,}1:L+1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}^{x}\bra{0\dots 0} \Big) \mathcal{U}_x^\dagger\Big(\id_{x-1}^\text{m}\otimes \id_{x}^{\text{in,}1:L\text{,out}}\otimes X_x^\text{m}\Big)\mathcal{U}_x\Big)\nonumber\\
=&\mathrm{tr}^{x}_{\text{in,}1:L+1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}^{x}\bra{0\dots 0} \Big) U_x^{1^\dagger}\dots U_x^{L+1^\dagger}\Big(\id_{x-1}^\text{m}\otimes \id_{x}^{\text{in,}1:L\text{,out}}\otimes X_x^\text{m}\Big)U_x^{L+1}\dots U_x^1\Big)\nonumber\\
=&\mathrm{tr}^{x}_{\text{in,}1:L+1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L}^{x}\bra{0\dots 0} \otimes \id_x^{L+1} \Big) \Big(\id_x^{0:L} \otimes \ket{0\dots 0}_{L+1}^{x}\bra{0\dots 0} \Big)\nonumber\\
&\cdot \Big(U_x^{1^\dagger}\dots U_x^{L^\dagger}\otimes \id_x^{L+1}\Big)\Big(\id_x^{0:L-1}\otimes U_x^{L+1^\dagger}\Big)\Big(\id_{x}^{0:L,\text{out}}\otimes X_x^\text{m}\Big) \Big(\id_x^{0:L-1}\otimes U_x^{L+1}\Big) \Big(U_x^{L}\dots U_x^1\otimes \id_x^{L+1}\Big)\Big)\nonumber\\
=&\mathrm{tr}^{x}_{\text{in,}1:L+1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L}^{x}\bra{0\dots 0} \otimes \id_x^{L+1} \Big)\Big(U_x^{1^\dagger}\dots U_x^{L^\dagger}\otimes \id_x^{L+1}\Big) \nonumber\\
&\cdot \Big(\id_x^{0:L} \otimes \ket{0\dots 0}_{L+1}^{x}\bra{0\dots 0} \Big) \Big(\id_x^{0:L-1}\otimes U_x^{L+1^\dagger}\Big)\Big(\id_{x}^{0:L,\text{out}}\otimes X_x^\text{m}\Big)\nonumber\\
&\cdot \Big(\id_x^{0:L-1}\otimes U_x^{L+1}\Big) \Big(U_x^{L}\dots U_x^1\otimes \id_x^{L+1}\Big)\Big)\nonumber\\
=&\mathrm{tr}^{x}_{\text{in,}1:L}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L}^{x}\bra{0\dots 0}  \Big)U_x^{1^\dagger}\dots U_x^{L^\dagger}  \Big(\id_x^{0:L-1}\otimes\mathrm{tr}^{x}_{L+1}\Big(\Big(\id_x^{L} \otimes \ket{0\dots 0}_{L+1}^{x}\bra{0\dots 0} \Big)\nonumber\\
&\cdot  U_x^{L+1^\dagger}\Big(\id_{x}^{L,\text{out}}\otimes X_x^\text{m}\Big) U_x^{L+1}\Big)\Big) U_x^{L}\dots U_x^1\Big)\nonumber\\
=&\mathrm{tr}^{x}_{\text{in,}1:L}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1:L}^{x}\bra{0\dots 0}  \Big)U_x^{1^\dagger}\dots U_x^{L^\dagger}  \Big(\id_x^{0:L-1}\otimes \mathcal{F}^{L+1}\Big(\id_x^\text{out}\otimes X_x^\text{m}\Big)\Big) U_x^{L}\dots U_x^1\Big)\nonumber\\
=&\mathrm{tr}^{x}_{\text{in,}1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1}^{x}\bra{0\dots 0}  \Big)U_x^{1^\dagger} \Big(\id_x^{0}\otimes \mathcal{F}^{2} \Big(\dots\mathcal{F}^{L+1}\Big(\id_x^\text{out}\otimes X_x^\text{m}\Big)\Big)\Big)  U_x^1\Big)\nonumber\\
=& \tilde{\mathcal{F}}^1_x\Big( \mathcal{F}^{2} \Big(\dots\mathcal{F}^{L+1}\Big(\id_x^\text{out}\otimes X_x^\text{m}\Big)\Big)\Big) ,
\end{align}

\begin{align}
\tilde{\mathcal{F}}^1_x:\B(\h_x^1)&\to \B(\h_{x-1}^\text{m})\\
X_x^1&\mapsto \mathrm{tr}^{x}_{\text{in,}1}\Big(\Big(\rho_x^\text{in}\otimes \id^{\text{m}}_{x-1}\otimes \ket{0\dots 0}_{1}^{x}\bra{0\dots 0}  \Big)U_x^{1^\dagger} \Big(\id_x^{0}\otimes X_x^1\Big)  U_x^1\Big).\nonumber
\end{align}

Analogously, we get

\begin{equation}
\tilde{\mathcal{S}}_x(X_x^\text{out})=\tilde{\mathcal{F}}^1_x\Big( \mathcal{F}^{2} \Big(\dots\mathcal{F}^{L+1}\Big(\id_x^\text{m}\otimes X_x^\text{out}\Big)\Big)\Big).
\end{equation}

In fact, \(\tilde{\mathcal{S}}_x\) is the adjoint channel to the channel

\begin{align*}
\mathcal{T}_x:\B(\h^\text{m})&\to \B(\h^\text{out})\\
X_{x-1}^\text{m}&\mapsto \mathrm{tr}^x_{0:L,\text{m}} \left( \U_{x}\left(  \rho_x^\text{in}\otimes X^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right) \U_{x}\right)
\end{align*}

as it is

\begin{align*}
\mathrm{tr}(X_x^{\text{out}^\dagger}\mathcal{T}_x(X_{x-1}^\text{m}))&=\mathrm{tr}\left(X_x^{\text{out}^\dagger}\mathrm{tr}^x_{0:L,\text{m}} \left( \U_{x}\left(  \rho_x^\text{in}\otimes X^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right) \U_{x}^\dagger\right)\right)\\
&=\sum_{k}\mathrm{tr}\left(X_x^{\text{out}^\dagger}  \left(\bra{k}\otimes \id^\text{out}\right) \U_{x}\left( \rho_x^\text{in}\otimes X^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right) \U_{x}^\dagger\left(\ket{k}\otimes \id^\text{out}\right)\right)\\
&=\sum_{k}\mathrm{tr}\left(\U_{x}^\dagger\left(\ket{k}\otimes \id^\text{out}\right)X_x^{\text{out}^\dagger}  \left(\bra{k}\otimes \id^\text{out}\right) \U_{x}\left( \rho_x^\text{in}\otimes X^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right) \right)\\
&=\mathrm{tr}\left(\U_{x}^\dagger\left(\id_x^{0:L,\text{m}}\otimes X_x^{\text{out}^\dagger}\right)  \U_{x}\left( \rho_x^\text{in}\otimes X^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right) \right)\\
&=\mathrm{tr}\left(\U_{x}^\dagger\left(\id_x^{0:L,\text{m}}\otimes X_x^{\text{out}^\dagger}\right)  \U_{x}\left(\id_x^{\text{in},1:L+1}\otimes X^\text{m}_{x-1}\right) \left( \rho_x^\text{in}\otimes \id^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right)\right)\\
&=\mathrm{tr}\left(\left( \rho_x^\text{in}\otimes \id^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right)\U_{x}^\dagger\left(\id_x^{0:L,\text{m}}\otimes X_x^{\text{out}^\dagger}\right)  \U_{x}\left(\id_x^{\text{in},1:L+1}\otimes X^\text{m}_{x-1}\right) \right)\\
&=\mathrm{tr}\left(\mathrm{tr}^x_{\text{in},1:L+1}\left(\left( \rho_x^\text{in}\otimes \id^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right)\U_{x}^\dagger\left(\id_x^{0:L,\text{m}}\otimes X_x^{\text{out}^\dagger}\right)  \U_{x}\right) X^\text{m}_{x-1} \right)\\
&=\mathrm{tr}\left(\left(\mathrm{tr}^x_{\text{in},1:L+1}\left(\U_{x}^\dagger\left(\id_x^{0:L,\text{m}}\otimes X_x^{\text{out}}\right)  \U_{x}\left( \rho_x^\text{in}\otimes \id^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right)\right)\right)^\dagger X^\text{m}_{x-1} \right)\\
&=\mathrm{tr}\left(\left(\mathrm{tr}^x_{\text{in},1:L+1}\left(\left( \rho_x^\text{in}\otimes \id^\text{m}_{x-1}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0} \right)\U_{x}^\dagger\left(\id_x^{0:L,\text{m}}\otimes X_x^{\text{out}}\right)  \U_{x}\right)\right)^\dagger X^\text{m}_{x-1} \right)\\
&=\mathrm{tr}\left(\left(\tilde{\mathcal{S}}_x(X_x^\text{out})\right)^\dagger X^\text{m}_{x-1} \right)
\end{align*}

for any orthonormal basis \(\{\ket{k}\}_k\) of \(\h^{0:L,\text{m}}\) and any \(X_x^\text{out}\in \B(\h_x^\text{out}),\ X_{x-1}^\text{m}\in \B(\h_{x-1}^\text{m})\).

Together, we get

\begin{align*}
\delta C
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \Big(\sigma_x^l \mathrm{tr}_{l-1}^x \Big(U_{m_{l}\,x}^{l}\dots U_{j+1\,x}^{l} \Big[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^l \Big( \rho_x^{l-1}\otimes \ket{0\dots 0}_l^x \bra{0\dots 0} \Big) U_{1\, x}^{l^\dagger} \dots U_{j\, x}^{l^\dagger}\Big]\\
&\cdot U_{j+1\, x}^{l^\dagger} \dots U_{m_{l}\, x}^{l^\dagger}\Big) \Big) +\sum_{z=1}^{x-1}\mathrm{tr}\Bigg(\omega_{zx}^l \mathrm{tr}_{l-1}^z \Big(U_{m_{l}\,z}^{l}\dots U_{j+1\,z}^{l} \Big[ K_{j\,z}^l, U_{j\,z}^l\dots U_{1\,z}^l \Big( \rho_z^{l-1}\otimes \ket{0\dots 0}_l^z \bra{0\dots 0} \Big) U_{1\, z}^{l^\dagger} \dots U_{j\, z}^{l^\dagger}\Big]\\
&\cdot  U_{j+1\, z}^{l^\dagger} \dots U_{m_{l}\, z}^{l^\dagger}\Big)\Bigg)\Bigg].
\end{align*}

As \(\h_z^\text{m}=\h_x^\text{m}\), \(\h_z^\text{in}=\h_x^\text{in}\), \(\h_z^\text{out}=\h_x^\text{out}\) and \(U_{j\,z}^l=U_{j\,x}^l\) for all \(x,z\in \{1,...,N\}\) and we only operate on \(x\) or \(z\), respectively, we can drop the subscripts \(x\) and \(z\) on the unitaries and traces. By also pulling out the traces \(\mathrm{tr}_{l-1}\) we get

\begin{align*}
\delta C
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \Big(\Big(\id^{l-1}\otimes \sigma_x^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l} \Big[ K_{j}^l, U_{j}^l\dots U_{1}^l \Big( \rho_x^{l-1}\otimes \ket{0\dots 0}_l\bra{0\dots 0} \Big) U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger}\Big]\\
&\cdot U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger}\Big)  +\sum_{z=1}^{x-1}\mathrm{tr}\Big(\Big(\id^{l-1}\otimes\omega_{zx}^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l} \Big[ K_{j}^l, U_{j}^l\dots U_{1}^l \Big( \rho_z^{l-1}\otimes \ket{0\dots 0}_l \bra{0\dots 0} \Big)\\
&\cdot  U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger}\Big] U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger}\Big)\Bigg].
\end{align*}

As it is

\begin{align*}
\mathrm{tr}(A[B,C]D)=\mathrm{tr}(ABCD-ACBD)=\mathrm{tr}(CDAB-DACB)=\mathrm{tr}([C,DA]B)
\end{align*}

for all \(A,B,C,D\in \B(\h)\) for some Hilbert space \(\h\), it is

\begin{align*}
\delta C
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \Big( \Big[ U_{j}^l\dots U_{1}^l \Big( \rho_x^{l-1}\otimes \ket{0\dots 0}_l\bra{0\dots 0} \Big) U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger}, U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger} \Big(\id^{l-1}\otimes \sigma_x^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l}  \Big]K_{j}^l \Big)\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\Big(\Big[U_{j}^l\dots U_{1}^l \Big( \rho_z^{l-1}\otimes \ket{0\dots 0}_l \bra{0\dots 0} \Big) U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger},U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger} \Big(\id^{l-1}\otimes\omega_{zx}^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l} \Big]K_{j}^l\Big)\Bigg]\\
=&\frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}\Big(M_{j\,x}^lK_j^l\Big)
\end{align*}

\begin{align*}
M_{j\,x}^l=&  \Big[ U_{j}^l\dots U_{1}^l \Big( \rho_x^{l-1}\otimes \ket{0\dots 0}_l\bra{0\dots 0} \Big) U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger}, U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger} \Big(\id^{l-1}\otimes \sigma_x^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l}  \Big]\\
&+\sum_{z=1}^{x-1}\Big[U_{j}^l\dots U_{1}^l \Big( \rho_z^{l-1}\otimes \ket{0\dots 0}_l \bra{0\dots 0} \Big) U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger},U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger} \Big(\id^{l-1}\otimes\omega_{zx}^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l} \Big]
\end{align*}

\begin{align*}
\rho_x^{l-1}&=\mathcal{E}^{l-1}\left(\dots \mathcal{E}^1\left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \right)\dots\right),\\
\sigma_x^l&= \mathcal{F}^{l+1}\left(\dots\mathcal{F}^{L+1}\left(\id_{x}^{\text{m}} \otimes \ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\right)\dots\right),\\
\omega_{zx}^l&=\mathcal{F}_z^{l+1}\big(\dots \mathcal{F}_z^{L+1}\big(\mathcal{S}_{z+1}\big(\dots \mathcal{S}_{x-1}\big(\tilde{\mathcal{S}}_x\big(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\big)\big)\dots\big)\big)\dots\big).
\end{align*}

In [79] it was shown that the matrices \(K_j^l\) that maximise a derivative of the cost function in the form

\begin{align*}
\delta C=\frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}\Big(M_{j\,x}^lK_j^l\Big)
\end{align*}

are the matrices

\begin{equation}
K_j^l=\frac{i2^{m_{l-1}}\eta}{N}\sum_{x=1}^N \mathrm{tr}_\text{rest}\Big(M_{j\,x}^l\Big),
\end{equation}

if one demands, that the sum over the matrix elements of \(K_j^l\) should be smaller than some constant. Thereby, \(\mathrm{tr}_\text{rest}\) denotes the trace over all qubits not affected by \(U_j^l\) and \(\eta\) the learning rate. Altogether, this leads to the algorithm in Box 1:

Box 1: Training algorithm local cost with pure output for \(M=1\)

1. Initialize:

Choose the perceptron unitaries \(U_j^l\) randomly.

2. Feedforward:
For \(x=1,...,N\) do the following:

2.1. Set
\(\rho_x^0=\rho_x^\text{in}\otimes\rho_{x-1}^\text{m}\) where \(\rho^\text{m}_0\) is given.

2.2. For \(l=1,...,L+1\) set
\(\rho_x^l=\mathrm{tr}_{l-1}\left( U^l(\rho_x^{l-1}\otimes \ket{0\dots 0}_{l}\bra{0\dots 0}) U^{l^\dagger} \right).\)

2.3. Set
\(\rho_x^\text{m}=\mathrm{tr}_\text{out}^x(\rho_x^{L+1})\)
\(\rho_x^\text{out}=\mathrm{tr}_\text{m}^x(\rho_x^{L+1}).\)

3. Feedbackward: For \(x=N,...,1\) do the following:

3.1. Set \(\sigma_x^{L+1}=\id ^\text{m}_x\otimes \ket{\phi_x^{\text{out}}}\bra{\phi_x^{\text{out}}}\).

3.2. For \(l=L,...,1\) set
\(\sigma_x^l=\mathrm{tr}_{l+1}\left( (\id ^l\otimes \ket{0\dots 0}_{l+1}\bra{0\dots 0}) U^{l+1^\dagger} (\id ^l\otimes \sigma_x^{l+1})  U^{l+1} \right)\).

3.3. Set
\(\omega_{xx}^\text{m}=\mathrm{tr}_{1,\text{in}}\left( (\rho_x^\text{in}\otimes \id _{x-1}^m\otimes \ket{0\dots 0}_{1}\bra{0\dots 0}) U^{1^\dagger} (\id _x^0\otimes \sigma_x^{1})  U^{1} \right)\).

3.4. For \(z=x-1,...,1\) do the following steps:

3.4.1. Set \(\omega_{zx}^{L+1}=\omega_{z+1\,x}^\text{m}\otimes \id ^\text{out}\).

3.4.2. For \(l=L,...,1\) set
\(\omega_{zx}^l=\mathrm{tr}_{l+1}\left( (\id ^l\otimes \ket{0\dots 0}_{l+1}\bra{0\dots 0}) U^{l+1^\dagger}(\id ^l\otimes \omega_{zx}^{l+1})   U^{l+1} \right).\)

3.4.3. Set
\(\omega_{zx}^\text{m}=\mathrm{tr}_{1,\text{in}}\left( (\rho_x^\text{in}\otimes \id _{x-1}^m\otimes \ket{0\dots 0}_{1}\bra{0\dots 0}) U^{1^\dagger}  (\id _x^0\otimes \omega_{zx}^{1})  U^{1} \right).\)

4. Update the network:  For \(l=1,...,L+1\) and \(j=1,...,m_l\) do the following:

4.1. For \(x=1,...,N\) set

\begin{align*}
M_{j\,x}^l=& \Bigg[\prod_{k=j}^{1} U_k^l\left( \rho_x^{l-1}\otimes\ket{0\dots 0}_l\bra{0\dots 0}\right) \prod_{k=1}^{j} U_k^{l^\dagger},\\
&\ \prod_{k=j+1}^{m_l} U_k^{l^\dagger} \left(\id ^{l-1}\otimes \sigma_x^l\right)  \prod_{k=m_l}^{j+1} U_k^l \Bigg]\\
&+ \sum_{z=1}^{x-1} \Bigg[ \prod_{k=j}^{1} U_k^l \left( \rho_z^{l-1}\otimes\ket{0\dots 0}_l\bra{0\dots 0}\right)  \prod_{k=1}^{j} U_k^{l^\dagger} ,\\
&\ \prod_{k=j+1}^{m_l} U_k^{l^\dagger} \left(\id ^{l-1}\otimes \omega_{zx}^l\right) \prod_{k=m_l}^{j+1} U_k^l \Bigg].
\end{align*}

4.2. Set

\[K_j^l=\frac{i2^{m_{l-1}} \eta}{N}\sum_{x=1}^N \mathrm{tr}_\text{rest}(M_{j\,x}^l),\]

where \(\mathrm{tr}_\text{rest}\) denotes that the trace is taken over all systems that are unaffected by \(U_j^l\).

4.3. Update each unitary $U_j^l$ according to $U_j^l\rightarrow e^{i\epsilon K_j^l} U_j^l$.

5. Repeat: Repeat steps 2. to 4. until the cost function reaches its maximum.

§ OPTIMISING THE LOCAL COST WITH MIXED OUTPUT STATES CLASSICALLY

If we want to learn from mixed output states, i.e. a training set of the form

\[S=\left(\rho_0^\text{m},\left(\rho_1^\text{in},\sigma_1^\text{out}\right),\dots ,  \left(\rho_N^\text{in},\sigma_N^\text{out}\right)  \right)\]

the first idea would be to use the fidelity for mixed states given by

\[F(\rho,\sigma )=\left(\mathrm{tr}\sqrt{\sqrt{\rho}\sigma \sqrt{\rho}}\right)^2\]

in the optimisation process, but the algorithm is not as easy to get with this, because it is not easy to determine the derivative of this anymore. But we can instead use the Hilbert-Schmidt norm for the cost function such that we get the cost function

\begin{equation}
C=C_\text{local,mixed}=\frac{1}{N}\sum_{x=1}^N \Vert \rho_x^\text{out}- \sigma_x^\text{out} \Vert _2 ^2= \frac{1}{N} \sum_{x=1}^N \mathrm{tr}\left(\rho_x^\text{out}- \sigma_x^\text{out}\right)^2.
\end{equation}

Because of the cyclic rule of trace, this is just

\[ 	C= \frac{1}{N} \sum_{x=1}^N \mathrm{tr}\left(\rho_x^{\text{out}^2}+ \sigma_x^{\text{out}^2} - 2 \sigma_x^\text{out}\rho_x^\text{out}\right). \]

We now want to get a similar algorithm to the one in Box 1. As before, we start the training by choosing the perceptron unitaries \(U_j^l\) randomly at first and updating the perceptron unitaries in each step according to

\[U_j^l\mapsto U_j^{l\prime}=e^{i\epsilon K_j^l}U_j^l,\]

where the \(K_j^l\) are chosen in a way so that they minimise the derivative of the cost function, as opposed to the case with pure output states as the norm distance has to be minimised.

Up to first order in \(\epsilon\) it is because of the cyclic rule of trace

\[ \Delta	C= \frac{1}{N} \sum_{x=1}^N 2\mathrm{tr}\left(\rho_x^{\text{out}}\Delta \rho_x^\text{out}\right) - 2\mathrm{tr}\left( \sigma_x^\text{out}\Delta\rho_x^\text{out}\right) +\mathcal{O}(\epsilon ^2) =\frac{2}{N} \sum_{x=1}^N \mathrm{tr}\left(\left(\rho_x^{\text{out}}- \sigma_x^\text{out}\right)\Delta\rho_x^\text{out}\right)+\mathcal{O}(\epsilon ^2). \]

Like before, it is

\[\rho_x^\text{out}=\mathrm{tr}_{\text{in,}1:L\text{,m}}^x(\mathrm{tr}_{\text{m}}^{x-1}(U^{L+1}_{m_{L+1}} \dots U^{l}_{j} \dots U^{1}_{1} (\rho_x^{\text{in}}\otimes \rho_{x-1}^{\text{m}} \otimes \ket{0\dots 0}_{ {1:L+1}}^x\bra{0\dots 0}) U^{1^\dagger}_{1} \dots U^{l^\dagger}_{j} \dots U^{L+1^\dagger}_{m_{L+1}}))\]

and hence

\begin{align*}
\Delta  \rho_x^\text{out}
=&  i\epsilon\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{\text{in,}1:L\text{,m}}^x\Big(\mathrm{tr}_{\text{m}}^{x-1}\Big( U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\Big[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^1 \big( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \\
& \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\big) U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\Big] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \Big) \\
&+\mathcal{U}_x(\rho_x^{\text{in}}\otimes \Delta \rho_{x-1}^{\text{m}} \otimes \ket{0\dots 0}_{ 1:L+1}^x\bra{0\dots 0})\mathcal{U}_x^\dagger +\mathcal{O}(\epsilon^2)\Big)\Big) .
\end{align*}

We already know, that

\begin{align*}
\Delta \rho_x^m=&i\epsilon \sum_{z=1}^{x}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x:z}\Bigg(\mathrm{tr}_{\text{m}}^{x-1:z-1} \Bigg( \mathcal{U}_x\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x}^z\rho_y^{\text{in}} \right.\\
&\otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ {1:L+1}}^{z:x}\bra{0\dots 0}\Bigg) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger\Bigg)\Bigg)+\mathcal{O}(\epsilon^2),
\end{align*}

\begin{align*}
\Delta  \rho_x^\text{out}
=&  i\epsilon\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \Bigg[ \mathrm{tr}_{\text{in,}1:L\text{,m}}^x\Big(\mathrm{tr}_{\text{m}}^{x-1}\Big( U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\left[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \right. \right.\\
& \left. \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right) \left. \left. U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \right) \Big)\Big) + \sum_{z=1}^x \mathrm{tr}_{\text{in,}1:L\text{,m}}^x\Bigg(\mathrm{tr}_{\text{m}}^{x-1}\Bigg( \mathcal{U}_x\\
&\cdot\Bigg(\rho_x^{\text{in}}\otimes \mathrm{tr}_{{\text{in,}1:L\text{,out}}}^{x-1:z}\Bigg(\mathrm{tr}_{\text{m}}^{x-2:z-1} \Bigg( \mathcal{U}_{x-1}\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1  \\
&\cdot   \left( \bigotimes_{y=x-1}^z\rho_y^{\text{in}} \otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ {1:L+1}}^{z:x-1}\bra{0\dots 0}\right) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x-1}^\dagger\Bigg)\Bigg) \\
& \otimes \ket{0\dots 0}_{ 1:L+1}^{x}\bra{0\dots 0}\Bigg)\mathcal{U}_x^\dagger \Bigg) \Bigg) \Bigg] +\mathcal{O}(\epsilon^2) .\\
=&  i\epsilon\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \Bigg[ \mathrm{tr}_{\text{in,}1:L\text{,m}}^{x}\Big(\mathrm{tr}_{\text{m}}^{x-1}\Big( U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\left[ K_{j\,x}^l, U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \right. \right.\\
& \left. \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right) \left. \left. U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \right) \Big)\Big) + \sum_{z=1}^x \mathrm{tr}_{\text{in,}1:L\text{,m}}^{x:z}\Bigg(\mathrm{tr}_{\text{m}}^{x-1:z-1} \\
&\cdot \Bigg( \mathcal{U}_{x}\dots \mathcal{U}_{z+1} U_{m_{L+1}\, z}^{L+1} \dots U_{j+1\, z}^{l}\Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1 \left( \bigotimes_{y=x}^z\rho_y^{\text{in}} \otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ {1:L+1}}^{z:x}\bra{0\dots 0}\right) \\
&\cdot    U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \U_{x}^\dagger\Bigg)\Bigg)   \Bigg) \Bigg] +\mathcal{O}(\epsilon^2) .
\end{align*}

Therefore, we get up to first order

\begin{align*}
\delta C
=& \frac{i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l}\Bigg[  \mathrm{tr} \left( \left( \id_x^{{\text{in,}1:L\text{,m}}} \otimes \id_{x-1}^\text{m} \otimes \left(\rho_x^{\text{out}}- \sigma_x^\text{out}\right)\right)U_{m_{L+1}\,x}^{L+1}\dots U_{j+1\,x}^{l}\left[ K_{j\,x}^l,\right.\right.\\
&\left. \left. U_{j\,x}^l\dots U_{1\,x}^1 \left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \otimes \ket{0\dots 0}_\text{1:L+1}^x \bra{0\dots 0}\right)U_{1\, x}^{1^\dagger}\dots U_{j\, x}^{l^\dagger}\right] U_{j+1\, x}^{l^\dagger} \dots U_{m_{L+1}\, x}^{L+1^\dagger} \right)\\
&+\sum_{z=1}^{x-1}\mathrm{tr}\Bigg(\left(\id_{x:z}^{\text{in,}1:L\text{,m}}\otimes \id_{z-1}^\text{m}\otimes \id_{x-1:z}^\text{out}\otimes\left(\rho_x^{\text{out}}- \sigma_x^\text{out}\right)\right)\mathcal{U}_x\dots \mathcal{U}_{z+1}U_{m_{L+1}\, z}^{L+1}\dots U_{j+1\, z}^{l} \nonumber\\
&\cdot  \Bigg[ K_{j\, z}^l, U_{j\, z}^{l} \dots U_{1\,z}^1   \left( \bigotimes_{y=x}^z\rho_y^{\text{in}} \otimes \rho_{z-1}^\text{m}  \otimes \ket{0\dots 0}_{ 1:L+1}^{z:x}\bra{0\dots 0}\right) U_{1\, z}^{1^\dagger }\dots U_{j\, z}^{l^\dagger }\Bigg]\\
&\cdot U_{j+1\, z}^{l^\dagger }\dots U_{m_{L+1}\, z}^{L+1^\dagger } \U_{z+1}^\dagger\dots \mathcal{U}_x^\dagger \Bigg) \Bigg].
\end{align*}

We can see, that this has the same form as before with pure states except that we have \(\rho_x^{\text{out}}- \sigma_x^\text{out}\) instead of \(\ket{\phi_x^\text{out}}\bra{\phi_x^\text{out}}\). The other difference is, that we have to minimise and not maximise the cost function, hence with

\begin{align*}
\chi_x^\text{out}=&\sigma_x^\text{out}-\rho_x^\text{out}\\
\rho_x^{l-1}=&\mathcal{E}^{l-1}\left(\dots \mathcal{E}^1\left( \rho_x^\text{in}\otimes \rho_{x-1}^\text{m} \right)\dots\right)\\
\chi_x^l=& \mathcal{F}^{l+1}\left(\dots\mathcal{F}^{L+1}\left(\id_{x}^{\text{m}} \otimes \chi_x^\text{out}\right)\dots\right)\\
\omega_{zx}^l=&\mathcal{F}_z^{l+1}\big(\dots \mathcal{F}_z^{L+1}\big(\mathcal{S}_{z+1}\big(\dots \mathcal{S}_{x-1}\big(\tilde{\mathcal{S}}_x\big(\chi_x^\text{out}\big)\big)\dots\big)\big)\dots\big)\\
M_{j\,x}^l=&  \Big[ U_{j}^l\dots U_{1}^l \Big( \rho_x^{l-1}\otimes \ket{0\dots 0}_l\bra{0\dots 0} \Big) U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger}, U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger} \Big(\id^{l-1}\otimes \chi_x^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l}  \Big]\\
&+\sum_{z=1}^{x-1}\Big[U_{j}^l\dots U_{1}^l \Big( \rho_z^{l-1}\otimes \ket{0\dots 0}_l \bra{0\dots 0} \Big) U_{1}^{l^\dagger} \dots U_{j}^{l^\dagger},U_{j+1}^{l^\dagger} \dots U_{m_{l}}^{l^\dagger} \Big(\id^{l-1}\otimes\omega_{zx}^l\Big) U_{m_{l}}^{l}\dots U_{j+1}^{l} \Big]
\end{align*}

we get

\begin{align*}
\delta C
=&-\frac{2i}{N}\sum_{x=1}^{N}\sum_{l=1}^{L+1}\sum_{j=1}^{m_l} \mathrm{tr}\Big(M_{j\,x}^lK_j^l\Big).
\end{align*}

The matrices \(K_j^l\) that minimize the derivative of the cost function are then given by

\begin{equation}
K_j^l=\frac{i2^{m_{l-1}+1}\eta}{N}\sum_{x=1}^N \mathrm{tr}_\text{rest}\Big(M_{j\,x}^l\Big).
\end{equation}

Therefore, the optimisation algorithm is the following one.

Box 2: Training algorithm local cost with mixed output for \(M= 1\)

1. Initialize:

Choose the perceptron unitaries \(U_j^l\) randomly.

2. Feedforward:
For \(x=1,...,N\) do the following:

2.1. Set
\(\rho_x^0=\rho_x^\text{in}\otimes\rho_{x-1}^\text{m}\) where \(\rho^\text{m}_0\) is given.

2.2. For \(l=1,...,L+1\) set
\(\rho_x^l=\mathrm{tr}_{l-1}\left( U^l(\rho_x^{l-1}\otimes \ket{0\dots 0}_{l}\bra{0\dots 0}) U^{l^\dagger} \right).\)

2.3. Set
\(\rho_x^\text{m}=\mathrm{tr}_\text{out}^x(\rho_x^{L+1})\)
\(\rho_x^\text{out}=\mathrm{tr}_\text{m}^x(\rho_x^{L+1}).\)

3. Feedbackward: For \(x=N,...,1\) do the following:

3.1. Set \(\chi_x^{L+1}=\id ^\text{m}_x\otimes  \left(\sigma_x^\text{out}- \rho_x^{\text{out}}\right)\).

3.2. For \(l=L,...,1\) set
\(\chi_x^l=\mathrm{tr}_{l+1}\left( (\id ^l\otimes \ket{0\dots 0}_{l+1}\bra{0\dots 0}) U^{l+1^\dagger} (\id ^l\otimes \chi_x^{l+1})  U^{l+1} \right)\).

3.3. Set
\(\omega_{xx}^\text{m}=\mathrm{tr}_{1,\text{in}}\left( (\rho_x^\text{in}\otimes \id _{x-1}^m\otimes \ket{0\dots 0}_{1}\bra{0\dots 0}) U^{1^\dagger} (\id _x^0\otimes \chi_x^{1})  U^{1} \right)\).

3.4. For \(z=x-1,...,1\) do the following steps:

3.4.1. Set \(\omega_{zx}^{L+1}=\omega_{z+1\,x}^\text{m}\otimes \id ^\text{out}\).

3.4.2. For \(l=L,...,1\) set
\(\omega_{zx}^l=\mathrm{tr}_{l+1}\left( (\id ^l\otimes \ket{0\dots 0}_{l+1}\bra{0\dots 0}) U^{l+1^\dagger}(\id ^l\otimes \omega_{zx}^{l+1})   U^{l+1} \right).\)

3.4.3. Set
\(\omega_{zx}^\text{m}=\mathrm{tr}_{1,\text{in}}\left( (\rho_x^\text{in}\otimes \id _{x-1}^m\otimes \ket{0\dots 0}_{1}\bra{0\dots 0}) U^{1^\dagger}  (\id _x^0\otimes \omega_{zx}^{1})  U^{1} \right).\)

4. Update the network:  For \(l=1,...,L+1\) and \(j=1,...,m_l\) do the following:

4.1. For \(x=1,...,N\) set

\begin{align*}
M_{j\,x}^l=& \Bigg[\prod_{k=j}^{1} U_k^l\left( \rho_x^{l-1}\otimes\ket{0\dots 0}_l\bra{0\dots 0}\right) \prod_{k=1}^{j} U_k^{l^\dagger},\\
&\ \prod_{k=j+1}^{m_l} U_k^{l^\dagger} \left(\id ^{l-1}\otimes \chi_x^l\right)  \prod_{k=m_l}^{j+1} U_k^l \Bigg]\\
&+ \sum_{z=1}^{x-1} \Bigg[ \prod_{k=j}^{1} U_k^l \left( \rho_z^{l-1}\otimes\ket{0\dots 0}_l\bra{0\dots 0}\right)  \prod_{k=1}^{j} U_k^{l^\dagger} ,\\
&\ \prod_{k=j+1}^{m_l} U_k^{l^\dagger} \left(\id ^{l-1}\otimes \omega_{zx}^l\right) \prod_{k=m_l}^{j+1} U_k^l \Bigg].
\end{align*}

4.2. Set

\[K_j^l=\frac{i2^{m_{l-1}} \eta}{N}\sum_{x=1}^N \mathrm{tr}_\text{rest}(M_{j\,x}^l),\]

where \(\mathrm{tr}_\text{rest}\) denotes that the trace is taken over all systems that are unaffected by \(U_j^l\).

4.3. Update each unitary $U_j^l$ according to $U_j^l\rightarrow e^{i\epsilon K_j^l} U_j^l$.

5. Repeat: Repeat steps 2. to 4. until the cost function reaches its maximum.

Note that the \(\chi_x^l\) and \(\omega_{zx}^l\) are not normed, so no states, anymore.

§ OPTIMISING THE GLOBAL COST WITH PURE OUTPUT CLASSICALLY

The global cost for pure output states is given by

\begin{align*}
C=C_\text{global,pure}&=\left( \bra{\phi_1^\text{out}}\otimes \dots\otimes \bra{\phi_N^\text{out}} \right) \mathrm{tr}_\text{m}^N \left(\rho^\text{OUT}\right)  \left( \ket{\phi_1^\text{out}}\otimes \dots\otimes \ket{\phi_N^\text{out}} \right)\\
&=\mathrm{tr}\left( \left( \ket{\phi_1^\text{out}}\bra{\phi_1^\text{out}}\otimes \dots \otimes \ket{\phi_N^\text{out}}\bra{\phi_N^\text{out}} \otimes \id_N^\text{m}   \right) \rho^\text{OUT} \right)\\
&=\mathrm{tr}\left( \left( \ket{\phi_1^\text{out}}\bra{\phi_1^\text{out}}\otimes \dots \otimes \ket{\phi_N^\text{out}}\bra{\phi_N^\text{out}} \otimes \id_N^\text{m}   \right) \left(\N_{\U_N}\circ \dots \circ \N_{\U_1}\right) \left(\rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes \rho_N^\text{in}\right) \right).
\end{align*}

With \(\sigma_i^\text{out}= \ket{\phi_i^\text{out}}\bra{\phi_i^\text{out}} \),

\[ \mathcal{N}_{\U}\left( \rho^{0} \right) =\mathrm{tr}^{0:L}(\mathcal{U}(\rho^{0}\otimes \ket{0\dots 0}_{1:L+1}\bra{0\dots 0})\mathcal{U}^\dagger) \]

and hence

\begin{align*}
\left(\N_{\U_N}\circ \dots \circ \N_{\U_1}\right) \left(\rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes \rho_N^\text{in}\right) =\mathrm{tr}^{0:L}_{1:N}(\mathcal{U}_N\dots \U_1( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{1:x}\bra{0\dots 0})\mathcal{U}_1^\dagger \dots \U_N^\dagger )
\end{align*}

we get

\begin{align*}
C=&\mathrm{tr}\left( \left( \id_{1:N}^{0:L} \otimes \sigma_1^\text{out}\otimes \dots \otimes \sigma_N^\text{out} \otimes \id_N^\text{m}   \right) \mathcal{U}_N\dots \U_1( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{1:x}\bra{0\dots 0})\mathcal{U}_1^\dagger \dots \U_N^\dagger \right)\\
=&\mathrm{tr}\Big( \left( \id_{1:N}^{0:L} \otimes \id_1^\text{out} \otimes \sigma_2^\text{out} \otimes \dots \otimes \sigma_N^\text{out} \otimes \id_N^\text{m}   \right)  \left( \id_{1:N}^{0:L} \otimes \sigma_1^\text{out}\otimes \id^\text{out}_{2:N}\otimes \id_N^\text{m}   \right) \left(\U_N\dots \U_2 \otimes \id_{0}^\text{m}\otimes \id_1^{\text{in},1:L,\text{out}}\right) \left(\U_1\otimes \id_{2:N}^{\text{in},1:L+1}\right)\\
&\cdot \left( \rho_0^\text{m}\otimes \rho_1^\text{in} \ket{0\dots 0}_{1:L+1}^{1}\bra{0\dots 0}\otimes \id_{2:N}^{\text{in},1:L+1}\right)  \left( \id_0^\text{m}\otimes \id_1^{\text{in},1:L+1}\otimes \rho_2^\text{in}\otimes \dots \otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{2:x}\bra{0\dots 0}\right) \\
&\cdot \left(\U_1^\dagger\otimes \id_{2:N}^{\text{in},1:L+1}\right)  \left(\U_2^\dagger\dots \U_N^\dagger \otimes \id_{0}^\text{m}\otimes \id_1^{\text{in},1:L,\text{out}}\right) \Big)\\
=& \mathrm{tr}\Big( \left( \id_{1:N}^{0:L} \otimes \id_1^\text{out} \otimes \sigma_2^\text{out} \otimes \dots \otimes \sigma_N^\text{out} \otimes \id_N^\text{m}   \right)   \left(\U_N\dots \U_2 \otimes \id_{0}^\text{m}\otimes \id_1^{\text{in},1:L,\text{out}}\right)  \left( \id_{1:N}^{0:L} \otimes \sigma_1^\text{out}\otimes \id^\text{out}_{2:N}\otimes \id_N^\text{m}   \right) \left(\U_1\otimes \id_{2:N}^{\text{in},1:L+1}\right)  \\
&\cdot \left( \rho_0^\text{m}\otimes \rho_1^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{1}\bra{0\dots 0}\otimes \id_{2:N}^{\text{in},1:L+1}\right) \left(\U_1^\dagger\otimes \id_{2:N}^{\text{in},1:L+1}\right)  \\
&\cdot   \left( \id_0^\text{m}\otimes \id_1^{\text{in},1:L+1}\otimes \rho_2^\text{in}\otimes \dots \otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{2:x}\bra{0\dots 0}\right) \left(\U_2^\dagger\dots \U_N^\dagger \otimes \id_{0}^\text{m}\otimes \id_1^{\text{in},1:L,\text{out}}\right) \Big)\\
=& \mathrm{tr}\Big( \left( \id_{2:N}^{0:L} \otimes \sigma_2^\text{out}\otimes \dots \otimes \sigma_N^\text{out} \otimes \id_N^\text{m}   \right)   \U_N\dots \U_2   \Big( \mathrm{tr}_0^\text{m}\Big( \mathrm{tr}_1^{\text{in}, 1:L,\text{out}}\Big(\left( \sigma_1^\text{out}\otimes \id_{0,1}^{\text{m}} \otimes \id^{\text{in},1:L}_{1}  \right) \U_1  \\
&\cdot \left( \rho_0^\text{m}\otimes \rho_1^\text{in} \otimes\ket{0\dots 0}_{1:L+1}^{1}\bra{0\dots 0}\right)   \U_1^\dagger\Big) \Big) \otimes \id_{2:N}^{\text{in},1:L+1}  \Big) \left( \id_1^\text{m}\otimes  \rho_2^\text{in}\otimes \dots \otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{2:x}\bra{0\dots 0}\right) \U_2^\dagger\dots \U_N^\dagger  \Big).
\end{align*}

By setting \(\tau_0^\text{m}=\rho_0^\text{m}\), \(\tau_x^0=\tau_{x-1}^\text{m}\otimes \rho_x^\text{in}\),

\[\tau_x^l=\mathcal{E}^l\left(\tau_x^{l-1}\right)\]

\begin{align}
\tau_x^\text{m}&=\mathrm{tr}_{x-1}^\text{m}\left(\mathrm{tr}_x^{\text{in},1:L,\text{out}}\left( \left(\sigma_x^\text{out}\otimes \id_{x-1,x}^\text{m}\otimes \id_x^{\text{in},1:L}\right) \U_x \left(\tau_{x-1}^\text{m}\otimes \rho_x^\text{in} \otimes \ket{0...0}_x^{1:L+1}\bra{0...0}\right) \U_x^\dagger\right)\right)\\
&=\mathrm{tr}_x^\text{out}\left(\left(\sigma_x^\text{out}\otimes \id_x^\text{m}\right)\tau_x^{L+1}\right)
\end{align}

for \(x=1,...,N\) we obtain

\begin{align}
C=& \mathrm{tr}\Big( \left( \id_{2:N}^{0:L} \otimes \sigma_2^\text{out}\otimes \dots \otimes \sigma_N^\text{out} \otimes \id_N^\text{m}   \right)   \U_N\dots \U_2   \Big( \tau_1^\text{m} \otimes \id_{2:N}^{\text{in},1:L+1}  \Big) \nonumber \\
&\cdot\left( \id_1^\text{m}\otimes  \rho_2^\text{in}\otimes \dots \otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{2:x}\bra{0\dots 0}\right) \U_2^\dagger\dots \U_N^\dagger  \Big) \nonumber \\
=& \mathrm{tr}\Big( \left( \id_{2:N}^{0:L} \otimes \sigma_2^\text{out}\otimes \dots \otimes \sigma_N^\text{out} \otimes \id_N^\text{m}   \right)   \U_N\dots \U_2   \left( \tau_1^\text{m}\otimes  \rho_2^\text{in}\otimes \dots \otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{1:L+1}^{2:x}\bra{0\dots 0}\right) \U_2^\dagger\dots \U_N^\dagger  \Big) \nonumber \\
=& \dots\nonumber\\
\end{align}

Note that the \(\tau_x^l\) do not have trace 1. By using the same update rule as for the local cost and the same argument we get

\begin{align*}
\Delta \rho^\text{OUT}=& i \epsilon \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{0:N-1}^{\text{m}}\Big(   \mathrm{tr}_{1:N}^{\text{in},1:L}\Big( \U_N \dots \U_{z+1} U^{L+1}_{m_{L+1}\, z}\dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \U_{z-1}\dots \U_1 \\
&\left( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{1:N}^{1:L+1}\bra{0\dots 0} \right)\U_1^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{L+1^\dagger}_{m_{L+1}\, z} \U_{z+1}^\dagger \dots \U_N^\dagger \Big) \Big).
\end{align*}

Therefore we get in the same way as before

\begin{align*}
\delta C=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{0:N}^{\text{m}}\Big(   \mathrm{tr}_{1:N}^{\text{in},1:L,\text{out}}\Big( \U_N \dots \U_{z+1} U^{L+1}_{m_{L+1}\, z}\dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \U_{z-1}\dots \U_1 \\
&\left( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{1:N}^{1:L+1}\bra{0\dots 0} \right)\U_1^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{L+1^\dagger}_{m_{L+1}\, z} \U_{z+1}^\dagger \dots \U_N^\dagger \\
&\left(\sigma_1^\text{out}\otimes\dots\otimes \sigma_N^\text{out} \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L}\right)\Big) \Big).
\end{align*}

With the cyclic rule of trace it follows that

\begin{align*}
\delta C=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{0:N}^{\text{m}}\Big(   \mathrm{tr}_{1:N}^{\text{in},1:L,\text{out}}\Big( U^{l+1^\dagger}_z\dots U^{L+1^\dagger}_z\U_{z+1}^\dagger \dots \U_N^\dagger \left(\sigma_1^\text{out}\otimes\dots \otimes \sigma_N^\text{out} \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L}\right)\\
& \U_N \dots \U_{z+1} U^{L+1}_z\dots U^{l+1}_z U^{l}_{m_{l}\, z} \dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \U_{z-1}\dots \U_1 \\
&\left( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{1:N}^{1:L+1}\bra{0\dots 0} \right)\U_1^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{l^\dagger}_{m_{l}\, z}  \Big)\Big)\\
=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{0:N}^{\text{m}}\Big(   \mathrm{tr}_{1:N}^{\text{in},1:L,\text{out}}\Big( \Big( \sigma_1^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out} \otimes U^{l+1^\dagger}_z\dots U^{L+1^\dagger}_z\U_{z+1}^\dagger \dots \U_N^\dagger \Big(\sigma_z^\text{out}\otimes\dots \otimes \sigma_N^\text{out} \\
&\otimes \id_{z-1:N}^\text{m}\otimes \id_{z:N}^{\text{in},1:L}\Big) \U_N \dots \U_{z+1} U^{L+1}_z\dots U^{l+1}_z  \otimes \id_{0:z-2}^\text{m}\otimes \id_{1:z-1}^{\text{in},1:L} \Big) \Big(U^{l}_{m_{l}\, z} \dots U^l_{j+1\,z} \\
&\Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \U_{z-1}\dots \U_1 \left( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes\rho_z^\text{in} \otimes \ket{0\dots 0}_{1:z}^{1:L+1}\bra{0\dots 0} \right)\U_1^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] \\
&U^{l^\dagger}_{j+1\, z}\dots U^{l^\dagger}_{m_{l}\, z}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)  \Big)\Big).
\end{align*}

By writing

\begin{align}
A=&U^{l+1^\dagger}_z\dots U^{L+1^\dagger}_z\U_{z+1}^\dagger \dots \U_N^\dagger \Big(\sigma_z^\text{out}\otimes\dots \otimes \sigma_N^\text{out}\otimes \id_{z-1:N}^\text{m}\otimes \id_{z:N}^{\text{in},1:L}\Big) \U_N \dots \U_{z+1} U^{L+1}_z\dots U^{l+1}_z,\\
B=&U^{l}_{m_{l}\, z} \dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \U_{z-1}\dots \U_1 \left( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes\rho_z^\text{in} \otimes \ket{0\dots 0}_{1:z}^{1:L+1}\bra{0\dots 0} \right)\nonumber\\
&\U_1^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{l^\dagger}_{m_{l}\, z}
\end{align}

in shorthand we obtain

\begin{align*}
\delta C	=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{0:N}^{\text{m}}\Big(   \mathrm{tr}_{1:N}^{\text{in},1:L,\text{out}}\Big( \Big( \sigma_1^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out} \otimes A  \otimes \id_{0:z-2}^\text{m}\otimes \id_{1:z-1}^{\text{in},1:L} \Big)\\
& \Big(B\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)  \Big)\Big)\\
=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{0:N}^{\text{m}}\Big(   \mathrm{tr}_{1:N}^{\text{in},1:L,\text{out}}\Big( \Big(  A  \otimes \id_{0:z-2}^\text{m}\otimes \id_{1:z-1}^{\text{in},1:L,\text{out}} \Big) \Big( \sigma_1^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out}   \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L} \otimes \id_{z:N}^\text{out} \Big) \\
& \Big(B\otimes \id_{z+1:N}^{\text{in},1:L+1}\Big) \Big(\id_{0:z}^\text{m}\otimes \id_{1:z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big) \Big)\Big)\\
=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{0:N}^{\text{m}}\Big(   \mathrm{tr}_{1:N}^{\text{in},1:L,\text{out}}\Big( \Big(\id_{0:z}^\text{m}\otimes \id_{1:z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)  \\
& \Big(  A  \otimes \id_{0:z-2}^\text{m}\otimes \id_{1:z-1}^{\text{in},1:L,\text{out}} \Big)  \Big( \sigma_1^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out}   \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L} \otimes \id_{z:N}^\text{out} \Big) \Big(B\otimes \id_{z+1:N}^{\text{in},1:L+1}\Big)  \Big)\Big)\\
=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{z-1:N}^{\text{m}}\Big(   \mathrm{tr}_{z:N}^{\text{in},1:L,\text{out}}\Big( \Big(\id_{z-1,z}^\text{m}\otimes \id_{z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)  A\\
& \Big(\mathrm{tr}_{0:z-2}^\text{m}\Big(\mathrm{tr}_{1:z-1}^{\text{in},1:L,\text{out}}\Big(\Big( \sigma_1^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out}   \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L} \otimes \id_{z:N}^\text{out} \Big) B\Big)\Big) \otimes \id_{z+1:N}^{\text{in},1:L+1} \Big) \Big)\Big).
\end{align*}

By substituting

\begin{equation}
E=\mathrm{tr}_{0:z-2}^\text{m}\Big(\mathrm{tr}_{1:z-1}^{\text{in},1:L,\text{out}}\Big(\Big( \sigma_1^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out}   \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L} \otimes \id_{z:N}^\text{out} \Big) B\Big)\Big)
\end{equation}

we get

\begin{align*}
\delta C	=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{z-1:N}^{\text{m}}\Big(   \mathrm{tr}_{z:N}^{\text{in},1:L,\text{out}}\Big( \Big(\id_{z-1,z}^\text{m}\otimes \id_{z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in} \otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)  A\\
& \Big(E \otimes \id_{z+1:N}^{\text{in},1:L+1} \Big) \Big)\Big)\\
=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{z-1:z}^{\text{m}}\Big(   \mathrm{tr}_{z}^{\text{in},1:L,\text{out}}\Big( \mathrm{tr}_{z+1:N}^{\text{in},1:L,\text{out,m}}\Big(\Big(\id_{z-1,z}^\text{m}\otimes \id_{z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in} \\
&\otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)  A\Big)E \Big)\Big).
\end{align*}

By further substituting

\begin{equation}
D=\mathrm{tr}_{z+1:N}^{\text{in},1:L,\text{out,m}}\Big(\Big(\id_{z-1,z}^\text{m}\otimes \id_{z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in}\otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)  A\Big)
\end{equation}

we obtain

\begin{align*}
\delta C	=& i  \sum_{z=1}^{N} \sum_{l=1}^{L+1} \sum_{j=1}^{m_l} \mathrm{tr}_{z-1:z}^{\text{m}}\Big(   \mathrm{tr}_{z}^{\text{in},1:L,\text{out}}\Big( DE \Big)\Big).
\end{align*}

By resubstituting \(B\) in \(E\) we get

\begin{align}
E=&\mathrm{tr}_{0:z-2}^\text{m}\Big(\mathrm{tr}_{1:z-1}^{\text{in},1:L,\text{out}}\Big(\Big( \sigma_1^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out}   \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L} \otimes \id_{z:N}^\text{out} \Big) U^{l}_{m_{l}\, z} \dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \nonumber\\
& \U_{z-1}\dots \U_1 \left( \rho_0^\text{m}\otimes \rho_1^\text{in}\otimes \dots \otimes\rho_z^\text{in} \otimes \ket{0\dots 0}_{1:z}^{1:L+1}\bra{0\dots 0} \right)\U_1^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{l^\dagger}_{m_{l}\, z}  \Big)\Big) \nonumber\\
=&\mathrm{tr}_{1:z-2}^\text{m}\Big(\mathrm{tr}_{2:z-1}^{\text{in},1:L,\text{out}}\Big(\Big( \sigma_2^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out}   \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L} \otimes \id_{z:N}^\text{out} \Big) U^{l}_{m_{l}\, z} \dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \U_{z-1}\dots \U_2 \nonumber\\
&  \Big( \mathrm{tr}_{0}^\text{m} \Big( \mathrm{tr}_{1:z-1}^{\text{in},1:L,\text{out}} \Big(\left(\sigma_1^\text{out}\otimes \id_1^{\text{in},1:L}\otimes \id_{0,1}^\text{m}\right)\U_1 \left( \rho_0^\text{m}\otimes \rho_1^\text{in} \otimes \ket{0\dots 0}_{1}^{1:L+1}\bra{0\dots 0} \right)\U_1^\dagger\Big)\Big) \otimes \rho_2^\text{in}\otimes \dots \otimes\rho_z^\text{in}\nonumber \\
&\otimes \ket{0\dots 0}_{2:z}^{1:L+1}\bra{0\dots 0} \Big)\U_2^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{l^\dagger}_{m_{l}\, z}  \Big)\Big) \nonumber\\
=&\mathrm{tr}_{1:z-2}^\text{m}\Big(\mathrm{tr}_{2:z-1}^{\text{in},1:L,\text{out}}\Big(\Big( \sigma_2^\text{out}\otimes\dots \otimes \sigma_{z-1}^\text{out}   \otimes \id_{0:N}^\text{m}\otimes \id_{1:N}^{\text{in},1:L} \otimes \id_{z:N}^\text{out} \Big) U^{l}_{m_{l}\, z} \dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \U_{z-1}\dots \U_2 \nonumber\\
&  \Big( \tau_1^\text{m} \otimes \rho_2^\text{in}\otimes \dots \otimes\rho_z^\text{in}\nonumber\otimes \ket{0\dots 0}_{2:z}^{1:L+1}\bra{0\dots 0} \Big)\U_2^\dagger \dots \U_{z-1}^\dagger U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{l^\dagger}_{m_{l}\, z}  \Big)\Big) \nonumber\\
=& U^{l}_{m_{l}\, z} \dots U^l_{j+1\,z} \Big[ K^l_{j\,z}, U^l_{j\,z} \dots U^1_{1\,z} \Big( \tau_{z-1}^\text{m} \otimes\rho_z^\text{in}\otimes \ket{0\dots 0}_{z}^{1:L+1}\bra{0\dots 0} \Big) U^{1^\dagger}_{1\, z}\dots U^{l^\dagger}_{j\, z}    \Big] U^{l^\dagger}_{j+1\, z}\dots U^{l^\dagger}_{m_{l}\, z} .
\end{align}

Also by resubstituting \(A\) in \(D\) it is

\begin{align}
D=&\mathrm{tr}_{z+1:N}^{\text{in},1:L,\text{out,m}}\Big(\Big(\id_{z-1,z}^\text{m}\otimes \id_{z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_N^\text{in}\otimes \ket{0\dots 0}_{z+1:N}^{1:L+1}\bra{0\dots 0} \Big)   U^{l+1^\dagger}_z\dots U^{L+1^\dagger}_z\nonumber\\
&\U_{z+1}^\dagger \dots \U_N^\dagger \Big(\sigma_z^\text{out}\otimes\dots \otimes \sigma_N^\text{out}\otimes \id_{z-1:N}^\text{m}\otimes \id_{z:N}^{\text{in},1:L}\Big) \U_N \dots \U_{z+1} U^{L+1}_z\dots U^{l+1}_z \Big) \nonumber \\
=&\mathrm{tr}_{z+1:N-1}^{\text{in},1:L,\text{out,m}}\Big(\Big(\id_{z-1,z}^\text{m}\otimes \id_{z}^{\text{in},1:L,\text{out}}\otimes \rho_{z+1}^\text{in}\otimes \dots \otimes\rho_{N-1}^\text{in}\otimes \ket{0\dots 0}_{z+1:N-1}^{1:L+1}\bra{0\dots 0} \Big)   U^{l+1^\dagger}_z\dots U^{L+1^\dagger}_z\nonumber\U_{z+1}^\dagger \dots \U_{N-1}^\dagger \\
&\Big( \sigma_z^\text{out}\otimes\dots \otimes \sigma_{N-1}^\text{out}\otimes \id_{z-1:N-2}^\text{m}\otimes \id_{z:N-1}^{\text{in},1:L} \otimes \mathrm{tr}_{N}^{\text{in},1:L,\text{out,m}} \Big(\Big(\id_{N-1}^\text{m}\otimes \rho_N^\text{in} \otimes \ket{0\dots 0}_{N}^{1:L+1}\bra{0\dots 0} \Big)\U_N^\dagger \nonumber \\
&\Big(\sigma_N^\text{out}\otimes \id_{N-1,N}^\text{m}\otimes \id_{N}^{\text{in},1:L}\Big) \U_N\Big)\Big) \U_{N-1}\dots \U_{z+1} U^{L+1}_z\dots U^{l+1}_z \Big) \nonumber .
\end{align}

Let \(\chi_N^\text{m}=\id_N^\text{m}\), \(\chi_x^{L+1}=\chi_{x}^\text{m}\otimes \sigma_x^\text{out}\),

\[\chi_x^l=\mathcal{F}^{l+1}\left(\chi_x^{l+1}\right)\]

\begin{align}
\chi_x^\text{m}=&\mathrm{tr}_{x+1}^{\text{in},1:L+1} \Big(\Big(\id_{x}^\text{m}\otimes \rho_{x+1}^\text{in} \otimes \ket{0\dots 0}_{x+1}^{1:L+1}\bra{0\dots 0} \Big)\U_{x+1}^\dagger\Big( \chi_{x+1}^\text{m}\otimes \sigma_{x+1}^\text{out}\otimes \id_{x}^\text{m}\otimes \id_{x+1}^{\text{in},1:L}\Big) \U_{x+1}\Big)\\
=&\mathrm{tr}_{x+1}^{\text{in}} \Big(\Big(\id_x^\text{m}\otimes \rho_x^\text{in}\Big)\chi_{y+1}^0\Big)
\end{align}

for \(x=1,...,N\) which, analogously to \(\tau_x^\text{m}\), do not have trace \(1\).
Hence, we get

\begin{align}