\section{Introduction}
\label{sec:intro}


For a given precoding support $\mathbbmss{X}\subset\mathbbmss{C}$ and penalty function $u(\cdot):\mathbbmss{X}\mapsto\mathbbmss{R}$, the \ac{glse} precoder constructs the transmit vector ${\boldsymbol{x}}\in\mathbbmss{X}^N$ from the data vector ${\boldsymbol{s}}\in \mathbbmss{C}^K$ and the channel matrix $\mathbf{H}\in\mathbbmss{C}^{N\times K}$ as ${\boldsymbol{x}}=\glse{{\boldsymbol{s}},\rho|\mathbf{H}}$ where $\rho$ is a power control factor and \cite{bereyhi2017wsa}
\begin{align}
\glse{{\boldsymbol{s}},\rho|\mathbf{H}}=\mathop{\mathrm{argmin}}_{{\boldsymbol{v}} \in {\mathbbmss{X}^N}} \norm{\mathbf{H} {\boldsymbol{v}}-\sqrt{\rho} \hspace*{.5mm}{\boldsymbol{s}}}^2 +u({\boldsymbol{v}}). \label{eq:glse}
\end{align}
The generality of $\mathbbmss{X}$ and $u(\cdot)$ allows for addressing various forms of constraints on the transmit vector. Compared to the classical approaches for imposing such constraints, the studies in \cite{sedaghat2016lse,sedaghat2017newclass,bereyhi2017wsa,bereyhi2017isit} have shown significant enhancements obtained via the \ac{glse} precoding scheme. Nevertheless, the computational complexity of this scheme has been remained as the main chal- lenge and is intended to be addressed in this paper.

The main motivation of this study comes from the great deal of interest being received recently by massive \ac{mimo} systems \cite{hoydis2013massive}. Form implementational points of view, however, these systems confront the problem of high \ac{rf}-cost which raises due to the vast number of \ac{rf}-chains needed in such setups. The initial approach to overcome this issue is to restrict the \ac{papr} of the transmit vector \cite{mohammed2013per,chen2017low}. In this case, nonlinear power amplifiers with lower dynamic ranges can be employed, and the total \ac{rf}-cost can be significantly reduced. Another approach is \ac{tas} \cite{li2014energy,asaad2017tas} in which a subset of transmit antennas is kept active at each transmission interval, and therefore, the number of required \ac{rf}-chains is reduced. Although such approaches combat the issue of high \ac{rf}-cost, the conventional algorithms significantly degrade the performance. In this case, \ac{glse} precoders reduce this degradation by finding the optimal transmit vector which satisfies the constraints imposed by these approaches.
In general, \ac{glse} precoders solve an optimization problem in each transmission interval. This task is not trivial for choices of $u(\cdot)$ and $\mathbbmss{X}$ which are non-convex. For cases with convex optimization problems, the precoder can be implemented via generic linear programming algorithms. The high computational complexity of these algorithms for large dimensions, however, leaves the implementation of \ac{glse} precoders as an issue in massive \ac{mimo} setups. \ac{gamp} \cite{rangan2011generalized} proposes a low complexity iterative approach for several estimation problems based on approximating the loopy belief propagation algorithm in the large limit \cite{donoho2009message}. The algorithm is known to considerably outperform other available iterative approaches. The underlying estimation problems, which are addressed by \ac{gamp}, are mathematically similar to the \ac{glse} precoding scheme, and therefore, the algorithm can be employed to design a class of iterative precoders based on the \ac{glse} scheme.

The main contribution of this paper is to adopt and tune the \ac{gamp} algorithm to address the \ac{glse} precoding scheme, recently proposed in \cite{sedaghat2016lse,sedaghat2017newclass,bereyhi2017wsa,bereyhi2017isit}. The developed iterative scheme is referred to as ``\ac{glse}-\ac{gamp}'' precoding and exhibits low complexity characteristic. Using the fact that the \ac{glse} and \ac{glse}-\ac{gamp} precoders consider same optimization problems, we further propose a tuning strategy based~on~the~asymptotic results in \cite{sedaghat2016lse,sedaghat2017newclass,bereyhi2017wsa,bereyhi2017isit} derived via the replica method. Our numerical investigations show that the performance of \ac{glse}-\ac{gamp} precoders tuned by the proposed strategy is accurately consistent with asymptotics of corresponding \ac{glse} precoders.


\subsection*{Notation}
Throughout the paper, scalars, vectors and matrices are represented with non-bold, bold lower case and bold upper case letters, respectively. $\mathbf{I}_K$ is a $K \times K$ identity matrix, and $\mathbf{H}^{\mathsf{H}}$ is the Hermitian of $\mathbf{H}$. The set of real and integer numbers are denoted by $\mathbbmss{R}$ and $\mathbbmss{Z}$, and $\mathbbmss{C}$ represents the complex plane. For $s\in\mathbbmss{C}$, $\re{s}$, $\img{s}$ and $\mathbf{s}\coloneqq \left[ \re{s} \hspace*{1mm} \img{s} \right]^\mathsf{T}$ identify the real part, imaginary part and augmented vector, respectively,
and the expression $\mathbf{s} \in \mathbbmss{S}$ indicates that $\mathbf{s}$ is the augmented version of $s\in\mathbbmss{S}$. For $\mathbf{f}({\boldsymbol{x}})=\left[ f_1({\boldsymbol{x}}), \ldots, f_n({\boldsymbol{x}}) \right]^\mathsf{T}$,
the gradient operator is defined as $\gred{{\boldsymbol{x}}} \mathbf{f}({\boldsymbol{x}}) \coloneqq [ \gred{{\boldsymbol{x}}} f_1({\boldsymbol{x}}), \ldots,\gred{{\boldsymbol{x}}} f_n({\boldsymbol{x}}) ]^\mathsf{T}$. $\norm{\cdot}$ and $\norm{\cdot}_1$ denote the Euclidean and $\ell_1$-norm, respectively. Considering the random variable $x$, $\mathrm{p}_x$ represents either the probability mass or density function. Moreover, $\hspace*{-.1mm}\mathbb{E}\hspace*{-.1mm}$~identifies~the expectation. For sake of compactness, $\set{1, \ldots , N}$ is abbreviated by $[N]$, and we define $\tilde{\phi}(x,\lambda)\coloneqq \exp(-{x^2}/{\lambda})$ and
$\tilde{\mathrm{Q}}(x,\lambda) \coloneqq \int_x^\infty \tilde{\phi}(u,\lambda) \mathrm{d} u / \lambda$
for a given non-negative real $\lambda$.


