\section{Introduction}
Massive MIMO is one of the key technologies used to boost capacity in currently under deployment 5G wireless networks \cite{redbook}; indeed, leveraging the joint coherent transmission/reception from several antenna elements, and the fully digital beamforming at the base station (BS), the massive MIMO technology permits an aggressive spatial multiplexing of a large number of user equipments (UEs) in the same time-frequency slot.

Recently, RISs have emerged as one of the most striking innovations for the evolution of 5G systems into 6G systems \cite{basar2019wireless,RuiZhang_TWC_2019}. A RIS is (usually) a passive planar structure, made of several sub-lambda spaced reflective elements, whose electromagnetic characteristics and in particular the phase offsets imposed when reflecting impinging waves can be controlled and adaptively changed. RISs can be thus employed to extend coverage, and/or to improve the signal-to-interference plus noise ratio (SINR) at the intended receiver location. Several studies have recently appeared showing the benefits that RISs can bring in different applications such as mobile edge computing networks\cite{Bai_MEC_JSAC2020}, physical layer security systems\cite{Hong_PLS_TCOM2020}, cognitive radio networks\cite{Zhang_Cognitive_TVT2020} and radar systems\cite{buzzi2021radar}.

In this paper, we consider the case in which a RIS is placed at short distance and in front of an active antenna array (see Fig. 1), with a non-large number of elements. The UEs to be served are placed in the backside of the active antenna array. The structure is to be intended as a low-complexity approximation of a massive MIMO system. A similar structure, but in a different setting and under different hypothesis, was also considered in \cite{jamali2020intelligent}. In the considered RIS-aided massive MIMO architecture, each active antenna is equipped with a dedicated RF chain and illuminates the RIS. Each passive element at the RIS receives a superposition of the signals transmitted by the active antennas and adds a desired phase shift to the overall signal. Leveraging the additional system optimization parameters introduced by the tunable reflecting elements of the RIS, a simpler transceiver architecture (i.e., less antennas and less RF chains) can be used than that usually employed in conventional massive MIMO systems with large-scale antenna array.

We firstly provide a signal and channel model suitable for the analysis of such a system, then we propose a channel estimation procedure which, exploiting only the active antennas and using different configurations of the RIS phase shifts, is capable to estimate the channel from the large number of RIS elements to all the UEs.  We derive a closed-form expression for an achievable downlink spectral efficiency per user, by using the popular hardening lower-bound, and formulate a generic optimization problem which can be used for two purposes: i) the minimization of the cross-correlation among channel vectors to reduce the multiuser interference, and ii) striking a balance between the conflicting requirements that the RIS phases should not lower the norms of the composite channels, while they should minimize the cross-correlation among channel vectors. In the numerical results, we show the effectiveness of the considered architecture assuming both omnidirectional and directional active antennas which allow us to overcome the traditional system with only active antennas.
\section{System Model}

Consider a single-cell system where a base station is equipped with a RIS-aided antenna array, and serves $K$ single-antenna UEs. We denote by $N_A$ the elements of the active antenna array, and by $N_R> N_A$ the number of configurable reflective elements of the RIS. The uplink channel between the $k$-th UE and the RIS is represented by the $N_R$-dimensional vector ${\mathbf{h}}_k$, while the matrix-valued channel from the RIS to the receive antenna array is represented through the $(N_A \times N_R)$-dimensional matrix $\mathbf{H}$.
The tunable phase shifts introduced by the reflective elements of the RIS are represented through a diagonal $(N_R \times N_R)$-dimensional matrix $\boldsymbol{\Phi}$; in particular, the $\ell$-th coefficient on the diagonal of $\boldsymbol{\Phi}$ is $e^{\i \phi_\ell}$, with $\phi_\ell$ the phase shift associated to the $\ell$-th element of the RIS.
Accordingly, the composite uplink channel from the $k$-th UE to the active array is the $N_A$-dimensional vector $\bar{\mathbf{h}}_k= \mathbf{H} \boldsymbol{\Phi} \mathbf{h}_k$.

Under the rich-scattering assumption, usually verified at sub-6 GHz frequencies, the channel vector $\mathbf{h}_k$ is modeled as a complex Gaussian distributed vector with zero mean and covariance matrix $\beta_k \mathbf{I}_{N_R}$, with $\mathbf{I}_{N_R}$ denoting the identity matrix of order $N_R$ and $\beta_k$ representing the large-scale fading coefficient for the $k$-th UE.

\subsection{Channel Model}
We model the RIS-to-active array $(N_A \times N_R)$-dimensional matrix $\mathbf{H}$. To this aim, and for the sake of simplicity, we consider the 2D scheme represented in Fig. \ref{fig:antenna}, where we have an active array that illuminates a large RIS placed at a short distance, that we denote by $D$; the UEs to be served are placed in the backside of the active antenna array. We make the following assumptions:
\begin{itemize}
	\item[i)] We neglect the blockage that the active array can introduce on the electromagnetic radiation reflected by the RIS. In practical 3D scenarios, the active array can be placed laterally with respect to the RIS so as to simply avoid the blockage problem.
	\item[ii)] The direct signal from the active antenna array to the mobile UEs will be neglected because it is considerably weaker than that from the RIS-reflected signal.
	\item[iii)] Given the short distance between the RIS and the active antenna array, the channel $\mathbf{H}$ is modeled as a deterministic quantity.
	\item[iv)] We neglect the mutual coupling between the elements of the active array and between the elements of the RIS. This is an assumption that is usually verified provided that the element spacing at the active array and at the RIS, i.e. $d_A$ and $d_R$, respectively, does not fall below the typical $\lambda/2$ value \cite{qian2021mutual}, with $\lambda$ being the wavelength corresponding to the radiated frequency.
	\item[v)] We assume that each element of the RIS is in the far field of each element of the active antenna array. However,  the whole RIS is not required to be in the far field of the whole active antenna array.
\end{itemize}
With regard to assumption v) we note that the far field condition among the generic $i$-th transmit antenna and $j$-th RIS element is $d_{i,j}>2 \max\{ \Delta_A^2, \Delta_R^2\}/\lambda$ and
$d_{i,j}\gg \max\{ \Delta_A, \Delta_R, \lambda\}$, with $\Delta_A$ and $\Delta_R$ the size of the active antenna and RIS element, respectively \cite{stutzman2012antenna}. The far-field condition among the whole RIS and the whole active array would have instead required that the fulfillement of the previous conditions with $d_{i,j}$,  $\Delta_A$ and $\Delta_R$ replaced by $D$, $N_A d_A$ and $N_R d_R$. The absence of the far field assumption between the whole active array and the RIS, although making the matrix $\mathbf{H}$ modelling somewhat more involved, is crucial in order to have a full-rank matrix channel that will permit communicating along many signal space directions.

Based on the above assumptions, and according to standard electromagnetic radiation theory, the $(i,j)$-th entry of $\mathbf{H}$, say $\mathbf{H}(i,j)$, representing the uplink channel coefficient from the $j$-th element of the RIS to the $i$-th active antenna, can be written as:
\begin{equation}
	\mathbf{H}(i,j)=\displaystyle \sqrt{\rho G_A(\theta_{i,j})G_R(\theta_{i,j})} \frac{\lambda}{4 \pi d_{i,j}} e^{-\i 2\pi d_{i,j}/\lambda} \; ,
	\label{eq:H_ij}
\end{equation}
where $G_A(\theta_{i,j})$ and $G_R(\theta_{i,j})$ represent the $i$-th active antenna and the $j$-th passive RIS element gains corresponding to the look angle $\theta_{i,j}$ (see Fig. \ref{fig:antenna}), while $\rho\leq1$ is a real-valued term modeling the RIS efficiency in reflecting the impinging waves.

\begin{figure}[t]
	\begin{center}	\includegraphics[scale=0.3]{fig_array_geometry}
	\end{center}
	\caption{Considered antenna scheme. $\theta_{i,j}$ is the look angle of the $i$-th element of the active array with respect to the $j$-th element of the RIS.}
	\label{fig:antenna}
\end{figure}

\subsection{Active antenna configuration}

In order to see the impact of the active antenna configuration, we consider two kinds of active antennas: \textit{omnidirectional} and \textit{directional}. In the first case, the active antennas are omnidirectional and radiate power in all the angular directions, i.e., $G_A(\theta)=3$ dB $\forall \theta \in [0, \pi]$. In the latter case, the active antennas radiate power only in one angular sector, i.e., $G_A(\theta)=g_{\text{dir}}$, $\forall \theta \in [-\alpha, \alpha]$. Note that, given the RIS and active array relative position depicted in Fig. \ref{fig:antenna}, if $N_R \gg N_A$ some RIS elements are external to the angular sector in which each active antenna radiate power. This situation strongly degrades the performance of the system because, actually, the ``effective'' elements of the RIS are only the ones that receive power from active antennas. To avoid this undesired behaviour, we can increase the active antenna spacing $d_A$ according to

\begin{equation}
	d_A=\displaystyle \frac{(N_R-1)d_R-2D\tan{\alpha}}{N_A - 2} \,
	\label{spacing_directional}
\end{equation}
where the above relationship can be obtained with simple geometric formulations and the notation refers to Fig. \ref{fig:antenna}.

\section{Uplink/Downlink signal processing}

\subsection{Uplink channel estimation}
During the uplink training phase, each UE transmits a known pilot sequence; we denote by $\boldsymbol{\varphi}_k$ the pilot sequence, of length $\tau_p$ and squared norm $\tau_p$, assigned to the $k$-th UE. Since, in this paper we assume $N_R>N_A$, in order to allow the active array to estimate the channels $\mathbf{h}_1, \mathbf{h}_K$, we use a similar approach as in \cite{Buzzi_TCCN_2021}. In particular, we assume that the pilot sequences are transmitted by the UEs $Q$ times, each one with a different RIS configuration. Without loss of generality, we assume $Q=N_R/N_A$ in order to guarantee generation of a number of observables that is not smaller than the number of unknown coefficients.  Pilots are drawn from a set of orthogonal sequences, which, for the sake of simplicity, are assumed to be real.
Denoting by $\eta_k^{(u,t)}$ the uplink transmit power from the $k$-th UE during the training phase, the data observed at the active antenna array when the RIS assumes the $q$-th configuration can be arranged in the following $(N_A \times \tau_p)$-dimensional matrix:
\begin{equation}
	\mathbf{Y}^{(q)}=\displaystyle \sum_{k=1}^K \sqrt{\eta_k^{(u,t)}}\mathbf{H} \boldsymbol{\Phi}_T^{(q)} \mathbf{h}_k \boldsymbol{\varphi}_k^T + \mathbf{N}^{(q)} \; ,
