\section{Introduction}

In the course of the increasing popularity of deep neural networks for applications in technical, ecological and many other fields, there has been huge progress in the research on understanding the mathematical properties of the mathematical mapping implemented by a deep neural network. The approximation properties, or expressiveness, of neural network functions have attracted intense interest. The central result in this field is the classic \emph{universal approximation theorem}, which states that any continuous function can be approximated with arbitrary accuracy (in terms of uniform approximation or $L_p$ norms) by neural network functions that have only one hidden layer for nearly every activation function \citep{cybenko1989approximation, hornik1991approximation}. In recent years, these kind of results have become more precise, for instance, in terms of estimates on the order of magnitude of parameters needed to achieve an approximation of a certain accuracy \citep{petersen2018optimal,yarotsky2017error,yarotsky2018optimal}, we refer to \citep{devore2020neural} for an overview on recent developments. Since empirical observations showed that depth has a significant impact on the performance of neural networks, a lot of research has been dedicated to the effect of depth on the expressive power of neural networks \citep{telgarsky2016benefits, mhaskar2016deep, montufar2014number,raghu2016expressive, rolnick2017power,lin2017does,cohen2016expressive}.
Besides this, the role of width in expressiveness of network functions has been investigated \citep{lu2017expressive, hanin2017approximating,hanin2017universal,johnson2018deep,kidger2019universal,park2020minimum}.
An important result states that a width larger than the input dimension is needed to allow universal approximation on arbitrary compact sets. It has also been shown that the capability to learn disconnected or bounded decision regions depends on whether the width is larger than the input dimension \citep{nguyen2018neural,beise2018decision}.

In a further line of research, the capabilities of neural networks to memorize finite sample sets is investigated  \citep{yun2018small,vershynin2020memory,bubeck2020network,park2020provable}, which provides mathematical foundations to the surprising findings that neural networks can interpolate the training data \citep{zhang2016understanding} in practise.

In this work, we investigate what kind of subsets $M\subset \mathbb{R}^n$ admit universal approximation by neural network functions with a maximum layer width less than or equal to the input dimension, i.e. network functions that violate the necessary condition for universal approximation given in \citet{hanin2017universal} and \citet{park2020minimum} for arbitrary compact sets.

Let us introduce the following notation: For a set $D\subset \mathbb{R}^n$ we denote by $D^\circ$ the set of interior points, by $\overline{D}$ its closure and by $\partial D=\overline{D}\setminus D^\circ$ the boundary of $D$. $\left\Vert \cdot\right\Vert$ is the Euclidean norm.
For a linear subspace $U\subset \mathbb{R}^n$, we denote by $\dim(U)$ its dimension and by $P_U:\mathbb{R}^n\rightarrow U$ the orthogonal linear projection on $U$. We denote by $e_j$ the unit vector of the $j$-te coordinate axis in $\mathbb{R}^n$. For some $v\in\mathbb{R}^n$ and $j\in\{1,...,n\}$, we write $v^{(j)}$ for the $j-$th component of $v$.

We consider neural network functions $F:\mathbb{R}^{n_{0}}\rightarrow \mathbb{R}^{n_{L}}$ where $n_0$ is called the input dimension and $n_L$ the output dimension. Our network functions have the following form
\begin{equation}\label{networkFun}
	F:=W_L(A_{L-1}\circ ... \circ A_1)+b_L
\end{equation}
where $A_j(x)=\sigma(W_jx+b_j)$ with $W_j\in \mathbb{R}^{n_{j}\times n_{j-1}}$ (weights), $b_j\in \mathbb{R}^{n_j}$ (bias), where $j=1,...,L$, and $\sigma:\mathbb{R} \rightarrow \mathbb{R}$ the activation function. The application of $\sigma$ and the inverse image $\sigma^{-1}$ are understood to be applied element wise to vectors or subsets of $\mathbb{R}^n$. The widely used activation function \emph{rectified linear unit}, shortly ReLU, is defined by $t\mapsto\max\{t,0\}$. As a shorthand notation of the function implemented by the first $k$ layers of a network function, we set
\begin{equation}
	F_k=W_{k}\left(A_{k-1}\circ ... \circ A_1\right)+b_k,
\end{equation}
for $k\in\{1,...,L-1\}$.
We call $n_j$  the \emph{width} of layer $j=1,...,L$. The width of the network is defined as $\omega(F)=\max\{n_j: \, j=1,...,L\}$ and $L$ is called the depth of the network. For $m,L \in \mathbb{N}$ and an activation function $\sigma$, $\mathcal{NN}_{\sigma}^{k}(m,L)$ is the set of network functions $F:\mathbb{R}^{n_0}\rightarrow \mathbb{R}^{k}$ of the form (\ref{networkFun}) with activation $\sigma$ and of maximum width $m$, i.e. $\omega(F)\leq m$, and depth $L$, i.e. $L$ layers including the final linear layer. In case that $ \mathbb{R}^{k}=\mathbb{R}$ we omit the latter and write $\mathcal{NN}_{\sigma}(m,L)$ instead, and in case the the depth is not specified, meaning that arbitrary depth is allowed, we write $\mathcal{NN}_{\sigma}^{k}(m)$, $\mathcal{NN}_{\sigma}(m)$, respectively.

We say that a compact subset $M\subset\mathbb{R}^{n_0}$ admits universal uniform approximation of functions in some function class of mappings from $M$ to $\mathbb{R}$ by network function of $\mathcal{NN}_{\sigma}(m,L)$ or $\mathcal{NN}_{\sigma}(m)$, if for every $f:M\rightarrow \mathbb{R}$ in that class and every $\varepsilon>0$, there exists some  $F\in\mathcal{NN}_{\sigma}(m,L)$, respectively $F\in\mathcal{NN}_{\sigma}(m)$, such that
\begin{equation}\label{conditionUniformApprox}
	\max_{x\in M}\lvert f(x)-F(x)\rvert<\varepsilon.
\end{equation}

\section{Related work}

Approximation properties of width bounded neural networks have been studied in several papers. A common goal in these lines of research is to provide upper and lower bounds on the minimum width of the network, i.e. the minimum number of neurons per layers, here denoted by $\omega_\text{min}$, needed for universal approximation in $C(M,\mathbb{R}^{n_L})$, $L_p(\mathbb{R}^{n_0},\mathbb{R}^{n_L})$ and $L_p(M,\mathbb{R}^{n_L})$ for  $M\subset\mathbb{R}^{n_0}$. Since in this work we are interested in universal approximation in $C(M,\mathbb{R})$, we refer to \citet{lu2017expressive}, \citet{kidger2019universal} and \citet{park2020minimum} for the case of $L_p$ approximation. For the case of $C(M,\mathbb{R}^{n_L})$, it was proven in \citet{hanin2017approximating} that
\begin{equation}
	\small{n_0+1\leq \omega_\text{min}\leq n_0+ n_L, \text{  for ReLU activation},}
\end{equation}
which has been tightened in \citet{park2020minimum} to
\begin{equation}
	\small{\omega_\text{min}\geq \max\{n_0+1,n_L\} \text{  for ReLU and STEP activation},}
\end{equation}
where STEP refers to the threshold activation that maps to $1$ on $\{x\geq 0\}$ and $0$ otherwise.
It has been shown at the same time in \citet{johnson2018deep} and \citet{beise2018decision} that
\begin{equation}\label{ourEstimate}
	\small{\omega_\text{min}\geq n_0+1 \text{  for ReLU + injective continuous activation}.}
\end{equation}
The activation functions allowed in \citet{johnson2018deep} are slightly more general, as it is only required that they admit arbitrary accurate uniform approximation by injective continuous functions on arbitrary compact subset of $\mathbb{R}$. In \citet{kidger2019universal} the above upper bounds are extended to further classes of activation functions. For these results and a general overview on recent developments in this field, we refer to \citet{park2020minimum}.

In this work, we only consider width less than or equal to $n_0$ and investigate uniform approximation on certain compacts subsets of the input space. According to the above result (\ref{ourEstimate}), we cannot expect universal approximation for all compact $M\subset \mathbb{R}^{n_0}$. However, a common assumption in applications is that approximation is only needed on a certain subset.
We derive some topological conditions on $M$ that allow or exclude a kind of universal approximation under these circumstances.
Although narrow neural networks, as they are considered in this work, are not common in practical applications, many networks exhibit a decaying layer width at the later layers. From this perspective, our results give theoretical insights on the kind feature extraction earlier layers need to implement in order to allow the later layers to solve a given machine learning task.
Our main contributions are as follows:
\begin{enumerate}

	\item We show that for the case of two disjoint compact sets, the existence of a cone-like sector that contains one of these sets and does not intersect with the other one, is sufficient to allow exact fit of functions that take constant values on each of these sets by network functions from $\mathcal{NN}_\mathrm{ReLU}(n_0,4)$. A weaker result is concluded for the case of multiply pairwise compact components.
	\item We show that any function on an arbitrary finite sample set can be approximated with arbitrary accuracy by network functions in $\mathcal{NN}_{\cos} (1,3)$.
	\item A maximum principle is given for $F\in\mathcal{NN}_{\sigma}(n_0)$ with $\sigma$ continuous and monotonic. This allows to conclude that the lower bound $\omega_\text{min}\geq n_0+1$ in (\ref{ourEstimate}) is sharp for a wide range of subsets $M$, e.g. when $M=\partial D\cup\{c\}$ for some $D$ with non-empty interior and $c\in D^{\circ}$.
\end{enumerate}

\section{Universal approximation on subsets}

The first result follows from existing results on approximation by network functions. The first point can be deduced from \citet[Theorem 1]{hanin2017approximating}. The second point can be derived from \citet[Proposition 2]{hanin2017approximating} or from \citet[proof of Theorem 3.2.]{hardt2016identity}.
\begin{proposition}\label{projAndFinite}
	\leavevmode
	\begin{enumerate}
		\item
		Let $M$ be some compact subset of $\mathbb{R}^{n_0}$.  If there exists a subspace $U\subset \mathbb{R}^{n_0}$ with $\dim(U)<n_0$ and such that for some $b\in\mathbb{R}^{n_0}$, $P_U$ is injective as a mapping from $M+b$ to $U$, then $M$ admits universal uniform approximation of continuous function by network functions in $\mathcal{NN}_\mathrm{ReLU}(\dim(U)+1)$.
		\item Let $M=\{x_1,...,x_m\}\subset \mathbb{R}^{n_0}$ be a finite set and $f:M\rightarrow \mathbb{R}$ some mapping. Then there exists a network function $F\in \mathcal{NN}_\mathrm{ReLU}(2)$ such that $F(x)=f(x)$ for all $x\in M$.
	\end{enumerate}
\end{proposition}
The proof is almost straight forward, given the result from the aforementioned works. A proof can be found in the supplementary material.

We next extend the second statement of the latter proposition to the case of a finite collection of disjoint $n_0$-dimensional compact sets. In Theorem \ref{maxPrinc} we will show that arbitrary uniform approximation is not admitted on such kind of sets. However, exact fit is possible for functions that are constant on each of these components under certain geometrical properties on their mutual location. The approximation of such functions arises naturally in classification tasks where each point in a particular disjoint compact set belongs to the same class, i.e. these sets constitute different disjoint clusters of the data.

In the case of functions that can only take two values, which corresponds to a two class classification, we obtain:
\begin{theorem}\label{twoClassApprox}
	Let $M=K_1\cup K_2$, with disjoint compact sets $K_1,K_2\subset \mathbb{R}^{n_0}$ such that for some $c\in\mathbb{R}^{n_0}$ and some linearly independent $v_1,...,v_{n_0}\in\mathbb{R}^{n_0}$ the set $K_1$ is contained in the sector
	\[ S:=\{x\in\mathbb{R}^{n_0}: x=c+ \sum\limits_{j=1}^{n_0} \lambda_j v_j, \, \lambda_j > 0 \}\]
	and $K_2\subset \mathbb{R}^{n_0}\setminus \overline{S}$.
	Then for every function $f:M\rightarrow \mathbb{R}$ that takes constant values on $K_1$ and $K_2$, respectively, there exists a network function $F\in\mathcal{NN}_\mathrm{ReLU}(n_0,4)$ such that
	\[f(x)=F(x)\]
	for every $x\in M$.