\section{Problem Formulation}
\label{sec:sys}
Consider a Gaussian broadcast \ac{mimo} setup in which a sequence of data symbols $\set{s_k}$ for $k\in [K]$ is transmitted to $K$ single-antenna users simultaneously. The transmitter is equipped with $N$ transmit antennas. The channel is considered to be quasi-static fading and perfectly known at the transmitter. By employing the \ac{glse} precoding scheme given in \eqref{eq:glse} with some penalty $u(\cdot)$ and precoding support $\mathbbmss{X}\subseteq \mathbbmss{C}$, the transmit vector is constructed as ${\boldsymbol{x}}_{N\times 1}=\glse{{\boldsymbol{s}},\rho|\mathbf{H}}$~where~${\boldsymbol{s}}_{K\times 1}\coloneqq \left[ s_1, \ldots, s_K\right]^\mathsf{T}$ 
and $\rho$ is a non-negative power control factor. For this setup, we assume that the following constraints hold.
\begin{enumerate}[label=(\alph*)]
\item ${\boldsymbol{s}}_{K \times 1}$ has \ac{iid} zero-mean complex Gaussian entries with unit variance
\item $u(\cdot)$ decouples meaning that $u({\boldsymbol{v}})=\sum_{j=1}^N u(v_j)$.
\item $N$ and $K$ grow large, such that the load factor $\alpha\coloneqq K/N$ is kept fixed in both $N$ and $K$.
\item $\mathbf{H}^{\mathsf{H}} \mathbf{H} = \mathbf{U} \mathbf{D} \mathbf{U}^{\mathsf{H}}$ in which $\mathbf{U}$ is an $N\times N$ unitary matrix, and $\mathbf{D}$ is a diagonal matrix with asymptotic eigenvalue distribution $\mathrm{p}_\mathbf{D}$. For $\mathrm{p}_\mathbf{D}$, we define the Stieltjes transform as $\mathrm{G}_\mathbf{D}(s)= \mathbb{E} \left\lbrace (d-s)^{-1} \right\rbrace $ with the expectation being taken over $d\sim \mathrm{p}_\mathbf{D}$ and the $\mathrm{R}$-transform as $\mathrm{R}_\mathbf{D} (\omega) = \mathrm{G}_\mathbf{D}^{-1} (-\omega)- \omega^{-1}$ where $\mathrm{G}_\mathbf{D}^{-1} (\cdot)$ denotes the inverse with respect to composition.
\end{enumerate}
By proper choices of the support $\mathbbmss{X}$ and penalty $u(\cdot)$, the \ac{glse} precoder can impose several constraints on the transmit vector.
\begin{itemize}
\item Setting $u({\boldsymbol{v}})=\lambda\norm{{\boldsymbol{v}}}^2$ and $\mathbbmss{X}=\set{x \in \mathbbmss{C}: \abs{x}^2<P}$, the transmit vector is restricted to have a limited \ac{papr}. In fact in this case, the peak power is set to $P$ and a desired constraint on the \ac{papr} is imposed by tuning $\lambda$ such that the average power is accordingly restricted.
\item Let $u({\boldsymbol{v}})\hspace*{-.5mm}=\hspace*{-.5mm}\lambda\norm{{\boldsymbol{v}}}^2+\mu\norm{{\boldsymbol{v}}}_1$ and $\mathbbmss{X} = \mathbbmss{C}$; then, the number of active transmit antennas is constrained.
\end{itemize}


