Low-Complexity QL-QR Decomposition Based Beamforming Design for Two-Way MIMO Relay Networks
Wei Duan
~\IEEEmembership{Student Member
~IEEE
%\boldmath In this paper, we investigate the optimization problem of joint source and relay beamforming matrices for a two-way amplify-and-forward (AF) multi-input multi-output (MIMO) relay system. The system consisting of two source nodes and two relay nodes is considered and the linear minimum mean-square-error (MMSE) is employed at both receivers. We assume individual relay power constraints and study an important design problem, a so-called determinant maximization (DM) problem. Since this DM problem is nonconvex, we consider an efficient iterative algorithm by using an MSE balancing result to obtain at least a locally optimal solution. The proposed algorithm is developed based on QL, QR and Choleskey decompositions which differ in the complexity and performance. Analytical and simulation results show that the proposed algorithm can significantly reduce computational complexity compared with their existing two-way relay systems and have equivalent bit-error-rate (BER) performance as the singular value decomposition (SVD) based on a regular block diagonal (RBD) scheme.
MSE_Duality.tex
\documentclass[journal, two column]{IEEEtran}
%\documentclass[journal, one column]{IEEEtran}

\usepackage{cite}
\usepackage{graphicx}
\usepackage{psfrag}
\usepackage{subfigure}
\usepackage{url}
\usepackage{stfloats}
\usepackage{amsmath}
\usepackage{array}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{algorithmicx}
\usepackage[ruled]{algorithm}
\usepackage{algpseudocode}
\usepackage{algpascal}
\usepackage{algc}


%\usepackage{anysize}


\hyphenation{op-tical net-works semi-conduc-tor}

%\renewcommand\baselinestretch{2}



\begin{document}


\title{Low-Complexity QL-QR Decomposition Based Beamforming Design for Two-Way MIMO Relay Networks}
\author{Wei Duan,~\IEEEmembership{Student Member,~IEEE} Wei Song and Moon Ho Lee,~\IEEEmembership{Senior Member, IEEE}
        % <-this % stops a space
%\thanks{The work presented in this paper was supported by
%MEST 2012-002521 NRF.}% <-this % stops a space~\IEEEmembership{Student Member,~IEEE},~\IEEEmembership{Senior Member, IEEE}
\thanks{Wei Duan,  and Moon Ho Lee are with the Institute of Information and Communication, Chonbuk National University, 664-14 Deokjin-dong, Jeonju 561-756, Republic of Korea.(e-mail:
sinder@live.cn, moonho@jbnu.ac.kr)}
%\thanks{Wei Song, Sang Seob Song, ssong@jbnu.ac.kr,  is with College of Information Technology, Eastern Liaoning University, Dandong, 118003, China. (e-mail:
%sw658@hotmail.com)}
\thanks{Wei Song is with College of
Information Technology, Eastern liaoning university,  Dandong,
118003, P.R.China. (e-mail: sw658@hotmail.com)}% <-this % stops a space
} \maketitle




\begin{abstract}
%\boldmath
In this paper, we investigate the optimization problem of joint
source and relay beamforming matrices for a two-way
amplify-and-forward (AF) multi-input multi-output (MIMO) relay
system. The system consisting of two source nodes and two relay
nodes is considered and the linear minimum mean-square-error
(MMSE) is employed at both receivers. We assume individual relay
power constraints and study an important design problem, a
so-called determinant maximization (DM) problem. Since this DM
problem is nonconvex, we consider an efficient iterative algorithm
by using an MSE balancing result to obtain at least a locally
optimal solution. The proposed algorithm is developed based on QL,
QR and Choleskey decompositions which differ in the complexity and
performance. Analytical and simulation results show that the
proposed algorithm can significantly reduce computational
complexity compared with their existing two-way relay systems and
have equivalent bit-error-rate (BER) performance as the singular
value decomposition (SVD) based on a regular block diagonal (RBD)
scheme.


\end{abstract}
\textbf{Index Terms:}  Two-way relay channel, MIMO,  QL-QR
decomposition, Choleskey decomposition, determinant maximization,
amplify-and-forward.





\section{Introduction}

Recently, wireless relay networks have been the focus of a lot of
research because the relaying transmission is a promising
technique which can be applied to extend the coverage or increase
the system capacity. There are various cooperative relaying
schemes have been proposed, such as amplify-and-forward (AF)
\cite{b1} and \cite{c2}, decode-and-forward (DF) \cite{b2},
denoise-and-forward (DNF) \cite{b3}, and compress-and-forward (CF)
\cite{b4} cooperative relaying protocols. Among these approaches,
AF is most widely used due to without detecting the transmitted
signal. Therefore, an AF relay scheme requires a less processing
power at the relays compared to other schemes.

In one-way relaying (OWR) approach, to completely exchange
information between two base stations, four time slots are
required in uplink (UL) and downlink (DL) communications, which
leads to a loss of one-half spectral resources \cite{b5}. In order
to solve this problem, a two-way relaying approach has been
considered in \cite{b6}, \cite{b7}, and \cite{b8}. In a typical
two-way relaying scheme, the communication is completed in two
steps. First, the transmitters send their symbols to two relays,
simultaneously. After receiving the signals, each relay processes
them based on an efficient relaying scheme to produce new signals.
Then the processed signals are broadcasted to both receiver nodes.

Multi-input multi-output (MIMO) relay systems have been
investigated in [10]--[13]. It is shown that, by employing
multiple antennas at the transmitter and/or the receiver, one can
significantly improve the transmission reliability by leveraging
spatial diversity. Relay precoder design methods have been
investigated in \cite{b13}\raisebox{0.1mm}, \cite{b14},
\cite{b15}. A problem in designing optimal beamforming vectors for
multicasting is challenging due to its nonconvex nature. In
\cite{b13}, the authors propose a transceiver precoding scheme at
the relay node by using zero-forcing (ZF) and MMSE criteria with
certain antenna configurations. The information theoretic capacity
of the multi-antenna multicasting is studied in \cite{b14}, along
with the achievable rates using lower complexity transmission
schemes, as the number of antennas or users goes to infinity. In
\cite{b15}, the authors propose an alternative method to
characterize the capacity region of two-way relay channel (TWRC)
by applying the idea of rate profile.

Joint optimization of the relay and source nodes for the MIMO TWRC
have been studied in \cite{b8}, \cite{b16}. In \cite{b8}, the
authors develop a unified framework for optimizing two-way linear
non-regenerative MIMO relay systems and show that the optimal
relay and source matrices have a general beamforming structure.
The joint source node and relay precoding design for minimizing
the mean squared error in a MIMO two-way relay (TWR) system is
studied in \cite{b16}.

Since singular value decomposition (SVD) and/or generalized SVD
(GSVD) are widely used to find the orthogonal complement to solve
an optimization problem \cite{c2}, \cite{b8}, \cite{b15},
\cite{b19}, but their computational complexities are extremely
high. In order to reduce the complexity, the SVD can be replaced
with a less complex QR decomposition \cite{f1} in this work.
However, this approach leads to degrading the BER performance. In
addition, it is difficult to realize in TWRC. In this paper, we
investigate the joint source and relay precoding matrix
optimization for a two way-relay amplify-and-forward relaying
system where two source nodes and two relay nodes are equipped
with multiple antennas. Also, in order to apply the QL/QR
decomposition to the TWRC, we design a three part relay filter.
Compared with existing works such as [9]--[14], the contributions
of this paper can be summarized as follows. Firstly, we
investigate a two-way MIMO relay system using the criteria which
minimizes an MSE of the signal waveform estimation for both two
source nodes. We prove an optimal sum-MSE solution can be obtained
as the Winer filter while signal-to-noise-ratio (SNR) at both
source nodes are equivalent \cite{b17} which leading to an MSE
balancing result. Secondly, we propose a new cooperative scenario,
i.e., the QL-QR compare with the Choleskey decomposition which
significantly reduces the computational complexity of the optimal
design. In this proposed design, the channels of its left side are
decomposed by the QL decomposition while those of its right side
factorized by the QR decomposition. And the equivalent noise
covariance is decomposed by the Choleskey decomposition. We also
design the three part relay filter, which is comprised by a left
filer, a middle filter, and a right filter, to efficiently combine
two source nodes and the relay nodes. By these approaches, the
received signals at both two source nodes are able to be redeemed
as either lower or upper triangular matrices. Stemming from one of
the properties of triangular matrices such that their determinant
is identical to the multiplication of their eigenvalues, we are
able to straightforwardly solve the optimization problem as a
determinant maximization problem. Also, we can obtain the BER
performance equivalent to that of SVD-RBD scheme.

The rest of this paper is organized as follows. Section II
describes a system model of the TWRC and raises a sum-MSE problem.
In Section III, we propose an iterative QL-QR algorithm and a
joint optimal beamforming design. In Sections IV, we discuss the
computational complexity of an efficient channel model. The
simulation results are presented to show the excellent performance
of our proposed algorithm for the TWRC in Section V. Section VI
concludes this paper.

\emph{Notations:} $\mathbf{A}^*$ and $\mathbf{A}^T$ denote the
conjugate and the transpose of a matrix $\mathbf{A}$,
respectively. $\mathbf{d}(.)$ denotes a diagonal matrix and an
$N\times N$ identity matrix is denoted by $\mathbf{I}_N$.
$\mathrm{E}(.)$ stands for the statistical expectation and
$(.)^{H}$ denotes a Hermitian transpose of a given matrix. $tr(.)$
and $R(.)$ denotes matrix trace and the range of a matrix.

\begin{figure}
\begin{center}
\includegraphics [width=60mm,height=34mm]{TWR3.eps}
\caption{\label{1} Proposed LQ-QR Amplify-and-Forward MIMO TWR
system.}
\end{center}
\end{figure}


\section{System Model and Sum-MSE}
 We consider a TWRC consisting of two source nodes
$\mathbf{S}_{1}$ and $\mathbf{S}_{2}$, and two relay nodes
$\mathbf{R}_{1}$ and $\mathbf{R}_{2}$  as shown in Fig. \ref{1}.
The source and relay nodes are equipped with \emph{M} and \emph{N}
antennas, respectively. We adopt the relay protocol with two time
slots introduced in \cite{b13}. In the first time slot, the
information vector $x_{i}\in{\mathbb{C}^{M\times{1}}}$ is linearly
processed by a precoding matrix
$\mathbf{V}_{i}\in{\mathbb{C}^{M\times{M}}}$ and then be
transmitted to the relay nodes. The received signals at
$\mathbf{R}_{i}$, $i\in{\{1,2\}}$, can be expressed as
\begin{eqnarray}\label{t1}
&&y_{R_{1}}=\textbf{H}_{1,1}s_{1}+\textbf{H}_{1,2}s_{2}+n_{R_{1}}\notag\\
&&y_{R_{2}}=\textbf{H}_{2,1}s_{1}+\textbf{H}_{2,2}s_{2}+n_{R_{2}},
\end{eqnarray}
where $y_{R_{i}}\in{\mathbb{C}^{N\times1}}, i\in\{1,2\}$,
indicates the received signal vector,
$\textbf{H}_{i,j}\in{\mathbb{C}^{N\times{M}}}, i,j\in\{1,2\}$,
represents the channel matrix from source $j$ to relay $i$, as
shown in Fig.1, $s_{i}\in{\mathbb{C}^{M\times1}}$ is the
transmitted symbol vector from $\mathbf{S}_{i}$, and
$n_{R_{i}}\sim{CN(0,\sigma^{2}_{R_{i}}\mathbf{I}_{N})}$ represents
the additive white Gaussian noise (AWGN) vector with zero mean and
variance $\sigma^{2}_{R_{i}}$ at relay node $i$. The term $s_{i}$
is subject to a power constraint,
$tr\{E(s_{i}s^{H}_{i})\}\leq{P_{i}}$ with
$tr\{E(x_{i}x^{H}_{i})\}\leq{\frac{P_{i}}{M}\mathbf{I}_{M}}$,
where $P_{i}$ is the transmit power at $\mathbf{S}_{i}$.

To find an appropriate power normalization vector $\rho_{R}$, we
express the total transmission power at the source node with
$\mathbf{V}_{i}$ as
\begin{eqnarray}\label{Power nomarlizeition}
tr\{\mathbf{V}_{1}\mathbf{V}^{H}_{1}\!+\!\mathbf{V}_{2}\mathbf{V}^{H}_{2}\}\!\!\!\!&=&\!\!\!\!tr\left\{\rho^{2}_{R}\left(\textbf{V}^{b}_{1}\left(\textbf{V}^{b}_{1}\right)^{H}+\textbf{V}^{b}_{2}\left(\textbf{V}^{b}_{2}\right)^{H}\right)\right\}\notag\\
&=&\!\!\!\!tr\left\{\rho^{2}_{R}\left(P_{1}+P_{2}\right)\right\}.
\end{eqnarray}
In this paper, we assume that each transmit antenna satisfies the
unity transmission power constraint. To satisfy the power
constraint, we propose the following power normalization vector
\begin{eqnarray}\label{Power nomarlizeition1}
\rho_{R}=1/\sqrt{P_{1}+P_{2}}.
\end{eqnarray}
In the second time slot, after power normalization, the relay node
$\mathbf{R}_{i}$ linearly amplifies $y_{R_{i}}$ with an
$N\times{N}$ matrix $\mathbf{F}_{i}$ and then broadcasts the
amplified signal vector $x_{R_{i}}$ to source nodes 1 and 2. The
signals transmitted from relay node $i$ can be expressed as
\begin{eqnarray}\label{R power}
x_{R_{i}}=\rho_{R}\textbf{F}_{i}y_{R_{i}}.
\end{eqnarray}
Using \eqref{t1} and \eqref{R power}, the received signal vectors
at $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ can be, respectively,
written as
\begin{eqnarray}\label{s2}
y_{1}&=&\textbf{H}^{T}_{1,1}\textbf{F}_{1}\textbf{H}_{1,2}s_{1}+\textbf{H}^{T}_{1,1}\textbf{F}_{1}\textbf{H}_{1,2}s_{2}
+\textbf{H}^{T}_{2,1}\textbf{F}_{2}\textbf{H}_{2,2}s_{1}\notag\\
&&+\textbf{H}^{T}_{2,1}\textbf{F}_{2}\textbf{H}_{2,2}s_{2}
+\textbf{H}^{T}_{1,1}\textbf{F}_{1}n_{R_{1}}
+\textbf{H}^{T}_{2,1}\textbf{F}_{2}n_{R_{2}}+n_{1}\notag\\
y_{2}&=&\textbf{H}^{T}_{1,2}\textbf{F}_{1}\textbf{H}_{1,1}s_{1}+\textbf{H}^{T}_{1,2}\textbf{F}_{1}\textbf{H}_{1,1}s_{2}
+\textbf{H}^{T}_{2,2}\textbf{F}_{2}\textbf{H}_{2,1}s_{1}\notag\\
&&+\textbf{H}^{T}_{2,2}\textbf{F}_{2}\textbf{H}_{2,1}s_{2}+\textbf{H}^{T}_{1,2}\textbf{F}_{1}n_{R_{1}}+\textbf{H}^{T}_{2,2}\textbf{F}_{2}n_{R_{2}}+n_{2},\notag\\
\end{eqnarray}
where ${\textbf{H}}^{T}_{i,j}$, $i,j\in\{1,2\}$, indicates the
$M\times{N}$ channel matrix from the relay node $i$ to the source
node $j$, and $n_{i}$, $i\in\{1,2\},$ is an $M\times 1$ noise
vector at $\mathbf{S}_{i}$.
%Here
%$\overline{i}=2$ for $i=1$, and $\overline{i}=1$ for $i=2$, $E(.)$
%stands for the statistical expectation, $I_{M}$ is and $M\times M$
%indentity matrix, and $(.)^{H}$ denotes matrix Hermitian
%transpose.
%In this paper, we assume that all MIMO channels are quasi-static,
%that is they remain constant over one time frame, but can change
%to another value in the next time frame.
We assume that the relay nodes perfectly know the channel state
information (CSI) of $\textbf{H}_{i,j}$. The relay node
$\mathbf{R}_{i}$ performs the optimizations of $\textbf{F}_{i}$
and $\mathbf{V}_{i}$, and then transmits the information to the
source nodes 1 and 2. Since source node $i$ knows its own
transmitted signal vector $s_{i}$ and full CSI, the
self-interference components in \eqref{s2} can be efficiently
cancelled. The effective received signal vectors are given by
\begin{eqnarray}\label{s3}
\widetilde{y}_{1}&=&\textbf{H}^{T}_{1,1}\textbf{F}_{1}\textbf{H}_{1,2}s_{2}+\textbf{H}^{T}_{2,1}\textbf{F}_{2}\textbf{H}_{2,2}s_{2}+\textbf{H}^{T}_{1,1}\textbf{F}_{1}n_{R_{1}}\notag\\
&&+\textbf{H}^{T}_{2,1}\textbf{F}_{2}n_{R_{2}}+n_{1}\notag\\
&=&\mathbf{\widetilde{H}}_{1}s_{2}+\widetilde{n}_{1},
\end{eqnarray}
\begin{eqnarray}\label{a}
\widetilde{y}_{2}&=&\textbf{H}^{T}_{1,2}\textbf{F}_{1}\textbf{H}_{1,1}s_{1}+\textbf{H}^{T}_{2,2}\textbf{F}_{2}\textbf{H}_{2,1}s_{1}+\textbf{H}^{T}_{1,2}\textbf{F}_{1}n_{R_{1}}\notag\\
&&+\textbf{H}^{T}_{2,2}\textbf{F}_{2}n_{R_{2}} +n_{2}\notag\\
&=&\mathbf{\widetilde{H}}_{2}s_{1}+\widetilde{n}_{2},
\end{eqnarray}
where
$\widetilde{\textbf{H}}_{1}=\textbf{H}^{T}_{1,1}\textbf{F}_{1}\textbf{H}_{1,2}+\textbf{H}^{T}_{2,1}\textbf{F}_{2}\textbf{H}_{2,2}$
and
$\widetilde{\textbf{H}}_{2}=\textbf{H}^{T}_{1,2}\textbf{F}_{1}\textbf{H}_{1,1}+\textbf{H}^{T}_{2,2}\textbf{F}_{2}\textbf{H}_{2,1}$
are the equivalent MIMO channels seen at source nodes
$\mathbf{S}_{1}$ and $\mathbf{S}_{2}$, respectively. The vectors
$\widetilde{n}_{1}=\textbf{H}^{T}_{1,1}\textbf{F}_{1}n_{R_{1}}+\textbf{H}^{T}_{2,1}\textbf{F}_{2}n_{R_{2}}+n_{1}$
and
$\widetilde{n}_{2}=\textbf{H}^{T}_{1,2}\textbf{F}_{1}n_{R_{1}}+\textbf{H}^{T}_{2,2}\textbf{F}_{2}n_{R_{2}}
+n_{2}$ are the equivalent noises at source node $\mathbf{S}_{1}$
and $\mathbf{S}_{2}$, respectively.

Due to the lower computational complexity, linear receivers are
applied at source node $i$ to retrieve the transmitted signals
sent from the other nodes. The estimated signal waveform vector is
given as
$\widehat{s}_{\overline{i}}=\mathbf{W}^{H}_{i}\widetilde{y}_{i}$,
where $\mathbf{W}_{i}$ is an $M\times{M}$ weight matrix, with
$\overline{i}=2$ for $i=1$ and $\overline{i}=1$ for $i=2$. From
\eqref{s3}, the MSE matrix of the signal waveform estimation
denoted by
$\mathbf{MSE}_{\overline{i}}=\mathrm{E}[(\widehat{s}_{i}-{s}_{i})(\widehat{s}_{i}-{s}_{i})^{H}]$,
which can be further written as
\begin{eqnarray}\label{mse}
\mathbf{MSE}_{i}&=&(\mathbf{W}^{H}_{i}\widetilde{\mathbf{H}}_{i}-\mathrm{I}_{M})(\mathbf{W}^{H}_{i}\widetilde{\mathbf{H}}_{i}-\mathrm{I}_{M})^{H}\notag\\
&&+\mathbf{W}^{H}_{i}\mathbf{C}_{n_{i}}\mathbf{W}_{i}
%\mathbf{MSE}_{2}&=&(\mathbf{W}^{H}_{2}\widetilde{\mathbf{H}}_{2}-\mathrm{I}_{N})(\mathbf{W}^{H}_{2}\widetilde{\mathbf{H}}_{2}-\mathrm{I}_{N})^{H}\notag\\
%&&+\mathbf{W}^{H}_{2}\mathbf{C}_{n_{2}}\mathbf{W}_{2},
\end{eqnarray}
where
$\mathbf{C}_{n_{i}}=\textbf{H}^{T}_{i,i}\textbf{F}_{i}\textbf{F}^{H}_{i}\textbf{H}^{*}_{i,i}+\textbf{H}^{T}_{\overline{i},i}\textbf{F}_{\overline{i}}\textbf{F}^{H}_{\overline{i}}\textbf{H}^{*}_{\overline{i},i}+\mathbf{I}_{M}$
%and
%$\mathbf{C}_{n_{2}}=\textbf{H}^{T}_{1,2}\textbf{F}_{1}\textbf{F}^{H}_{1}\textbf{H}^{*}_{1,2}+\textbf{H}^{T}_{2,2}\textbf{F}_{2}\textbf{F}^{H}_{2}\textbf{H}^{*}_{2,2}+\mathbf{I}_{N}$
is the equivalent noise covariance. The sum-MSE of the two source
nodes in the proposed system model can be written as:
\begin{eqnarray}\label{s4}
\mathbf{MSE}_{sum}=\mathbf{MSE}_{1}+\mathbf{MSE}_{2}.
\end{eqnarray}
Note that the sum-MSE minimization criterion measures the overall
transmission performance of both the DL and the UL. Since the two
data streams are transmitted at different directions during the
two time slots are considered in the TWR network.


\section{Joint Source and relay Beamforming Design}

In this section, we develop an iterative QL-QR algorithm by using
the MSE balancing result. The QL-QR algorithm involves two steps,
i.e., the linear receiver matrix optimization and the joint source
and relay beamformer design.

%\subsection{Overview GSVD-BD Precoding Algorithm}
%The design of GSVD-BD precoding algorithm is performed in two
%steps. The received signals at source $i$ is given as
%\begin{eqnarray}\label{overview+}
%y^{G}_{i}=\mathbf{H}^{T}_{i}\mathbf{F}^{G}_{R}\mathbf{H}_{\overline{i}}x_{\overline{i}}+\widetilde{n}^{G}_{i},
%\end{eqnarray}
%where
%$\widetilde{n}^{G}_{i}=\mathbf{H}^{T}_{i}\mathbf{F}^{G}_{R}n_{R}+n_{i}$
%with $n_{R}$ and $n_{i}$ are noise vectors at relay node and
%source node $i$. The GSVD of channel $\mathbf{H}$ can be surmised
%as follow
%\begin{eqnarray}\label{overview0}
%[\mathbf{H}^{G}_{1},
%\mathbf{H}^{G}_{2}]=\mathbf{K}^{G}[\mathbf{D}^{G}_{1}\mathbf{V}^{G}_{1},
%\mathbf{D}^{G}_{2}\mathbf{V}^{G}_{2}],
%\end{eqnarray}
%where $\mathbf{K}^{G}$ and $V^{G}_{i}$ are unitary matrices,
%$\mathbf{D}^{G}_{i}$ are eigenvalues of channel
%$\mathbf{H}^{G}_{i}$ with
%$(\mathbf{D}^{G})^{2}_{1}+(\mathbf{D}^{G})^{2}_{2}=I_{N}$. The
%first precoding filter is used to completely the MU interference
%with noise, then approximate parallel SU-MIMO channels are
%obtained where GSVD is applied at MAC phase. The source precoders
%design are surmised as
%\begin{eqnarray}\label{overview1}
%\mathbf{B}^{G}_{1}=\mathbf{V}^{G}_{1}\mathbf{\Lambda}^{G}_{1}
%~~\mathrm{and}~~
%\mathbf{B}^{G}_{2}=\mathbf{V}^{G}_{2}\mathbf{\Lambda}^{G}_{2},
%\end{eqnarray}
%where $\Lambda^{G}_{i}$ is diagnoal matrix. The second precoding
%filter is used to constraint relay power and implement block
%diagonal signals at source which can be written as
%\begin{eqnarray}\label{overview2}
%\mathbf{F}^{G}_{R}=(\mathbf{K}^{T})^{+}(\mathbf{\Pi}\mathbf{\Lambda}_{R})\mathbf{K}^{+},
%\end{eqnarray}
%where $(.)^{+}$ denotes pseudo-inverse, $\mathbf{\Pi}$ is
%permutation matrix and $\mathbf{\Lambda}_{R}$ is diagonal matrix.
%Applying \eqref{overview0}, \eqref{overview1} and
%\eqref{overview2} to \eqref{overview+}, we have
%\begin{eqnarray}\label{overview3}
%y_{i}=\mathbf{V}^{*}_{i}\mathbf{D}^{T}_{i}\mathbf{\Pi}\mathbf{\Lambda}_{R}\mathbf{D}_{\overline{i}}\mathbf{\Lambda}_{\overline{i}}x_{\overline{i}}+\widetilde{n}^{G}_{i}.
%\end{eqnarray}
%It can be verified that both $\mathbf{B}^{G}_{1}$ and
%$\mathbf{B}^{G}_{2}$ have at most one non-zero entry per row and
%per column. Thus orthogonal transmission can be realized by
%rotating $y_{i}$ using $\mathbf{V}^{G}_{1}$ and
%$\mathbf{V}^{G}_{2}$ , respectively. This orthogonality is the
%main feature of GSVD and it allows performance optimization using
%power allocation matrixes $\mathbf{\Lambda}^{G}_{1}$,
%$\mathbf{\Lambda}^{G}_{2}$ and $\mathbf{\Lambda}^{G}_{R}$ by a
%low-cost iterative water-filling algorithm after $\mathbf{\Pi}$ is
%fixed.

\begin{figure}
\begin{center}
\includegraphics [width=115mm,height=55mm]{DM.eps}
\caption{\label{2} Examples of the SNR regions achieved in a TWRC
with two relays.}
\end{center}
\end{figure}

\subsection{Proposed Optimal Detector and Optimization Problem}
We would like to find the jointly optimal beamforming vectors
$\mathbf{W}_{i}$, $\mathbf{V}_{i}$ and $\mathbf{F}_{i}$ such as
the following sum-MSE is minimized
\begin{eqnarray}\label{MSEc}
\min\limits_{\mathbf{W}_{1},\mathbf{W}_{2},\mathbf{F}_{1},\mathbf{F}_{2},\mathbf{V}_{1},\mathbf{V}_{2}}~~\mathbf{MSE}_{sum}.
\end{eqnarray}
According to \eqref{R power}, we consider the following
transmission power constraint at relay node
\begin{eqnarray}\label{total R}
tr{\left(\textbf{F}_{1}\textbf{D}_{1}\textbf{F}^{H}_{1}+\textbf{F}_{2}\textbf{D}_{2}\textbf{F}^{H}_{2}\right)}\leq{P_{R_{1}}+P_{R_{2}}=P_{R}},
\end{eqnarray}
where
$\textbf{D}_{1}=\rho^{2}_{R}(\textbf{H}_{1,1}\mathbf{V}_{1}\mathbf{V}^{H}_{1}\textbf{H}^{H}_{1,1}+\textbf{H}_{1,2}\mathbf{V}_{2}\mathbf{V}^{H}_{2}\textbf{H}^{H}_{1,2}+\mathbf{I}_{N})$
and
$\textbf{D}_{2}=\rho^{2}_{R}(\textbf{H}_{2,1}\mathbf{V}_{1}\mathbf{V}^{H}_{1}\textbf{H}^{H}_{2,1}+\textbf{H}_{2,2}\mathbf{V}_{2}\mathbf{V}^{H}_{2}\textbf{H}^{H}_{2,2}+\mathbf{I}_{N})$.
The $P_{R_{i}}$ denotes the power constraint at the relay node
$\mathbf{R}_{i}$, and $P_{R}$ is the total relay power. The
transmission power constraint at two source nodes can be written
as
\begin{eqnarray}\label{source power}
tr{(\mathbf{V}_{i}\mathbf{V}^{H}_{i})}\leq{P_{i}},~~~~~ i=1,2
\end{eqnarray}
where $P_{i}$ is the available power at the $i$th source node.
According to \eqref{MSEc}, \eqref{total R} and \eqref{source
power}, the joint optimization problem of the sum-MSE can be
formulated as follows:
\begin{eqnarray}\label{MSEduality}
&\min\limits_{\mathbf{W}_{1},\mathbf{W}_{2},\mathbf{F}_{1},\mathbf{F}_{2},\mathbf{V}_{1},\mathbf{V}_{2}}&~~~~~~~~~~~~\mathbf{MSE}_{sum}\notag\\
&s.t.&~tr{\left(\textbf{F}_{1}\textbf{D}_{1}\textbf{F}^{H}_{1}+\textbf{F}_{2}\textbf{D}_{2}\textbf{F}^{H}_{2}\right)}\leq{P_{R}}\notag\\
&s.t.&~~~~~~~~tr{(\mathbf{V}_{i}\mathbf{V}^{H}_{i})}\leq{P_{i}}.
\end{eqnarray}
It is shown in \cite{b17} that at the optimum, $SNR_{1}=SNR_{2}$
holds true, thus leading to an SNR balancing result. Otherwise, if
$SNR_{1}>SNR_{2}$, then $P_{2}$ can be reduced to retain
$SNR_{1}=SNR_{2}$, and this reduction of $P_{2}$ will not violate
the power constraint, i.e.,
\begin{eqnarray}\label{MSEbalance}
P_{1}\cdot SNR_{1}=P_{2}\cdot SNR_{2}.
\end{eqnarray}
In Fig. \ref{2}, we show two examples of the SNR regions with
$\alpha_{1}=0.5$ and $\alpha_{2}=0.3$, where $\omega_{i}\in[0,1]$
is a Lagrange multiplier weight value and $\alpha_{i}\in[0,1]$ is
an SNR weight value. We have assumed the sum of SNR is a constant
value. It is clear that the SNR region of $\alpha_{1}$ is larger
than that of $\alpha_{2}$. For further details, see \cite{b17}.
%To simplify  optimization problem, we obtain the necessary and
%suifficient conditions of $\mathbf{MSE}_{1}=\mathbf{MSE}_{2}$ or
%$\varepsilon_{1}=\varepsilon_{2}$. It is intuitively clear that
%source precoders and the relay beamformers lie in the subspace
%spanned by the channel vectors \cite{a4}, we have
%\begin{eqnarray}\label{total R power}
%\mathbf{F}_{1}&=&\left(\mathbf{H}^{T}_{4}-\mathbf{H}_{2}^{T}\right)\left(\mathbf{H}^{T}_{1}-\mathbf{H}_{3}^{T}\right)^{\dag}\mathbf{F}_{2}\notag\\
%&=&\aleph \mathbf{F}_{2}
%\end{eqnarray}
%and
%\begin{eqnarray}\label{total R power1}
%\mathbf{V}_{1}\!\!\!\!\!&=&\!\!\!\!\!\!\left(\mathbf{H}^{T}_{1}\mathbf{F}_{1}\mathbf{H}_{2}\!+\!\mathbf{H}^{T}_{3}\mathbf{F}_{2}\mathbf{H}_{4}\right)\!
%\left(\mathbf{H}^{T}_{2}\mathbf{F}_{1}\mathbf{H}_{1}\!+\!\mathbf{H}^{T}_{4}\mathbf{F}_{2}\mathbf{H}_{3}\right)\!\!^{\dag}\mathbf{V}\!_{2}\notag\\
%&=&\!\!\!\!\!\Theta \mathbf{V}_{2},
%\end{eqnarray}
%where $(.)^{\dag}$ stands for pseudo-inverse.
As discussed in \cite{a5}, the optimization problems have the
performance matrix that are functions of SNR , namely the MSE at
the output of a linear-MMSE (LMMSE) filter of each user
\begin{eqnarray}\label{MSEbalance}
\mathbf{MSE}=\frac{1}{1+SNR}.
\end{eqnarray}
By these two approaches, the max-min optimization problem in
\eqref{MSEduality} can be efficiently written as
\begin{eqnarray}\label{MSEequal}
&\min\limits_{\mathbf{W}_{1},\mathbf{F}_{i},\mathbf{V}_{2}}&~~~~~~~~~~~~\mathbf{MSE}_{1}\\
&s.t.&~tr{(\textbf{F}_{1}\textbf{D}_{1}\textbf{F}^{H}_{1}+\textbf{F}_{2}\textbf{D}_{2}\textbf{F}^{H}_{2}})\leq{P_{R}}\label{MSE22}\\
&s.t.&~~~~~~\mathbf{MSE}_{1}=\mathbf{MSE}_{2}\label{MSE33},
\end{eqnarray}
where $i\in{1,2}$. Since the optimization problem \eqref{MSEequal}
is nonconvex, it is difficult to obtain the globally optimal
solution. In this paper, we present a locally optimal solution of
the joint optimization problem over $\mathbf{W}_{i}$,
$\mathbf{V}_{i}$ and $\mathbf{F}_{i}$ where $i=1,2$, which can be
solved by three stages, i.e., $1:$ The linear receiver weighted
matrices are optimized with the fixed source precoding matrix
$\mathbf{V}_{i}$ and relay amplifying matrices $\mathbf{F}_{i}$
($\mathbf{W}_{i}$ is not in constraints \eqref{MSE22} and
\eqref{MSE33}). $2:$ With given $\mathbf{W}_{i}$ and fixed
$\mathbf{F}_{i}$, update $\mathbf{V}_{i}$. $3:$ With given
$\mathbf{W}_{i}$ and $\mathbf{V}_{i}$, obtain suboptimal
$\mathbf{F}_{i}$ to solve \eqref{MSEequal}.

\emph{Lemma 1 }: For any fixed $\mathbf{V}_{i}$ and
$\mathbf{F}_{i}$, the minimization problems in \eqref{MSEequal}
are convex quadratic problems and the optimal $\mathbf{W}_{i}$ can
be obtained as the Wiener filter which is used to decode $s_{i}$
shown as follows
\begin{eqnarray}\label{MMSE}
\mathbf{W}^{o}_{i}=(\mathbf{\widetilde{H}}_{i}\mathbf{\widetilde{H}}^{H}_{i}+\mathbf{C}_{n_{i}})^{-1}\mathbf{\widetilde{H}}_{i},
\end{eqnarray}
\emph{Proof}: For source node $i$, the MSE can be further
expressed as:
\begin{eqnarray}\label{msep}
\mathbf{MSE}_{i}&=&\mathbf{W}^{H}_{i}\widetilde{\mathbf{H}}_{i}\widetilde{\mathbf{H}}^{H}_{i}\mathbf{W}_{i}-\mathbf{W}^{H}_{i}\widetilde{\mathbf{H}}_{i}-\widetilde{\mathbf{H}}^{H}_{i}\mathbf{W}_{i}\notag\\
&&+\mathbf{I}_{M}+\mathbf{W}^{H}_{i}\mathbf{C}_{n_{i}}\mathbf{W}_{i}
\end{eqnarray}
Based on \eqref{msep}, the derivation of an optimal MSE detection
matrix $\mathbf{W}^{opt}_{i}$ is equivalent to solving the
following equation:
\begin{eqnarray}\label{msep1}
\frac{\partial\mathbf{MSE}_{i}}{\partial\mathbf{W}^{H}_{i}}=2\widetilde{\mathbf{H}}_{i}\widetilde{\mathbf{H}}^{H}_{i}\mathbf{W}_{i}-2\widetilde{\mathbf{H}}_{i}+2\mathbf{C}_{n_{i}}\mathbf{W}_{i}=0.
\end{eqnarray}
Then, we may obtain closed-form solution of $\mathbf{W}_{i}$,
which is
\begin{eqnarray}\label{msep1}
\mathbf{W}^{o}_{i}=(\mathbf{\widetilde{H}}_{i}\mathbf{\widetilde{H}}^{H}_{i}+\mathbf{C}_{n_{i}})^{-1}\mathbf{\widetilde{H}}_{i}.
\end{eqnarray}
This completes the proof.

%\begin{eqnarray}\label{s6}
%\varepsilon_{1}&=&\textrm{arg}~min~\mathrm{E}\{\|x_{\overline{i}}-\mathbf{W}^{H}_{i}\widetilde{y}_{i}\|^{2}\}\notag\\
%&=&\textrm{tr}\{(\mathbf{W}^{H}_{i}\widetilde{y}_{i}-x_{\overline{i}})\times(\mathbf{W}^{H}_{i}\widetilde{y}_{i}-x_{\overline{i}})^{H}\}\notag\\
%&=&P_{2}\cdot\textrm{tr}\{\mathbf{W}^{H}_{1}\Phi_{1}\mathbf{W}_{1}+\mathbf{W}^{H}_{1}\Phi_{2}\mathbf{W}_{1}+\mathbf{W}^{H}_{1}\Phi_{3}\mathbf{W}_{1}\notag\\
%&&+\mathbf{W}^{H}_{1}\Phi_{4}\mathbf{W}_{1}+\mathbf{W}^{H}_{1}\Phi_{5}\mathbf{W}_{1}-\mathbf{W}^{H}_{1}\Omega_{1}-\mathbf{W}^{H}_{1}\Omega_{2}\notag\\
%%&&-\Omega_{3}\mathbf{W}_{1}-\Omega_{4}\mathbf{W}_{1}-I_{M}\}
%\end{eqnarray}
%where $\Phi_{1}=\mathbf{A}\mathbf{A}^{H}$,
%$\Phi_{2}=\mathbf{A}\mathbf{B}^{H}$,
%$\Phi_{3}=\mathbf{B}\mathbf{A}^{H}$,
%$\Phi_{4}=\mathbf{B}\mathbf{A}^{H}$ and
%$\Phi_{5}=\mathbf{H}^{T}_{1}{\mathbf{F}}_{1}{\mathbf{F}}^{H}_{1}\mathbf{H}^{*}_{1}\sigma^{2}_{R_{1}}
%+\mathbf{H}^{T}_{3}{\mathbf{F}}_{2}{\mathbf{F}}^{H}_{2}\mathbf{H}^{*}_{3}\sigma^{2}_{R_{2}}+\sigma^{2}_{i}$
%with
%$\mathbf{A}=\textbf{H}^{T}_{1}\mathbf{F}_{1}\textbf{H}_{2}\mathbf{V}_{2}$
%and
%$\mathbf{B}=\textbf{H}^{T}_{3}\mathbf{F}_{2}\textbf{H}_{4}\mathbf{V}_{2}$,
%respectively. And
%$\Omega_{1}=\mathbf{H}^{T}_{1}\mathbf{F}_{1}\mathbf{H}_{2}\mathbf{V}_{2}$,
%$\Omega_{2}=\mathbf{H}^{T}_{3}\mathbf{F}_{2}\mathbf{H}_{4}\mathbf{V}_{2}$,
%$\Omega_{3}=\mathbf{V}^{H}_{2}\mathbf{H}^{H}_{2}\mathbf{F}_{1}\mathbf{H}^{*}_{1}$
%and
%$\Omega_{4}=\mathbf{V}^{H}_{2}\mathbf{H}^{H}_{4}\mathbf{F}^{H}_{2}\mathbf{H}^{*}_{3}$.
%We have assumed that the signals and noise are independent with
%each other. Based on \eqref{s6}, the derivation of optimal MSE
%detection matrix $W^{opt}_{i}$ is equivalent to solving the
%following equation:
%\begin{eqnarray}\label{derivative}
%\frac{\sigma{\varepsilon_{1}}}{\sigma{W_{1}}}&=&W_{1}(\Phi_{1}+\Phi_{2}+\Phi_{3}+\Phi_{4}+\Phi_{5})\notag\\
%&&-\Omega_{1}-\Omega_{2}-\Omega_{3}-\Omega_{4}=0.
%\end{eqnarray}
%Obviously, the closed-form expression of $W^{opt}_{i}$ can be
%obtained, which can be expressed as
%\begin{eqnarray}\label{optimal W}
%W^{opt}_{1}&=&(\Omega_{1}+\Omega_{2}+\Omega_{3}+\Omega_{4})\times\notag\\
%&&(\Phi_{1}+\Phi_{2}+\Phi_{3}+\Phi_{4}+\Phi_{5})^{-1}.
%\end{eqnarray}
%It is simply that
%$\Phi_{1}+\Phi_{2}+\Phi_{3}+\Phi_{4}=\widetilde{H}_{i}\widetilde{H}^{H}_{i}$,
%$\Phi_{5}=\Phi_{i}$ and
%$\Omega_{1}+\Omega_{2}+\Omega_{3}+\Omega_{4}=\widetilde{H}_{i}$
%which are shown in \emph{lemma} 1, similarly $\varepsilon_{2}$ can
%be calculated as same way, this completes the proof.

With the optimal $\mathbf{W}^{o}_{1}$ fixed, the outer
minimization problem in \eqref{MSEequal} can be rewritten as
\begin{eqnarray}\label{MSEopt}
&\min\limits_{\mathbf{F}_{1},\mathbf{F}_{2},\mathbf{V}_{1},\mathbf{V}_{2}}&~~~~~~~~~~~~\mathbf{MSE}^{o}_{1}\notag\\
&s.t.&~tr{(\textbf{F}_{1}\textbf{D}_{1}\textbf{F}^{H}_{1}+\textbf{F}_{2}\textbf{D}_{2}\textbf{F}^{H}_{2})}\leq {P_{R}}\notag\\
%&s.t.&~~~tr{(\textbf{F}_{2}\textbf{D}_{2}\textbf{F}^{H}_{2})}\leq (1-a){P_{R}}\notag\\
&s.t.&~~~~~~~\mathbf{MSE}_{1}=\mathbf{MSE}_{2},
\end{eqnarray}
where $\mathbf{MSE}^{o}_{1}$ is the MSE matrix using
$\mathbf{W}^{o}_{1}$. By substituting \eqref{MMSE} into
\eqref{mse}, we have
\begin{eqnarray}\label{MSEequ}
\mathbf{MSE}^{o}_{1}=[\mathbf{I}_{M}+\mathbf{\widetilde{H}}^{H}_{1}\mathbf{C}^{-1}_{n_{1}}\mathbf{\widetilde{H}}_{1}]^{-1}.~~~
i=1,2.
\end{eqnarray}
Note that the matrix inversion lemma is used to obtain
\eqref{MSEequ}.
%According to the \emph{Theorem 1} in [14], the optimal duality MSE
%solution is given as
%\begin{eqnarray}\label{Optimal.duality}
%MSE_{S_{1}}=MSE_{S_{2}}.
%\end{eqnarray}
%The problem \eqref{MSEduality} is highly nonconvex and a
%closed-form expression of the optimal $\textbf{F}_{i}$ and
%$\textbf{V}_{i}$ is intractable. From \eqref{R power}, the total
%relay power can be expressed as
%\begin{eqnarray}\label{total R power1}
%P_{R_{1}}+P_{R_{2}}=\textbf{F}^{H}_{1}\textbf{D}_{1}\textbf{F}_{1}+\textbf{F}^{H}_{2}\textbf{D}_{2}\textbf{F}_{2},
%\end{eqnarray}
%where
%$\textbf{D}_{1}=\rho^{2}_{R}(P_{1}\textbf{H}_{1}\textbf{H}^{H}_{1}+P_{2}\textbf{H}_{2}\textbf{H}^{H}_{2}+I_{M})$
%and
%$\textbf{D}_{2}=\rho^{2}_{R}(P_{1}\textbf{H}_{3}\textbf{H}^{H}_{3}+P_{2}\textbf{H}_{4}\textbf{H}^{H}_{4}+I_{M})$.
%Therefore, based on \emph{Lemma 1}, \eqref{Optimal.duality} and
%\eqref{total R power}, the MSE duality may be optimized by solving
%the following problem
%\begin{eqnarray}\label{total R power2}
%&min&~~~~~~~~~2MSE_{S_{1}}\notag\\
%&s.t.&~\textbf{F}^{H}_{1}\textbf{D}_{1}\textbf{F}_{1}+\textbf{F}^{H}_{2}\textbf{D}_{2}\textbf{F}_{2}\leq{P_{R}}\notag\\
%&s.t.&~~~MSE_{S_{1}}=MSE_{S_{2}},
%\end{eqnarray}

\subsection{Joint Optimal Source and Relay Beamforming Matrices Design and Iterative Algorithm}


In this section, we focus on the source and relay beamforming
matrices design and develop an iterative algorithm which is
suboptimal for the general case, but has a much lower
computational complexity. For the fixed $\mathbf{F}_{i}$, the
source precoding matrix $\mathbf{V}_{i}$ is optimized by solving
the following problem
\begin{eqnarray}\label{MSEV}
&\min\limits_{\mathbf{V}_{1},\mathbf{V}_{2}}&~~~~~~~~~~~~tr\left[\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}\right]^{-1}\notag\\
&s.t.&~~~tr\{\rho^{2}_{R}(\textbf{V}^{H}_{1}\mathbf{\Psi}_{1}\textbf{V}_{1}+\textbf{V}^{H}_{2}\mathbf{\Psi}_{2}\textbf{V}_{2})\}\leq{P_{R}}\notag\\
 &s.t.&~~~~~~~~~~~~~~~~tr\{\mathbf{V}^{H}_{i}\mathbf{V}_{i}\}\leq P_{i},
