\section{Introduction} \label{sec_introduction}

The Shannon capacity of a channel characterizes the ultimate transmission efficiency limit of reliable communication over a channel \cite{Shannon1948}. It has played a fundamental role in the research of communication theory and directed the development of communication systems for decades.

In many scenarios of practical interest, however, the perfect knowledge about a channel may not be available or may not be fully utilized for implementing the transceivers. Important examples include channels with uncertainty (like fading in wireless communication systems) \cite{2002Fading}, with non-ideal transceiver hardware \cite{zhang2011general}, or with constrained receiver structure \cite{salz1995com}. A commonly adopted practice of the receiver under such circumstances is to use a prescribed decoding metric, which may not be matched to the actual channel transition probability, and the mismatch capacity has been introduced to characterize the highest information rate under a prescribed decoding metric; see, e.g., \cite{1998Reliable} \cite{2020Information} and references therein.

Unfortunately, the mismatch capacity is still an open problem to date \cite{1995Channel}. So instead of it, researchers have been focusing on deriving its achievable lower bounds. The generalized mutual information (GMI) is a relatively simple lower bound \cite{1993Information}, which has found extensive applications in various setups; see, e.g., \cite{2002Fading} \cite{zhang2011general} \cite{2012Nearest}. However, the LM rate  \cite{Merhav1994On} is a tighter lower bound, by replacing the independent and identically distributed (i.i.d.) codebook ensemble for the GMI by constant-composition codebook ensemble.\footnote{The LM rate is \textcolor{black}{still not} the tightest lower bound of the mismatch capacity. There are several ways of improving the GMI and the LM rate by considering more structured codebook ensembles \cite{2020Information}. These are beyond the scope of this paper.
}

From the perspective of optimization, the primal forms of the GMI and the LM rate are highly similar, except that for the LM rate there is an additional constraint on the marginal distribution of the sought-for maximizing joint probability distribution over the input and output alphabets; see, e.g., \cite[Thm. 1]{2000Mismatched}. Consequently, the dual form of the GMI is easy to compute, since it can be written as a maximization over a real number \cite[Eqn. (12)]{2000Mismatched}, which is readily solved by a one-dimensional line search. On the other hand, the dual form of the LM rate further involves a maximization over a function of the channel input \cite[Eqn. (11)]{2000Mismatched}. Consequently the computation of the LM rate is much more challenging, and there is rare work on the numerical computation of the LM rate. Interior point methods (e.g., \cite{sdpt3,2010cvx}) can be applied for computing the LM rate, but it is typically memory intensive and computationally expensive. It is therefore desirable to develop effective computation models and efficient numerical algorithms for the LM rate.




In this paper, we formulate the LM rate computation as an optimal transport (OT) problem with an extra constraint, motivated by the similarity between these two optimization problems. 
Our contribution consists of three parts.
First, we propose an OT model with an extra constraint for the LM rate. 
Second, \textcolor{black}{we show that} this model can be solved directly using the well-known Sinkhorn algorithm \cite{2013sinkhorn}. 
Since our OT model is strongly convex, no additional regularization is required. 
This guarantees that we only need \textcolor{black}{several hundred of} Sinkhorn iterations to converge rather than tens of thousands of iterations for classical OT problems.
Last, we show that the extra constraint \textcolor{black}{in the OT problem formulation} can be efficiently handled by finding the root of a one-dimensional monotone function, and the feasibility of the extra constraint is guaranteed at each iteration. Numerical experiments show that for the case of QPSK and 16-QAM, our proposed algorithm is efficient and accurate. We also point out that our method is highly scalable and can be easily generalized to large-scale cases, such as 256-QAM.



The remaining \textcolor{black}{part} of this paper is organized as follows. 
In Section \ref{sec_formulation_2}, after briefly reviewing the basic definition of the LM rate, we write it in the form of OT model with an extra constraint.
Next, we present the numerical methods for this problem, including the Sinkhorn algorithm and the treatment of the extra constraint in Section \ref{sec_Sinkhorn_3}. 
In Section \ref{sec_num_4}, the simulation results demonstrate the advantages of our approach. We finally conclude the paper in Section \ref{sec_conclu_5}.