\end{equation}
where $\boldsymbol{\Phi}_T^{(q)}$ is the $q$-th configuration of the RIS assumed
$
	\boldsymbol{\Phi}_T^{(q)}=\text{diag}\left(\boldsymbol{\phi}^{(q)}\right),
$
and $\boldsymbol{\phi}^{(q)}$ is a $N_R$-dimensional vector whose $\ell$-th entry is $e^{\i\widetilde{\phi}_{\ell,q}}$ with $\widetilde{\phi}_{\ell,q}$ a random angle chosen from a quantized set of phase shifts available at the RIS.
In order to estimate the $k$-th UE channel $\mathbf{h}_k$, the data $\mathbf{Y}^{(q)}$ is projected along the $k$-th UE pilot sequence, leading to
\begin{equation}
	\mathbf{y}_{k,q}=\sqrt{\eta_k^{(u,t)}}\tau_p \mathbf{H}\boldsymbol{\Phi}_T^{(q)} \mathbf{h}_k + \sum_{\substack{j=1\\ j \neq k}}^K
	\sqrt{\eta_j^{(u,t)}} \mathbf{H} \boldsymbol{\Phi}_T^{(q)} \mathbf{h}_j \rho_{j,k} + \mathbf{n}_{k,q} \; ,
	\label{eq:y_k}
\end{equation}
with $\rho_{j,k}=\boldsymbol{\varphi}_j^T \boldsymbol{\varphi}_k$, and $\mathbf{n}_{k,q} = \mathbf{N}^{(q)} \boldsymbol{\varphi}_k \sim {\cal CN}(0,\sigma^2 \tau_p \mathbf{I}_{N_A})$. These observables are collected for all $q=1,\ldots,Q$ as the following $N_R \times 1$-dimensional vector
$
	\widetilde{\mathbf{y}}_k=\left[ \mathbf{y}_{k,1}^T, \; \ldots , \; \mathbf{y}_{k,Q}^T\right]^T.
$
After some manipulations, $\widetilde{\mathbf{y}}_k$ can be written as
\begin{equation}
	\widetilde{\mathbf{y}}_k=\sqrt{\eta_k^{(u,t)}}\tau_p \widetilde{\mathbf{H}}_T \mathbf{h}_k + \sum_{\substack{j=1\\ j \neq k}}^K
	\sqrt{\eta_j^{(u,t)}} \widetilde{\mathbf{H}}_T \mathbf{h}_j \rho_{j,k} + \widetilde{\mathbf{n}}_k \; ,
	\label{y_k_tilde2}
\end{equation}
where $\widetilde{\mathbf{n}}_k=\left[ \mathbf{n}_{k,1}^T, \; \ldots , \; \mathbf{n}_{k,Q}^T\right]^T$ and
\begin{equation}
	\widetilde{\mathbf{H}}_T=\left[ \left(\mathbf{H}\boldsymbol{\Phi}_T^{(1)}\right)^T, \; \ldots , \; \left(\mathbf{H}\boldsymbol{\Phi}_T^{(Q)}\right)^T\right]^T.
	\label{H_tilde_T}
\end{equation}

In order to estimate the \textit{essential directions} of the channels $\left\lbrace \mathbf{h}_1, \ldots, \mathbf{h}_K  \right\rbrace$, we use the economy size singular valued decomposition (SVD) of the matrix $\widetilde{\mathbf{H}}_T$
\begin{equation}
	\widetilde{\mathbf{H}}_T=\mathbf{U}\bm{\Lambda}\mathbf{V}^H \approx \widetilde{\mathbf{U}}\widetilde{\bm{\Lambda}}\widetilde{\mathbf{V}}^H,
	\label{H_tilde_T_svd}
\end{equation}
where $\mathbf{U},\bm{\Lambda},\mathbf{V}$ are the $N_R$-dimensional matrices corresponding to the regular SVD and, say $q$ the  number of significant eigenvalues of the matrix $\widetilde{\mathbf{H}}_T$, $\widetilde{\mathbf{U}}$ and $\widetilde{\mathbf{V}}$ contains the columns of $\mathbf{U}$ and $\mathbf{V}$ corresponding to the $q$ highest eigenvalues and $\widetilde{\bm{\Lambda}}$ is the diagonal matrix containing these eigenvalues. Note that the approximation in Eq. \eqref{H_tilde_T_svd} holds if the matrix $\widetilde{\mathbf{H}}_T$ has $q$ eigenvalues significantly higher with respect to the remaining ones. In this paper, it holds because of the structure of the matrix $\mathbf{H}$. Given the properties of the SVD, we have
$
\mathbf{h}_k= \mathbf{V} \mathbf{V}^H \mathbf{h}_k = \mathbf{V} \mathbf{v}_k,
$
where $\mathbf{v}_k=\mathbf{V}^H \mathbf{h}_k= \left[ \widetilde{\mathbf{v}}_k^T, \overline{\mathbf{v}}_k^T\right]^T$ and using the approximation in Eq. \eqref{H_tilde_T_svd} the following holds
\begin{equation}
	\mathbf{h}_k \approx \widetilde{\mathbf{V}} \widetilde{\mathbf{v}}_k.
	\label{v_k_approx}
\end{equation}
We perform a channel estimation technique aimed at determining the strongest directions of the channels, corresponding to the highest eigenvalues, which we can estimate with highest power, and thereby recover the channels exploiting Eq. \eqref{v_k_approx}. Using the above considerations, we focus on the $q$-dimensional vector
\begin{equation}
	\begin{array}{lll}
		\overline{\mathbf{y}}_k= &\!\!\!\!\!\widetilde{\mathbf{U}}^H\widetilde{\mathbf{y}}_k\!\!\approx\!\!\sqrt{\eta_k^{(u,t)}}\tau_p \widetilde{\bm{\Lambda}}\widetilde{\mathbf{v}}_k\!\! +\!\!\!  \displaystyle \sum_{\substack{j=1\\ j \neq k}}^K
		\!\!\sqrt{\eta_j^{(u,t)}} \widetilde{\bm{\Lambda}}\widetilde{\mathbf{v}}_j\rho_{j,k} \!+\! \widetilde{\mathbf{U}}^H\widetilde{\mathbf{n}}_k \; .
	\end{array}
	\label{y_k_bar2}
\end{equation}
Thus, the linear MMSE estimate of $\widetilde{\mathbf{v}}_k$ can be expressed as
$
	\widehat{\widetilde{\mathbf{v}}}_k=\mathbf{R}_{\widetilde{\mathbf{v}}_k \overline{\mathbf{y}}_k}\mathbf{R}_{\overline{\mathbf{y}}_k \overline{\mathbf{y}}_k}^{-1}\overline{\mathbf{y}}_k \; ,
$
with
\begin{equation}
	\mathbf{R}_{\widetilde{\mathbf{v}}_k \overline{\mathbf{y}}_k} \triangleq E\left\{ \widetilde{\mathbf{v}}_k \overline{\mathbf{y}}_k^H \right\}=\sqrt{\eta_k^{(u,t)}} \tau_p \beta_k \widetilde{\bm{\Lambda}}^H,
\end{equation}
\begin{equation}
	\mathbf{R}_{\overline{\mathbf{y}}_k \overline{\mathbf{y}}_k}\!\! \triangleq \!E\left\{ \overline{\mathbf{y}}_k \overline{\mathbf{y}}_k^H \right\}\!\!=\!\!\Bigg[ \displaystyle \sum_{j=1}^K \eta_j^{(u,t)} \beta_j  \rho_{j,k}^2\Bigg]\!  \widetilde{\bm{\Lambda}} \widetilde{\bm{\Lambda}}^H\!\! + \sigma^2 \tau_p \mathbf{I}_{q}.
\end{equation}
The channel estimate $\widehat{\mathbf{h}}_k$ can be thus written as
\begin{equation}
	\widehat{\mathbf{h}}_k \approx \widetilde{\mathbf{V}} \mathbf{R}_{\widetilde{\mathbf{v}}_k \overline{\mathbf{y}}_k}\mathbf{R}_{\overline{\mathbf{y}}_k \overline{\mathbf{y}}_k}^{-1}\overline{\mathbf{y}}_k.
	\label{v_k_approx_estimate}
\end{equation}
It is easy to realize that the LMMSE estimation of the vector containing the strongest directions of the channel $\widehat{\mathbf{h}}_k$ is Gaussian distributed with zero mean and covariance matrix
\begin{equation}
	\mathbf{R}_{\widehat{\mathbf{h}}_k \widehat{\mathbf{h}}_k}=\widetilde{\mathbf{V}}\mathbf{R}_{\widetilde{\mathbf{v}}_k \overline{\mathbf{y}}_k}\mathbf{R}_{\overline{\mathbf{y}}_k \overline{\mathbf{y}}_k}^{-1}
\mathbf{R}_{\widetilde{\mathbf{v}}_k \overline{\mathbf{y}}_k}^H \widetilde{\mathbf{V}}^H.
\end{equation}

\subsection{Downlink data transmission}
During the downlink data transmission phase, the complex envelope of the signal received by the $k$-th UE in the $n$-th symbol interval can be shown to be expressed as
\begin{equation}
	\begin{array}{lll}
		r_k(n)=&\sqrt{\eta_k^{(d)}} \mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T \mathbf{w}_k x_k(n)  \\ &\displaystyle +\sum_{\substack{j=1\\ j\neq k}}^K
		\sqrt{\eta_j^{(d)}} \mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T \mathbf{w}_j x_j(n) + z_k(n) \, ,
	\end{array} \label{eq:downlink_data}