\section{\ac{glse}-\ac{gamp} Precoders}
\label{sec:result}
The \ac{glse} scheme can be considered as a max-sum problem which can be addressed via the \ac{gamp} algorithm~\cite{rangan2011generalized}.
\subsection{\ac{gamp} Algorithm}
The \ac{gamp} algorithm, proposed in \cite{rangan2011generalized}, intends to estimate ${\boldsymbol{v}}_{N\times 1}$ from ${\boldsymbol{s}}_{K\times 1}$ iteratively considering the following setup.
\begin{enumerate}[label=(\alph*)]
\item Each entry of ${\boldsymbol{v}}$ is generated from the corresponding entry of some ${\boldsymbol{a}}\in\mathbbmss{A}^{N}$ via $\mathrm{p}_{v|a}$.
\item The entries of ${\boldsymbol{s}}$ are obtained form the entries of the vector ${\boldsymbol{z}}_{K\times 1}$ through identical scalar channels with $\mathrm{p}_{s|z}$.
\item ${\boldsymbol{z}}$ is a random linear transform of ${\boldsymbol{v}}$, i.e., ${\boldsymbol{z}}=\mathbf{H} {\boldsymbol{v}}$ for some random $K\times N$ matrix $\mathbf{H}$.
\end{enumerate}
Depending on the estimation scheme, the \ac{gamp} algorithm~is developed to address the ``max-sum'' or ``sum-product'' problems. The max-sum \ac{gamp} algorithm iteratively determines the \ac{map} estimation
\begin{align}
{\boldsymbol{x}} = \argmax_{{\boldsymbol{v}}} \sum_{n=1}^N f_\mathrm{in}(v_n,a_n)+\sum_{k=1}^K f_\mathrm{out}(z_k,y_k) \label{eq:map}
\end{align}
for some scalar functions $f_\mathrm{in}(\cdot,\cdot)$ and $f_\mathrm{out}(\cdot,\cdot)$ which represent the conditional distributions $\mathrm{p}_{v|a}$ and $\mathrm{p}_{s|z}$. The sum-product \ac{gamp} algorithm, moreover, addresses the \ac{mmse} estimation where ${\boldsymbol{x}} = \mathbb{E}\set{{\boldsymbol{v}}|{\boldsymbol{s}},{\boldsymbol{a}}}$.
\subsection{The \ac{gamp}-\ac{glse} Algorithm}
By comparing \ac{glse} precoding with \eqref{eq:map}, it is observed that the precoding scheme solves a max-sum problem in which ${\boldsymbol{z}}\coloneqq\mathbf{H} {\boldsymbol{v}}$ with $\mathbf{H}$ being the channel matrix, $v_k \in \mathbbmss{X}$ for $k\in [K]$, and
$f_\mathrm{in}(v_n,a_n) = - u(v_n)$ and
$f_\mathrm{out}(z_k,s_k) = - \abs{z_k-\sqrt{\rho} s_k}^2$.
As the result, the \ac{gamp} algorithm can be applied to iteratively construct the transmit vector ${\boldsymbol{x}}$. By some lines of derivations, the max-sum \ac{gamp} algorithm can be adopted to the \ac{glse} scheme in \eqref{eq:glse}. The resulting algorithm is referred to as ``\ac{glse}-\ac{gamp}'' algorithm and is represented in Algorithm~\ref{A-GAMP} for the precoding support $\mathbbmss{X} \subseteq \mathbbmss{C}$ and the complex-valued matrix $\mathbf{H}$. The variables and functions in the algorithm, for $k\in [K]$ and $n\in [N]$, are defined as follows.
\begin{itemize}
\item The real two-dimensional vectors $\mathbf{w}_k$, ${\mathbf{z}}_k$, $\mathbf{y}_k$, $\mathbf{s}_k$, $\mathbf{u}_n$ and ${\mathbf{x}}_n$ are the augmented forms of the complex scalars $w_k$, $z_k$, $y_k$, $s_k$, $u_n$ and $x_n$, respectively.
\item The matrices $\mathbf{R}^\mathrm{w}_k$, $\mathbf{R}^\mathrm{y}_k$, $\mathbf{R}^\mathrm{u}_n$ and $\mathbf{R}^\mathrm{x}_n$ are real $2\times 2$ matrices, and $\mathbf{Q}_{kn}$ is defined as
\begin{align}
\mathbf{Q}_{kn} \coloneqq 
\begin{bmatrix}
\re{h_{kn}} &-\img{h_{kn}}\\
\img{h_{kn}} &\hphantom{-}\re{h_{kn}}
\end{bmatrix}
\end{align}
with $h_{kn}$ representing the entry $(k,n)$ of $\mathbf{H}$.\vspace*{1mm}
\item $\mathbf{g}_{\mathrm{out}}\left(\cdot \right)$ is the output thresholding function defined as
\begin{align}
\mathbf{g}_{\mathrm{out}}\left(\mathbf{w},\mathbf{s},\mathbf{R}\right) &\coloneqq \gred{\mathbf{w}} \min_{{\mathbf{z}}\in\mathbbmss{C}} \mathcal{E}_{\mathrm{out}} \pn{{\mathbf{z}},\mathbf{w},\mathbf{s},\mathbf{R}} \label{eq:g_out}
\end{align}
where the function $\mathcal{E}_{\mathrm{out}} \pn{\cdot}$ is determined by
\begin{align}
\mathcal{E}_{\mathrm{out}} \pn{{\mathbf{z}},\mathbf{w},\mathbf{s},\mathbf{R}} = \frac{1}{2} &\pn{{\mathbf{z}} - \mathbf{w}}^\mathsf{T}  \mathbf{R}^{-1} \pn{{\mathbf{z}} - \mathbf{w}} \nonumber \\
 &+  \norm{{\mathbf{z}}-\sqrt{\rho} \hspace*{.5mm} \mathbf{s}}^2
\end{align}
\item $\mathbf{g}_{\mathrm{in}}\left(\cdot \right)$ is the input thresholding function~being defined as
\begin{align}
\mathbf{g}_{\mathrm{in}}\left( \mathbf{u},\mathbf{R}\right) \coloneqq \mathop{\mathrm{argmin}}_{{\mathbf{x}}\in\mathbbmss{X}} \mathcal{E}_{\mathrm{in}} \pn{{\mathbf{x}},\mathbf{u}, \mathbf{R}}. \label{eq:g_In}
\end{align}
where the function $\mathcal{E}_{\mathrm{in}} \pn{\cdot}$ is evaluated by
\begin{align}
\mathcal{E}_{\mathrm{in}} \pn{{\mathbf{x}},\mathbf{u},\mathbf{R}} =  \frac{1}{2} \pn{\mathbf{u}-{\mathbf{x}}}^\mathsf{T} \mathbf{R}^{-1} \pn{\mathbf{u}-{\mathbf{x}}} + u({\mathbf{x}}). \label{eq:E_in}
\end{align}
\item The initial conditions are
${\mathbf{x}}_n\pn{1}= \arg \min_{{\mathbf{x}}\in\mathbbmss{X}} u({\mathbf{x}})$ and
$\mathbf{R}_n^{\mathrm{x}}\pn{1}= \left[\gredSq{{\mathbf{x}}} u\pn{\mathbf{u}_n\pn{1}}\right]^{-1}$.
\end{itemize}