\section{Problem Formulation} \label{sec_formulation_2}
\subsection{The LM Rate}
We consider a discrete memoryless communication channel model with transition law $W(y|x)$ over the (finite, discrete) channel input alphabet $\mathcal{X}= \{x_1,\cdots,x_M\}\subset\mbbR^2$ and the channel output alphabet $\mathcal{Y} \subset\mbbR^2$.
Thus, given a probability measure $P_{X}$ on $\mathcal{X}$, we are able to define the joint probability distribution $P_{XY}$ on $\mathcal{X}\times\mathcal{Y}$ and the output distribution $P_{Y}$ on $\mathcal{Y}$ by the following:
\begin{align*}
    & P_{XY}(x_i,B)= W(B|x_i)P_{X}(x_i), ~ \forall B\subset\mathcal{Y}, \\
    & P_{Y}(B) = \sum_{i=1}^{M} W(B|x_{i})P_{X}(x_i), ~ \forall B\subset\mathcal{Y}.
\end{align*}
This discrete memoryless channel is in fact a mapping from the input alphabet $\mathcal{X}$ to the output alphabet $\mathcal{Y}$ with the channel transition law $W(y|x)$.

\textcolor{black}{For transmission of rate $R$, a block length-$n$ codebook $\mathcal{C}$ consists of $2^{nR}$ vectors $\mathbf{x}(m) = (x^{(1)}(m),\cdots,x^{(n)}(m))\in\mathcal{X}^{n}$. The encoder maps the message $m$ uniformly randomly selected from the set $\mathcal{M} = \{1,\cdots,2^{nR}\}$ to the corresponding codeword $\mathbf{x}(m)$. After receiving the output sequence $\mathbf{y} = (y^{(1)},\cdots,y^{(n)})$, the decoder forms the estimation on the message $m$ following the decoding rule
\begin{equation*}
    \hat{m} = \argmin_{j\in\mathcal{M}}\sum_{k=1}^{n}d(x^{(k)}(j),y^{(k)}),
\end{equation*}
where $d : \mathcal{X}\times\mathcal{Y}\to\mbbR$ is  a non-negative function called the \textit{decoding metric}.}


\textcolor{black}{With the notations defined above, the LM rate, as one of the achievable lower bound of the mismatch capacity, is defined as:}
\begin{subequations}\label{IM}
	\begin{align}
	I_{LM}(P_{X}) & := \min_{\textcolor{black}{\gamma\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})}} D(\gamma\|P_{X}P_{Y}) \label{IM_obj} \\
	\textrm{s.t.}\quad \int_{y\in\mathcal{Y}}\gamma(x_i,y) \mathrm{d} y &= P_X(x_i),\ \forall x_i\in\mathcal{X}, \label{IM_Px} \\
	\sum_{x_i\in\mathcal{X}}\gamma(x_i,y) &= P_Y(y),\ \forall y\in\mathcal{Y}, \label{IM_PY} \\
	\int_{y\in\mathcal{Y}}\sum_{x_i\in\mathcal{X}}\gamma d(x_i,y) \mathrm{d} y &
	\le \int_{y\in\mathcal{Y}}\sum_{x_i\in\mathcal{X}}P_{XY}d(x_i,y) \mathrm{d} y. \label{IM_ineq}
	\end{align}
\end{subequations}
Here $D(\cdot\|\cdot)$ denotes the relative entropy function and
\textcolor{black}{$\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ denotes the set of all joint probability distributions on $\mathcal{X}\times\mathcal{Y}$.}






\subsection{The OT Problem}
In \eqref{IM}, we can see \textcolor{black}{that} the LM rate is similar to the OT problem \cite{peyre2019computational} in terms of optimization objective and constraints. Naturally, we would like to derive an equivalent OT-type problem to \eqref{IM}, which creates a prerequisite for computing the LM rate by using the Sinkhorn algorithm \cite{2013sinkhorn}.