\end{equation}
where $\eta_k^{(d)}$ and $\mathbf{w}_k$ represent the downlink transmit power and the downlink beamformer for UE $k$, respectively, while $x_k(n)$ is the data symbol, drawn from a constellation with unit average energy, intended for UE $k$ in the $n$-th signaling interval. The term $z_k(n)$ is the AWGN term, modeled as ${\cal CN}(0, \sigma^2_k)$.

\section{Performance measures} \label{Performance_measures}
For benchmarking purposes, we start considering the ideal case in which perfect CSI is available both at the active array and at the UEs.
Given \eqref{eq:downlink_data}, the $k$-th UE downlink SINR, for perfect CSI, is written as
\begin{equation}
	\gamma_k^{(d)}= \displaystyle \frac{\eta_k^{(d)} \left| \mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T   \mathbf{w}_k\right|^2}{
		\displaystyle \sum_{\substack{j=1\\ j\neq k}}^K \eta_j^{(d)} 	\left| \mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T \mathbf{w}_j \right|^2 + \sigma^2_k  \; .
	}
	\label{eq:DL_SINR}
\end{equation}
Consequently, assuming Gaussian-distributed codewords, the spectral efficiency (SE) for the $k$-th UE  can be obtained by using the classical Shannon expression, i.e.:
\begin{equation}
	\text{SE}_k^{(d)}= \xi \log_2 (1+ \gamma_k^{(d)}) \; , \quad [\mbox{bit/s/Hz}] \; .
	\label{DL_SE}
\end{equation}
where $\xi$ accounts for the fraction of the channel coherence interval used for the downlink data transmission.
In general, when the channel realizations are not perfectly known, the UE cannot clearly tell what is ``signal'' and what is
``interference'' in Eq. \eqref{eq:downlink_data}. With imperfect CSI, the intuitive notion of SINR in Eq. \eqref{eq:DL_SINR} is in general not rigorously related to a corresponding notion of information theoretic achievable rate \cite{Caire_Bounds2018}.

\subsection{The upper-bound (UB) to the system performance}
An upper bound of the downlink SE, can be obtained by optimistically assuming perfect knowledge of the channel estimates at the UEs. This bound gives us a first sense of the system performance, and is given by~\cite{Caire_Bounds2018}
\begin{equation}
	\text{SE}_{k,\rm{UB}}^{(d)}= \bar{\xi}~E\left\{\log_2 (1+ \gamma_k^{(d)})\right\} \; , \quad [\mbox{bit/s/Hz}] \; ,
	\label{rate_Dl_UB}
\end{equation}
where the pre-log factor $\bar{\xi}<1$ accounts for the fraction of the channel coherence interval used for downlink data transmission, and for the pilot overhead. In Eq. \eqref{rate_Dl_UB} the expectation is computed over the small-scale fading quantities, and the beamforming vectors $\mathbf{w}_1,\ldots, \mathbf{w}_K$ in Eq. \eqref{eq:DL_SINR} are designed upon the channel estimates.

\subsection{The hardening lower-bound (LB) to system performance} \label{LB_performance}
To obtain an achievable SE expression in closed form, we derive the so-called \textit{hardening} lower-bound on the capacity, as detailed in~\cite{redbook}.
Under the assumption of Rayleigh fading, and using a conjugate beamformer, i.e. letting
$\mathbf{w}_k=\mathbf{H}^* \boldsymbol{\Phi}^* \widehat{\mathbf{h}}_k^*$, an achievable SE for the $k$-th UE can be obtained by treating all the interference and noise contributions as uncorrelated effective noise. Letting,
\begin{equation}
	\begin{array}{lll}
		&r_k(n)=\underbrace{\sqrt{\eta_k^{(d)}} E\left\{\mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T \mathbf{w}_k\right\}}_{\text{DS}_k} x_k(n)
		\\ & + \displaystyle \underbrace{\sqrt{\eta_k^{(d)}}\left[\mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T \mathbf{w}_k - E\left\{\mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T \mathbf{w}_k\right\}\right]}_{\text{BU}_k} x_k(n)
		\\ &\displaystyle +\sum_{\substack{j=1\\ j\neq k}}^K
		\underbrace{\sqrt{\eta_j^{(d)}} \mathbf{h}_k^T \boldsymbol{\Phi}^T \mathbf{H}^T \mathbf{w}_j}_{\text{UI}_{k,j}} x_j(n) + z_k(n) \, ,
	\end{array}
	\label{eq:downlink_data_hardening}
\end{equation}
then the hardening lower-bound is given by
\begin{equation}
	\text{SE}_{k,\rm{LB}}^{(d)} =  \bar{\xi} \log_2(1+ \gamma_{k, \text{LB}}^{(d)}) \, ,
	\label{eq:hardening}
\end{equation} \vspace*{-2mm}
with
\begin{equation}
	\gamma_{k, \text{LB}}^{(d)}= \displaystyle\frac{\left|\text{DS}_k\right|^2}{
		E\{\left|\text{BU}_k\right|^2\} + \displaystyle \sum_{\substack{j=1\\ j\neq k}}^K E
		\{\left| \text{UI}_{k,j} \right|^2 \}+ \sigma^2_k }.
\end{equation}

The closed-form expression of $\gamma_{k, \text{LB}}^{(d)}$ is reported in Eq. \eqref{eq:SINR_D_LB} at the top of next page, where $\mathcal{P}_k$ includes the indices of the UEs using the same pilot as user $k$.
Due to space constraints, we omit the proof including the derivation of~\eqref{eq:SINR_D_LB}.
\begin{figure*}[t]
	\begin{equation}
		\gamma_{k, \text{LB}}^{(d)}= \displaystyle \frac{{\eta_k^{(d)}} \text{tr}^2\left[   \mathbf{A}_{\boldsymbol{\Phi}} \mathbf{R}_{\widehat{\mathbf{h}}_k \widehat{\mathbf{h}}_k}^*\right]}
		{\left\lbrace\begin{array}{lr}
				\eta_k^{(d)} \text{tr}\left( \beta_k \mathbf{A}_{\boldsymbol{\Phi}} \mathbf{R}_{\widehat{\mathbf{h}}_k \widehat{\mathbf{h}}_k}^* \mathbf{A}_{\boldsymbol{\Phi}}\right) + \displaystyle \sum_{j \in {\cal P}_k} \eta_j^{(d)}\left(
				\displaystyle \frac{\beta_j \sqrt{\eta_j^{(u,t)}}}{\beta_k \sqrt{\eta_k^{(u,t)}}}\right)^2 \left[
				\text{tr}^2(\mathbf{A}_{\boldsymbol{\Phi}} \mathbf{R}_{\widehat{\mathbf{h}}_k \widehat{\mathbf{h}}_k}^*)+ \beta_k\text{tr}(\mathbf{A}_{\boldsymbol{\Phi}} \mathbf{R}_{\widehat{\mathbf{h}}_k \widehat{\mathbf{h}}_k}^* \mathbf{A}_{\boldsymbol{\Phi}})
				\right]    +  \\
				\hfill \displaystyle \sum_{j \notin {\cal P}_k}    \eta_j^{(d)} \beta_k\text{tr}\left[ \mathbf{A}_{\boldsymbol{\Phi}}
				\mathbf{R}_{\widehat{\mathbf{h}}_j \widehat{\mathbf{h}}_j}^*\mathbf{A}_{\boldsymbol{\Phi}}\right]
				+ \sigma^2_k \end{array}\right \rbrace}
	\label{eq:SINR_D_LB}
	\end{equation}
	\hrulefill
\end{figure*}

Finally, we can note that the analysis of the considered structure is, from a mathematical point of view, different from that required to analyze conventional massive MIMO systems with one active antenna array with a large number of elements.
As it is well known, in a conventional massive MIMO system, when the number of antenna elements grows without bound we can observe two phenomena: the \textit{channel hardening} (the channel norm tends to be deterministic) and the \textit{favorable propagation} (channel vectors of two different UEs become almost orthogonal).
In the case of RIS-aided massive MIMO, letting $N_R$ grow without bound and keeping $N_A$ fixed does not provide the above described effects. It can be shown that these properties hold \textit{only} if the number of active antennas increases without bound, we omit this proof here due to the lack of space.

\section{Optimization of the RIS phase shifts} \label{RIS_optimization}
In order to optimize the system performance for finite $N_A$, we define the matrix
$
\mathbf{Q}\left( \boldsymbol{\Phi}\right)=\mathbf{S}\left( \boldsymbol{\Phi}\right)\mathbf{S}\left( \boldsymbol{\Phi}\right)^H,
$
with
$
\mathbf{S}\left( \boldsymbol{\Phi}\right)=\mathbf{H} \boldsymbol{\Phi} \overline{\mathbf{H}},
$
and $\overline{\mathbf{H}}=\left[ \mathbf{h}_1, \ldots, \mathbf{h}_K\right ]$.
The entries of matrix $\mathbf{Q}\left( \boldsymbol{\Phi}\right)$ give us an idea on the potential interference that the UEs suffer for a given RIS configuration $\boldsymbol{\Phi}$.
We thus focus on the matrix $\mathbf{Q}\left( \boldsymbol{\Phi}\right)$ in order to formulate two problems aimed at the optimal configuration of the RIS phase offsets. We define two objective functions:
\begin{equation*}
f_1\left( \boldsymbol{\Phi}\right)= \displaystyle \sum_{k=1}^K \sum_{k'=k+1}^K \!\!\! |\bar{\mathbf{h}}_k^H \bar{\mathbf{h}}_{k'}|, \; \text{and} \; f_2\left( \boldsymbol{\Phi}\right)= \displaystyle \frac{f_1\left( \boldsymbol{\Phi}\right)}{\displaystyle \sum\nolimits_{k=1}^K \|\bar{\mathbf{h}}_k\|^2} \, ,
\end{equation*}
where
$
\bar{\mathbf{h}}_k=\mathbf{H} \boldsymbol{\Phi} \mathbf{h}_k \,
$
and the functional dependence on $\boldsymbol{\Phi}$ is omitted for ease of notation.
The minimization problem with objective function $f_1\left( \boldsymbol{\Phi}\right)$ is aimed at the minimization of the cross-correlation among channel vectors to reduce interference. The one with objective function $f_2\left( \boldsymbol{\Phi}\right)$ strikes a balance between the conflicting requirements that the RIS phases should not lower the norms of the composite channels, while at the same time should minimize cross-correlation among channel vectors.
The generic optimization problem can be written as
\begin{subequations}\label{Prob:RIS_opt}
\begin{align}
&\displaystyle\min_{\boldsymbol{\Phi}}\;  f_{(\cdot)}\left( \boldsymbol{\Phi}\right)\\
&\;\textrm{s.t.}\; |\phi_i|=1,\; i=1,\ldots, N_R
\end{align}
\end{subequations}
where $(\cdot) \, \in \left\lbrace 1,2 \right\rbrace$.
The optimization problem \eqref{Prob:RIS_opt} is not convex, hence we propose a solution based on the alternating maximization theory\cite{BertsekasNonLinear}.
We consider the $N_R$ optimization problems, $i=1,\ldots, N_R$,
\begin{subequations}\label{Prob:RIS_opt_ith}
	\begin{align}
		&\displaystyle\min_{\phi_i}\;  f_{(\cdot)}\left( \phi_i\right)\\
		&\;\textrm{s.t.}\; |\phi_i|=1,
	\end{align}