The update rules in Algorithm~\ref{A-GAMP} are derived by extending the sum-max \ac{gamp} algorithm to the case with a complex-valued matrix $\mathbf{H}$ and an arbitrary input support $\mathbbmss{X}\subseteq \mathbbmss{C}$. The extension is followed by determining the update rules for the corresponding loopy belief propagation algorithm and then taking some steps similar to \cite[Appendix~C]{rangan2011generalized}. The detailed derivations are skipped due to the page limit and is represented in the extended version of the manuscript.
\begin{remark}
One should distinguish between the  \ac{glse} scheme and the \ac{glse}-\ac{gamp} algorithm. In fact, the former is a least square based scheme to design transmit signals which fulfill some desired constraints. The \ac{glse}-\ac{gamp} algorithm, on the other hand, proposes an iterative approach based on \ac{gamp} to address the \ac{glse} scheme. For some choices of the penalty function, precoding support and channel matrix, the \ac{glse}-\ac{gamp} algorithm converges to the transmit signal given by the \ac{glse} scheme. There are however some particular cases in which the \ac{glse}-\ac{gamp} algorithm does not converge. For these cases, Algorithm~\ref{A-GAMP} does not give the desired transmit signal. To avoid the divergence in such cases, we need to modify the algorithm. This issue is briefly discussed in Section~\ref{sec:num}. 
\end{remark}
In contrast to \ac{glse} precoders, \ac{glse}-\ac{gamp} precoders exhibit low complexity characteristic. Considering Algorithm~\ref{A-GAMP} and noting that the matrices in \eqref{eq:A1}-\eqref{eq:A9} are fixed $2\times 2$ matrices, it is straightforward to show that the total worst-case complexity of \ac{glse}-\ac{gamp} precoders per iteration is $\mathcal{O}(KN)$. The number of iterations, moreover, does not grow with the dimensions. Therefore, one can conclude that the overall complexity of the precoding scheme is $\mathcal{O}(KN)$~as~well.

\begin{figure}
\vspace*{-2.5mm}
\end{figure}
\begin{algorithm}[t]
\caption{ \ac{glse}-\ac{gamp} Precoding Algorithm}
\label{A-GAMP}
\begin{algorithmic}[0]
\Initiate Start from $t=1$ and for $ k \in[K]$ let $\mathbf{y}_k \pn{0} =\boldsymbol{0}$. Set $\mathbf{x}_n\pn{1}$ and $\mathbf{R}_n^{\mathrm{x}}\pn{1}$ for $ n\in[N]$ to their initial conditions. \vspace*{2mm}
\While\NoDo $t< T$
\For\NoDo $k\in[K]$
\begin{subequations}
\begin{align}
\mathbf{R}_k^{\mathrm{w}} \pn{t} &= \sum_{n=1}^N \mathbf{Q}_{k n} \mathbf{R}_n^{\mathrm{x}} \pn{t} \mathbf{Q}_{k n}^\mathsf{T} \label{eq:A1} \\
{\mathbf{z}}_k \pn{t} &=\sum_{n=1}^N \mathbf{Q}_{k n} \mathbf{x}_n \pn{t} \label{eq:A2} \\
\mathbf{w}_k \pn{t} &={\mathbf{z}}_k\pn{t} -\mathbf{R}_k^{\mathrm{w}} \pn{t} \mathbf{y}_k \pn{t-1} \label{eq:A3} \\
\mathbf{y}_k \pn{t} &=\mathbf{g}_{\mathrm{out}}( \mathbf{w}_k \pn{t} ,\mathbf{s}_k,\mathbf{R}_k^{\mathrm{w}} \pn{t} ) \label{eq:A4} \\
\mathbf{R}_k^{\mathrm{y}} \pn{t} &=-\gred{\mathbf{w}} \mathbf{g}_{\mathrm{out}} ( \mathbf{w}_k \pn{t} ,\mathbf{s}_k,\mathbf{R}_k^{\mathrm{w}} \pn{t} )\label{eq:A5}
\end{align}
\end{subequations}
\EndFor
\For\NoDo $n \in[N]$
\begin{subequations}
\begin{align}
\mathbf{R}_n^{\mathrm{u}} \pn{t} &= \left[ \sum_{k=1}^K \mathbf{Q}_{k n}^\mathsf{T} \mathbf{R}_{k}^{\mathrm{y}} \pn{t} \mathbf{Q}_{k n}\right]^{-1} \label{eq:A6} \\
\mathbf{u}_n \pn{t} &=\mathbf{x}_n \pn{t} + \mathbf{R}_n^{\mathrm{u}} \pn{t}\left[ \sum_{k=1}^K \mathbf{Q}_{k n}^\mathsf{T} \mathbf{y}_k \pn{t} \right] \label{eq:A7} \\
\mathbf{x}_n \pn{t+1} &=\mathbf{g}_{\mathrm{in}} ( \mathbf{u}_n \pn{t} ,\mathbf{R}_n^{\mathrm{u}} \pn{t} ) \label{eq:A8} \\
\mathbf{R}_n^{\mathrm{x}} \pn{t+1} &= \left[\gred{\mathbf{u}} \mathbf{g}_{\mathrm{in}} ( \mathbf{u}_n \pn{t} ,\mathbf{R}_n^{\mathrm{u}} \pn{t} ) \right] \mathbf{R}_n^{\mathrm{u}} \pn{t} \label{eq:A9}
\end{align}
\end{subequations}
\EndFor
\EndWhile \vspace*{1mm}
\Out ${\mathbf{x}}_n (T)$ for $n \in [N]$.
\end{algorithmic}
\end{algorithm}