Consider a discrete probability distribution $\mu_i=P_X(x_i)$ on $\mathcal{X}$. It satisfies the following normalization and \textcolor{black}{average} power constraints:
\begin{equation}\label{pc}
    \sum_{i=1}^{M}\mu_{i} = 1,\quad\sum_{i=1}^{M}\mu_{i}\|x_{i}\|_{2}^{2} = 1.
\end{equation}
\textcolor{black}{For given channel transition law $W(y|x)$ and decoding metric $d(x,y)$, the right-hand side of \eqref{IM_ineq} is a constant:}
\begin{equation}\label{const_T}
\begin{aligned}
    T&=\int_{\mathcal{Y}}\sum_{i=1}^{M} P_{XY}(x_i,y)d(x_i,y) \mathrm{d} y \\
    &= \int_{\mathcal{Y}}\sum_{i=1}^{M} \mu_{i}W(y|x_i)d(x_i,y) \mathrm{d} y.
\end{aligned}
\end{equation}
Moreover, it is easy to show that
\begin{align*}
    & D(\gamma\|P_{X}P_{Y})=\int_{\mathcal{Y}}\sum_{i=1}^{M}
        \gamma(x_{i},y)\log\gamma(x_{i},y)\mathrm{d} y-C, \\
    & C=\sum_{i=1}^MP_{X}(x_i)\log P_{X}(x_i)
        +\int_{\mathcal{Y}}P_{Y}(y)\log P_{Y}(y) \mathrm{d} y.
\end{align*}
By neglecting the constant $C$, minimizing the objective function $D(\gamma\|P_{X}P_{Y})$ is equivalent to maximizing the entropy
\begin{equation*}
    -\int_{\mathcal{Y}}\sum_{i=1}^{M}\gamma(x_{i},y)
        \log\gamma(x_{i},y)\mathrm{d} y.
\end{equation*}
Thus, we have obtained the following optimal transport problem \eqref{LM_obj}-\eqref{LM_y} with an extra constraint \eqref{LM_ineq}:
\begin{small}
\begin{subequations} \label{LM}
	\begin{align} 
	\min_{\textcolor{black}{\gamma\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})}}&\int_{\mathcal{Y}}\sum_{i=1}^{M}\gamma(x_{i},y)\log\gamma(x_{i},y)\mathrm{d} y \label{LM_obj} \\
	\textbf{s.t.}~&\int_{\mathcal{Y}}\gamma(x_{i},y) \mathrm{d} y = \mu_i,\quad i = 1,\cdots,M, \label{LM_x} \\
	&\sum_{i=1}^{M}\gamma(x_{i},y) = P_{Y}(y),\quad \forall y\in \mathcal{Y}, \label{LM_y} \\
	&\int_{\mathcal{Y}}\sum_{i=1}^{M}\gamma(x_{i},y) d(x_i,y)\mathrm{d} y \leq T. \label{LM_ineq}
	\end{align}
\end{subequations}
\end{small}
\begin{remark}
   
    \textcolor{black}{In the classical OT theory, the cost function is only related to distance and should be independent of variables. Thus, \eqref{LM_obj}-\eqref{LM_ineq} is not a standard OT form, since the ``cost function'' $\log\gamma(x_i,y)$ in \eqref{LM_obj} is related to the variables $\gamma(x_i,y)$. However, the optimization objective function \eqref{LM_obj} is obtained by complicated and lengthy inequality scaling of the original problem (2.29) in \cite{2020Information}, which is in standard OT form. For detailed derivation, the interested readers are referred to Chapter 2 in \cite{csiszar2011information}. On the other hand, from the numerical point of view, the optimization problem \eqref{LM_obj}-\eqref{LM_ineq} can be efficiently solved by the Sinkhorn algorithm, which is well-known for its high efficiency in solving OT problems. In these senses, we might as well call \eqref{LM_obj}-\eqref{LM_ineq} an OT problem.}
    
   
   
   
\end{remark}