\end{subequations}
with $\phi_{\ell}, \ell=1,\ldots, i-1,i+1,\ldots, N_R$ remain fixed. This optimization problem can be solved using an exhaustive search. It is proved in \cite[Proposition 2.7.1]{BertsekasNonLinear} that iteratively solving \eqref{Prob:RIS_opt_ith} monotonically decreases the value of the objective of \eqref{Prob:RIS_opt}, and converges to a first-order optimal point if the solution of \eqref{Prob:RIS_opt_ith} is unique for any $i$. The generic optimization algorithm can be summarized as in Algorithm \ref{Alg:RIS_opt}.
\begin{algorithm}[t]
	\caption{Optimization of the RIS phase shifts}
	\begin{algorithmic}[1]
		\label{Alg:RIS_opt}
		\STATE Choose a starting random RIS configuration of the phase shifts $\boldsymbol{\Phi}$;
		\REPEAT
		\FOR{$i=1\to N_R$}
		\STATE Solve Problem \eqref{Prob:RIS_opt_ith};
		\STATE Set the $(i,i)$-th entry of the matrix $\boldsymbol{\Phi}$ as the solution of Problem \eqref{Prob:RIS_opt_ith}, $\phi_i^*$ say.
		\ENDFOR
		\UNTIL{convergence of the objective function}
	\end{algorithmic}
\end{algorithm}

\begin{figure*}[t]
	\begin{center}	\includegraphics[scale=0.41]{Fig_RIS_UatF_Directional}
	\end{center}
	\caption{CDF of the Spectral efficiency (SE) per user in the case of omnidirectional (Omni.) and directional (Dir.) active antennas with $N_R=64$ and $N_A=16$. Comparison between the case of perfect CSI (PCSI), and imperfect CSI (ICSI) with upper-bound (UB) and lower-bound (LB) of the performance.}
	\label{fig:RIS_bounds}
\end{figure*}

\section{Numerical Results}
In our simulation setup, we consider a communication bandwidth of $W = 20$ MHz centered over the carrier frequency
$f_0 = 1.9$ GHz. The additive thermal noise is assumed to have a power spectral density of $-174$ dBm/Hz,
while the front-end receiver at the active array and at the UEs is assumed to have a noise figure of 5 dB. The heights of active array and RIS are 10 meters while the heigh of the UEs is 1.5 meters. The $K=8$ UEs are uniformly distributed in an angular region of $120 \deg$ in front of the RIS at distances in $[10,400]$ meters. We model the large scale fading coefficients $\beta_k$ from the UEs to the RIS according to \cite{3GPP_38901_model} and the matrix $\mathbf{H}$ according to Eq. \eqref{eq:H_ij} with $D=5\lambda$, $\rho=1$, $d_A=\lambda/2$ and in the case of omnidirectional active antennas $G_R\left( \theta \right) = 3$ dB $\forall \theta$, $d_R=\lambda/2$; in the case of directional active antennas  $d_R$ follows Eq. \eqref{spacing_directional} with $g_{\text{dir}}$ and $\alpha$, specified in the figures.  With regard to the channel estimation procedure, $q$ is chosen as the number of eigenvalues which collects the 98\% of the total sum of eigenvalues. The pilot sequences have length $\tau_p = 16$, and the uplink transmit power for the pilot is assumed to be $\eta_k=800$ mW. The coherence interval is $\tau_c = 200$ samples long, and $\tau_c-\tau_p$ samples are used for downlink data. This implies $\bar{\xi}=1-\tau_p/\tau_c$. The maximum transmit power of the active array is $P_{\rm max}^{{\rm A}}=7$ dBW and it is equally distributed between the UEs, i.e., $\eta_k^{(d)}= P_{\rm max}^{{\rm A}}/K \norm{\mathbf{w}_k}^2$.

\begin{figure*}[t]
	\begin{center}	\includegraphics[scale=0.41]{Fig_RIS_optimization}
	\end{center}
	\caption{CDF of the Spectral efficiency (SE) per user in the case of omnidirectional (Omni.) and directional (Dir.) active antennas with $N_R=64$ and $N_A=16$. Comparison between the case of perfect CSI (PCSI), and imperfect CSI (ICSI) with upper-bound (UB) of the performance with RIS phase shifts optimization with objective functions $f_1(\cdot)$ and $f_2(\cdot)$.}
	\label{fig:RIS_optimization}
\end{figure*}

In Fig. \ref{fig:RIS_bounds} we report the cumulative distribution function (CDF) of the spectral efficiency (SE) per user using the expressions discussed in Section \ref{Performance_measures} for perfect CSI (PCSI) and imperfect CSI (ICSI). In this figure, we assume a random configuration of the RIS phase shifts. We can see that in the case of directional active antennas the performance in terms of SE increases with respect to the omnidirectional case due to the higher values in the matrix $\mathbf{H}$ which models the communication between the RIS and the active array. With the considered system, we can note that the hardening LB is not as tight as in the case of classical massive MIMO because of the reasons discussed in Section \ref{LB_performance}.
In Fig. \ref{fig:RIS_optimization} we report the CDF of the SE per user in the case of PCSI and ICSI using the UB of the performance. We report the cases of RIS with random configuration of phase shifts (Random RIS) and in the cases of phase shift optimization (Opt. RIS) detailed in Section \ref{RIS_optimization} where we denote by Obj. 1 the Problem \eqref{Prob:RIS_opt} with objective function $f_1(\cdot)$ and Obj. 2 the one with objective function $f_2(\cdot)$. We compare the performance with the traditional Massive MIMO with $N_A$ antennas. We can see that the optimized RIS with directional active antennas strongly increases the performance of the system compared with the traditional massive MIMO, especially in the PCSI case. Additionally, we can also see that the proposed channel estimation procedure is effective and that the phase shift optimization works well also in this case, especially when directional active antennas are assumed.

\section{Conclusions}
We considered a RIS-aided massive MIMO system in which a RIS is placed at short distance and in front of an active antenna array, with a non-large number of elements. The UEs to be served are placed in the backside of the active antenna array. We provided a signal and channel model suitable for the analysis of such a system and proposed a channel estimation procedure which, using different configurations of the RIS phase shifts, is capable to estimate the channel from the large number of RIS elements and all the UEs exploiting only the active antennas. We derived a closed-form expression for for an achievable downlink SE per user, by using the popular hardening lower-bound, and formulated a generic optimization problem which can be used for two purposes. In the numerical results, we showed the effectiveness of the considered architecture assuming both omnidirectional and directional active antennas, and of the phase shifts optimal configuration algorithm.

\bibliographystyle{IEEEtran}


\section{Introduction}
\label{intro}
It is known that the classical least square regression models achieve the optimal efficiency when the noises are Gaussian, however, they always underperform if the data is contaminated by non-Gaussian noises or outliers. Some robust regression models have been well developed in the past decades such as the median regression, the modal regression, the Huber regression and the least trimmed squares regression, etc.
Moreover, a new robust regression model named the maximum correntropy criterion regression (MCCR) has been theoretically studied within the frame of statistical learning in \cite{2015Learning}. Correntropy is constructed based on a kernel function and it is a generalized similarity measure between two random variables (see \cite{2006Generalized,gunduz2009correntropy,2007Correntropy,He2011Robust,2012Maximum} and \cite{2016Correntropy}). Correntropy induced loss function $\sigma^{2}(1-e^{-\frac{t^2}{\sigma^2}})$ and its scale parameter $\sigma$ are used to control the trade-off between the convergence rate and the robustness of the estimated regression function in MCCR models. While the scale parameter goes to infinity, the convergence rate of estimator has been established under some moment conditions or some absolute bounded conditions on the response variable (see Theorem 4 and 5 in \cite{2015Learning}), while the scale parameter is a fixed constant large enough, the convergence rate has been presented under the symmetric and uniformly bounded noises (see Theorem 6 in \cite{2015Learning}). However, the case that the scale parameter goes to zero has not been considered.

Some other works of MCCR have been presented. The model with the mixed symmetric stable noises and a fixed scale parameter was discussed in \cite{2018Learning}. The situation that the response variable with conditions of \begin{math} (1+\varepsilon)\end{math}-moment and the scale parameter greater than 1 was studied in \cite{feng2020learning}. Particularly, \cite{feng2020statistical} declared that the model with the correntropy induced loss function $\sigma^{-1}(1-e^{-\frac{t^2}{\sigma^2}})$ and a tending-to-zero scale parameter was modal regression. It should be noted that the loss functions for MCCR and the modal regression are essentially different when the scale parameters go to zero.