\end{theorem}
The assumptions on the mutual location of the compact components in Theorem \ref{twoClassApprox} is essential in the proof, but we do not know whether they can be weakened significantly. Increasing the depth will certainly allow more complex configurations. However, from Theorem \ref{maxPrinc} it will be clear that $K_2$ cannot completely enclose $K_1$ when the maximum width is upper bounded by $n_0$, no matter how many layers the network contains. The conditions of Theorem \ref{twoClassApprox} are verified and inspected by numerical experiments in Section \ref{sec:experiment}.

For the proof of the latter result we need the following.
\begin{proposition}\label{sepaCompressReLU}
	Let $K $, $M$ be two compact sets in $\mathbb{R}^{n_0}$ that are strictly separable by a linear hyper surface. Then there exists an $F\in\mathcal{NN}_\mathrm{ReLU}^{n_0}(n_0,2)$ such that $F(K)$ is a single vector in $\mathbb{R}^{n_0}\setminus M$ and $F(x)=x$ for all $x\in M$. Moreover, $F$ can be arranged in a way that given $\varepsilon>0$
	\[\min_{x\in K}\lVert F(K)-x\rVert <\varepsilon.\]
\end{proposition}
As a consequence of Proposition \ref{sepaCompressReLU} and Proposition \ref{projAndFinite}, we obtain the following for the case of more than two distinct components.
\begin{corollary}\label{finiteKompaktApprox}
	Let $M=\bigcup_{j=1}^d K_j$ with pairwise disjoint compact sets $K_j\subset\mathbb{R}^{n_0}$, $j=1,...,d$ such that for every $K_j$ there exist $v_j\in\mathbb{R}^{n_0}$ and $q_j\in\mathbb{R}$ such that
	\[v_j^T x>q_j \text{ for all } x\in K_j,\ \text{and }\  v_j^T x<q_j \text{ for all } x\in M\setminus K_j,\]
	i.e. there are linear hypersurfaces separating each $K_j$ from the remaining $K_l$, $l\neq j$. Then for every function $f:M\rightarrow \mathbb{R}$ that takes constant values on every $K_j$, $j=1,...,d$, there exists a network function $F\in\mathcal{NN}_\mathrm{ReLU}(n_0)$ with
	\[f(x)=F(x)\]
	for every $x\in M$.
\end{corollary}
A proof of this corollary is given in the supplementary material.
\begin{proof}(Proposition \ref{sepaCompressReLU})
	Let us remind that $e_1,...,e_{n_0}$ are the standard basis vectors in $\mathbb{R}^{n_0}$.
	By assumption, there is a $v\in \mathbb{R}^{n_0}$, with $\lVert v \rVert=1$ and a $q\in\mathbb{R}$ such that $v^T x>q$ for $x\in K$ and $v^T x<q$ for $x\in M$.
	For a given $\varepsilon>0$ we may assume that the linear hyper surface defined by $v,q$, $H:=\{x\in \mathbb{R}^{n_0}: v^Tx=q\}$
	is so near to $K$ that for some $a\in K$ that realizes the minimum distance of $K$ to $H$ we have
	\begin{equation}\label{KnearHyper}
		v^Ta=q+\varepsilon/2.
	\end{equation}
	Indeed, otherwise we can shift the hyper surface accordingly by increasing $q$.
	Then, by the fact that $\lVert v \rVert=1$, (\ref{KnearHyper}) means that $a$ has a distance equal to $\varepsilon/2$ to $H$ and $\tilde{a}:=a-\varepsilon/2\, v\in H$ is the unique vector in $H$ that realizes this distance.
	Let $V_1\in\mathbb{R}^{n_0\times n_0}$ such that $x\mapsto V_1x$ implements a length preserving rotation that maps $v$ to $-e_1$. Then $x\mapsto V_1x-V_1\tilde{a}$
	maps $H$ to $\mathrm{span}(e_2,...,e_{n_0})$ and $K_1:=V_1K-V_1\tilde{a}$ is contained in the half space of vectors with negative first component and with $a_1:=V_1a-V_1\tilde{a}=-\varepsilon/2\, e_1$, and  $M_1:=V_1M-V_1\tilde{a}$ is contained in the half space of vectors with positive first component. By compactness of $K_1$ and $M_1$, we can
	find $u_1,...,u_{n_0}\in \mathrm{span}(e_2,...,e_{n_0})$, such that for a  $\delta>0$, $-u_1-\delta e_1,...,-u_{n_0}-\delta e_1$ are linearly independent and such that with $\delta$ sufficiently small
	\begin{align*}
	K_1\subset S^-&:=\{x\in\mathbb{R}^{n_0}: x=\sum\limits_{j=1}^{n_0} \lambda_j (-u_j-\delta e_1), \, \lambda_j > 0 \} \\
	M_1\subset S^+&:=\{x\in\mathbb{R}^{n_0}: x=\sum\limits_{j=1}^{n_0} \lambda_j (u_j+\delta e_1), \, \lambda_j > 0 \}.
	\end{align*}
	Let $V_2\in\mathbb{R}^{n_0\times n_0}$ such that $x\mapsto V_2x$ maps $-u_j-\delta e_1$ to $-e_j$ for $j=1,...,n_0$. Then, by linearity, every $u_j+\delta e_1$ is mapped to $e_j$ for $j=1,...,n_0$. Thus $S^-$ is mapped to the cone of vectors having negative components only, and $S^+$ is mapped to the cone of vectors having positive components only. Hence, the application of $\mathrm{ReLU}$ maps all of $V_2K_1$ to $0$ and maps $V_2M_1$ identically. We can now apply the inverse of the linear affine mappings which will map $V_2M_1$ back to $M$, i.e. for
	\[F(x):=V_1^{-1}V_2^{-1}\mathrm{ReLU}(V_2(V_1x-V_1\tilde{a})) +\tilde{a}.\]
	With this $F$ we have $F(M)=M$ and $F(K)=F(a)=\tilde{a}$. Note that by (\ref{KnearHyper}) and the definition of $\tilde{a}$
	\[\min\limits_{x\in K}\lVert x - F(a)\rVert=\varepsilon/2<\varepsilon\]
	which concludes the proof.
\end{proof}

\begin{proof}(Theorem \ref{twoClassApprox})
	Let $V$ be the matrix that results from the concatenation of the columns $v_1,...,v_{n_0}$ and set $W_1:=-V^{-1}$ and $b_1:=W_1(-c)$. Then the mapping $x\mapsto W_1x $ maps the $v_j$ to $-e_j$, $j=1,...,n_0$ and $W_1K_1+b_1$ is a subset of
	\[ S^-:=\{x\in\mathbb{R}^{n_0}: x=\sum\limits_{j=1}^{n_0} \lambda_j (-e_j), \, \lambda_j > 0 \}\]
	and $W_1K_2+b_1\subset \mathbb{R}^{n_0}\setminus \overline{S^-}$. Thus, the one layer network function $F_1(x):=\mathrm{ReLU}(W_1x+b_1)$ maps all $x\in K_1$ to $0$ and $x\in K_2$ are mapped to
	\[ \{x\in\mathbb{R}^{n_0}: x=\sum\limits_{j=1}^{n_0} \lambda_j e_j, \, \lambda_j \geq 0 \}\setminus \{0\}.\]
	Hence, with a sufficiently small $q>0$, $F_1(K_1)$ and $F_1(K_2)$ are strictly separated by the linear hyper surface
	\[\{x\in\mathbb{R}^{n_0}: (1,1,....,1)x=q \}.\]
    The twofold application of Proposition \ref{sepaCompressReLU} yields two network functions that can be concatenated to a three layer network function $F_2\in \mathcal{NN}_\mathrm{ReLU}^{n_0}(n_0,3)$ with $F_2(x)=W_4\,\mathrm{ReLU}(W_3\,\mathrm{ReLU}(W_2x+b_2)+b_3)+b_4$ and such that $F_2(F_1(K_1))$ and $F_2(F_1(K_2))$ are two distinct vectors $u_1,u_2$ in $\mathbb{R}^{n_0}$, respectively. For a given $f:M\rightarrow \mathbb{R}$ that takes constant values on $K_1$ and $K_2$, say $a_1,a_2\in \mathbb{R}$, respectively, we attach a linear affine layer $x\mapsto w^T x+ b_5$, where $w\in\mathbb{R}^{n_0},\, b_5\in\mathbb{R}$ such that $w^Tu_1+b_5=a_1$ and $w^Tu_2+b_5=a_2$. This mapping is integrated in the final layer of $F_2\circ F_1$ and finally gives
	\[F(x)=w^T W_4\,\mathrm{ReLU}(W_3\,\mathrm{ReLU}(W_2\,\mathrm{ReLU}(W_1x+b_1)+b_2)+b_3)+w^Tb_4+b_5.\]
\end{proof}
In some applications it might be natural to apply periodic activation functions \citep{sitzmann2020implicit} instead of common monotonic activation functions like sigmoid or rectified linear unit. As an application of Kronecker's theorem, we obtain that approximation on finite sets can be accomplished by very narrow and shallow neural networks.
\begin{theorem}\label{cosineNet}
	Let $M=\{x_1,...,x_m\}\subset \mathbb{R}^{n_0}$, then $M$ admits universal uniform approximation of all functions $f:M\rightarrow \mathbb{R}$ by network functions in $\mathcal{NN}_{\cos}(1,3)$.
\end{theorem}
It is interesting that, in contrast to our latter result, it is shown in \citep{sontag1997shattering} that for a very wide class of parametric families, including most usual neural networks, there are at least $n$ parameters needed to shatter arbitrary sets of $2n+3$ samples. Even only $2n+1$ samples are theoretically possible if one is restricted to analytic function classes. In Theorem \ref{cosineNet} we stay below this threshold since the number of parameters we need in the proof is $n_0+3$ and this quantity is independent of the number of samples $m$. This is possible since the application of $\cos$ as activation function excludes our neural network from the function classes covered by the results in \citep{sontag1997shattering}.We further mention that there exist infinite sets, and even sets that have certain size in terms of Hausdorff measure, where Kronecker's theorem applies \citep{wik1964some}. This might give means to extent Theorem \ref{cosineNet} to an infinite sample approximation version. But it is an open question if this is possible.

\begin{proof}(Theorem \ref{cosineNet})
	The main tool in this proof is Kronecker's theorem  \citep[Th. 9.1, page 217] {katznelson2004introduction}. In order to apply this result, we have to map the $x_1,...,x_m$ to real numbers that are linearly independent over the field of rationals. To this end, let $x\mapsto W_1x+b_1$, with $W_1\in\mathbb{R}^{1\times n_0}$, $b_1\in\mathbb{R}^1$, be a suitable linear affine mapping that maps the $x_j$ to distinct $y_j:=W_1x_j+b_1$, $\{y_1,...,y_m\}\subset (0,\infty)$. With a common parameter $\alpha\in\mathbb{R}$, to be defined next, the action of the first layer is given by $\cos(\alpha (W_1 x+b_1))$.
	It is now  shown that there exists an $\alpha$ such that $\cos(\alpha y_j)$, $j=1,...m$, are linearly independent over the field of rationals. For a fixed $c=(c_1,...,c_m)^T\in\mathbb{Q}^m\setminus\{0\}$, consider
	\[T_c(\alpha)=\sum_{j=1}^m c_j\, \cos(\alpha y_j).\]
	Then $T_c$ is holomorphic in the whole complex plane, i.e. $T_c$ is an entire function, that is not constantly zero since the $y_1,...,y_m$ are non-negative and distinct. Hence, by the identity theorem for holomorphic functions, $T_c$ can only have finitely many zeros on a bounded interval, say $[0,1]$. As the set of all possible $c$ in $\mathbb{Q}^m\setminus\{0\}$ is countable, there must be an $\alpha\in[0,1]$ such that $T_c(\alpha)\neq 0$ for all $c\in\mathbb{Q}^m\setminus\{0\}$. We fix such an $\alpha$ for the determination of the first layer
	and obtain  $z_j:=\cos(\alpha (W_1 x_j+b_1))$, $j=1,...,m$ that, by construction, are linearly independent over the field of rationals. Now, let $f:M\rightarrow \mathbb{R}$ be some mapping, $\varepsilon>0$, and $\delta>0$ such that $\delta\lvert f(x_j)\rvert\leq  1$ for all $j$. Kronecker's theorem implies that we can find a $W_2\in \mathbb{R}$ such that
	\[\max_{j=1,...,m}\lvert \cos(W_2 z_j)-\delta f(x_j)\rvert <\varepsilon \delta.\]
	Thus $F(x)=1/\delta(\cos(W_2(\cos(\alpha W_1 x+ \alpha b_1))))$
	approximates $f$ with an accuracy of $\varepsilon$.