\section{THE Sinkhorn ALGORITHM} \label{sec_Sinkhorn_3}
In this section, we turn to the Sinkhorn algorithm for the equivalent OT problem to the LM rate. First, we need to discretize the integral in \eqref{LM}. We might as well consider a set of uniform grid points $\{y_j\}_{j=1}^N$ in $\mathcal{Y}$, and the corresponding rectangular discretization of the integral. This leads to the discrete form of \eqref{LM}:
\begin{subequations}\label{OT}
\begin{align}
    \min_{\textcolor{black}{\gamma_{ij}}}\quad&\sum_{i=1}^{M}\sum_{j=1}^{N}\gamma_{ij}\log \gamma_{ij} \label{OT_obj} \\
    \textbf{s.t.}\quad&\sum_{j=1}^{N}\gamma_{ij} = \mu_i ,\quad i = 1,\cdots,M, \label{OT_Px} \\
    &\sum_{i=1}^{M}\gamma_{ij} = \nu_j,\quad j = 1,\cdots,N, \label{OT_Py} \\
    & \sum_{i=1}^{M}\sum_{j=1}^{N} d_{ij}\gamma_{ij}\le T. \label{OT_ineq}
\end{align}
\end{subequations}
Here $\gamma_{ij}=\gamma(x_i,y_j),\; \nu_j=P_Y(y_j)$ and $d_{ij}=d(x_i,y_j)$. 

\begin{remark}
   In this work, we illustrate our key idea through the rectangular \textcolor{black}{formula} with only first-order accuracy for numerical integration. The idea can be easily generalized to a higher-order numerical discretization, e.g., the trapezoid formula. \textcolor{black}{Not only that, the Sinkhorn algorithm, which we will discuss later in this section, only requires minor modifications to be suitable for higher-order discretization.}
  
  
\end{remark}





Introducing the dual variables $\bda\in\mbbR^{M}, \bdb\in\mbbR^{N}$ and $\lambda\in\mbbR^{+}$, the Lagrangian of \eqref{OT} \textcolor{black}{can be written as}:
\begin{equation}
\begin{aligned}
	&\mathcal{L}(\gamma; \bda, \bdb, \lambda) = \sum_{i=1}^{M}\sum_{j=1}^{N}\gamma_{ij}\log\gamma_{ij} +\sum_{i=1}^{M}\alpha_{i}(\sum_{j=1}^{N}\gamma_{ij} - \mu_i)\\
	& + \sum_{j=1}^{N}\beta_{j}(\sum_{i=1}^{M}\gamma_{ij} - \nu_j)+\lambda(\sum_{i=1}^{M}\sum_{j=1}^{N}d_{ij}\gamma_{ij} - T).
\end{aligned}
\end{equation}
Taking the derivative of $\mathcal{L}(\gamma; \bda, \bdb, \lambda)$ with respect to $\gamma_{ij}$ \textcolor{black}{leads to}
\begin{equation*}
    \gamma_{ij} = \phi_{i}\Lambda_{ij}\psi_{j},
\end{equation*}
in which
\begin{equation*}
    \phi_{i} = \mathrm{e}^{-\alpha_{i}-1/2},\;
    \psi_{j} = \mathrm{e}^{-\beta_{j}-1/2},\;
    \Lambda_{ij} = \mathrm{e}^{-\lambda d_{ij}}.
\end{equation*}
Substituting the above formula into \eqref{OT_Px} and \eqref{OT_Py} yields
\begin{equation*}
\begin{aligned}
&\phi_{i}\sum_{j=1}^{N}\Lambda_{ij}\psi_{j} = \mu_i,\quad  i = 1,\cdots,M,   \\
&\psi_{j}\sum_{i=1}^{M}\Lambda_{ij}\phi_{i} = \nu_j,\quad  j = 1,\cdots,N.  
\end{aligned}
\end{equation*}
Since $\Lambda_{ij}>0$, we can alternatively update $\phi_i$ and $\psi_j$ as follows:
\begin{subequations}\label{sink_iter}
\begin{align}
&\phi_{i}^{(\ell + 1)} = \mu_{i}/\sum_{j=1}^{N}\Lambda_{ij}^{(\ell)}\psi_{j}^{(\ell)},\quad i = 1,\cdots,M, \label{sink_px} \\
&\psi_{j}^{(\ell + 1)} = \nu_{j}/\sum_{i=1}^{M}\Lambda_{ij}^{(\ell)}\phi_{i}^{(\ell+1)},\quad j = 1,\cdots,N. \label{sink_py}
\end{align}
\end{subequations}
This iterative formula is the well-known Sinkhorn algorithm \cite{sinkhorn1967diagonal}. Different from the classical OT problem, we need to deal with the extra constraint \eqref{OT_ineq}. By taking the derivative of $\mathcal{L}(\gamma; \bda, \bdb, \lambda)$ with respect to $\lambda$, we have
\begin{equation}\label{G_lambda}
    G(\lambda)\triangleq \sum_{i=1}^{M}\sum_{j=1}^{N} \phi_{i}\psi_{j} d_{ij} \mathrm{e}^{-\lambda d_{ij}} - T=0.