MCCR models with small or tending-to-zero scale parameters are worth to be considered. First, the scale parameters of MCCR play the same roles as the bandwidths of kernel density estimations, when the sample size $n$ is large or where the sample points are dense, the small or tending-to-zero scale parameters are expected to decrease the fitting error. Second, the scale parameters are similar to the smoothing parameters of local polynomial estimations, the small ones are needed where the regression functions have fluctuated. Last, \cite{2015Learning} has noted that the smaller scale parameters provide MCCR models with better robustness, it is meaningful to study the consistency of the estimates with nice robustness.

In this paper, we first give the theoretical study on MCCR models with tending-to-0 scale parameters. The convergence rates of estimates are discussed under the assumption of symmetric noises when the sample size goes to infinity. Then, the simulations and applications are presented to compare MCCR, Huber, and least square (LS) regression models as the sample size is finite. The rest of this paper is organized as follows. In Section 2, we give a short introduction to MCCR. The convergence rates of MCCR estimates with tending to 0 scale parameters under mixed symmetric noise are presented in Section 3. Section 4 and 5 display the numerical experiments and applications. The last is Section 6 for some conclusions.

\section{ Definition of MCCR}
We first assume that the regression model is given as
\begin{equation}
Y=f^{*}(X)+\varepsilon,
\end{equation}
where $X$ is the independent variable takes values in a compact metric space $\mathcal{X}$, the dependent variable $Y\in\mathcal{Y}=\mathbb{R}$, and $\varepsilon$ is the noise subject to $\mathbb{E}(\varepsilon|X)=0$. The purpose of a regression model is to learn the unknown regression function $f^*$ based on the given observations $z=\left\{(x_i,y_i)\right\}_{i=1}^n$, which are  independently sampled from the distribution $\rho$ of $(X, Y)$ on  $\mathcal{X} \times \mathcal{Y}$.

Given two scalar random variables $U$ and $V$, the correntropy $\mathcal{V}_\sigma$ between $U$ and $V$ is defined as $\mathcal{V}_{\sigma}(U,V)=\mathbb{E}{\mathcal{K}_{\sigma}(U,V)}$, where $\mathcal{K}_{\sigma}(\cdot,\cdot)$ is the Gaussian kernel function with the scale parameter $\sigma>0$. Theoretically, the idea of MCCR is to infer $f^*$  based on the maximization of $\mathcal{V}_\sigma$. Empirically, MCCR obtains the estimates of $f^*$ by maximizing the mean version of $\mathcal{V}_{\sigma}$ as follows

$$\hat{\mathcal{V}}_{\sigma,z}(f)=\frac{1}{n}{\sum_{i=1}^{n}{\mathcal{K}}_\sigma(y_i,f(x_i))}.$$
Let the output function of MCCR as
$$f_z=\arg\mathop{\max}\limits_{f\in \mathcal{H}}\hat{\mathcal{V}}_{\sigma,z}(f),$$
where the hypothesis space $\mathcal{H}$ is assumed to be a compact subset of a continuous functions space.
Furthermore, let the correntropy induced loss function $l_\sigma:\mathbb{R} \to [0,+\infty)$  as (see Definition 1 in \cite{2015Learning})
\begin{equation}
l_\sigma(t)=\sigma^{2}(1-e^{-\frac{t^2}{\sigma^2}}),
\end{equation}
then the estimator $f_z$ of MCCR can be equivalently expressed as
\begin{equation}
f_z=\arg\mathop{\min}\limits_{f\in \mathcal{H}}\frac{1}{n}{\sum_{i=1}^{n}l_\sigma(y_i-f(x_i))}.
\end{equation}
The discussions on the loss functions in \cite{2015Learning} and \cite{2013Robust} show that the scale parameter $\sigma$ tunes the tradeoff between robustness and convergence rate of MCCR model, a decreasing $\sigma$ enhances the robustness and reduces the convergence rate and vice versa.

\section{MCCR with mixture of symmetric noises}
In this section, we give the theoretical research on MCCR with a mixture of symmetric noises.
\subsection{Definitions and assumptions}
First, let $\rho_\mathcal{X}$, $\rho_\mathcal{Y}$ be the marginal distributions of $X, Y$ respectively. For any $f \in \mathcal{H}$, the empirical error in (3) is defined as
$${\varepsilon}_{z}^{\sigma}(f)=\frac{1}{n}{\sum_{i=1}^{n}l_\sigma(y_i-f(x_i))},$$
and its population version ${\varepsilon}^{\sigma}(f)$ is defined as
$${\varepsilon}^{\sigma}(f)=\int\limits_{\mathcal{X} \times \mathcal{Y}}l_{\sigma}(y-f(x))d{\rho}.$$
Let the distance between $f$ and $f^*$ under $L_{\rho_\mathcal{X}}^2$ be $\Vert f-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2$. Moreover, let notation $a_1 \lesssim a_2$ for $a_1, a_2\in\mathbb{R} $ means that there exists a positive constant $a_3$ such that $a_1 \leq a_{3}a_2$.
\begin{definition}[Mixture of symmetric distributions]
The univariate distribution function $P(t)$ (density function $p(t)$) is a mixed symmetric distribution if it is a convex combination of symmetric distributions $P_{i}(t)$ (density function $p_{i}(t)$), where $i=1, 2, \cdots, K$ and $K$ is a positive integer, that is, there exists $w_1, w_2,\cdots, w_K$ with $w_{i}>0, i=1, 2, \cdots, K$ and $\sum_{i=1}^{K}w_i=1$, such that for any $t\in \mathbb{R}$, there are
$$P(t)={\sum_{i=1}^{K}w_{i}P_{i}(t)},\quad p(t)={\sum_{i=1}^{K}w_{i}p_{i}(t)}.$$
\end{definition}
\begin{definition}[$l^{2}$-empirical covering number]
Let $x=\left\{x_1,x_2,\cdots,x_n\right\}\subset {\mathcal{X}}^n$ and  $\eta>0$. The $l^{2}$-empirical covering number $\mathcal{N}_{2}(\mathcal{H},\eta)$ of the hypothesis space $\mathcal{H}$ with radius $\eta$ is defined as follows (See \cite{anthony2009neural, guo2013concentration} and \cite{2007Learnability})
\[
\begin{aligned}
\mathcal{N}_{2}(\mathcal{H},\eta)&:=\sup\limits_{n\in \mathbb{N}} \sup\limits_{x\in {\chi}^n} \inf \left\{ l\in \mathbb{N}:\exists \left\{f_i\right\}_{i=1}^l \subset \mathcal{H}\, \mbox{such that for each}\, f\in \mathcal{H}\,\mbox{there exists some}\right.\\
&\left. i\in \left\{i=1,2,\cdots,l\right\}\, \mbox{with}\; \frac{1}{n}\sum_{j=1}^{n}\vert{f(x_{j})-f_{i}(x_j)\vert}^2\leq {\eta}^2 \right\}.
\end{aligned}
\]
\end{definition}
\begin{assumption}[Mixture of symmetric noises]
The noise $\varepsilon$ is a mixture of symmetric noises with 0 as the center, that is, for any $x\in \mathcal{X}$, the density $p_{\varepsilon|x}$ of the noise  $\varepsilon$ holds the following forms
$$p_{\varepsilon|x}(t)={\sum_{i=1}^{K}w_{i}p_{\varepsilon,i|x}(t)},\quad t\in \mathbb{R},$$
where $K$ is a positive integer, $w_{i}>0,i=1,2,\cdots,K$ with $\sum_{i=1}^{K}w_i=1$, and $p_{\varepsilon,i|x}$  is the density function of a symmetric distribution $P_{\varepsilon,i|x}$  with 0 as the center.
\end{assumption}

\begin{assumption}
 For any $x\in \mathcal{X}$ and $i=1,2,\cdots,K$, the Fourier transform $\widehat{p_{\varepsilon,i|x}}$  of $p_{\varepsilon,i|x}$  is positive, and there exist two positive constants $c_{0,i},C_{0,i}$, such that $\widehat{p_{\varepsilon,i|x}}(\xi)\geq C_{0,i}$ for $\xi \in [-c_{0,i},c_{0,i}]$.
\end{assumption}

\begin{assumption}[Complexity Assumption]
There exist positive constants $s$ and $c$ with $0<s<2$, such that
$$\log \mathcal{N}_2(\mathcal{H},\eta)\leq {c\eta}^{-s}, \forall \eta >0.$$
\end{assumption}

Moreover, we assume that $\sup\limits_{f\in \mathcal{H}}\Vert f\Vert_{\infty}\leq M$ and $\Vert f^*\Vert_{\infty}\leq M$, where $M$ is a positive constant. Assumptions 1 and 2 allow the noises to be kinds of random variables, such as Gaussian distribution, Cauchy distribution, Laplace distribution, Linnik distribution, symmetric stable distribution and the convex combinations of the above distributions (See \cite{1998Tails, fama1968some} and \cite{miller1978properties}). Assumption 3 is a common representation  for the complexity of  hypothesis space in statistical learning (See \cite{zhou2002covering} and \cite{cucker2007learning}).

\subsection{ Theoretical results on convergence rates}
We give two theorems on the learning efficiency of MCCR in this subsection.

\begin{theorem}
Suppose that Assumptions 1, 2 and $f^{*}\in \mathcal{H}$ hold. For a fixed scale parameter $\sigma>0$, we have
$$f^*=\mathop{\arg\min}_{f\in \mathcal{H}}{\varepsilon}^{\sigma}(f),$$
and for any $f\in \mathcal{H}$,
$$c_{\sigma}\Vert f-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2\leq {\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)\leq \Vert f-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2,$$
where $c_{\sigma}$ is a positive constant which will be given explicitly in the  proof.
\end{theorem}

Theorem 1 shows that the population version estimator of MCCR can accurately represent the unknown function $f^*$ when $f^{*}\in \mathcal{H}$ and the noise is a mixture of symmetric distributions with positive Fourier transform and the Fourier transform has a lower bound of positive constant in a finite region. In this sense, $f_z$ can be regarded as an unbiased estimator of $f^*$. In addition, if let $f=f_z$, Theorem 1 shows that the excess risk of MCCR can be bounded by $L_{\rho_\mathcal{X}}^2$ distance between the MCCR estimator $f_z$ and the real unknown function $f^*$, which helps us to get Theorem 2.