\end{proof}
\section{Maximum principle}
We will prove a maximum principle for network functions of $\mathcal{NN}_{\sigma}(n_0)$ for a wide class of activation functions. This principle can be viewed as a root cause why universal approximation with functions from $\mathcal{NN}_{\sigma}(n_0)$  on arbitrary compact sets is impossible: as shown in \citet{johnson2018deep} and \citet{beise2018decision}, see (\ref{ourEstimate}). This also immediately leads to a topological condition on subsets of $\mathbb{R}^{n_0}$ that do not admit universal approximation by functions in $\mathcal{NN}_{\sigma}(n_0)$.
\begin{theorem}[Maximum Principle]\label{maxPrinc}
	Let $M$ be some compact subset of $\mathbb{R}^{n_0}$ and $\sigma$ a continuous, monotonic  activation function. Then every $F\in\mathcal{NN}_{\sigma}(n_0)$ takes its maximum value on $\partial M$.
\end{theorem}
Note that, considering $-F$ instead of $F$, the latter result also implies that the minimum value is taken on $\partial M$.
This implies that universal uniform approximation of continuous functions is impossible on compact sets with non-empty interior. Even more, it can be concluded that:
\begin{corollary}
	Let $M$ be some compact subset of $\mathbb{R}^{n_0}$ such that $\partial D\subset M$ and $M\cap D^\circ\neq \emptyset$ for some non-empty open set $D\subset \mathbb{R}^n$. Then $M$ does not admit universal uniform approximation of continuous function by network functions from $\mathcal{NN}_{\sigma}(n_0)$.
\end{corollary}

The following proposition provides the main observation for the proof of Theorem \ref{maxPrinc}.
\begin{proposition}\label{outNotInje}
	Let $M$ be a compact subset of $\mathbb{R}^{n_0}$, $\sigma$ a continuous, monotonic activation function and $F\in\mathcal{NN}_{\sigma}^{n_0}(n_0)$ such that all weight matrices $W_j$, $j=1,...,L$, are square and have full rank. If for some $x\in M^{\circ}$ the image $F(x)$ is a boundary point of $F(M)$, i.e. $F(x)\in\partial F(M)$, then there is an $\tilde{x}\in \partial M$ such that $F(\tilde{x})=F(x)$.
\end{proposition}
\begin{proof}
	In case that $\sigma$ is injective, the whole network function $F$ is injective since the weight matrices are assumed to be square and of full rank. Hence, the invariance of the domain theorem gives that no $x\in M^\circ$ can be mapped to $\partial F(M)$.
	The remaining case is that $\sigma$ is partially constant. Then, given an $x\in M^\circ$ with $F(x)\in\partial F(M)$, let $\tilde{k}$ be the smallest $k\in\{1,...,L-1\}$ such that $\sigma(F_{\tilde{k}}(x))\in\partial \sigma(F_{\tilde{k}}(M))$. Notice that such a $\tilde{k}$ must exist since, according to the invariance of the domain theorem, the last injective, linear affine mapping $x\mapsto W_Lx$ cannot map inner points to boundary points.
	Then, necessarily for at least one component of $F_{\tilde{k}}(x)$, a partially constant part of $\sigma$ is active. More explicitly, say $(F_{\tilde{k}}(x))^{(j)}\in [a_j,b_j]\subset \mathbb{R}$ with $\sigma(t)=c_j$, where $c_j$ is the constant value taken for all $t\in [a_j,b_j]$, for some indices $j\in I\subset\{1,...,n_0\}$. In case that one of those intervals can be enlarged arbitrarily to either $[a_j,\infty)$ or $(-\infty,b_j]$, which would be the case for $\mathrm{ReLU}$ activation, we can find a $t\in\mathbb{R}$ such that $x_{\tilde{k}}=F_{\tilde{k}}(x)+t e_j\in \partial F_{\tilde{k}}(M)$ and with $\sigma(x_{\tilde{k}})=\sigma(F_{\tilde{k}}(x))$.
	In a next step we show that such an $x_{\tilde{k}}$ can also be found in cases where the length of the intervals $[a_j,b_j]$ are upper bounded. In those cases we may assume that the intervals $[a_j,b_j]$ are maximum in the sense that at their left and right end $\sigma$ is not constant any more, and $I$ is maximum in the sense that for every other index, the activation function $\sigma$ is injective in a neighbourhood. Indeed, if such components would not exist, the component wise application of $\sigma$ would be injective in a small neighbourhood of $F_{\tilde{k}}(x)$ and the invariance of the domain theorem would exclude that $\sigma(F_{\tilde{k}}(x))\in\partial \sigma(F_{\tilde{k}}(M))$. By the selection of $\tilde{k}$, $F_{\tilde{k}}(x)$ is still an inner point of $F_{\tilde{k}}(M)$. The fact that the inner point $F_{\tilde{k}}(x)\in (F_{\tilde{k}}(M))^\circ$ is mapped to a boundary point by the element wise application of $\sigma$ implies that there is some non-empty $\tilde{I}\subset I$ with corresponding $t_{j}\in\mathbb{R}$ such that
	\begin{align}
	    t_{j}\in [-(F_{\tilde{k}}(x))^{(j)}+a_j,&-(F_{\tilde{k}}(x))^{(j)}+b_j],\  j \in \tilde{I} \label{inAB} \\
	    x_{\tilde{k}}:=F_{\tilde{k}}(x)+\sum\limits_{j\in \tilde{I}} &t_j e_j\in \partial F_{\tilde{k}}(M) \label{innerOfM}.
	\end{align}
	Indeed, otherwise,
	\[F_{\tilde{k}}(x)+\sum\limits_{j=1}^{n_0}t_j e_j\in  (F_{\tilde{k}}(M))^{\circ}\]
	where for some small $\varepsilon>0$ and all $j\in I $, $t_j$ may be chosen arbitrarily in $(-(F_{\tilde{k}}(x))^{(j)}+a_j-\varepsilon,-(F_{\tilde{k}}(x))^{(j)}+b_j+\varepsilon)$ and $\varepsilon>\lvert t_j\rvert >0$ for $j\notin I$. By the definition of $I$ and the corresponding intervals $[a_j,b_j]$, this would imply that $\sigma\left(F_{\tilde{k}}(x)\right)$ is an inner point of $\sigma\left(F_{\tilde{k}}(M)\right)$. Now, (\ref{innerOfM}) implies that there is a $\tilde{x}\in M$ such that
	\[F_{\tilde{k}}(\tilde{x})=F_{\tilde{k}}(x)+\sum\limits_{j=1}^{n_0} t_j e_j \in \partial F_{\tilde{k}}(M),\]
	with $t_j$ as in (\ref{inAB}).
	Considering that $W_{\tilde{k}}$ is square and has full rank, the inverse image of $F_{\tilde{k}}(\tilde{x})$ under $x\mapsto W_{\tilde{k}}x+b_{\tilde{k}}$ is a boundary point of $\sigma (F_{\tilde{k}-1}(M))$ in case that $\tilde{k}>1$ and a boundary point of $M$ in case  when $\tilde{k}=1$. In the last case we are done. In the first case, the same reasoning can be applied to $\tilde{x}$ instead of $x$, and the iterative application finally leads to a boundary point of $M$ that is mapped to $\sigma\left(F_{\tilde{k}}(\tilde{x})\right)=\sigma\left(F_{\tilde{k}}(x)\right)$ and hence follows the same trajectory when passed through the last layers from $\tilde{k}$ to $L$.
\end{proof}
\begin{proof}(Theorem \ref{maxPrinc})
	We may assume that $M^{\circ}$ is non-empty and that every weight matrix $W_j$, $j=1,...,L-1$, is square and of full rank. Otherwise the result is trivial, because $F(M)=\partial F(M)$.
	Hence the assumption of Proposition \ref{outNotInje} holds for $\tilde{F}=A_{L-1}\circ ... \circ A_1$.
	Since the linear mapping $x\mapsto W_L x$ from $\mathbb{R}^{n_{L-1}} $ to $\mathbb{R}$ takes its maximum and minimum values on $\partial \tilde{F}(M)$, the latter proposition also implies that $F$ takes its maximum value on $\partial M$.
\end{proof}
It is not clear to us whether the conditions of Theorem \ref{maxPrinc} can be weakened in a way that it also applies in some situations to network function of width larger than the input dimensions. The following might be a natural question in this context:
Let $F:\mathbb{R}^{n_0}\rightarrow \mathbb{R}$ be a neural network function of arbitrary width but with weight matrices having rank less than or equal to $n_0$. Does a maximum principle similar to Theorem \ref{maxPrinc} apply in this case? In the supplementary material we give an example that shows that this is not case.

We conclude this section with a comment on implications of our results on generalisation.
In \citet{zhang2016understanding} the researchers made the observation that deep neural networks can interpolate randomly annotated data, which shows the enormous capabilities of large networks to memorize data. In Proposition \ref{projAndFinite}, Theorem \ref{twoClassApprox} and Corollary \ref{sepaCompressReLU} we show that network functions of width less than or equal than to input dimension in many configurations also implement an exact fit of the training data. At the same time these network functions are restricted by the maximum principle given in Theorem \ref{maxPrinc}, which may contradict the underlying distribution of the data at hand.

\section{Uniqueness}

A maximum principle similar to Theorem \ref{maxPrinc} also applies for harmonic and holomorphic functions, that is, the null space of the Laplace operator and the null space of $\partial_{\overline{z}}:=1/2(\partial_x+i\partial_y)$ ($x$ real part, $y$ imaginary part), respectively \citep{rudinRealComplex}. As a consequence, a uniqueness theorem applies which states that two harmonic, respectively holomorphic, functions coincide on the interior of a given set when they coincide on its boundary, respectively.  In view of Theorem \ref{maxPrinc}, it is hence natural to ask whether a similar result applies in the case of network functions in $\mathcal{NN}_{\sigma}(n_0)$.  However, the proof for the case of harmonic and holomorphic functions relies on the fact that these functions form a vector space, which is not the case for the set $\mathcal{NN}_{\sigma}(n_0)$.
The following examples show that a uniqueness theorem does not hold, in general, for the network functions with (strictly) monotonic activation functions considered in this work.
\begin{example}
	Let $x\mapsto\mathrm{ReLU}(-\mathrm{ReLU}(x)+1)$ and $x\mapsto\mathrm{ReLU}(-\mathrm{ReLU}(x-1)+1)$: both map $0$ to $1$ and $2$ to $0$, but they do not coincide at $1$.