\end{equation}
Thus, we can update $\lambda$ by finding the roots of
\begin{equation*}
    G(\lambda^{(\ell+1)})\triangleq \sum_{i=1}^{M}\sum_{j=1}^{N} \phi_{i}^{(\ell+1)}\psi_{j}^{(\ell+1)} d_{ij} \mathrm{e}^{-\lambda^{(\ell+1)} d_{ij}} - T.
\end{equation*}
\textcolor{black}{Noticing} that
\begin{equation*}
    G'(\lambda)= -\sum_{i=1}^{M}\sum_{j=1}^{N}\phi_{i}\psi_{j}d_{ij}^{2}\mathrm{e}^{-\lambda d_{ij}} < 0,
\end{equation*}
the function $G(\lambda)$ is monotonic. Thus, The problem of finding the root of $G(\lambda)$ can be easily solved by the Newton's method. The pseudo-code \textcolor{black}{of the proposed algorithm} is present in Algorithm \ref{alg:OT}.


\begin{algorithm}[H]
	\caption{The Sinkhorn Algorithm}
	\label{alg:OT}
	\begin{algorithmic}[1]
		\REQUIRE Decoding metric $d_{ij}$; Marginal distributions $\mu_i, \nu_j$; Iteration number $K$.
	
		\ENSURE Minimal value of $\sum_{i=1}^{M}\sum_{j=1}^{N}\gamma_{ij}\log\gamma_{ij}$.
		\STATE \textbf{Initialization:} $\bm{\phi} = \mathbf{1}_{M}, \bm{\psi} = \mathbf{1}_{N}, \lambda = 1$;
		\FOR{$\ell = 1 : K$}
	
		\STATE $\Lambda_{ij} \gets \mathrm{e}^{-\lambda d_{ij}}$,\quad $i = 1,\cdots,M,\ j = 1,\cdots,N$
		\FOR{$i = 1 : M$}
		\STATE $\phi_{i} \gets \mu_{i}/\sum_{j=1}^{N}\Lambda_{ij}\psi_{j}$
		\ENDFOR
		\FOR{$j = 1 : N$}
		\STATE $\psi_{j} \gets \nu_{j}/\sum_{i=1}^{M}\Lambda_{ij}\phi_{i}$
		\ENDFOR
	
		\STATE Solve $G(\lambda) = 0$ for $\lambda$ with Newton's method
	
		\ENDFOR
		\RETURN $\sum_{i=1}^{M}\sum_{j=1}^{N} \phi_{i}\psi_{j}\Lambda_{ij}\log\left(\phi_{i}\psi_{j} \Lambda_{ij}\right)$
	\end{algorithmic}
\end{algorithm}
Finally, we need to discuss the feasibility of the Sinkhorn iteration, especially \textcolor{black}{under} the extra constraint \eqref{OT_ineq}. Depending on the value of $G(0)$, there are two cases:
\begin{itemize}
    \item $G(0)>0$: In this case, $G(\lambda)=0$ has a unique solution on $(0,+\infty)$ since $G'(\lambda)<0$. Thus the extra constraint \eqref{OT_ineq} is obviously satisfied.
    
    \item $G(0)\le 0$: In this case, the extra constraint \eqref{OT_ineq} is already satisfied. We only need to set $\lambda=0$ instead of solving $G(\lambda)=0$ at  line 8 of Algorithm \ref{alg:OT}.