\begin{theorem}
Suppose that Assumptions 1, 2, 3 and $f^{*}\in \mathcal{H}$ hold. Let $f_z$ be given in equation (3) and scale parameter $\sigma=n^{\theta}$. For any $0<{\delta}<1$, with probability at least $1-\delta$, it holds that
$$\Vert f_z-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2 \lesssim n^{-\frac{8\theta+2}{2+s}}\log\left(\frac{1}{\delta}\right),$$
where $\theta \in (-\frac{1}{4},0)$.
\end{theorem}

Theorem 2 shows that, when $n$ goes to infinity, the estimator $f_z$ can learn the conditional mean function $f^*$ by the rate ${\mathcal{O}}(n^{-\frac{8\theta+2}{2+s}})$ as the scale parameter $\sigma$ tends to 0 via ${\mathcal{O}}(n^{\theta})$  with $\theta \in (-\frac{1}{4},0)$, and the optimal rate is \begin{math} {\mathcal{O}}(n^{-1})\end{math} in the sense of approximation. Together with \cite{2015Learning}, where the scale parameter goes to infinity or is a constant large enough are discussed, this work completes the theoretical analysis for MCCR. And, the convergence rate is established when the noise is the mixture of symmetric distributions with some common conditions, which ensures the model to learn the unknown mean regression functions under mild conditions. Last, this result reveals that MCCR models hold favorable learning rate as well as adequate robustness via tending to 0 scale parameters.  All these lead to the merits of MCCR models in dealing with kinds of noises and outliers, some cases will be shown in Section 4.

\subsection{ Proofs of Theorem 1 and 2}
First, we display the following lemma which is needed to prove Theorem 2. It is given in \cite{2007Multi} and adopted in kinds of statistical learning works such as \cite{guo2013concentration, 2018Learning} and \cite{feng2020statistical}.