\end{eqnarray}
where
$\mathbf{\Phi}=\mathbf{\widehat{H}}^{H}_{1}\mathbf{\widehat{C}}^{-1}_{1}\mathbf{\widehat{H}}_{1}$
,
$\mathbf{\Psi}_{1}=\mathbf{H}^{H}_{1,1}\mathbf{F}^{H}_{1}\mathbf{F}_{1}\mathbf{H}_{1,1}+\mathbf{H}^{H}_{2,1}\mathbf{F}^{H}_{2}\mathbf{F}_{2}\mathbf{H}_{2,1}$,
and
$\mathbf{\Psi}_{2}=\mathbf{H}^{H}_{1,2}\mathbf{F}^{H}_{1}\mathbf{F}_{1}\mathbf{H}_{1,2}+\mathbf{H}^{H}_{2,2}\mathbf{F}^{H}_{2}\mathbf{F}_{2}\mathbf{H}_{2,2}$.
The Lagrangian function associated with the problem \eqref{MSEV}
is given by
\begin{eqnarray}\label{MSEL}
L_{\mathbf{V}}&=&tr\left[\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}\right]^{-1}+\sum^{2}_{i=1}\mu_{i}\left(tr\{\mathbf{V}^{H}_{i}\mathbf{V}_{i}\}-P_{i}\right)\notag\\
&&+\mu_{3}\left\{tr\{\rho^{2}_{R}(\textbf{V}^{H}_{1}\mathbf{\Psi}_{1}\textbf{V}_{1}+\textbf{V}^{H}_{2}\mathbf{\Psi}_{2}\textbf{V}_{2})\}-P_{R}\right\},
\end{eqnarray}
where $\mu_{i}\geq0$ is the Lagrange multiplier.