\subsection{Tuning \ac{glse}-\ac{gamp} precoders}
\label{sec:tune_gen}
In order to impose a given set of constraint on the transmit signal, the corresponding \ac{glse}-\ac{gamp} precoder should be tuned. As an example, consider the case in which the number of active transmit antennas, as well as the average transmit power, is desired to be restricted via a \ac{glse}-\ac{gamp} precoder. In this case, one may set $\mathbbmss{X} = \mathbbmss{C}$ and $u({\boldsymbol{v}})=\lambda\norm{{\boldsymbol{v}}}^2+\mu\norm{{\boldsymbol{v}}}_1$. The factors $\lambda$ and $\mu$ in this case control the average transmit power and the fraction of active antennas, respectively. Consequently for given constraints, these factors need to be tuned. Nevertheless, the derivation of an exact tuning strategy is not a trivial problem as the constrained parameters, i.e.,~the~ave- rage power or fraction of active antennas, cannot be derived in terms of the tuning factors straightforwardly. We therefore propose a tuning strategy based on the asymptotics of the \ac{glse}-\ac{gamp} algorithm and its connection to the \ac{glse} scheme. The large-system performance of \ac{glse}-\ac{gamp} precoders is studied through asymptotic analyses of ``state evolution'' equations; see \cite{javanmard2013state} and the references therein. Following the results in the literature, e.g. \cite{rangan2014convergence,rangan2016fixed}, it is shown that for choices of $\mathbf{H}$, $\mathbbmss{X}$ and $u(\cdot)$, in which the \ac{glse}-\ac{gamp} algorithm converges, the asymptotic performance of the algorithm coincides with the large-system performance of \ac{glse} precoders investigated in \cite{bereyhi2017wsa,bereyhi2017isit}. This result indicates that in the large-system limit, the tuning factors for \ac{glse}-\ac{gamp} and \ac{glse} precoders are the same. Therefore, for a given set of constraint, we derive the tuning factors of the \ac{glse}-\ac{gamp} precoders by tuning the corresponding \ac{glse} precoders.
\begin{tun}
Assume that the constraints $f_j({\boldsymbol{x}})/N=C_j$ are desired to be satisfied via a \ac{glse}-\ac{gamp} precoder with penalty $u(\cdot)$ and support $\mathbbmss{X}$ which are controlled by $\lambda_j$~for~$j\in [J]$. Here, $f_j(\cdot)$ are decoupling functions meaning that $f_j({\boldsymbol{x}})=\sum_{n=1}^N f_j(x_n)$. To tune $\lambda_j$ accordingly, we define
\begin{align}
\mathrm{x}=\mathop{\mathrm{argmin}}_{v\in\mathbbmss{X}} \abs{v- s_0}^2+ \xi \hspace*{1mm} u(v) \label{eq:snigle}
\end{align}
where $s_0\sim \mathcal{CN}\left( 0, \sigma^2 \right)$ with
\begin{align}
\sigma^2=\left[\mathrm{R}_\mathbf{D}(-\chi)\right]^{-2}\frac{\partial}{\partial \chi} \left[ ( \lambda_s \chi- \mathsf{p} ) \mathrm{R}_\mathbf{D}(-\chi) \right].
\end{align}
and $\xi=\left[\mathrm{R}_\mathbf{D}(-\chi)\right]^{-1}$ for $\chi$ and $\mathsf{p}$ which satisfy  $\mathsf{p}= \mathbb{E} \abs{\mathrm{x}}^2$~and
\begin{align}
 \frac{\sigma^2 \chi}{\xi} &=  \mathbb{E} \re{\mathrm{x}^* s_0} . \label{eq:fix}
\end{align}
The precoder is then accordingly tuned by choosing $\lambda_j$ for $j\in [J]$ such that the equations  $\mathbb{E} f_j(\mathrm{x})=C_j$ are satisfied.
\end{tun}
\begin{derivation}
The derivation follows~the~marginal decoupling property of the \ac{glse} precoders presented in \cite{bereyhi2017wsa,bereyhi2017isit}. In fact, using the property, it is concluded that $f_j({\boldsymbol{x}})/N$ asymptotically converges to $\mathbb{E} f_j(\mathrm{x})$. By taking the approach illustrated at the beginning of the section, the tuning strategy is obtained.
\end{derivation}
The proposed tuning strategy evaluate the decoupled \ac{glse} precoder\footnote{See Proposition~2 in \cite{bereyhi2017wsa} for the decoupling property of \ac{glse} precoders. A more general version of the property is represented in \cite[Section II-A]{bereyhi2017isit}.} by finding $\chi$ and $\mathsf{p}$ form the fixed-point equations. The asymptotic constrained parameters are then determined by taking the expectation $\mathbb{E} f_j(\mathrm{x})$ and set it equal to $C_j$. One should note that the strategy in general is heuristic, since it tunes the precoders for the large-system limit. Nevertheless, the numerical investigations show that for several cases, the \ac{glse}-\ac{gamp} precoders are well tuned via this strategy.








\section{Applications of \ac{glse}-\ac{gamp} Precoders}
In this section, we investigate two special cases of \ac{glse}-\ac{gamp} precoders with \ac{tas} and limited \ac{papr}. Throughout the analyses, we assume that $\mathbf{H}$ represents an \ac{iid} Rayleigh fading channel with variance $1/N$, i.e., $\mathrm{R}_{\mathbf{J}}(\omega)=\alpha (1-\omega)^{-1}$.

\subsection{\ac{glse}-\ac{gamp} Precoder with \ac{tas}}
\label{sec:tas-gamp}