\begin{lemma}
Let $\mathcal{F}$ be a set of measurable functions on $\mathcal{Z}$, and $B,c>0$, $\tau \in [0,1]$ be constants such that each of $f\in \mathcal{F}$ satisfies $\Vert f\Vert_{\infty}\leq B$ and $\mathbb{E}f^2\leq c(\mathbb{E}f)^\tau$. If for some $a>0$ and $s\in (0,2)$,
$$\log \mathcal{N}_{2}(\mathcal{F},\eta)\leq a{\eta}^{-s},\forall \eta>0,$$
then there exists a constant $c_s$ such that for any $t>0$, with probability at least $1-e^{-t}$, for all $f\in \mathcal{F}$, there holds
$$\mathbb{E}f-\frac{1}{n}\sum_{i=1}^{n}f(z_i)\leq \frac{1}{2}{\zeta}^{1-\tau}(\mathbb{E}f)^{\tau}+c_{s}\zeta+2\left(\frac{ct}{n}\right)^{\frac{1}{2-\tau}}+\frac{18Bt}{n},$$
where
$$\zeta:=max\left\{c^{\frac{2-s}{4-2\tau+s\tau}}\left(\frac{a}{n}\right)^{\frac{2}{4-2\tau+s\tau}},
B^{\frac{2-s}{2+s}}\left(\frac{a}{n}\right)^{\frac{2}{2+s}}\right\}.$$
\end{lemma}
\
\\
\textbf{Proof of Theorem 1.} The techniques we used here is similar with \cite{fan2016consistency}. First, we prove the left inequality. From the definition of ${\varepsilon}^{\sigma}(f)$, we have
\[
\begin{split}
{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)&=
\int\limits_{\mathcal{X} \times \mathcal{Y}}l_{\sigma}(y-f(x))d\rho-\int\limits_{\mathcal{X} \times \mathcal{Y}}l_{\sigma}(y-f^*(x))d\rho  \\
&={\sigma}^2\int\limits_\mathcal{X}\int\limits_\mathcal{Y}
\left[-\exp\left(-\frac{(y-f(x))^2}{{\sigma}^2}\right)+\exp\left(-\frac{(y-f^*(x))^2}{{\sigma}^2}\right)\right]p_{y|x}(y)dyd{\rho}_\mathcal{X}(x)\\
&={\sigma}^2\int\limits_\mathcal{X}\int_{-\infty}^{+\infty}
\left[\exp\left(-\frac{t^2}{{\sigma}^2}\right)-\exp\left(-\frac{(t-(f(x)-f^*(x)))^2}{{\sigma}^2}\right)\right]p_{\varepsilon|x}(t)dtd{\rho}_\mathcal{X}(x)\\
\end{split}
\]
By Plancherel formula, we have
\[
\begin{split}
{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)&=
\frac{{\sigma}^3}{2\sqrt{\pi}}\int\limits_\mathcal{X}\int_{-\infty}^{+\infty}
\exp\left(-\frac{{\sigma}^2{\xi}^2}{4}\right)[1-\exp(-i{\xi}(f(x)-f^*(x)))]\widehat{p_{\varepsilon|x}}(\xi)d{\xi}d{\rho}_\mathcal{X}(x)\\
&=\frac{{\sigma}^3}{2\sqrt{\pi}}\int\limits_\mathcal{X}\int_{-\infty}^{+\infty}\exp\left(-\frac{{\sigma}^2{\xi}^2}{4}\right)
2{\sin}^2\left(\frac{\xi(f(x)-f^*(x))}{2}\right)\widehat{p_{\varepsilon|x}}(\xi)d{\xi}d{\rho}_\mathcal{X}(x)
\end{split}
\]
where $\widehat{p_{\varepsilon|x}}$ is the Fourier transform of $p_{\varepsilon|x}$, the second equation holds because ${\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)$ is real for $\forall f\in \mathcal{H}$. From the linearity of the Fourier transformation and Assumption 1, we get
$$\widehat{p_{\varepsilon|x}}(\xi)=\sum_{i=1}^{K}w_{i}\widehat{p_{\varepsilon,i|x}}(\xi).$$
where $\widehat{p_{\varepsilon,i|x}}$ is the Fourier transformation of $p_{\varepsilon,i|x}$ ($i=1,2,\cdots,K$). Furthermore, from Assumption 2, it is known that exists positive constants $c_{0,i}, C_{0,i}$ with $i=1,2,\cdots,K$, such that
\[
\begin{split}
{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)
&=\frac{{\sigma}^3}{\sqrt{\pi}}\int\limits_\mathcal{X}\sum_{i=1}^{K}w_i\int_{-\infty}^{+\infty}\exp\left(-\frac{{\sigma}^2{\xi}^2}{4}\right)
{\sin}^2\left(\frac{\xi(f(x)-f^*(x))}{2}\right)\widehat{p_{\varepsilon,i|x}}(\xi)d{\xi}d{\rho}_\mathcal{X}(x)\\
&\geq\frac{{\sigma}^3}{\sqrt{\pi}}\int\limits_\mathcal{X}\sum_{i=1}^{K}w_i\int_{-c_{0,i}}^{c_{0,i}}\exp\left(-\frac{{\sigma}^2{\xi}^2}{4}\right)
{\sin}^2\left(\frac{\xi(f(x)-f^*(x))}{2}\right)C_{0,i}d{\xi}d{\rho}_\mathcal{X}(x)
\end{split}
\]
And, for $\forall x\in \mathcal{X}$, $|f(x)-f^*(x)| \leq 2M$. When $|\xi| \leq \frac{\pi}{2M}$, from Jordan's inequality, it holds that
$$\frac{{\xi}^{2}(f(x)-f^*(x))^2}{{\pi}^2} \leq {\sin}^2\left(\frac{\xi(f(x)-f^*(x))}{2}\right).$$
Then, let $c^{'}=\min\{c_{0,1},c_{0,2},\cdots,c_{0,K},\frac{\pi}{2M}\}$, so we have
\begin{equation}
\begin{split}
{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)& \geq
\frac{{\sigma}^3}{{\pi}^{\frac{5}{2}}}\int\limits_\mathcal{X}\sum_{i=1}^{K}w_i \int_{-c^{'}}^{c^{'}}
{\xi}^{2}\exp\left(-\frac{{\sigma}^2{\xi}^2}{4}\right)C_{0,i}(f(x)-f^*(x))^2d{\xi}d{\rho}_\mathcal{X}(x)\\
&:=c_{\sigma}\int\limits_\mathcal{X}(f(x)-f^*(x))^2d{\rho}_\mathcal{X}(x)
\end{split}
\end{equation}
where
\begin{equation}
\begin{split}
c_{\sigma}=\frac{{\sigma}^3}{{\pi}^{\frac{5}{2}}}\sum_{i=1}^{K}w_i C_{0,i}\int_{-c^{'}}^{c^{'}}{\xi}^2
\exp\left(-\frac{{\sigma}^2{\xi}^2}{4}\right)d{\xi}.
\end{split}
\end{equation}
is a positive constant. This implies that ${\varepsilon}^{\sigma}(f)\geq {\varepsilon}^{\sigma}(f^*)$ for any $f\in \mathcal{H}$. In other words,
$$f^*=\mathop{\arg\min}\limits_{f\in \mathcal{H}}{\varepsilon}^{\sigma}(f).$$
Second,
\[
\begin{aligned}
&{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)={\sigma}^2\int\limits_\mathcal{X}\int_{-\infty}^{+\infty} \left[\exp\left(-\frac{t^2}{{\sigma}^2}\right)-\exp\left(-\frac{(t-(f(x)-f^*(x)))^2}{{\sigma}^2}\right)\right]p_{\varepsilon|x}(t)dtd{\rho}_\mathcal{X}(x)\\
&:={\sigma}^2\int\limits_\mathcal{X}[F_{x}(f(x)-f^*(x))-F_{x}(0)]d{\rho}_\mathcal{X}(x)
\end{aligned}
\]
where $F_{x}:\mathbb{R} \to \mathbb{R}$ is defined as
$$F_{x}(u):=1-\int_{-\infty}^{+\infty}\exp\left(-\frac{(t-u)^2}{{\sigma}^2}\right)p_{\varepsilon|x}(t)dt,x\in \mathcal{X}.$$
From Taylor expansions, we know that
$$F_{x}(f(x)-f^*(x))-F_{x}(0)=F_{x}^{'}(0)(f(x)-f^*(x))+\frac{F_{x}^{''}({\zeta}_x)}{2}(f(x)-f^*(x))^2,$$
where for $\forall x\in \mathcal{X}$, ${\zeta}_{x}$ is between 0 and $f(x)-f^*(x)$. Due to the noise is symmetric with 0 as the center, for $\forall x\in \mathcal{X}$, we have
$$F_{x}^{'}(0)=-\int_{-\infty}^{+\infty}\exp\left(-\frac{t^2}{{\sigma}^2}\right)\left(\frac{2t}{{\sigma}^2}\right)p_{\varepsilon|x}(t)dt=0.$$
In addition, for $\forall x\in \mathcal{X}$,
\[
\begin{split}
F_{x}^{''}({\zeta}_x)&=\int_{-\infty}^{+\infty}2\exp\left(-\frac{(t-{\zeta}_x)^2}{{\sigma}^2}\right)\left(\frac{{\sigma}^2-2(t-{\zeta}_x)^2}{{\sigma}^4}\right)p_{\varepsilon|x}(t)dt \leq \frac{2}{{\sigma}^2}.
\end{split}
\]
Therefore,
\begin{equation}
\begin{split}
{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)&={\sigma}^2\int\limits_\mathcal{X}0+\frac{F_{x}^{''}({\zeta}_x)}{2}(f(x)-f^*(x))^2d{\rho}_\mathcal{X}(x)\\
&\leq \int\limits_\mathcal{X}(f(x)-f^*(x))^2d{\rho}_\mathcal{X}(x)=\Vert f-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2.
\end{split}
\end{equation}
Combine (4) and (6), we obtain
$$c_{\sigma}\Vert f-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2\leq {\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)\leq \Vert f-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2,$$
where $c_{\sigma}$ is a positive constant given in (5). This completes the proof of Theorem 1.

\begin{remark}
When $\sigma$ changes with the sample size $n$, we will make further analysis here. From (5), it is known that there exists a constant $A_i\in [-c^{'},c^{'}], i=1,2,\cdots,K$, such that,
\[
\begin{split}
c_{\sigma}&=\frac{{\sigma}^3}{{\pi}^{\frac{5}{2}}}\sum_{i=1}^{K}w_i{A_i}^2
\exp\left(-\frac{{\sigma}^2{A_i}^2}{4}\right)C_{0,i}(c^{'}+c^{'})\\
&={\sigma}^3\sum_{i=1}^{K}\exp\left(-\frac{{\sigma}^2{A_i}^2}{4}\right)2{\pi}^{-\frac{5}{2}}c^{'}w_i{A_i}^2C_{0,i}\\
&:={\sigma}^3\sum_{i=1}^{K}\exp\left(-\frac{{\sigma}^2{A_i}^2}{4}\right)c_1,
\end{split}
\]
where
$$c_1=2{\pi}^{-\frac{5}{2}}c^{'}w_i{A_i}^2C_{0,i},$$
is a positive constant.
\end{remark}
\
\\ \textbf{Proof of Theorem 2.}
First, we prove that Theorem 2 satisfies the conditions of Lemma 1. We assume that the definition of the function set $\mathcal{F}_{\mathcal{H}}$ is as follows
$$\mathcal{F}_{\mathcal{H}}=\left\{g|g(z)=-{\sigma}^2\exp\left\{-\frac{(y-f(x))^2}{{\sigma}^2}\right\}+{\sigma}^2\exp\left\{-\frac{(y-f^*(x))^2}{{\sigma}^2}\right\}
,f\in \mathcal{H},z\in \mathcal{Z}\right\}.$$
Then, for $\forall g\in \mathcal{F}_{\mathcal{H}}$, we have
$$\Vert g\Vert_{\infty}\leq {\sigma}^2+{\sigma}^2\leq 2{\sigma}^2,$$
and
$$\mathbb{E}g={\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*).$$
Moreover, introducing the auxiliary function $h(t)=-{\sigma}^2\exp\{-\frac{t^2}{\sigma^2}\}$, $t \in \mathbb{R}$, it is easy to see that $\| h^{'} \|_{\infty}=\sqrt{2/e} \sigma$. By taking $t_1=y-f(x)$, $t_2=y-f^{*}(x)$ and applying the mean value theorem to $h$, we see that
\[
\begin{split}
\mathbb{E}g^2&=\int\limits_\mathcal{Z}\left(-{\sigma}^2\exp\left\{-\frac{(y-f(x))^2}{{\sigma}^2}\right\}+{\sigma}^2\exp\left\{-\frac{(y-f^*(x))^2}{{\sigma}^2}\right\}\right)^2d{\rho}\\
&\leq \int\limits_\mathcal{Z}\|h^{'} \|_{\infty}^2(f(x)-f^*(x))^2d{\rho}\\
&\leq 2e^{-1}{\sigma}^2\int\limits_\mathcal{Z}(f(x)-f^*(x))^2d{\rho}\\
&\lesssim {\sigma}^{-1}\left(\sum_{i=1}^{K}\exp\left\{-\frac{\sigma^2A_i^2}{4}\right\}\right)^{-1}\mathbb{E}g.
\end{split}
\]
On the other hand, for $\forall g_1,g_2\in \mathcal{F}_{\mathcal{H}},$ $\exists f_1,f_2\in \mathcal{H}$, such that
$$g_1(z)=-{\sigma}^2\exp\left\{-\frac{(y-f_1(x))^2}{{\sigma}^2}\right\}+{\sigma}^2\exp\left\{-\frac{(y-f^*(x))^2}{{\sigma}^2}\right\},$$
and
$$g_2(z)=-{\sigma}^2\exp\left\{-\frac{(y-f_2(x))^2}{{\sigma}^2}\right\}+{\sigma}^2\exp\left\{-\frac{(y-f^*(x))^2}{{\sigma}^2}\right\}.$$
So, we have
$$\| g_1-g_2 \|_{\infty} \leq \sqrt{2/e} \sigma \cdot \| f_1-f_2 \|_{\infty}.$$
Then, under the assumption of complexity and $0<s<2$, the $l^{2}$-empirical covering numbers of $\mathcal{F}_{\mathcal{H}}$ and $\mathcal{H}$ have the following relation
$$\log \mathcal{N}_2(\mathcal{F}_{\mathcal{H}},\eta)\leq \log \mathcal{N}_2(\mathcal{H},\frac{\eta}{\sqrt{2/e}\sigma})\lesssim \sigma^s \eta^{-s}.$$
Now, applying Lemma 1 to the random variable $g$ with $B=2\sigma^2$, $c={\sigma}^{-1}\left(\sum_{i=1}^{K}\exp\left\{-\frac{\sigma^2A_i^2}{4}\right\}\right)^{-1}$, $\tau=1$, $a=\sigma^s$, then for any $0<{\delta}<1$, with probability at least $1-\delta$, it holds that
$$[{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)]-[{\varepsilon}_{z}^{\sigma}(f)-{\varepsilon}_{z}^{\sigma}(f^*)]-\frac{1}{2}[{\varepsilon}^{\sigma}(f)-{\varepsilon}^{\sigma}(f^*)]
\lesssim (\sigma^{\frac{3s-2}{2+s}}n^{-\frac{2}{2+s}}+\sigma^{-1}n^{-1})\log\left(\frac{1}{\delta}\right).$$
Because $f_z=\mathop{\arg\min}\limits_{f\in \mathcal{H}}{\varepsilon}^{\sigma}_{z}(f)$, we have
\[
\begin{split}
{\varepsilon}^{\sigma}(f_z)-{\varepsilon}^{\sigma}(f^*) & \leq 2[{\varepsilon}^{\sigma}(f_z)-{\varepsilon}^{\sigma}(f^*)]-2[{\varepsilon}_{z}^{\sigma}(f_z)-{\varepsilon}_{z}^{\sigma}(f^*)]-[{\varepsilon}^{\sigma}(f_z)-{\varepsilon}^{\sigma}(f^*)]\\
& \lesssim (\sigma^{\frac{3s-2}{2+s}}n^{-\frac{2}{2+s}}+\sigma^{-1}n^{-1})\log\left(\frac{1}{\delta}\right).
\end{split}
\]
Therefore, according to Theorem 1, for any $0<{\delta}<1$, with probability at least $1-\delta$, it holds that
\[
\begin{split}
\Vert f_{z}-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2 & \leq \frac{1}{c_\sigma}[{\varepsilon}^{\sigma}(f_z)-{\varepsilon}^{\sigma}(f^*)]\\
& \lesssim (\sigma^{\frac{3s-2}{2+s}}n^{-\frac{2}{2+s}}+\sigma^{-1}n^{-1})\sigma^{-3}\left(\sum_{i=1}^{K}\exp\left\{-\frac{\sigma^2A_i^2}{4}\right\}\right)^{-1}\log\left(\frac{1}{\delta}\right)\\
& \lesssim (\sigma^{-\frac{8}{2+s}}n^{-\frac{2}{2+s}}+\sigma^{-4}n^{-1})\log\left(\frac{1}{\delta}\right).
\end{split}
\]
Let $\sigma=n^{\theta}$ with $\theta \in (-\frac{1}{4},0)$, then we have
$$\Vert f_{z}-f^*\Vert_{L_{\rho_\mathcal{X}}^2}^2 \lesssim n^{-\frac{8\theta+2}{2+s}}\log\left(\frac{1}{\delta}\right).$$
This completes the proof of Theorem 2.

\vspace*{+0.1cm}
\section{Simulations }
In this section, we use synthetic data of finite samples to show the effectiveness of MCCR model under mixed symmetric noises and compare the robustness of MCCR to Huber and least squares (LS) regression models. Assume further that the space $\mathcal{H}$ is a bounded subset of a reproducing kernel Hilbert space $\mathcal{H}_{\mathcal{K}}$, where $\mathcal{K}$ is a Mercer kernel. At this time, the solution (3) of the MCCR model is equivalent to
\begin{equation}
f_z=\mathop{\arg\min}\limits_{f\in \mathcal{H}_{\mathcal{K}}}\frac{1}{n}{\sum_{i=1}^{n}l_\sigma(y_i-f(x_i))}+\lambda\Vert f \Vert_{\mathcal{K}}^2,
\end{equation}
where $\lambda$ is a positive regularization parameter. The representor theorem ensures that $f_z$ can be modeled by
$$f_{z}(x)=\sum_{i=1}^{n}\alpha_{i}\mathcal{K}(x,x_i)+b_0,$$
where $\bm{\alpha}=(\alpha_1,\cdots,\alpha_n)^\top \in \mathbb{R}^n$ and $b_0 \in \mathbb{R}$ are learned from (7). We use Gaussian kernel $\mathcal{K}(x,x')=\exp\{-\frac{\Vert x-x'\Vert^2}{h^2}\}$ and the iterative weighted least squares to perform the algorithms. The initial values of the iterations  for MCCR regression and Huber regression are $\bm{\alpha}=(0,\cdots,0)^\top$ and $b_0=0$. The scale parameter $\sigma$,  bandwidth parameter $h$ and regularization parameter $\lambda$ in $l_{\sigma}$ loss and Huber's loss, the bandwidth parameter $h$ and  regularization parameter $\lambda$ in LS regression are all selected by the rule of five-fold cross validation.

\vspace*{+0.2cm}
\noindent \textbf{Example 1.} The real function is $f(x)=e^{-7.5x}\cos(10\pi x)$, where $X_{i}$ are independent sample points from a uniform distribution $U(0,1)$, and $Y_{i}=f(X_{i})+\varepsilon_{i}$ with $\varepsilon_{i}\sim 0.8N(0,0.1^2)+0.2N(0,0.5^2)$. $N(0,0.1^2)$ is the background noise and $N(0,0.5^2)$ is used to generate outliers, $i=1, 2, \cdots, n$ and $n=200$. The estimates  are shown in Fig.1.\\

\begin{figure*}[htbp]
\centering
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{1-1-eps-converted-to}}
\hspace*{1cm}
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{1-2-eps-converted-to}}
\caption{Real function, estimates and scatter plots of Example 1. (a) Real function (black solid curve) and scatter plots. Circles are sample points with background noise and red crosses are outliers. (b) Red, blue and green solid curves are the estimates of MCCR,  Huber and LS regressions, respectively.}
\label{fig:1}
\end{figure*}

\vspace*{+0.2cm}
\noindent \textbf{Example 2.} The real function is $f(x)= -1+1.5x+0.2\phi(x-0.6)$, where $\phi(x)$ is the density function of normal distribution $N(0,0.04^2)$, $X_{i}$ are independent sample points from a uniform distribution $U(0,1)$, and $Y_{i}=f(X_{i})+\varepsilon_{i}$ with $\varepsilon_{i}\sim 0.8N(0,0.1^2)+0.2N(0,1)$. $N(0,0.1^2)$ is the background noises and $N(0,1)$ is used to generate outliers, $i=1, 2, \cdots, n$ and $n=200$. The estimates  are shown in Fig.2.

\begin{figure*}[htbp]
\centering
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{2-1-eps-converted-to}}
\hspace*{1cm}
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{2-2-eps-converted-to}}
\caption{Real function, estimates and scatter plots of Example 2. (a) Real function (black solid curve) and scatter plots. Circles are sample points with background noise and red crosses are outliers. (b) Red, blue and green solid curves are the estimates of  MCCR,  Huber  and LS  regressions, respectively.}
\label{fig:2}
\end{figure*}

\vspace*{+0.2cm}
\noindent \textbf{Example 3.} The real function is $f(x)=\frac{\sin(x)}{x}$, where $X_{i}$ are independent sample points from a uniform distribution $U(-10,10)$, and $Y_{i}=f(X_{i})+\varepsilon_{i}$ with $\varepsilon_{i}\sim 0.8N(0,0.1^2)+0.2Cauchy(0,0.2)$. $N(0,0.1^2)$ is the background noise and $Cauchy(0,0.2)$ is used to generate outliers, $i=1, 2, \cdots, n$ and $n=200$. The estimates are shown in Fig.3.

\begin{figure*}[htbp]
\centering
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{3-1-eps-converted-to.pdf}}
\hspace*{1cm}
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{3-2-eps-converted-to}}
\caption{Real function, estimates and scatter plots of Example 3. (a) Real function (black solid curve) and scatter plots. Circles are sample points with background noise and red crosses are outliers. (b) Red, blue and green solid curves are the estimates of MCCR,  Huber and LS regressions, respectively.}
\label{fig:3}
\end{figure*}