\emph{Case 1}: When $\mu_{i}=0$, making the derivative of
$L_{\mathbf{V}}$ with respect to $\mathbf{V}_{2}$ be zero, we
obtain
\begin{eqnarray}\label{MSEL2}
\frac{\partial L_{\mathbf{V}}}{\partial
\mathbf{V}_{2}}=-\left[\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}\right]^{-2}\mathbf{V}^{H}_{2}\mathbf{\Phi}=0.
\end{eqnarray}
Since $\mathbf{V}_{2}$ and $\mathbf{\Phi}$ are nonsingular
matrices, \eqref{MSEL2} can be represented as
\begin{eqnarray}\label{MSEL3}
\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}=0.
\end{eqnarray}
Simplifying \eqref{MSEL3}, $\mathbf{I}_{M}>0$ and
$\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}\geq0$.
Consequently, in \emph{Case 1}, the optimal solution is not
existent.

\emph{Case 2}: When $\mu_{i}>0$, we rewrite the Lagrangian
function as
\begin{eqnarray}\label{MSEL4}
L_{\mathbf{V}}&=&\!\!\!\!\left[\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}\right]^{-1}-\mu_{1}P_{1}-\mu_{2}P_{2}-\mu_{3}P_{R}\notag\\
&&\!\!\!\!+\mathbf{V}^{H}_{2}\Upsilon\Upsilon^{H}\mathbf{V}_{2}\!+\!\mu_{3}\rho^{2}_{R}\mathbf{V}_{1}^{H}\mathbf{\Psi}_{1}\mathbf{V}_{1}\!+\!\mu_{1}\mathbf{V}_{1}^{H}\mathbf{V}_{1}\!,
\end{eqnarray}
where
$\Upsilon\Upsilon^{H}=\mu_{2}\mathbf{I}_{M}+\mu_{3}\rho^{2}_{R}\mathbf{\Psi}_{2}$.
We obtain the derivative of $L_{\mathbf{V}}$ as
\begin{eqnarray}\label{MSEL5}
\frac{\partial L_{\mathbf{V}}}{\partial
\mathbf{V}_{2}}=-\left[\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}\right]^{-2}\mathbf{V}^{H}_{2}\mathbf{\Phi}+\mathbf{V}^{H}_{2}\Upsilon\Upsilon^{H}=0.
\end{eqnarray}
Since $\mathbf{V}^{H}_{2}$ and $\mathbf{\Phi}$ are nonsingular
matrices, multiply both sides by
$\left(\mathbf{V}^{H}_{2}\right)^{-1}$ and $\mathbf{\Phi}^{-1}$,
we have
\begin{eqnarray}\label{MSEL6}
\left(\mathbf{V}^{H}_{2}\right)^{-1}\left[\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\mathbf{\Phi}\mathbf{V}_{2}\right]^{-2}\mathbf{V}^{H}_{2}=\Upsilon\Upsilon^{H}\mathbf{\Phi}^{-1}.
\end{eqnarray}
Due to $\mathbf{\Phi}$ is Hermitian and positive definite, we
apply the Choleskey decomposition of
$\mathbf{\Phi}=\Omega^{H}\Omega$, where $\Omega$ is a lower
triangular matrix. Consequently, we represent \eqref{MSEL6} as
\begin{eqnarray}\label{MSEL7}
&&\left(\Omega^{H}\right)^{-1}\left(\mathbf{V}^{H}_{2}\right)^{-1}\left[\mathbf{I}_{M}+\mathbf{V}^{H}_{2}\Omega^{H}\Omega\mathbf{V}_{2}\right]^{-2}\mathbf{V}^{H}_{2}\Omega^{H}\notag\\
&=&\left(\Omega^{H}\right)^{-1}\Upsilon\Upsilon^{H}\left(\Omega^{H}\Omega\right)^{-1}\Omega^{H}.
\end{eqnarray}
By the definition of the matrix identity as
\begin{eqnarray}\label{ID}
\left[\mathbf{I}_{M}+\mathbf{X}\mathbf{X}^{H}\right]^{-1}\mathbf{X}=\mathbf{X}\left[\mathbf{I}_{N}+\mathbf{X}^{H}\mathbf{X}\right]^{-1},
\end{eqnarray}
for any $M\times N$ matrix $\mathbf{X}$, we can rewrite
\eqref{MSEL7} as
\begin{eqnarray}\label{MSEL8}
\left[\mathbf{I}_{M}+\Omega\mathbf{V}_{2}\mathbf{V}^{H}_{2}\Omega^{H}\right]^{-2}=\left(\Omega^{H}\right)^{-1}\Upsilon\Upsilon^{H}\Omega^{-1}.
\end{eqnarray}
Solving \eqref{MSEL8} for $\mathbf{V}_{2}$, we obtain
\eqref{MSEL7} as
\begin{eqnarray}\label{MSEL8}
\mathbf{V}_{2}=\left(\nabla\nabla^{H}-\mathbf{\Phi}^{-1}\right)^{\frac{1}{2}},
\end{eqnarray}
where $\nabla=\left(\Upsilon^{H}\Omega\right)^{-\frac{1}{2}}$.
Obviously, the precoding matrix ${\textbf{V}}_1$ can be obtained
in the same way.

%Using decompositions \eqref{LQ} and \eqref{QR}, the received
%signals at receivers can be written as,
%\begin{eqnarray}\label{QR+SVD}
%y_{R_{\overline{i}}}&=&(\mathbf{H}^{T}_{R_{i},\overline{i}}\mathbf{F}_{i}\mathbf{H}_{R_{i},i}+\mathbf{H}^{T}_{R_{\overline{i}},\overline{i}}\mathbf{F}_{i}\mathbf{H}_{R_{\overline{i}},i})s_{\overline{i}}+\notag\\
%&&\mathbf{H}^{T}_{R_{i},\overline{i}}\mathbf{F}_{i}n_{R_{i}}+\mathbf{H}^{T}_{R_{\overline{i}},\overline{i}}n_{R_{\overline{i}}}+n_{\overline{i}}.
%\end{eqnarray}
\begin{figure}
\begin{center}
\includegraphics [width=70mm,height=30mm]{RF.eps}
\caption{\label{3} The relay filter design of the proposed QL-QR
technique.}
\end{center}
\end{figure}