\end{example}
\begin{example}
	Let $\sigma_\alpha(x)$ be the leaky ReLU function defined by $\sigma_\alpha(x)=x$ for $x\geq 0$ and $\sigma_\alpha(x)=\alpha x$ for $x<0$ with $\alpha>0$. Then
	\[x\mapsto \sigma_\alpha \left(\frac{1}{(1+\alpha)}\sigma_\alpha(x)+\frac{ \alpha }{(1+\alpha)}\right) \]
	and
	\[x\mapsto \sigma_\alpha (\frac{1}{2}\sigma_\alpha(x+1)) \]
	both map $-1$ to $0$ and $1$ to $1$, but the first one maps $0$ to $\alpha/ (1+\alpha)$, while the second one maps $0$ to $1/2$ and hence do not coincide unless $\alpha=1$.
\end{example}
A simple uniqueness result follows from the following observation.
\begin{proposition}\label{innerInj}
	Let $M$ be a compact subset of  $\mathbb{R}^{n_0}$, $F\in\mathcal{NN}_\mathrm{ReLU}^{n_0}(n_0)$ and $B=F^{-1}\left( (F(M))^{\circ}\right)$ the inverse image of the inner points of $F(M)$. Then $F$ is linear affine and bijective as a mapping from $B$ to $(F(M))^\circ$.
\end{proposition}
A proof of Proposition \ref{innerInj} is given in the supplementary material.
As an immediately consequence we have.
\begin{corollary}
	Let $F,G\in\mathcal{NN}_\mathrm{ReLU}^{n_0}(n_0)$ and $M\subset \mathbb{R}^{n_0}$ be some compact set. If $F$ and $G$ coincide on a set of vectors $v_1,...,v_{n_0+1}$  in
	\[K:=F^{-1}\left( (F(M))^{\circ}\right)\cap G^{-1}\left( (G(M))^{\circ}\right),\]
	such that $v_j-v_{n_0+1}$, $j=1,...,n_0$ are linearly independent, then $F(x)=G(x)$ for all $x\in K$.
\end{corollary}
It should be noted that in the latter corollary, the fact that $F,G$ coincide on $K$ does not imply that $F$ and $G$ are the same network functions as functions on $\mathbb{R}^{n_0}$.

\section{Numerical experiment}\label{sec:experiment}
We provide numerical results to illustrate that neural networks learn transformations similar to the construction of our theoretical derivations in the proof of Theorem \ref{twoClassApprox}. We formulate a toy example consisting of several (6 or 8) pair-wise disjoint balls of the same class encasing a center ball of a different class. The border balls use a radius of $0.125$ and a subset of the centres $(0.25, 0.25)$, $(0.5, 0.25)$, $(0.75, 0.25)$, $(0.25, 0.5)$, $(0.75, 0.5)$, $(0.25, 0.75)$, $(0.5, 0.75)$ and $(0.75, 0.75)$ while the center ball uses a radius of $ 0.01$ and the center $(0.5, 0.5)$. We generated randomly and uniformly 2000 data points for each border ball, however, as the classes would hence be unbalanced, we generated as many data points for the center ball as there are in total for the border balls. We used a multilayer perceptron (MLP) consisting of 4 layers, all of width 2, and the $\mathrm{ReLU}$ activation function (as stated in the theorem). For illustration purposes, we used an input dimension of 2 such that the neural network functions are of the form $F\in\mathcal{NN}_\mathrm{ReLU}(2)$. We used a batch size of 16, the Adam optimizer with a learning rate of 0.001 and we trained each model for 500 epochs using the mean-squared error (MSE) as our cost function. The different datasets and results are reported in Fig. \ref{fig:experiment}.  We rate an experiment as successful when the universal uniform approximation condition (\ref{conditionUniformApprox}) tends to zero. It can be observed in Fig. \ref{fig:experiment} (a) that, when the conditions of Theorem \ref{twoClassApprox} are fulfilled, the neural network can approximate the function correctly and the different layers are learning transformations similar to the constructions in the proof. When the conditions are violated, as is shown in Fig. \ref{fig:experiment} (b), the universal uniform approximation condition does not tend to zero and the intermediate layer transformations need to violate the constructions in the proof. Additional results are reported in Fig. A1 in the supplementary material.
\begin{figure}[ht]
	\vskip 0.2in
	\begin{center}
		\begin{minipage}{\textwidth}
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size/layer_0.png}\hfill
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size/layer_1.png}\hfill
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size/layer_2.png}\hfill
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size/layer_3.png}\hfill
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size/layer_4.png}\hfill
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size/layer_5.png}\hfill
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size/layer_6.png}\hfill
			\includegraphics[width=0.23\linewidth]{images/6_balls_tiny_size.png}\hfill
			\caption*{(a) 6 border balls - center ball radius: 0.01 - MSE: 0.0 - UUAC: 0.0}
		\end{minipage}
		\vskip 0.5 \baselineskip
		\begin{minipage}{\textwidth}
		    \includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size/layer_0.png}\hfill
        	\includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size/layer_1.png}\hfill
        	\includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size/layer_2.png}\hfill
        	\includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size/layer_3.png}\hfill
        	\includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size/layer_4.png}\hfill
        	\includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size/layer_5.png}\hfill
        	\includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size/layer_6.png}\hfill
        	\includegraphics[width=0.23\linewidth]{images/8_balls_tiny_size.png}\hfill
  			\caption*{(b) 8 border balls - center ball radius: 0.01 - MSE: 0.01135 - UUAC: 1.08494}
		\end{minipage}
		\caption{Datasets (a) and (b) together with the learned transformations by each activation function and layer. When the conditions of Theorem \ref{twoClassApprox} are fulfilled (a), the transformations are similar to the construction in the proof. Otherwise (b), the construction in the proof is violated. We report the final decision regions, the mean-square error (MSE) and the universal uniform approximation condition (UUAC). Each subplot's title describes the function which was applied to the previous subplot's data.}
		\label{fig:experiment}
	\end{center}
	\vskip -0.2in
\end{figure}

\section{Conclusion}

We proved that, although being not sufficient for universal approximation in the general sense, $\mathrm{ReLU}$ network functions of width less than or equal to the input dimension are sufficient to implement functions that give zero or arbitrary small training error in many machine learning tasks. In particular, we haven shown that for the case of two disjoint compact sets, exact fit in a two class classification task is feasible with 4 layers, provided the sets can be separated by a cone like sector.
However, it is not clear to what extent our conditions are sharp. We further proved that with $\cos$ as activation function, a three layer network of width 1 is sufficient to approximate any mapping on an arbitrary finite set.
We identified a maximum principle that holds in the case of width less than or equal to the input dimension for all common activation functions, which can be interpreted as a root cause why universal approximation is not possible under the aforementioned conditions. Our numerical experiments show that neural networks learn transformations similar to our theoretical derivations.

\begin{ack}
This work was partially supported by the MECO, Luxembourg project \textit{Artificial Intelligence for Safety Critical Complex Systems}.
The second author is supported by the Luxembourg National Research Fund (FNR) under the grant number 13043281
\end{ack}


\section{Introduction}\label{sect_intro}
The rapid development of communication technology makes wireless spectrum resources very tight \cite{8275026,7886242,Liu2012:Energy}.
Therefore, in recent years, cooperative communication has been paid much attention, for it helps improving spectrum utilization and system throughput in multi-user scenario.

Cooperative communication systems usually use outage probability as a metric, to measure the Quality-of-Service (QoS) and the robustness of the system.
An outage occurs when the received signal-to-noise ratio (SNR) falls below a certain threshold \cite{8435942, 2006.00664}.
In order to minimize the outage probability of the cooperative relay network, there are usually two intuitive approaches, that is, optimize relay selection scheme or power allocation scheme.

\textbf{Relay Selection Schemes:}
For a scenario with multi-relay, it is usually possible to select multiple relays to coordinate data transmission and assign orthogonal channels to avoid interference.
Jedrzejczak \textit{et al.} \cite{7438886} studied the relay selection problem by calculating of harmonic mean of channel gains.
Islam \textit{et al.} \cite{7778750} demonstrated the influence of network coverage capability on relay selection.
Das and Mehta \cite{7105886} proposed an approach for relay selection by analyzing the outage probability.
The draw back of employing too many relays is that, it may lead to extensive time consumption and frequency resource waste when forwarding signal.
To solve this issue, Bletsas \textit{et al.} \cite{1603719} proposed an opportunistic relay section scheme by only choosing the best relay according to channel state, which can obtain the full set gain.
However, these methods all assume an exact channel state information (CSI), which is not practical because of the existence of inevitable noise.

\textbf{Power Allocation Schemes:}
Based on the given relay, reasonable power allocation can further improve the received SNR and reduce the outage probability.
Given partial channel information, Wang and Chen \cite{7313006} derived a closed-form formula for optimal power allocation based on maximizing a tight capacity lower bound.
In \cite{5722051} and \cite{6108302}, the authors considered optimal power allocation scheme in different situations of convention relay and opportunistic relay.
Tabataba \textit{et al.} \cite{5618893} deduced the expression of system outage probability under the condition of high SNR, and studied power allocation using AF protocol.
However, these researches still require prior knowledge of the channel, which can not be further applied to other situations.

Recently, some researchers have successfully applied reinforcement learning (RL) methods to cooperative communication.
RL is one of the three paradigms of machine learning, which can achieve high-precision function fitting through powerful computing capability.
Unlike traditional methods, RL methods do not need prior knowledge of the environment, that is, we do not have to add any assumptions to the learning process.

In RL approaches, the source node is empowered with the learning ability to determine the optimal relay or power allocation for the current moment, based on previous observation of system state and rewards.
Shams \textit{et al.} \cite{6954557} employed Q-learning algorithm to solve the power control problem,
and Wang \textit{et al.} \cite{9072416} proposed a Q-learning based relay selection scheme in relay-aided communication scenarios.
The drawback of these studies is obvious, as they can be only suitable for simple problems with a low dimension.
Su \textit{et al.} \cite{8750861} proposed a deep Q network (DQN) based relay selection scheme with detailed mutual information (MI) as reward, but did not take power allocation into consideration.
In order to solve the joint optimization problem, Su \textit{et al.} \cite{9137340} employed convex optimization and DQN to deal with relay selection and power allocation, respectively.
However, these RL methods and frameworks can only be used in some simple environments, and do not work well in the high-dimensional search space.

In this paper, we propose a hierarchical reinforcement learning (HRL) approach for relay selection and power allocation, to minimize the outage probability of the two-hop cooperative communication system.
Unlike traditional optimization methods, our method can learn behavior policy without assuming any prior knowledge of channel state.
It is also different from existing RL-based methods that, we design an outage-based reward function, which uses a binary signal to represent success or failure of communication.
Furthermore, we propose a novel hierarchical framework to reduce searching space and improve learning efficiency.
Specifically, the contributions of this paper can be summarized as follows.
\begin{itemize}
    \item
    In our two-hop cooperative communication model, we transform the traditional outage probability optimization problem into a statistical problem, so that the RL method can be used to solve the problem.
    By employing RL methods, we no longer need to add any assumptions to channel distributions, and only rely on interaction with communication environment.

    \item
    We propose an outage-based reward function.
    Compared with other existing RL methods, our method needs less information fed back from the environment.
    Rewards are only determined by binary signals of success or failure, and do not include other concrete representations of information.
    It is practical because other additional feedback may be not available in certain situations.

    \item
    We further design an HRL framework with two levels for cooperative relay network, where relay selection and power allocation are disassembled into two optimization objectives.
    Traditional deep reinforcement learning (DRL) methods considers relay selection and power allocation together, which leads to a more complex action space and may affected the learning performance.
    By decomposing different optimization objectives into different levels, complex action space is therefore simplified in our HRL framework.

\end{itemize}

The rest of this paper is organized as follows.
Section \ref{sect preli} introduces the preliminaries of DRL.
Section \ref{sect model} analyzes our system model and Section \ref{sect problem} formulates the outage minimization problem.
Section \ref{sect DRL} describes our outage-based method using DQN framework.
Section \ref{sect HRL} describes our proposed HRL framework and learning algorithm, and presents our pre-training algorithm in detail.
Section \ref{sect Result} presents simulation results.
Finally, Section \ref{sect Conclusion} concludes this paper and outlines future works.

