\begin{align}
\label{eq:def:qba:msg}
\sigma_a &\propto \rho_a \star \bigstar_{i\in\nb{a}} m_{i\to a},\\
\label{eq:def:qbi:msg}
\sigma_i &\propto \bigstar_{a\in\nb{a}} m_{a\to i},
\end{align}
where $\{m_{i\to a},  m_{a\to i}\in \DensOp(\hilbert_i)\}_{(i,a)\in\set{E}}$ are (straightly) PD messages such that
\begin{align}
\label{eq:msg:QFG:update:2:ia}
m_{i\to a} &\propto \bigstar_{c\in\nb{i}\xk{a}} m_{c\to i}, \\
\label{eq:msg:QFG:update:2:ai}
m_{a\to i} &\propto \tr_{\nb{a}\xk i}\left( \rho_a \star \Tensor_{k\in\nb{a}} m_{k\to a} \right) \star m_{i\to a}^{-1}.
\end{align}
Conversely,~\eqref{eq:msg:QFG:update:2:ia} and~\eqref{eq:msg:QFG:update:2:ai} hold if $\left(\{\sigma_a\}_{a\in\set{F}},\{\sigma_i\}_{i\in\set{V}}\right)$ defined according to~\eqref{eq:def:qba:msg} and~\eqref{eq:def:qbi:msg} form an internal stationary point of $\bethe$.
\end{theorem}
\begin{proof}
This theorem is a quantum analogy of Theorem~\ref{thm:bethe:BPA}.
Part of the ideas of this proof originated from~\cite{yedidia2005constructing}.
\par
Suppose $\left(\{\sigma_a\}_{a\in\set{F}},\{\sigma_i\}_{i\in\set{V}}\right)\in\set{L}(\set{G})$ is an \emph{interior} stationary point of the constrained Bethe minimization problem.
In this case, the Lagrangian of this problem can be expressed as
\[
L \defeq \bethe + \sum_{a\in\set{F}} \gamma_a \cdot \left( \tr\left( \sigma_a \right)-1 \right) + \sum_{i\in\set{V}} \gamma_i \cdot \left( \tr\left( \sigma_i \right) - 1 \right) + \sum_{(i,a)\in\set{E}} \tr\left( \lambda_{i,a} \cdot \left(\sigma_i - \tr_{\partial{a}\xk i}\left( \sigma_a \right) \right) \right),
\]
where $\{\gamma_a\in\Reals\}_{a\in\set{F}}$, $\{\gamma_i\in\Reals\}_{i\in\set{V}}$, and $\{\lambda_{i,a}\in \LinearOp(\hilbert_i)\}_{(i,a)\in\set{E}}$ are Lagrangian dual variables.
At the stationary point, it must hold that
\begin{align}
\label{eq:kkt:1}
\frac{\partial{L}}{\partial{\gamma_a}} &=0 &&\forall a\in\set{F},\\
\label{eq:kkt:2}
\frac{\partial{L}}{\partial{\gamma_i}} &=0 &&\forall i\in\set{V},\\
\label{eq:kkt:3}
\left.\frac{\D}{\D t}\right|_{t=0} L(\lambda_{i,a} + t\cdot H_{i,a}) &= 0 &&\forall H_{i,a}\in\HermitianOp(\hilbert_i),\, \forall (i,a)\in\set{E},\\
\label{eq:kkt:4}
\left.\frac{\D}{\D t}\right|_{t=0} L(\sigma_a + t\cdot H_a) &= 0 &&\forall H_a\in\HermitianOp(\hilbert_\nb{a}),\, \forall a\in\set{F},\\
\label{eq:kkt:5}
\left.\frac{\D}{\D t}\right|_{t=0} L(\sigma_i + t\cdot H_i) &= 0, &&\forall H_i\in\HermitianOp(\hilbert_i),\, \forall i\in\set{V}.
\end{align}
Eqs.~\eqref{eq:kkt:1},~\eqref{eq:kkt:2},~\eqref{eq:kkt:3} are equivalent to the constraints of the problem.
By Spectral theorem and the first-order perturbation theory~\cite{Kato1966},~\eqref{eq:kkt:4} and~\eqref{eq:kkt:5} can be expanded as
\begin{align*}
-\tr\left( H_a \cdot \log{\rho_a} \right) + \tr\left( H_a \cdot \left( I+\log{\sigma_a} \right) \right) + \tr(H_a) \cdot \gamma_a - \sum_{i\in\nb{a}} \tr\left( \lambda_{i,a} \cdot \tr_{\partial{a}\xk i} H_a \right) & = 0,\\
(1-d_i) \cdot \tr\left( H_i \cdot \left( I+\log{\sigma_i} \right) \right) + \tr(H_i) \cdot \gamma_i + \sum_{a\in\nb{i}} \tr\left( H_i \cdot \lambda_{i,a} \right) & = 0.
\end{align*}
Solving the above equations for $\{\sigma_a\}_{a\in\set{F}}$ and $\{\sigma_i\}_{i\in\set{V}}$, respectively, we have
\begin{align}
\label{eq:kkt:4b}
\sigma_a &= \exp\left( \log{\rho_a} + \sum_{i\in\nb{a}} \lambda_{i,a} - (1+\gamma_a) \cdot I \right) && \forall a\in\set{F},\\
\label{eq:kkt:5b}
\sigma_i &= \exp\left( \frac{1}{d_i-1} \cdot \left( \sum_{a\in\nb{i}} \lambda_{i,a} + \left(1+\gamma_i\right) \cdot I \right) \right) && \forall i\in\set{V}.
\end{align}
Eqs.~\eqref{eq:def:qba:msg} and~\eqref{eq:def:qbi:msg} can be shown by making the identification that
\begin{equation}\label{eq:QBP:lambda}
\lambda_{i,a} = \log{m_{i\to a}} = \sum_{c\in\nb{i}\xk a}\log{m_{c\to i}} \qquad \forall (i,a)\in\set{E}.
\end{equation}
In this case,~\eqref{eq:msg:QFG:update:2:ai} can be derived from the constraints that $\tr_{\nb{a}\xk{i}} \sigma_a = \sigma_i$ for each $(i,a){\in\set{E}}$.
\par
Conversely, suppose there exist some PD messages $\{m_{i\to a},  m_{a\to i}\in \DensOp(\hilbert_i)\}_{(i,a)\in\set{E}}$ satisfying~\eqref{eq:msg:QFG:update:2:ia} and~\eqref{eq:msg:QFG:update:2:ai}.
Let $\left(\{\sigma_a\}_{a\in\set{F}},\{\sigma_i\}_{i\in\set{V}}\right)\in\set{L}(\set{G})$ be defined according to~\eqref{eq:def:qba:msg} and~\eqref{eq:def:qbi:msg}.
By choosing $\{\gamma_a\}_{a\in\set{F}}$, $\{\gamma_i\}_{i\in\set{V}}$, and $\{\lambda_{i,a}\}_{(i,a)\in\set{E}}$ by such that~\eqref{eq:kkt:4b}, \eqref{eq:kkt:5b}, and~\eqref{eq:QBP:lambda} hold simultaneously, we have satisfied~\eqref{eq:kkt:4} and~\eqref{eq:kkt:5}.
Eqs.~\eqref{eq:kkt:1} and~\eqref{eq:kkt:2} also hold since $\{\sigma_a\}_{a\in\set{F}}$ and $\{\sigma_i\}_{i\in\set{V}}$ are density operators.
Finally,~\eqref{eq:kkt:3} holds since $\tr_{\nb{a}\xk{i}} \sigma_a = \sigma_i$ for each $(i,a){\in\set{E}}$, which in turn can be derived from~\eqref{eq:msg:QFG:update:2:ai}.
Thus, we have verified $\left(\{\sigma_a\}_{a\in\set{F}},\{\sigma_i\}_{i\in\set{V}}\right)$ to be a stationary point.
\end{proof}
\par
Under suitable conditions, by Proposition~\ref{prop:approx:distri},~\eqref{eq:msg:QFG:update:2:ai} can be approximated as
\begin{equation}\label{eq:msg:QFG:update:3:ai}
m_{a\to i} \propto \tr_{\nb{a}\xk i}\left( \rho_a \star \Tensor_{k\in\nb{a}} m_{k\to a} \right) \star m_{i\to a}^{-1}
\approx \tr_{\nb{a}\xk i}\left( \rho_a \star \Tensor_{k\in\nb{a}\xk{i}} m_{k\to a} \right),
\end{equation}
which is exactly the second half of the BP fixed point condition (recall Definition~\ref{def:QFG:BP:fixed:points}).
Therefore, Theorem~\ref{thm:bethe:QBPA} provides an interpretation of the quantum BP algorithm as an iterative method for finding a stationary point of the constrained Bethe minimization problem.
Unfortunately, there is no guarantee of convergence of such an iterative method.
Moreover, even for acyclic QFGs, we have made multiple approximations to link the BP fixed points to the partition sum (including Theorem~\ref{thm:quantum:Bethe:Gibbs} and Eq.~\eqref{eq:msg:QFG:update:3:ai}).
Despite such concerns, as illustrated in the next section, the quantum BP algorithm still shows a rather promising performance in some numerical applications.
\par
The following corollary generalizes Corollary~\ref{cor:Zb:decompose} and is a direct result of Theorem~\ref{thm:bethe:QBPA}.
\begin{corollary}
Consider a QFG $\set{G}$ representing the factorization $\rho=\bigstar_{a\in\set{F}} \rho_a$.
For a collection of positive messages $\{m_{i\to a}, m_{a\to i}\}_{i,a}$ satisfying~\eqref{eq:msg:QFG:update:2:ia} and~\eqref{eq:msg:QFG:update:2:ai}, we have
\begin{align}
&\defeq \exp\left(-\bethe(\{\sigma_a\}_{a\in\set{F}}, \{\sigma_i\}_{i\in\set{V}})\right)\\
&= Z_{\mathsf{induced}}(\{m_{i\to a},m_{a\to i}\}_{(i,a)\in\set{E}}),
\end{align}
where the density operators $\{\sigma_a\}_{a\in\set{F}}$ and $\{\sigma_i\}_{i\in\set{V}}$ are defined according to~\eqref{eq:def:qba:msg} and~\eqref{eq:def:qbi:msg}.
\end{corollary}
\section{Numerical Example}\label{sec:QFG:numerical}
\begin{figure}\centering
\begin{tikzpicture}[factor/.style={rectangle, minimum size=1cm, draw}]
\node[draw=none,inner sep = 0pt, outer sep = 0pt] (O) {};
\foreach \t in {1,...,6}{
\pgfmathsetmacro\posX{2*cos(60*\t)};
\pgfmathsetmacro\posY{2*sin(60*\t)};
\node[factor] (P\t) at ([xshift=-\posX cm, yshift=\posY cm]O) {$\rho_\t$};
\foreach \t in {1,...,6}{
\pgfmathtruncatemacro\a{mod(\t,6)+1};
\pgfmathtruncatemacro\b{mod(\t+1,6)+1};
\path (P\t) edge [line width = 2pt] (P\a);
\path (P\t) edge [line width = 2pt] (P\b);
\end{tikzpicture}
\caption{QFG in Section~\ref{sec:QFG:numerical}.}
\label{fig:QFG:numerical}
\end{figure}
\begin{figure}\centering
\begin{tikzpicture}[nodes = {draw=none, minimum size = 0pt, inner sep = 0pt, outer sep = 0pt}]
% Scatter Plot
\node (plot) {\includegraphics[width=6cm]{1155001913-cao_michael_xuan-202104-phd-Figure_QFG_Numerical_Example_1.png}};
\node (o) at (plot.south west) {};
\node (xend) at (plot.south east) {};
\node (yend) at (plot.north west) {};
\path (o) edge[draw=none]
node[anchor=north, yshift=-15pt] {$Z$}
node[pos=0, anchor=north, yshift=-2pt, font=\scriptsize] {$10^{1}$}
node[pos=.7505, anchor=north, yshift=-2pt, font=\scriptsize] {$10^{2}$}
\path (o) edge[draw=none]
node[anchor=south, yshift=20pt, sloped] {$Z_{\mathsf{induced}}$}
node[pos=0, anchor=east, xshift=-2pt, font=\scriptsize] {$10^{1}$}
node[pos=.7505, anchor=east, xshift=-2pt, font=\scriptsize] {$10^{2}$}
\draw[red,dashed] (o) -- (plot.north east);
% Histogram
\node[right = 35pt of plot.east, anchor=west] (plotH) {\includegraphics[width=6cm]{1155001913-cao_michael_xuan-202104-phd-Figure_QFG_Numerical_Example_1h.png}};
\node (o) at (plotH.south west) {};
\node (xend) at (plotH.south east) {};
\node (yend) at (plotH.north west) {};
\path (o) edge[draw=none]
node[anchor=north, yshift=-15pt] {Relative Error}
node[pos=0, anchor=north, yshift=-2pt, font=\scriptsize] {$0$}
node[pos=.1111, anchor=north, yshift=-2pt, font=\scriptsize] {$.05$}
node[pos=.2222, anchor=north, yshift=-2pt, font=\scriptsize] {$.10$}
node[pos=.3333, anchor=north, yshift=-2pt, font=\scriptsize] {$.15$}
node[pos=.4444, anchor=north, yshift=-2pt, font=\scriptsize] {$.20$}
node[pos=.5555, anchor=north, yshift=-2pt, font=\scriptsize] {$.25$}
node[pos=.6666, anchor=north, yshift=-2pt, font=\scriptsize] {$.30$}
node[pos=.7777, anchor=north, yshift=-2pt, font=\scriptsize] {$.35$}
node[pos=.8888, anchor=north, yshift=-2pt, font=\scriptsize] {$.40$}
\path (o) edge[draw=none]
node[anchor=south, yshift=20pt, sloped] {Number of Appearances}
node[pos=0, anchor=east, xshift=-2pt, font=\scriptsize] {$0$}
node[pos=.1, anchor=east, xshift=-2pt, font=\scriptsize] {$100$}
node[pos=.2, anchor=east, xshift=-2pt, font=\scriptsize] {$200$}
node[pos=.3, anchor=east, xshift=-2pt, font=\scriptsize] {$300$}
node[pos=.4, anchor=east, xshift=-2pt, font=\scriptsize] {$400$}
node[pos=.5, anchor=east, xshift=-2pt, font=\scriptsize] {$500$}
node[pos=.6, anchor=east, xshift=-2pt, font=\scriptsize] {$600$}
node[pos=.7, anchor=east, xshift=-2pt, font=\scriptsize] {$700$}
node[pos=.8, anchor=east, xshift=-2pt, font=\scriptsize] {$800$}
node[pos=.9, anchor=east, xshift=-2pt, font=\scriptsize] {$900$}
node[pos=1, anchor=east, xshift=-2pt, font=\scriptsize] {$1000$}
\end{tikzpicture}
\caption{The partition sum $Z$ and the induced partition sum $Z_{\mathsf{induced}}$ produced by Algorithm~\ref{alg:BP:QFG}. The plot consists of 10,000 instances based on the QFG in Figure~\ref{fig:QFG:numerical}, where each local factors $\{\rho_i\}_{i=1}^6$ are generated independently as $\rho_i\gets U\cdot \mbox{diag}(\boldsymbol\lambda_1^{16}) \cdot U$, and where (for each $\rho_i$) $U$ is some $16$-by-$16$ unitary matrix (independently) randomly generated according to Haar measure, and where each $\{\lambda_k\}_{k=1}^{16}$ are independently uniformly distributed on $[0,1]$.}
\label{fig:QFG:numerical:1}
\end{figure}
\begin{figure}\centering
\begin{tikzpicture}[nodes = {draw=none, minimum size = 0pt, inner sep = 0pt, outer sep = 0pt}]
% Scatter Plot
\node (plot) {\includegraphics[width=6cm]{1155001913-cao_michael_xuan-202104-phd-Figure_QFG_Numerical_Example_2.png}};
\node (o) at (plot.south west) {};
\node (xend) at (plot.south east) {};
\node (yend) at (plot.north west) {};
\path (o) edge[draw=none]
node[anchor=north, yshift=-15pt] {$Z$}
node[pos=0, anchor=north, yshift=-2pt, font=\scriptsize] {$10^{3}$}
node[pos=.6225, anchor=north, yshift=-2pt, font=\scriptsize] {$10^{4}$}
\path (o) edge[draw=none]
node[anchor=south, yshift=20pt, sloped] {$Z_{\mathsf{induced}}$}
node[pos=0, anchor=east, xshift=-2pt, font=\scriptsize] {$10^{3}$}
node[pos=.6225, anchor=east, xshift=-2pt, font=\scriptsize] {$10^{4}$}
\draw[red,dashed] (o) -- (plot.north east);
% Histogram
\node[right = 35pt of plot.east, anchor=west] (plotH) {\includegraphics[width=6cm]{1155001913-cao_michael_xuan-202104-phd-Figure_QFG_Numerical_Example_2h.png}};
\node (o) at (plotH.south west) {};
\node (xend) at (plotH.south east) {};
\node (yend) at (plotH.north west) {};
\path (o) edge[draw=none]
node[anchor=north, yshift=-15pt] {Relative Error}
node[pos=0, anchor=north, yshift=-2pt, font=\scriptsize] {$0$}
node[pos=.0833, anchor=north, yshift=-2pt, font=\scriptsize] {$.05$}
node[pos=.1667, anchor=north, yshift=-2pt, font=\scriptsize] {$.10$}
node[pos=.2500, anchor=north, yshift=-2pt, font=\scriptsize] {$.15$}
node[pos=.3333, anchor=north, yshift=-2pt, font=\scriptsize] {$.20$}
node[pos=.4167, anchor=north, yshift=-2pt, font=\scriptsize] {$.25$}
node[pos=.5000, anchor=north, yshift=-2pt, font=\scriptsize] {$.30$}
node[pos=.5833, anchor=north, yshift=-2pt, font=\scriptsize] {$.35$}
node[pos=.6667, anchor=north, yshift=-2pt, font=\scriptsize] {$.40$}
node[pos=.7500, anchor=north, yshift=-2pt, font=\scriptsize] {$.45$}
node[pos=.8333, anchor=north, yshift=-2pt, font=\scriptsize] {$.50$}
node[pos=.9167, anchor=north, yshift=-2pt, font=\scriptsize] {$.55$}
\path (o) edge[draw=none]
node[anchor=south, yshift=20pt, sloped] {Number of Appearances}
node[pos=0, anchor=east, xshift=-2pt, font=\scriptsize] {$0$}
node[pos=.125, anchor=east, xshift=-2pt, font=\scriptsize] {$100$}
node[pos=.250, anchor=east, xshift=-2pt, font=\scriptsize] {$200$}
node[pos=.375, anchor=east, xshift=-2pt, font=\scriptsize] {$300$}
node[pos=.500, anchor=east, xshift=-2pt, font=\scriptsize] {$400$}
node[pos=.625, anchor=east, xshift=-2pt, font=\scriptsize] {$500$}
node[pos=.750, anchor=east, xshift=-2pt, font=\scriptsize] {$600$}
node[pos=.875, anchor=east, xshift=-2pt, font=\scriptsize] {$700$}
node[pos=1, anchor=east, xshift=-2pt, font=\scriptsize] {$800$}
\end{tikzpicture}
\caption{The partition sum $Z$ and the induced partition sum $Z_{\mathsf{induced}}$ produced by Algorithm~\ref{alg:BP:QFG}. The plot consists of 10,000 instances based on the QFG in Figure~\ref{fig:QFG:numerical}, where each local factors $\{\rho_i\}_{i=1}^6$ are generated independently as $\rho_i\gets U\cdot \mbox{diag}(\boldsymbol\lambda_1^{16}) \cdot U$, and where (for each $\rho_i$) $U$ is some $16$-by-$16$ unitary matrix (independently) randomly generated according to Haar measure, and where each $\{\lambda_k\}_{k=1}^{16}$ are independently distributed according to $\abs{\mathcal{N}(1,1)}$. ($\mathcal{N}(\mu,\sigma)$ denotes the Gaussian distribution with mean $\mu$ and standard variation $\sigma$.)}
\label{fig:QFG:numerical:2}
\end{figure}
In this section, we consider the QFG as in Figure.~\ref{fig:QFG:numerical}, where each local operator $\{\rho_a\}_{a=1,\ldots,6}$ is generated in the same fashion as $\rho_a$ and $\rho_b$ for Figure~\ref{fig:approx:distributive:star}.
We are interested in comparing the induced partition sum $Z_{\mathsf{induced}}$ at BP fixed points with the actual partition sum $Z(\set{G})$.
In particular, in Figure~\ref{fig:QFG:numerical:1}~and~\ref{fig:QFG:numerical:2}, we have plotted the statistics of the relative error $\eta\defeq\abs{(Z_{\mathsf{induced}}-Z)/Z)}$ \wrt different distributions of the eigenvalues.
The plots are based on $10,000$ simulations.
% Chapter 4 Case Study**********************************************************
\chapter[IRs of Quantum Channels with Memory]{Bounding and Estimating the Classical Information Rate of Quantum Channels with Memory}\label{chapter:QCwM}
In this chapter, we consider the transmission rate of classical information over a finite-dimensional quantum channel with memory.
A \emph{quantum channel with memory}~\cite{bowen2004quantum, kretschmann2005quantum, caruso2014quantum} is a quantum channel (from input system $\system{A}$ to output system $\system{B}$) equipped with a memory system $\system{S}$; namely it is a CPTP map from the set of density operators on $\hilbert_\system{A}\tensor\hilbert_\system{S}$ to the set of  density operators on $\hilbert_\system{B}\tensor\hilbert_{\system{S}'}$.
The system $\system{S}$ can be understood either as a state of the channel (as illustrated in Figure~\ref{fig:interpret:1}), or as a part of the environment that does not decay between consecutive channel uses (as illustrated in Figure~\ref{fig:interpret:2}).
Interesting examples of quantum channels with  memory include spin chains~\cite{bose2003quantum} and fiber optic links~\cite{ball2004exploiting}.
\begin{figure}[t]\centering
\begin{subfigure}[b]{0.41\textwidth}\centering
\begin{tikzpicture}[nodes/.style={draw=none. inner sep= 0pt, outer sep= 0pt},
smallcircle/.style={draw, minimum size = 12pt, circle}]
\node[draw, minimum width = 30pt, minimum height = 22.5pt] (C) {$\operator{N}$};
\node[smallcircle, below = 18pt of C] (M) {};
\node[left = 3pt of M.center, anchor=center, smallcircle] {};
\node[right = 3pt of M.center, anchor=center, smallcircle] {};
\draw (C.west|-M.south) -- (C.east|-M.south);
\draw (C.west|-M.south) arc (270:90:18pt);
\draw (C.east|-M.south) arc (-90:90:18pt);
\draw ([yshift=5pt]C.west) -- ([yshift=5pt, xshift=-16pt]C.west) node[pos=1, left](A) {$\system{A}$};
\draw ([yshift=5pt]C.east) -- ([yshift=5pt, xshift=16pt]C.east) node[pos=1, right](B) {$\system{B}$};
\node[yshift=16pt] at (A|-M.south) {$\system{S}$};
\node[yshift=16pt] at (B|-M.south) {$\system{S}'$};
\node[below=0pt of M, font=\small] {quantum memory};
\end{tikzpicture}
\caption{Memory as the state of the channel.}
\label{fig:interpret:1}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}\centering
\begin{tikzpicture}[nodes/.style={draw=none, inner sep=0pt, outer sep=0pt}]
\node[draw, minimum width = 30pt, minimum height = 75pt] (U) {$\operator{U}$};
\draw ([yshift=27pt]U.west) -- ([yshift=27pt,xshift=-22.5pt]U.west) node[left, pos=1] (A) {$\system{A}$};
\draw (U.west) -- ([xshift=-22.5pt]U.west) node[left, pos=1] (E) {$\bra{0}$};
\draw[latex-] ([yshift=-27pt]U.west) -- ([yshift=-27pt,xshift=-22.5pt]U.west) -- ([yshift=-48pt,xshift=-22.5pt]U.west) -- ([yshift=-48pt,xshift=22.5pt]U.east) -- ([yshift=-27pt,xshift=22.5pt]U.east) -- ([yshift=-27pt]U.east);
\draw (U.east) -- ([xshift=13.5pt]U.east) node [pos=1, right, draw, minimum size=18pt] (M) {};
\node[minimum size = 2.25pt, fill = black, inner sep = 0pt, outer sep = 0pt, below = 4.5pt of M.center, circle, anchor = center] (m) {};
\draw[decoration={markings, mark = at position 1 with {\arrow[>=latex]{>}}}, postaction={decorate}, draw = none] (m.center)--([yshift=10.5pt,xshift=4.5pt]m.center);
\draw(m.center)--([yshift=9.45pt, xshift=4.05pt]m.center);
\draw ([xshift=6pt]m.center) arc (0:180:6pt);
\draw ([yshift=27pt]U.east) -- ([yshift=27pt,xshift=42pt]U.east) node[right, pos=1] (B) {$\system{B}$};
\draw[line width = 2pt] (M.east) -- ([xshift=20.25pt]M.east) node[pos=1, minimum size = 7.5pt, inner sep = 0pt, outer sep = 0pt] (Eend) {};
\draw[line width = 2pt] (Eend.south west) -- (Eend.north east);
\draw[line width = 2pt] (Eend.north west) -- (Eend.south east);

\node[yshift=-27pt] at (A|-U.west) (S) {$\system{S}$};
\node[yshift=-27pt] at (B|-U.east) {$\system{S}^\prime$};

\path[draw=none] (E) edge[draw=none] node[xshift=-18pt](env) {$\Bigg\{$} (S);
\node[left=0pt of env, anchor = center, rotate=90, font = \small] {environment};
\node[right=123pt of env] {$\Bigg\}$};
\end{tikzpicture}
\caption{Memory as undecayed partial environment.}
\label{fig:interpret:2}
\end{subfigure}
\caption{Interpretations of quantum channels with memory.}
\end{figure}
\par
Classical communication over such channels is accomplished by encoding classical data into some density operators before the transmission and applying measurements to the outputs of the channel.
In the most generic case, an ensemble and a measurement on the \emph{joint} input and output systems of multiple channels uses can be used for encoding and decoding, respectively.
The scenario involving a $k$-channel joint ensemble and a $k$-channel joint measurement is depicted in Figure~\ref{fig:generic_memory}, where
\begin{itemize}
\item the encoding process $\mathcal{E}$ is described by some ensemble $\{\prob_\rv{X}(x), \rho_{\system{A}_1^k}^{(x)}\}_{x\in\set{X}}$ on the joint input system $(\system{A}_1, \ldots, \system{A}_k)$, with $\set{X}$ being the input alphabet, $\prob_\rv{X}(x)$ being the input distribution, and $\rho_{\system{A}_1^k}^{(x)}$ being the density operator on the input systems $\system{A}_1^k$ corresponding to the classical input $x$;
\item the decoding process $\mathcal{D}$ is described by some positive-operator valued measurement (POVM) $\{\Lambda_{\system{B}_1^k}^{(y)}\}_{y\in\set{Y}}$ on the joint output system $(\system{B}_1,\ldots,\system{B}_k)$, with $\set{Y}$ being the output alphabet;
\item the classical input and output are represented by some random variables $\rv{X}$ and $\rv{Y}$, respectively.
\end{itemize}
For comparison, Figure~\ref{fig:generic_memoryless} shows the corresponding memoryless setup.
The above arrangement results in a (classical) channel from $\rv{X}$ to $\rv{Y}$, whose rate of transmission is given by
\begin{equation}\label{eq:capacity:jointk}
\infoRate(\mathcal{E},\operator{N}^{\boxtimes k},\mathcal{D})= \limsup_{n\to\infty}\frac{1}{n} \mutualInfo(\rv{X}_1^{n};\rv{Y}_1^{n}),
\end{equation}
where we use the above transmission scheme $n$ times consecutively (as depicted in Figure~\ref{fig:generic_memory_multiple}) and where
\[
\begin{aligned}
\operator{N}^{\boxtimes k} \defeq
&\left( \operator{N}_{\system{A}_k\system{S}_{k-1} \to \system{B}_k\system{S}_k} \otimes \id_{\system{B}_1^{k-1} \to \system{B}_1^{k-1}} \right) \circ
\left(\id_{\system{A}_k \to \system{A}_k} \otimes \operator{N}_{\system{A}_{k-1}\system{S}_{k-2} \to \system{B}_{k-1}\system{S}_{k-1}} \otimes \id_{\system{B}_1^{k-2} \to \system{B}_1^{k-2}} \right)\\
&\circ \cdots \circ
\Big( \id_{\system{A}_2^k \to \system{A}_2^k} \otimes \operator{N}_{\system{A}_1\system{S}_0 \to \system{B}_1\system{S}_1} \Big).
\end{aligned}
\]
Here, $\mutualInfo$ stands for the mutual information.
As a fundamental result, this quantity can be simplified to $\mutualInfo(\rv{X};\rv{Y})$ for the memoryless case~\cite{shannon1948mathematical, cover2012elements}.
Optimizing $\infoRate(\mathcal{E},\operator{N}^{\boxtimes k},\mathcal{D})$ over $\mathcal{E}$ and $\mathcal{D}$ (with $k\to\infty$) yields the classical capacity of the quantum channel $\operator{N}$, namely
\begin{equation}
\capacity(\operator{N}) = \limsup_{k}\frac{1}{k}\sup_{\mathcal{E},\mathcal{D}}
\infoRate(\mathcal{E},\operator{N}^{\boxtimes k},\mathcal{D}).
\end{equation}
\begin{figure}[t]\centering
\begin{subfigure}[b]{.49\textwidth}\centering
\begin{tikzpicture}[nodes/.style={draw=none, inner sep=0pt, outer sep= 0pt},
bignode/.style={minimum height = 100pt, minimum width = 20pt, draw},
smallnode/.style={minimum size = 15pt, draw}]
\node[bignode] (cq) {$\mathcal{E}$};
\node[bignode, right = 70pt of cq] (qc) {$\mathcal{D}$};
\path (cq.north) edge[draw=none] node[smallnode, anchor=north] (T1) {$\operator{N}$} (qc.north);
\path (cq.east|-T1) edge node[above] {$\system{A}_1$} (T1);
\path (T1) edge node[above] {$\system{B}_1$} (qc.west|-T1);
\path (cq.south) edge[draw=none] node[smallnode, anchor=south] (Tk) {$\operator{N}$} (qc.south);
\path (cq.east|-Tk) edge node[below] {$\system{A}_k$} (Tk);
\path (Tk) edge node[below] {$\system{B}_k$} (qc.west|-Tk);
\path (T1.center) edge[draw=none] node[smallnode, pos=.3333, anchor=center] (T2) {$\operator{N}$} node[pos=.6667, anchor=center] (Tdots) {$\vdots$} (Tk.center);
\path (cq.east|-T2) edge node[above] {$\system{A}_2$} (T2);
\path (T2) edge node[above] {$\system{B}_2$} (qc.west|-T2);
\path (cq.west) edge[line width=2pt] node[pos=1, above] {$\rv{X}$} ([xshift=-20pt]cq.west);
\path (qc.east) edge[line width=2pt] node[pos=1, above] {$\rv{Y}$} ([xshift=20pt]qc.east);
\draw[draw=none] ([yshift=10pt]T1.north) -- (T1.north) node[pos=0, above=0pt, color=white] {$\system{S}_0$};
\draw[draw=none] ([yshift=-10pt]Tk.south) -- (Tk.south) node[pos=0, below=0pt, color=white] {$\system{S}_k$};
\end{tikzpicture}
\settowidth{\somelength}{(a)}
\setlength{\somelength}{\dimexpr\textwidth-\somelength\relax}
\caption{\parbox[t]{\somelength}{Using the memoryless quantum channel $\operator{N}$ consecutively.}}
\label{fig:generic_memoryless}
\end{subfigure}
\begin{subfigure}[b]{.49\textwidth}\centering
\begin{tikzpicture}[nodes/.style={draw=none, inner sep=0pt, outer sep= 0pt},
bignode/.style={minimum height = 100pt, minimum width = 20pt, draw},
smallnode/.style={minimum size = 15pt, draw}]
\node[bignode] (cq) {$\mathcal{E}$};
\node[bignode, right = 70pt of cq] (qc) {$\mathcal{D}$};
\path (cq.north) edge[draw=none] node[smallnode, anchor=north] (T1) {$\operator{N}$} (qc.north);
\path (cq.east|-T1) edge node[above] {$\system{A}_1$} (T1);
\path (T1) edge node[above] {$\system{B}_1$} (qc.west|-T1);
\path (cq.south) edge[draw=none] node[smallnode, anchor=south] (Tk) {$\operator{N}$} (qc.south);
\path (cq.east|-Tk) edge node[below] {$\system{A}_k$} (Tk);
\path (Tk) edge node[below] {$\system{B}_k$} (qc.west|-Tk);
\path (T1.center) edge[draw=none] node[smallnode, pos=.3333, anchor=center] (T2) {$\operator{N}$} node[pos=.6667, anchor=center] (Tdots) {$\vdots$} (Tk.center);
\path (cq.east|-T2) edge node[above] {$\system{A}_2$} (T2);
\path (T2) edge node[above] {$\system{B}_2$} (qc.west|-T2);
\path (cq.west) edge[line width=2pt] node[pos=1, above] {$\rv{X}$} ([xshift=-20pt]cq.west);
\path (qc.east) edge[line width=2pt] node[pos=1, above] {$\rv{Y}$} ([xshift=20pt]qc.east);

\draw([yshift=10pt]T1.north) -- (T1.north) node[pos=0, above=0pt] {$\system{S}_0$};
\path (T1) edge node[right] {$\system{S}_1$} (T2);
\path (T2) edge node[pos=.75, right] {$\system{S}_2$} (Tdots);
\path (Tdots) edge node[pos=.25, right] {$\system{S}_{k-1}$} (Tk);
\draw ([yshift=-10pt]Tk.south) -- (Tk.south) node[pos=0, below=0pt] {$\system{S}_k$};
\end{tikzpicture}
\settowidth{\somelength}{(b)}
\setlength{\somelength}{\dimexpr\textwidth-\somelength\relax}
\caption{\parbox[t]{\somelength}{Using the quantum channel with memory $\operator{N}$ consecutively.}}
\label{fig:generic_memory}
\end{subfigure}
\begin{subfigure}[b]{0.96\columnwidth}\centering
\begin{tikzpicture}[nodes/.style={draw=none, inner sep=0pt, outer sep=0pt},
smallnode/.style={minimum size = 15pt, draw}]
\node[smallnode] (cq1) {$\mathcal{E}$};
\node[smallnode, below =40pt of cq1.center, anchor=center] (cq2)
\node[below =40pt of cq2.center, anchor=center] (cqdots)
\node[smallnode, below =40pt of cqdots.center, anchor = center] (cqn)
\node[smallnode, right = 150pt of cq1] (qc1) {$\mathcal{D}$};
\node[smallnode, right = 150pt of cq2] (qc2) {$\mathcal{D}$};
\node[below = 40pt of qc2.center, anchor = center] (qcdots) {$\vdots$};
\node[smallnode, right = 150pt of cqn] (qcn) {$\mathcal{D}$};

\path (cq1) edge[draw=none] node[smallnode] (T1) {$\operator{N}^{\boxtimes k}$} (qc1);
\path (cq1) edge node[above] {$\system{A}_1^k$} (T1);
\path (T1) edge node[above] {$\system{B}_1^k$} (qc1);

\path (cq2) edge[draw=none] node[smallnode] (T2) {$\operator{N}^{\boxtimes k}$} (qc2);
\path (cq2) edge node[above] {$\system{A}_{k+1}^{2k}$} (T2);
\path (T2) edge node[above] {$\system{B}_{k+1}^{2k}$} (qc2);
\path (cqdots) edge[draw=none] node (Tdots) {$\vdots$} (qcdots);
\path (cqn) edge[draw=none] node[smallnode] (Tn) {$\operator{N}^{\boxtimes k}$} (qcn);
\path (cqn) edge node[below] {$\system{A}_{(n\!-\!1)\cdot k+1}^{n\cdot k}$} (Tn);
\path (Tn) edge node[below] {$\system{B}_{(n\!-\!1)\cdot k+1}^{n\cdot k}$} (qcn);

\path (cq1.west) edge[line width=2pt] node[pos=1, above] {$\rv{X}_1$} ([xshift=-20pt]cq1.west);
\path (cq2.west) edge[line width=2pt] node[pos=1, above] {$\rv{X}_2$} ([xshift=-20pt]cq2.west);
\path (cqn.west) edge[line width=2pt] node[pos=1, above] {$\rv{X}_n$} ([xshift=-20pt]cqn.west);
\path (qc1.east) edge[line width=2pt] node[pos=1, above] {$\rv{Y}_1$} ([xshift=20pt]qc1.east);
\path (qc2.east) edge[line width=2pt] node[pos=1, above] {$\rv{Y}_2$} ([xshift=20pt]qc2.east);
\path (qcn.east) edge[line width=2pt] node[pos=1, above] {$\rv{Y}_n$} ([xshift=20pt]qcn.east);

\draw ([yshift=10pt]T1.north) -- (T1.north) node[pos=0, above=0pt] {$\system{S}_0$};
\path (T1) edge node[right] {$\system{S}_k$} (T2);
\path (T2) edge node[right] {$\system{S}_{2k}$} (Tdots);
\path (Tdots) edge node[right] {$\system{S}_{(n-1)\cdot k}$} (Tn);
\draw ([yshift=-10pt]Tn.south) -- (Tn.south) node[pos=0, below=0pt] {$\system{S}_{n\cdot k}$};
\end{tikzpicture}
\caption{Generic classical communication corresponding to~\eqref{eq:capacity:jointk}.}
\label{fig:generic_memory_multiple}
\end{subfigure}
\caption{Classical communications over quantum channels.}
\label{fig:ccq}
\end{figure}
\par
In this chapter, we are interested in computing and bounding the information rate as in~\eqref{eq:capacity:jointk} for finite-dimensional quantum channels with memory using only separable input ensembles and local output measurements, \ie, the case $k=1$, which is depicted in Figure~\ref{fig:ccq2}.
This restriction is equivalent to the scenario where no quantum computing device is present at the sending or receiving end or where our manipulation of the channel is limited to a single-channel use.
% The above sentence is exactly the same as a sentence in the preface.
The difficulty of the problem lies with the presence of the quantum memory.
In the simplest situation, the memory system exhibits classical properties under certain ensembles and measurements.
In this case, the resulting classical communication setup is equivalent to a finite-state-machine channel (FSMC)~\cite{gallager1968information}.
Though the evaluation of the information rate of an FSMC is nontrivial in general, efficient stochastic methods for estimating and bounding this quantity have been developed~\cite{arnold2006simulation, sadeghi2009optimization}.
\begin{figure}[t]\centering
\begin{tikzpicture}[nodes/.style={draw=none, inner sep=0pt, outer sep=0pt},
smallnode/.style={minimum size = 15pt, draw},
bignode/.style={minimum height = 150pt, minimum width = 20pt, draw}]
\node[bignode] (M) {}; \node[rotate=90] at (M) {Classical Encoder};
\node[bignode, right=200pt of M] (Mp) {}; \node[rotate=90] at (Mp) {Classical Decoder};
\path (M.north) edge[draw=none]
node[smallnode, anchor=north, pos=.25] (cq1) {$\mathcal{E}$}
node[smallnode, anchor=north, pos=.5] (T1) {$\operator{N}$}
node[smallnode, anchor=north, pos=.75] (qc1) {$\mathcal{D}$} (Mp.north);
\path (M.south) edge[draw=none]
node[smallnode, anchor=south, pos=.25] (cqn) {$\mathcal{E}$}
node[smallnode, anchor=south, pos=.5] (Tn) {$\operator{N}$}
node[smallnode, anchor=south, pos=.75] (qcn) {$\mathcal{D}$} (Mp.south);
\path (cq1.center) edge[draw=none]
node[smallnode, anchor=center, pos=.3333] (cq2) {$\mathcal{E}$}
node[anchor=center, pos=.6667] (cqdots) {$\vdots$} (cqn.center);
\path (T1.center) edge[draw=none]
node[smallnode, anchor=center, pos=.3333] (T2) {$\operator{N}$}
node[anchor=center, pos=.6667] (Tdots) {$\vdots$} (Tn.center);
\path (qc1.center) edge[draw=none]
node[smallnode, anchor=center, pos=.3333] (qc2) {$\mathcal{D}$}
node[anchor=center, pos=.6667] (qcdots) {$\vdots$} (qcn.center);
\path (M.east|-cq1) edge[line width=2pt] node[above] {$\rv{X}_1$} (cq1.west);
\path (M.east|-cq2) edge[line width=2pt] node[above] {$\rv{X}_2$} (cq2.west);
\path (M.east|-cqn) edge[line width=2pt] node[above] {$\rv{X}_n$} (cqn.west);
\path (cq1) edge node[above] {$\system{A}_1$} (T1);
\path (cq2) edge node[above] {$\system{A}_2$} (T2);
\path (cqn) edge node[above] {$\system{A}_n$} (Tn);
\path (T1) edge node[above] {$\system{B}_1$} (qc1);
\path (T2) edge node[above] {$\system{B}_2$} (qc2);
\path (Tn) edge node[above] {$\system{B}_n$} (qcn);
\path (qc1.east) edge[line width=2pt] node[above] {$\rv{Y}_1$} (Mp.west|-qc1);
\path (qc2.east) edge[line width=2pt] node[above] {$\rv{Y}_2$} (Mp.west|-qc2);
\path (qcn.east) edge[line width=2pt] node[above] {$\rv{Y}_n$} (Mp.west|-qcn);
\path (M.west) edge[line width=2pt] node[above, pos=1] {$\rv{U}$} ([xshift=-20pt]M.west);
\path (Mp.east) edge[line width=2pt] node[above, pos=1] {$\hat{\rv{U}}$} ([xshift=20pt]Mp.east);
\path ([yshift=10pt]T1.north) edge node[pos=0, above=0pt] {$\system{S}_0$} (T1.north);
\path (T1) edge node[right] {$\system{S}_1$} (T2);
\path (T2) edge node[right] {$\system{S}_2$} (Tdots);
\path (Tdots) edge node[right] {$\system{S}_n$} (Tn);
\path ([yshift=-10pt]Tn.south) edge node[pos=0, below=0pt] {$\system{S}_{n+1}$} (Tn.south);
\end{tikzpicture}
\caption{Classical communication over a quantum channel with memory using a separable ensemble and local measurements.}
\label{fig:ccq2}
\end{figure}
\par
Our work is highly inspired by~\cite{arnold2006simulation}, where the authors considered the information rate of FSMCs.
In particular, for an indecomposable FSMC~\cite{gallager1968information} with its channel law described by \pmf $W$, its information rate, which is independent of the initial channel state, is given by
\begin{equation}\label{eq:def:fsmc:ir:1}
\infoRate_W(Q) = \lim_{n\to\infty} \frac{1}{n} \mutualInfo(\rv{X}_1^n;\rv{Y}_1^n),
\end{equation}
where $\rv{X}_1^n = \left(\rv{X}_1,\ldots,\rv{X}_n\right)$ is the channel input process characterized by some sequence of distributions $\{Q^{(n)}\}_{n}$ and where $\rv{Y}_1^n = \left(\rv{Y}_1,\ldots,\rv{Y}_n\right)$ is the channel output process.
Although, except for very special cases, there are no single-letter or other simple expressions for information rate available, efficient stochastic techniques have been developed for estimating the information rate for \emph{stationary} and \emph{ergodic} input processes $\{Q^{(n)}\}_{n}$~\cite{arnold2006simulation, sharma2001entropy, pfister2001achievable}.
(For these techniques, under mild conditions, the numerical estimate of the information rate converges with probability one to the true value when the length of the channel input sequence goes to infinity.)
In this chapter, we extend such techniques to quantum channels with memory; in particular, we use \emph{factor graphs} for quantum probabilities~\cite{loeliger2017factor} (see Section~\ref{subsec:FG:QP}) for estimating quantities of interest.
These graphical models are useful for visualizing the relevant computations and for providing a clear comparison between the setup considered in this chapter and its classical counterparts in~\cite{arnold2006simulation} and~\cite{sadeghi2009optimization}.\footnote{
Clearly, the graphical models that we use are very similar to tensor networks (see, for example, the discussion in Appendix~\ref{app:figures}~A of~\cite{loeliger2017factor}).
A benefit of the graphical models that we use (including the corresponding terminology), is that they are compatible with the graphical models that are being used in classical information processing.}
\par
Our work is also partially inspired by~\cite{sadeghi2009optimization}, where the authors proposed upper and lower bounds based on some so-called auxiliary FSMCs, which are often lower-complexity approximations of the original FSMC.
They also provided efficient methods for optimizing these bounds.
Such techniques have been proven useful for FSMCs with large state spaces, when the above-mentioned information rate estimation techniques can be overly time-consuming.
Interestingly enough, the lower bounds represent achievable rates under mismatched decoding, where the decoder bases its computations not on the original FSMC but on the auxiliary FSMC~\cite{ganti2000mismatched}.
(See~\cite{sadeghi2009optimization} for a more detailed discussion of this topic and for further references.)
In this chapter, we also consider auxiliary channels and their induced bounds.
However, the auxiliary channels of our interest are chosen from a larger set of channels called \emph{quantum-state channels}, which will be defined in Section~\ref{QCwM:sec:3:QCM}.
We also propose a method for optimizing these bounds.
In particular, our method for optimizing the lower bound is ``data-driven'' in the sense that only  the input/output sequences of the original channel are needed, but not the mathematical model of the original channel.
\par
One must note that even if we can efficiently compute or bound the information rate, it is still a long way to go to compute the classical capacity of a quantum channel with memory.
On the one hand, maximizing $\infoRate(\mathcal{E},\operator{N},\mathcal{D})$ is a difficult problem.
(The analogous classical problems have been addressed in~\cite{arimoto1972algorithm},~\cite{blahut1972computation}, and~\cite{vontobel2008generalization}.)
On the other hand, due to the superadditivity property~\cite{hastings2009superadditivity} of quantum channels, which happens to be more common for quantum channels with memory~\cite{macchiavello2004transition, karimipour2006entanglement, lupo2010transitional} (compared with memoryless quantum channels), it is inevitable to consider joint ensembles on input systems and joint measurements on output systems across multiple channel uses.
\par
The rest of this chapter is organized as follows.
Section~\ref{sec:2:FSMC} reviews the method of estimating the information rate of an FSMC.
Section~\ref{QCwM:sec:3:QCM} models the classical communication scheme over a quantum channel with memory, and defines the notion of quantum-state channels as an equivalent description.
A graphical notation for representing such channels is also presented in this section.
Section~\ref{QCwM:sec:4:IR} estimates the information rate of such channels.
Section~\ref{QCwM:sec:5:UBLB} considers the upper and lower bounds induced by auxiliary quantum-state channels, and presents methods for optimizing them.
Section~\ref{QCwM:sec:6:example} contains numerical examples.
\section[Review of FSMCs]{Review of (Classical) Finite-State Machine Channels:
Information Rate, its Estimation, and Bounds}
\label{sec:2:FSMC}
In this section, we review the methods developed in~\cite{arnold2006simulation} for estimating the information rate of a (classical) FSMC, and the auxiliary-channel-induced upper and lower bounds studied in~\cite{sadeghi2009optimization}.
As we will see, the development in later sections about quantum channels will have many similarities, but also some important differences.
We emphasize that this section is a \emph{brief review} of~\cite{arnold2006simulation} and~\cite{sadeghi2009optimization} for the purpose of introducing necessary tools and ideas for later sections.
\subsection{Finite-State Machine Channels (FSMCs) and their Graphical Representation}
A (time-invariant) finite-state machine channel (FSMC) consists of an input alphabet $\set{X}$, an output alphabet $\set{Y}$, a state alphabet $\set{S}$, all of which are finite, and a channel law $W(y,s'|x,s)$, where the latter equals the probability of receiving $y\in\set{Y}$ and ending up in state $s'\in\set{S}$ given the channel input $x\in\set{X}$ and the previous channel state $s\in\set{S}$.
The relationship among the input, output, and state processes $\rv{X}_1^n,\rv{Y}_1^n,\rv{S}_0^n$ of $n$-channel uses can be described by the conditional \pmf
\begin{equation}\label{eq:FSMC:n:conditional}
W(\vy_1^n,\vs_1^n|\vx_1^n,s_0)  \defeq \prob_{\rv{Y}_1^n,\rv{S}_1^n|\rv{X}_1^n,\rv{S}_0}(\vy_1^n,\vs_1^n|\vx_1^n,s_0) = \prod_{\ell=1}^{n} W(y_\ell,s_\ell|x_\ell,s_{\ell-1}),
\end{equation}
where $x_\ell\in\set{X}$, $y_\ell\in\set{Y}$, and $s_\ell\in\set{S}$ for each $\ell$.
\begin{example}[Gilbert--Elliott channels]\label{example:GEC}
A notable class of examples of FSMCs are the Gilbert--Elliott channels~\cite{mushkin1989capacity}, which behave like a binary symmetric channel (BSC) with cross-over probability $p_s$ controlled by the channel state $s\in\{``\mathrm{b}",``\mathrm{g}"\}$, where usually $\bigl| \pbad - \frac{1}{2} \bigr| < \bigl| \pgood - \frac{1}{2} \bigr|$.
The state process itself is a first-order stationary ergodic Markov process that is independent of the input process.\footnote{
The independence of the state process on the input process is a particular feature of the Gilbert--Elliott channel.
In general, the state process of a finite-state channel can depend on the input process.}
(For more details, see, \eg, the discussions in~\cite{sadeghi2009optimization}.)
\end{example}
Given an input process $\{Q^{(n)}\}_{n}$ and an initial state \pmf $\prob_{\rv{S}_0}(s_0)$, we can write down the joint \pmf of  $(\rv{X}_1^n,\rv{Y}_1^n,\rv{S}_0^n)$ as
\begin{equation} \label{eq:FSMC:glabal:distribution}
g(\vx_1^n,\vy_1^n,\vs_0^n) \defeq \prob_{\rv{X}_1^n,\rv{Y}_1^n,\rv{S}_0^n}(\vx_1^n,\vy_1^n,\vs_0^n) = \prob_{\rv{S}_0}(s_0) \!\cdot\! Q^{(n)}(\vx_1^n) \!\cdot\!
\prod_{\ell=1}^{n} W(y_\ell,s_\ell|x_\ell,s_{\ell-1}).
\end{equation}
The factorization of $g(\vx_1^n,\vy_1^n,\vs_0^n)$ as shown in~\eqref{eq:FSMC:glabal:distribution} can be visualized with the help of an NFG as in Figure~\ref{fig:FMSC:high:level:1}.
\begin{figure}\centering
\begin{tikzpicture}[
factor/.style={rectangle, minimum width=1cm, minimum height=.7cm, draw},
sfactor/.style={rectangle, minimum size=.4cm, draw}]
\node[sfactor] (S) {$\prob_{\rv{S}_0}$};
\node[factor] (E1) [right=.7cm of S] {$W$};
\draw (S) -- (E1) node[above=-.05cm,midway] {$s_0$};
\node[sfactor] (X1) [above=1.1cm of E1] {$Q$};
\draw (X1) -- (E1) node[right=-.05cm, midway] {$x_1$};
\draw (E1.south) -- ([yshift=-.6cm]E1.south) node[right=-.05cm] {$y_1$};

\node[factor] (E2) [right=1cm of E1] {$W$};
\draw (E1) -- (E2) node[above=-.05cm,midway] {$s_1$};
\node[sfactor] (X2) [above=1.1cm of E2] {$Q$};
\draw (X2) -- (E2) node[right=-.05cm, midway] {$x_2$};
\draw (E2.south) -- ([yshift=-.6cm]E2.south) node[right=-.05cm] {$y_2$};

\node[factor,draw=none] (Edummy) [right=1cm of E2] {$\cdots$};
\node[sfactor,draw=none] (Xdummy) [above=1.1cm of Edummy] {$\cdots$};

\node[factor] (En) [right=1cm of Edummy] {$W$};
\draw (E2) -- (Edummy) node[above=-.05cm,midway] {$s_2$};

\draw (Edummy) -- (En) node[above=-.05cm,midway] {$s_{n-1}$};
\node[sfactor] (Xn) [above=1.1cm of En] {$Q$};
\draw (Xn) -- (En) node[right=-.05cm, midway] {$x_n$};
\draw (En.south) -- ([yshift=-.6cm]En.south) node[right=-.05cm] {$y_n$};

\node[sfactor, inner sep=0pt] (ee) [right=.5cm of En] {$1$};
\draw (En) -- (ee) node[above=-.05cm, midway] {$s_n$};

\begin{pgfonlayer}{bg}
\draw[dashed, blue, line width=1.5pt, fill=blue!10] ([xshift=-.8cm,yshift=0.45cm]X1) rectangle ([xshift=.7cm,yshift=-0.45cm]Xn);
\draw[dashed, red, line width=1.5pt, fill= red!10] ([xshift=-.7cm,yshift=0.55cm]S) rectangle ([xshift=.6cm,yshift=-0.55cm]ee);
\end{pgfonlayer}
\end{tikzpicture}
\caption[Channel with a classical state: closing the top box yields the input process $Q^{(n)}$, closing the bottom box yields the joint channel law $W(\vy_1^n|\vx_1^n)$.]{Channel with a classical state: closing the {\color{blue}top} box yields the input process $Q^{(n)}$, closing the {\color{red}bottom} box yields the joint channel law $W(\vy_1^n|\vx_1^n)$.}
\label{fig:FMSC:high:level:1}
\end{figure}
In particular:
\begin{enumerate}[label=\alph*)]
\item \label{closing:the:box} The part of the factor graph inside the {\color{red}bottom} box represents $W\left(\vy_1^n,\vs_1^n|\vx_1^n,s_0\right)$, \ie, the probability of obtaining $\vy_1^n$ and $\vs_1^n$ given $\vx_1^n$ and $s_0$.
After applying the ``closing-the-box'' operation (see Section~\ref{subsec:marginal:acyclic:FGs}) \wrt the {\color{red}bottom} box, we obtain the joint channel law $W(\vy_1^n|\vx_1^n) \defeq \sum_{\vs_0^n} \prob_{\rv{S}_0}(s_0) \cdot W(\vy_1^n,\vs_1^n|\vx_1^n,s_0) $.
\item The part of the factor graph inside the {\color{blue}top} box represents the input process $Q^{(n)}(\vx_1^n)$.
Here, for simplicity, the input process is an {i.i.d.}~process characterized by the \pmf $Q$, \ie, $Q^{(n)}(\vx_1^n) = \prod_{\ell=1}^{n} Q(x_{\ell})$.
\item The marginal function $g(\vx_1^n,\vy_1^n) \defeq \sum_{\vs_0^n} g(\vx_1^n,\vy_1^n,\vs_0^n)$is the marginal \pmf of $\vx_1^n$ and $\vy_1^n$.
The marginal function $g(\vs_0^n) \defeq \sum_{\vx_1^n,\vy_1^n} g(\vx_1^n,\vy_1^n, \vs_0^n)$ is the marginal \pmf of $\vs_0^n$.
Other marginal \pmfs can be obtained similarly.
\end{enumerate}
Using the ``closing-the-box'' operations, such factor graph representations can be useful in computing a number of quantities of interests.
For example, to prove that~\eqref{eq:FSMC:n:conditional} is indeed a valid conditional \pmf, it suffices to show that
\begin{equation}\label{eq:verfy:FSMC:conditional:distribution}
\sum_{\vs_1^n,\vy_1^n} W(\vy_1^n,\vs_1^n|\cvx_1^n,\cs_0) = 1
\quad \forall \cvx_1^n\in\set{X}^{n},\,\cs_0\in\set{S},
\end{equation}
which can be verified via a sequence of ``closing-the-box'' operations as shown in Figure~\ref{fig:FMSC:high:level:2} in the Appendix~\ref{app:figures}.
Such techniques are at the heart of the information-rate-estimation methods as in~\cite{arnold2006simulation}.
The details are reviewed in the next subsection.
\subsection{Information Rate Estimation} \label{QCwM:sec:FSMC:IR}
The approach of~\cite{arnold2006simulation} for estimating information rate of FSMCs, as reviewed in this section, is based on the Shannon--McMillan--Breiman theorem (see \eg,~\cite{cover2012elements}) and suitable generalizations.
We make the following assumptions.
\begin{itemize}
\item As already mentioned, the derivations in this chapter are for the case where the input process $\rv{X} = (\rv{X}_1,\rv{X}_2,\ldots)$ is an {i.i.d.}~process.
The results can be generalized to other stationary ergodic input processes that can be represented by a finite-state-machine source (FSMS).
Technically, this is done by defining a new state that combines the source state and the channel state.
\item We assume that the FSMC is indecomposable, which roughly means that in the long term the behavior of the channel is independent of the initial channel state distribution $\prob_{\rv{S}_0}$ (see~\cite[Section~4.6]{gallager1968information} for the exact definition).
For such channels and stationary ergodic input processes, the information rate $\infoRate_W$ in~\eqref{eq:def:fsmc:ir:1} is well defined.
\end{itemize}
Let $W(\vy_1^n|\vx_1^n)$ be the joint channel law of an FSMC satisfying the assumptions above.
As aforementioned, the information rate of such a channel using the {i.i.d.}~input distribution $\{Q^{(n)}\defeq Q^{\tensor{n}}\}_{n}$ is given by~\eqref{eq:def:fsmc:ir:1}, where the input process $\rv{X}_1^n$ and the output process $\rv{Y}_1^n$ are jointly distributed according to
\begin{equation} \label{eq:joint:fsmc}
\prob_{\rv{X}_1^n,\rv{Y}_1^n}(\vx_1^n,\vy_1^n) = \prod_{\ell=1}^{n} Q(x_\ell) \cdot W(\vy_1^n|\vx_1^n).
\end{equation}
One can rewrite~\eqref{eq:def:fsmc:ir:1} as
\begin{equation}\label{eq:def:fsmc:ir:2}
\infoRate_W(Q) = \entropicRate(\rv{X}) + \entropicRate(\rv{Y}) - \entropicRate(\rv{X},\rv{Y}),
\end{equation}
where the \emph{entropic rates} $\entropicRate(\rv{X})$, $\entropicRate(\rv{Y})$ and $\entropicRate(\rv{X},\rv{Y})$ are defined as
\begin{align}
\entropicRate(\rv{X}) &\defeq \lim_{n\to\infty} \frac{1}{n} \entropy(\rv{X}_1^n), \\
\entropicRate(\rv{Y}) &\defeq \lim_{n\to\infty} \frac{1}{n} \entropy(\rv{Y}_1^n), \\
\entropicRate(\rv{X},\rv{Y}) &\defeq \lim_{n\to\infty} \frac{1}{n} \entropy(\rv{X}_1^n,\rv{Y}_1^n).
\end{align}
\par
We proceed as in~\cite{arnold2006simulation}.
(For more background information, see the references in~\cite{arnold2006simulation}, in particular~\cite{ephraim2002hidden}.)
Namely, because of~\eqref{eq:def:fsmc:ir:2} and
\begin{align}
\label{eq:converge:x:1}
-\frac{1}{n}\log{\prob_{\rv{X}_1^n}(\rv{X}_1^n)} &\stackrel{n\to\infty}{\longrightarrow} \entropicRate(\rv{X})
&& \quad \text{in probability},\\
\label{eq:converge:y:1}
-\frac{1}{n}\log{\prob_{\rv{Y}_1^n}(\rv{Y}_1^n)} &\stackrel{n\to\infty}{\longrightarrow} \entropicRate(\rv{Y})
&& \quad \text{in probability},\\
\label{eq:converge:xy:1}
-\frac{1}{n}\log{\prob_{\rv{X}_1^n,\rv{Y}_1^n}(\rv{X}_1^n,\rv{Y}_1^n)} &\stackrel{n\to\infty}{\longrightarrow} \entropicRate(\rv{X},\rv{Y})
&& \quad \text{in probability},
\end{align}
by choosing some large number $n$, we have the approximation
\begin{equation}\label{eq:FSMC:ir:estimate:1}
\infoRate_W(Q) \approx -\frac{1}{n}\log{\prob_{\rv{X}_1^n}(\cvx_1^n)} - \frac{1}{n}\log{\prob_{\rv{Y}_1^n}(\cvy_1^n)} + \frac{1}{n}\log{\prob_{\rv{X}_1^n,\rv{Y}_1^n}(\cvx_1^n,\cvy_1^n)}
\end{equation}
where $\cvx_1^n$ and $\cvy_1^n$ are some input and output sequences, respectively, randomly generated according to
\begin{equation}\label{eq:fsmc:joint:xy:distribution}
\prob_{\rv{X}_1^n,\rv{Y}_1^n}(\cvx_1^n,\cvy_1^n) = \sum_{\vs_0^n} \prob_{\rv{S}_0}(s_0)\cdot Q^{(n)}(\cvx_1^n) \cdot W(\cvy_1^n,\vs_1^n|\cvx_1^n,s_0),
\end{equation}
where $W(\cvy_1^n,\vs_1^n|\cvx_1^n,s_0)$ is defined in~\eqref{eq:FSMC:n:conditional}.
Note that $\cvx_1^n$ can be obtained by simulating the input process, and $\cvy_1^n$ can be obtained by simulating the channel for the given input string $\cvx_1^n$.
The latter can be done by keeping track of $\prob_{\rv{Y}_\ell|\rv{X}_1^\ell,\rv{Y}_1^{\ell\!-\!1}}(y_\ell|\cvx_1^\ell,\cvy_1^{\ell\!-\!1})$, which is proportional to $\prob_{\rv{Y}_\ell,\rv{Y}_1^{\ell\!-\!1}|\rv{X}_1^\ell}(y_\ell,\cvy_1^{\ell\!-\!1}|\cvx_1^\ell)$, and can be efficiently calculated by applying suitable ``closing-the-box'' operations as in Figure~\ref{fig:CFSM:channel:simulation:Y} in the Appendix~\ref{app:figures}.
\par
We continue by showing how the three terms appearing on the right-hand side of~\eqref{eq:FSMC:ir:estimate:1} can be computed efficiently.
We show it explicitly for the second term, and then outline it for the first and the third term. \par
In order to efficiently compute the second term on the right-hand side of~\eqref{eq:FSMC:ir:estimate:1}, \ie, $-\frac{1}{n}\log{\prob_{\rv{Y}_1^n}(\cvy_1^n)}$, we consider the \emph{state metric} defined in~\cite{arnold2006simulation} as
\begin{equation}\label{eq:def:classical:channel:state:metric:Y:1}
\muY_{\ell}(s_{\ell}) \defeq \sum_{\vx_1^{\ell}} \sum_{\vs_0^{\ell\!-\!1}} \prob_{\rv{S}_0}(s_0) \cdot Q^{(\ell)}(\vx_1^\ell) \cdot W(\cvy_1^\ell,\vs_1^\ell|\vx_1^\ell,s_0).
\end{equation}
In this case,
\begin{equation} \label{eq:calculate:p:y:1}
\prob_{\rv{Y}_1^n}(\cvy_1^n) = \sum_{s_n} \muY_n(s_n),
\end{equation}
and the calculation of $\muY_\ell(s_{\ell})$ can be done iteratively as
\begin{align}
\muY_\ell(s_{\ell}) & = \sum_{x_{\ell}} \sum_{s_{\ell\!-\!1}} \muY_{\ell\!-\!1}(s_{\ell\!-\!1}) \cdot Q(x_{\ell}|\vx_1^{\ell\!-\!1}) \cdot W(\cy_{\ell},s_{\ell}|x_{\ell},s_{\ell\!-\!1}) \\
\label{eq:recursive:state:metric:Y:1}
& = \sum_{x_{\ell}} \sum_{s_{\ell\!-\!1}} \muY_{\ell\!-\!1}(s_{\ell\!-\!1}) \cdot Q(x_{\ell}) \cdot W(\cy_{\ell},s_{\ell}|x_{\ell},s_{\ell\!-\!1}).
\end{align}
Eq.~\eqref{eq:recursive:state:metric:Y:1} is visualized in Figure~\ref{fig:CFSM:estimate:hY} as applying suitable ``closing-the-box'' operations to the NFG in Figure~\ref{fig:FMSC:high:level:1}. \par
However, since the value of $\muY_{\ell}(s_\ell)$ tends to zero as $\ell$ grows, such recursive calculations are numerically inconvenient.
A solution is to \emph{normalize} $\muY_{\ell}(s_\ell)$ after each use of~\eqref{eq:recursive:state:metric:Y:1} and to keep track of the scaling coefficients.
\begin{equation}\label{eq:recursive:state:metric:Y:2}
\bmuY_{\ell}(s_\ell) \defeq \frac{1}{\lambdaY_\ell} \sum_{x_\ell}\sum_{s_{\ell\!-\!1}} \bmuY_{\ell\!-\!1}(s_\ell) \cdot Q(x_\ell) \cdot W(\cy_\ell,s_\ell|x_\ell,s_{\ell\!-\!1}),
\end{equation}
where the scaling factor $\lambdaY_\ell > 0$ is chosen such that $\sum_{s_\ell} \bmuY_\ell(s_\ell) = 1$.
With this, Eq.~\eqref{eq:calculate:p:y:1} can be rewritten as
\begin{equation} \label{eq:calculate:p:y:2}
\prob_{\rv{Y}_1^n}(\cvy_1^n) = \prod_{\ell=1}^n \lambdaY_{\ell}.
\end{equation}
Finally, we arrive at the following efficient procedure for computing $-\frac{1}{n}\log{\prob_{\rv{Y}_1^n}(\vy_1^n)}$:
\begin{itemize}
\item For $\ell = 1,\ldots, n$, iteratively compute the normalized state
metric and with that the scaling factors $\lambdaY_{\ell}$.
\item Conclude with the result
\begin{equation} \label{eq:calculate:p:y:3}
-\frac{1}{n}\log{\prob_{\rv{Y}_1^n}(\vy_1^n)} = \frac{1}{n} \sum_{\ell=1}^n
\log(\lambdaY_\ell).
\end{equation}
\end{itemize} \par
The third term on the right-hand side of~\eqref{eq:FSMC:ir:estimate:1} can be evaluated by an analogous procedure, where the state metric $\muY_\ell(s_\ell)$ is replaced by the state metric
\begin{equation} \label{eq:def:classical:channel:state:metric:XY:1}
\muXY_\ell(s_\ell) \defeq \sum_{\vs_0^{\ell\!-\!1}} \prob_{\rv{S}_0}(s_0) \cdot Q^{(\ell)}(\cvx_1^\ell) \cdot W(\cvy_1^\ell,\vs_1^\ell|\cvx_1^\ell).
\end{equation}
The iterative calculation of $\muXY_{\ell}(s_{\ell})$ is visualized in Figure~\ref{fig:CFSM:estimate:hXY}. \par
Finally, the first term on the right-hand side of~\eqref{eq:FSMC:ir:estimate:1} can be trivially evaluated if $\rv{X}$ is an {i.i.d.}~process, and with a similar approach as above if it is described by an FSMS.
\par
The above discussion is summarized as Algorithm~\ref{alg:SPA}.
On the side, note that for each $\ell=2,\ldots,n$, the quantities $\lambdaY_\ell$ and $\lambdaXY_\ell$ in the algorithm are the conditional probabilities $\prob_{\rv{Y}_\ell|\rv{Y}_{1}^{\ell\!-\!1}}(\cy_\ell|\cvy_1^{\ell\!-\!1})$ and $\prob_{\rv{X}_\ell\rv{Y}_\ell|\rv{X}_{1}^{\ell\!-\!1}\rv{Y}_{1}^{\ell\!-\!1}}(\cx_{\ell},\cy_\ell|\cvx_1^{\ell\!-\!1},\cvy_1^{\ell\!-\!1})$, respectively.
\begin{algorithm}
\caption{Estimating the Information Rate of an FSMC}
\begin{algorithmic}[1]
\Require{An indecomposable FSMC channel law $W$, a input distribution $Q$, a positive integer $n$ large enough.}
\Ensure{$\infoRate_W(Q) \approx \entropy(\rv{X}) + \hat\entropicRate(\rv{Y}) - \hat\entropicRate(\rv{X},\rv{Y})$.}
\State Initialize the channel state distribution $\prob_{\rv{S}_0}$ as a uniform distribution over $\set{S}$
\State Generate an input sequence $\cvx_1^n \sim Q^{\tensor n}$
\State Generate a corresponding output sequence $\cvy_1^n$
\State $\bmuY_{0}\gets \prob_{\rv{S}_0}$
\ForEach{$\ell=1,\ldots,n$}
\State $\muY_\ell(s_\ell) \gets \sum_{x_\ell,s_{\ell\!-\!1}} \bmuY_{\ell\!-\!1}(s_{\ell\!-\!1}) \cdot Q(x_\ell) \cdot W(\cy_\ell,s_\ell|x_\ell,s_{\ell\!-\!1})$
\State $\lambdaY_\ell \gets \sum_{s_\ell} \muY_\ell(s_\ell)$;
\State $\bmuY_\ell \gets \muY_\ell / \lambdaY_\ell$
\EndFor
\State $\hat\entropicRate(\rv{Y}) \gets -\frac{1}{n} \sum_{\ell=1}^n \log(\lambdaY_\ell)$
\State $\bmuXY_{0} \gets \prob_{\rv{S}_0}$
\ForEach{$\ell=1,\ldots,n$}
\State $\muXY_\ell(s_\ell) \gets \sum_{s_{\ell\!-\!1}} \bmuXY_{\ell\!-\!1}(s_{\ell\!-\!1}) \cdot Q(\cx_\ell) \cdot W(\cy_\ell,s_\ell|\cx_\ell,s_{\ell\!-\!1})$
\State $\lambdaXY_\ell \gets \sum_{s_\ell} \muXY_\ell(s_\ell)$
\State $\bmuXY_\ell \gets \muXY_\ell/\lambdaXY_\ell$
\EndFor
\State $\hat\entropicRate(\rv{X},\rv{Y})\gets - \frac{1}{n} \sum_{\ell=1}^n \log(\lambdaXY_\ell)$
\State $\entropy(\rv{X}) \gets -\sum_{x} Q(x) \log{Q(x)}$
\State Estimate $\infoRate_W(Q)$ as $\entropy(\rv{X}) + \hat\entropicRate(\rv{Y}) - \hat\entropicRate(\rv{X},\rv{Y})$.
\end{algorithmic}
\label{alg:SPA}
\end{algorithm}
\subsection{Auxiliary Channels and Bounds on the Information Rate}
\label{QCwM:sec:aux}
As mentioned earlier in this chapter, auxiliary channels\footnote{Technically speaking, an auxiliary channel can be defined as \emph{any} channel with the same input/output alphabet. For example, an auxiliary channel for an FSMC can be just another FSMC with smaller state space; in contrast, in Section~\ref{QCwM:sec:5:UBLB}, an auxiliary channel can also be a quantum-state channel.} are introduced when the state space of the FSMC is too large, making the calculation in Algorithm~\ref{alg:SPA} (practically) intractable.
More precisely, given an auxiliary forward FSMC (AF-FSMC) $\hat{W}(y_\ell,\hat{s}_\ell|x_\ell,\hat{s}_{\ell\!-\!1})$ and an auxiliary backward FSMC (AB-FSMC) $\hat{V}(x_\ell,\hat{s}_\ell|y_\ell,\hat{s}_{\ell\!-\!1})$,  a pair of upper and lower bounds of the information rate is given in~\cite{arnold2006simulation, sadeghi2009optimization} as
\begin{align}
\label{eq:def:IRUB}
\IRUB^{(n)}_{W}(\hat{W}) &\defeq \frac{1}{n} \sum_{\vx_1^n,\vy_1^n}
Q(\vx_1^n) W(\vy_1^n|\vx_1^n) \log{\frac{W(\vy_1^n|\vx_1^n)}{\QWaux(\vy_1^n)}},\\
\label{eq:def:IRLB}
\IRLB_{W}^{(n)}(\hat{V})  &\defeq \frac{1}{n} \sum_{\vx_1^n,\vy_1^n} Q(\vx_1^n) W(\vy_1^n|\vx_1^n) \log{\frac{\hat{V}(\vx_1^n|\vy_1^n)}{Q(\vx_1^n)}},
\end{align}
where $\QWaux(\vy_1^n) \defeq \sum_{\vx_1^n} Q(\vx_1^n) \cdot \hat{W}(\vy_1^n|\vx_1^n)$.
To see that~\eqref{eq:def:IRUB} and~\eqref{eq:def:IRLB} are, respectively, upper and lower bounds, one can verify that
\begin{align}
\label{eq:IRUB:minus:IR}
\IRUB_{W}(\hat{W}) - \infoRate_{W} &= \frac{1}{n} \infdiv{\QW(\rv{Y}_1^n)}{\QWaux(\rv{y}_1^n)},\\
\label{eq:IR:minus:IRLB}
\infoRate_{W} - \IRLB_{W}(\hat{V}) &= \frac{1}{n} \sum_{\vy_1^n} \QW(\vy_1^n) \cdot \infdiv{V(\rv{X}_1^n|\vy_1^n)}{\hat{V}(\rv{X}_1^n|\vy_1^n)},
\end{align}
where the backward channel $V(\vx|\vy)$ is defined as $V(\vx|\vy)\defeq Q(\vx)W(\vy|\vx)/\QW(\vy)$.
In particular, given an AF-FSMC $\hat{W}$,~\cite{sadeghi2009optimization} considered the induced AB-FSMC $\hat{V}(\vx|\vy)\defeq Q(\vx)\hat{W}(\vy|\vx)/\QWaux(\vy)$.
In this case,
\begin{equation}
\IRLB_{W}^{(n)}(\hat{V}) = \frac{1}{n} \sum_{\vx_1^n,\vy_1^n} Q(\vx_1^n) W(\vy_1^n|\vx_1^n) \log\frac{\hat{W}(\vy_1^n|\vx_1^n)}{\QWaux(\vy_1^n)}.
\end{equation}
The difference function $\Delta_{W}^{(n)}(\hat{W})$ is defined as
\begin{equation} \label{eq:delta} \begin{aligned}
\Delta_{W}^{(n)}(\hat{W}) &\defeq \IRUB^{(n)}_W(\hat{W}) - \IRLB_{W}^{(n)}(\hat{V})\\
&=\frac{1}{n} \sum_{\vx_1^n,\vy_1^n} Q(\vx_1^n) W(\vy_1^n|\vx_1^n) \log{\frac{W(\vy_1^n|\vx_1^n)}{\hat{W}(\vy_1^n|\vx_1^n)}}\\
&=\frac{1}{n} \infdiv{Q(\rv{X}_1^n) W(\rv{Y}_1^n|\rv{X}_1^n)}{Q(\rv{X}_1^n) \hat{W}(\rv{Y}_1^n|\rv{X}_1^n)}.
\end{aligned}\end{equation}
Apparently, $\Delta_{W}^{(n)}(\hat{W})\geqslant 0$, and equality holds if and only if $\hat{W}(\vy_1^n|\vx_1^n) = W(\vy_1^n|\vx_1^n)$ for all $\vx_1^n$ and $\vy_1^n$ with positive support \wrt $\prob_{\rv{X}_1^n,\rv{Y}_1^n}$ as in~\eqref{eq:joint:fsmc}.
An efficient algorithm for finding a local minimum of the difference function was proposed in~\cite{sadeghi2009optimization}; we refer to~\cite{sadeghi2009optimization} for further details.
\section[QCwM \& their Graphical Representation]{Quantum Channel with Memory and their Graphical Representation} \label{QCwM:sec:3:QCM}
In this section, we formalize our notations and modeling of quantum channels with memory~\cite{bowen2004quantum, kretschmann2005quantum, caruso2014quantum} and of classical communications over such channels.
In particular, we will define a class of channels named \emph{quantum-state channels},  which is an alternative description of the classical communications over quantum channels with memory.
In addition, we will introduce several NFGs for representing these channels and processes.
\subsection{Classical Communication over a Quantum Channel with Memory} \label{QCwM:sec:classical:comm:quantum}
As mentioned earlier in this chapter, we define a quantum channel with memory as follows.
\begin{definition} \label{def:quantum:channel:memory}
A \emph{quantum channel with memory} is a CPTP map:
\begin{equation} \label{eq:def:qcm}
\operator{N}:\DensOp(\hilbert_\system{A}\tensor\hilbert_\system{S}) \to \DensOp(\hilbert_\system{B}\tensor\hilbert_{\system{S}'}),
\end{equation}
where $\system{A}$ is the input system, $\system{B}$ is the output system, $\system{S}$ and $\system{S}'$ are, respectively, the memory systems before and after the channel use.
The Hilbert spaces $\hilbert_\system{A}$, $\hilbert_\system{B}$, and $\hilbert_\system{S} = \hilbert_{\system{S}'}$ are the state spaces corresponding to those systems.
\end{definition}
\par
We consider classical communication over such channels using some separable input ensemble and local output measurements; namely, the encoder and decoder are, respectively, some classical-to-quantum and quantum-to-classical channels involving a single input or output system.
In particular, given an ensemble $\{\rho^{(x)}_\system{A}\}_{x\in\set{X}}$ and a measurement $\{\Lambda_\system{B}^{(y)}\}_{y\in\set{Y}}$, we define the encoding and decoding function, respectively, as
\begin{align}
\text{Encoding }\mathcal{E}: \prob_\rv{X} &\mapsto \sum_{x\in\set{X}} \prob_\rv{X}(x) \rho^{(x)}_\system{A}
&& \forall\ \prob_\rv{X} \text{ \pmf over }\set{X},\\
\text{Decoding }\mathcal{D}: \sigma_\system{B} &\mapsto \left\{\tr(\Lambda_\system{B}^{(y)} \cdot \sigma_\system{B})\right\}_{y\in\set{Y}}
&& \forall\ \sigma_\system{B}\in\PositiveOp(\hilbert_\system{B}).
\end{align}
We emphasize that in our setup, the ensemble $\{\rho^{(x)}_\system{A}\}_{x\in\set{X}}$ and measurements $\{\Lambda_\system{B}^{(y)}\}_{y\in\set{Y}}$ are given and fixed.
Furthermore, we assume that one does not have access to the memory systems of the channel.
For the case of {i.i.d.}~inputs, the memory system $\system{S}$ before each channel use shall be independent of the input system $\system{A}$, namely, the joint memory-input operator shall take the form of $\rho_\system{A}\tensor\rho_\system{S}$ at each channel input.\footnote{
More generally, for FSMSs, this statement also holds by conditioning on all previous inputs.}
\par
With this, the probability of receiving $y\in\set{Y}$, given that $x\in\set{X}$  was sent and given that the density operator of the memory system \emph{before} the usage of the channel was $\rho_\system{S}$, equals
\begin{equation}\label{eq:channel:law:1}
\prob_{\rv{Y}|\rv{X};\system{S}}(y|x;\rho_\system{S}) = \tr\left( \Lambda_\system{B}^{(y)} \cdot \tr_{\system{S}'}\left( \operator{N}(\rho^{(x)}_\system{A}\tensor \rho_\system{S}) \right) \right),
\end{equation}
which can also be written as
\begin{equation}\label{eq:channel:law:2}
\prob_{\rv{Y}|\rv{X};\system{S}}(y|x;\rho_\system{S}) = \tr\left( (\Lambda_\system{B}^{(y)}\tensor I_\system{S}) \cdot \operator{N}( \rho^{(x)}_\system{A}\tensor\rho_\system{S}) \right).
\end{equation}
Moreover, assuming that $y$ was observed, the density operator of the
memory system \emph{after} the channel use is given by
\begin{equation} \label{eq:channel:evolution:1}
\rho_{\system{S}'} = \frac{\tr_\system{B}\left( (\Lambda_\system{B}^{(y)} \tensor I_\system{S}) \cdot \operator{N}(\rho^{(x)}_\system{A} \tensor \rho_\system{S})\right)}{\tr\left( (\Lambda_\system{B}^{(y)} \tensor I_\system{S}) \cdot \operator{N}(\rho^{(x)}_\system{A} \tensor \rho_\system{S}) \right)}.
\end{equation}
Notice that the denominator in~\eqref{eq:channel:evolution:1} equals the expressions in~\eqref{eq:channel:law:1} and~\eqref{eq:channel:law:2}.
One should note that, though the input and the memory systems are independent before each channel use (given {i.i.d.}~inputs), the output and the memory systems after each channel use can be correlated or even entangled.
In particular, this translates to the fact that the measurement outcome $y$ can have an influence on the memory system as indicated in~\eqref{eq:channel:evolution:1}.
\par
Consider using the channel $n$ times consecutively with the above scheme.
The joint channel law, namely the conditional \pmf of the channel outputs $\rv{Y}_1^n$ given the channel inputs $\rv{X}_1^n$ and the initial channel state $\rho_{\system{S}_0}$, can be computed iteratively using~\eqref{eq:channel:law:2} and~\eqref{eq:channel:evolution:1}.
In particular, the joint conditional \pmf can be computed as
\begin{equation}\label{eq:joint:1}
\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0}) = \prod_{\ell=1}^n \prob_{\rv{Y}_\ell|\rv{X}_\ell;\system{S}_{\ell\!-\!1}}
\end{equation}
where we compute the density operators $\{\rho_{\system{S}_\ell}\}_{\ell=1}^n$ iteratively using~\eqref{eq:channel:evolution:1} as
\begin{equation}\label{eq:channel:evolution:2}
\rho_{\system{S}_\ell} = \frac{\tr_\system{B}\left( (\Lambda_\system{B}^{(y_\ell)}\tensor I_\system{S}) \cdot \operator{N}(\rho^{(x_\ell)}_\system{A}\tensor\rho_{\system{S}_{\ell\!-\!1}}) \right)}{\tr\left( (\Lambda_\system{B}^{(y_\ell)}\tensor I_\system{S}) \cdot \operator{N}(\rho^{(x_\ell)}_\system{A} \tensor \rho_{\system{S}_{\ell\!-\!1}}) \right)}.
\end{equation}
\subsection{Quantum-State Channels} \label{QCwM:sec:QSC}
For each channel-ensemble-measurement configuration ($\operator{N}$, $\{\rho^{(x)}_\system{A}\}_{x\in\set{X}}$, $\{\Lambda_\system{B}^{(y)}\}_{y\in\set{Y}}$) as introduced above, one ends up with a joint conditional \pmf, as  in~\eqref{eq:joint:1}.
However, this relationship is not bijective.
In particular, consider some unitary operators $U_\system{A}$ and $U_\system{B}$ acting on $\hilbert_\system{A}$ and $\hilbert_\system{B}$, respectively.
The following setup induces exactly the same joint conditional \pmf:
\begin{align*}
&\tilde{\operator{N}}: \tilde\rho_\system{AS} \mapsto (U_\system{B} \tensor I_\system{S}) \cdot \operator{N}\left( (U_\system{A} \tensor I_\system{S}) \tilde\rho_\system{AS} (U_\system{A}^\Herm \tensor I_\system{S}) \right) \cdot (U_\system{B}^\Herm \tensor I_\system{S}),\\
&\tilde\rho^{(x)}_\system{A} \defeq U_\system{A}^\Herm \cdot \rho^{(x)}_\system{A} \cdot U_\system{A} &&\quad\forall x\in\set{X},\\
&\tilde\Lambda_\system{B}^{(y)} \defeq U_\system{B}^\Herm \cdot \Lambda_\system{B}^{(y)} \cdot U_\system{B} &&\quad\forall y\in\set{Y}.
\end{align*}
Such redundancy is not only tedious, but also detrimental when we try to compare different channels; in particular, when we try to introduce proper auxiliary channels to approximate the original communication scheme.
\par
In this subsection, we introduce a class of channels called~\emph{quantum-state channels} to eliminate such redundancies.
In particular, notice that the statistical behavior of the aforementioned communication scheme is fully specified via~\eqref{eq:channel:law:2} and~\eqref{eq:channel:evolution:1}; which are in turn determined by the set of completely positive mappings $\{\operator{N}^{y|x}\}_{x\in\set{X},y\in\set{Y}}$ defined as
\begin{equation}\label{eq:def:qsc}
\operator{N}^{y|x}: \rho_\system{S} \mapsto \tr_\system{B}\left( (\Lambda_\system{B}^{(y)} \tensor I_\system{S}) \cdot \operator{N}(\rho^{(x)}_\system{A} \tensor \rho_\system{S}) \right).
\end{equation}
In this case,~\eqref{eq:channel:law:2},~\eqref{eq:channel:evolution:1}, and~\eqref{eq:joint:1} can be rewritten, respectively, as
\begin{align}
\label{eq:channel:law:3}
\prob_{\rv{Y}|\rv{X};\system{S}}(y|x;\rho_\system{S}) &= \tr\left( \operator{N}^{y|x}(\rho_\system{S}) \right),\\
\label{eq:channel:evolution:3}
\rho_{\system{S}'} &= \operator{N}^{y|x}(\rho_\system{S}) \big/ \tr\left( \operator{N}^{y|x}(\rho_\system{S}) \right),\\
\label{eq:joint:2}
\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0} &= \tr\left( \operator{N}^{y_n|x_n} \circ \cdots \circ \operator{N}^{y_1|x_1}(\rho_{\system{S}_0}) \right).
\end{align}
Thus, the operators $\{\operator{N}^{y|x}\}_{x\in\set{X},y\in\set{Y}}$ fully specify the joint conditional \pmf as in~\eqref{eq:joint:2}.
Moreover, such specification is also \emph{unique}; namely, any two sets of channel-ensemble-measurement configuration shall end up with the same joint channel law if and only if the mappings defined in~\eqref{eq:def:qsc} are identical.
This inspires us to make the following definition.
\begin{definition}[Quantum-State Channel] \index{quantum-state channel}
A (finite indexed) set of completely positive operators $\{\operator{N}^{y|x}\}_{x\in\set{X},y\in\set{Y}}$ (acting on the same Hilbert space) is said to be a  \emph{(classical-input classical-output) quantum-state channel} (CC-QSC) if $\sum_{y\in\set{Y}} \operator{N}^{y|x}$ is trace-preserving for each $x \in \set{X}$.
\end{definition}
Given any channel-ensemble-measurement configuration as described in Section~\ref{QCwM:sec:classical:comm:quantum}, one can always define a corresponding CC-QSC by~\eqref{eq:def:qsc}.
On the other hand, as stated in the proposition below, the converse is also true.
\begin{proposition}\label{prop:quantum:state:channel}
For any CC-QSC $\{\operator{N}^{y|x}\}_{x\in\set{X},y\in\set{Y}}$, there exists some quantum channel with memory $\operator{N}$ as in~\eqref{eq:def:qcm} such that~\eqref{eq:def:qsc} holds with the ensemble $\{\rho_\system{A}^{(x)} = \braket{x}\}_{x\in\set{X}}$ and the measurement $\{\Lambda_\system{B}^{(y)} = \braket{y}\}_{y\in\set{Y}}$.
Here, $\hilbert_\system{A}$ and $\hilbert_\system{B}$ are defined such that $\{\bra{x}\}_x$ and $\{\bra{y}\}_y$ are orthonormal bases of $\hilbert_\system{A}$ and $\hilbert_\system{B}$, respectively.
\end{proposition}
\begin{proof}
It suffices to show that there exists a CPTP map $\operator{N}: \DensOp(\hilbert_\system{A} \tensor \hilbert_\system{S}) \to \DensOp(\hilbert_\system{B} \tensor \hilbert_\system{S})$ such that for all $\rho_\system{S} \in \DensOp(\hilbert_\system{S})$, and $x\in\set{X}$,
\[
\operator{N}: \braket{x}\tensor\rho_\system{S} \mapsto \sum_{y\in\set{Y}} \braket{y}\tensor\operator{N}^{y|x}(\rho_\system{S}).
\]
Such an $\operator{N}$ can be constructed as
\[
\operator{N}: \rho \mapsto \sum_{x,y,k}\left(\bra{y}\!\ket{x} \otimes E^{y|x}_k\right) \cdot \rho \cdot \left( \bra{y}\!\ket{x} \otimes E^{y|x}_k\right)^{\Herm},
\]
where $\left\{E^{y|x}_k\right\}_{k}$ is a Kraus representation of $\operator{N}^{y|x}$.
It remains to check if $\operator{N}$ is a CPTP, which is indeed the case:
\begin{align*}
\sum_{x,y,k} \left(\bra{y}\!\ket{x} \otimes E^{y|x}_k\right)^{\Herm} \cdot \left(\bra{y}\!\ket{x} \otimes E^{y|x}_k\right)
&= \sum_x \sum_{y,k} \braket{x} \otimes (E^{y|x}_k)^\Herm E^{y|x}_k \\
&= \sum_x \braket{x} \otimes I = I. \qedhere
\end{align*}
\end{proof}
\subsection{Visualization using Normal Factor Graphs} \label{QCwM:subsec:NFGs}
In this subsection, we focus on the computations of~\eqref{eq:channel:law:3},~\eqref{eq:channel:evolution:3}, and~\eqref{eq:joint:2} for the situation where the involved channel $\operator{N}$ is of finite dimension.
In analogy to the FSMCs, we demonstrate how to use NFGs to facilitate and visualize the relevant computations (see Section~\ref{subsec:FG:QP} and~\cite{loeliger2017factor}).
\par
By Proposition~\ref{prop:quantum:state:channel}, let us consider a CC-QSC $\{\operator{N}^{y|x}\}_{x\in\set{X},y\in\set{Y}}$ acting on $\hilbert_\system{S}$, where $d=\dim(\hilbert_\system{S})$ is finite, and $\{\bra{s}\}_{s\in\set{S}}$ is an orthonormal basis of $\hilbert_\system{S}$.
(Apparently, $\size{\set{S}}=d$.)
Since for each $x$ and $y$, $\operator{N}^{y|x}$ is a completely positive map, there must exist finitely many (\emph{not} necessarily unique) matrices $\{ F_k^{y|x} \in \mathbb{C}^{\set{S}\times\set{S}} \}_{k}$ such that
\begin{equation} \label{eq:Kraus:quantum:state:channel}
\bigl[\operator{N}^{y|x}(\rho_\system{S})\bigr] = \sum_{k} F_k^{y|x} \cdot [\rho_\system{S}] \cdot (F_k^{y|x})^\Herm
\quad \forall \rho_\system{S} \in \DensOp(\hilbert_\system{S}),
\end{equation}
where $\bigl[\operator{N}^{y|x}(\rho_\system{S})\bigr]$ and $[\rho_\system{S}]$ are, respectively, the matrix representation of the operator $\operator{N}^{y|x}(\rho_\system{S})$ and $\rho_\system{S}$ under $\{\bra{s}\}_{s\in\set{S}}$.
The reason for such matrices $\{F_k^{y|x}\}_{k}$ to exist is the same as for the Kraus operators of CPTP maps (see Theorem~\ref{thm:Kraus}).
Also note that $\sum_{y\in\set{Y}} \mathcal{E}^{y|x}$ is trace-preserving, thus it must hold that
\begin{equation} \label{eq:operator:sum:representation:condition:1}
\sum_{y\in\set{Y}}\sum_{k} (F_k^{y|x})^\Herm F_k^{y|x} = I
\quad \forall x \in\set{X}.
\end{equation}
Now, define a set of functions $\{W^{y|x}\}_{x\in\set{X},y\in\set{Y}}$ as
\begin{equation}\label{eq:def:channel:function:representation}
W^{y|x}:(s',s,\tilde{s}',\tilde{s}) \mapsto \sum_k F_k^{y|x}(s',s) \conj{F_k^{y|x}(\tilde{s}',\tilde{s})},
\end{equation}
where $s',s,\tilde{s}',\tilde{s}\in\set{S}$ are indices of the corresponding matrices, namely, $F_k^{y|x}(s',s)$ is the $(s',s)$-th entry of the matrix $F_k^{y|x}$.
In this case, one can rewrite~\eqref{eq:channel:law:3},~\eqref{eq:channel:evolution:3} and~\eqref{eq:joint:2}, respectively, into
\begin{align}
\label{eq:channel:law:4}
\prob_{\rv{Y}|\rv{X};\system{S}}(y|x;\rho_\system{S})
&= \sum_{s',\tilde{s}':\:\atop s'=\tilde{s}'} \sum_{s,\tilde{s}} W^{y|x}(s',s,\tilde{s}',\tilde{s})\cdot[\rho_\system{S}]_{s,\tilde{s}}, \\
\label{eq:channel:evolution:4}
&= \frac{\sum_{s,\tilde{s}} W^{y|x}(s',s,\tilde{s}',\tilde{s}) \cdot [\rho_\system{S}]_{s,\tilde{s}}}{\sum_{s',\tilde{s}':\:\atop s'=\tilde{s}'}\sum_{s,\tilde{s}} W^{y|x}(s',s,\tilde{s}',\tilde{s}) \cdot [\rho_\system{S}]_{s,\tilde{s}}}, \\
\label{eq:joint:3}
\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})
&= \sum_{s_n,\tilde{s}_n:\:\atop s_n=\tilde{s}_n} \sum_{\vs_0^{n\!-\!1},\tilde{\vs}_0^{n\!-\!1}} [\rho_{\system{S}_0}]_{s_0,\tilde{s}_0} \cdot \prod_{\ell=1}^n W^{y_\ell|x_\ell}(s_\ell,s_{\ell\!-\!1},\tilde{s}_\ell,\tilde{s}_{\ell\!-\!1}).
\end{align}
By rearranging the entries of $W^{y|x}$ (for each $x,y$) into a matrix $[W^{y|x}]\in\mathbb{C}^{\set{S}^2\times\set{S}^2}$ as
\begin{equation}\label{eq:qsc:matrix:1}
[W^{y|x}]_{(s',\tilde{s}'),(s,\tilde{s})} \defeq W^{y|x}(s',s,\tilde{s}',\tilde{s}),
\end{equation}
where $(s',\tilde{s}')\in\set{S}^2$ is the first index, and $(s,\tilde{s})\in\set{S}^2$ is the second index of $[W^{y|x}]$, we can simplify~\eqref{eq:channel:law:4},~\eqref{eq:channel:evolution:4}, and~\eqref{eq:joint:3} as
\begin{align}
\label{eq:channel:law:5}
\prob_{\rv{Y}|\rv{X};\system{S}}(y|x;\rho_\system{S})
&= \tr([W^{y|x}] \cdot [\rho_\system{S}]),\\
\label{eq:channel:evolution:5}
&= \frac{[W^{y|x}] \cdot [\rho_\system{S}]}{\tr([W^{y|x}] \cdot [\rho_\system{S}])},\\
\label{eq:joint:4}
\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})
&= \tr\left( [W^{y_{n}|x_{n}}] \cdots [W^{y_{1}|x_{1}}] \cdot [\rho_{\system{S}_0}] \right),
\end{align}
Here we treat $[\rho_\system{S}]$ as a length-$d^2$ vector indexed by $(s,\tilde{s})\in\set{S}^2$ in the above equations.
\par
By considering $\{W^{y|x}\}_{x,y}$ as a function of six variables, we can represent it using a factor vertex of degree six in an NFG as in Figure~\ref{fig:NFG:QSC:single}.
In this case, Eqs.~\eqref{eq:channel:law:4} and~\eqref{eq:channel:law:5} can be visualized as ``closing the {\color{red} outer} box'' in the factor graph.
Similarly,~\eqref{eq:channel:evolution:4} and~\eqref{eq:channel:evolution:5} can be visualized as ``closing the {\color{blue}inner} box''.
The factor graph corresponding to using the channel $n$~times consecutively is depicted in Figure~\ref{fig:NFG:QSC:multiple}, where~\eqref{eq:joint:3} and~\eqref{eq:joint:4} are visualized as closing the {\color{red}outermost} box.
Interestingly, this ``closing-the-box'' operation can be carried out by a sequence of simpler ``closing-the-box'' operations as shown in the figure.
\begin{figure}\centering
\begin{tikzpicture}[node/.style={draw=none},
factor/.style={rectangle, minimum width=1cm, minimum height=1.7cm, draw},
sfactor/.style={rectangle, minimum size=.4cm, draw}]
\node[factor] (W) {$W$};
\node[above left=-.5cm and 1.2cm of W] (sp) {};
\node[below left=-.5cm and 1.2cm of W] (sp') {};
\node[above right=-.5cm and 1.2cm of W] (s) {};
\node[below right=-.5cm and 1.2cm of W] (s') {};
\node[above=.7cm of W] (X) {$x$};
\node[below=.7cm of W] (Y) {$y$};
\node[sfactor, inner sep=0pt, midway] (EW) [right=1.6cm of W] {$=$};
\node[sfactor, inner sep=0pt, midway] (rho) [left=1.6cm of W] {};
\node at (rho) {$\rho_\system{S}$};
\draw (sp.east-|W.west) -| (rho) node[above=-0.05cm,pos=0.25] {$s$};
\draw (sp'.east-|W.west) -| (rho) node[below,pos=0.25] {$\tilde{s}$};
\draw (sp.east-|W.east) -| (EW) node[above=-0.05cm,pos=0.25] {$s'$};
\draw (sp'.east-|W.east) -| (EW) node[below,pos=0.25] {$\tilde{s}'$};
\draw (X) -- (W);
\draw (Y) -- (W);
\begin{pgfonlayer}{bg}
\draw[dashed, red, line width=1.5pt,fill= red!30]
([xshift=-.7cm,yshift=1.2cm]rho) rectangle
\draw[dashed, blue, line width=1.5pt,fill= blue!30]
([xshift=-.4cm,yshift=1cm]rho) rectangle
\end{pgfonlayer}
\end{tikzpicture}
\caption{Representation of $\{W^{y|x}\}_{x,y}$ using an NFG.}
\label{fig:NFG:QSC:single}
\end{figure}
\begin{figure}\centering
\begin{tikzpicture}[
factor/.style ={rectangle, minimum width=1cm, minimum height=1.7cm, draw},
sfactor/.style={rectangle, minimum size=.4cm, draw},
label/.style={anchor=south east, circle, draw, inner sep=.5pt,
outer sep=5pt, font=\scriptsize}]
\node[sfactor] (S) {}; \node at (S) {\resizebox{.375cm}{!}{$\rho_{\system{S}_0}$}};
\node[factor] (E1) [right=.5cm of S] {$W$};
\draw (S.north) |- ([yshift=.6cm]E1.west) node[above=-0.05cm,pos=.75] {$s_0$};
\draw (S.south) |- ([yshift=-.6cm]E1.west) node[below,pos=.75] {$\tilde{s}_0$};
\draw (E1.north) -- ([yshift=1cm]E1.north) node[right] {$x_1$};
\draw (E1.south) -- ([yshift=-1cm]E1.south) node[right] {$y_1$};
\node[factor] (E2) [right=.6cm of E1] {$W$};
\draw ([yshift=.6cm]E1.east) |- ([yshift=.6cm]E2.west) node[above=-0.05cm,pos=.75] {$s_1$};
\draw ([yshift=-.6cm]E1.east) |- ([yshift=-.6cm]E2.west) node[below,pos=.75] {$\tilde{s}_1$};
\draw (E2.north) -- ([yshift=1cm]E2.north) node[right] {$x_2$};
\draw (E2.south) -- ([yshift=-1cm]E2.south) node[right] {$y_2$};
\node[factor, draw=none] (Edummy) [right=.6cm of E2] {$\cdots$};
\draw ([yshift=.6cm]E2.east) |- ([yshift=.6cm]Edummy.west) node[above=-0.05cm,pos=.75] {$s_2$};
\draw ([yshift=-.6cm]E2.east) |- ([yshift=-.6cm]Edummy.west) node[below,pos=.75] {$\tilde{s}_2$};
\node[factor] (En) [right=.6cm of Edummy] {$W$};
\draw ([yshift=.6cm]Edummy.east) |- ([yshift=.6cm]En.west) node[above=-0.05cm,pos=.75] {$s_{\!n-\!1}$};
\draw ([yshift=-.6cm]Edummy.east) |- ([yshift=-.6cm]En.west) node[below,pos=.75] {$\tilde{s}_{n\!-\!1}$};
\draw (En.north) -- ([yshift=1cm]En.north) node[right] {$x_n$};
\draw (En.south) -- ([yshift=-1cm]En.south) node[right] {$y_n$};
\node[sfactor, draw=none] (Edummy1) [above=-0.50cm of Edummy] {$\cdots$};
\node[sfactor, draw=none] (Edummy2) [below=-0.45cm of Edummy] {$\cdots$};
\node[sfactor,draw=none] (Xdummy) at ([yshift=1cm]E1.north-|Edummy) {$\cdots$};
\node[sfactor,draw=none] (Ydummy) at ([yshift=-1cm]E1.south-|Edummy) {$\cdots$};
\node[sfactor, inner sep=0pt] (ee) [right=.5cm of En] {$=$};
\draw ([yshift=.6cm]En.east) -| (ee.north) node[above=-0.05cm,pos=.25] {$s_n$};
\draw ([yshift=-.6cm]En.east) -| (ee.south) node[below,pos=.25] {$\tilde{s}_n$};
\begin{pgfonlayer}{bg}
\draw[dashed, red,line width=1.5pt,fill= red!10] ([xshift=-.7cm,yshift=1.5cm]S) rectangle ([xshift=1.5cm,yshift=-1.6cm]En);
\draw[dashed, red,line width=1.5pt,fill= red!20] ([xshift=-.6cm,yshift=1.4cm]S) rectangle	([xshift=.6cm,yshift=-1.5cm]En);
\node[label] at ([xshift=.6cm,yshift=-1.5cm]En) {$n$};
\draw[dashed, red,line width=1.5pt,fill= red!30] ([xshift=-.5cm,yshift=1.3cm]S) rectangle ([xshift=.6cm,yshift=-1.4cm]E2);
\node[label] at ([xshift=.6cm,yshift=-1.4cm]E2) {2};
\draw[dashed, red,line width=1.5pt,fill= red!50] ([xshift=-.4cm,yshift=1.2cm]S) rectangle ([xshift=.6cm,yshift=-1.3cm]E1);
\node[label] at ([xshift=.6cm,yshift=-1.3cm]E1) {1};
\end{pgfonlayer}
\end{tikzpicture}
\caption[The joint channel law~\eqref{eq:joint:3} and~\eqref{eq:joint:4} can be visualized as the result of the ``closing of the outermost box'' above, which can in turn be carried out by a sequence of ``closing-the-box'' operations as indicated.]{The joint channel law~\eqref{eq:joint:3} and~\eqref{eq:joint:4} can be visualized as the result of the ``closing of the {\color{blue}outermost} box'' above, which can in turn be carried out by a sequence of ``closing-the-box'' operations as indicated.}
\label{fig:NFG:QSC:multiple}
\end{figure}
\par
A number of statistical quantities and density operators of interest can be computed and visualized as ``closing-the-box'' operations on suitable NFGs similar to that of Figure~\ref{fig:NFG:QSC:multiple}.
The following example highlights how quantities of this kind can be computed in such a manner.
\begin{figure}[t]\centering
\resizebox{\columnwidth}{!}{\begin{tikzpicture}[
factor/.style={rectangle, minimum width=1cm, minimum height=1.7cm, draw},
sfactor/.style={rectangle, minimum size=.5cm, draw},
darksolid/.style={rectangle, minimum size=.15cm, draw, fill = black,
inner sep=0pt, outer sep = 0pt}]
\node[sfactor] (S) {}; \node at (S) {$\rho_{\system{S}_0}$};
\node[factor] (E1) [right=.7cm of S] {$W$};
\draw (S.north) |- ([yshift=.6cm]E1.west) node[above=-0.05cm,pos=.75] {$s_0$};
\draw (S.south) |- ([yshift=-.6cm]E1.west) node[below,pos=.75] {$\tilde{s}_0$};
\node[sfactor] (X1) [above=.7cm of E1] {}; \node at (X1) {$Q$};
\draw (X1) -- (E1) node[right, midway] {$x_1$};
\draw (E1.south) -- ([yshift=-.8cm]E1.south) node[darksolid]{} node (Y1) [right] {$\cy_1$};
\node[factor, draw=none] (Edummy1) [right=1cm of E1] {};
\node at ([yshift=.6cm]E1.east-|Edummy1) {$\cdots$};
\node at ([yshift=-.6cm]E1.east-|Edummy1) {$\cdots$};
\draw ([yshift=.6cm]E1.east) |- ([yshift=.6cm]Edummy1.west) node[above=-0.05cm,pos=.75] {$s_1$};
\draw ([yshift=-.6cm]E1.east) |- ([yshift=-.6cm]Edummy1.west) node[below,pos=.75] {$\tilde{s}_1$};
\node at (X1-|Edummy1) {$\cdots$};
\node at (Y1-|Edummy1) {$\cdots$};
\node[factor] (El0) [right=1cm of Edummy1] {$W$};
\draw ([yshift=.6cm]Edummy1.east) |- ([yshift=.6cm]El0.west) node[above=-0.05cm,pos=.75] {$s_{\ell-2}$};
\draw ([yshift=-.6cm]Edummy1.east) |- ([yshift=-.6cm]El0.west) node[below,pos=.75] {$\tilde{s}_{\ell-2}$};
\node[sfactor] (Xl0) [above=.7cm of El0] {}; \node at (Xl0) {$Q$};
\draw (Xl0) -- (El0) node[right, midway] {$x_{\ell\!-\!1}$};
\draw (El0.south) -- ([yshift=-.8cm]El0.south) node[darksolid]{} node[right] {$\cy_{\ell\!-\!1}$};
\node[factor] (El) [right=1cm of El0] {$W$};
\draw ([yshift=.6cm]El0.east) |- ([yshift=.6cm]El.west) node[above=-0.05cm,pos=.85] {$s_{\ell\!-\!1}$};
\draw ([yshift=-.6cm]El0.east) |- ([yshift=-.6cm]El.west) node[below,pos=.85] {$\tilde{s}_{\ell\!-\!1}$};
\node[sfactor] (Xl) [above=.7cm of El] {}; \node at (Xl) {$Q$};
\draw (El.north) -- (Xl.south) node[darksolid,midway,red,fill=red]{} node[right,midway] {$x_{\ell}$};
\draw (El.south) -- ([yshift=-.8cm]El.south) node[darksolid]{} node[right] {$\cy_{\ell}$};
\node[factor] (El2) [right=1cm of El] {$W$};
\draw ([yshift=.6cm]El.east) |- ([yshift=.6cm]El2.west) node[above=-0.05cm,pos=.75] {$s_\ell$};
\draw ([yshift=-.6cm]El.east) |- ([yshift=-.6cm]El2.west) node[below,pos=.75] {$s_\ell'$};
\node[sfactor] (Xl2) [above=.7cm of El2] {}; \node at (Xl2) {$Q$};
\draw (Xl2) -- (El2) node[right, midway] {$x_{\ell\!+\!1}$};
\draw (El2.south) -- ([yshift=-.8cm]El2.south) node[darksolid]{} node[right] {$\cy_{\ell\!+\!1}$};
\node[factor, draw=none] (Edummy2) [right=1cm of El2] {};
\node at ([yshift=.6cm]El2.east-|Edummy2) {$\cdots$};
\node at ([yshift=-.6cm]El2.east-|Edummy2) {$\cdots$};
\draw ([yshift=.6cm]El2.east) |- ([yshift=.6cm]Edummy2.west) node[above=-0.05cm,pos=.75] {$s_{\ell\!+\!1}$};
\draw ([yshift=-.6cm]El2.east) |- ([yshift=-.6cm]Edummy2.west) node[below,pos=.75] {$s_{\ell\!+\!1}'$};
\node at (X1-|Edummy2) {$\cdots$};
\node at (Y1-|Edummy2) {$\cdots$};
\node[factor] (En) [right=1cm of Edummy2] {$W$};
\draw ([yshift=.6cm]Edummy2.east) |- ([yshift=.6cm]En.west) node[above=-0.05cm,pos=.75] {$s_{n\!-\!1}$};
\draw ([yshift=-.6cm]Edummy2.east) |- ([yshift=-.6cm]En.west) node[below,pos=.75] {$s_{n\!-\!1}'$};
\node[sfactor] (Xn) [above=.7cm of En] {}; \node at (Xn) {$Q$};
\draw (Xn) -- (En) node[right, midway] {$x_n$};
\draw (En.south) -- ([yshift=-.8cm]En.south) node[darksolid]{} node[right] {$\cy_n$};
\node[sfactor, inner sep=0pt] (ee) [right=.5cm of En] {$=$};
\draw ([yshift=.6cm]En.east) -| (ee.north) node[above=-0.05cm,pos=.25] {$s_n$};
\draw ([yshift=-.6cm]En.east) -| (ee.south) node[below,pos=.25] {$s_n'$};
\begin{pgfonlayer}{bg}
\draw[dashed, black, line width=1.5pt, fill=yellow!20] ([xshift=-.8cm,yshift=2.4cm]S) rectangle ([xshift=.8cm,yshift=-2.3cm]ee);
\draw[dashed, blue, line width=1.5pt, fill=blue!20] ([xshift=-.7cm,yshift=2.2cm]El2) rectangle ([xshift=.6cm,yshift=-2.1cm]ee);
\node[anchor=south east, fill= blue!50] at ([xshift=.6cm,yshift=-2.1cm]ee) {$\lvec{\sigma}_{\system{S}_{\ell}}^{(\cvy_{\ell\!+\!1}^n)}$};
\draw[dashed, red, line width=1.5pt, fill=red!20] ([xshift=-.6cm,yshift=2.2cm]S) rectangle ([xshift=.8cm,yshift=-2.1cm]El0);
\node at ([xshift=-.6cm,yshift=2.2cm]S) (tempA) {};
\node at ([xshift=.8cm,yshift=-2.1cm]El0) (tempB) {};
\node[anchor=south west,fill=red!50] at (tempA |- tempB) {$\rvec{\sigma}_{\system{S}_{\ell\!-\!1}}^{(\cvy_1^{\ell\!-\!1})}$};
\end{pgfonlayer}
\end{tikzpicture}}
\caption{Computation of the marginal \pmf $\prob_{\rv{X}_\ell|\rv{Y}_1^n;\system{S}_0}$ using a sequence of ``closing-the-box'' operations.}
\label{fig:qsc:bcjr}
\end{figure}
\begin{example}[BCJR~\cite{bahl1974optimal} decoding for CC-QSCs]
For fixed $\cvy_1^n\in\set{Y}^n$ and a given initial density operator $\rho_{\system{S}_0}$, the conditional probability $\prob_{\rv{X}_\ell|\rv{Y}_1^n;\system{S}_0}(x_\ell|\cvy_1^n;\rho_{\system{S}_0})$ can be computed via
\begin{equation}\label{eq:qsc:bcjr:1}
\prob_{\rv{X}_\ell|\rv{Y}_1^n;\system{S}_0}(\cdot|\cvy_1^n;\rho_{\system{S}_0}) \propto \prob_{\rv{X}_\ell,\rv{Y}_1^n|\system{S}_0}(\cdot,\cvy_1^n|\rho_{\system{S}_0}),
\end{equation}
where the right-hand side of~\eqref{eq:qsc:bcjr:1} is a marginal \pmf defined as
\begin{equation}\label{eq:qsc:bcjr:2}
\prob_{\rv{X}_\ell,\rv{Y}_1^n|\system{S}_0}(x_\ell,\cvy_1^n|\rho_{\system{S}_0}) = \sum_{\vx_{1}^{\ell\!-\!1},\vx_{\ell\!+\!1}^{n}} \sum_{\vs_0^n,\tilde{\vs}_0^n} [\rho_{\system{S}_0}]_{s_0,\tilde{s}_0} \cdot \prod_{i=1}^n Q(x_i) \cdot \prod_{j=1}^n W^{\cy_j|x_j}(s_j,s_{j\!-\!1},\tilde{s}_j,\tilde{s}_{j\!-\!1}),
\end{equation}
where we have assumed that the input process $\rv{X}_1^n$ is {i.i.d.}~characterized by some \pmf $Q$.
The evaluation of~\eqref{eq:qsc:bcjr:2} can be carried out efficiently using a sequence of ``closing-the-box'' operations as visualized in Figure~\ref{fig:qsc:bcjr}.
These operations can be roughly divided into the following three steps:
\begin{enumerate}
\item Closing the {\color{red}left inner} box: this results in an operator $\rvec{\sigma}_{\system{S}_{\ell\!-\!1}}^{(\cvy_1^{\ell\!-\!1})}$ on $\hilbert_{\system{S}_{\ell\!-\!1}}$.
\item Closing the {\color{blue}right inner} box: this results in another operator $\lvec{\sigma}_{\system{S}_{\ell}}^{(\cvy_{\ell\!+\!1}^n)}$ on $\hilbert_{\system{S}_{\ell}}$.
\item Applying the ``closing-the-box'' operation to the yellow box: the result is the marginal \pmf $\prob_{\rv{X}_\ell,\rv{Y}_1^n|\system{S}_0}(x_\ell,\cvy_1^n|\rho_{\system{S}_0})$, from which the desired conditional probability $\prob_{\rv{X}_\ell|\rv{Y}_1^n;\system{S}_0}$ $(x_\ell|\cvy_1^n;\rho_{\system{S}_0})$ can be easily obtained by normalization.
\end{enumerate}
The operators mentioned in 1) and 2) can be computed recursively, using a sequence of ``closing-the-box'' operations.
Namely, one can carry out the computations in 1) consecutively with $\ell=1,2,\ldots,n$; and the computations in 2) consecutively with $\ell=n,n\!-\!1,\ldots,1$.
This provides an efficient way to evaluate $\prob_{\rv{X}_\ell|\rv{Y}_1^n;\system{S}_0}(x_\ell|\cvy_1^n;\rho_{\system{S}_0})$ for each $\ell=1,\ldots,n$; and thus provides an efficient symbol-wise decoding algorithm.
The idea in this example is conceptually identical to that of the BCJR decoding algorithm for an FSMC.
\end{example}
As shown in the above example, very often the desired functions or quantities are based on the same partial results.
The NFG framework is very helpful to visualize these partial results and to show how they can be combined to obtain the desired functions and quantities.
\par
We emphasize that the functions $\{W^{y|x}\}_{x,y}$ defined in~\eqref{eq:def:channel:function:representation} are unique for a given finite-dimensional CC-QSC $\{\operator{N}^{y|x}\}_{x,y}$; even though such uniqueness does not apply to the Kraus operators $\{F^{y|x}\}_k$ being used to define $\{W^{y|x}\}_{x,y}$.
This can be proven by making the identification that
\begin{equation}\label{eq:qsc:matrix:identification}
\bigl[\operator{N}^{y|x}(\rho_\system{S})\bigr] = [W^{y|x}]\cdot [\rho_\system{S}]
\qquad \forall \rho_\system{S}\in\DensOp(\hilbert_\system{S})
\end{equation}
for all $x$ and $y$.
Moreover, we argue that the functions $\{W^{y|x}\}_{x,y}$, are an \emph{equivalent} way to specify a CC-QSC, or classical communication over a quantum channel with memory as described at the beginning of this section.
Namely, for any set of complex-valued functions $\{W^{y|x}\}_{x,y}$ on $\set{S}^{4}$ satisfying some constraints to be clarified later, there must exist a unique CC-QSC $\{\operator{N}^{y|x}\}_{x,y}$ such that~\eqref{eq:qsc:matrix:identification} holds; and thus, there must exist some corresponding channel-ensemble-measurement configuration, unique up to its channel law.
As for such constraints, we rearrange the entries of $W^{y|x}$ (for each $x,y$) into another matrix $\llbracket W^{y|x}\rrbracket \in \mathbb{C}^{\set{S}^2 \times \set{S}^2}$ (\aka Choi--Jamio{\l}kowski matrix~\cite{jamiolkowski1972linear}), whose entries are defined as
\begin{equation}\label{eq:qsc:matrix:2}
\llbracket W^{y|x}\rrbracket_{(s',s),(\tilde{s}',\tilde{s})} \defeq W^{y|x}(s',s,\tilde{s}',\tilde{s}),
\end{equation}
where $(s',s)\in\set{S}^2$ is the first index, and $(\tilde{s}',\tilde{s})\in\set{S}^2$ is the second index of $\llbracket W^{y|x}\rrbracket$.
Notice that, $\llbracket W^{y|x}\rrbracket$ is a PSD matrix, and satisfies the equation
\begin{equation} \label{eq:operator:sum:representation:condition:2}
\sum_{y\in\set{Y}} \sum_{s',\tilde{s}':\: s'=\tilde{s}'} \llbracket W^{y|x}\rrbracket_{(s',s),(\tilde{s}',\tilde{s})} = \delta_{s,\tilde{s}}
\quad \forall x\in\set{X}.
\end{equation}
In this case, the ``equivalence'' can be shown by the following proposition.
\begin{proposition} \label{prop:CFR}
Let $\set{X}$, $\set{Y}$ be finite sets, and $\hilbert_\system{S}$ be a finite-dimensional Hilbert space with an orthonormal basis $\{\bra{s}\}_{s\in\set{S}}$.
For any set of functions
\[
\{W^{y|x}:\set{S}\times\set{S}\times\set{S}\times\set{S} \to \mathbb{C}\}_{x\in\set{X},y\in\set{Y}}
\]
such that their matrix form $\{\llbracket W^{y|x}\rrbracket\}_{x,y}$ consists of PSD matrices and satisfies~\eqref{eq:operator:sum:representation:condition:2}, there must exist a unique CC-QSC $\{\operator{N}^{y|x}\}_{x,y}$ acting on $\hilbert_\system{S}$ such that~\eqref{eq:qsc:matrix:identification} holds.
\end{proposition}
\begin{proof}
The idea of the proof is to consider the eigenvalue decomposition of $\llbracket W^{y|x}\rrbracket$, and reconstruct $\operator{N}^{y|x}$ by following the equations~\eqref{eq:def:channel:function:representation} and~\eqref{eq:Kraus:quantum:state:channel} backwardly.
We omit the details here.
\end{proof}
\par
Let us conclude this section by pointing out that the functions $\{W^{y|x}\}_{x,y}$, particularly the corresponding NFG, can be constructed from the channel-ensemble-measurement configuration ($\operator{N}$, $\{\rho^{(x)}_\system{A}\}_{x\in\set{X}}$, $\{\Lambda_\system{B}^{(y)}\}_{y\in\set{Y}}$) as in Figure~\ref{fig:CFR}.
This can be justified by checking~\eqref{eq:def:qsc} and~\eqref{eq:qsc:matrix:identification}.
\begin{figure}\centering
\begin{tikzpicture}[
node/.style={draw=none, inner sep=0pt, outer sep=0pt},
factor/.style={rectangle, minimum width=1cm, minimum height=3cm, draw},
sfactor/.style={rectangle, minimum size=.75cm, draw, inner sep=1pt}]
\node[sfactor] (rhoX) {};
\node at (rhoX) {$\rho_\system{A}^{(x)}$};
\node[factor, right=1cm of rhoX] (N) {$\operator{N}$};
\draw[*-] (rhoX.north) |- ([yshift=1cm]N.west);
\draw[-*] (rhoX.south) |- ([yshift=-1cm]N.west);
\node at ([xshift=-.5cm,yshift=2.3cm]rhoX) (X) {$x$};
\draw (X) |- (rhoX);
\draw ([xshift=-3.50cm,yshift=1.25cm]N.west) -- ([yshift=1.25cm]N.west) node[above=1pt,pos=0,anchor=south west] {$s$};
\draw[-*] ([xshift=-3.50cm,yshift=-1.25cm]N.west) -- ([yshift=-1.25cm]N.west) node[below=1pt,pos=0,anchor=north west] {$\tilde{s}$};
\draw[-*] ([xshift=3.50cm,yshift=1.25cm]N.east) -- ([yshift=1.25cm]N.east) node[above=1pt,pos=0,anchor=south east] {$s'$};
\draw ([xshift=3.50cm,yshift=-1.25cm]N.east) -- ([yshift=-1.25cm]N.east) node[below=1pt,pos=0,anchor=north east] {$\tilde{s}'$};
\node[sfactor, right=1cm of N] (M) {};
\node at (M) {$\Lambda^{(y)}_{\system{B}}$};
\draw[-*] (M.north) |- ([yshift=1cm]N.east);
\draw[*-] (M.south) |- ([yshift=-1cm]N.east);
\node at ([xshift=.5cm,yshift=-2.3cm]M) (Y) {$y$};
\draw (Y) |- (M);
\begin{pgfonlayer}{bg}
\draw[dashed, red, line width=1.5pt, fill= red!20] ([xshift=-1.0cm,yshift=1.7cm]rhoX) rectangle ([xshift=1cm,yshift=-1.7cm]M);
\node[red,anchor=north west] at ([xshift=-1.0cm,yshift=-1.7cm]rhoX) {$W^{y|x}(s',s,\tilde{s}',\tilde{s})$};
\end{pgfonlayer}
\end{tikzpicture}
\caption{NFG representation of the channel-ensemble-measurement configuration ($\operator{N}$, $\{\rho^{(x)}_\system{A}\}_{x\in\set{X}}$, $\{\Lambda_\system{B}^{(y)}\}_{y\in\set{Y}}$).}
\label{fig:CFR}
\end{figure}
\section{Information Rate and its Estimation}\label{QCwM:sec:4:IR}
In this section, we focus on the information rate of the communication scheme described in Section~\ref{QCwM:sec:3:QCM}.
As defined in~\eqref{eq:capacity:jointk}, the information rate is the limit superior of the average mutual information $\frac{1}{n}\mutualInfo\left(\rv{X}_1^n;\rv{Y}_1^n\right)$ between the input and output processes $\rv{X}_1^n$ and $\rv{Y}_1^n$ as $n$ tends to infinity.
We assume that $\rv{X}_1^n$ is distributed according to some {i.i.d.}~process\footnote{
For more general type of sources, like a finite-state-machine source (FSMS), one can consider ``merging'' the memory of the source into that of the channel, and thus obtaining an equivalent memoryless input process.}
characterized by the \pmf $Q$, \ie, $Q^{(n)}(\vx_1^n) = \prod_{\ell=1}^{n} Q(x_{\ell})$.
In this case, the joint distribution of $(\rv{X}_1^n,\rv{Y}_1^n)$ is given by
\begin{equation}\label{eq:joint:5}
\prob_{\rv{X}_1^n,\rv{Y}_1^n|\system{S}_0}(\vx_1^n,\vy_1^n|\rho_{\system{S}_0}) =  \prod_{\ell=1}^n Q(x_\ell) \cdot \prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0}),
\end{equation}
where $\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}$ is specified in~\eqref{eq:joint:1},~\eqref{eq:joint:2},~\eqref{eq:joint:3} or~\eqref{eq:joint:4}, depending on which notation we use to specify the channel (see Propositions~\ref{prop:quantum:state:channel} and~\ref{prop:CFR}).
It is obvious that the value of~\eqref{eq:joint:5}, and thus the information rate, depends on the initial density operator $\rho_{\system{S}_0}$.
In this sense, we denote the information rate as a function of the input \pmf $Q$, the CC-QSC $\{\operator{N}^{y|x}\}_{x,y}$ describing the channel, and the initial density operator $\rho_{\system{S}_0}$, namely
\begin{align}
\label{eq:def:qsc:IR:2}
\infoRate(Q,\{\operator{N}^{y|x}\}_{x,y},\rho_{\system{S}_0})
&\defeq \limsup_{n\to\infty} \infoRate^{(n)}(Q,\{\operator{N}^{y|x}\}_{x,y},\rho_{\system{S}_0}),\\
\label{eq:def:qsc:IR:1}
\infoRate^{(n)}(Q,\{\operator{N}^{y|x}\}_{x,y},\rho_{\system{S}_0})
&\defeq \frac{1}{n}\mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)(\rho_{\system{S}_0}).
\end{align}
Here, $\mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)(\rho_{\system{S}_0})$ is the mutual information between $\rv{X}_1^n$ and $\rv{Y}_1^n$; and the latter are jointly distributed according to~\eqref{eq:joint:5}.
The argument $\rho_{\system{S}_0}$ emphasizes the dependency of $\frac{1}{n}\mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)$ on $\rho_{\system{S}_0}$.
\par
Similar to the case of an FSMC, the dependency of the information rate on the initial density operator usually cannot be ignored.
However, as already mentioned in~Section~\ref{QCwM:sec:FSMC:IR}, for a class of FSMCs, namely the indecomposable FSMCs, it is known that the information rate is independent of the initial channel state~\cite{gallager1968information}.
An indecomposable FSMC, intuitively speaking, is an FSMC whose state distribution, given different initial states, tends to be indistinguishable as $n\to\infty$, independently of the input sequence realized.
A quantum analogy was proposed by Bowen, Devetak, and Mancini~\cite{bowen2005bounds}, where they defined the indecomposable quantum channels with memory, and proved that the quantum entropic bound for such channels is independent of the initial density operator.
\par
In the remainder of this section we firstly define the indecomposability of CC-QSCs, and prove the independence of the information rate as in~\eqref{eq:def:qsc:IR:2} from the initial density operator.
Secondly, we generalize the methods in Algorithm~\ref{alg:SPA} for estimating such information rate efficiently.
\par
The definition of an indecomposable CC-QSC in this chapter is similar (but different) and closely related to that of an indecomposable (quantum) channel with memory in~\cite{bowen2005bounds}.
Namely, an indecomposable channel with memory equipped with separable input ensemble and local output measurement will always induce an indecomposable CC-QSC, but not necessarily vice versa.
Moreover, in~\cite{bowen2005bounds} the classical capacity of quantum channels with finite memory was considered, where the capacity is essentially the Holevo bound and where the latter was proven to be achievable~\cite{bowen2004quantum}.
However, in our work, we focus on the situation where the ensemble and the measurement are fixed.
\subsection{Indecomposable Quantum-State Channel}
\begin{definition}
A CC-QSC $\{\operator{N}^{y|x}\}_{x,y}$ is said to be \emph{indecomposable}\index{indecomposable} if for any initial density operators $\alpha_{\system{S}_0}$ and $\beta_{\system{S}_0}$, the following statement holds: for any $\epsilon>0$, there exists some positive integer $N$ s.t.
\begin{equation}
\norm{\alpha_{\system{S}_n}^{(\vx_1^n)} - \beta_{\system{S}_n}^{(\vx_1^n)}}_1 < \varepsilon \quad \forall n\geqslant N,\ \forall \vx_1^n\in\set{X}^n,
\end{equation}
\begin{align}
\alpha_{\system{S}_n}^{(\vx_1^n)}
&\defeq \sum_{\vy_1^n} \operator{N}^{y_n|x_n} \circ \cdots \circ \operator{N}^{y_1|x_1}(\alpha_{\system{S}_0}),\\
\beta_{\system{S}_n}^{(\vx_1^n)}
&\defeq \sum_{\vy_1^n} \operator{N}^{y_n|x_n} \circ \cdots \circ \operator{N}^{y_1|x_1}(\beta_{\system{S}_0}),
\end{align}
and where $\norm{A}_1$ is the trace distance for an operator $A$ on $\hilbert_\system{S}$, \ie, $\norm{A}_1\defeq\frac{1}{2}\tr\sqrt{A^\Herm A}$.
\end{definition}
\begin{theorem}\label{thm:quantum:indecom}\hspace{-5pt}\footnote{
A similar result regarding indecomposable/forgetful quantum channel with memory can be found in~\cite{kretschmann2005quantum} and~\cite{bowen2005bounds}.}
The information rate of an indecomposable CC-QSC with an {i.i.d.}~input process is independent of the initial density operator.
Namely, if $\{\operator{N}^{y|x}\}_{x,y}$ is indecomposable, then
\begin{equation}
\infoRate^{(n)}(Q,\{\operator{N}^{y|x}\}_{x,y},\alpha_{\system{S}_0}) - \infoRate^{(n)}(Q,\{\operator{N}^{y|x}\}_{x,y},\beta_{\system{S}_0}) \stackrel{n\to\infty}{\longrightarrow} 0
\end{equation}
for any initial density operators $\alpha_{\system{S}_0}$, $\beta_{\system{S}_0}
\in\DensOp(\hilbert_{\system{S}_0})$.
\end{theorem}
In the proof below, we follow a similar idea as in~\cite{gallager1968information} for indecomposable FSMCs, and as that in~\cite{bowen2005bounds} for indecomposable quantum channels with memory.
\begin{proof}
Let $\system{A}$ and $\system{B}$ be quantum systems described by Hilbert spaces $\hilbert_\system{A}$ and $\hilbert_\system{B}$, respectively, where $\{\bra{x}\}_{x\in\set{X}}$ and $\{\bra{y}\}_{y\in\set{Y}}$ are orthonormal bases of $\hilbert_\system{A}$ and $\hilbert_\system{B}$, respectively.
Let $\system{A}_1^n$ and $\system{B}_1^n$ be $n$ copies of $\system{A}$ and $\system{B}$, respectively.
Let $\rho_{\system{S}_0}$ be some initial density operator; and let the joint density operator on the system $\system{A}_1^n\system{B}_1^n$ be
\[
\rho_{\system{A}_1^n\system{B}_1^n} \defeq \sum_{\vx_1^n} Q(\vx_1^n) \cdot \braket{\vx_1^n} \tensor \sum_{\vy_1^n} \tr\left( \operator{N}^{\vy_1^n|\vx_1^n}(\rho_{\system{S}_0}) \right) \cdot \braket{\vy_1^n},
\]
where $\operator{N}^{\vy_1^n|\vx_1^n} \defeq \operator{N}^{y_n|x_n} \circ \cdots \circ \operator{N}^{y_1|x_1}$.
In this case, it is not hard to see that
\[
\mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)[\rho_{\system{S}_0}] = \qmutualInfo(\system{A}_1^n;\system{B}_1^n)[\rho_{\system{S}_0}].
\]
In fact, one can easily check that
\begin{align*}
\qEntropy(\system{A}_1^n) &= \entropy(\rv{X}_1^n),\\
\qEntropy(\system{B}_1^n) &= \entropy(\rv{Y}_1^n),\\
\qEntropy(\system{A}_1^n,\system{B}_1^n) &= \entropy(\rv{X}_1^n,\rv{Y}_1^n).
\end{align*}
In particular, $\qEntropy(\system{A}_1^n)$ is independent of the initial density operator $\rho_{\system{S}_0}$.
We also claim that, for each $\rho_{\system{S}_0} \in \DensOp(\hilbert_{\system{S}_0})$ and a positive integer $N < n$,
\begin{align}
\label{eq:indecomposable:AB}
\qmutualInfo(\system{A}_1^N\system{B}_1^N;\system{A}_{N+1}^n\system{B}_{N+1}^n)
&\leqslant 2 \qEntropy(\system{S}_N),\\
\label{eq:indecomposable:B}
\qmutualInfo(\system{B}_1^N;\system{B}_{N+1}^n)
&\leqslant 2 \qEntropy(\system{S}_N),
\end{align}
where the density operator for $\system{S}_N$ is defined as (depending on $\rho_{\system{S}_0}$)
\[
\rho_{\system{S}_N} \defeq \sum_{\vx_1^N} Q(\vx_1^N) \cdot \sum_{\vy_1^N} \operator{N}^{\vy_1^N|\vx_1^N}(\rho_{\system{S}_0}).
\]
\par
\textbf{Proof of~\eqref{eq:indecomposable:AB}: }
We define a class of CPTP maps $\{\Phi_a^b: \DensOp(\hilbert_{\system{S}_a}) \to \DensOp(\hilbert_{\system{A}_a^b \system{B}_a^b})\}_{a<b\in\mathbb{N}}$ as
\[
\Phi_a^b: \rho_{\system{S}_a} \mapsto \sum_{\vx_a^b} Q(\vx_a^b) \cdot \braket{\vx_a^b} \tensor \sum_{\vy_a^b} \tr\left( \operator{N}^{\vy_a^b|\vx_a^b}(\rho_{\system{S}_a}) \right) \cdot \braket{\vy_a^b}.
\]
Since the input process $Q$ is {i.i.d.}, we can rewrite $\rho_{\system{A}_1^n\system{B}_1^n}$ for each positive integer $N<n$, as
\[
\rho_{\system{A}_1^n\system{B}_1^n} = \left( I_{\system{A}_1^N\system{B}_1^N} \otimes \Phi_{N+1}^n \right) \left( \rho_{\system{A}_1^N\system{B}_1^N\system{S}_N} \right),
\]
\[
\rho_{\system{A}_1^N\system{B}_1^N\system{S}_N} \defeq \sum_{\vx_1^N} Q(\vx_1^N) \braket{\vx_1^N} \tensor \sum_{\vy_1^N} \operator{N}^{\vy_1^N|\vx_1^N}(\rho_{\system{S}_0}) \braket{\vy_1^N}.
\]
Hence, by data processing inequality for quantum mutual information (see \eg,~\cite[Theorem~11.9.4]{wilde2017quantum}), one must have
\[
\qmutualInfo(\system{A}_1^N\system{B}_1^N;\system{A}_{N+1}^n\system{B}_{N+1}^n) \leqslant \qmutualInfo(\system{A}_1^N\system{B}_1^N;\system{S}_N).
\]
Additionally, by subadditivity of joint entropy, we have
\begin{align*}
\qmutualInfo(\system{A}_1^N\system{B}_1^N;\system{S}_N)
&\defeq \qEntropy(\system{A}_1^N\system{B}_1^N) + \qEntropy(\system{S}_N) - \qEntropy(\system{A}_1^N\system{B}_1^N\system{S}_N)\\
&\leqslant \qEntropy(\system{A}_1^N\system{B}_1^N) + \qEntropy(\system{S}_N) - \abs{\qEntropy(\system{A}_1^N\system{B}_1^N) - \qEntropy(\system{S}_N)}\\
&\leqslant 2\qEntropy(\system{S}_N).
\end{align*}
Combining the above two inequalities, we have proven~\eqref{eq:indecomposable:AB}.
\par
\textbf{Proof of~\eqref{eq:indecomposable:B}: }
\eqref{eq:indecomposable:B} can be shown via the same approach above by considering another class of CPTP maps $\{\Psi_a^b: \DensOp(\hilbert_{\system{S}_a}) \to \DensOp(\hilbert_{\system{B}_a^b})\}_{a<b\in\mathbb{N}}$ as
\[
\Psi_a^b: \rho_{\system{S}_a} \mapsto \sum_{\vx_a^b} Q(\vx_a^b) \cdot \sum_{\vy_a^b} \tr\left( \operator{N}^{\vy_a^b|\vx_a^b}(\rho_{\system{S}_a})\right) \cdot \braket{\vy_a^b}.
\]
We omit the details.
\par
Now return to the main proof.
Given the initial density operators $\alpha_{\system{S}_0}$, and $\beta_{\system{S}_0}$, we define $\alpha_{\system{A}_1^n\system{B}_1^n}$, $\beta_{\system{A}_1^n\system{B}_1^n}$ and $\alpha_{\system{S}_N}$, $\beta_{\system{S}_N}$ in a similar fashion as we have defined $\rho_{\system{A}_1^n\system{B}_1^n}$ and $\rho_{\system{S}_N}$ based on $\rho_{\system{S}_0}$.
In this case, one obtains
\begin{align}
\nonumber
&\bigabs{\qEntropy(\alpha_{\system{A}_1^n\system{B}_1^n}) - \qEntropy(\beta_{\system{A}_1^n\system{B}_1^n})} - \bigabs{\qEntropy(\alpha_{\system{A}_{N\!+\!1}^n\system{B}_{N\!+\!1}^n}) - \qEntropy(\beta_{\system{A}_{N\!+\!1}^n\system{B}_{N\!+\!1}^n})}\\
\nonumber
\overset{\text{(a)}}{\leqslant}
&\biggabs{\qEntropy(\alpha_{\system{A}_1^N\system{B}_1^N}) \!-\! \qEntropy(\beta_{\system{A}_1^N\system{B}_1^N})} \!+\! \biggabs{\qmutualInfo(\system{A}_1^N\system{B}_1^N;\system{A}_{N\!+\!1}^n\system{B}_{N\!+\!1}^n)[\alpha_{\system{A}_1^n\system{B}_1^n}] \!-\! \qmutualInfo(\system{A}_1^N\system{B}_1^N;\system{A}_{N\!+\!1}^n\system{B}_{N\!+\!1}^n)[\beta_{\system{A}_1^n\system{B}_1^n}]}\\
\label{eq:tail:1}
\overset{\text{(b)}}{\leqslant}
& N\cdot\log(\dim{\hilbert_\system{AB}}) + 2\cdot\max\left\{\qEntropy(\alpha_{\system{S}_N}),\qEntropy(\beta_{\system{S}_N})\right\},
\end{align}
where we have used the triangle inequality in step (a), and Corollary~\ref{cor:qEntropy:basics} and~\eqref{eq:indecomposable:AB} in step (b).
Similarly, using~\eqref{eq:indecomposable:B}, one can prove
\begin{equation}\label{eq:tail:2}
\bigabs{\qEntropy(\alpha_{\system{B}_1^n}) - \qEntropy(\beta_{\system{B}_1^n})} - \bigabs{\qEntropy(\alpha_{\system{B}_{N+1}^n}) - \qEntropy(\beta_{\system{B}_{N+1}^n})}
\leqslant N \cdot \log(\dim{\hilbert_\system{B}}) + 2\cdot\max\left\{\qEntropy(\alpha_{\system{S}_N}),\qEntropy(\beta_{\system{S}_N})\right\}.
\end{equation}
By assumption, there exists some positive integer $d$ such that $\max\{\dim{\hilbert_\system{A}},\dim{\hilbert_\system{B}},\dim{\hilbert_\system{S}}\} \leqslant d$.
Thus, we have
\begin{fleqn}\begin{equation*}\phantom{=}
\frac{1}{n}\bigabs{\mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)[\alpha_{\system{S}_0}] - \mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)[\beta_{\system{S}_0}]}
\end{equation*}\end{fleqn}
\begin{fleqn}\begin{equation*}
=\frac{1}{n}\bigabs{\qmutualInfo(\system{A}_1^n;\system{B}_1^n)[\alpha_{\system{S}_0}] - \qmutualInfo(\system{A}_1^n;\system{B}_1^n)[\beta_{\system{S}_0}]}
\end{equation*}\end{fleqn}
\begin{fleqn}\begin{equation*}
=\frac{1}{n}\bigabs{\left( \qEntropy(\alpha_{\system{B}_1^n}) - \qEntropy(\alpha_{\system{A}_1^n\system{B}_1^n}) \right) - \left(\qEntropy(\beta_{\system{B}_1^n}) - \qEntropy(\beta_{\system{A}_1^n\system{B}_1^n})\right)}
\end{equation*}\end{fleqn}
\begin{fleqn}\begin{equation*}\overset{(\mathrm{c})}{\leqslant}
\frac{1}{n} \bigabs{\qEntropy(\alpha_{\system{B}_1^n}) - \qEntropy(\beta_{\system{B}_1^n})} + \frac{1}{n}\bigabs{\qEntropy(\alpha_{\system{A}_1^n\system{B}_1^n}) - \qEntropy(\beta_{\system{A}_1^n\system{B}_1^n})}
\end{equation*}\end{fleqn}
\begin{fleqn}\begin{equation*}\overset{(\mathrm{d})}{\leqslant}
\frac{3N\!+\!4}{n}\cdot\log{d} + \frac{1}{n}\bigabs{\qEntropy(\alpha_{\system{B}_{N+1}^n}) - \qEntropy(\beta_{\system{B}_{N+1}^n})} + \frac{1}{n}\bigabs{\qEntropy(\alpha_{\system{A}_{N+1}^n\system{B}_{N+1}^n}) - \qEntropy(\beta_{\system{A}_{N+1}^n\system{B}_{N+1}^n})}
\end{equation*}\end{fleqn}
\begin{fleqn}\begin{equation*}
=\frac{3N\!+\!4}{n}\cdot\log{d} + \frac{1}{n}\bigabs{\qEntropy(\Psi_{N+1}^n(\alpha_{\system{S}_N})) - \qEntropy(\Psi_{N+1}^n(\beta_{\system{S}_N}))} + \frac{1}{n}\bigabs{\qEntropy(\Phi_{N+1}^n(\alpha_{\system{S}_N})) - \qEntropy(\Phi_{N+1}^n(\beta_{\system{S}_N}))},
\end{equation*}\end{fleqn}
where we have used the triangle inequality in step (c), and~\cite[Theorem~11.8]{nielsen2011quantum}, \eqref{eq:tail:1}, \eqref{eq:tail:2} in step (d).
Using a loose variant of Fannes' inequality~\cite{fannes1973continuity}\footnote{
Namely, we used the inequality $\abs{\qEntropy(\rho)-\qEntropy(\sigma)} \leqslant \log{\dim}\cdot\norm{\rho-\sigma}_1+e^{-1}$.
Note that tighter variants of Fannes' inequality exist, but the above inequality is good enough to prove the desired result.}, we have
\begin{align*}
\bigabs{\qEntropy(\Psi_{N+1}^n(\alpha_{\system{S}_N})) \!-\! \qEntropy(\Psi_{N+1}^n(\beta_{\system{S}_N}))}
&\leqslant (n\!-\!N) \cdot \log{d}\ \cdot \norm{\Psi_{N+1}^n(\alpha_{\system{S}_N}) \!-\! \Psi_{N+1}^n(\beta_{\system{S}_N})}_1 + e^{-1},\\
\bigabs{\qEntropy(\Phi_{N+1}^n(\alpha_{\system{S}_N})) \!-\! \qEntropy(\Phi_{N+1}^n(\beta_{\system{S}_N}))}
&\leqslant 2\!\cdot\!(n\!-\!N)\cdot\log{d} \cdot \norm{\Phi_{N+1}^n(\alpha_{\system{S}_N}) \!-\! \Phi_{N+1}^n(\beta_{\system{S}_N})}_1 + e^{-1}.
\end{align*}
Moreover, by the contractivity of the trace distance, we have,
\begin{align*}
\norm{\Psi_{N+1}^n(\alpha_{\system{S}_N})-\Psi_{N+1}^n(\beta_{\system{S}_N})}_1
& \leqslant \norm{\alpha_{\system{S}_N}-\beta_{\system{S}_N}}_1,\\
\norm{\Phi_{N+1}^n(\alpha_{\system{S}_N})-\Phi_{N+1}^n(\beta_{\system{S}_N})}_1
& \leqslant \norm{\alpha_{\system{S}_N}-\beta_{\system{S}_N}}_1.
\end{align*}
This allows us to bound the difference of the information rate by
\[
\frac{1}{n}\bigabs{\mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)[\alpha_{\system{S}_0}] - \mutualInfo(\rv{X}_1^n;\rv{Y}_1^n)[\beta_{\system{S}_0}]} \leqslant \frac{3N\!+\!4}{n}\cdot\log{d} + \frac{3(n-N)}{n} \cdot \log{d} \cdot \norm{\alpha_{\system{S}_n} - \beta_{\system{S}_n}}_1 + \frac{2}{n\cdot e}.
\]
Finally, because the CC-QSC is indecomposable, for any $\varepsilon>0$, we can choose $N$ large enough such that
\[
\norm{\alpha_{\system{S}_N}-\beta_{\system{S}_N}}_1 < \frac{\varepsilon}{6\cdot\log{d}}
\]
and then choose an integer $M>N$ such that
\[
\frac{3N+4}{M}\cdot\log{d} + \frac{2}{M\cdot e}< \frac{\varepsilon}{2}.
\]
This will ensure that for any $n>M$, we have
\[
\frac{3N+4}{n} \cdot \log{d} + \frac{3(n-N)}{n} \cdot \log{d} \cdot \norm{\alpha_{\system{S}_n} - \beta_{\system{S}_n}}_1 + \frac{2}{n\cdot e} < \varepsilon,
\]
which concludes the proof.
\end{proof}
\subsection{Estimation of the Information Rate}
The development in this section is very similar to the development in Subsection~\ref{QCwM:sec:FSMC:IR}.
In particular, we follow the same approach as in~\eqref{eq:def:fsmc:ir:2}--\eqref{eq:FSMC:ir:estimate:1}.
This similarity stems from the similarity of the factor graphs in Figures~\ref{fig:FMSC:high:level:1} and~\ref{fig:qsc:bcjr}, and highlights one of the benefits of the factor-graph approach that we take to estimate information rate of quantum channels with memory.\par
We make the following assumptions.
\begin{itemize}
\item As already mentioned, the derivations in this chapter are for the case where the input process $\rv{X}_1^n=(\rv{X}_1,\ldots,\rv{X}_n)$ is an {i.i.d.}~process.
The results can be generalized to other stationary ergodic input processes that can be represented by a finite-state-machine source (FSMS).
Technically, this is done by defining a new state that combines the FSMS state and the channel state.
\item We assume that the corresponding quantum-state channel $\{\operator{N}^{y|x}\}_{x\in\set{X},y\in\set{Y}}$ is finite-dimensional and indecomposable.
We also assume it can be represented by some functions $\{W^{y|x}\}_{x,y}$ as defined in~\eqref{eq:def:channel:function:representation}.
\end{itemize}
The major difference compared with Section~\ref{QCwM:sec:FSMC:IR} is the conditional \pmf $\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}$, and thus the joint \pmf $\prob_{\rv{Y}_1^n,\rv{X}_1^n|\system{S}_0}$ as specified in~\eqref{eq:joint:4} and~\eqref{eq:joint:5}, respectively.
In this case, in order to compute $-\frac{1}{n}\log{\prob_{\rv{Y}_1^n}(\cvy_1^n)}$ and $-\frac{1}{n}\log{\prob_{\rv{X}_1^n\rv{Y}_1^n}(\cvx_1^n,\cvy_1^n)}$ using a similar method as in Section~\ref{QCwM:sec:FSMC:IR}, we consider the state metrics $\{\sigmaY_\ell\}_{\ell=1}^n$ and $\{\sigmaXY_\ell\}_{\ell=1}^n$ (which are operators on $\hilbert_{\system{S}_\ell}$ for each $\ell$) defined \wrt $\cvy_1^n$ and \wrt $\cvx_1^n$ and $\cvy_1^n$, respectively, as
\begin{align}
\sigmaY_\ell
&\defeq \sum_{\vx_1^\ell} Q^{(\ell)}(\vx_1^\ell) \cdot \operator{N}^{\cy_n|x_n} \circ \cdots \circ \operator{N}^{\cy_1|x_1}(\rho_{\system{S}_0}),\\
\sigmaXY_{\ell}
&\defeq \operator{N}^{\cy_n|\cx_n} \circ \cdots \circ \operator{N}^{\cy_1|\cx_1}(\rho_{\system{S}_0}).
\end{align}
In this case, we have $\prob_{\rv{Y}_1^n}(\cvy_1^n) = \tr(\sigmaY_n)$, and $\prob_{\rv{X}_1^n\rv{Y}_1^n}(\cvx_1^n,\cvy_1^n) = \tr(\sigmaXY_n)$.
Notice that $\{\sigmaY_\ell\}_{\ell}$ and $\{\sigmaXY_\ell\}_{\ell}$ can be computed iteratively as
\begin{align}
\label{eq:recursive:quantum:state:metric:Y:1}
&= \sum_{x_\ell} Q(x_\ell) \cdot [W^{y|x}] \cdot [\sigmaY_{\ell-1}],\\
\label{eq:def:quantum:channel:state:metric:XY:1}
&= [W^{y|x}] \cdot [\sigmaXY_{\ell-1}],
\end{align}
where we treat $[\sigma^\rv{Y}_{\ell}]$ and $[\sigma^\rv{XY}_{\ell}]$ as length-$d^2$ vectors indexed by $(s,\tilde{s})\in\set{S}^2$ in the above two equations.
(See~\eqref{eq:qsc:matrix:1} and~\eqref{eq:joint:4} for notations.)
Moreover, we can also introduce normalizing coefficients $\{\lambda_\ell^\rv{Y}\}_{\ell}$ and $\{\lambda_\ell^\rv{XY}\}_{\ell}$, similar to~\eqref{eq:recursive:state:metric:Y:2}, for the sake of numerical stability.
In the latter case, we have iterative updating rules
\begin{align}
\label{eq:recursive:quantum:state:metric:Y:2}
[\bsigmaY_\ell] &= \frac{1}{\lambdaY_\ell} \cdot \sum_{x_\ell} Q(x_\ell) \cdot
[W^{y|x}] \cdot [\bsigmaY_{\ell-1}],\\
\label{eq:def:quantum:channel:state:metric:XY:2}
[\bsigmaXY_\ell] &= \frac{1}{\lambdaXY_\ell} \cdot [W^{y|x}]
\cdot [\bsigmaXY_{\ell-1}],
\end{align}
where the scaling factors $\lambdaY_\ell>0$ and $\lambdaXY_\ell>0$ are chosen such that $\tr(\bsigmaY_\ell)=1$ and $\tr(\bsigmaXY_\ell)=1$, respectively.
In addition, one can verify that $\prob_{\rv{Y}_1^n}(\cvy_1^n) = \prod_{\ell=1}^n \lambdaY_\ell$, and $\prob_{\rv{X}_1^n\rv{Y}_1^n}(\cvx_1^n,\cvy_1^n) = \prod_{\ell=1}^n \lambdaXY_\ell$.
\par
The above discussion is summarized as Algorithm~\ref{alg:SPA:2}.
The computations corresponding to Line 3, 5--9 and 12--16 are visualized in Figures~\ref{fig:QFSM:channel:simulation:Y},~\ref{fig:QFSM:estimate:hY}, and~\ref{fig:QFSM:estimate:hXY} in the Appendix~\ref{app:figures}, respectively.
\begin{algorithm}
\caption{Estimating the Information Rate of a CC-QSC}
\begin{algorithmic}[1]
\Require{An indecomposable CC-QSC $\{\operator{N}^{y|x}\}_{x\in\set{X},y\in\set{Y}}$, which can be represented by functions $\{W^{y|x}\}_{x,y}$, a input~distribution~$Q$, a positive~integer~$n$ large enough.}
\Ensure{$\infoRate^{(n)}(Q,\{\operator{N}^{y|x}\}_{x,y}) \approx \entropy(\rv{X})+\hat\entropicRate(\rv{Y}) - \hat\entropicRate(\rv{X},\rv{Y})$.}
\State Initialize the memory density operator $\rho_{\system{S}_0} \gets \braket{0_\system{S}}$
\State Generate an input sequence $\cvx_1^n \sim Q^{\tensor n}$
\State Generate a corresponding output sequence $\cvy_1^n$
\State $\bsigmaY_0 \gets \rho_{\system{S}_0}$
\ForEach{$\ell=1,\ldots,n$}
\State $[\sigmaY_\ell] \gets \sum_{x_\ell} Q(x_\ell)\cdot [W^{\cy_{\ell}|x}] \cdot [\bsigmaY_{\ell-1}]$
\State $\lambdaY_\ell \gets \tr(\sigmaY_\ell)$
\State $\bsigmaY_{\ell} \gets \sigmaY_\ell/\lambdaY_\ell$
\EndFor
\State $\hat\entropicRate(\rv{Y}) \gets -\frac{1}{n} \sum_{\ell=1}^n \log(\lambdaY_{\ell})$
\State $\bsigmaXY_0 \gets \rho_{\system{S}_0}$
\ForEach{$\ell=1,\ldots,n$}
\State $[\sigmaXY_\ell] \gets [W^{\cy_\ell|\cx_\ell}] \cdot [\bsigmaXY_{\ell-1}]$
\State $\lambdaXY_\ell \gets \tr(\sigmaXY_\ell)$
\State $\bsigmaXY_{\ell} \gets \sigmaXY_\ell/\lambdaXY_\ell$
\EndFor
\State $\hat\entropicRate(\rv{X},\rv{Y}) \gets -\frac{1}{n} \sum_{\ell=1}^n \log(\lambdaXY_{\ell})$
\State $\entropy(\rv{X}) \gets -\sum_{x} Q(x) \log{Q(x)}$
\State Estimate $\infoRate^{(n)}(Q,\{\operator{N}^{y|x}\}_{x,y})$ as $\entropy(\rv{X}) +\hat\entropicRate(\rv{Y}) - \hat\entropicRate(\rv{X},\rv{Y})$.
\end{algorithmic}
\label{alg:SPA:2}
\end{algorithm}
\section[IRUB/IRLB and their Optimization]{Information rate upper/lower bounds and their Optimization} \label{QCwM:sec:5:UBLB}
In this section, we consider auxiliary channels and their induced upper and lower bounds on the information rate.
As mentioned in the beginning of this chapter auxiliary channels are often introduced as a low-complexity approximation of the original channel, which are useful in mismatch decoding.
The techniques developed in this section only require the channel input/output data, but not the channel model itself.
This is particularly useful when the channel is only made physically, but not mathematically, available.
In this case, the task of minimizing the difference between the upper and lower bound is equivalent to finding the channel model (within a specified class of channel models) best fitting the \emph{empirical} channel law.
Similarly, minimizing the upper bound corresponds to finding the channel model best fitting the \emph{empirical} channel output distribution, and maximizing the lower bound corresponds to finding the channel model best fitting the \emph{empirical} reverse channel law.
Motivated by the above scenarios, we particularly consider the auxiliary channels chosen from the domain of all CC-QSCs with the same input and output alphabet as the original channel, and acting on a memory system of a certain dimension (which can be different from the memory dimension of the original channel).
Throughout this section, we assume the original channel as described in Section~\ref{QCwM:sec:3:QCM} is indecomposable, and that all the involved Hilbert spaces are of finite dimension, and that the alphabets $\set{X}$ and $\set{Y}$ are finite.
\par
Suppose we have some auxiliary CC-QSC $\{\hat{\operator{N}}^{y|x}\}_{x,y}$, describable by some functions $\{\hat{W}^{y|x}\}_{x,y}$ as in~\eqref{eq:def:channel:function:representation}.
Let $\hat{P}_{\rv{Y}_1^n|\rv{X}_1^n,\hat{\system{S}}_0}$ denote its joint channel law, similar to~\eqref{eq:joint:2},~\eqref{eq:joint:3}, or~\eqref{eq:joint:4}.
\begin{equation}
\hat{P}_{\rv{Y}_1^n|\rv{X}_1^n,\hat{\system{S}}_0}(\vy_1^n|\vx_1^n;\hat{\rho}_{\system{S}_0}) \defeq \tr\left( [\What^{y_{n}|x_{n}}] \cdots [\What^{y_{1}|x_{1}}] \cdot [\hat{\rho}_{\system{S}_0}] \right).
\end{equation}
We follow a similar approach as in~\cite{arnold2006simulation, sadeghi2009optimization}, and define the quantities
\begin{align}
\label{eq:IRUB}\begin{aligned}
\IRUB_{W}^{(n)}(\hat{W}) \defeq &
\frac{1}{n} \sum_{\vx_1^n,\vy_1^n} Q^{(n)}(\vx_1^n) \cdot \prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})\\
&\cdot \log\frac{\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})}{\sum_{\cvx_1^n} Q^{(n)}(\cvx_1^n) \hat{P}_{\rv{Y}_1^n|\rv{X}_1^n,\hat{\system{S}}_0}(\vy_1^n|\cvx_1^n;\rho_{\hat{\system{S}}_0})},
\end{aligned}\\
\label{eq:IRLB}\begin{aligned}
\IRLB_{W}^{(n)}(\hat{W}) \defeq &
\frac{1}{n} \sum_{\vx_1^n,\vy_1^n} Q^{(n)}(\vx_1^n) \cdot \prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})\\
&\cdot \log\frac{\hat \prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})}{\sum_{\cvx_1^n} Q^{(n)}(\cvx_1^n) \hat{P}_{\rv{Y}_1^n|\rv{X}_1^n,\hat{\system{S}}_0}(\vy_1^n|\cvx_1^n;\rho_{\hat{\system{S}}_0})},
\end{aligned}
\end{align}
where $\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}$ is defined in~\eqref{eq:joint:2},~\eqref{eq:joint:3} or~\eqref{eq:joint:4}.
By following similar arguments like those in~\eqref{eq:IRUB:minus:IR} and~\eqref{eq:IR:minus:IRLB}, one can verify that
\begin{equation}
\IRLB_{W}^{(n)}(\hat{W}) \leqslant \infoRate^{(n)}_{W} \leqslant \IRUB_{W}^{(n)}(\hat{W}),
\end{equation}
where the first inequality holds with equality if and only if $\hat{P}_{\rv{Y}_1^n|\rv{X}_1^n,\hat{\system{S}}_0}(\vy_1^n|\vx_1^n;\rho_{\hat{\system{S}}_0})$ and $\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})$ coincide for all $\vx_1^n$ and $\vy_1^n$ with positive support of $\prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}$ and where the second inequalities holds with equality if and only if $\hat{P}_{\rv{Y}_1^n|\hat{\system{S}}_0}(\vy_1^n|\rho_{\hat{\system{S}}_0})$ and $\prob_{\rv{Y}_1^n|\system{S}_0}(\vy_1^n|\rho_{\system{S}_0})$ coincide for all $\vy_1^n$ with positive support of $\prob_{\rv{Y}_1^n|\system{S}_0}$.
Another quantity of interest is the \emph{difference function} defined as
\begin{equation} \label{eq:DIFF}
\Delta_{W}^{(n)}(\hat{W}) \defeq \IRUB_{W}^{(n)}(\hat{W}) - \IRLB_{W}^{(n)}(\hat{W}).
\end{equation}
Explicit expressions of~\eqref{eq:IRUB},~\eqref{eq:IRLB}, and~\eqref{eq:DIFF} are given by
\begin{align}
\label{eq:IRUB:explicit}
\IRUB^{(n)}_W(\hat{W}) &= \frac{1}{n}\expectationwrt{
\log\frac{\tr\left([W^{\rv{Y}_n|\rv{X}_n}] \cdots [W^{\rv{Y}_1|\rv{X}_1}] \cdot [\rho_{\system{S}_0}]\right)}{\sum_{\vx_1^n} Q^{(n)}(\vx_1^n) \cdot \tr\left([\hat{W}^{\rv{Y}_n|x_n}] \cdots [\hat{W}^{\rv{Y}_1|x_1}] \cdot [\rho_{\hat{\system{S}}_0}]\right)}}{\rv{X}_1^n\rv{Y}_1^n},\\
\label{eq:IRLB:explicit}
\IRLB_{W}^{(n)}(\hat{W}) &= \frac{1}{n}\expectationwrt{
\log\frac{\tr\left([\hat{W}^{\rv{Y}_n|\rv{X}_n}] \cdots [\hat{W}^{\rv{Y}_1|\rv{X}_1}]\cdot [\rho_{\hat{\system{S}}_0}]\right)}{\sum_{\vx_1^n} Q^{(n)}(\vx_1^n) \cdot \tr\left([\hat{W}^{\rv{Y}_n|x_n}]\cdots[\hat{W}^{\rv{Y}_1|x_1}] \cdot [\rho_{\hat{\system{S}}_0}]\right)}}{\rv{X}_1^n\rv{Y}_1^n},\\
\label{eq:DIFF:explicit}
\Delta_{W}^{(n)}(\hat{W}) &= \frac{1}{n}\expectationwrt{
\log\frac{\tr\left([W^{\rv{Y}_n|\rv{X}_n}] \cdots [W^{\rv{Y}_1|\rv{X}_1}] \cdot [\rho_{\system{S}_0}]\right)}{\tr\left([\hat{W}^{\rv{Y}_n|\rv{X}_n}] \cdots [\hat{W}^{\rv{Y}_1|\rv{X}_1}] \cdot [\rho_{\hat{\system{S}}_0}]\right)}}{\rv{X}_1^n\rv{Y}_1^n},
\end{align}
respectively, where $\rv{X}_1^n$ and $\rv{Y}_1^n$ are random variables distributed according to the joint distribution $Q^{(n)}(\vx_1^n)\cdot \prob_{\rv{Y}_1^n|\rv{X}_1^n;\system{S}_0}(\vy_1^n|\vx_1^n;\rho_{\system{S}_0})$ and where $\left\langle\cdot\right\rangle$ stands for the expectation function.
\par
In the remainder of this section, we propose an algorithm based on the gradient-descent method and the techniques described in Section~\ref{QCwM:subsec:NFGs} and~\ref{QCwM:sec:4:IR} for optimizing the quantities in~\eqref{eq:IRUB},~\eqref{eq:IRLB}, and~\eqref{eq:DIFF}.
In particular, we consider $\{\hat{W}^{y|x}\}_{x,y}$ to be an \emph{interior} point in the domain of CC-QSCs, namely
\begin{itemize}
\item The Choi--Jamio{\l}kowski matrices $\llbracket \hat{W}^{y|x}\rrbracket$, defined similarly as~\eqref{eq:qsc:matrix:2}, are strictly positive definite for each $x$ and $y$,
\item Eq.~\eqref{eq:operator:sum:representation:condition:2} holds by replacing $W^{y|x}$ with $\hat{W}^{y|x}$, namely $\sum_{y\in\set{Y}} \sum_{s',\tilde{s}':\: s'=\tilde{s}'} \llbracket \hat{W}^{y|x}\rrbracket_{(s',s),(\tilde{s}',\tilde{s})} = \delta_{s,\tilde{s}}$ for all $x\in\set{X}$.
\end{itemize}
For any set of functions $\{H^{y|x}:\set{S}^4\to \mathbb{C}\}_{x,y}$ such that $\llbracket H^{y|x}\rrbracket$ (again, defined similarly as~\eqref{eq:qsc:matrix:2}) is Hermitian for each $x$ and $y$ and such that
\begin{equation}\label{eq:tangent:linear:constraint}
\sum_{y\in\set{Y}}\sum_{s',\tilde{s}':\: s'=\tilde{s}'} \llbracket H^{y|x}\rrbracket_{(s',s),(\tilde{s}',\tilde{s})} = 0
\quad \forall x\in\set{X},
\end{equation}
the functions $\{\hat{W}^{y|x}+t\cdot H^{y|x}\}_{x,y}$ describe a valid CC-QSC, for all $t$ in some neighborhood of $0$.
In this case, the directional derivatives of the functions $\IRLB_{W}^{(n)}$, $\IRUB_{W}^{(n)}$, and $\Delta_{W}^{(n)}$ at $\{\hat{W}^{y|x}\}_{x,y}$ \emph{along} $\{H^{y|x}\}_{x,y}$ is well defined, and can be expressed as~\eqref{eq:IRUB:dev:2},~\eqref{eq:IRLB:dev:2}, and~\eqref{eq:DIFF:dev:2} on next page,
\begin{sidewaysfigure}
\begin{align}
\left.\frac{\D}{\D t}\right|_{t=0} \IRUB_{W}^{(n)}(\hat{W}+tH)
&\propto -\frac{1}{n} \expectationwrt{
\sum_{k=1}^n \sum_{\vx_1^n} Q^{(n)}(\vx_1^n) \cdot \tr\left( [\hat{W}^{\rv{Y}_n|x_n}] \cdots [\hat{W}^{\rv{Y}_{k\!+\!1}|x_{k\!+\!1}}] [H^{\rv{Y}_k|x_k}] [\hat{W}^{\rv{Y}_{k\!-\!1}|x_{k\!-\!1}}] \cdots [\hat{W}^{\rv{Y}_1|x_1}] \cdot [\rho_{\hat{\system{S}}_0}] \right)}{\rv{Y}_1^n} \nonumber\\
\label{eq:IRUB:dev:2}
&= -\frac{1}{n} \sum_{\vx_1^n,\vy_1^n} \prob_{\rv{X}_1^n,\rv{Y}_1^n|\system{S}_0}(\vx_1^n,\vy_1^n|\rho_{\system{S}_0}) \cdot \sum_{k} \sum_{s',s,\tilde{s}',\tilde{s}} \rvec{\varrho}_{\hat{\system{S}}_{k\!-\!1}}^{(\vy_1^{k\!-\!1})}(s,\tilde{s}) \cdot H^{y_k|x_k}(s',s,\tilde{s}',\tilde{s}) \cdot \lvec{\varrho}_{\hat{\system{S}}_k}^{(\vy_{k\!+\!1}^n)}(s',\tilde{s}')\\
\left.\frac{\D}{\D t}\right|_{t=0} \IRLB_{W}^{(n)}(\hat{W}+tH)
\sum_{k=1}^n \sum_{\vx_1^n} Q^{(n)}(\vx_1^n) \cdot \tr\left([\hat{W}^{\rv{Y}_n|x_n}] \cdots [\hat{W}^{\rv{Y}_{k\!+\!1}|x_{k\!+\!1}}] [H^{\rv{Y}_k|x_k}] [\hat{W}^{\rv{Y}_{k\!-\!1}|x_{k\!-\!1}}] \cdots [\hat{W}^{\rv{Y}_1|x_1}] \cdot [\rho_{\hat{\system{S}}_0}] \right)}{\rv{Y}_1^n}\\
\sum_{k=1}^n \tr\left([\hat{W}^{\rv{Y}_n|\rv{X}_n}] \cdots [\hat{W}^{\rv{Y}_{k\!+\!1}|\rv{X}_{k\!+\!1}}] [H^{\rv{Y}_k|\rv{X}_k}] [\hat{W}^{\rv{Y}_{k\!-\!1}|\rv{X}_{k\!-\!1}}] \cdots [\hat{W}^{\rv{Y}_1|\rv{X}_1}] \cdot [\rho_{\hat{\system{S}}_0}] \right)}{\rv{X}_1^n\rv{Y}_1^n}
\end{aligned}\nonumber\\
\label{eq:IRLB:dev:2}
&-\frac{1}{n} \sum_{\vx_1^n,\vy_1^n} \prob_{\rv{X}_1^n,\rv{Y}_1^n|\system{S}_0}(\vx_1^n,\vy_1^n|\rho_{\system{S}_0}) \cdot \sum_{k}\sum_{s',s,\tilde{s}',\tilde{s}} \rvec{\varrho}_{\hat{\system{S}}_{k\!-\!1}}^{\,(\vy_1^{k\!-\!1})}(s,\tilde{s}) \cdot H^{y_k|x_k}(s',s,\tilde{s}',\tilde{s}) \cdot \lvec{\varrho}_{\hat{\system{S}}_k}^{\,(\vy_{k\!+\!1}^n)}(s',\tilde{s}')\\
&+\frac{1}{n} \sum_{\vx_1^n,\vy_1^n} \prob_{\rv{X}_1^n,\rv{Y}_1^n|\system{S}_0}(\vx_1^n,\vy_1^n|\rho_{\system{S}_0}) \cdot \sum_{k}\sum_{s',s,\tilde{s}',\tilde{s}} \rvec{\varrho}_{\hat{\system{S}}_{k\!-\!1}}^{(\vx_1^{k\!-\!1},\vy_1^{k\!-\!1})}(s,\tilde{s}) \cdot H^{y_k|x_k}(s',s,\tilde{s}',\tilde{s})\cdot \lvec{\varrho}_{\hat{\system{S}}_k}^{(\vx_{k\!+\!1}^n,\vy_{k\!+\!1}^n)}(s',\tilde{s}')
\end{aligned}\\
\left.\frac{\D}{\D t}\right|_{t=0} \Delta_{W}^{(n)}(\hat{W}+tH)
\sum_{k=1}^n \tr\left( [\hat{W}^{\rv{Y}_n|\rv{X}_n}] \cdots [\hat{W}^{\rv{Y}_{k\!+\!1}|\rv{X}_{k\!+\!1}}] [H^{\rv{Y}_k|\rv{X}_k}] [\hat{W}^{\rv{Y}_{k\!-\!1}|\rv{X}_{k\!-\!1}}] \cdots [\hat{W}^{\rv{Y}_1|\rv{X}_1}] \cdot [\rho_{\hat{\system{S}}_0}]\right)}{\rv{X}_1^n\rv{Y}_1^n}\nonumber\\
\label{eq:DIFF:dev:2}
&=-\frac{1}{n}\sum_{\vx_1^n,\vy_1^n} \prob_{\rv{X}_1^n,\rv{Y}_1^n|\system{S}_0}(\vx_1^n,\vy_1^n|\rho_{\system{S}_0}) \cdot \sum_{k}\sum_{s',s,\tilde{s}',\tilde{s}} \rvec{\varrho}_{\hat{\system{S}}_{k\!-\!1}}^{\,(\vx_1^{k\!-\!1},\vy_1^{k\!-\!1})}(s,\tilde{s}) \cdot H^{y_k|x_k}(s',s,\tilde{s}',\tilde{s}) \cdot \lvec{\varrho}_{\hat{\system{S}}_k}^{\,(\vx_{k\!+\!1}^n,\vy_{k\!+\!1}^n)}(s',\tilde{s}')
\end{align}
\end{sidewaysfigure}
where we define the messages $\{\rvec{\varrho}_{\system{S}_\ell}^{(\cvy_1^\ell)}\}_\ell$, $\{\lvec{\varrho}_{\system{S}_\ell}^{(\cvy_{\ell\!+\!1}^n)}\}_\ell$, $\{\rvec{\varrho}_{\system{S}_\ell}^{(\cvx_1^\ell,\cvy_1^\ell)}\}_\ell$, and $\{\lvec{\varrho}_{\system{S}_\ell}^{(\cvx_{\ell\!+\!1}^n,\cvy_{\ell\!+\!1}^n)}\}_\ell$ in a recursive manner as
\begin{align}
[\rvec{\varrho}_{\system{S}_\ell}^{(\cvy_1^\ell)}] &\defeq
\sum_{\vx_1^\ell} Q(\vx_1^\ell) \cdot [\hat{W}^{\cy_\ell|x_\ell}] \cdots [\hat{W}^{\cy_1|x_1}] \cdot [\rho_{\system{S}_0}],\\
[\lvec{\varrho}_{\system{S}_\ell}^{(\cvy_{\ell\!+\!1}^n)}] &\defeq
\sum_{\vx_{\ell\!+\!1}^n} Q(\vx_{\ell+1}^n) \cdot [I_{\system{S}_n}] \cdot [\hat{W}^{\cy_n|x_n}] \cdots [\hat{W}^{\cy_{\ell\!+\!1}|x_{\ell\!+\!1}}],\\
[\rvec{\varrho}_{\system{S}_\ell}^{(\cvx_1^\ell,\cvy_1^\ell)}] &\defeq
[\hat{W}^{\cy_\ell|\cx_\ell}] \cdots [\hat{W}^{\cy_1|\cx_1}] \cdot [\rho_{\system{S}_0}],\\
[\lvec{\varrho}_{\system{S}_\ell}^{(\cvx_{\ell\!+\!1}^n,\cvy_{\ell\!+\!1}^n)}] &\defeq
[I_{\system{S}_n}] \cdot [\hat{W}^{\cy_n|\cx_n}] \cdots [\hat{W}^{\cy_{\ell\!+\!1}|\cx_{\ell\!+\!1}}].
\end{align}
Recall that, in above equations, $[I_{\system{S}_n}]$ is a row vector, and $[\rho_{\system{S}_0}]$ is a column vector.
\par