\end{itemize}











	%






\section{NUMERICAL SIMULATIONS} \label{sec_num_4}
In this section, we use the model and algorithm developed in \textcolor{black}{previous sections} to \textcolor{black}{compute} the LM rate for different modulation schemes over the additive white Gaussian noise (AWGN) channels \textcolor{black}{subject to rotation and scaling} with the following transition law 
\begin{equation}\label{awgneq}
    Y = H X + Z, \ \text{with } Z\sim \mathcal{N}(0,\sigma_{Z}^{2}).
\end{equation}
The matrix $H\in\mbbR^{2\times 2}$ is a combination of rotation and scaling transformations as
\textcolor{black}{
\begin{equation*}
    H = \begin{pmatrix} \eta_1 & 0 \\ 0 & \eta_2 \end{pmatrix}\begin{pmatrix} \cos \theta & \sin \theta\\ -\sin \theta & \cos \theta \end{pmatrix}.
\end{equation*}}
Here the parameters $\eta_1,\;\eta_2$ indicate the scaling of the signal, and the parameter $\theta$ indicates the degree of rotation on the signal. In this work, we consider the case of $\eta_1=1$ and $\eta_2=\eta$. \textcolor{black}{Note that}, the specific values of $\eta$ and $\theta$ can not be priorly known.

The decoding metric is $d(x,y) := \|y-\hat{H}x\|_{2}^{2}$, where $\hat{H}$ is an approximation that is based on part of  the information we know about $H$. We consider the scenario in which the decoder is unaware of the mismatch effect. Thus, in the decoding rule, \textcolor{black}{we have} $\hat{H}=I$.

In the following simulations, we consider two classical modulation schemes: QPSK and 16-QAM. Under the power constraint \eqref{pc}, their constellation points (also known as the alphabet) are illustrated in Fig \ref{Conspoint}. For both cases, all the constellation points in $\mathcal{X}$ are restricted in the region $[-1,1]\times[-1,1]$.


\begin{figure}[H]
	\centerline{\includegraphics[width=0.5\textwidth]{data/Constellations.png}}
	\caption{Constellation points under the power constraint. Left: The QPSK modulation scheme. Right: The 16-QAM modulation scheme.}
	\label{Conspoint}
\end{figure}
Meanwhile, according to the channel law \eqref{awgneq}, the channel output \textcolor{black}{alphabet} $\mathcal{Y}$ \textcolor{black}{is} $\mbbR^2$ itself. Under some reasonable assumptions, e.g. $\eta\le 1,\sigma_z^2\le \frac{1}{2}$, we are able to truncate the \textcolor{black}{alphabet} $\mathcal{Y}$ with a \textcolor{black}{sufficiently large region, e.g.} $[-8, 8]\times[-8, 8]$. Correspondingly, we can discretize the \textcolor{black}{region} by a set of uniform grid points $\{y_j\}_{j=1}^{N}$:
\begin{equation*}
\begin{aligned}
    y_{r\sqrt{N}+s} &= (-8 + r\Delta y, -8 + (s-1)\Delta y), \Delta y = \frac{16}{\sqrt{N}-1},\\
    r &= 0,1,\cdots,\sqrt{N}-1, s = 1,2,\cdots,\sqrt{N},
\end{aligned}
\end{equation*}


In Algorithm \ref{alg:OT}, the termination condition of the Newton iteration is that the update step size of $\lambda$ is less than $10^{-14}$. Moreover, the $\text{SNR} \triangleq 1/(2\sigma_{Z}^{2})$ is used in the sequel. 




All the experiments are conducted on a platform with 128G RAM, and one Intel(R) Xeon(R) Gold 5117 CPU @2.00GHz with 14 cores.