\section{Preliminaries}\label{sect preli}
In the field of cooperative communication, recent studies have proposed several machine learning methods for relay selection or power allocation.
Traditional communication methods make assumptions about CSI, while these methods use data for learning and then making channel predicting.

As one of the three paradigms of machine learning, RL is an emerging tool to solve decision-making problems such as resource management in communication \cite{8931561,8761525,8771176}.
RL methods use an agent, which can be regarded as an intelligent robot, to interact with the environment.
The agent in RL methods has no access to prior knowledge of the environment, but can only get familiar with the environment through an interactive process called Markov Decision Process (MDP).
In order to minimize outage probability, the agent will repeatedly interact with the communication environment, and choose a suitable relay and allocate transmission power according to current system state.
In addition, it continuously adjusts its behavior policy according to the feedback from the communication environment.
In this section, we introduce these preliminaries of reinforcement learning.

\subsection{Markov Decision Process}\label{subsect_preliA}
An MDP consists of an environment $\mathcal{E}$, a state space $\mathcal{S}$, an action space $\mathcal{A}$, and a reward function $\mathcal{S}\times\mathcal{A}\to\mathcal{R}$. At each discrete time step $t$, the agent observes the current state $s_t\in\mathcal{S}$, and selects an action $a_t\in\mathcal{A}$ according to a policy $\pi$: $\mathcal{S}\to \mathcal{P}(\mathcal{A})$, which maps states to a probability distribution over actions. After executing action $a_t$, the agent receives a scalar reward $r(s_t,a_t)$ from the environment $\mathcal{E}$ and observes the next state $s_{t+1}$ according to the transition probability $p(s_{t+1}|s_t,a_t)$. This process will continue until a terminal state is reached.

The goal of the agent is to find the optimal policy to maximize the expected long-term discounted reward, \emph{i.e.}, maximize the expected accumulated return $R_t=\sum\nolimits_{i=t}^T \gamma^{i-t}r(s_i,a_i)$ from each state $s_t$, where $T$ denotes the total step, and $\gamma\in[0,1]$ denotes the discount factor that trades off the importance of immediate and future rewards.

Action-value function is usually used to describe the expected return after selecting action $a_t$ in state $s_t$ according to policy $\pi$.
\begin{equation}\label{eq2_1} Q^{\pi}(s_t,a_t)=\mathbb{E}_{s_t\sim\mathcal{E},a_t\sim\pi(s_t)}(R_t|s_t,a_t). \end{equation}
And we can obtain the preceding action-value function via recursive relationship known as Bellman function.
\begin{equation}\label{eq2_2} \begin{aligned}
Q^{\pi}(s_t,a_t)=\mathbb{E}_{s_{t+1}\sim\mathcal{E}}\Big[r(s_t,a_t)+\gamma\mathbb{E}_{a_{t+1}\sim\pi(s_{t+1})}\big[Q^{\pi}(s_{t+1},a_{t+1})\big]\Big].
\end{aligned}\end{equation}
Moreover, the optimal action-value function $Q^{\ast}(s_t,a_t)=\max_{\pi\in\Pi}Q^{\pi}(s,a)$ gives the maximum action value under state $s$ and action $a$, and it also obeys Bellman function.
\begin{equation}\label{eq2_3} Q^{\ast}(s,a)=\mathbb{E}_{s^{\prime}\sim\mathcal{E}}[r+\gamma\max\limits_{a^{\prime}}Q^{\pi}(s^{\prime},a^{\prime})]. \end{equation}

\subsection{Reinforcement Learning}\label{subsect_preliB}
In practice, we usually do not know the underlying state transition probability, \textit{i.e.}, in a model-free situation. It requires the agent to interact with the environment and learn from the feedback, constantly adjust its behavior to maximize the expected reward.

In reinforcement learning, temporal difference (TD) methods \cite{RL:intro2} are proposed via combining Monte Carlo methods and dynamic programming methods, which enable the agent to learn directly from raw experience. Therefore, we have the following well-known Q-learning algorithm.
\begin{equation}\label{eq2_5} Q(s_t,a_t)\leftarrow Q(s_t,a_t)+\alpha\delta_t \end{equation}
with
\begin{equation}\label{eq2_4} \delta_t=r(s_t,a_t)+\gamma\max\limits_{a_{t+1}\in\mathcal{A}}Q(s_{t+1},a_{t+1})-Q(s_t,a_t), \end{equation}
where $\delta_{t}$ denotes TD error and $\alpha\in[0,1]$ denotes learning rate.

Through continuous iterative updating, the Q value of different actions selected in each state finally tends to be stable, which can then provide a policy for the subsequent action selection.

\section{System Model}\label{sect model}
Consider a wireless network where exists an $N_S$-antenna source $S$, an $N_D$-antenna destination $D$, and a group of single-antenna relays $R=\{R_1,R_2,\dots,R_K\}$, as shown in Fig. \ref{model}.
We assume that the source is far from the destination, and the help of relay nodes is needed.
Due to the limitation of equipment of relays, we consider a half-duplex signaling mode where the communication from $S$ to $D$ via the selected relay $R_i$ will take two time slots.
In the first time slot, $S$ broadcasts its signal, and all the other nodes, include the destination, listen to this transmission.
In the second time slot, the selected relay forwards the decoded signal to destination.
  \begin{figure}[ht]
    \centering

    \includegraphics[scale=0.75]{fig2/Relay_new_2.pdf}
    \caption{Relay network model.}
    \label{model}
  \end{figure}

Depending on how the cooperative relay processes the received signal, the relay mode can be mainly divided into amplify-and-forward (AF) and decode-and-forward (DF).
Next, we will analyze the MI obtained by using AF protocol and DF protocol, respectively.

\subsection{Amplify-and-Forward Relaying}\label{subsect_AF_Mode}
In the case that all relays can only scale the received signal and send it to the destination, we employ AF protocol to realize cooperative communication. In the first phase, the received signal at $R_i$ can be written as
\begin{equation}\label{eq3_1} y_{si}(t)=\sqrt{P_s}\boldsymbol h_{si}^{\dagger}(t) \boldsymbol x(t)+n_i(t), \end{equation}
where $P_s \in [0,P_{max}]$ is the transmission power at source with $P_{max}$ being the maximum value, $\boldsymbol x(t)$ is a $N_s \times 1$ data symbol vector and $\|\boldsymbol x(t)\|=1$,
$\boldsymbol h_{si}(t)=[h_{si}^1(t),\dots,h_{si}^{N_S}(t)]^T$ is an $N_S\times 1$ channel vector between source and relay,
and each element is a complex Gaussian random variable with zero mean and variance $\sigma_{si}^2$, $n_i(t)\sim\mathcal{CN}(0,\sigma_n^2)$ is the complex Gaussian noise at relay. Similarly, we have the received signal at $D$ which can be written as
\begin{equation}\label{eq3_2} \boldsymbol y_{sd}(t)=\sqrt{P_s}\boldsymbol h_{sd}^{\dagger}(t)\boldsymbol x(t)+\boldsymbol n_d(t), \end{equation}
where $\boldsymbol h_{sd}(t)$ denotes a $N_S \times N_D$ channel matrix, and $\boldsymbol n_d(t)\sim \mathcal{CN}(\boldsymbol 0,\sigma_n^2\boldsymbol I_{N_D})$ denotes the complex Gaussian noise at destination, where $\boldsymbol I$ is identity matrix.

In the second time slot, the selected relay amplifies the signal and transmits it to the destination.
The received signal at the destination from the relay can be written as
\begin{equation}\label{eq3_3} \begin{aligned}
&\boldsymbol y_{id,AF}(t)=\sqrt{P_r}\boldsymbol h_{id}^{\dagger}(t)\beta y_{si}(t)+\boldsymbol n_d(t), \\
\end{aligned}\end{equation}
where $P_r\in [0,P_{max}]$ is the transmission power at relay with $P_{max}$ being the maximum value,
$\boldsymbol h_{id}(t)=[h_{id}^1(t),\dots,h_{id}^{N_D}(t)]$ is a $1 \times N_D$ channel vector between relay and destination,
and similarly, each element is a complex Gaussian random variable with zero mean and variance $\sigma_{id}^2$.
$\beta$ is the amplification factor, which can be written as follows \cite{5710995, 1427716}.
\begin{equation}\label{eq3_4} \beta =\sqrt{\frac{1}{P_s\|\boldsymbol h_{si}\|^2+\sigma^2_n}}. \end{equation}

The destination combines the data from the source and the relay using maximal ratio combining (MRC), and after some manipulations according to \cite{1350931,4107949}, we have the following final end-to-end SNR.
\begin{equation}\label{eq3_5} \varphi_{z} =\frac{\varphi_{si}\varphi_{id}}{\varphi_{si}+\varphi_{id}+1}, \end{equation}
where $\varphi_{si}=P_s{\|\boldsymbol h_{si}\|}^2/\sigma_n^2$ and $\varphi_{id}=P_r{\|\boldsymbol h_{id}\|}^2/\sigma_n^2$.

Similarly, we can obtain the SNR of direct transmission from source to destination, which can be represented as $\varphi_{sd}=P_s{\|\boldsymbol h_{sd}\|}^2/\sigma_n^2$. Then we have the MI between the source and the destination using AF protocol.
\begin{equation}\label{eq3_7} \begin{aligned}
I_{AF}=\frac{1}{2}\log_2(1+\varphi_{AF})=\frac{1}{2}\log_2(1+\varphi_{sd}+\frac{\varphi_{si}\varphi_{id}}{\varphi_{si}+\varphi_{id}+1}).
\end{aligned}\end{equation}

\subsection{Decode-and-Forward Relaying}\label{subsect_DF_Mode}
Assume that all relays are able to decode the signal from the source, and then re-encode and transmit the signal to the destination. The first time slot in DF mode is the same as that in AF mode, and (\ref{eq3_1}) and (\ref{eq3_2}) have given the received signal at $R_i$ and $D$ in this time slot.

In the second time slot, different from that in AF mode, the selected relay decodes and forwards the signal to the destination, and the received signal at the destination from the relay can be written as
\begin{equation}\label{eq3_8} \boldsymbol y_{id,DF}(t)=\sqrt{P_r}\boldsymbol {h_{id}}^{\dagger}(t) y_{si}(t)+\boldsymbol n_d(t), \end{equation}

When employing DF protocol, the relays need to first successfully decode the signal from source on the condition that the MI is above the required transmission rate. Then, the signals received in both two time slots are combined at the destination using MRC \cite{4107949,1603719}, and we have the following instantaneous MI between the source and the destination.
\begin{equation}\label{eq3_9} \begin{aligned}
I_{DF}=\frac{1}{2}\log_2(1+\varphi_{DF})=\frac{1}{2}\log_2(1+\varphi_{sd}+\varphi_{id}).
\end{aligned}\end{equation}

\section{Problem Formulation}\label{sect problem}
Suppose that there is an agent in the communication environment, which has access to channel state in previous time slots.
The agent estimates current channel state based on historical CSI, and accordingly selects relay and allocates transmission power.
Afterwards, it receives a reward from the environment, which indicates whether the communication is successful.

In this section, we model this process as an MDP, where historical channel state is regarded as system state, and relay selection along with power allocation are considered as system action.
Then we describe the variables in our two-hop cooperative communication scenario and formulate our optimization problem.

\subsection{State Space}\label{subsect state}
Full observation of our two-hop communication system consists of the channel states between any two nodes in the previous time slot. Therefore, the state space in current time slot is a union of different wireless channel states, which can be denoted as
\begin{equation}\label{eq4_1} \mathcal{S}_t \triangleq [\boldsymbol h_{si}(t-1), \boldsymbol h_{id}(t-1), \boldsymbol h_{sd}(t-1)], \end{equation}
where the integer $i$ satisfies $i \in [1,K]$.