Figure 3 shows our proposed relay filter design, which forwards
the received signal (input) from $\mathbf{S}_{1}$ amplified by a
Left Filter (LF) matrix $\mathbf{F}_{L,i}$ and the signal from
$\mathbf{S}_{2}$ amplified by a Right Filter (RF) matrix
$\mathbf{F}_{R,i}$ to the Center Filter (CF) $\mathbf{F}_{D,i}$
that amplifies the outputs from the Left Filter (LF) matrix
$\mathbf{F}_{L,i}$ and the Right Filter (RF) matrix
$\mathbf{F}_{R,i}$, and forward them to $\mathbf{S}_{1}$ and
$\mathbf{S}_{2}$ (output).\footnote{For example: For
$\mathbf{S}_{1}$, the equivalent channel can be written as
$\widetilde{\mathbf{H}}_{1}=\textbf{H}^{T}_{1,1}\textbf{F}_{1}\textbf{H}_{1,2}+\textbf{H}^{T}_{2,1}\textbf{F}_{2}\textbf{H}_{2,2}=\textbf{H}^{T}_{1,1}\textbf{F}_{L,1}\textbf{F}_{D,1}\textbf{F}_{R,1}\textbf{H}_{1,2}+\textbf{H}^{T}_{2,1}\textbf{F}_{L,2}\textbf{F}_{D,2}\textbf{F}_{R,2}\textbf{H}_{2,2}$.
For $\mathbf{S}_{2}$, the equivalent channel can be written as
$\widetilde{\mathbf{H}}_{2}=\textbf{H}^{T}_{1,2}\textbf{F}_{1}\textbf{H}_{1,1}+\textbf{H}^{T}_{2,2}\textbf{F}_{2}\textbf{H}_{2,1}=\textbf{H}^{T}_{1,2}\textbf{F}_{R,1}\textbf{F}_{D,1}\textbf{F}_{L,1}\textbf{H}_{1,1}+\textbf{H}^{T}_{2,1}\textbf{F}_{R,2}\textbf{F}_{D,2}\textbf{F}_{L,2}\textbf{H}_{2,2}$.}


\emph{Lemma 2}: The optimal relay filter constructive of
$\mathbf{F}_{L}$, $\mathbf{F}_{R}$ and $\mathbf{F}_{D}$ matrices
i.e., for $\mathbf{R}_{1}$ and $\mathbf{R}_{2}$ can be designed as
\begin{eqnarray}\label{RL+RF}
&&\mathbf{F}_{L,1}=\mathbf{Q}^{*}_{L,1},~\mathbf{F}_{L,2}=\mathbf{Q}^{*}_{L,2},\notag\\
&&\mathbf{F}_{R,1}=\mathbf{Q}^{H}_{R,1},~\mathbf{F}_{R,2}=\mathbf{Q}^{H}_{R,2},\notag\\
&&\mathbf{F}_{D,1}~\mathrm{and}~ \mathbf{F}_{D,2}~\mathrm{ are~
diagonal~ matrix}.
\end{eqnarray}
\emph{Proof}: For $\mathbf{F}_{L,i}$ and $\mathbf{F}_{R,i}$, the
proof is similar as \emph{Theorem 3.1} in \cite{a6}. For
$\mathbf{F}_{D,i}$, using \emph{Theorem 2} in \cite{b8}, the
structure of $\mathbf{F}_{D,i}$ is optimal for the two cases of
\begin{eqnarray}\label{OP1}
&&(a):~~~~~R(\mathbf{H}_{1}\mathbf{V}_{1})\bot R(\mathbf{H}_{2}\mathbf{V}_{2}); (\mathbf{H}^{*}_{1})\bot R(\mathbf{H}^{*}_{2})\notag\\
&& ~~\mathrm{and}~\, ~~~R(\mathbf{H}_{3}\mathbf{V}_{1})\bot R(\mathbf{H}_{4}\mathbf{V}_{2}); (\mathbf{H}^{*}_{3})\bot R(\mathbf{H}^{*}_{4})\notag\\
 &&(b):~~~~~R(\mathbf{H}_{1}\mathbf{V}_{1})\| R(\mathbf{H}_{2}\mathbf{V}_{2}); (\mathbf{H}^{*}_{1})\bot R(\mathbf{H}^{*}_{2})\notag\\
&& ~~\mathrm{and}~\, ~~~R(\mathbf{H}_{1}\mathbf{V}_{1})\|
R(\mathbf{H}_{2}\mathbf{V}_{2}); (\mathbf{H}^{*}_{1})\bot
R(\mathbf{H}^{*}_{2}).
\end{eqnarray}
\emph{Case a}: If $N=2$, the optimal $\mathbf{F}_{D,i}$ is a
diagonal matrix given as
\begin{eqnarray}\label{OP4}
\mathbf{F}_{D,i} &\triangleq& {\left[\begin{array}{cc}
              \mathbf{f}_{d,i,1} & \mathbf{0} \\
             \mathbf{0} & \mathbf{f}_{d,i,2} \\
             \end{array}\right]}.
             \end{eqnarray}
If $N=2a$, $a={2,3,...}$, the optimal $\mathbf{F}_{D,i}$ is a
$2\times2$ block diagonal matrix given as
\begin{eqnarray}\label{OP2}
\mathbf{F}_{D,i} &\triangleq& {\left[\begin{array}{cc}
              \mathbf{F}_{D,i,1} & \mathbf{0} \\
             \mathbf{0} & \mathbf{F}_{D,i,2} \\
             \end{array}\right]},
             \end{eqnarray}
where $\mathbf{F}_{D,i,1}$ and $\mathbf{F}_{D,i,2}$ are
$\frac{N}{2}\times \frac{N}{2}$ matrices. \\
\emph{Case b}: The optimal $\mathbf{F}_{D,i}$ is defined as
\begin{eqnarray}\label{OP6}
\mathbf{F}^{\star}_{D,i} &\triangleq& {\left[\begin{array}{cc}
              \mathbf{f}_{d,i,1} \\
             \mathbf{f}_{d,i,2} \\
             \end{array}\right]}.
             \end{eqnarray}
\emph{Discussion}: In \emph{Case b}, since
$\mathbf{F}^{\star}_{D,i}$ is optimal, but the computational
complexity will be considerably increased compared with \emph{Case
a}, so we exclude it.
%Applying the previous per user MSE duality to these virtual users,
%we obtain the per stream MSE duality. The rows of
%$\mathbf{\widetilde{H}}_{i}$ determine the effective transmit
%power amplification. Then, the precoding filter can be obtained
%\begin{equation}\label{Equation_P1}
%{\textbf{V}}^{b}_{i}=(\mathbf{\widetilde{H}}_{i}\mathbf{\widetilde{H}}^{H}_{i})^{-1}\mathbf{\widetilde{H}}_{i}.
%\end{equation}
%Finally, the second precoding matrix ${\textbf{V}}^{b}_i$ for
%total system is
%\begin{equation}\label{P2}
%{\textbf{V}}^{b}_i = {\left[\begin{array}{ccc}
%              {v}_{1}^b & \ldots & 0 \\
%              \vdots & \ddots &  \vdots \\
%              0      &  \ldots & {v}_{n}^b \\
%              \end{array}\right]}_{M\times M}.
%\end{equation}
%The observation at the source node $\mathbf{S}_{1}$ can be given
%as
%\begin{eqnarray}\label{equ1}
%\widehat{y}_{1}&=&\!\!\!\!(\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}\mathbf{R}_{2}+\mathbf{L}^{T}_{3}\mathbf{F}_{D,2}\mathbf{R}_{4})\mathbf{V}_{2}x_{2}\notag\\
%&&\!\!\!\!+\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}\mathbf{Q}^{H}_{R,1}n_{R_{1}}\!+\!\mathbf{L}^{T}_{3}\mathbf{F}_{D,2}\mathbf{Q}^{H}_{R,2}n_{R_{2}}\!+\!n_{2}\notag\\
%&=&\!\!\!\!\mathbf{\widehat{H}}_{1,2}x_{2}+\widehat{n}_{1},
%\end{eqnarray}
%where
%\begin{eqnarray}\label{Hnew}
%\mathbf{\widehat{H}}_{1,2}=(\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}\mathbf{R}_{2}+\mathbf{L}^{T}_{3}\mathbf{F}_{D,2}\mathbf{R}_{4})\mathbf{V}_{2}.
%\end{eqnarray}

For \emph{Case a}, Before we develop a numerical method to solve
vector $\mathbf{F}_{D,i}$, let us have some insights into the
structure of this suboptimal relay beamforming matrix. To simplify
relay beamforming matrix $\mathbf{F}_{D,i}$, we introduce a
following property:

\emph{Property 1}: The statistical behavior of a unitary matrix
$\mathbf{U}$ remains unchanged when multiplied by any unitary
matrix $\mathbf{T}$ independent of $\mathbf{U}$. In other worlds,
$\mathbf{TU}$ has the same distribution as $\mathbf{U}$, i.e., in
\eqref{RL+RF},
\begin{eqnarray}\label{Prop}
&&|\mathbf{F}_{1}|=|\mathbf{F}_{L,1}\mathbf{F}_{D,1}\mathbf{F}_{R,1}|=|\mathbf{F}_{D,1}|\notag\\
&&|\mathbf{F}_{2}|=|\mathbf{F}_{L,2}\mathbf{F}_{D,2}\mathbf{F}_{R,2}|=|\mathbf{F}_{D,2}|.
\end{eqnarray}
Now, let us introduce the
following $\mathbf{QL}$ decompositions%$\mathbf{GQL}$,
%specifically, if $\mathbf{H}_{1}$ and $\mathbf{H}_{3}$ are two
%$N\times M$ matrices and satisfy $N\leq{M}\leq{2N}$, then
%$\mathbf{GQL}$ decomposition can be applied on them as
\begin{eqnarray}\label{LQ}
\left[\mathbf{H}_{1,1}\mathbf{V}_{1},
\mathbf{H}_{2,1}\mathbf{V}_{1}\right]&=&\left[\mathbf{Q}_{L,1}\mathbf{L}_{1},
\mathbf{Q}_{L,2}\mathbf{L}_{2} \right],
%&=&\mathbf{\Pi}_{L}\left[\mathbf{Q}_{1}\mathbf{L}_{1},
%\mathbf{Q}_{3}\mathbf{L}_{3}\right]
\end{eqnarray}
where $\mathbf{Q}_{L,i}$ for $i=1,2$, is a unitary matrix with a
dimension ${\mathbb{C}^{N\times{N}}}$,
%$\mathbf{\Pi}_{L}=\left[\mathbf{Q}_{1},\mathbf{Q}_{3}\right]^{T}\mathbf{Q}_{L}$
%is permutation matrix
and $\left\{\mathbf{L}_{1}, \mathbf{L}_{2}\right\}\in
{\mathbb{C}^{N\times{M}}}$ are lower triangular matrices. %By briefly
%outline the $\mathbf{QL}$ decomposition principle, we have
%\begin{eqnarray}\label{LQ1}
%\mathbf{Q}_{L}\mathbf{Q}^{H}_{L}&=&\left(\left[\mathbf{Q}_{1},\mathbf{Q}_{3}\right]^{*}\mathbf{\Pi}_{L}\right)\left(\left[\mathbf{Q}_{1},\mathbf{Q}_{3}\right]^{*}\mathbf{\Pi}_{L}\right)^{H}\notag\\
%&=&\mathbf{Q}^{2}_{1}+\mathbf{Q}^{2}_{3}\notag\\
%&=&\mathbf{I}_{N}.
%\end{eqnarray}
%$\mathbf{\Sigma}_{1}=[\mathbf{0}^{T}_{(N-M)\times{M}},
%\mathbf{\Lambda}^{T}_{1}]^{T}$ and
%$\mathbf{\Sigma}_{2}=[\mathbf{\Lambda}^{T}_{2},
%\mathbf{0}^{T}_{(N-M)\times{M}}]^{T}$ satisfy
%$\mathbf{\Sigma}_{i}=\mathbf{d}[\lambda(1),\lambda(2),...,\lambda(n)]$
%$\mathbf{\Sigma}^{2}_{1}+\mathbf{\Sigma}^{2}_{3}=\mathbf{I}_{M}$
%and $\mathbf{\Lambda}_{1}$ and $\mathbf{\Lambda}_{2}$ are two
%nonnegative diagonal matrices with with elements
%$\lambda(1)\geq{\lambda(i)}\geq{\lambda(n)}$, $i\in{[1,n]}$.
Similarly, let us introduce another decomposition, namely,
$\mathbf{QR}$ decomposition as
%of $\mathbf{H}_{1,2}$ and $\mathbf{H}_{2,2}$ shown as follow
\begin{eqnarray}\label{QR}
\left[\mathbf{H}_{1,2}\mathbf{V}_{2},
\mathbf{H}_{2,2}\mathbf{V}_{2}\right]&=&\left[\mathbf{Q}_{R,1}\mathbf{R}_{1},
\mathbf{Q}_{R,2}\mathbf{R}_{2} \right],
%&=&\mathbf{\Pi}_{R}\left[\mathbf{Q}_{2}\mathbf{R}_{1},
%\mathbf{Q}_{2}\mathbf{R}_{4}\right]
\end{eqnarray}
where $\mathbf{Q}_{R,i}\in{\mathbb{C}^{N\times{M}}}$ for $i=1,2$,
is a unitary matrix, and $\left\{\mathbf{R}_{1},
\mathbf{R}_{2}\right\} \in{\mathbb{C}^{M\times{M}}}$ are upper
triangular matrices. Substituting \eqref{LQ}, \eqref{QR}, and
\eqref{RL+RF} back into \eqref{s3}, i.e,. for $\mathbf{S}_{1}$, we
may get equivalent received signals shown as,
\begin{eqnarray}\label{equ}
\widehat{y}_{1}&=&\!\!\!\!(\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}\mathbf{R}_{1}+\mathbf{L}^{T}_{2}\mathbf{F}_{D,2}\mathbf{R}_{2})x_{2}\notag\\
&&\!\!\!\!+\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}n_{R_{1}}\!+\!\mathbf{L}^{T}_{2}\mathbf{F}_{D,2}n_{R_{2}}\!+\!n_{2}\notag\\
&=&\!\!\!\!\mathbf{\widehat{H}}_{1,2}x_{2}+\widehat{n}_{1},
\end{eqnarray}
where
$\mathbf{\widehat{H}}_{1}=\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}\mathbf{R}_{1}+\mathbf{L}^{T}_{2}\mathbf{F}_{D,2}\mathbf{R}_{2}$
and
$\widehat{n}_{1}=\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}n_{R_{1}}\!+\!\mathbf{L}^{T}_{2}\mathbf{F}_{D,2}n_{R_{2}}\!+\!n_{2}$
are efficient channel and noise coefficients, obtained from the
covariance of $\widehat{n}_{1}$, we have
\begin{eqnarray}\label{equa}
\mathbf{\widehat{C}}_{1}&=&\widehat{n}_{1}\widehat{n}^{H}_{1}\notag\\
&=&\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}\mathbf{F}^{H}_{D,1}\mathbf{L}^{*}_{1}+\mathbf{L}^{T}_{2}\mathbf{F}_{D,2}\mathbf{F}^{H}_{D,2}\mathbf{L}^{*}_{2}+\mathbf{I}_{N}.
\end{eqnarray}

For fixed $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$, using \eqref{equ}
and \emph{Property 1}, the optimal problem \eqref{MSEopt} becomes
\begin{eqnarray}\label{MSEnew}
&\max\limits_{\mathbf{F}_{D,1},
\mathbf{F}_{D,2}}&~~~tr\left(\mathbf{I}_{N}+\mathbf{\widehat{H}}^{H}_{1}\mathbf{\widehat{C}}_{1}^{-1}\mathbf{\widehat{H}}_{1}\right)\label{MSEnew}\\
&s.t.&tr{(\textbf{F}^{H}_{D,1}\mathbf{D}_{1}\textbf{F}_{D,1}\!+\!\textbf{F}^{H}_{D,2}\mathbf{D}_{2}\textbf{F}_{D,2})}\!\leq\!{P_{R}}.
%&s.t.&~~~tr{(\textbf{F}^{H}_{D,2}\mathbf{D}_{2}\textbf{F}_{D,2})}\leq{(1-a)P_{R}},
\end{eqnarray}
%where
%$\widehat{\mathbf{D}}_{1}=P_{1}\mathbf{L}^{H}_{1}\mathbf{L}_{1}+P_{2}\mathbf{R}^{H}_{2}\mathbf{R}_{2}+\mathbf{I}_{M}$
%and
%$\widehat{\mathbf{D}}_{2}=P_{1}\mathbf{L}^{H}_{3}\mathbf{L}_{3}+P_{2}\mathbf{R}^{H}_{4}\mathbf{R}_{4}+\mathbf{I}_{M}$.
Then, \eqref{MSEnew} can be represented as
\begin{eqnarray}\label{MSEnew1}
tr\left((\mathbf{\widehat{H}}^{H}_{1}\mathbf{\widehat{C}}_{1}^{-1}\mathbf{\widehat{H}}_{1})+n\right),
\end{eqnarray}
where the lemma $tr(\mathbf{A+B})=tr(\mathbf{A})+tr(\mathbf{B})$
has been used. Since the matrix $\mathbf{\widehat{C}}_{1}$ is
Hermitian and positive definite, we can decompose this matrix
using Cholesky factorization as
\begin{eqnarray}\label{Choleskey}
\mathbf{\widehat{C}}_{1}&=&\mathbf{\Xi}^{H}_{1}\mathbf{\Xi}_{1}
%&=&{\left[\begin{array}{cc}
             % \mathbf{R_{4}\Xi^{-1}_{1}} & \mathbf{*} \\
             % \mathbf{0} & \mathbf{R_{2}\Xi^{-1}_{2}} \\
              %\end{array}\right]}{\left[\begin{array}{cc}
             % \mathbf{R_{2}\Xi^{-1}_{1}} & \mathbf{0} \\
              %\mathbf{*} & \mathbf{R_{4}\Xi^{-1}_{2}} \\
             % \end{array}\right]},
\end{eqnarray}
where $\mathbf{\Xi}_{1}$ denote a lower triangular matrix.
%also $\widetilde{\mathbf{H}}$ can be represented as
%\begin{eqnarray}\label{BD}
%\widetilde{\mathbf{H}}={\left[\begin{array}{cc}
 %             {\mathbf{R}^{H}_{2}} {\mathbf{R}^{H}_{4}} \\
 %             \end{array}\right]}{\left[\begin{array}{cc}
 %             {\mathbf{A}_{1}} & \mathbf{0} \\
 %             \mathbf{0} & {\mathbf{A}_{2}} \\
 %             \end{array}\right]}{\left[\begin{array}{cc}
 %             {\mathbf{\Sigma}_{1}} \\
 %             {\mathbf{\Sigma}_{3}} \\
 %            \end{array}\right]}.