As it was discussed, \ac{tas} can be directly addressed at the transmit side by using \ac{glse} scheme with $u(v)=\lambda \abs{v}^2 + \mu \abs{v}$.
The corresponding \ac{glse}-\ac{gamp} precoder is therefore given by Algorithm~\ref{A-GAMP}
where $\mathbbmss{X}=\mathbbmss{C}$,
$\mathbf{g}_{\mathrm{out}}\left(\mathbf{w},\mathbf{s},\mathbf{R}\right) = \mathbf{G}_{\mathrm{w}} \mathbf{w} + \mathbf{G}_{\mathrm{s}} \mathbf{s}$ and 
$\gred{\mathbf{w}}\mathbf{g}_{\mathrm{out}}\left(\mathbf{w},\mathbf{s},\mathbf{R}\right) = \mathbf{G}_{\mathrm{w}}$, respectively 
with
\begin{subequations}
\begin{align}
\mathbf{G}_{\mathrm{w}} &\coloneqq-2\hspace*{.5mm}\mathbf{A}^{\mathsf{T}}\mathbf{A}- \pn{\mathbf{A}-\mathbf{I}_2}^\mathsf{T} \mathbf{R}^{-1} \pn{\mathbf{A}-\mathbf{I}_2}  \\
\mathbf{G}_{\mathrm{s}} \hspace*{1mm} &\coloneqq-2\left[ 2\mathbf{A}^{\mathsf{T}}\mathbf{A}\mathbf{R}-\mathbf{A}^{\mathsf{T}}+\pn{\mathbf{A}-\mathbf{I}_2}^\mathsf{T} \mathbf{R}^{-1} \mathbf{A} \mathbf{R}\right]
\end{align}
\end{subequations}
and $\mathbf{A}\coloneqq\pn{\mathbf{I}_2+2\mathbf{R}}^{-1}$. For the input thresholding function, the analytic evaluation of the function from the augmented form in \eqref{eq:g_In} is not trivial. We thus employ the complex scalar form of the equation which results in $\mathbf{g}_{\mathrm{in}}\left( \mathbf{u},\mathbf{R}\right) = \mathbf{G} (\mathbf{u}) \mathbf{f}(\mathbf{u})$ and $\gred{\mathbf{u}} \mathbf{g}_{\mathrm{in}}\left( \mathbf{u},\mathbf{R}\right) = \mathbf{G} (\mathbf{u}) \mathbf{F}(\mathbf{u})$ where
\begin{align}
\mathbf{G} (\mathbf{u}) =
\begin{cases}
    \mathbf{G}_{\mathrm{u}} \hspace*{.3mm}  \vspace*{.8mm}  \qquad & {\norm{\mathbf{u}}} \geq \tau \\
    0             &  {\norm{\mathbf{u}}} < \tau \label{eq:ave_in}
\end{cases}
\end{align}
with $\tau \coloneqq 2 \mu \left[\tr{\mathbf{R}^{-1}}\right]^{-1}$ and $\mathbf{G}_{\mathrm{u}} \coloneqq \pn{\mathbf{I}_2+2\lambda\mathbf{R}}^{-1}$, and
\begin{subequations}
\begin{align}
\mathbf{f}(\mathbf{u}) &\coloneqq \left[1-\frac{\tau}{\norm{\mathbf{u}}}\right] \mathbf{u}, \label{eq:bf_u}\\
\mathbf{F}(\mathbf{u}) &\coloneqq \frac{\tau}{\norm{\mathbf{u}}^3} \mathbf{u}\buu^\mathsf{T} + \left[1-\frac{\tau}{\norm{\mathbf{u}}}\right] \mathbf{I}_2 . \label{eq:F_u}
\end{align}
\end{subequations}
By setting $\mu=0$, the \ac{glse} scheme reduces to \ac{rzf} precoding, and thus, the \ac{glse}-\ac{gamp} algorithm iteratively constructs the output of the \ac{rzf} precoder.
\subsubsection*{Tuning Strategy}
We employ the strategy in Section~\ref{sec:tune_gen} to tune $\mu$ and $\lambda$ such that the fraction of active antennas and the average transmit power are $\eta$ and $P$, respectively. For this case, $J=2$ and $f_1({\boldsymbol{x}})=\norm{{\boldsymbol{x}}}^2$ and $f_2({\boldsymbol{x}})=\norm{{\boldsymbol{x}}}_0$. Consequently, $\lambda$ and $\mu$ are determined from the fixed-point equations $\tilde{\phi}(\xi \mu  ;\theta) = \eta$ for $\theta=(\rho + P)/{\alpha}$ and 
\begin{align}
(1+2\xi \lambda)^2 &=\frac{\theta}{P} \left[\eta- 2 \hspace*{.2mm} \xi \mu \hspace*{.3mm} \tilde{\mathrm{Q}} ({ \xi \mu};\theta ) \right]
\end{align}
and $\xi$ is determined in terms of $\lambda$ and $\mu$ through
\begin{align}
\alpha \xi &= \frac{1}{2} + \frac{\xi}{1+2\xi \lambda} \left[\eta-  \xi \mu \hspace*{.3mm} \tilde{\mathrm{Q}} ({ \xi \mu};\theta ) \right].
\end{align}




\subsection{\ac{glse}-\ac{gamp} Precoder with \ac{papr} Constraint}
\label{sec:papr-gamp}
The precoder in Section~\ref{sec:tas-gamp} can further take the \ac{papr} constraint into account by setting $\mathbbmss{X} \hspace*{-.5mm} = \hspace*{-.5mm} \set{ x\in\mathbbmss{C}:  \abs{x}^2 < P_{\max} }$. The support in this case imposes a peak power constraint on the transmit signal which along with the penalty function restricts both the \ac{papr} and the number of active antennas\footnote{See \cite[Section IV-B]{bereyhi2017wsa} for further illustrations.}. Considering Algorithm~\ref{A-GAMP}, the output function for this setup remains unchanged , and the input function read
\begin{align}
\mathbf{g}_{\mathrm{in}}\left(\mathbf{u},\mathbf{R}\right) =
\begin{cases}
  \dfrac{{\mathbf{u}}}{{\norm{\mathbf{u}}}} \sqrt{P_{\max}}  \vspace*{1mm}   & \tilde{\tau} \leq {\norm{\mathbf{u}}} \\
    \mathbf{G}_{\mathrm{u}} \hspace*{.3mm} \mathbf{f}(\mathbf{u}) \vspace*{.8mm}  \qquad &\tau \leq {\norm{\mathbf{u}}} < \tilde{\tau} \\
    0             & 0 \leq {\norm{\mathbf{u}}} < \tau \label{eq:peak_in}