In order to characterize the temporal correlation between time slots for each channel, we employ the following widely adopted Gaussian Markov block fading autoregressive model \cite{5710995, ref1}.
\begin{equation}\label{eq4_2} \boldsymbol h_{ij}(t)=\rho\boldsymbol h_{ij}(t-1)+\sqrt{1-\rho^2}\boldsymbol e(t), \end{equation}
where $\rho$ denotes the normalized channel correlation coefficient between corresponding elements in $\boldsymbol h(t)$ and $\boldsymbol h(t-1)$, $\boldsymbol e(t)\sim\mathcal{CN}(\boldsymbol 0,\sigma^2\boldsymbol I)$ denotes the error variable and is uncorrelated with $\boldsymbol h_{ij}(t)$. According to Jake's fading spectrum, we have $\rho=J_0(2\pi f_d\tau)$ where $J_0(\cdot)$ denotes zeroth-order Bessel function of the first kind, $f_d$ and $\tau$ denote Doppler frequency and the length of time slot, respectively.

\subsection{Action Space}\label{subsect action}
Full action space includes relay selection $\boldsymbol a^{R}(t)$, source power allocation $\boldsymbol a^{P_s}(t)$, and relay power allocation $\boldsymbol a^{P_r}(t)$.

Considered that the total power is constraint, \textit{i.e.}, $P_{s}+P_{r} \leq P_{max}$, we can assume that the sum of power used by the source and its selected optimal relay is $P_{max}$. So $P_r$ can be directly represented by the difference between $P_{max}$ and $P_s$. Then, we can reduce the number of actions that need to be optimized and derive the following reduced action space.
\begin{equation}\label{eq4_3 action space} \mathcal{A}_t \triangleq [\boldsymbol a^{R}(t), \boldsymbol a^{P_s}(t)]. \end{equation}

The first part of action space is relay selection, which is denoted by
\begin{equation}\label{eq4_3 relay selection} \boldsymbol a^{R}(t)=[a_1^R(t),a_2^R(t),\dots,a_K^R(t)], \end{equation}
where $a_k^R(t)=1$ means relay $R_k$ is selected in time slot $t$ and $a_k^R(t)=0$ otherwise.

The second part of is power allocation for the source node. Similarly, it is denoted by
\begin{equation}\label{eq4_3_power allocation} \boldsymbol a^{P_s}(t)=[a_1^{P_s}(t),a_2^{P_s}(t),\dots,a_{L-1}^{P_s}(t)], \end{equation}
where $P_{max}$ is divided into $L$ power-levels, and $a_l^{P_s}(t)=1$ means the $l$-th power-level is selected for source node transmission in time slot $t$, and $a_l^{P_s}(t)=0$ otherwise.

  \begin{figure*}[t]
    \centering

    \includegraphics[scale=0.8]{fig2/DQN.pdf}
    \caption{DRL framework for relay selection and power allocation.}
    \label{DQN}
  \end{figure*}

\subsection{Reward and Optimization Problem}\label{subsect reward}
The outage probability minimizing problem that jointly optimizes relay selection and power allocation can be intuitively formulated as
\begin{equation}\label{eq_problem_1}
    \min\limits_{\mathcal{A}_t} \ {\rm \textbf{Prob}}(I<\lambda),
\end{equation}
where the positive scalar $\lambda$ is denoted to be the outage threshold, and $I$ can be calculated according to (\ref{eq3_7}) and (\ref{eq3_9}).

Traditional methods establish a probabilistic model, where the distribution employed to describe channel uncertainty is assumed artificially.
However, we do not rely on underlying channel distributions in this paper, and thus traditional probabilistic analysis methods are not applicable.