%\end{eqnarray}
By substituting \eqref{Choleskey} back into \eqref{MSEnew1}, we
can simply rewrite the optimal problem as
\begin{eqnarray}\label{MSEnew2}
&&\mathrm{max}~\left(\mathbf{MSE}^{o}_{1}\right)^{-1}\!\!\notag\\
&=&\mathrm{max}~tr\left((\mathbf{\widehat{H}}^{H}_{1}\left(\mathbf{\Xi}^{H}_{1}\mathbf{\Xi}_{1}\right)^{-1}\mathbf{\widehat{H}}_{1})+n\right)\notag\\
&\overset{a}{=}&\mathrm{max}~tr\left(\left(\mathbf{\widehat{H}}^{H}_{1}\mathbf{\Xi}^{-1}_{1}\right)\left(\mathbf{\widehat{H}}^{H}_{1}\mathbf{\Xi}^{-1}_{1}\right)^{H}\right)\notag\\
&=&\mathrm{max}~tr\left(\mathbf{B}_{1}\mathbf{B}^{H}_{1}\right),
\end{eqnarray}
where $\overset{a}{=}$ denotes $n$ has nothing to do with the
maximum solution and
$\mathbf{B}_{1}=\mathbf{\widehat{H}}^{H}_{1}\mathbf{\Xi}^{-1}_{1}$.
Thus, the optimal problem can be represented as the determinant
maximization of $|\mathbf{B}_{1}|^{2}$.

In \emph{Case a}, since $\mathbf{F}_{D,i}$ is the block diagonal
matrix, its determinant can be written as
\begin{eqnarray}\label{target}
\mathbf{det}{\mathbf{F}_{D,i}}=\mathbf{det}{\mathbf{F}_{D,i,1}}\cdot\mathbf{det}{\mathbf{F}_{D,i,2}}.
\end{eqnarray}
Let $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, and $\mathbf{D}$ be
an $\frac{N}{4}\times\frac{N}{4}$ matrix. We can define
$\mathbf{det}{\mathbf{F}_{D,i,i}}$, for $i\in{1,2}$ as,
\begin{eqnarray}\label{target1}
\mathbf{det}{\mathbf{F}_{D,i,i}} &=& {\left|\begin{array}{cc}
              \mathbf{A} & \mathbf{D} \\
             \mathbf{B} & \mathbf{C} \\
             \end{array}\right|}\notag\\
&=& {\left|\begin{array}{cc}
              \mathbf{A} & \mathbf{0} \\
             \mathbf{B} & \mathbf{I} \\
             \end{array}\right|}
             {\left|\begin{array}{cc}
              \mathbf{I} & \mathbf{A}^{-1} \\
             \mathbf{0} & \mathbf{C}-\mathbf{BA^{-1}D} \\
             \end{array}\right|}\notag\\
&=& |\mathbf{A}|\left|\mathbf{C}-\mathbf{BA^{-1}D}\right|,
             \end{eqnarray}
where $\mathbf{I}$ stands for an $\frac{N}{4}\times\frac{N}{4}$
identity matrix. In \eqref{target1}, to obtain maximum
$\mathbf{det}{\mathbf{F}_{D,i,i}}$, we should minimize
$\mathbf{BA^{-1}D}$. Let us introduce the SVD of $\mathbf{B}$,
$\mathbf{A}$, and $\mathbf{D}$ as
\begin{eqnarray}\label{target3}
\mathbf{B}=\mathbf{U}_{B}\mathbf{\Sigma}_{B}\mathbf{\Lambda}^{H}_{B},~
\mathbf{A}=\mathbf{U}_{A}\mathbf{\Sigma}_{A}\mathbf{\Lambda}^{H}_{A},~
\mathbf{D}=\mathbf{U}_{D}\mathbf{\Sigma}_{D}\mathbf{\Lambda}^{H}_{D},
\end{eqnarray}
where $\mathbf{U}_{i}$, $\mathbf{\Lambda}_{i}$, $i\in{\left\{A, B,
D\right\}}$, are the unitary matrices, and $\mathbf{\Sigma}_i$ is
an $\frac{N}{4}\times\frac{N}{4}$ diagonal matrix. Substituting
\eqref{target3} back into $\mathbf{BA^{-1}D}$, we have
\begin{eqnarray}\label{target4}
&&\min~tr\left(\mathbf{U}_{B}\mathbf{\Sigma}_{B}\mathbf{\Lambda}^{H}_{B}\left(\mathbf{U}_{A}\mathbf{\Sigma}_{A}\mathbf{\Lambda}^{H}_{A}\right)^{-1}\mathbf{U}_{D}\mathbf{\Sigma}_{D}\mathbf{\Lambda}^{H}_{D}\right)\notag\\
&\triangleq&\min~tr\left(\mathbf{\Sigma}_{B}\left(\mathbf{\Sigma}_{A}\right)^{-1}\mathbf{\Sigma}_{D}\right)\notag\\
&\triangleq&\min_{b_{i},d_{i}}~\sum^{n/4}_{i}\frac{b_{i}d_{i}}{a_{i}},
\end{eqnarray}
where $b_{i}$, $d_{i}$, and $a_{i}$ are the diagonal elements of
$\mathbf{\Sigma}_{B}$, $\mathbf{\Sigma}_{D}$ and
$\mathbf{\Sigma}_{A}$, respectively. To simplify our discussion,
we assume $\mathbf{F}_{D,i,i}$ is a semi-positive matrix, thus, we
have the minimum solution as $b_{i}d_{i}=0$. Interestingly, if
both $b_{i}$ and $d_{i}$ are $0$, $\mathbf{F}_{D,i,i}$ is a
diagonal matrix. Otherwise, it is a lower/upper triangular matrix.
In addition, for $\mathbf{S}_{1}$, the equivalent channel
$\widehat{\mathbf{H}}_{1}$, since the terms $\mathbf{L}^{T}_{1}$,
$\mathbf{L}^{T}_{2}$, $\mathbf{R}_{1}$, and $\mathbf{R}_{2}$ are
upper triangular matrices, the optimal $\mathbf{F}_{D,i}$ should
be an upper triangular matrix. Since the equivalent channel
$\widehat{\mathbf{H}}_{2}$, $\mathbf{L}_{1}$, $\mathbf{L}_{2}$,
$\mathbf{R}^{T}_{1}$, and $\mathbf{R}^{T}_{2}$ are lower
triangular matrices for $\mathbf{S}_{2}$, the optimal
$\mathbf{F}_{D,i}$ is a lower triangular matrix. Therefore, if and
only if $\mathbf{F}_{D,i}$ is a diagonal matrix,
the sum-MSE is optimal in our proposed method. %Then, the MMSE
%based $\mathbf{F}_{D,i}$ is given by
%\begin{eqnarray}\label{MSEfilter}
%\mathbf{F}^{mmse}_{D,i}=\left(\mathbf{RL}^{T}\right)^{H}\left(\left(\mathbf{RL}^{T}\right)\left(\mathbf{RL}^{T}\right)^{H}\right)^{-1},
%\end{eqnarray}
This completes the proof for \emph{Lemma 2}.
%$\mathbf{F}_{D,1}$ and $\mathbf{F}_{D,2}$ are block diagonal
%matrix, it's simply that, for $\widehat{y}_{1}$, the equivalent
%channel $\widehat{\mathbf{H}}_{1,2}$ can be obtained as
%\begin{eqnarray}\label{OP3}
%\widehat{\mathbf{H}}_{1,2} &\triangleq& {\left[\begin{array}{cc}
%              \mathbf{U}_{D,i,1} & \mathbf{0} \\
%             \mathbf{0} & \mathbf{U}_{D,i,2} \\
%             \end{array}\right]},
%             \end{eqnarray}
%where $\mathbf{U}_{D,i,1}$ and $\mathbf{U}_{D,i,2}$ are
%$\frac{N}{2}\times \frac{N}{2}$ up triangular matrices. For
%$\widehat{y}_{2}$, the equivalent channel is
%\begin{eqnarray}\label{OP5}
%\widehat{\mathbf{H}}^{T}_{1,2} &\triangleq&
%{\left[\begin{array}{cc}
%              \mathbf{U}^{T}_{D,i,1} & \mathbf{0} \\
%             \mathbf{0} & \mathbf{U}^{T}_{D,i,2} \\
%             \end{array}\right]}.
%             \end{eqnarray}
%this solution can efficient help us to solve determinant
%maximization problem in the rest of the paper.
%
%\emph{Proof}: Since the noise power have been assumed with zero
%mean and unit variance, the determinant of effective channel to
%noise coefficient ratio is greater, the signal power will be
%higher, it means the SNR will be increase\footnote{The discussion
%is
%$\frac{J_{1}J^{H}_{1}}{J_{2}J^{H}_{2}}=\left|\frac{J_{1}}{J_{2}}\right|^{2},$
%where $J_{1}$ and $J_{2}$ are nonsingular matrix, it's simply that
%if $\left|\frac{J_{1}}{J_{2}}\right|$ is maximum,
%$\frac{J_{1}J^{H}_{1}}{J_{2}J^{H}_{2}}$ will achieve optimal.}.
%For simplicity of discussion, we only consider the case of $M=N$,
%the extension of our discussions to the case of $M\neq N$ is
%straightforward. The equivalent channel can be written as
%\begin{equation}\label{equivaleng signal}
%{\mathbf{\widehat{H}}}_{1,2} = {\left[\begin{array}{cc}
%              \mathbf{A} & \mathbf{0} \\
%              \mathbf{B} & \mathbf{C} \\
 %             \end{array}\right]}_{N\times N}
%              {\left[\begin{array}{cc}
%              \mathbf{M} & \mathbf{N} \\
%              \mathbf{P} & \mathbf{Q} \\
%              \end{array}\right]}_{N\times N}
 %              {\left[\begin{array}{cc}
 %             \mathbf{E} & 0 \\
%              \mathbf{F} & \mathbf{G} \\
%              \end{array}\right]}_{N\times N},
%              \end{equation}
%              and
%\begin{equation}\label{equivaleng signal}
%{\mathbf{\widehat{n}}}_1 = {\left[\begin{array}{cc}
%              \mathbf{A} & \mathbf{0} \\
%              \mathbf{B} & \mathbf{C} \\
 %             \end{array}\right]}_{N\times N}
%              {\left[\begin{array}{cc}
%              \mathbf{M} & \mathbf{N} \\
%              \mathbf{P} & \mathbf{Q} \\
%              \end{array}\right]}_{N\times N}
%               {\left[\begin{array}{cc}
%              \mathbf{U} & \mathbf{K} \\
%              \mathbf{T} & \mathbf{Z} \\
%              \end{array}\right]}_{N\times N},
%              \end{equation}
%obviously, with inequality of $\mathrm{Det}|\mathbf{ABC}|\leq
%\mathrm{Det}{|\mathbf{A}|\cdot \mathrm{Det}|\mathbf{B}|\cdot
%\mathrm{Det}|\mathbf{C}|}$, the determinant of efficient channel
%can be rewritten as
%\begin{eqnarray}\label{equivalent}
%\mathrm{det}|{\mathbf{\widetilde{H}}}_i| &\leq&
%\mathrm{det}{\left|\begin{array}{cc}
%              \mathbf{A} & \mathbf{0} \\
%              \mathbf{B} & \mathbf{C} \\
% %             \end{array}\right|}
%              \mathrm{det}{\left|\begin{array}{cc}
%              \mathbf{M} & \mathbf{N} \\
%              \mathbf{P} & \mathbf{Q} \\
%              \end{array}\right|}
%               \mathrm{det}{\left|\begin{array}{cc}
%              \mathbf{E} & 0 \\
%              \mathbf{F} & \mathbf{G} \\
%              \end{array}\right|}\notag\\
%              &=&
%\left|\mathbf{AC}\right|\left|\mathbf{E}\mathbf{G}\right|\left|\mathbf{MQ}-\mathbf{NP})\right|\notag\\
%&\leq&
%|\mathbf{AC}|\left|\mathbf{E}\mathbf{G}\right||\mathbf{MQ}|.
%{\mathbf{\widetilde{H}}}_i &\triangleq& {\left[\begin{array}{cc}
              %\mathbf{A} & \mathbf{0} \\
              %\mathbf{B} & \mathbf{C} \\
              %\end{array}\right]}_{N\times N}
              %{\left[\begin{array}{cc}
             % \mathbf{M} & \mathbf{N} \\
              %\mathbf{P} & \mathbf{Q} \\
             % \end{array}\right]}_{N\times M}\notag\\
            %  &=&{\left[\begin{array}{cc}
            %  \mathbf{AM} & \mathbf{AN} \\
            %  \mathbf{BM+PC} & \mathbf{BN+CQ} \\
            %  \end{array}\right]}_{N\times M}.
 %             \end{eqnarray}
%The determinate of \eqref{equivalent} is
%\begin{eqnarray}\label{equivalent signal}
%{\mathbf{\widetilde{H}}}_i &\triangleq& \left|
%\begin{array}{cc}
%\mathbf{AM} & \mathbf{AN} \\
%\mathbf{BM+PC} & \mathbf{BN+CQ}\notag\\
%\end{array}
%\right|\notag\\
%&=&|\mathbf{AM(BN+CQ)-AN(BM+PC)}|\notag\\
%&\leq& \mathbf{|A||M(BN+CQ)-N(BM+PC)|}\notag\\
%&=& |\mathbf{A}||\mathbf{MCQ}-\mathbf{NPC})|\notag\\
%&\leq& |\mathbf{A}||\mathbf{C}||\mathbf{MQ}-\mathbf{NP})|\notag\\
%&\leq& |\mathbf{A}||\mathbf{C}||\mathbf{MQ}|.
 %\end{eqnarray}
%It's interestingly \eqref{equivalent} can achieve upper bound is
%if and only if $\mathbf{NP}=0$, otherwise matrix $\mathbf{A}$
%should be formed as
%\begin{eqnarray}\label{A}
%\mathbf{A} &\triangleq& {\left[\begin{array}{cc}
%              \mathbf{M} & \mathbf{0} \\
%              \mathbf{P} & \mathbf{Q} \\
%              \end{array}\right]}~,~
%{\left[\begin{array}{cc}
%              \mathbf{M} & \mathbf{0} \\
%              \mathbf{0} & \mathbf{Q} \\
%              \end{array}\right]}~\mathrm{or}~
%{\left[\begin{array}{cc}
%              \mathbf{M} & \mathbf{N} \\
%              \mathbf{0} & \mathbf{Q} \\
%              \end{array}\right]}.
%              \end{eqnarray}
%\emph{Remark}: The first term is low triangular matrix with
%eigenvalues distributed in diagonal part, the third term is
%diagonal matrix, thus, only matrix $\mathbf{A}$ is low triangular
%or diagonal matrix, we may fast compute eigenvalues of
%$\mathbf{S}_{i}$, this completes proof.

\emph{Property 2}: For any $M\times N$ rectangular matrices
$\mathbf{G}$ and $\mathbf{J}$, matrices $\mathbf{A}$ and
$\mathbf{B}$ are lower/upper triangular matrices based on QR or QL
decomposition of $\mathbf{G}$ and $\mathbf{J}$. If
$a_{i,i}+b_{i,i}\neq0$, where $a_{i,i}$ and $b_{i,i}$ are diagonal
elements of matrices $\mathbf{A}$ and $\mathbf{B}$, respectively,
we can easily obtain
\begin{eqnarray}\label{det}
\mathrm{det}\left(\mathbf{A}+\mathbf{B}\right)=\prod^{m}_{i=1}(a_{i,i}+b_{i,i})
\geq\mathrm{det}\mathbf{A}+\mathrm{det}\mathbf{B},
\end{eqnarray}
Consequently, we have
\begin{eqnarray}\label{det1}
\mathrm{det}\mathbf{\widehat{H}}^{H}_{1}&=&\prod^{m}_{i=1}\left(l_{1,i,i}f_{D,1,i}r_{1,i,i}+l_{2,i,i}f_{D,2,i}r_{2,i,i}\right)\notag\\
%&=&\prod^{m}_{i=1}\left(l_{1,i,i}\tau f_{D,2,i}r_{2,i,i}+l_{3,i,i}f_{D,2,i}r_{4,i,i}\right)\notag\\
&=&\prod^{m}_{i=1}(\varsigma_{1}+\varsigma_{2}),
\end{eqnarray}
where $\varsigma_{1}=l_{1,i,i}f_{D,1,i}r_{1,i,i}$,
$\varsigma_{2}=l_{2,i,i}f_{D,2,i}r_{2,i,i}$, $l_{1,i,i}$,
$f_{D,1,i}$, $r_{1,i,i}$, $l_{2,i,i}$, $f_{D,2,i}$, and
$r_{2,i,i}$ are diagonal elements of $\mathbf{L}_{1}$,
$\mathbf{F}_{D,1}$, $\mathbf{R}_{1}$, $\mathbf{L}_{2}$,
$\mathbf{F}_{D,2}$, and $\mathbf{R}_{2}$, respectively. Since
$\mathbf{\Xi}_{1}$ and $\mathbf{\Xi}_{1}^{-1}$ are also lower
triangular matrices, we have
\begin{eqnarray}\label{det2}
\mathrm{det}\mathbf{B}_{1}&=&\mathrm{det}\mathbf{\widehat{H}}^{H}_{1}\mathrm{det}\mathbf{\Xi}_{1}^{-1}\notag\\
&=&\prod^{m}_{i=1}(\varsigma_{1}+\varsigma_{2})\xi_{i}
\end{eqnarray}
where $\xi_{i}$ is the diagonal element of
$\mathbf{\Xi}^{-1}_{1}$.
%\begin{eqnarray}\label{MSEnew2}
%tr\left(\sum^{n}_{i=1}({v^{b}_{i,1}}^{2}a^{2}_{1}\lambda^{2}_{1}\varsigma^{2}_{1}+a^{2}_{2}\lambda^{2}_{3}\varsigma^{2}_{2})+n\right),
%\end{eqnarray}
%where $\varsigma_{i}$ denote diagonal element of low triangular
%matrix $\Xi_{i}$.
%\emph{Discussion 2}: $\mathbf{A}$ is low triangular matrix,
%obviously $\mathbf{A}$ may equivalent yield as
%\begin{eqnarray}\label{LTA}
%\mathbf{A}=\mathbf{Q}_{2}{\left[\begin{array}{cc}
%              {\Xi_{1}}\Sigma^{-1}_{1} & \mathbf{0} \\
%               *& {\Xi_{2}}\Sigma^{-1}_{3} \\
%              \end{array}\right]}\mathbf{U}^{H}_{1},
%\end{eqnarray}
%where $*$ is any matrix which satisfy
%$|*|=1/2\sqrt{\Xi_{1}\Xi_{2}\Sigma^{-1}_{1}\Sigma^{-1}_{3}}$.
%
%If $\mathbf{A}$ is diagonal matrix, it may be written as
%\begin{eqnarray}\label{LTA1}
%\mathbf{A}=\mathbf{Q}_{2}{\left[\begin{array}{cc}
%              {\Xi_{1}}\Sigma^{-1}_{1} & \mathbf{0} \\
%              \mathbf{0} & {\Xi_{2}}\Sigma^{-1}_{3} \\
%              \end{array}\right]}\mathbf{U}^{H}_{1}.
%\end{eqnarray}
Now, the optimization problem can be reformulated as
\begin{eqnarray}
&\max\limits_{\mathbf{F}_{D,1},\mathbf{F}_{D,2}}&~~~~~~~~~~~~~~~~~~~~|\mathbf{B}_{1}|^{2}\label{optnew}\\
&s.t.&{\left(\!\!\!\begin{array}{cc}
              |\mathbf{B}_{1}|^{2} & \mathbf{I}_{N} \\
              \mathbf{I}_{N} & \!\!\!\mathbf{I}_{N}+\mathbf{\Xi}^{-1}_{1}\widehat{\mathbf{H}}^{H}_{1}\widehat{\mathbf{H}}_{1}\left(\mathbf{\Xi}^{-1}_{1}\right)^{H} \\
              \end{array}\!\!\!\!\!\right)}\label{optnew1}\\