\end{cases}
\end{align}
with the corresponding gradient
\begin{align}
\gred{\mathbf{u}}\mathbf{g}_{\mathrm{in}}\left(\mathbf{u},\mathbf{R}\right) =
\begin{cases}
  \dfrac{\tilde{\mathbf{u}}\tilde{\mathbf{u}}^\mathsf{T}}{{\norm{\mathbf{u}}}^3} \sqrt{P_{\max}} \vspace*{.8mm}    & \tilde{\tau} \leq {\norm{\mathbf{u}}} \\
    \mathbf{G}_{\mathrm{u}} \hspace*{.3mm}\mathbf{F}(\mathbf{u})\vspace*{.8mm}   \qquad &\tau \leq {\norm{\mathbf{u}}} < \tilde{\tau} \\
    0             &0\leq {\norm{\mathbf{u}}} < \tau, \label{eq:peak_grad_in}
\end{cases}
\end{align}
where $\tilde{\mathbf{u}}\coloneqq \left[\mathrm{u}_2,-\mathrm{u}_1\right]^\mathsf{T}$, $\tau \coloneqq 2 \mu \left[\tr{\mathbf{R}^{-1}}\right]^{-1}$, and
\begin{align}
\tilde{\tau} \coloneqq \left(1+\frac{4\lambda}{\tr{\mathbf{R}^{-1}}} \right) \sqrt{P_{\max}} + \frac{2\mu}{\tr{\mathbf{R}^{-1}}}.
\end{align}
$\mathbf{G}_\mathrm{u}$, $\mathbf{f}(\mathbf{u})$ and $\mathbf{F}(\mathbf{u})$ are moreover given as in Section~\ref{sec:tas-gamp}.
By setting $\mu=0$, the~precoder~employs all the transmit antennas and restricts only the \ac{papr}. In this case, $\mathbf{F}(\mathbf{u})=\mathbf{I}_2$, $\mathbf{f}(\mathbf{u})=\mathbf{u}$, and $\tau$ reduces to zero.
\subsubsection*{Tuning Strategy}
Consider the same constraints as for the case without the \ac{papr} restriction. From Section~\ref{sec:tune_gen}, $\lambda$ and $\mu$ for the average power $P$ and the fraction of active antennas~$\eta$ are given by the fixed-point equations $\tilde{\phi}(\xi \mu  ;\theta) = \eta$ and
\begin{align}
(1+2\xi \lambda)^2 &=\frac{\theta}{P} \left[\Delta_1 (\xi \mu) - 2 \hspace*{.2mm} \xi \mu \hspace*{.3mm} \Delta_2 (\xi \mu)\right].
\end{align}
Here, $\xi$ is a function of $\lambda$ and $\mu$ which satisfies
\begin{align}
\alpha \xi &= \frac{1}{2} + \frac{\xi}{1+2\xi \lambda} \left[\Delta_1 (\xi \mu) - 2 \hspace*{.2mm} \xi \mu \hspace*{.3mm} \Delta_2 (\xi \mu)\right].
\end{align}
Moreover, $\theta=(\rho + P)/{\alpha}$ and we have defined
\begin{subequations}
\begin{align}
\Delta_1(\xi \mu) &\coloneqq \tilde{\phi}(\xi \mu  ;\theta) - \tilde{\phi}(\xi \mu + (1+2\xi \lambda)\sqrt{P_{\max}}  ;\theta), \\
\Delta_2(\xi \mu) &\coloneqq \tilde{\mathrm{Q}} (\xi \mu  ;\theta) - \tilde{\mathrm{Q}} (\xi \mu + (1+2\xi \lambda)\sqrt{P_{\max}}  ;\theta) .
\end{align}
\end{subequations}

\section{Numerical Investigations}
\label{sec:num}
To investigate the performance of \ac{glse}-\ac{gamp} precoders, we define the distortion measure for a given $\rho$ as
\begin{align}
\mathsf{D}(\rho) \coloneqq \frac{1}{K} \mathbb{E} \norm{\mathbf{H}{\boldsymbol{x}}-\sqrt{\rho}\hspace*{.2mm} {\boldsymbol{s}}}^2 \label{eq:D_rho}
\end{align}
which determines the average distortion caused by the multiuser interference at receive terminals. It is moreover shown that the achievable ergodic rate per user can be bounded from below in terms of $\mathsf{D}(\rho)$ as proved in \cite{sedaghat2016lse}.

\begin{figure}[t]
\hspace*{-.8cm}  
\resizebox{1.085\linewidth}{!}{
\pstool[width=.35\linewidth]{Figs/Fig1.eps}{
\psfrag{D}[c][t][0.23]{$\mathsf{D}(\rho=1)$ in [dB]}
\psfrag{alpha-inv}[c][c][0.25]{$\alpha^{-1}$}
\psfrag{PP=3-070}[l][l][0.24]{$\eta=0.7$}
\psfrag{PP=3-100}[l][l][0.24]{$\eta=1$}
\psfrag{PP=3-030}[l][l][0.24]{$\eta=0.3$}
\psfrag{PP=3-050}[l][l][0.24]{$\eta=0.5$}
\psfrag{-4}[r][c][0.18]{$-4$}
\psfrag{-8}[r][c][0.18]{$-8$}
\psfrag{-12}[r][c][0.18]{$-12$}
\psfrag{-16}[r][c][0.18]{$-16$}
\psfrag{-20}[r][c][0.18]{$-20$}

\psfrag{2}[c][b][0.18]{$2$}
\psfrag{2.5}[c][b][0.18]{$2.5$}
\psfrag{3}[c][b][0.18]{$3$}
\psfrag{3.5}[c][b][0.18]{$3.5$}
\psfrag{4}[c][b][0.18]{$4$}
\psfrag{4.5}[c][b][0.18]{$4.5$}
\psfrag{5}[c][b][0.18]{$5$}
}}
\caption{Distortion at $\rho=1$ vs. $\alpha^{-1}$ for $P=0.3$ and various $\eta$. Circles depict the performance of the \ac{glse}-\ac{gamp} precoder for $N=64$ and $T=20$. Solid lines denote the asymptotic performance of the corresponding \ac{glse} scheme determined by the replica method.
\vspace*{-3.5mm}}
\label{fig:1}
\end{figure}