\subsection{Algorithm Verfication}
We first study the convergence of our Sinkhorn algorithm by considering the residual errors of \eqref{sink_iter} and \eqref{G_lambda}:
\begin{equation*}
\begin{aligned}
r_{\phi} &= \sum_{i=1}^{M}|\phi_{i}\sum_{j=1}^{N}\Lambda_{ij}\psi_{j} - \mu_{i}|,\\
r_{\psi} &= \sum_{j=1}^{N}|\psi_{j}\sum_{i=1}^{M}\Lambda_{ij}\phi_{i} - \nu_{j}|,\\
r_{\lambda} &= |G(\lambda)|.
\end{aligned}
\end{equation*}
The parameters are set to
\begin{equation} \label{expe01}
    (\eta,\theta) = (0.9, \pi/18), \; N = 250,000, \;
    \text{SNR} = 0 \textrm{dB}.
\end{equation}
\begin{figure}[H]
	\centerline{\includegraphics[width=0.5\textwidth]{data/Res.png}}
	\caption{The convergent trajectories of the residual error for $r_{\phi}$ (Orange), $r_{\psi}$ (Blue) and $r_{\lambda}$ (Purple). Left: The QPSK modulation scheme. Right: The 16-QAM modulation scheme.}
	\label{Res}
\end{figure}

\begin{table*}[ht]\label{compare}
	\renewcommand\arraystretch{1.25}

	\centering 

	\caption{Comparison between the Sinkhorn algorithm and CVX. Columns 3-5 are the averaged computational time and the speed-up ratio of the Sinkhorn algorithm. Column 6 is the averaged difference of the LM rate computed by two methods.} 
	\label{Table_cvx} 



	\setlength{\tabcolsep}{6.5mm}{
		\begin{tabular}{c|c|c|c|c|c}
		
		
			\toprule
			\multirow{2}{*}{} & \multirow{2}{*}{$N$} & \multicolumn{2}{c|}{Computational time (s)} &
			\multirow{2}{*}{Speed-up ratio} & \multirow{2}{*}{Average difference}\\
			\cline{3-4}
			 &  & Sinkhorn & CVX & &\\%[3pt]只改一行   
			\hline 
		
			\multirow{3}{*}{QPSK} & $100$ & $0.37\times10^{0}$ & $2.61\times10^{1}$ & $7.05\times10^{1}$ & $1.99\times10^{-7}$\\
			 & $225$ & $0.64\times10^{0}$ & $8.66\times10^{1}$ & $1.35\times10^{2}$ & $6.45\times10^{-7}$\\
			 & $400$ & $1.28\times10^{0}$ & - & - & -\\
		
			\hline
		
			\multirow{3}{*}{16-QAM} & $100$ & $0.96\times10^{0}$ & $1.49\times10^{2}$ & $1.55\times10^{2}$ & $4.13\times10^{-7}$\\
			 & $225$ & $1.96\times10^{0}$ & $9.27\times10^{2}$ & $4.73\times10^{2}$ & $1.37\times10^{-6}$\\
			 & $400$ & $3.81\times10^{0}$ & - & - & -\\
		
			\hline
			\multirow{2}{*}{256-QAM} & $100$ & $1.49\times10^{1}$ & - & - & -\\
		
		
			 & $1600$ & $2.24\times10^{2}$ & - & - & -\\
			\bottomrule
			\multicolumn{6}{}{}\\[1pt]
		
		
		
	\end{tabular}}
\end{table*}

In Fig. \ref{Res}, we output the convergent trajectories of the residual errors with respect to iteration steps. We can see that the three curves all decrease rapidly and reach the machine accuracy near $250$ iterations. It is worth mentioning that classical OT problems usually require tens of thousands of Sinkhorn iterations to converge. This significantly illustrates the \textcolor{black}{dramatic} efficiency advantage of our OT model and Sinkhorn algorithm. 

To further illustrate the accuracy and efficiency of the Sinkhorn algorithm for the OT problem \eqref{OT}. We use CVX \cite{2010cvx} as the baseline. The averaged computational time and the averaged difference of the optimal values between the two methods are \textcolor{black}{listed} in Table \ref{Table_cvx}. To reduce the influence of noise, we repeat each experiment for $100$ times. Except for $N$, other parameters are the same as \eqref{expe01}. \textcolor{black}{Since it is difficult for CVX to handle large-scale problems, we restrict $N$ to small scales, e.g. $100,\;225,\;400$ and $1600$.} From the table, we can see the optimal values obtained by the two methods are almost the same. But our Sinkhorn algorithm has a significant advantage (one or two orders of magnitude) in computational speed. Moreover, for slightly larger scale problems, CVX has failed to output convergent results. 