&s.t.&tr{(\textbf{F}_{D,1}\mathbf{D}_{1}}\textbf{F}^{H}_{D,1}+\textbf{F}_{D,2}\mathbf{D}_{2}\textbf{F}^{H}_{D,2})\!\leq\!{P_{R}}\label{optnew2},
%&s.t.&~~~tr{(\textbf{F}^{H}_{D,2}\widehat{\mathbf{D}}_{2}\textbf{F}_{D,2})}\leq{(1-a)P_{R}}\label{optnew3},
\end{eqnarray}
It is clear that \eqref{optnew}-\eqref{optnew2} is a convex
problem for beamformer vectors $\mathbf{F}_{D,1}$ and
$\mathbf{F}_{D,2}$, which can be efficiently solved by the
interior-point method \cite{a7}.
%thus the KKT conditions are necessary and sufficient to find the
%optimal solution. The Lagrangian of problem \eqref{optnew} can be
%form as
%\begin{eqnarray}\label{mmm}
%\textit{L}&=&\!\!\!\!\!\mathbf{F}_{D,2}\Delta
%\mathbf{F}_{D,2}^{H}\!+\!\!\lambda_{1}\!\!
%\left(\mathrm{diag}(\tau)\mathbf{F}_{D,2}\mathbf{\widetilde{D}}_{1}
%\mathbf{F}_{D,2}^{H}\mathrm{diag}(\tau)-aP_{R}\right)\notag\\
%&&\!\!\!\!\!+\lambda_{2}
%%\left(\mathbf{F}_{D,2}\mathbf{\widetilde{D}}_{2}
%\mathbf{F}_{D,2}^{H}-(1-a)P_{R}\right),
%\end{eqnarray}
%where $\mathbf{F}_{D,2}=[f_{D,2,1},...f_{D,2,n}]$,
%\begin{equation}\label{mmm1}
%\Delta={\left[\begin{array}{ccc}
%              \left|\varsigma_{1}\xi_{1}\right|^{2} & \ldots & 0 \\
%              \vdots & \ddots &  \vdots \\
%              0      &  \ldots &   \left|\varsigma_{n}\xi_{n}\right|^{2} \\
%              \end{array}\right]},
%\end{equation}
%\begin{equation}\label{mmm2}
%\mathbf{\widetilde{D}}_{1}=\!\!{\left[\begin{array}{ccc}
%              \!\!|l_{1,1,1}|^{2}\!+\!|r_{2,1,1}|^{2}\!+\!1 & \ldots & 0 \\
%              \vdots & \ddots &  \vdots \\
%              0      &  \ldots &  |l_{1,n,n}|^{2}\!+\!|r_{2,n,n}|^{2}\!+\!1 \\
%              \end{array}\!\!\right]}
%\end{equation}
%and
%\begin{equation}\label{mmm2}
%\mathbf{\widetilde{D}}_{2}=\!\!{\left[\begin{array}{ccc}
%              \!\!|l_{3,1,1}|^{2}\!+\!|r_{4,1,1}|^{2}\!+\!1 & \ldots & 0 \\
%              \vdots & \ddots &  \vdots \\
%              0      &  \ldots &  |l_{3,n,n}|^{2}\!+\!|r_{4,n,n}|^{2}\!+\!1 \\
%              \end{array}\!\!\right]}.
%\end{equation}
%The corresponding KKT conditions are given by
%\begin{eqnarray}\label{mmm3}
%&&\sum^{i=n}_{i=1}\left(\mathbf{f}^{2}_{D,2,i}\Delta_{i}+\lambda_{1}
%\mathbf{f}^{2}_{D,2,i}\mathbf{d}_{1,i}+\lambda_{2}\tau^{2}\mathbf{f}^{2}_{D,2,i}\mathbf{d}_{2,i}\right)=\mathbf{P}_{R}\notag\\
%&&\lambda_{1}\left(\mathrm{diag}(\tau)\mathbf{F}_{D,2}\mathbf{\widetilde{D}}_{1,i}
%\mathbf{F}_{D,2}^{H}\mathrm{diag}(\tau)-aP_{R}\right)=0\notag\\
%&&\lambda_{2} \left(\mathbf{F}_{D,2}\mathbf{\widehat{D}}_{2,i}
%\mathbf{F}_{D,2}^{H}-(1-a)P_{R}\right)=0,
%\end{eqnarray}
%where $\Delta_{i}$ is element of matrices $\Delta$. By further
%recalling optimal phase in the first term of \eqref{mmm3}, the
%explicit solution structure for $\mathbf{f}_{D,2,i}$, which is
%given as
%\begin{eqnarray}\label{mmm4}
%\mathbf{f}_{D,2,i}=\sqrt{\frac{\mathbf{P}_{R}}{N}\left(\Delta+\lambda_{1}
%\mathbf{d}_{1,i}+\lambda_{2}\tau^{2}\mathbf{d}_{2,i}\right)^{-1}},
%\end{eqnarray}
%where we have assumed each channel is with uncial power. The above
%KKT conditions \eqref{mmm4} are easy to solve since, for any given
%$\lambda$, \eqref{mmm4} is a univariate cubic equation of
%$\mathbf{f}_{D,2,i}$.

In summary, we outline the iterative beamforming design algorithm
as follows ($\mathbf{QL}-\mathbf{QR}~\mathbf{Algorithm}$):

\alglanguage{pseudocode}
\begin{algorithm}[h]
\caption{QL-QR Algorithm} $\mathbf{1}$.~~$\mathbf{Initialize}$:
$\mathbf{F}^{(n)}_{i}$, $\mathbf{W}^{(n)}_{i}$,
$\mathbf{V}^{(n)}_{i}$, $\mathbf{H}^{(n)}_{i,j}$,
for $i,j=1,2$, set $n=0$;\\
$\mathbf{2}$.~~$\mathbf{Repeat}$:
\begin{algorithmic}[1]
\For {$n\leftarrow n+1$}\\
for given $\mathbf{F}^{(n-1)}_{i}$, $\mathbf{H}^{(n-1)}_{i,j}$ and
$\mathbf{V}^{(n-1)}_{i}$ update $\mathbf{W}^{(n)}_{i}$ using
(22);\\
for fixed $\mathbf{F}^{(n-1)}_{i}$, $\mathbf{H}^{(n-1)}_{i,j}$,
$\mathbf{W}^{(n)}_{i}$ update $\mathbf{V}^{(n)}_{i}$ using
(35);\\
decompose $(\mathbf{H}^{(n)}_{1,1}\mathbf{V}^{(n)}_{1},
\mathbf{H}^{(n)}_{2,1}\mathbf{V}^{(n)}_{1})$ using $\mathbf{QL}$
decomposition as (42);\\
decompose $(\mathbf{H}^{(n)}_{1,2}\mathbf{V}_{2}^{(n)},
\mathbf{H}^{(n)}_{2,2}\mathbf{V}^{(n)}_{2})$ using $\mathbf{QR}$
decomposition as
(43);\\
compute
$\widehat{\mathbf{C}}^{(n)}_{1}$;\\
decompose $\widehat{\mathbf{C}}^{(n)}_{1}$ using Cholesky decomposition as (49);\\
for fixed $\mathbf{\widehat{H}}_{1}$, compute $\mathrm{det}\mathbf{B}^{(n)}_{1}$ using (57);\\
$\mathbf{Until}$ $\mathbf{MSE}_{1}$ converges.
      \EndFor
\end{algorithmic}
\end{algorithm}

Since in the QL-QR algorithm, the solution of each subproblem is
optimal, we conclude that the total MSE value is decreased as the
number of iterations increases. Meanwhile, the total MSE is lower
bounded.

\begin{figure}
\begin{center}
\includegraphics [width=70mm,height=24mm]{Gneral.eps}
\caption{\label{4} The extended system model. (a): The $K$ pair
source nodes scenario. (b): The $T$ relay nodes scenario.}
\end{center}
\end{figure}

\emph{Discussion}: The extended two system models are shown in
Fig. \ref{4}, which are the multipair scenario with two relay
nodes and $K$ pair source nodes, and the $Z$ ($Z$ should be even
number) relay nodes scenario with two source nodes. In Figure
4(a), each pair sources and two relay nodes can be seen as a
group. Since each pair source nodes are independent with each
other, we can design that there are $K$ RFs , LFs and one CF
equipped at each relay nodes. Therefor, the extended system model
(a) can be seen as $K$ parallel of our proposed system model. In
Figure 4(b), two source nodes and every two relay nodes can be
seen as one group. Obviously, the extended system model (b) can be
seen as $\frac{Z}{2}$ parallel of our proposed system model.