Estimates shown in these three Figures tell us that the MCCR has a comparable performance with Huber regression on fitness and robustness, and both of them outperform LS regression when data contains outliers. Moreover, for each example, the simulations are repeated 100 times and the mean square error(MSE) between the estimator $f_z$ and the real function $f^*$ are given in Table 1 ( in brackets are the standard deviations), which supports the similar conclusion. Furthermore, let $\theta=-0.2$ and the sample size  $n$ goes from $100$ to $200, 300, \cdots, 600$, MSEs and their standard deviations are shown in Table 2 and Fig.4, which provide some evidence to Theorem 2.

\begin{table}[htbp]
\centering
\caption{The MSE between $f_z$ and $f^*$ of MCCR, Huber, and LS regression}
\label{tab:1}
\begin{tabular}{llll}
\hline\noalign{\smallskip}
  & MCCR & Huber & LS  \\
\noalign{\smallskip}\hline\noalign{\smallskip}
Example 1 & $\mathbf{0.0014}(0.0006)$ & 0.0017(0.0006) & 0.0046(0.0015)\\
Example 2 & $\mathbf{0.0027}(0.0008)$ & 0.0045(0.0018) & 0.0165(0.0062)\\
Example 3 & $\mathbf{0.0006}(0.0004)$ & 0.0007(0.0004) & 0.2595(1.5240)\\
\noalign{\smallskip}\hline
\end{tabular}
\end{table}

\begin{table}[htbp]
\centering
\caption{The MSE of MCCR for sample size $n$}
\label{tab:2}
\begin{tabular}{llll}
\hline\noalign{\smallskip}
  &Example 1 & Example 2 & Example 3 \\
\noalign{\smallskip}\hline\noalign{\smallskip}
 n=100 & 0.0034(0.0025) & 0.0066(0.0275) & 0.0013(0.0006) \\
 n=200 & 0.0017(0.0025) & 0.0035(0.0081) & 0.0006(0.0003) \\
 n=300 & 0.0010(0.0003) & 0.0031(0.0095) & 0.0004(0.0002) \\
 n=400 & 0.0007(0.0002) & 0.0026(0.0078) & 0.0003(0.0001) \\
 n=500 & 0.0006(0.0002) & 0.0013(0.0002) & 0.0003(0.0001) \\
 n=600 & 0.0005(0.0001) & 0.0011(0.0002) & 0.0002(0.0001) \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table}

\begin{figure}[htbp]
\centering
  \includegraphics[width=10cm,height=7cm]{mse-eps-converted-to}
\caption{MSE of MCCR  for  sample size $n$ .}
\label{fig:4}
\end{figure}

\vspace*{+1cm}
\section{Application}
In this section, we use MCCR, Huber and LS regression to learn from the speed-flow data, which was firstly studied in \cite{petty1996freeway} and is publicly available from the \texttt{R}-package \texttt{hdrcde}. Let $X$ be the traffic flow, evaluated by the number of vehicles per hour per lane,  $Y$ be the speed of the vehicle, measured in miles per hour, scatter plots of two data sets collected on two separate lanes (lane 2 and lane 3) of the 4-lane Californian highway I-880 in 1993 are displayed in Fig.5. Red, blue and green solid curves are the estimates of MCCR, Huber and LS regressions, respectively. The algorithm is the iterative weighted least squares and the hyperparameter selection method is five-fold cross-validation. The results show that MCCR holds the best robustness, LS regression is the worst and Huber regression is between the two.

\begin{figure*}
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{4-2-eps-converted-to}}
\hspace*{1cm}
\subfigure[]{
\includegraphics[width=2.5in,height=2.5in]{5-2-eps-converted-to}}
\caption{Estimates and scatter plots of  the speed-flow data. Red, blue and green solid curves are the estimates of  MCCR,  Huber  and LS  regressions, respectively. (a) Speed-Flow data for lane $2$ ; (b) Speed-Flow data for lane 3.}
\label{fig:5}
\end{figure*}

\section{Conclusion}
Maximum Correntropy Criterion Regression models with tending-to-zero scale parameters are theoretically studied within the frame of statistical learning. The model can learn the conditional mean function $f^*$ with the rate ${\mathcal{O}}(n^{-\frac{8\theta+2}{2+s}})$ as the scale parameter $\sigma$ tends to 0 via ${\mathcal{O}}(n^{\theta})$, where $\theta \in (-\frac{1}{4},0)$ and $s$ is the complexity parameter of the hypothesis space. The optimal convergence rate of estimator is \begin{math} {\mathcal{O}}(n^{-1})\end{math} in the sense of approximation. This study implies that MCCR models can attain both accurate approximation and adequate robustness when the sample size is large enough.  Simulations and applications of MCCR, Huber and LS regression models are evaluated, both of them show some evidences for that the MCCR model holds the best robustness.\par
\ \\
\ \\
\textbf{Acknowledgements}
\ \\
This research was supported by the National Science Foundation of China (Grant NO. 11671012), the Key University Science Research Project of Anhui, China (Grant NO. KJ2017A028) and  the Open Project Program of School of Mathematical Sciences of Anhui University (Grant NO. Y01002431).
\ \\
\ \\

\bibliographystyle{apalike}