\subsection{Results and Discussions}

Below, we present the computational results of the LM rate under different modulation schemes, different parameters $(\eta,\theta)$, and different SNRs. For comparison, we also present the computational results of GMI \cite{zhang2011general} under the same setup. As we know, GMI is \textcolor{black}{generally lower} than the LM rate \cite{2020Information}. These results not only help us quantitatively understand the relationship between the two rates, but also help verify the correctness of our model and algorithm. In the numerical \textcolor{black}{experiments}, we consider four sets of parameters $(\eta,\theta)$, namely,
\begin{equation*}
   (0.9, \frac{\pi}{18}),\quad
   (0.9, \frac{\pi}{12}),\quad
   (0.8, \frac{\pi}{18}), \quad
   (0.8, \frac{\pi}{12}).
\end{equation*}
We set $N=250,000$ to \textcolor{black}{ensure} the discretization accuracy. And we use $500$ iterations for each \textcolor{black}{experiment} to ensure the algorithm \textcolor{black}{convergence}. We also repeat each \textcolor{black}{experiment} for $100$ times to reduce the influence of noise.
\begin{figure}[H]
	\centerline{\includegraphics[width=0.40\textwidth]{data/QPSK_Fin.png}}
	\caption{LM rate (solid) and GMI (dashed) versus SNR for the QPSK modulation scheme under different mismatched cases, including $(\eta,\theta) = (0.9, \pi/18)$ (Purple), $(\eta,\theta) = (0.8, \pi/18)$ (Blue), $(\eta,\theta) = (0.9, \pi/12)$ (Red), and $(\eta,\theta) = (0.8, \pi/12)$ (Orange).}
	\label{QPSK}
\end{figure}

\begin{figure}[H]
	\centerline{\includegraphics[width=0.40\textwidth]{data/16QAM_Fin.png}}
	\caption{LM rate (solid) and GMI (dashed) versus SNR for the 16-QAM modulation scheme under different mismatched cases, including $(\eta,\theta) = (0.9, \pi/18)$ (Purple), $(\eta,\theta) = (0.8, \pi/18)$ (Blue), $(\eta,\theta) = (0.9, \pi/12)$ (Red), and $(\eta,\theta) = (0.8, \pi/12)$ (Orange).}
	\label{16QAM}
\end{figure}
In Fig. \ref{QPSK}, we can see the comparison of the LM rate and GMI verse SNR for the QPSK modulation scheme. We observe that the LM rate is \textcolor{black}{higher} than GMI with the same parameters. Especially when $\text{SNR} > 8$dB, there is about $2\%$ gain in the LM rate. In \textcolor{black}{addition}, as $\eta$ decreases (from $0.9$ to $0.8$) or $\theta$ increases (from $\pi/18$ to $\pi/12$), both LM rate and GMI decrease accordingly. These results agree with intuition. In Fig. \ref{16QAM}, we \textcolor{black}{display} the comparison for the 16-QAM modulation scheme. From this, we can draw the same conclusions as those for Fig. \ref{QPSK}. 







\section{CONCLUSION} \label{sec_conclu_5}

In this paper, we \textcolor{black}{studied} the computation problem of the LM rate, which is a lower bound for mismatch capacity. Our contributions are twofold. First, we showed that the computation of the LM rate \textcolor{black}{can} be reformulated into the Optimal Transport problem with an extra constraint. Second, we proposed a Sinkhorn-type algorithm to solve the above problem. For the extra constraint, we show that it is equivalent to seeking the root of a one-dimensional monotonic \textcolor{black}{function}. Numerical \textcolor{black}{experiments} show that our approach to computing the LM rate is efficient and accurate. Moreover, we can observe a \textcolor{black}{noticeable} gain in the LM rate compared to GMI.






\bibliographystyle{bibliography/IEEEtran}