%\subsection{Sum of Rate Discussion}
%
%In this subsection, we carry out an analysis of the sum of rate.
%The achievable sum of rate (ASR) of an AF MIMO TWRC, which is
%given by
%\begin{eqnarray}\label{s7}
%\mathbf{R}_{sum}&=&\frac{a}{2}~\mathbf{R}_{S_{1}}+\frac{1-a}{2}~\mathbf{R}_{S_{2}}\notag\\
%&=&\frac{a}{2}~\mathbf{E}~\mathrm{log}_{2}\left(1+\mathbf{SNR}_{1}\right)\notag\\
%&&+\frac{1-a}{2}~\mathbf{E}~\mathrm{log}_{2}\left(1+\mathbf{SNR}_{2}\right),
%\end{eqnarray}
%where
%\begin{eqnarray}\label{s9}
%\mathbf{SNR}_{1}&=&\frac{P_{1}(\mathbf{L}^{T}_{1}\mathbf{F}_{D,1}\mathbf{R}_{2}+\mathbf{L}^{T}_{3}\mathbf{F}_{D,2}\mathbf{R}_{4})(\Phi
%_{2}\Sigma _{1}+\Phi _{4}\Sigma _{3})^{H}}{\Phi _{2}\Phi
%^{H}_{2}\sigma^{2}_{n_{1}}+\Phi _{4}\Phi
%^{H}_{4}\sigma^{2}_{n_{1}}+\mathrm{I}_{N}\sigma^{2}_{1}}\notag\\
%&=&\frac{\sum^{m}_{i=1}p_{1,i}\left(r^{2}_{2,i}\lambda^{2}_{1,i}a^{2}_{1,i}+r^{2}_{4,i}\lambda^{2}_{2,i}a^{2}_{3,i}\right)}{\sum^{m}_{i=1}(r^{2}_{2,i}a^{2}_{1,i}\sigma^{2}_{r,1}+r^{2}_{4,i}a^{2}_{3,i}\sigma^{2}_{r,2})+I_{N}\sigma^{2}_{1}}\notag\\
%&\triangleq&\frac{\sum^{m}_{i=1}p_{1,i}\left(r^{2}_{2,i}\lambda^{2}_{1,i}a^{2}_{1,i}+r^{2}_{4,i}\lambda^{2}_{2,i}a^{2}_{3,i}\right)}{\sum^{m}_{i=1}\left(r^{2}_{2,i}a^{2}_{1,i}+r^{2}_{4,i}a^{2}_{3,i}+1\right)\sigma^{2}}
%\end{eqnarray}
%and
%\begin{eqnarray}\label{s10}
%\mathbf{SNR}_{2}&=&\frac{P_{2}(\Phi _{1}\Sigma _{2}+\Phi
%_{3}\Sigma _{4})(\Phi _{1}\Sigma _{2}+\Phi _{3}\Sigma
%_{4})^{H}}{\Phi _{2}\Phi ^{H}_{2}\sigma^{2}_{n_{1}}+\Phi _{3}\Phi
%^{H}_{3}\sigma^{2}_{n_{1}}+\mathrm{I}_{N}\sigma^{2}_{2}}\notag\\
%&=&\frac{\sum^{m}_{i=1}p_{2,i}\left(r^{2}_{1,i}\lambda^{2}_{2,i}a^{2}_{1,i}+r^{2}_{2,i}\lambda^{2}_{4,i}a^{2}_{3,i}\right)}{\sum^{m}_{i=1}(r^{2}_{2,i}a^{2}_{1,i}\sigma^{2}_{r,1}+r^{2}_{4,i}a^{2}_{3,i}\sigma^{2}_{r,2})+I_{N}\sigma^{2}_{2}}\notag\\
%&\triangleq&\frac{\sum^{m}_{i=1}p_{2,i}\left(r^{2}_{1,i}\lambda^{2}_{2,i}a^{2}_{1,i}+r^{2}_{2,i}\lambda^{2}_{4,i}a^{2}_{3,i}\right)}{\sum^{m}_{i=1}\left(r^{2}_{1,i}a^{2}_{2,i}+r^{2}_{3,i}a^{2}_{4,i}+1\right)\sigma^{2}}
%\end{eqnarray}
%where $\Phi_{i}=\mathbf{RA}$ and elements $\lambda, r, a$ are
%diagonal part of matrix $\mathbf{\Sigma}, \mathbf{R}, \mathbf{A}$.
%In \eqref{s10}, $\triangleq$ means we have assumed relay noise and
%source noise are with same variance $\sigma$. By using Jensen's
%inequality $\mathrm{E}~\mathrm{log}_{2}(1+x)\leq
%\mathrm{log}_{2}\left(1+\mathrm{E}(x)\right)$ to drive a upper
%bound of sum of rate, which is given by
%\begin{eqnarray}\label{s1}
%\mathbf{R}_{sum}&\leq&\frac{a}{2}~\mathrm{log}_{2}\left(1+~\mathbf{E}~(\mathbf{SNR}_{1})\right)\notag\\
%&&+\frac{1-a}{2}~\mathrm{log}_{2}\left(1+\mathbf{E}~(\mathbf{SNR}_{2})\right).
%\end{eqnarray}
%It's shown in \cite{a1} that at the optimal, $SNR_{1}=SNR_{2}$
%holds true, thus leading to an $SNR$ balacing result. Hence, the
%max-min optimization problem of \eqref{s7} is equivalent to
%\begin{eqnarray}\label{snr1=snr2}
%\mathbf{R}^{o}_{sum}\leq
%\frac{1}{2}\mathrm{log}_{2}\left(1+~\mathbf{E}~(\mathbf{SNR}_{1})\right).
%\end{eqnarray}
%By using expectation definition, $\mathbf{E}({a})=\int af(a)da$,
%where $f(a)$ is probability density function (PDF) of value $a$,
%we represent \eqref{snr1=snr2} as
%\begin{eqnarray}\label{s1}
%\mathbf{R}^{o}_{sum}&\leq&
%\frac{1}{2}\mathrm{log}_{2}\left(1+\int^{\infty}_{0}\mathbf{SNR}_{1}f(\mathbf{SNR}_{1})dx \right)\notag\\
%&=&\frac{1}{2}\mathrm{log}_{2}\left(1+\sum^{m}_{i=1}\int^{\infty}_{0}\frac{\rho
%x(\lambda^{2}_{2}+\lambda^{2}_{4})}{\sigma^{2}(x+1)}f_{\Gamma}(x)dx
%\right)\notag\\
%&=&\frac{1}{2}\mathrm{log}_{2}\left(1+\sum^{m}_{i=1}\psi(x-ln(1+x))f^{'}_{\Gamma}(x)\right),
%\end{eqnarray}
%where $\psi =\frac{\rho
%(\lambda^{2}_{2}+\lambda^{2}_{4})}{\sigma^{2}}$,
%$f^{'}_{\Gamma}(x)=f_{\Gamma}(x-1)$ and $x$ denote effective
%channel covariance.
%
%\begin{figure}
%\begin{center}
%\includegraphics [width=70mm,height=40mm]{TWR2.eps}
%\caption{\label{123} Proposed GLQ-GQR Amplify-and-Forward MIMO
%two-way relay system.}
%\end{center}
%\end{figure}

\section{Computational Complexity Analysis}
In this section, we measure the performance of the proposed QL-QR
scheme in terms of the computational complexity compared with
existing algorithms by using the total number of floating point
operations (FLOPs). A flop is defined as a real floating
operation, i.e., a real addition, multiplication, division, and so
on. In \cite{c1}, the authors show the computational complexity of
the real Choleskey decomposition. For complex numbers, a
multiplication followed by an addition needs 8 FLOPs, which leads
to 4 times its real computation. According to \cite{a3}, the
required number of FLOPs of each matrix is described as follows:\\\\
1. Multiplication of $m\times{n}$ and $n\times{p}$ complex
matrices: $8mnp-2mp$;\\
2. Multiplication of $m\times{n}$ and $n\times{m}$ complex
matrices: $4nm\times{(m+1)}$;\\
3. SVD of an $m\times{n}(m\leq{n})$ complex matrix where only
$\Sigma$ is obtained: $32(mn^{2}-n^{3}/3)$;\\
4. SVD of an $m\times{n}(m\leq{n})$ complex matrix where only
$\Sigma$ and $\Lambda$ are obtained: $32(nm^{2}+2m^{3})$,\\
5. SVD of an $m\times{n}(m\leq{n})$ complex matrix where
$U$,$\Sigma$, and $\Lambda$ are obtained:
$8(4n^{2}m+8nm^{2}+9m^{3})$;\\
6. Inversion of an $m\times{m}$ real matrix using Gauss-Jordan
elimination: $2m^{3}-2m^{2}+m$;\\
7. Cholesky factorization of an $m\times{m}$ complex matrix:
$8m^{3}/3$.\\
8. QR or QL decomposition of an $m\times n$ conplex matrix
$16\left(n^{2}m-nm^{2}+\frac{1}{3}m^{3}\right)$.\\

For the conventional RBD method \cite{b18}, the authors consider a
linear MU-MIMO precoding scheme for DL MIMO systems. For the
non-regenerative MIMO relay systems \cite{b19}, the authors
investigate a precoding design for a 3-node MIMO relay network. In
\cite{c2}, a relay-aided system based on a quasi-EVD channel is
proposed. We compare the required number of FOLPs of our proposed
method with conventional precoding algorithm, such as the
conventional RBD, the non-regenerative MIMO relay system, and the
CD-BD algorithm as shown in Tables I, II, III, and IV,
respectively, under the assumption that $N_{T}=N_{R}$ and
$\overline{N}_{i}=N_{T}-N_{i}$.

\begin{figure}
\begin{center}
\includegraphics [width=92mm,height=75mm]{Computational1.eps}\\
\caption{\label{5} The complexity comparisons for required FLOPs
versus the number of the users $K$.}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics [width=92mm,height=75mm]{Computational2.eps}\\
\caption{\label{6} The complexity comparisons for required FLOPs
versus the number of the receive antennas $N_{i}$ for each user.}
\end{center}
\end{figure}

For instance, the $(2,2,2)\times6$ case denotes a system with
three users $(K=3)$, where each user is equipped with two antennas
$(N_{i}=2)$ and the total number of transmit antennas is six
$(N_{T}=2\times3=6)$. The required number FLOPs of the QL-QR
algorithm, the conventional RBD, the non-regenerative MIMO relay
system, and the CD-BD algorithm are counted as 33530, 40824,
45306, 34638, respectively. From these results, we can see that
the reduction in the number of FLOPs of our proposed precoding
method is $17.87\%$, $25.99\%$, and $3.20\%$ on an individual
basis compared to the conventional RBD, the non-regenerative MIMO
relay systems, and the CD-BD algorithm. Thus, our proposed QL-QR
algorithm exhibits lower complexity than conventional algorithms.
In addition, the complexity reduces as $N_{i}$ and $N_{T}$
increase with fixed $K$.

We summarize our calculation results of the required number of
FLOPs of the alternative methods in Tables I, II, III, and IV and
show them in Figures \ref{5} and \ref{6}. Figure \ref{5} shows the
computational complexity where $N_{i}=2$ and a value of $K$
varies. And Figure \ref{6} shows the computational complexity
where $K=4$ and a value of $N_{i}$ varies. For the conventional
RBD method, the orthogonal complementary vector $\mathbf{V}_{k,0}$
requires $K$ times SVD operations. If only $\mathbf{V}_{k,0}$ is
obtained, it is not computationally efficient. In step 5, the
efficient channel
$\mathbf{H}_{eff}=\mathbf{H}_{i}\mathbf{P}^{a}_{i}$ is decomposed
by the SVD with a dimension $\mathbf{R}_{eff}\times N_{T}$, where
$\mathbf{R}_{eff}$ is the rank of $\mathbf{H}_{eff}$. In the
nonregenerative MIMO relay method and the CD-BD algorithm, two SVD
operations are performed for the channels from the source to relay
and from relay to the destination, and then the efficient channel
covariance matrix is measured. In the nonregenerative MIMO relay
method, the authors compute $\mathbf{A}$ using the EVD an then
they diagonalize $\mathbf{G}$. In the CD-BD algorithm, the authors
calculate $\mathbf{V}^{a}_{i}$ by the SVD of
$\mathbf{H}^{\dag}_{mse}$ and then they structure
$\mathbf{V}^{b}_{i}$ by using the Choleskey decomposition.

In our proposed QL-QR algorithm, we take advantage of QL and QR
decompositions instead of the SVD operation, and then we compute
an efficient channel as well as decompose a noise covariance
matrix by the Choleskey decomposition. Finally, we calculate the
determinant of $\mathbf{B}^{2}_{i}$ to solve an optimization
problem. Obviously, our proposed QL-QR algorithm outperforms
conventional algorithms in the light of the computational
complexity.


\begin{table*}[tbp]
\centering  % ������
\caption{Computational complexity of the proposed QL-QR
Algorithm.}
\begin{tabular}{| c | c | c | c |} \hline % {lccc} ��ʾ����Ԫ�ض��뷽ʽ��left-l,right-r,center-c
Step   & Operations  & FLOPS & Case: $(2,2,2)\times 6$ \\ \hline  % \hline �ڴ������滭һ����
1 &$\mathbf{V}_{1}, \mathbf{V}_{2}$ &$2\times
K(40N^{3}_{i}-24N^{2}_{i}+17N_{i})$ & 1560
\\ \hline        % \\ ��ʾ���¿�ʼһ��
2 &$\mathbf{Q}_{L,1}\mathbf{L}_{1},
\mathbf{Q}_{L,2}\mathbf{L}_{2}$
&$2\times16K(N^{2}_{T}N_{i}-N_{T}N^{2}_{i}+\frac{1}{3}N^{3}_{i})$
& 4864
\\ \hline        % \\ ��ʾ���¿�ʼһ��
3 &$\mathbf{Q}_{R,1}\mathbf{R}_{1},
\mathbf{Q}_{R,2}\mathbf{R}_{2}$
&$2\times16K(N^{2}_{T}N_{i}-N_{T}N^{2}_{i}+\frac{1}{3}N^{3}_{i})$
& 4864
\\ \hline       % & ��ʾ�еķָ���
4 &$\mathbf{H}^{T}_{1,1}\mathbf{F}_{1}\mathbf{H}_{1,2}$
&$8N^{2}_{T}N_{i}+4N_{T}N^{2}_{i}+2N_{T}N_{i}$&
696\\
\hline 5 &$\mathbf{H}^{T}_{2,1}\mathbf{F}_{2}\mathbf{H}_{2,2}$
&$8N^{2}_{T}N_{i}+4N_{T}N^{2}_{i}+2N_{T}N_{i}$ & 696\\ \hline 6 &
$\mathbf{\widehat{C}}_{1}$
 &$2K(32N^{2}_{T}N_{i}+8N_{T}N_{i}+2N^{2}_{T}-4N_{i}+3N_{T})$ & 14856 \\ \hline 7 &$(\mathbf{\Xi}^{H}_{i}\mathbf{\Xi}_{i})^{-1}$ &$K(\frac{14}{3}N^{3}_{T}-2N^{2}_{T}+N_{T})$ &
 2826
 \\
\hline 8 & $\mathrm{det}\mathbf{B}^{2}_{1}$
 &$4K(N^{3}_{T}+N^{2}_{T}+2N_{T})$ & 3168
\\ \hline Total &  & & 33530 \\
\hline
\end{tabular}
\end{table*}

\begin{table*}[tbp]
\centering  % ������
\caption{Computational complexity of the nonregenerative MIMO
relay system \cite{b19}.}
\begin{tabular}{| c | c | c | c |} \hline % {lccc} ��ʾ����Ԫ�ض��뷽ʽ��left-l,right-r,center-c
Step   & Operations  & FLOPS & Case: $(2,2,2)\times 6$ \\ \hline  % \hline �ڴ������滭һ����
1
&$\mathbf{U}^{a}_{i}\mathbf{\Sigma}^{a}_{i}\mathbf{\Lambda}^{aH}_{i}$
&$8K(4N^{2}_{T}N_{i}+8N_{T}N^{2}_{i}+9N^{3}_{i})$ & 13248
\\ \hline        % \\ ��ʾ���¿�ʼһ��
2
&$\mathbf{U}^{a}_{j}\mathbf{\Sigma}^{a}_{j}\mathbf{\Lambda}^{aH}_{i}$
&$ 8K(4N^{2}_{T}N_{i}+8N_{T}N^{2}_{i}+9N^{3}_{i})$ & 13248
\\ \hline       % & ��ʾ�еķָ���
3 &$\mathbf{H}^{H}_{i}\mathbf{H}_{i}$ &$4KN_{i}N_{T}(N_{i}+1)$&
432\\
\hline 4 &$\mathbf{H}^{H}_{j}\mathbf{H}_{j}$
&$4KN_{i}N_{T}(N_{i}+1) $ & 432\\ \hline
5 &$\mathbf{H}^{H}_{i}[\sigma^{2}_{1}\sigma^{2}_{2}(\mathbf{H}_{j}\mathbf{F})^{H}\mathbf{H}_{j}\mathbf{F}+\mathbf{I}]^{-1}\mathbf{H}_{i}$ &$2K(N^{3}_{i}+8N_{i}N^{2}_{T}+4N^{2}_{i}N_{T}+2N_{i}N_{T}-N^{2}_{i}+N_{i})$ & 4212 \\
\hline 6 &$\mathbf{V}_{A}\mathbf{\Lambda}_{A}\mathbf{V}^{H}_{A}$
 &$8K(4N^{2}_{T}N_{i}+8N_{T}N^{2}_{i}+9N^{3}_{i}+\frac{1}{2}N_{i}) $ & 13272 \\ \hline 7 &$\mathrm{diag}(\widetilde{\mathbf{G}})$ &$K[4N_{i}N_{T}(N_{i}+1)+2N^{3}_{i}-2N^{2}_{i}+N_{i}]$ & 462

\\ \hline Total &  & & 45306 \\
\hline
\end{tabular}
\end{table*}

\begin{table*}[tbp]
\centering  % ������
\caption{Computational complexity of the conventional RBD
\cite{b18}.}
\begin{tabular}{| c | c | c | c |} \hline % {lccc} ��ʾ����Ԫ�ض��뷽ʽ��left-l,right-r,center-c
Step   & Operations  & FLOPS & Case: $(2,2,2)\times 6$ \\ \hline  % \hline �ڴ������滭һ����
1
&$\mathbf{U}^{a}_{i}\mathbf{\Sigma}^{a}_{i}\mathbf{\Lambda}^{aH}_{i}$
&$32K(N_{T}\overline{N}^{2}_{i}+2\overline{N}^{3}_{i})$ & 21504
\\ \hline        % \\ ��ʾ���¿�ʼһ��
2
&$\left((\mathbf{\Sigma}^{a}_{i})^{T}\mathbf{\Sigma}^{a}_{i}+\rho^{2}\mathbf{I}\right)^{-1/2}$
&$ K(18N_{T}N^{2}_{i}-2N^{2}_{i})$ & 336
\\ \hline       % & ��ʾ�еķָ���
3 &$\mathbf{V}^{a}_{i}\mathbf{D}^{a}_{i}$ &$8KN^{3}_{T}$&
5184\\
\hline 4 &$\mathbf{H}_{i}\mathbf{P}^{a}_{i}$ &$K(8N_{T}N^{2}_{i}-2N^{2}_{i}) $ & 552\\
\hline
5 &$\mathbf{U}^{b}_{i}\mathbf{\Sigma}^{b}_{i}\mathbf{V}^{bH}_{i}$ &$64K(\frac{9}{8}N^{3}_{i}+N_{T}N^{2}_{i}+\frac{1}{2}N^{2}_{T}N_{i})$ & 13248 \\
\hline Total &  & & 40824  \\
\hline
\end{tabular}
\end{table*}

\begin{table*}[tbp]
\centering  % ������
\caption{Computational complexity of the CD-BD algorithm
\cite{c2}.}
\begin{tabular}{| c | c | c | c |} \hline % {lccc} ��ʾ����Ԫ�ض��뷽ʽ��left-l,right-r,center-c
Step   & Operations  & FLOPS & Case: $(2,2,2)\times 6$ \\ \hline  % \hline �ڴ������滭һ����
1
&$\mathbf{U}^{H}_{i,1}\mathbf{\Sigma}_{i,1}\mathbf{\Lambda}_{i,1}$
&$8K(4N^{2}_{T}N_{i}+8N_{T}N^{2}_{i}+9N^{3}_{i})$ & 13248
\\ \hline        % \\ ��ʾ���¿�ʼһ��
2
&$\mathbf{\Lambda}^{H}_{i,2}\mathbf{\Sigma}_{i,2}\mathbf{U}_{i,2}$
&$8K(4N^{2}_{T}N_{i}+8N_{T}N^{2}_{i}+9N^{3}_{i})$ & 13248
\\ \hline       % & ��ʾ�еķָ���
3
&$\mathbf{H}_{i,2}\mathbf{W}\mathbf{H}_{i,1}$&$K[8N_{i}N^{2}_{T}-2N_{i}N_{T}+4N_{i}N_{T}\times(N_{i}+1)]$&
2088\\
\hline 4 &$\mathbf{L}^{H}_{i}\mathbf{L}_{i}$ &$2K(N_{i}+2N_{T}N_{i}\times(N_{i}+1)+4N^{3}_{i}/3) $ & 508\\
\hline 5 &$\mathbf{H}^{\dag}_{mse}$
&$4N^{3}_{R}/3+12N^{2}_{R}N_{T}-2N^{2}_{R}-2N_{T}N_{R}$ & 2736 \\
\hline 6 &$\mathbf{H}_{i,i}\mathbf{V}^{a}_{i}\mathbf{V}^{b}_{i}$
&$8K[4N_{T}N^{2}_{i}-4N^{3}_{i}/3+N^{2}_{i}(N_{i}+1)]$ & 2336
\\\hline
7 &$(\mathbf{Q}_{i}\mathbf{Q}^{H}_{i}+\sigma^{2}_{i}\Psi_{i})^{-1}$ &$K[4N_{R}N_{i}\times(N_{i}+1)+3N_{i}+2N^{3}_{i}-2N^{2}_{i}]$ & 474 \\
\hline Total &  & & 34638  \\
\hline
\end{tabular}
\end{table*}

\section{Simulation Results}
In this section, we study the performance of the proposed QL-QR
algorithm for two-way MIMO relay networks. All the simulations are
performed on the assumption that all the channels are the Rayleigh
fading channel and they are independently generated following
$\sim{CN(0,1)}$. The noise variances $\sigma^{2}_{i}$ are equally
given as $\sigma^{2}$. The total relay power constraint can be
written as
\begin{eqnarray}\label{Constraint}
P_{R_{1}}+P_{R_{2}}=aP_{R}+(1-a)P_{R}=P_{R},
\end{eqnarray}
where $a\in[0,1]$ is an auxiliary value as well as a power
allocation coefficient between two relay nodes. All the simulation
results are averaged over 1000 channel trials.

\begin{figure}
\begin{center}
\includegraphics [width=92mm,height=75mm]{SMI.eps}\\
\caption{\label{7} The achieved SMI for $N=4, 2$.}
\end{center}
\end{figure}

In Fig. \ref{7}, we compare the sum mutual information (SMI) of
various MU-MIMO schemes where full CSI is known at each node. We
set $P_{1}=P_{2}=10~\mathrm{dB}$, $M=1$, and an equal power budget
for the two relays ($a=0.5$ is assumed). The negative SMI is
adopted in \cite{b15} which can be defined as
\begin{eqnarray}\label{SMI}
\mathbf{MI}_{sum}=\mathrm{log}_{2}\left|\mathbf{MSE}_{1}\right|+\mathrm{log}_{2}\left|\mathbf{MSE}_{2}\right|.
\end{eqnarray}
In our proposed method, the SMI  shown in the simulation results
is calculated as $-2~\mathrm{log}_{2}|\mathbf{B}^{2}_{i}|$ by
using \eqref{Choleskey}, \eqref{det1}, and \eqref{det2}. It can be
observed that the proposed QL-QR algorithm has the same SMI
performance as an optimal solution in \cite{b15}.


\begin{figure}
\begin{center}
\includegraphics [width=92mm,height=75mm]{11.eps}\\
\caption{\label{8} The SMI versus the number of the relays $T$.}
\end{center}
\end{figure}

Figure \ref{8} shows the performance of our proposed SMI
performance versus the number of the relays, $T$ which is even. We
consider a practical scenario with different relay power
constraints. We set $P_{R}=30~\mathrm{dB}$ and $a=0.5$. It is
clear that, for different values of $P_{1}$ and $P_{2}$, a
solution of our proposed QL-QR algorithm shows better performance
than a max-power solution.

Figure \ref{9} exhibits the BER performance of the BD water
filling, the RBD, the SVD-RBD, and our proposed QL-QR method,
where the QPSK modulation is made use of. As pointed out in
\cite{e1}, the BER performance for a MIMO precoding system is
actually determined by the energy of the transmitted signal. To
simplify our discussion, we assume $a=0$. In the RBD,
$\mathrm{det}\left(\overline{\mathbf{H}}\overline{\mathbf{H}}^{H}\right)=\prod^{m}_{i=1}\lambda^{2}_{i}$,
where $\overline{\mathbf{H}}\in{\mathbb{C}^{N\times M}}$, for
$M<N$, is an equivalent channel matrix with its eigenvalues
$\lambda_{i}$. In our proposed QL-QR method, for source node
$\mathbf{S}_{1}$, we have
$\mathrm{det}\left(\mathbf{\widehat{H}}_{1}\mathbf{\widehat{H}}^{H}_{1}\right)=\prod^{m}_{i=1}\varsigma^{2}_{1}$.
Under the stipulaton that $\mathrm{det}\mathbf{F}_{D,i}=1$, we are
able to easily obtain $\lambda_{i}=\varsigma_{1}$. Therefore, our
proposed QL-QR method has the same BER performance as the SVD-RBD
method.

\begin{figure}
\begin{center}
\includegraphics [width=92mm,height=75mm]{BER1.eps}\\
\caption{\label{9} BER performance on the Rayleigh fading channel
.}
\end{center}
\end{figure}


\section{Conclusions}

This paper studies joint optimization problem of an AF based on
the MIMO TWRC, where two source nodes exchange their messages with
two relay nodes. A relay filter has been designed, which is able
to efficiently join the source and the relay nodes. Our main
contribution is that the optimal beamforming vectors can
efficiently be computed using determinant maximization techniques
through an iterative QL-QR algorithm based on a MSE balancing
method. Our proposed QL-QR algorithm can significantly reduce the
computational complexity and has the equivalent BER performance to
the SVD-BD algorithm.


\section*{Acknowledgment}
The authors would like to thank the anonymous reviewer for their
great constructive comments that improved this paper.

\begin{thebibliography}{99}
\bibitem{b1}
J. N. Laneman, D. N. C. Tse, and G. W. Wornell, ��Cooperative
diversity in wireless networks: Efficient protocols and outage
behavior,�� IEEE Trans. Inf. Theory, vol. 50, no. 12, pp.
3062�C3080, Dec. 2004.

\bibitem{c2}
Wei Duan, Wei Song, Sang Seob Song, and Moon Ho Lee, "Precoding
Method Interference Management for Quasi-EVD Channel," The
Scientific World Journal, 28th. Aug. 2014.

\bibitem{b2}
T. M. Cover and A. A. E. Gamal, ��Capacity theorems for the relay
channel,�� IEEE Trans. Inf. Theory, vol. 25, no. 5, pp. 572�C584,
Sep. 1979.

\bibitem{b3}
Zhongyuan Zhao, Mugen Peng, Zhiguo Ding, Wenbo Wang, Hsiao-Hwa
Chen, "Denoise-and-Forward Network Coding for Two-Way Relay MIMO
Systems," IEEE Trans. Veh. Technol., vol. 63, no. 2, pp. 775-788,
Feb. 2014.

\bibitem{b4}
Sebastien Simoens, Olga Munoz-Medina, Josep Vidal, Aitor Del Coso,
"Compress-and-Forward Cooperative MIMO Relaying With Full Channel
State Information," IEEE Trans. Signal Pross., vol. 58, no. 2, pp.
781-791, Feb. 2010.

\bibitem{b5}
B. Rankov and A. Wittneben, ��Spectral efficient protocols for
halfduplex fading relay channels," IEEE J. Sel. Areas Commun.,
vol. 25, no. 2, pp. 379-389, Feb. 2007.

\bibitem{b6}
T. Koike-Akino, P. Popovski, and V. Tarokh, ��Optimized
constellations for two-way wireless relaying with physical network
coding,�� IEEE J. Sel. Areas Commun., vol. 27, no. 5, pp.
773�C787, Jun. 2009.

\bibitem{b7}
S. Xu and Y. Hua, ��Optimal design of spatial source-and-relay
matrices for a non-regenerative two-way MIMO relay system,�� IEEE
Trans. Wireless Commun., vol. 5, no. 10, pp. 1645�C1655, May 2011.

\bibitem{b8}
Y. Rong, ��Joint source and relay optimization for two-way linear
nonregenerative MIMO relay communications,�� IEEE Trans. Signal
Process., vol. 60, no. 12, pp. 6533�C6546, Dec. 2012.

\bibitem{b9}
K.-J. Lee, H. Sung, E. Park, and I. Lee, ��Joint optimization for
one and two-way MIMOAFmultiple-relay systems,�� IEEE
Trans.Wireless Commun., vol. 9, no. 12, pp. 3671�C3681, Dec. 2010.

\bibitem{b10}
H. Bolcskei, R. U. Nabar, O. Oyman, and A. J. Paulraj, ��Capacity
scaling laws in MIMO relay networks," IEEE Trans. Wireless
Commun., vol. 5, no. 6, pp. 1433-1444, June 2006.

\bibitem{b11}
K.-J. Lee, K. W. Lee, H. Sung, and I. Lee, ��Sum-rate maximization
for two-way MIMO amplify-and-forward relaying systems," in Proc.
IEEE Veh. Technol. Conf., pp. 1-5, Apr. 2009.

\bibitem{b12}
F. Roemer and M. Haardt, ��Algebraic norm-maximizing (ANOMAX)
transmit strategy for two-way relaying with MIMO amplify and
forward relays,�� IEEE Signal Process. Lett., vol. 16, no. 10, pp.
909�C912, Oct. 2009.

\bibitem{b13}
T. Unger and A. Klein, ��Duplex schemes in multiple antenna
two-hop relaying,�� EURASIP Journal on Advances in Signal
Processing, vol. 2008,Article ID 128592, 2008.

\bibitem{b14}
N. Jindal and Z. Luo, ��Capacity limits of multiple antenna
multicast,�� in Proceedings of the IEEE International Symposium on
Information Theory (ISIT ��06), pp. 1841�C1845, Seattle,Wash, USA,
July 2006.

\bibitem{b15}
R. Zhang, Y.-C. Liang, C. C. Chai, and S. Cui, ��Optimal
beamforming for two-way multi-antenna relay channel with analogue
network coding,�� IEEE Journal on Selected Areas in
Communications, vol. 27, no. 5, pp. 699�C712, 2009.

\bibitem{b16}
R. Wang and M. Tao, ��Joint source and relay precoding designs for
MIMO two-way relaying based on MSE criterion,�� IEEE. Trans.
Signal Process., vol. 60, no. 3, pp. 1352�C1365, Mar. 2012.

\bibitem{f1}
H. Wang, L. Li, L. Song, and X. Gao, ��A linear precoding scheme
for downlink multiuser MIMO precoding systems,�� IEEE Commun.
Lett., vol. 15, no. 6, pp. 653�C655, June 2011.

\bibitem{a1}
H. Park, H. J. Yang, J. Chun and R. Adve, ��A closed-form power
allocation and signal alignment for a diagonalized MIMO two-way
relay channel with linear receivers,�� IEEE Trans. Signal
Process., vol. 60, no. 11, pp. 5948-5962, Nov. 2012.

\bibitem{b17}
Veria Havary-Nassab, Shahram Shahbazpanahi and Ali Grami, "Optimal
Distributed Beamforming for Two-Way Relay Networks," IEEE Trans.
Signal Process., vol. 58, no. 3, pp. 1238-1250, Mar. 2010.

\bibitem{a5}
Pramod Viswanath, Venkat Anantharam and David N. C. Tse, "Optimal
Sequences, Power Control, and User Capacity of Synchronous CDMA
Systems with Linear MMSE Multiuser Receivers," IEEE Trans. Inf.
Theory, vol. 45, no. 6, pp. 1968�C1983, 1999.

\bibitem{a6}
Rui Zhang, Ying-Chang Liang, Chin Choy Chai and Shuguang Cui, "
Optimal Beamforming for Two-Way Multi-Antenna Relay Channel with
Analogue Network Coding," IEEE J. Sel. Areas Commun., vol. 27, no.
5, pp. 699 - 712, June 2009

\bibitem{a4}
Cristian Lameiro, Javier V��a and Ignacio Santamar��a,
"Amplify-and-Forward Strategies in the Two-Way Relay Channel With
Analog Tx�CRx Beamforming," IEEE Trans. Veh. Technol., vol. 62,
no.2, pp. 642-654, Feb. 2013.


\bibitem{a2}
Yue Rong, "Simplified Algorithms for Optimizing Multiuser
Multi-Hop MIMO Relay Systems," IEEE Trans. Commun. vol. 59, no.
10, pp. 2896-2904, Oct. 2011.

\bibitem{a7}
S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge, U.K.:
Cambridge Univ. Press, 2004.

\bibitem{c1}
X. Gao, L. Li, P. Zhang, L. Song, G. Zhang, Q. Wang and H. Tian,
"Low-complexity RBD precoding method for downlink of multi-user
MIMO system," Electronics Letters, wol. 46, no. 23, pp. 1574 -
1576, Nov. 2010.


\bibitem{a3}
G. H. Golub and C. F. Van Loan, Matrix Computations, Johns Hopkins
Studies in the Mathematical Sciences, Johns Hopkins University
Press, Baltimore, Md, USA, 3rd edition, 1996.

\bibitem{b18}
H. Wang, L. Li, L. Song, and X. Gao, ��A linear precoding scheme
for downlink multiuser mimo precoding systems,�� IEEE
Communications Letters, vol. 15, no. 6, pp. 653�C655, 2011.

\bibitem{b19}
R. Mo and Y. H. Chew, ��Precoder design for non-regenerative MIMO
relay systems,�� IEEE Transactions on Wireless Communications,
vol. 8, no. 10, pp. 5041�C5049, 2009.

\bibitem{d1}
S. Vishwanath, N. Jindal, and A. Goldsmith, ��On the capacity of
multiple input multiple output broadcast channels,�� in
Proceedings of the International Conference on Communications (ICC
��02), pp. 1444�C1450, May 2002.

\bibitem{e1}
C. B. Peel, B. M. Hochwald, and A. L. Swindlehurst, ��A
vectorperturbation technique for near capacity multiantenna
multiuser communication��part I: channel inversion and
regularization,�� IEEE Trans. Commun., vol. 52, no. 1, pp.
195��202, Jan. 2005.


\end{thebibliography}
\end{document}

1502.04860.extracted_bib
{b1 J. N. Laneman, D. N. C. Tse, and G. W. Wornell, ��Cooperative diversity in wireless networks: Efficient protocols and outage behavior,�� IEEE Trans. Inf. Theory, vol. 50, no. 12, pp. 3062�C3080, Dec. 2004.
1502.04860.extracted_bib
{c2 Wei Duan, Wei Song, Sang Seob Song, and Moon Ho Lee, "Precoding Method Interference Management for Quasi-EVD Channel," The Scientific World Journal, 28th. Aug. 2014.
1502.04860.extracted_bib
{b2 T. M. Cover and A. A. E. Gamal, ��Capacity theorems for the relay channel,�� IEEE Trans. Inf. Theory, vol. 25, no. 5, pp. 572�C584, Sep. 1979.
1502.04860.extracted_bib
{b3 Zhongyuan Zhao, Mugen Peng, Zhiguo Ding, Wenbo Wang, Hsiao-Hwa Chen, "Denoise-and-Forward Network Coding for Two-Way Relay MIMO Systems," IEEE Trans. Veh. Technol., vol. 63, no. 2, pp. 775-788, Feb. 2014.
1502.04860.extracted_bib
{b4 Sebastien Simoens, Olga Munoz-Medina, Josep Vidal, Aitor Del Coso, "Compress-and-Forward Cooperative MIMO Relaying With Full Channel State Information," IEEE Trans. Signal Pross., vol. 58, no. 2, pp. 781-791, Feb. 2010.
1502.04860.extracted_bib
{b5 B. Rankov and A. Wittneben, ��Spectral efficient protocols for halfduplex fading relay channels," IEEE J. Sel. Areas Commun., vol. 25, no. 2, pp. 379-389, Feb. 2007.
1502.04860.extracted_bib
{b6 T. Koike-Akino, P. Popovski, and V. Tarokh, ��Optimized constellations for two-way wireless relaying with physical network coding,�� IEEE J. Sel. Areas Commun., vol. 27, no. 5, pp. 773�C787, Jun. 2009.
1502.04860.extracted_bib
{b7 S. Xu and Y. Hua, ��Optimal design of spatial source-and-relay matrices for a non-regenerative two-way MIMO relay system,�� IEEE Trans. Wireless Commun., vol. 5, no. 10, pp. 1645�C1655, May 2011.
1502.04860.extracted_bib
{b8 Y. Rong, ��Joint source and relay optimization for two-way linear nonregenerative MIMO relay communications,�� IEEE Trans. Signal Process., vol. 60, no. 12, pp. 6533�C6546, Dec. 2012.
1502.04860.extracted_bib
{b9 K.-J. Lee, H. Sung, E. Park, and I. Lee, ��Joint optimization for one and two-way MIMOAFmultiple-relay systems,�� IEEE Trans.Wireless Commun., vol. 9, no. 12, pp. 3671�C3681, Dec. 2010.
1502.04860.extracted_bib
{b10 H. Bolcskei, R. U. Nabar, O. Oyman, and A. J. Paulraj, ��Capacity scaling laws in MIMO relay networks," IEEE Trans. Wireless Commun., vol. 5, no. 6, pp. 1433-1444, June 2006.
1502.04860.extracted_bib
{b11 K.-J. Lee, K. W. Lee, H. Sung, and I. Lee, ��Sum-rate maximization for two-way MIMO amplify-and-forward relaying systems," in Proc. IEEE Veh. Technol. Conf., pp. 1-5, Apr. 2009.
1502.04860.extracted_bib
{b12 F. Roemer and M. Haardt, ��Algebraic norm-maximizing (ANOMAX) transmit strategy for two-way relaying with MIMO amplify and forward relays,�� IEEE Signal Process. Lett., vol. 16, no. 10, pp. 909�C912, Oct. 2009.
1502.04860.extracted_bib
{b13 T. Unger and A. Klein, ��Duplex schemes in multiple antenna two-hop relaying,�� EURASIP Journal on Advances in Signal Processing, vol. 2008,Article ID 128592, 2008.
1502.04860.extracted_bib
{b14 N. Jindal and Z. Luo, ��Capacity limits of multiple antenna multicast,�� in Proceedings of the IEEE International Symposium on Information Theory (ISIT ��06), pp. 1841�C1845, Seattle,Wash, USA, July 2006.
1502.04860.extracted_bib
{b15 R. Zhang, Y.-C. Liang, C. C. Chai, and S. Cui, ��Optimal beamforming for two-way multi-antenna relay channel with analogue network coding,�� IEEE Journal on Selected Areas in Communications, vol. 27, no. 5, pp. 699�C712, 2009.
1502.04860.extracted_bib
{b16 R. Wang and M. Tao, ��Joint source and relay precoding designs for MIMO two-way relaying based on MSE criterion,�� IEEE. Trans. Signal Process., vol. 60, no. 3, pp. 1352�C1365, Mar. 2012.
1502.04860.extracted_bib
{f1 H. Wang, L. Li, L. Song, and X. Gao, ��A linear precoding scheme for downlink multiuser MIMO precoding systems,�� IEEE Commun. Lett., vol. 15, no. 6, pp. 653�C655, June 2011.
1502.04860.extracted_bib
{a1 H. Park, H. J. Yang, J. Chun and R. Adve, ��A closed-form power allocation and signal alignment for a diagonalized MIMO two-way relay channel with linear receivers,�� IEEE Trans. Signal Process., vol. 60, no. 11, pp. 5948-5962, Nov. 2012.
1502.04860.extracted_bib
{b17 Veria Havary-Nassab, Shahram Shahbazpanahi and Ali Grami, "Optimal Distributed Beamforming for Two-Way Relay Networks," IEEE Trans. Signal Process., vol. 58, no. 3, pp. 1238-1250, Mar. 2010.
1502.04860.extracted_bib
{a5 Pramod Viswanath, Venkat Anantharam and David N. C. Tse, "Optimal Sequences, Power Control, and User Capacity of Synchronous CDMA Systems with Linear MMSE Multiuser Receivers," IEEE Trans. Inf. Theory, vol. 45, no. 6, pp. 1968�C1983, 1999.
1502.04860.extracted_bib
{a6 Rui Zhang, Ying-Chang Liang, Chin Choy Chai and Shuguang Cui, " Optimal Beamforming for Two-Way Multi-Antenna Relay Channel with Analogue Network Coding," IEEE J. Sel. Areas Commun., vol. 27, no. 5, pp. 699 - 712, June 2009
1502.04860.extracted_bib
{a4 Cristian Lameiro, Javier V��a and Ignacio Santamar��a, "Amplify-and-Forward Strategies in the Two-Way Relay Channel With Analog Tx�CRx Beamforming," IEEE Trans. Veh. Technol., vol. 62, no.2, pp. 642-654, Feb. 2013.
1502.04860.extracted_bib
{a2 Yue Rong, "Simplified Algorithms for Optimizing Multiuser Multi-Hop MIMO Relay Systems," IEEE Trans. Commun. vol. 59, no. 10, pp. 2896-2904, Oct. 2011.
1502.04860.extracted_bib
{a7 S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge, U.K.: Cambridge Univ. Press, 2004.
1502.04860.extracted_bib
{c1 X. Gao, L. Li, P. Zhang, L. Song, G. Zhang, Q. Wang and H. Tian, "Low-complexity RBD precoding method for downlink of multi-user MIMO system," Electronics Letters, wol. 46, no. 23, pp. 1574 - 1576, Nov. 2010.
1502.04860.extracted_bib
{a3 G. H. Golub and C. F. Van Loan, Matrix Computations, Johns Hopkins Studies in the Mathematical Sciences, Johns Hopkins University Press, Baltimore, Md, USA, 3rd edition, 1996.
1502.04860.extracted_bib
{b18 H. Wang, L. Li, L. Song, and X. Gao, ��A linear precoding scheme for downlink multiuser mimo precoding systems,�� IEEE Communications Letters, vol. 15, no. 6, pp. 653�C655, 2011.
1502.04860.extracted_bib
{b19 R. Mo and Y. H. Chew, ��Precoder design for non-regenerative MIMO relay systems,�� IEEE Transactions on Wireless Communications, vol. 8, no. 10, pp. 5041�C5049, 2009.
1502.04860.extracted_bib
{d1 S. Vishwanath, N. Jindal, and A. Goldsmith, ��On the capacity of multiple input multiple output broadcast channels,�� in Proceedings of the International Conference on Communications (ICC ��02), pp. 1444�C1450, May 2002.
1502.04860.extracted_bib
{e1 C. B. Peel, B. M. Hochwald, and A. L. Swindlehurst, ��A vectorperturbation technique for near capacity multiantenna multiuser communication��part I: channel inversion and regularization,�� IEEE Trans. Commun., vol. 52, no. 1, pp. 195��202, Jan. 2005.