In our problem, the agent can only make use of the communication result that denotes success or failure from the cooperative communication environment.
Therefore, we define the following indicator function of event $I<\lambda$, which represents the result after each selection.
\begin{equation}\label{eq4_4_indicator}
    f(\boldsymbol a^{R}, \boldsymbol a^{P_s}; \boldsymbol h) \triangleq \mathbbm{1}_{I<\lambda}=\left\{
        \begin{aligned}
            &1, &if\ I<\lambda\ \ \\
            &0, &otherwise
        \end{aligned}
    \right.
\end{equation}

Consider the fact that, when an indicator function is employed to represent each occurrence of an event, then the expectation of the indicator function can be used to calculate the probability of the original event.
Therefore, we can reformulate the optimization problem for minimizing outage probability of our communication system in the form of statistics.
Then, this problem can be solved using RL methods, which is formulated as follows.
\begin{equation}\label{eq4_4_problem} \begin{aligned}
\min\limits_{\mathcal{A}_t} \ &\mathbb{E} \left[ \frac{1}{T}\sum\limits_{t=1}^T f\big(\boldsymbol a^{R}(t), \boldsymbol a^{P_s}(t);\boldsymbol h(t)\big) \right] \\
s.t.\
&\textbf C_1: \sum_{k=1}^K a_k^R(t)=1, \\
&\textbf C_2: \sum_{l=1}^{L-1} a_l^{P_s}(t)=1, \\
&\textbf C_3: a_k^R(t),a_l^{P_s}(t) \in \{0,1\}.
\end{aligned}\end{equation}

In MDP, the reward is fed back to the agent to evaluate the selected action under current system state.
In this paper, we design an outage-based reward function, which only consists of the communication result.
Note that, the goal of the agent is to find the optimal behavior policy to maximize the expected long-term discounted reward, so our binary reward function is denoted as
\begin{equation}\label{eq4_4_reward}
r_t = 1-f\big(\boldsymbol a^{R}(t), \boldsymbol a^{P_s}(t);\boldsymbol h(t)\big).
\end{equation}

\section{DRL Based Solution}\label{sect DRL}
In reinforcement learning, we often estimate the action-value function by using Bellman function as an iterative update to converge to the optimal.
Unfortunately, traditional RL method which only employs Bellman function has no generalization ability.
Since channel state is uncountable, traditional RL method often fails to make decisions when faced with channel state that has never appeared.
DRL is a combination of deep neural networks (DNN) and traditional RL method, which is proposed to solve generalization problem in large state or action space \cite{Mnih2015, IMPALA}.
In this section, we adopt DQN framework to design a DRL-based solution for relay selection and power allocation.

The DRL framework for relay selection and power allocation is shown in Fig. \ref{DQN}.
Note that the source node has no prior knowledge of the communication system, which means the distributions of wireless channels between any two nodes are all unknown to it.
The agent can only observe the current state from the communication environment and get state $s_t$.
Then a deep neural network is employed as the nonlinear function approximator to deal with the input data, which can be applied to estimate the action-value function for high dimensional state space.
According to (\ref{eq2_3}), we have
\begin{equation}\label{eq4.5_5}
Q^{\pi}(s_t, a_t; \theta)=\mathbb{E}\Big[r_{e,t}+\gamma_l\max\limits_{a_{t+1}}Q^{\pi}(s_{t+1}, a_{t+1}; \theta)\Big].
\end{equation}

After calculation, the agent chooses the best action which can induce the maximal value of $Q^{\pi}(s_t, a_t; \theta)$, and then the environment will give the corresponding reward and update system state.
So far, we have obtained a complete experience tuple $e_t=(s_t,a_t,r_{t},s_{t+1})$, which will be stored in experience replay buffer $\mathcal{B}=\{e_1,e_2,\dots,e_t\}$.
When training, a batch of experience will be sampled and used to optimize a set of loss functions below.
\begin{equation}\label{eq4.5_6} L_i(\theta_i)=\mathbb{E}_{e_t\sim\mathcal{B}}\Big[\big(y_i-Q(s_t, a_t; \theta_i)\big)^2\Big], \end{equation}
with
\begin{equation}\label{eq4.5_7} y_i=r + \gamma_l \max\limits_{a_{t+1}}Q^{\pi}(s_{t+1}, a_{t+1}; \theta^{-}), \end{equation}
where $\theta^{-}$ is parameters from previous iteration in a separate target network which are held fixed when optimizing, and will be replaced by $\theta_{i-1}$ from the evaluate network after a period of time.

Then differentiate operation on these loss functions will be carried out, and we can yield the following expression.
\begin{equation}\label{eq4.5_8} \begin{aligned}
\nabla_{\theta_i}L_i(\theta_i)=\mathbb{E}\Big[\big(y_i-Q(s_t,a_t;\theta_i)\big) \nabla_{\theta_i}Q(s_t,a_t;\theta_i)\Big].
\end{aligned}\end{equation}

The following standard non-centered RMSProp optimization algorithm \cite{A3C,lecture} is then adopted to minimize the loss function and update parameters in Q network .
\begin{equation}\label{eq4.5_9} \theta\leftarrow\theta-\eta\frac{\Delta\theta}{\sqrt{\upsilon+\epsilon}}, \end{equation}
with
\begin{equation}\label{eq4.5_10} \upsilon=\kappa \upsilon +(1-\kappa)\Delta\theta^2, \end{equation}
where $\kappa$ is a momentum and $\Delta\theta$ is the accumulated gradients.

Note that, this is a model-free approach, as the agent using state and reward sampled from the environment rather than estimating transition probability. And this is an offline policy, because an epsilon greedy method will be employed as the behavior policy.
For a two-hop cooperative relay network, Algorithm \ref{algoDQN} employs DQN framework to make dynamic relay selection and power allocation.
In the Evaluation part of this paper, we will test the algorithm and use it for comparison.
Pseudocode of the algorithm can be found in Algorithm \ref{algoDQN}.

\begin{algorithm}[htb]
    \caption{DRL Based Relay Selection and Power Allocation}
    \label{algoDQN}
    \begin{algorithmic}[1]
        \STATE Initialize experience replay buffer $\mathcal{B}$.
        \STATE Initialize Q network with random weights $\theta$.
        \STATE Initialize target Q network with $\theta^{-}=\theta$.
        \FOR{episode $u=1,2,\dots,u_{max}$}
            \STATE Initialize the environment, get state $s_1$.
            \FOR{time slot $t=1,2,\dots,t_{max}$}
                \STATE Choose action $a_t$ using epsilon-greedy method with a fix parameter $\epsilon$.
                \STATE Execute action $a_t$, and observe reward $r_{t}$ and next state $s_{t+1}$.
                \STATE Collect and save the tuple $e_t$ in $\mathcal{B}$.
                \STATE Sample a batch of transitions $(s_j,a_j,r_{j},s_{j+1})$ from $\mathcal{B}$.
                \IF {episode terminates at time slot $j+1$}
                    \STATE Calculate $y_j$ according (\ref{eq4.5_7});
                \ELSE
                    \STATE Set $y_j=r_{j}$.
                \ENDIF
                \STATE Perform gradient descent and update Q-network according to (\ref{eq4.5_8}).
                \STATE Update current state, and every $C$ steps reset $\theta^{-}=\theta$.
            \ENDFOR
        \ENDFOR
    \end{algorithmic}
\end{algorithm}

\section{HRL Based Solution}\label{sect HRL}
Traditional DRL-based approaches put all variables together in its action, which results in a complex search space.
HRL is a recent technology based on DRL, which has developed rapidly in recent years and is considered a promising method to solve problems with sparse rewards in complex environment.
HRL enables more efficient exploration of the environment by abstracting complex tasks into different levels \cite{HDQN, HIRO, 8629360}.
In this section, we propose a novel two-level HRL framework for cooperative communication, to learn relay selection policy and power allocation policy in different levels.

  \begin{figure*}[ht]
    \centering

    \includegraphics[scale=0.85]{fig2/Framework.pdf}
    \caption{HRL framework for decomposing relay selection and power allocation into different levels.}
    \label{Framework}
  \end{figure*}

\subsection{Proposed Framework}\label{subsect_HRLA}

As shown in Fig. \ref{Framework}, the communication agent has two levels. In the higher level, the meta-controller receives observation of state from the external communication environment, and outputs a goal $g_t$. The controller in the lower level is supervised with the goals that are learned and proposed by meta-controller, which observes the state of the external communication environment and selects an action $a_t$. Note that the meta-controller gives a goal every $n$ steps, and the goal will remain until the low-level controller reaches the terminal. We employ a standard experience replay buffer, and it is worth to note that, experience tuples $(s_t,g_t,r_{e,t},s_{t+n})$ for meta controller and $(s_t,g_t,a_t,r_{i,t},s_{t+1})$ for controller are stored in disjoint spaces for training.

Hierarchical environment includes state, high-level goal, low-level action and reward, which are described specifically below.

\textbf{State:}
State in our hierarchical framework is the same as that in DRL environment, which consists of the channel states between any two nodes in the previous time slot. The expression for the state space can be referred to (\ref{eq4_1}).

\textbf{High-Level Goal:}
In cooperative communication systems, we can intuitively find that relay selection plays a major role.
Therefore, we separate the different action components to make different levels, and extract relay selection as high-level goals for overall planning.
Denote $g_t$ to be the goal in higher level, we then have
\begin{equation}\label{eq5_1} g_t \in \mathcal{G} \triangleq [\boldsymbol g^R(t)]=[g_1^R(t),g_2^R(t),\dots,g_K^R(t)]. \end{equation}

In fact, the goal selection in high-level is similar to relay selection action in the previous DRL method. Therefore, $\boldsymbol g^R(t)$ should meet the same constraint that $\sum_{k=1}^K g_k^R(t)=1,\ g_k^R(t) \in \{0,1\}$.

\textbf{Low-Level Action}:
By decomposing relay selection and power allocation into different levels, we can further reduce the action space. Then the low-level action space only has one variable $\boldsymbol a^{P_s}(t)$, which satisfies $C_2$ and $C_3$ in (\ref{eq4_4_problem}).

\textbf{Reward}:
Note that the higher level and the lower level are working in different time scales. Meta-controller first proposes a temporarily fixed goal for the lower level, and then controller performs actions in a period of time according to both system state and high-level goal and receive feedbacks from the environment.
Therefore, we can denote the internal reward for low-level controller as
\begin{equation}\label{eq5_2} r_{i,t}=1-f\big(\boldsymbol g^{R}(t), \boldsymbol a^{P_s}(t);\boldsymbol h(t)\big). \end{equation}

On the other hand, we use communication success rate of a given relay over a period of time $n$ to measure the quality of current relay selection. Therefore, the external reward for high-level meta-controller can be represented as
\begin{equation}\label{eq5_3} r_{e,t}=\frac{1}{n}\sum\limits_{t=1}^{n} r_{i,t}, \end{equation}
which the agent aims to maximize its expectation.

\newcounter{TempEqCnt}

\subsection{Hierarchical Learning Policy}\label{subsect_HRLB}
For meta-controller in higher level, we use gradient bandit method to learn goal-policy for dynamically proposing goals according to a given system state.
Recall that, we have $K$ relays to choose from.
Therefore, we first establish the following probability distribution.
\begin{equation}\label{eq5_4} \pi^h_t(R_i) \triangleq Pr\{g_t=R_i\} \triangleq \frac{e^{M_t(R_i)}}{\sum_{b=1}^{K} e^{M_t(R_b)}}, \end{equation}
where $\pi^h$ is the high-level policy, and $\pi^h_t(R_i)$ denotes the probability that relay $R_i$ is selected as the goal in time slot $t$. $M_t(R_b)$ denotes the preference value for choosing relay $R_b$, which will be updated every $n$ steps.

Then, we employ stochastic gradient descent to update the preference values.
\begin{equation}\label{eqappendix_1}
M_{t+n}(R_i) \triangleq M_t(R_i) + \zeta\frac{\partial\mathbb{E}[r_{e,t}]}{\partial{M_t(R_i)}},
\end{equation}
where $\zeta>0$ denotes learning step size, and the expectation of $r_{e,t}$ can be equally calculated by $\sum_{b=1}^K\pi^h_t(R_b)Q^{\ast}_h(R_b)$.
Replace the expectation form in equation (\ref{eqappendix_1}), and we then have
\begin{equation}\label{eqappendix_3} \begin{aligned}
\frac{\partial\mathbb{E}[r_{e,t}]}{\partial{M_t(R_i)}} &= \frac{\partial}{\partial{M_t(R_i)}}\left[\sum\limits_{b=1}^K\pi^h_t(R_b)Q^{\ast}_h(R_b)\right] \\
&= \sum\limits_{b=1}^K Q^{\ast}_h(R_b)\frac{\partial\pi^h_t(R_b)}{\partial{M_t(R_i)}} \\
&= \sum\limits_{b=1}^K (Q^{\ast}_h(R_b)-\bar{r}_{e,t})\frac{\partial\pi^h_t(R_b)}{\partial{M_t(R_i)}} \\
&= \sum\limits_{b=1}^K \pi^h_t(R_b)\frac{(Q^{\ast}_h(R_b)-\bar{r}_{e,t})}{\pi^h_t(R_b)}\frac{\partial\pi^h_t(R_b)}{\partial{M_t(R_i)}}, \\
\end{aligned}\end{equation}
where the newly introduced scalar $\bar{r}_{e,t}$ is independent of $b$. It denotes the average of all external rewards, \textit{i.e.}, the average success rate of our communication system. Further, the partial derivative part can be further written as
\begin{equation}\label{eqappendix_4} \begin{aligned}
\frac{\partial\pi^h_t(R_b)}{\partial{M_t(R_i)}} &= \frac{\partial}{\partial{M_t(R_i)}}\left[\frac{e^{M_t(R_i)}}{\sum_{b=1}^{K} e^{M_t(R_b)}}\right] \\
&= \frac{\mathbbm{1}_{i=b}e^{M_t(R_i)}\sum_{b=1}^{K} e^{M_t(R_b)}-e^{M_t(R_i)}e^{M_t(R_b)}}{\Big(\sum_{b=1}^{K} e^{M_t(R_b)}\Big)^2} \\
&= \mathbbm{1}_{i=b}\pi^h_t(R_b)-\pi^h_t(R_b)\pi^h_t(R_i).
\end{aligned}\end{equation}

Note that, $\mathbb{E}[r_{e,t}|g_t]=Q^{\ast}_h(g_t)$ and $r_{e,t}$ is independent of the other variables.
Then we can derive the following equation by employing the form of expectation.
\begin{equation}\label{eqappendix_5} \begin{aligned}
\frac{\partial\mathbb{E}[r_{e,t}]}{\partial{M_t(R_i)}} &= \sum\limits_{b=1}^K \pi^h_t(R_b)\big(Q^{\ast}_h(R_b)-\bar{r}_{e,t}\big)\big(\mathbbm{1}_{i=b}-\pi^h_t(R_i)\big) \\
&= \mathbb{E}\Big[(r_{e,t}-\bar{r}_{e,t})\big(\mathbbm{1}_{i=b}-\pi^h_t(R_i)\big) \Big].
\end{aligned}\end{equation}

In the training process, sampling is conducted every $n$ time steps, and the gradient in (\ref{eqappendix_1}) is replaced by the expectation value of the single sample. Therefore, we can finally obtain the following update expression of preference value.
\begin{equation}\label{eqappendix_6}
M_{t+n}(R_i) \triangleq M_t(R_i) + \zeta(r_{e,t}-\bar{r}_{e,t})\big(\mathbbm{1}_{i=b}-\pi^h_t(R_i)\big).
\end{equation}

For controller in lower level, it learns action-policy for selecting actions according to both state and goal, which aims to maximize the long-term discounted expected internal reward.

In order to reflect the difference in values of the different actions, we perform some changes to the architecture of traditional deep Q network.
Inspired by \cite{Duel}, we further employ a dueling network, which can enhance the stability of DRL algorithm by ignoring subtle changes in the environment and focusing on key states.

Schematic illustration of dueling architecture is shown in Fig. \ref{DuelingDNN}.
The input layer and hidden layers are the same as that in traditional deep Q network.
The key difference is that, there is a sub-output layer in our dueling network, where the traditional Q function output is separated into a state-goal valuation function $V(s_t,g_t)$ and an advantage evaluation function $A(s_t,g_t,a_t)$.
In state-goal valuation part, there is only one neuron which represents the assessment of the current state and goal.
In advantage evaluation part, the number of neurons is equal to that in output layer, representing advantage of choosing each optional action.
  \begin{figure}[ht]
    \centering

    \includegraphics[scale=1.0]{fig2/DNN.pdf}
    \caption{Dueling network in controller for low-level learning.}
    \label{DuelingDNN}
  \end{figure}

Consider the fact that $\mathbb{E}_{a_t\sim\pi_l}[Q(s_t,g_t,a_t)]=V(s_t,g_t)$, we have $\mathbb{E}_{a_t\sim\pi_l}[A(s_t,g_t,a_t)]=0$. Therefore, to meet this property, we can rewrite the advantage part as
\begin{equation}\label{eq_5_5_adv}
A(s_t,g_t,a_t) \triangleq A(s_t,g_t,a_t;\theta,w_2)-\frac{1}{|A|}\sum\limits_{a^{\prime}}A(s_t,a^{\prime};\theta,w_2),
\end{equation}
and thus have the following expression of Q function by combining the two parts in sub-output layer.
\begin{equation}\label{eq_5_5} \begin{aligned}
Q_l(s_t,g_t,a_t) \triangleq V(s_t,g_t;\theta,w_1)+ A(s_t,g_t,a_t;\theta,w_2),
\end{aligned}\end{equation}
where $\theta$ denotes parameters in common part of DNN (\textit{i.e.}, three columns on the left in Fig. \ref{DuelingDNN}), and $w_1$ and $w_2$ are parameters in separated fully connected  sub-output layer for valuation function and advantage function, respectively. Note that, the output of our dueling network is still the same as that of the traditional network, that is, the estimated expected return for each action $a \in \mathcal{A}$ under the current state $s_t$ and goal $g_t$.

By employing temporal difference method, the optimal value of Q function can be written as
\begin{equation}\label{eq5_6} \begin{aligned}
Q_l^{\ast}(s_t,g_t,a_t;\theta,w_1,w_2)
&=\max\limits_{\pi^l_t}\mathbb{E}\left[\sum\limits_{j=t}^{\infty}\gamma_l^{j-t}r_{i,j}\right] \\
&=\max\limits_{\pi^l_t}\mathbb{E}\Big[r_{i,j} + \gamma_l \max\limits_{a_{t+1}}Q_l^{\ast}(s_{t+1},g_t,a_{t+1};\theta,w_1,w_2)\Big],
\end{aligned}\end{equation}
where
$\pi_t^l$ denotes low-level policy for power allocation in time slot $t$. Note that meta-controller and controller work on different timescales. The controller operates at each time step, while the meta-controller operates on a longer timescale of $n$ time steps.

When training, a batch of memories are sampled from experience replay buffer $B_l$. We calculate a set of loss functions, and derive a batch of gradients as
\begin{equation}\label{eq5_7} \begin{aligned}
\nabla_{\theta_l}L_l(\theta_l)=\mathbb{E}\Big[\big(r_{i,t}+\gamma_l&\max\limits_{a_{t+1}}Q_l(s_{t+1},g_t,a_{t+1};\theta_l^{\prime},w_1,w_2) \\
&-Q_l(s_t,g_t,a_t;\theta_l,w_1,w_2)\big)\nabla_{\theta_l}Q_l(s_t,g_t,a_t;\theta_l,w_1,w_2)\Big].
\end{aligned}\end{equation}
Then the RMSProp optimizer in (\ref{eq4.5_9}) is employed to update network parameters.
Please refer to Algorithm \ref{algo1} for detailed procedure of the hierarchical algorithm.
\begin{algorithm}[htb]
    \caption{HRL Based Relay Selection and Power Allocation}
    \label{algo1}
    \begin{algorithmic}[1]
        \STATE Initialize experience replay buffer: $\mathcal{B}_h$ for higher level, and $\mathcal{B}_l$ for lower level.
        \STATE Initialize deep neural network parameters $\theta_l$ for low-level Q-network.
        \STATE Initialize exploration probability $\epsilon=1$ and anneal factor $\sigma$ for controller.

        \FOR{episode $u=1,2,\dots,u_{max}$}
            \STATE Initialize the environment, obtain initial state $s_0$.
            \STATE Choose goal $g_u$ according to policy $\pi^h$.
            \FOR{time slot $t=1,2,\dots,t_{max}$}

                        \STATE Set goal $g_t=g$ for current time slot.
                        \STATE Choose action $a_t$ using epsilon-greedy method with parameter $\epsilon_{l,g_t}$.
                        \STATE Execute action $a_t$, receive internal reward $r_{i,t}$ from the environment and observe next state $s_{t+1}$.

                        \STATE Collect and save the tuple $(s_t,g_t,a_t,r_{i,t},s_{t+1})$ in $B_l$.
                        \STATE Randomly choose a batch of index, and sample transitions from $B_l$.
                        \STATE Perform gradient descent and update low-level Q-network according to (\ref{eq5_7}).
                        \STATE Update current state.

                \STATE Update internal exploration probability $\epsilon\leftarrow\epsilon-\sigma$.
            \ENDFOR
            \STATE Calculate external reward $r_{e}$ according to (\ref{eq5_3}).
            \STATE Collect and save the tuple $(s_0,g_u,r_{e},s_{t_{max}})$ in $B_h$.
            \STATE Read data from $B_h$, update preference values according to (\ref{eqappendix_6}).
            \STATE Update probability distribution according to (\ref{eq5_4}).
        \ENDFOR
    \end{algorithmic}
\end{algorithm}

\section{Evaluation}\label{sect Result}
In this section, we first introduce the setup of simulation environment. We then carry out experiments to evaluate the proposed algorithms.

\subsection{Experiment Setup}\label{subsect_setup}
Similar to \cite{ref1}, in the two-hop cooperative relay network, channel vectors between any two nodes in each time slot are calculated according to formula (\ref{eq4_2}), where the correlation coefficient $\rho$ is set to be 0.95.
The maximum total power for transmission $P_{max}$, which is the sum of power for source and relay, is limited to 4W.
On the other hand, we set the outage threshold as $\lambda=2.0$ during training process, which means an outage will occur when the MI is lower than 2.0.

To implement our proposed framework, learning step size $\zeta$ in high-level meta-controller is set to be 0.1. Elements in preference value vector are initialized to be 0, and the vector will be updated after each inner loop is completed.

In low-level controller, we use two separate dueling deep Q networks, which share the same structure.
Both dueling deep Q networks include two hidden layers, each of which has 50 neurons, and we employ ReLU function for all hidden layers as activation function.
The number of neurons in input layer is equal to the sum of numbers of states and goals, and the number of neurons in output layer corresponds to the dimensions of low-level action.

Note that, deep Q network requires the action space to be discrete, so we have discretized the power allocation in the environment and set $L$ different power-levels for the agent to choose from. For comparison, we use the following methods as baseline in our experiments.

\textbf{Random Selection}:
For each time slot, the agent randomly selects a relay to perform cooperative communication with random transmission power.

\textbf{DRL Based Approach}:
Our DRL method for minimizing outage probability is proposed in Section \ref{sect DRL}.
We employ traditional DQN framework to make relay selection and power allocation at the same time, and it is now used as one of baseline methods.

\subsection{Numerical Results}\label{subsect_results}
Consider that our DRL method and HRL method both use the structure of deep Q network, therefore we first study the influence of different hyper-parameters on the convergence performance, to obtain the optimal network structure.
Note that, average success rate with different hyper-parameter values is tested 10 times, and mean curves and ranges are then recorded.

First of all, the learning rate for updating network parameter should have an appropriate value.
If the learning rate is too small, such as 0.1 (orange solid line), it will lead to local optimum.
If the learning rate is too large, such as 0.0001 (pink dotted line), it will then take much more time to converge.
As shown in Fig. \ref{HyperParameters}(a), we finally set the learning rate as 0.001 for the following simulations.

\begin{figure}[h]
	\centering

    \includegraphics[scale=0.41]{fig2/hyper_parameters.pdf}
	\caption{Average success rate under different parameters: (a) learning rate, (b) memory size, (c) batch size, (d) training interval.}
	\label{HyperParameters}
\end{figure}

Experience replay buffer stores experience tuples obtained by the agent.
In Fig. \ref{HyperParameters}(b), we study the effect of replay buffer size, \textit{i.e.}, memory size, on the performance of convergency.
However, unlike Fig. \ref{HyperParameters}(a), different memory sizes have little influence on the final value that average success rate converges to.
Therefore, we directly select the replay buffer size as 8000.

When training, a batch of data is sampled from the experience buffer to improve the DNN.
In Fig. \ref{HyperParameters}(c), we further fix the memory size, and study the effects of different batch sizes during training on the convergence performance.
It can be found that training with a small batch size cannot take advantage of all data stored in experience buffer, and converges slowly.
While training with a large batch size, such as 256 (pink dotted line), has the fastest convergence speed, although it will consume much more time during training process.

Finally, we investigate the convergence performance under different training intervals, as shown in Fig. \ref{HyperParameters}(d).
Theoretically, the shorter the training interval, the faster the convergence speed.
However, shorter training interval also means more training times, which will result in a certain waste of computing resources.
On the other hand, we find that the final convergence values with time intervals of 5 and 10 are very close.
Considering the above reasons, we finally set training interval as 10.

Set the above hyper-parameters to optimal values and apply them to all deep Q networks, then we carry out the following experiments.

In training process, we evaluate the performance of different methods in AF environment and DF environment, with relay number $K=10$ and power level number $L=10$. The result is depicted in Fig. \ref{TrainReward}.
\begin{figure}[ht]
	\centering

    \includegraphics[scale=0.6]{fig2/iteration_vs_probability_combine.pdf}
	\caption{Average success rate using different protocols.}
	\label{TrainReward}
\end{figure}

It can be observed that when using the method of random selection, the performance is always very poor.
On the other hand, both DRL method and our hierarchical algorithm can be effectively trained, and their average success rate curves eventually converge to a stable value with slight changes.

However, our HRL method can achieve a lower outage probability. Take employing DF protocol as an example, with DRL method, the average success rate is only about 0.82, which means the outage probability of communication system is about $18\%$.
When employing our hierarchical method, the average success rate is closer to 0.9, with an improvement of about $5\%$.
On the other hand, our HRL method has an obvious faster learning speed, which converges after about 10 iterations, while DRL method needs about 40 iterations to reach the convergence value.

When comparing results under with different protocols, Fig. \ref{TrainReward} shows that there is little difference in convergence speed of the same method, but the fluctuation of convergence value using DRL method is larger.
We also observe that the average success rate using AF protocol is generally lower than that using DF protocol, which leads to less successful experiences for agents to learn from.
The performance of DRL method is obviously affected by this factor, but our HRL method can still converge to a value with slighter changes.
What's more, in the case of using AF protocol, the average communication success rate obtained by our HRL method is approximately $15\%$ higher than that obtained by DRL method.
Compared with DF protocol, the performance gap between the two methods under AF protocol is larger.
In all, our HRL agent can learn a better strategy faster for dynamic relay selection and power allocation, in both AF communication environment and DF communication environment.

Then, we evaluate the performance of our proposed hierarchical method and DRL method under different search space scales.
We conduct this experiment in DF communication environment. As shown in Fig. \ref{TestKL}, we study two scenarios where the number of relays $K$ and power levels $L$ are both set to be 10 or 20.
\begin{figure}[ht]
	\centering

    \includegraphics[scale=0.6]{fig2/different_amount.pdf}
    \caption{Outage probability under different relay number $K$ and power level number $L$.}
	\label{TestKL}
\end{figure}

As $K$ and $L$ increased, there is a small increase in the average success rate of our HRL method, which increases from 0.87 to 0.91.
With optional power level $L$ increased, transmission power can be allocated more efficiently at source and relay node, resulting in an improved communication success rate.
In addition, Fig. \ref{TestKL} also vividly shows that the result curve obtained by our HRL method is much smoother than that obtained by DRL method.

One the other hand, DRL method performs worse under larger $K$ and $L$, where the average success rate drops by about $5\%$ after convergence.
In a larger search space, it becomes more difficult for the agent to select appropriate relay and power simultaneously.
As a result, there are fewer successful explorations, which leads to a problem of sparse reward.
Therefore, it takes more training iterations for DRL agent to converge, and the fluctuation becomes larger.
It is worth noting that the number of successful explorations is also small when using the AF protocol (in Fig. \ref{TrainReward}), but this situation is actually different from that shows in Fig. \ref{TestKL}.
In the previous experiment, even if the optimal action is taken, there is still a high probability that the communication will fail due to the uncertainty of the channel.
However, in this experiment, the optimal action with good return exists, but the agent may not be able to find this action policy as the search space is too large.
Therefore, the problem of sparse reward usually refers to the latter situation.

Traditional DRL methods usually perform poorly in the environment with sparse rewards, due to the lack of positive experience to learn from.
However, by making different hierarchies, our method can reduce the complexity of search space, which ensures the efficiency of exploration and learning.
Therefore, when employing our proposed hierarchical method, we can still obtain a more stable behavior policy for relay selection and power allocation.

After 100 iterations of training, we obtain dynamic relay selection and power allocation policies by applying both HRL method and DRL based method.
To further evaluate the robustness of different methods, we evaluate the performance by using these well-trained policies under different outage thresholds, and the result is depicted in Fig. \ref{TestReward}.
\begin{figure}[ht]
	\centering

	\includegraphics[scale=0.75]{fig2/rate_vs_probability.pdf}
    \caption{Outage probability under different outage threshold $\lambda$.}
	\label{TestReward}
\end{figure}

This experiment is conducted in DF environment, and the communication outage threshold ranges from 1.6 to 2.4 in testing process.
The only difference between testing and training is that, the parameters of all networks in testing process are fixed, which means DNN is only used to provide a best action rather than executing further learning.

As we can see from Fig. \ref{TestReward}, both HRL policy and DRL policy trained in a smaller search space can be applied to other situations.
However, DRL policy trained in a larger search space performs poorly when testing, while we can still obtain proper actions in different environments by following our HRL policy.

Take $\lambda=1.8$ as an example, the outage probability under different policies is different.
In terms of the small search space, the outage probability using HRL policy is lower than 0.03, and that using DRL policy and random selection are about 0.07 and 0.8.
In terms of the large search space, outage probability using HRL policy is about 0.04, and that using DRL policy and random selection are 0.8 and nearly 1.0.
It is obvious that our HRL method is more robust and can greatly reduce the outage probability, which means that HRL agent can perform better relay selection and adjust power allocation more reasonably according to the current state after training.

\section{Conclusion and Future Works}\label{sect Conclusion}
In this paper, we propose an HRL method to dynamically select relay and allocate transmission power in a two-hop cooperative communication model, in order to minimize outage probability under a total transmission power constraint.
Unlike traditional studies, our method does not require any assumptions about channel distribution, but relies on the interaction between the agent and the communication environment.
Compared with existing RL-based methods, we propose an outage-based reward function.
Our reward function uses only binary reward that indicates the result of communication, while other RL-based methods require concrete representations of feedback information, such as instantaneous SNR or MI.
We further design an HRL framework by decomposing relay selection and power allocation into two sub-tasks, which reduces the search space.
Simulation results show that our HRL method can reduce outage probability by $5\%$, and reach convergence 30 iterations earlier than DRL methods in both AF and DF communication environment.
In addition, the hierarchical method can effectively solves the problem of sparse reward, while other methods can hardly deal with it.

Our HRL method provides a novel way for the research of resource allocation and optimization in the field of communication. However, the total transmission power is discretized into enumerable power levels in our framework, which can be improved further. In future works, we would like to explore new methods applicable to continuous action space.

\bibliographystyle{IEEEtran}