The circles in Fig.~\ref{fig:1} show the distortion given by the \ac{glse}-\ac{gamp} precoder presented in Section~\ref{sec:tas-gamp} for various inverse load factors $\alpha^{-1}=N/K$ considering several constraints on the number of active antennas. The results have been given for $N=64$ antennas and $T=20$ iterations. The asymptotic performances of the corresponding \ac{glse} precoders, derived via the replica method in \cite{bereyhi2017isit}, have been also sketched with solid lines.
Here, $\rho=1$ and $\lambda$ is set such that $P=0.3$. As the figure shows, the \ac{glse}-\ac{gamp} precoder tracks accurately the performance of the \ac{glse} scheme, even for a practically moderate number of antennas.
For the \ac{papr}-limited precoder in Section~\ref{sec:papr-gamp}, the distortion at $\rho=1$ has been plotted in terms of $\alpha^{-1}$ in Fig.~\ref{fig:2}. The curves have been sketched for multiple \ac{papr} constraints. Similar to Fig.~\ref{fig:1}, solid lines correspond to the \ac{glse} scheme and circles denote the simulation results for the \ac{glse}-\ac{gamp} precoder with $N=64$ and $T=20$ for \ac{papr} $=3$ dB. Here, we have considered $P=0.5$, and $P_{\max}$ is tuned via the proposed strategy assuming all the antennas being active.
The figure depicts that by increasing the \ac{papr} up to $5$ dB, the performance of the precoder is sufficiently close to the case without \ac{papr} restriction. This observation suggests for employing the \ac{glse}-\ac{gamp} precoder, in order to reduce the transmit \ac{papr} without any significant performance loss. In this case, low efficiency power amplifiers can be utilized which can significantly reduce the \ac{rf}-cost

\begin{remark}
\label{remark:1}
It is known that the \ac{gamp} algorithm converges for \ac{iid} Gaussian matrices \cite{rangan2014convergence,rangan2016fixed}. However, by deviating from this assumption, the algorithm may diverge. This issue was recently addressed in \cite{rangan2017vector} via the \ac{vamp} algorithm. Consequently, for channel models with ill-conditioned matrices, one can develop a precoding algorithm based on the \ac{glse} scheme by taking a same approach while employing \ac{vamp}.
\end{remark}


\section{Conclusion}
\label{conclusion}
This paper has proposed a class of low complexity precoders based on the \ac{glse} scheme using the \ac{gamp} algorithm.
The numerical investigations have been consistent with the replica results for the \ac{glse} scheme given in \cite{sedaghat2016lse,sedaghat2017newclass,bereyhi2017wsa,bereyhi2017isit}. This consistency demonstrates that various implementational limitations in massive \ac{mimo} systems can be effectively overcome using some low-complexity, but effective, algorithms. As indicated in Remark~\ref{remark:1}, the \ac{glse}-\ac{gamp} precoders may fail in converging for channel models with ill-conditioned channel matrices, and therefore, an alternative algorithm can be proposed via \ac{vamp}. The extension under \ac{vamp} is however skipped and left as a possible future work.

\begin{figure}[t]
\hspace*{-.8cm}  
\resizebox{1.085\linewidth}{!}{
\pstool[width=.35\linewidth]{Figs/Fig2_copy_2.eps}{
\psfrag{D}[c][t][0.23]{$\mathsf{D}(\rho=1)$ in [dB]}
\psfrag{alpha-inv}[c][c][0.25]{$\alpha^{-1}$}
\psfrag{PAPR=3-AAA-BBB}[l][l][0.24]{${\mathsf{PAPR}}=3$ dB}
\psfrag{PAPR=0-AAA-BBB}[l][l][0.24]{${\mathsf{PAPR}} \downarrow 0$ dB}
\psfrag{PAPR=5-AAA-BBB}[l][l][0.24]{${\mathsf{PAPR}}=5$ dB}
\psfrag{PAPR=9-AAA-BBB}[l][l][0.24]{${\mathsf{PAPR}}\uparrow\infty$}
\psfrag{-4}[r][c][0.18]{$-4$}
\psfrag{-8}[r][c][0.18]{$-8$}
\psfrag{-12}[r][c][0.18]{$-12$}
\psfrag{-16}[r][c][0.18]{$-16$}
\psfrag{-20}[r][c][0.18]{$-20$}

\psfrag{2}[c][b][0.18]{$2$}
\psfrag{2.5}[c][b][0.18]{$2.5$}
\psfrag{3}[c][b][0.18]{$3$}
\psfrag{3.5}[c][b][0.18]{$3.5$}
\psfrag{4}[c][b][0.18]{$4$}
\psfrag{1.5}[c][b][0.18]{$1.5$}
\psfrag{1}[c][b][0.18]{$1$}
}}
\caption{Distortion at $\rho=1$ vs. $\alpha^{-1}$ for several \ac{papr}s. $P=0.5$ and $\eta=1$. Solid lines and circles respectively denote the results for \ac{glse} and \ac{glse}-\ac{gamp} algorithm with $N=64$ and $T=20$.
\vspace*{-4mm}}

\label{fig:2}
\end{figure}


