text:
\section{Introduction}

In recent times, the emergence and widespread use of large-scale data centers with
massive power requirements has elevated the problem of energy-efficient scheduling
to one of paramount importance (see e.g.~\cite{BirmanCR09}). A natural strategy
for achieving energy savings is that of {\em partial shutdown}, i.e. only
a subset of machines/processors are active at any point of time. This immediately leads to
the following scheduling question: {\em which set of machines should be activated
to serve a given set of jobs?} Note that such a schedule must address twin
objectives:
\begin{itemize}
	\item The total cost (e.g. in terms of energy consumption) of all
	machines used in the schedule must be small (thereby achieving energy
	efficiency).
	\item The sum of processing times of all jobs assigned to any machine
	must be small (thereby satisfying throughput requirements).
\end{itemize}
Motivated by this application, Khuller, Li, and Saha~\cite{KhullerLS10}
introduced the {\em machine activation} problem that involves scheduling jobs
to machines so as to minimize the cost of the machines used in the schedule, while
ensuring that the ``load'' on any machine is small.
Observe that treating each of these objectives individually leads to classical
problems in combinatorial optimization, namely {\em minimum set cover} and
{\em minimum makespan scheduling} on parallel machines,
that have been extensively studied over
the last thirty years. The novelty of the algorithm proposed in
\cite{KhullerLS10} for the machine activation problem lies in being able to
handle both objectives simultaneously.

More formally, let $M$ be a set of $m$ machines and $J$ be a set of $n$ jobs,
where the {\em processing time} of job $j$ on machine $i$ is $p_{ij} > 0$.
Further, machine $i$ has {\em startup cost} $c_i$. A {\em schedule}
is defined as an assignment $S: J\rightarrow M$ of jobs to machines; we denote
the set of jobs assigned to machine $i$ in schedule $S$ by $J^{(S)}_i$.
The set of {\em active machines} $M^{(S)}_A$ in schedule $S$ are the machines
to which at least one job has been assigned, i.e.
$M^{(S)}_A = \{i\in M: J^{(S)}_i \not= \emptyset\}$, and the cost of schedule $S$
is defined as $\sum_{i\in A^{(S)}} c_i$. The {\em load} $\ell^{(S)}_i$ on machine $i$
is the sum of processing times of all jobs assigned to machine $i$,
i.e. $\ell^{(S)}_i = \sum_{j\in J^{(S)}_i} p_{ij}$, and the {\em makespan}
$\ell^{(S)}_{\max}$ is the maximum load on a machine, i.e.
$\ell^{(S)}_{\max} = \max_{i\in M} \ell^{(S)}_i$. (Often, we will drop the
superscript $(S)$ in the above notation if the schedule is clear from the context.)
The objective of the machine activation problem is to obtain a schedule of minimum
cost, subject to the constraint that its makespan is at most some given value $\bf L$.

In real-life scheduling tasks, the set of jobs is often not known in advance.
This has motivated extensive algorithmic research in online scheduling problems,
where the  set of machines are available offline but the jobs appear online
and have to be scheduled to a machine when they arrive. A natural and
important question left open in \cite{KhullerLS10} was to obtain an algorithm for
the machine activation problem in the online model. Here, the set of machines,
their individual startup costs, and the budget on the total startup cost of the
machines activated by the schedule are known offline, but the jobs arrive online.
The processing time of a job on each machine is also revealed on arrival of the
job. The goal is to assign the arriving job to a machine such that the cost
of the resulting schedule is minimized subject to the constraint that its makespan
is at most $\bf L$. We call this the {\em online machine activation problem}.

\noindent
\paragraph{Our Contributions.}
Our main contribution is a randomized online algorithm
for the machine activation problem with a bicriteria competitive ratio of
$(O(\log (mn)\log m), O(\log m))$: suppose an offline
optimal schedule for an instance of our problem has cost $B$ and makespan at most $\bf L$;
then our online algorithm produces a schedule of expected
cost $O(B\log (mn)\log m)$ and makespan $O({\bf L}\log m)$
with high probability.
\begin{theorem}
\label{thm:main}
There is a randomized online algorithm for the machine activation problem that
has a bicriteria competitive ratio of $(O(\log (mn)\log m), O(\log m))$.
\end{theorem}

In the minimum set cover problem, we are given a collection of subsets defined on
a universe of elements. The goal is to select a minimum cost sub-collection
such that every element of the universe is in at least
one selected subset. If the elements appear online
and the current selection of subsets must cover every element that has
appeared thus far, a lower bound of $\Omega(\log m\log n)$ is
known~\cite{Korman05}
for $m$ sets and $n$ elements, under standard complexity assumptions.
Since the set cover problem is a special case of the machine activation problem
where the limit $\bf L$ on the makespan of the schedule is $\infty$,
the competitive ratio in the cost of the schedule for any online algorithm
for the machine activation problem must be $\Omega(\log m \log n)$.

On the other hand, the (minimum makespan) scheduling problem
for unrelated parallel machines is
defined as that of distributing $n$ jobs among $m$ machines so as to minimize the
makespan of the schedule (machines do not have cost). It can be shown using standard
techniques that an online algorithm that produces a schedule of makespan at most
$\alpha{\bf L}$ for the online machine activation
problem can be used to obtain an $O(\alpha)$-competitive algorithm for the online
scheduling problem on unrelated parallel machines.
It is well-known~\cite{AzarNR95} that the competitive ratio of any algorithm
for the latter problem is $\Omega(\log m)$; therefore, this
lower bound also holds for the online machine activation problem.
\begin{theorem}
\label{thm:lowerbound}
No algorithm for the online machine activation problem can have a competitive
ratio of $o(\log m)$ in the makespan. Further, under standard complexity
assumptions, no algorithm for the online machine activation problem can have
an approximation factor of $o(\log m \log n)$ in the cost of the schedule.
\end{theorem}
\noindent
\paragraph{Our Techniques.}
Our algorithm draws inspiration from the techniques used to solve the online
versions of the set cover problem and the scheduling problem for unrelated
parallel machines. So, let us first summarize the key ideas involved in these
two algorithms. For the latter problem, Aspnes {\em et al}~\cite{AspnesAFPW97}
gave an elegant solution based on the following exponential potential function:
if the current load on machine $i$ is
$\ell_i$, then its potential is  $a^{\ell_i}$ for some constant $a$.
The algorithm assigns the arriving job to a machine that suffers the minimum
increase of potential, i.e. to machine
$i = \arg\min_{i\in M} (a^{\ell_i+p_{ij}} - a^{\ell_i})$. Observe that if all
processing times are scaled down sufficiently and $a$ is small enough, then
$a^{\ell_i + p_{ij}} - a^{\ell_i} \simeq a^{\ell_i} (a-1) p_{ij}$, i.e. the
increase in potential function is linear in the processing time but exponential
in the current load on the machine. Therefore, the algorithm favors
lightly loaded machines in preference to those offering low processing times.
It was shown in \cite{AspnesAFPW97} that this strong ``bias'' for lightly loaded
machines ensures that the sum of the potentials of all machines is at most
$O(m)$ times that in an optimal offline schedule, thereby leading to a competitive
ratio of $O(\log m)$ on the makespan of the schedule.

For the online set cover problem, Alon {\em et al}~\cite{AlonAABN09} introduced
a two-phase online primal dual framework that works as follows. In the first
phase, the goal is to obtain a feasible fractional solution to the online
instance of the problem. In each step of this phase, in response to a new
constraint (i.e. a new element in the online set cover problem)
that arrives online, the fractional solution is updated to preserve
feasibility while ensuring that the cost incurred can be accounted for by a
suitably updated dual solution.\footnote{The dual was not used explicitly in Alon
{\em et al}'s original analysis. In fact, we will also not use the dual explicitly
even though our algorithm can also be analyzed via the dual.}
In the second phase, the fractional solution is rounded {\em online} to obtain
an integer solution. It is important to note that while the two phases are presented
sequentially for clarity, the algorithm needs to operate both phases (the fractional
updates followed by the rounding) in response
to response to the arrival of a new constraint.
These two phases aim to notionally distinguish between
the information-theoretic aspect of the online problem and the computational aspect
of rounding a fractional solution. In fact, Alon~{\em et al} showed
for $m$ sets and $n$ elements, this two-phase framework can be used to
design an algorithm that has a competitive ratio of $O(\log m\log n)$, where the
$O(\log m)$ factor arises in the first phase due to information-theoretic
limitations of the online algorithm, and the $O(\log n)$ factor arises in the
second phase due to computational limitations encapsulated by the integrality
gap of the linear program (LP) for set cover.

\begin{figure}[!htb]
	\centering
	\small
	Minimize $\quad\sum_{i\in M} c_i x_i\quad$ subject to
	\begin{eqnarray}
		\label{eqn:packing}		\sum_{j\in J} p_{ij} y_{ij} & \leq &  x_i {\bf L} \quad \forall i\in M \\
		\label{eqn:fraction}		y_{ij} & \leq & x_i\quad \forall~i\in M,~j\in J \\
		\label{eqn:covering}		\sum_{i\in M} y_{ij} & \geq & 1\quad \forall~j\in J \\
		\label{eqn:range1}		x_i & \in & \{0, 1\} \quad \forall i\in M \\
		\label{eqn:range2}		y_{ij} & \in & \{0, 1\}\quad \forall i\in M,~j\in J
	\end{eqnarray}
	\caption{\small The integer scheduling linear program (or ISLP). In the fractional scheduling linear
		program (or FSLP), Eqns.~\ref{eqn:range1} and \ref{eqn:range2} are relaxed to $0 \leq x_i \leq 1$
		for all machines $i\in M$, and $0\leq y_{ij}\leq 1$ for all machines $i\in M$ and jobs $j\in J$,
		respectively.}
	\label{fig:lp}
\end{figure}

This two-phase primal dual framework has since been extensively
used for various online problems (see \cite{BuchbinderN09b} for a survey);
our algorithm also uses this framework. However, to the best of our knowledge,
whereas all previous applications of the framework have been to LPs that have
exclusively covering or packing constraints, we apply the framework to a mixed LP.
Consider the integer LP formulation of our problem
given in Fig.~\ref{fig:lp} (we call this the {\em integer scheduling LP} or ISLP).
The variable $x_i$ is 1 iff machine $i$ is active,
and $y_{ij}$ is 1 iff job $j$ is assigned to machine $i$. In the fractional
relaxation (which we call the {\em fractional scheduling LP} or FSLP),
these variables are constrained to be in the range $[0, 1]$
instead of Eqns.~\ref{eqn:range1} and \ref{eqn:range2}.
Note that we have both covering
(Eqn.~\ref{eqn:covering}) and packing (Eqn.~\ref{eqn:packing}) constraints
in the ISLP/FSLP.

We interpret the online set cover algorithm as one that
maintains a bound on a potential function that is linear in the cost of the
active machines, whereas the online scheduling algorithm translates its
objective of minimizing makespan into maintaining a bound on the value of
a potential function that is exponential in the load on each machine. We
design a potential function that combines both objectives: it is linear in
the startup cost and exponential in the load on a machine. Our fractional
algorithm preferentially assigns (fractions of) jobs to machines in a way that
leads to a small increase in this potential. It should be noted that even
if we use this potential function, we cannot afford to simply assign
the entire job  to the machine that would suffer the minimum increase in
potential---such a
greedy strategy can be easily shown to fail even for the special case of the
online set cover problem. Instead, our algorithm creates a dynamic list of
machines in increasing order of its change in potential if the current job
were assigned to it, and then assigns larger fractions of the job to
machines that appear earlier in the order.
This allocation is also somewhat unusual in that
the increase in load on a machine is inverse in the
value of this potential function itself, i.e. inverse exponential in the
current load on the machine. For some of the machines in the
order, this might also involve increasing the fraction to which the machine
is active, i.e. increasing the value of $x_i$, if the assignment violates
Eqn.~\ref{eqn:packing} or \ref{eqn:fraction}. The analysis of the fractional
algorithm involves proving a bound on the value of the cumulative
potential function over all the machines. Depending on the behavior of a
fixed optimal offline solution (which is unknown to the algorithm and is
used only for analysis), we classify jobs into three different categories, and
prove a bound on the total increase in potential due to jobs in each individual
category using three different techniques:
\begin{itemize}
\item The increase of potential for jobs in category 1 is charged globally to
the offline optimal solution using a primal dual argument (without introducing
the dual solution explicitly) similar in spirit to that used by Alon {\em et al}
in~\cite{AlonAABN09} for the online set cover problem.
\item We give a bound on the increase in potential for each individual job in
category 2 by showing that the ratio of increase in potential to the fraction
of job assigned is bounded.
\item We give a global bound on the increase of potential for jobs in category
3 by using a recursive argument which is similar in spirit
to the one used in the online scheduling
algorithm for unrelated machines by Aspnes {\em et al}~\cite{AspnesAFPW97}.
\end{itemize}
The novelty of our analysis lies in being able to seamlessly combine, and
non-trivially extend, the disparate
techniques from \cite{AlonAABN09} and \cite{AspnesAFPW97}.
Ultimately, we show that the potential of the fractional schedule
produced by our algorithm is $O(m\log m)$, where the cost and makespan
of the offline optimal solution are respectively $\Omega(m)$ and 1 by
a suitable initial scaling. We hope that the algorithmic
techniques developed in this paper to simultaneously handle packing and
covering constraints will be useful for solving other online optimization
problems that can be expressed as mixed LPs as well.

In the second phase of the algorithm, we use an online randomized rounding scheme
to obtain an integer solution. To ensure that the expected cost of
the integer schedule is bounded by that of the fractional schedule, each machine $i$
is activated with probability proportional to $x_i$. This is implemented online using
standard techniques. The more challenging aspect of the rounding is the actual
scheduling of jobs to active machines. The natural approach would be to assign
job $j$ to machine $i$ with probability $y_{ij}$.
Translated to conditional probabilities, this implies
that job $j$ should be assigned to machine $i$ with probability $z_{ij} = y_{ij}/x_{ij}$
if machine $i$ is active, where $x_{ij}$ is the value of $x_i$ at the end of the
update to the fractional solution for job $j$. This immediately implies that the expected
load on a machine in the integer schedule is at most that in the fractional schedule.
However, our goal is to obtain a bound
on the {\bf makespan} of the integer solution; in fact, since the events of
jobs being assigned to a fixed machine are positively
correlated, a bound on the expected load does not immediately
yield a concentration bound on the load. Instead, we show that even if job $j$
were to be assigned to machine $i$ {\em unconditionally} with probability $z_{ij}$,
the expected load on machine $i$ given by $\sum_{j\in J} z_{ij}$ would be small.
This overcomes the problem of positive correlation mentioned above since the events are
no longer conditioned on machine $i$ being active. We now derive
concentration bounds on the load on a machine, which translates to a bound on the
makespan of the integer schedule thereby proving Theorem~\ref{thm:main}.

\noindent
\paragraph{Previous Work.}
Many variants of the machine scheduling (or {\em load balancing}) problem have been
extensively studied in the literature. Perhaps the most celebrated result in
scheduling theory is a 2-approximation for the offline minimum makespan scheduling problem
for unrelated machines due to Lenstra, Shmoys and Tardos~\cite{LenstraST90},
which was later simplified by Shmoys and Tardos~\cite{ShmoysT93}.
Rather surprisingly, this algorithm continues to offer the best competitive ratio for
this problem (and even for several natural special cases such as the restricted
assignment problem) even after more than two decades
of research.
In the online setting, Graham~\cite{Graham66, Graham69} showed that the natural greedy
heuristic achieves
a competitive ratio of $2-1/m$ for $m$ identical machines. The competitive ratio of
this problem has been subsequently improved in a series of results
(see e.g.~\cite{BartalFKV95} and subsequent improvements).
For the more general restricted assignment problem
where the processing time of each job $j$ on any machine is either some value $p_j$
or $\infty$, an online algorithm having competitive ratio $O(\log m)$ was designed
by Azar, Naor and Rom~\cite{AzarNR95}. This algorithm was later generalized to
the unrelated machines scenario by Aspnes {\em et al}~\cite{AspnesAFPW97}
with the same competitive ratio.
Various other models and objectives have been considered for the load balancing
problem; for a comprehensive survey, see \cite{Azar96} and \cite{Sgall96}.
In particular, the machine activation problem was introduced by
Khuller~{\em et al} in \cite{KhullerLS10}, where they gave an
$O(2(1+1/\epsilon)(1 + \ln (n/OPT)), 2+\epsilon)$-approximation algorithm for
any $\epsilon > 0$. Recently, this result was extended by Khuller and
Li~\cite{LiK11} to a more general set of cost functions.

\section{The Fractional Algorithm}
\label{sec:fractional}

\begin{figure}[!htb]
	\small
	\centering
	Minimize $\quad\sum_{i\in M} c_i x_i\quad$ subject to
	\begin{eqnarray}
		\label{eqn:packingnew}	\sum_{j\in J} p_{ij} y_{ij} & \leq & 6 x_i {\bf L} \\
		\label{eqn:fractionnew}		y_{ij} & \leq & 2 x_i\quad \forall~i\in M,~j\in J \\
		\label{eqn:coveringnew}		\sum_{i\in M} y_{ij} & \geq & 1\quad \forall~j\in J \\
		\label{eqn:range1new}		0\quad \leq & x_i & \leq \quad 1 \quad \forall i\in M \\
		\label{eqn:range2new}		0\quad \leq & y_{ij} & \leq \quad 1 \quad \forall i\in M,~j\in J
	\end{eqnarray}
	\caption{\small The relaxed fractional scheduling linear program (or RFSLP).
		Eqn~\ref{eqn:packingnew} is enforced only for partially active machines $i$.
		(Note that for inactive machines, Eqns.~\ref{eqn:packingnew} and
		\ref{eqn:fractionnew} are identical.)}
	\label{fig:lp-relaxed}
\end{figure}

In this section, we will describe the online updates to the fractional
solution to maintain feasibility for FSLP on receiving a new job $j$.
This involves updating the values of
$y_{ij}$ (denoting the fraction of job $j$ assigned to machine $i$)
so as to satisfy Eqn.~\ref{eqn:covering}, and corresponding updates
to the values of $x_i$ if Eqns.~\ref{eqn:packing} or \ref{eqn:fraction}
is violated.
In fact, we relax the constraints in FSLP in two ways.
Let machine $i$ be said to be {\em inactive}, {\em partially active} or
{\em fully active} depending on whether $x_i = 0$, $0 < x_i < 1$ or
$x_i = 1$ respectively. First, for technical reasons,
we relax Eqns.~\ref{eqn:packing} and \ref{eqn:fraction}
to Eqns.~\ref{eqn:packingnew} and \ref{eqn:fractionnew} respectively
(see Fig.~\ref{fig:lp-relaxed}). Further, we
enforce Eqn.~\ref{eqn:packing} only if $x_i < 1$, i.e. if machine
$i$ is not fully active. The load on a fully active machine will
be bounded separately in the analysis.
We call this the {\em relaxed fractional scheduling LP} or RFSLP.

Before describing these updates, let us set up some
conventions that we will use throughout the paper.
We divide all processing times by $\bf L$ at the outset;
this allows us to assume that the makespan of the optimal solution is 1.
We also assume that we know the value $\alpha$ of the optimal
offline (integer) solution, i.e. the minimum startup cost of an offline
assignment of jobs to machines that has makespan at most 1.
This is also without loss of generality because we can guess the
value of the optimal solution, doubling our guess whenever the
current algorithmic solution exceeds the cost bounds that we are
going to prove (thereby implying that our current guess is too small).
In the following discussion, it is sufficient to
know the value of $\alpha$ up to a multiplicative
factor of 2, but for simplicity
of presentation, we will assume that we know it exactly.

Our algorithm uses the value of $n$, which is the total number of
jobs that arrive online. If this value is not known offline, each job
estimates $n$ by assuming that it is the last job. We can show that
using such estimates for $n$ incurs a small {\em additive} factor of
$O(\log \log n)$ in the makespan, and an {\em additive} factor of
$O(\log n\log (mn))$ in the cost of the schedule.
However, for simplicity, the rest of the paper assumes that the value
of $n$ is known offline.

We define the {\em virtual cost} of job $j$ on machine $i$ as
\begin{equation*}
\eta_i(j) = \left\{ \begin{array}{ll}
c_i a^{\ell_i-1} p_{ij}, & \text{if machine $i$ is fully active, i.e. } x_i = 1 \\
c_i p_{ij}, & \text{otherwise}
\end{array}\right.
\end{equation*}
where $a$ is a constant that we will fix later.
(Recall that $\ell_i$ represents the load on machine $i$,
i.e. $\ell_i = \sum_{j\in J} p_{ij} y_{ij}$.)
Let $M(j)$ denote
an ordering of machines in non-decreasing order of virtual cost
$\eta_i(j)$ for job $j$.
Let $P(j)$ denote the maximal prefix of $M(j)$ such that
$\sum_{i\in P(j)} x_i < 1$. (Note that $P(j)$ may be empty.)
If $P(j)\not= M(j)$, then $k(j)$
denotes the first machine in $M(j)$ that is not in $P(j)$;
$k(j)$ is undefined if $P(j) = M(j)$.

If $x_i$ is increased to $x_i + \Delta x_i$ for a partially
active machine $i$, then
we say that the {\em effective capacity} created by this increase
for job $j$ is $\min(2 x_i, 6 \Delta x_i/p_{ij})$.
Note that the effective capacity created by an increase in $x_i$
is a feasible increase in the value of $y_{ij}$ independent of the
current load on machine $i$.

\subsection{The Algorithm}
The algorithm has two phases---an offline pre-processing phase, and
an online phase that (fractionally) schedules the arriving jobs.

\paragraph{Pre-processing.} We multiply the startup cost of every machine
by $\alpha/m$, and discard machines with startup cost greater than $m$
(after the scaling) at the outset. Further, for every machine $i$
whose startup cost is at most 1, we increase its
cost to 1, and initialize $x_i$ to 1. For all other
machines with $1 < c_i \leq m$, we initialize $x_i$ to $1/m$.
At the end of the pre-processing phase, we have the following properties:
\begin{itemize}
	\item The cost of an optimal solution is between $m$ and $2m$.
	\item The cost of every machine is between 1 and $m$.
	\item Every machine whose cost is 1 is fully active; all other
	machines have $x_i = 1/m$.
\end{itemize}

\paragraph{Online Algorithm.}
Suppose job $j$ arrives online.
We increase $y_{ij}$s using the following rules repeatedly until
Eqn.~\ref{eqn:coveringnew} is satisfied.
\begin{itemize}
\item {\bf Type A: $k(j)$ is undefined (i.e. $P(j) = M(j)$) or
$x_{k(j)} < 1$ (i.e. machine $k(j)$ is not fully active).}
We increase $x_i$ to
$\min(x_i(1+1/c_i n), 1)$ for each machine
$i\in P(j)$ (and also for $i= k(j)$ if it is defined), and correspondingly
increase $y_{ij}$ by the effective capacity created in each machine.
\item {\bf Type B: $x_{k(j)} = 1$ (i.e. machine $k(j)$ is fully active).}
We increase $x_i$ to
$\min(x_i(1 + 1/c_i n), 1)$ for each machine
$i\in P(j)$, and correspondingly increase $y_{ij}$ by the effective
capacity created in each machine. Further,
we increase $y_{k(j)j}$ by
$6/\eta_{k(j)}(j) n$.
\end{itemize}

\subsection{Analysis}
Our goal is to show the following bounds
on the makespan and cost of the fractional schedule.
\begin{lemma}\label{lma:fractional-main}
The fractional schedule produced by the online algorithm satisfies
$\sum_{j\in J} y_{ij} p_{ij} = O(\log m)$
for each machine $i$, and
$\sum_{i\in M} c_i x_i = O(m\log m)$.
\end{lemma}
\noindent
We introduce a potential function $\phi_i$ for machine $i$ defined as
\begin{equation*}
\phi_i = \left\{ \begin{array}{ll}
c_i a^{\ell_i-1}, & \text{if machine $i$ is fully active, i.e. $x_i = 1$} \\
c_i x_i, & \text{otherwise.}
\end{array}\right.
\end{equation*}
The cumulative potential function $\phi = \sum_{i\in M} \phi_i$.
Our goal will be to show that $\phi = O(m\log m)$.
This will immediately imply Lemma~\ref{lma:fractional-main}.

We prove the bound on $\phi$ in three steps. First, we bound the increase of
$\phi$ in the pre-processing phase; next, we bound the increase of $\phi$ in
each algorithmic step (of either type A or type B); and finally, we bound the
total number of algorithmic steps.

\noindent
{\bf Pre-processing.}
The next lemma bounds the increase of $\phi$ in the pre-processing phase.
\begin{lemma}
\label{lma:initial}
At the end of the pre-processing phase, $\phi \leq m$.
\end{lemma}
\begin{proof}
The startup cost of each machine that is fully active after pre-processing is 1;
on the other, every partially active machine $i$ has $c_i\leq m$ and $x_i = 1/m$
after pre-processing.
\end{proof}

\noindent
{\bf Single Algorithmic Step.}
Now, we bound the increase in $\phi$ due to a single algorithmic step of either type.
\begin{lemma}
\label{lma:singleA}
The increase in potential in a single algorithmic step of type A
is at most $2/n$.
\end{lemma}
\begin{proof}
The total increase in potential in an algorithmic step of type A (due to increase
in the value of $x_i$ for machines $i\in P(j)\cup \{k(j)\}$) is
at most
$\sum_{i\in P(j)\cup \{k(j)\}} c_i (x_i/c_i n) = (\sum_{i\in P(j)\cup \{k(j)\}} x_i)/n < 2/n$.
\end{proof}
\noindent
\begin{lemma}
\label{lma:singleB}
For any constant $1 < a < 13/12$, the increase in potential in a single algorithmic
step of type B is at most $2/n$.
\end{lemma}
\begin{proof}
The total increase in potential for machines $i\in P(j)$ due to an
algorithmic step of type B is at most
$\sum_{i\in P(j)} c_i (x_i/c_i n) = (\sum_{i\in P(j)} x_i)/n < 1/n$.
The other source of increase in potential is the scheduling of
a fraction of job $j$ to machine $k(j)$, due of which
the load on machine
$k(j)$ increases by $6/c_{k(j)}a^{\ell_{k(j)}-1} n$. The resulting
increase of potential $\phi_{k(j)}$ is
\begin{eqnarray*}
c_{k(j)} (a^{\ell_{k(j)} - 1 + 6/c_{k(j)}a^{\ell_{k(j)} - 1} n} - a^{\ell_{k(j)} - 1})
& = & c_{k(j)} a^{\ell_{k(j)} - 1} (a^{6/c_{k(j)}a^{\ell_{k(j) - 1}} n} - 1) \\
\quad = \quad c_{k(j)} a^{\ell_{k(j)} - 1} \left(\left(1 + (a-1)\right)^{6/c_{k(j)} a^{\ell_{k(j)} - 1} n} -  1 \right)
& < & c_{k(j)} a^{\ell_{k(j)} - 1} \cdot \frac{12(a-1)}{c_{k(j)} a^{\ell_{k(j)} - 1} n}
\quad < \quad \frac{1}{n}.
\end{eqnarray*}
The penultimate inequality follows from the property that
$(1 + x)^{1/y} < e^{x/y} < 1 + 2x/y$,
for any $y \geq x > 0$.
\end{proof}

\noindent
{\bf Number of Algorithmic Steps.}
We classify the algorithmic steps according to a fixed
optimal offline (integer) schedule that we call {\sc opt}. Suppose {\sc opt}
assigns job $j$ to machine ${\text{\sc opt}} (j)$, and let $M_{{\text{\sc opt}}}$ denote the machines
that are active in the optimal offline schedule. The three categories are:
\begin{enumerate}
\item ${\text{\sc opt}} (j) \in P(j)$.
\item ${\text{\sc opt}} (j) \notin P(j)$ and ${\text{\sc opt}} (j)$ is partially active.
\item ${\text{\sc opt}} (j)$ is fully active.
\end{enumerate}
The next lemma bounds the total increase in potential due to algorithmic
steps in the first category above.
\begin{lemma}
\label{lma:category1}
The total increase in potential due to all algorithmic steps in the first
category is $O(m \log m)$.
\end{lemma}
\begin{proof}
In any algorithmic step of the first category, the value of $x_{{\text{\sc opt}} (j)}$
either increases to
$x_{{\text{\sc opt}} (j)}\left(1 + \frac{1}{c_{{\text{\sc opt}} (j)} n}\right)$ or to 1.
Since $x_i$ is initialized
to at least $1/m$ for each machine $i$ in the pre-processing phase,
there are at most
$m + \sum_{i\in M_{{\text{\sc opt}}}} c_i n\log m = O(m n\log m)$
algorithmic steps of the first category.
The lemma now follows from Lemmas~\ref{lma:singleA} and \ref{lma:singleB}.
\end{proof}
\noindent
The next lemma bounds the total increase in potential due to algorithmic
steps in the second category.
\begin{lemma}
\label{lma:category2}
The total increase in potential due to all algorithmic steps in the second
category is $O(m)$.
\end{lemma}
\begin{proof}
Consider the first
$2 c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j} n$ algorithmic steps in the second category for any
particular job $j$. We have two cases:
\begin{itemize}
\item {\bf Case~1}: these
algorithmic steps contain at least $c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j} n$ steps of
type B, or
\item {\bf Case~2}: these algorithmic steps contain at least
$c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j} n$ steps of type A.
\end{itemize}
In case~1, each algorithmic step of type B creates an effective capacity
of $6/\eta_{k(j)}(j) n$ in machine $k(j)$. Since in each such algorithmic
step, $\eta_{k(j)}(j) \geq \eta_{{\text{\sc opt}} (j)}(j) = c_{{\text{\sc opt}} (j)}p_{{\text{\sc opt}} (j)j}$,
the total effective capacity created by these algorithmic steps is
at least 1.

In case~2, let $R(j)$ denote the set $P(j)\cup \{k(j)\}$
for the last of these algorithmic steps. (Note that since
${\text{\sc opt}} (j)\notin P(j)$, $k(j)$ is defined.) Further, let $x^{(1)}_i$ and
$x^{(2)}_i$ respectively denote
the value of $x_i$ for machine $i$ before the first
algorithmic step for job $j$, and after the last
algorithmic step.
For each machine $i\in R(j)$, $x_i$ has been increased in
each of at least
$c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j} n$ algorithmic steps of type A. Thus, for each
machine $i\in R(j)$,
\begin{equation}
\label{eqn:diff}
x^{(2)}_i
\geq x^{(1)}_i \left(1 + \frac{1}{c_i n}\right)^{c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j} n}
= x^{(1)}_i \left(\left(1 + \frac{1}{c_i n}\right)^{c_i n}\right)^{\frac{c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j}}{c_i}}
\geq x^{(1)}_i 2^{\frac{c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j}}{c_i}}
\geq x^{(1)}_i 2^{p_{ij}}
\end{equation}
since $c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j} \geq c_i p_{ij}$ for all machines $i\in R(j)$. The
total effective capacity created in these steps is
\begin{equation*}
\frac{6 (x^{(2)}_i - x^{(1)}_i)}{p_{ij}}
\geq \frac{6 x^{(2)}_i (1 - 2^{-p_{ij}})}{p_{ij}}
\geq \frac{6 x^{(2)}_i}{2}
> x^{(2)}_i.
\end{equation*}
The first inequality follows from Eqn.~\ref{eqn:diff} while the
second inequality follows from the observation that for any $z \leq 1$,
we have
$(1 - 2^{-z})/z \geq [1 - (1 - z + z^2/2)]/z = 1 - z/2 \geq 1/2$.
Hence, the effective capacity created by these algorithmic steps is at least
$\sum_{i\in R(j)} x^{(2)}_i \geq 1$. The lemma now follows from
Lemmas~\ref{lma:singleA} and \ref{lma:singleB} coupled with the fact that
$\sum_{j\in J} c_{{\text{\sc opt}} (j)} p_{{\text{\sc opt}} (j)j} \leq 2m$.
\end{proof}
\noindent
Finally, we bound the total increase of potential due to algorithmic
steps in the third category.
\begin{lemma}
\label{lma:category3}
For any $1 < a < 13/12$,
the total increase in potential due to all algorithmic steps in the third
category is at most $2 + (2/3)\sum_{i\in M_A} c_i a^{L_i - 1}$,
where $L_i$ is the final load on machine $i$ in the schedule produced
by the algorithm and $M_A$ is the set of machines that are fully
activated by the algorithm.
\end{lemma}
\begin{proof}
First,
we consider an algorithmic step of type A. In such a step, $k(j)$ must
be defined since ${\text{\sc opt}} (j)\notin~P(j)$; thus,
$\sum_{i\in P(j)\cup \{k(j)\}} x_i \geq 1$. Further, for each machine
$i\in P(j)\cup \{k(j)\}$, we have
$c_i p_{ij} = \eta_i(j) \leq \eta_{{\text{\sc opt}} (j)}(j)$.
Thus, the fraction of job $j$ assigned in this algorithmic step is at least
\begin{equation*}
\sum_{i\in P(j)\cup \{k(j)\}} \min\left(\frac{6 x_i}{c_i p_{ij} n}, 2x_i \right) \geq \min\left(\frac{3}{\eta_{{\text{\sc opt}} (j)}(j) n}, 1\right),
\end{equation*}
since $\sum_{i\in P(j)\cup\{k(j)\}} x_i \geq 1$.
We first consider the situation where the whole of job $j$ was assigned
in this algorithmic step. In this case, the increase in potential due
to job $j$ is
$\sum_{i\in P(j)\cup \{k(j)\}} (c_i x_i)/(c_i n) < 2/n$; cumulatively for
all jobs, such increases in potential add up to at most 2.
Otherwise, the sum of increase in $y_{ij}$ over all machines in
this algorithmic step is at least $3/\eta_{{\text{\sc opt}} (j)j} n$.
Now, consider an algorithmic step of type B.
Then the increase in $y_{k(j)j}$ is $6/\eta_{k(j)j} n \geq 6/\eta_{{\text{\sc opt}} (j)j} n$
since $\eta_{{\text{\sc opt}} (j)j} \geq \eta_{k(j)j}$. In either case,
the total number of algorithmic steps in the third category
for job $j$ is at most $\eta_{{\text{\sc opt}} (j)j} n/3$.

Note that $\eta_{{\text{\sc opt}} (j)j} \leq c_{{\text{\sc opt}} (j)} a^{L_{{\text{\sc opt}} (j)}-1} p_{{\text{\sc opt}} (j)j}$.
Therefore, summing over all jobs, the total increase in potential
due to algorithmic steps in the third category is bounded by
(using Lemmas~\ref{lma:singleA} and \ref{lma:singleB})
\begin{eqnarray*}
2 + \frac{2}{3} \sum_{j\in J} c_{{\text{\sc opt}} (j)} a^{L_{{\text{\sc opt}} (j)}-1} p_{{\text{\sc opt}} (j)j}
& = & 2 + \frac{2}{3} \sum_{i\in M_{{\text{\sc opt}}} \cap M_A} c_i a^{L_i-1} \left(\sum_{j: {\text{\sc opt}} (j) = i} p_{ij} \right) \\
\quad \leq \quad 2 + \frac{2}{3} \sum_{i\in M_{{\text{\sc opt}}} \cap M_A} c_i a^{L_i-1}
& \leq & 2 + \frac{2}{3} \sum_{i\in M_A} c_i a^{L_i-1},
\end{eqnarray*}
where the penultimate inequality follows from the fact that the
makespan of the optimal offline schedule is at most 1.
\end{proof}
\noindent
Finally, we bound the total potential $\phi$; this immediately yields
Lemma~\ref{lma:fractional-main}.
\begin{lemma}
The online fractional algorithm produces a schedule that satisfies
$\phi = O(m\log m)$.
\end{lemma}
\begin{proof}
Lemmas~\ref{lma:initial}, \ref{lma:category1},
\ref{lma:category2} and \ref{lma:category3} imply
$\phi \leq O(m\log m) + (2/3)\phi$,
which proves the lemma.
\end{proof}

\section{The Online Randomized Rounding Procedure}
\label{sec:rounding}

In this section, we give an online randomized rounding scheme for the fractional solution produced by
the algorithm in the previous section.

\subsection{The Algorithm}
For each machine $i$, we select (offline) a number $r_i$ uniformly at random and
independently from $[0, 1]$. In response to a new job $j$ arriving online, the
algorithm updates the schedule in three steps:
\begin{itemize}
\item {\bf Fractional step.} The fractional schedule is updated according to the algorithm
described in the previous section.
\item {\bf Activation step.} Each inactive machine $i$ that satisfies $r_i \leq 5 x_i(j)\ln (mn)$
is activated.
\item {\bf Assignment step.} Let $M_A(j)$ be the set of active machines in the integer solution. Let
$z_{ij} = y_{ij}/2 x_i(j)$ if $x_i(j) < 1/5\ln (mn)$ and $z_{ij} = y_{ij}$ otherwise.
Let $q_{ij}$ be the normalized probability proportional to $z_{ij}$ in a distribution defined on the
set of machines $M_A(j)$. We assign job $j$ to machine $i$ with probability $q_{ij}$.
\end{itemize}

\subsection{Analysis}
The next lemma is an immediate consequence of the fact that machine
$i$ is active in the integer solution with probability
$\min(5x_i\ln (mn), 1)$.
\begin{lemma}
\label{lma:cost-integer}
The total startup cost of all machines activated in the integer
schedule is $O(m \log m\log (mn))$ in expectation.
\end{lemma}
\begin{proof}
The expected startup cost of machine $i$ in the integer schedule
is at most $5c_i x_i\ln (mn)$;
the lemma now follows from Lemma~\ref{lma:fractional-main} and
linearity of expectation.
\end{proof}
\noindent
To bound the makespan of the integer solution, we first bound the probabilities $q_{ij}$.
(This lemma also shows that with high probability, at least one machine is active for even
one job, thereby proving correctness of the algorithm.)
\begin{lemma}
\label{lma:open}
With probability at least $1 - 1/m$, $q_{ij} \leq z_{ij}$ for all machines $i$ and jobs $j$.
\end{lemma}
\begin{proof}
We show that $\sum_{i\in M_A(j)} z_{ij} \geq 1$ with probability at least $1 - 1/mn$
for any job $j$; the lemma then follows using the union bound over all jobs.
We classify machines into $M_1(j)$ and $M_2(j)$ depending on whether or not
$x_i(j) \geq 1/5 \ln (mn)$. Every machine $i\in M_1(j)$ is also in $M_A(j)$. Therefore,
the contribution of machines in $M_1(j)$
to $\sum_{i\in M_A(j)} z_{ij}$ is exactly $\sum_{i\in M_1(j)} y_{ij}$. On
the other hand, the contribution of each machine $i\in M_2(j)$ to this sum
is $y_{ij}/2x_i(j)$ with probability $5 x_i(j)\ln (mn)$, and 0 otherwise. Define
a random variable $Z_{ij} = 1$ with probability $5 y_{ij} \ln (mn)/2$, and 0 otherwise.
Since $y_{ij}/2 x_i(j)\leq 1$,
\begin{equation*}
{\mathbb{P}}\left[\sum_{i\in M_2(j) \cap M_A(j)} z_{ij} < \sum_{i\in M_2(j)} y_{ij}\right]
\leq {\mathbb{P}}\left[\sum_{i\in M_2(j)} Z_{ij} < \sum_{i\in M_2(j)} y_{ij}\right]
< \frac{1}{mn},
\end{equation*}
by Chernoff bounds (cf. e.g.~\cite{MotwaniR97}).
\end{proof}
\noindent
\begin{lemma}
\label{lma:load-whp}
The makespan of the integer schedule is $O(\log m)$ with probability $1 - 2/m$.
\end{lemma}
\begin{proof}
First, we prove that the load on any machine in the integer schedule is
$O(\log m + \sum_{j\in J} y_{ij} p_{ij})$ with probability at least $1 - 1/m^2$
conditioned on the following:
\begin{itemize}
\item Lemma~\ref{lma:open} holds, i.e. $q_{ij}\leq z_{ij}$ for all machines $i$ and
jobs $j$, and
\item Machine $i$ is active from the outset in the integer algorithm.
\end{itemize}
Given that machine $i$ is active from the outset and $q_{ij} \leq z_{ij}$,
the load on machine $i$ due to job $j$ is $p_{ij}$ with probability at most $z_{ij}$.
Now,
\begin{equation*}
	\sum_{j\in J} z_{ij} p_{ij}
	= \sum_{j\in J: x_i(j) < 1/5 \ln (mn)} \frac{y_{ij} p_{ij}}{2x_i(j)} + \sum_{j\in J: x_i(j) \geq 1/5 \ln (mn)} y_{ij} p_{ij}
	\leq \sum_{j\in J: x_i(j) < 1} \frac{y_{ij} p_{ij}}{2x_i(j)} + \sum_{j\in J} y_{ij} p_{ij}.
\end{equation*}
Let job $j'$ immediately precede job $j$ in the online order; if $j$ is the first job,
$x_i(j') = 1/m$. Then,
\begin{equation*}
	\sum_{j\in J: x_i(j) < 1} \frac{y_{ij} p_{ij}}{2x_{ij}}
	\leq 3 \sum_{j\in J: x_i(j) < 1} \frac{x_i(j) - x_i(j')}{x_i(j)}
	\leq 3 \sum_{j\in J: x_i(j) < 1} \int_{w = x_i(j')}^{x_i(j)} \frac{dw}{w}
	\leq 3 \int_{1/m}^1 \frac{dw}{w}
	\leq 3 \ln m.
\end{equation*}
The first inequality follows from the fractional algorithm which assigns a
fraction $y_{ij} \leq 6(x_i(j) - x_i(j'))/p_{ij}$ of job $j$ to machine $i$.
Since $y_{ij} p_{ij} \leq p_{ij} \leq 1$ for all jobs $j$, it follows using
Chernoff bounds that the load on machine $i$ is
$O(\log m + \sum_{j\in J} y_{ij} p_{ij})$ with probability at least $1 - 1/m^2$.

Using the union bound over all machines and Lemma~\ref{lma:open},
we can now claim that the load on machine $i$ is
$O(\log m + \sum_{j\in J} y_{ij} p_{ij})$ (unconditionally) for all machines $i$,
with probability at least $1-2/m$.
The lemma now follows using Lemma~\ref{lma:fractional-main}.
\end{proof}
\noindent
Finally, we note that Lemmas~\ref{lma:cost-integer} and \ref{lma:load-whp} imply
Theorem~\ref{thm:main}.

\newpage

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\IEEEPARstart{M}{assive} machine type communication (mMTC) is a core use case in 5G and later generations wireless communication systems, where a large number of users are expected to be served \cite{ituvision2015,Bockelmann2018,Mikhaylov2019,Dang2020,Saad2020,Chen2020}. These users are expected to intermittently send small amounts of data, primarily in the uplink, i.e., from the users to the base station. As the mMTC users send a limited amount data, the impact of overhead signaling such as random access or scheduling requests is large. Therefore, a grant-free random access scheme is beneficial for these types of users. In a grant-free random access scheme all the active users will transmit their data at the same time without waiting for a grant from the base station. To be able to be distinguished by the base station, the active users also need to transmit a unique identifier.

Grant-free random access with massive MIMO, i.e., base stations with a large number of antennas, has been widely studied in many papers \cite{deCarvalho2017,Chen2018,Liu2018a,Liu2018b,Senel2018,Haghighatshoar2018,Becirovic2019,Chen2019,Shao2019,Ahn2019,Chen2020}. Conventionally in massive MIMO, the users transmit a pilot in each coherence interval that is used by the base station to estimate the channel. Often these pilots are mutually orthogonal such that the channel estimates will be good. When a grant-free random access scheme is applied to massive MIMO, the same pilots are used to detect the active users. If the number of users exceeds the dimension (in samples) of the channel coherence interval, the users cannot be allocated mutually orthogonal pilots, which leads to pilot collisions. There are two fundamentally different approaches for solving this problem. In the first,  pilots are drawn at random and allowed to be non-orthogonal \cite{Liu2018a,Chen2018,Senel2018,Haghighatshoar2018,Chen2019,Shao2019,Ahn2019}. In the other, the pilots are orthogonal, but the data transmission is spread over many coherence intervals \cite{deCarvalho2017,Becirovic2019}, i.e., pilot-hopping sequences are detected instead of pilots. To reduce the complexity of the detection problem, a common approach for grant-free user detection with massive MIMO is to utilize the sparsity of the solution, i.e., that only a small subset of users are active at the same time. Hence, many of the existing solutions use methods developed in the field of compressed sensing \cite{Chen2018,Liu2018a,Liu2018b,Senel2018,Haghighatshoar2018,Becirovic2019,Ahn2019}.

The scenario for detecting pilot-hopping sequences for grant-free random access was introduced in \cite{deCarvalho2017}. However, in \cite{deCarvalho2017} perfect detection of the pilot-hopping sequences was assumed. In \cite{Becirovic2019}, it was suggested that the non-negative least-squares (NNLS) is an appropriate problem formulation for the pilot-hopping sequence detection problem.

The active users transmit one pilot per coherence interval (the time-frequency interval in which the channel can be assumed to be constant), according to a predetermined sequence. The pilot-hopping sequences are unique for every user but more than one user can transmit the same pilot in a coherence interval. The pilots in each coherence interval are accompanied by an uplink data part. An example of this scheme can be found in Fig. \ref{fig:coherence-interval-RA-MaMi}, with three users, two pilots and a pilot-hopping sequence of length four.
\begin{figure}
	\centering
	\begin{tikzpicture}[every node/.style={rounded corners, rectangle,fill=white,inner sep=1pt,scale=0.8},scale=0.5]
	\node at (-1.3,0.5) {User 3};
	\node at (-1.3,1.5) {User 2};
	\node at (-1.3,2.5) {User 1};

	\draw [pattern=north west lines, pattern color=black]
	(0,0) rectangle node{$\bm{\phi}_1$} (1,1);
	\draw [fill=white]
	(1,0) rectangle node{$D_3(1)$}(3,1);

	\draw [pattern=dots, pattern color=black]
	(3,0) rectangle node{$\bm{\phi}_2$}(4,1);
	\draw [fill=white]
	(4,0) rectangle node{$D_3(2)$}(6,1);

	\draw [pattern=north west lines, pattern color=black]
	(6,0) rectangle node{$\bm{\phi}_1$}(7,1);
	\draw [fill=white]
	(7,0) rectangle node{$D_3(3)$}(9,1);

	\draw [pattern=north west lines, pattern color=black]
	(9,0) rectangle node{$\bm{\phi}_1$}(10,1);
	\draw [fill=white]
	(10,0)rectangle node{$D_3(4)$}(12,1);

	\draw [pattern=north west lines, pattern color=black]
	(0,1) rectangle node{$\bm{\phi}_1$}(1,2);
	\draw [fill=white]
	(1,1) rectangle node{$D_2(1)$}(3,2);

	\draw [pattern=north west lines, pattern color=black]
	(3,1) rectangle node{$\bm{\phi}_1$}(4,2);
	\draw [fill=white]
	(4,1) rectangle node{$D_2(2)$}(6,2);

	\draw [pattern=dots, pattern color=black]
	(6,1) rectangle node{$\bm{\phi}_2$}(7,2);
	\draw [fill=white]
	(7,1) rectangle node{$D_2(3)$}(9,2);

	\draw [pattern=dots, pattern color=black]
	(9,1) rectangle node{$\bm{\phi}_2$}(10,2);
	\draw [fill=white]
	(10,1)rectangle node{$D_2(4)$}(12,2);

	\draw [pattern=dots, pattern color=black]
	(0,2) rectangle node{$\bm{\phi}_2$}(1,3);
	\draw [fill=white]
	(1,2) rectangle node{$D_1(1)$}(3,3);

	\draw [pattern=north west lines, pattern color=black]
	(3,2) rectangle node{$\bm{\phi}_1$}(4,3);
	\draw [fill=white]
	(4,2) rectangle node{$D_1(2)$}(6,3);

	\draw [pattern=north west lines, pattern color=black]
	(6,2) rectangle node{$\bm{\phi}_1$}(7,3);
	\draw [fill=white]
	(7,2) rectangle node{$D_1(3)$}(9,3);

	\draw [pattern=dots, pattern color=black]
	(9,2) rectangle node{$\bm{\phi}_2$}(10,3);
	\draw [fill=white]
	(10,2)rectangle node{$D_1(4)$}(12,3);
	\end{tikzpicture}
	\caption{Three users transmitting pilot-hopping sequences of length four using two pilots, $\bm{\phi}_1$ and $\bm{\phi}_2$. Additionally, in coherence interval $t$ the user $k$ also transmits uplink data $D_k(t)$.  }
	\label{fig:coherence-interval-RA-MaMi}
\end{figure}

In this paper, two architectures for solving this problem using the approach in \cite{Becirovic2019} is proposed. To the best of the authors' knowledge, the only other implementations of similar detectors are \cite{Tran2019}, implementing the algorithm from \cite{Senel2018}, and \cite{Henriksson2020}, implementing the algorithm from \cite{Haghighatshoar2018}. However, these are based on the other principle: having non-orthogonal pilot sequences drawn from random sources.

A preliminary version of this work was presented in \cite{Sarband2020}. In addition to a more detailed presentation and more detailed results, an alternative NNLS algorithm is also considered here. The implementation of the alternative algorithm is shown to require smaller chip area and less energy, although it has slightly worse detection performance for hard\footnote{Unfavorable channel hardening and unfavorable propagation such that the linearization in (\ref{eq:asymptotic-limit}) is not a good approximation.} channel conditions.

In the next section, the detection algorithm is reviewed, algorithms for solving NNLS problems discussed, and the general proposed architectures introduced. Then, in Section~\ref{sec:mainresults}, implementation results and decisions are presented for a considered scenario. Finally, some concluding remarks are given in Section~\ref{sec:conclusions}.

\section{Methods and Procedures}
\subsection{Algorithm to Detect Pilot-Hoppping Sequence}
We consider a single-cell massive MIMO system where the base station is equipped with $M$ antennas and serves $K$ mMTC users, each with a single antenna, see Fig. \ref{fig:sm}. We assume that only a subset of these users is active at once. We consider a block fading channel model where the channel of user $k$ in coherence interval $t$ is denoted by $\mathbf{g}_k^t \in \mathbb{C}^M$. The users transmit $T$ coherence interval long pilot-hopping sequences. There are $\tau_\text{p}$ mutually orthogonal pilots. Each pilot consists of  $\tau_\text{p}$ symbols and have unit norm.
At the pilot phase of coherence interval $t$, the base station receives
\begin{equation}
\mathbf{Y}^t = \sum_{k=1}^{K}\sum_{j=1}^{\tau_\text{p}}\alpha_k S_{j,k}^t\sqrt{\tau_\text{p}p_k}\mathbf{g}_k^t\bm{\phi}_j^H + \mathbf{N}^t,
\end{equation}
where
\begin{equation}
\alpha_k = \begin{cases}
1 & \text{if user }k\text{ is active,} \\
0 & \text{otherwise,}
\end{cases}
\end{equation}
\begin{equation}
S_{j,k}^t = \begin{cases}
1, & \text{if user }k\text{ sends pilot }j\text{ at pilot phase }t\text{,} \\
0, & \text{otherwise,}
\end{cases}\label{eq:Sdef}
\end{equation}
$p_k$ is the transmit power of user $k$, $\bm{\phi}_j \in \mathbb{C}^{\tau_\text{p}} $ is the $j$th pilot, and $\mathbf{N}^t \in \mathbb{C}^{M \times \tau_\text{p}}$ is noise with i.i.d. $\mathcal{CN}(0,\sigma^2)$ elements.
The model uses the assumption that the users start their pilot-hopping sequences at the same coherence interval.
\begin{figure}
	\centering
	\begin{tikzpicture}[scale=0.7]
	\newdimen\R
	\R=1.2cm
	\newdimen\H
	\H=2.3cm
	\draw (0:\H) \foreach \x in {60,120,...,360} {  -- (\x:\H) };
	\node[user,fill=green] (k1) at (-0.5*\H,0.1*\H) {};
	\node[user] (k2) at (0*\H,-0.6*\H) {};
	\node[user] (k3) at (-0.4*\H,-0.3*\H) {};
	\node[user,label={\footnotesize user $k$}] (k) at (0.4025*\H,0.5898*\H) {};

	\path (-2\R,-1.5\R) node(K) {$K$ users};
	\draw[->] (K) .. controls +(up:0.7\R) and +(left:0.7\R) .. (k1);
	\draw[->] (K) .. controls +(right:0.7\R) and +(down:0.3\R)..(k2);
	\draw[->] (K)--(k3);
	\node[user,fill=green] (k4) at(0.5*\H,-0.6*\H) {};
	\node[user] (k5) at(0.6*\H,0.1*\H) {};
	\path (2.1\R,-\R) node(A) {active};
	\path (2.3\R,0.7\R) node(I) {inactive};
	\draw[->] (A)--(k4);
	\draw[->] (I)--(k5);
	\DTLforeach*{coordinates}{\x=x, \y=y}{
		\node[user] at (\x*\H,\y*\H) {};
	}
	\DTLforeach*{coordinatesgreen}{\x=x, \y=y}{
		\node[user,fill=green] at (\x*\H,\y*\H) {};
	}
	\pic at (0,0) {bs};
	\node (bs-point) at (-0.2\R,0.42\R) {};
	\node (antenna-point) at (0\R,0.47\R) {};
	\path (-2.3\R,1.3\R) node (M) {$M$ antennas};
	\draw[->] (M) .. controls +(right:1.3\R) and +(left:0.7\R) .. (bs-point);
	\draw[dashed,<->] (antenna-point) .. controls +(up:0.3\R) and +(left:0.3\R) .. node[above,sloped] {\footnotesize$\mathbf{g}_{k}$} (k) ;
	\end{tikzpicture}
	\caption{System model: a cell with a base station using $M$ antennas serving $K$ users, where only a subset is active. The channel between user $k$ and the base station is denoted by $\mathbf{g}_k$.}
	\label{fig:sm}
\end{figure}

In each coherence interval the base station computes an estimate of the received signal energy over each pilot as
\begin{equation}
E_{i,t} = \frac{(\mathbf{Y}^t\bm{\phi}_i)^H(\mathbf{Y}^t\bm{\phi}_i)}{M} - \sigma^2 = \frac{\lVert \mathbf{Y}^t\bm{\phi}_i \rVert^2}{M} - \sigma^2.
\end{equation}
With block Rayleigh fading,
\begin{equation}
\mathbf{g}_k\sim\mathcal{CN}\left(\mathbf{0},\beta_k\mathbf{I}_M\right),
\end{equation}
we have the channel hardening and favorable propagation properties of massive MIMO \cite{redbook}. The channel harding comes from
\begin{equation}
\frac{\lVert\mathbf{g}\rVert^2}{M} \rightarrow \beta_k, \text{ as }M\rightarrow \infty,
\end{equation}
while the favorable propagation comes from
\begin{equation}
\frac{\mathbf{g}_k^H\mathbf{g}_{k'}}{M} \rightarrow 0,\text{ as }M\rightarrow\infty, k\neq{k'}.
\end{equation}
Using these properties, we can state the ``asymptotic energies'' as
\begin{equation}
E_{i,t} = \frac{\lVert \mathbf{Y}^t\bm{\phi}_i \rVert^2}{M} - \sigma^2 \rightarrow \sum_{k=1}^{K}\alpha_k S_{i,k}^t\tau_\text{p}p_k\beta_k \text{ as }M\rightarrow \infty.
\label{eq:asymptotic-limit}
\end{equation}
From this limit, we can see that the asymptotic energies are linear in the user activity parameters. We assume that the users are using statistical channel inversion power control \cite{massivemimobook},
\begin{equation}
p_k = p \frac{\min_{k'}\beta_{k'}}{\beta_k} = p\frac{\beta_{\text{min}}}{\beta_k},
\end{equation}
where $\beta_{\text{min}} = \min_{k'} \beta_{k'}$ is the large-scale fading to the weakest user and $p$ is a system wide power parameter. With statistical channel information, the received signal power at the base station from each user is the same.

We use the asymptotic energies and the fact that the users perform statistical channel inversion to form the column vector
\begin{equation}
\mathbf{b} = \frac{1}{\tau_\text{p}p\beta_{\text{min}}}
\begin{pmatrix}
E_{1,1} \\ \vdots \\ E_{\tau_\text{p},1} \\ \vdots \\ E_{i,t} \\ \vdots \\ E_{1,T}\\ \vdots \\ E_{\tau_\text{p},T}\label{eq:bvector}
\end{pmatrix}
\end{equation}
and the $T\tau_p \times K$-matrix
\begin{equation}
\renewcommand\arraystretch{0.15}
\mathbf{A} = \begin{pmatrix}
S_{1,1}^1 & & \dots & &S_{1,K}^1 \\
\vdots & & \ddots & & \vdots \\
S_{\tau_\text{p},1}^1 & & \dots & &S_{\tau_\text{p},K}^1 \\
\vdots & & \ddots & &\vdots\\
S_{1,1}^T & & \dots & &S_{1,K}^T \\
\vdots & & \ddots & & \vdots \\
S_{\tau_\text{p},1}^T & & \dots & &S_{\tau_\text{p},K}^T \\
\end{pmatrix}.\label{eq:amatrix}
\end{equation}
One can note that $\sum_{j=1}^{\tau_p} S_{j, k}^t = 1,\ \forall k, t$, as user $k$ uses exactly one pilot sequence in coherence time interval $t$.

The $\mathbf{A}$-matrix for the three users in Fig.~\ref{fig:coherence-interval-RA-MaMi} is
\begin{equation}
\newcommand{\timezero}{%
	\left.\vphantom{\begin{matrix}0\\0\end{matrix}}\right\}%
	t = 1%
}
\newcommand{\timeone}{%
	\left.\vphantom{\begin{matrix}0\\0\end{matrix}}\right\}%
	t = 2%
}
\newcommand{\timetwo}{%
	\left.\vphantom{\begin{matrix}0\\0\end{matrix}}\right\}%
	t = 3%
}
\newcommand{\timethree}{%
	\left.\vphantom{\begin{matrix}0\\0\end{matrix}}\right\}%
	t = 4%
}
\mathbf{A} =
\begin{matrix}
\text{User}
 & \\
 \begin{matrix}
 1 & 2 & 3
 \end{matrix} & \\
\begin{pmatrix}\hspace*{0.2em}
\begin{matrix}
0 & 1 & 1\\
1 & 0 & 0
\end{matrix} \\
\hspace*{0.2em} \begin{matrix}
1 & \tikzmarknode{m1}{1} & 0\\
0 & \tikzmarknode{m2}{0} & 1\\
\end{matrix} \\
\hspace*{0.2em} \begin{matrix}
1 & 0 & 1\\
0 & 1 & 0\\
\end{matrix} \\
\hspace*{0.2em} \begin{matrix}
0 & 0 & 1\\
1 & 1 & 0
\end{matrix}
\end{pmatrix}
& \hspace*{-1em}\begin{matrix}
\timezero  \\ \timeone \\ \timetwo \\ \timethree
\end{matrix}
\end{matrix}
\begin{tikzpicture}[overlay,remember picture]
\node[fill=black,rounded corners,opacity=0.1,fit=(m1)(m2)]{};
\end{tikzpicture}\label{eq:aexample}
\end{equation}
The gray area corresponds to $S_{1,2}^2 = 1$ and $S_{2,2}^2 = 0$, i.e., user 2 transmits pilot sequence 1 in coherence interval 2.

With this vector and matrix, we can state the asymptotic energy equation as
\begin{equation}
\mathbf{Ax} \rightarrow \mathbf{b} \text{ as }M \rightarrow \infty,\label{eq:asymptotic-linear}
\end{equation}
where
$$\mathbf{x} = \begin{pmatrix}
\alpha_1 &\dots &\alpha_k &\dots &\alpha_K
\end{pmatrix}^T.$$

We use the fact that not all the users are active at the same time, i.e., that $\mathbf{x}$ is sparse, to solve the pilot-hopping sequence detection problem. In \cite{Becirovic2019}, it was shown that formulating this as an NNLS problem produce sparse solutions that solve the pilot-hopping sequence detection problem. It was also shown that using additional sparsity promoting techniques, such as NNLS combined with $L_1$-minimization, do not improve the results.

\subsection{Algorithms for Solving NNLS}
The non-negative least squares (NNLS) problem can be stated as
\begin{equation}
\begin{array}{ll}
\text{minimize} & \| \mathbf{A}\mathbf{x} - \mathbf{b}\|_2 \\
\text{subject to} & \mathbf{x} \geq \mathbf{0},
\end{array}
\end{equation}
where the design matrix, $\mathbf{A}$ and an observation vector, $\mathbf{b}$, are given. $\mathbf{0}$ denotes a vector with all zeros. The solution vector, $\mathbf{x}$, minimizes the error of $\mathbf{A}\mathbf{x} - \mathbf{b}$ in the least squares sense subject to all elements in $\mathbf{x}$ being non-negative.
Over the years, many different algorithms have been proposed to solve NNLS problems. For an overview of different NNLS algorithms, we refer the reader to \cite{Slawski2013b}.

The probably most well-known approach is the active set method by Lawson and Hanson \cite{Lawson1987}.
However, an active set-based algorithm relies heavily on branching making it unsuitable for a parallel high-speed implementation since parts of the execution will be sequential in nature.
Besides, the execution is data dependent making it non-deterministic.

Instead, we here use two algorithms with a much higher degree of computational parallelism. As they are iterative, there are still limitations, but as later seen, a very high degree of parallelism is still obtained. The execution in each iteration is deterministic simplifying the design of any control. In addition, since these algorithms rely heavily on matrix-vector multiplications, the structure of $\mathbf{A}$, being a binary matrix consisting of only 0/1 entries as illustrated in (\ref{eq:aexample}), can be readily utilized.

\subsubsection{Fast Projected Gradient for Solving NNLS -- Fast}
The fast projected gradient algorithm \cite{Polyak2015}, is described as
\begin{IEEEeqnarray}{rCl}
	\mathbf{x}^{(0)} & = &  \mathbf{0}, \ \mathbf{p}^{(0)} = \mathbf{0} \label{eq:initialization} \\
	\mathbf{x}^{(k+1)} & = & \max\left\{\mathbf{0},\ \mathbf{p}^{(k)}-\frac{1}{L}\mathbf{A}^T\left(\mathbf{A} \mathbf{p}^{(k)} -\mathbf{b}\right)\right\}  \label{eq:xupdate}\\
	\mathbf{p}^{(k+1)} & = & \mathbf{x}^{(k+1)} + s^{(k)}\left(\mathbf{x}^{(k+1)}- \mathbf{x}^{(k)}\right), \label{eq:pupdate}
\end{IEEEeqnarray}
where
\begin{equation}
c^{(k)} = \frac{1+\sqrt{1+4\left(c^{(k-1)}\right)^2}}{2}, \ s^{(k)} = \frac{c^{(k-1)}-1}{c^{(k)}},
\label{eq:sk}
\end{equation}
with $c^{(-1)} =  0$. $L$ is the Lipschitz constant which controls the step size in each iteration. It should be larger or equal to the spectral radius of $\mathbf{A}^T\mathbf{A}$, i.e., the largest absolute eigenvalue. For a given $\mathbf{A}$, this can be computed off-line. Replacing $\mathbf{p}^{(k)}$ with $\mathbf{x}^{(k)}$ in (\ref{eq:xupdate}) and not using (\ref{eq:pupdate}) and  (\ref{eq:sk})  gives the standard projected gradient algorithm for NNLS. The fast part comes from (\ref{eq:pupdate}), where the previous iteration value is used to accelerate the convergence. $s^{(k)}$ is the step length of the fast part and $\mathbf{p}^{(k)}$ is now denoted the predictor \cite{Polyak2015}.

Note here that the core of the algorithm, the matrix-vector multiplications, will be simple additions since the elements of $\mathbf{A}$ are 0 or 1 as seen from (\ref{eq:aexample}). As such, the architectures can not be directly used for NNLS problems without this property.

\subsubsection{Multiplicative Updates for Solving NNLS -- Mult}
In the multiplicative updates algorithm, the updates in each iteration are instead performed using multiplications and divisions and is described as \cite{Daube1986,Lee2001,Sha2007}:
\begin{IEEEeqnarray}{rCl}
\mathbf{x}^{(0)} & = & \mathbf{1} \\
\mathbf{d} & = & \mathbf{A}^T\mathbf{b}\label{eq:mul2} \\
\mathbf{e}^{(k+1)} & = & \mathbf{A}^T\mathbf{A}\mathbf{x}^{(k)}\label{eq:mul3} \\
\mathbf{x}^{(k+1)} & = & \mathbf{d}\odot \mathbf{x}^{(k)} \oslash {\mathbf{e}^{(k+1)}}, \label{eq:mul5}
\end{IEEEeqnarray}
where $\mathbf{1}$ denotes a vector with all ones and $\odot$ and $\oslash$ represent component-wise (Hadamard) vector multiplication and division, respectively.

Similar to the Fast algorithm, the matrix-vector multiplications consist of additions.

\subsection{Proposed Architectures\label{sec:proposedarchitecture}}
The proposed architectures are iso-morphic mappings of a single iteration of the respective algorithms. Hence, each clock cycle one complete iteration of the algorithms is performed.
The primary reason to do this is to obtain a high degree of parallelism and to easily benefit from the fact that $\mathbf{A}$ only consists of 0 and 1 entries, leading to additions only.
Since the matrices $\mathbf{A}$ and $\mathbf{A}^T$ are constant in all iterations, these constant matrix-vector multiplications can be converted to summation structures where sub-expression sharing can be used to
reduce the area and power consumption of the implementation.

\subsubsection{Fast }
The resulting architecture is shown in Fig.~\ref{fig:architecture}.
It should be noted that most signals are vectors.
There are registers to store $\mathbf{b}$, which is given at the start of a detection process, and $\mathbf{p}^{(k)}$ and $\mathbf{x}^{(k)}$, which are initialized to zero at the start as seen in (\ref{eq:initialization}).
\begin{figure}
	\centering
	\includegraphics[scale=1.2]{architecture.pdf}
	\caption{Proposed architecture of Fast algorithm with number of parallel signals annotated.}
	\label{fig:architecture}
\end{figure}

To implement multiplication by $s^{(k)}$, the value of it should be calculated in each iteration and then multiplied. As visible from (\ref{eq:sk}), the calculation of $s^{(k)}$ is costly in hardware, but it is possible to compute the $s^{(k)}$ values off-line and save them in a memory. Furthermore, $s^{(k)}$ converges to one after some time, so for $k$ and all later iterations, the values will saturate to the largest representable number less than or equal to one.

\subsubsection{Mult}
To avoid the relatively costly multiplications and divisions of the Mult algorithm, a logarithmic number system (LNS) \cite{Swartzlander1975} is used for those parts. This reduces the multiplications and divisions to additions and subtractions, respectively, at the expense of converting the numbers between the linear and logarithmic domains.

The multiplicative updates algorithm, partly in the logarithmic domain, is described as:
\begin{IEEEeqnarray}{rCl}
\hat{\mathbf{x}}^{(0)} & = & \mathbf{0} \\
\hat{\mathbf{d}} & = & \log \left(\mathbf{A}^T\mathbf{b}\right) \\
\hat{\mathbf{e}}^{(k+1)} & = & \log \left(\mathbf{A}^T\mathbf{A}\mathbf{x}^{(k)} \right) \\
\hat{\mathbf{x}}^{(k+1)} & = & \hat{\mathbf{d}} + \hat{\mathbf{x}}^{(k)} - {\hat{\mathbf{e}}^{(k+1)}} \\
\mathbf{x}^{(k+1)} & = & \exp \left(\hat{\mathbf{x}}^{(k+1)}\right).
\end{IEEEeqnarray}

The resulting architecture is shown in Fig.~\ref{fig:MULarchitecture}. Note that $\hat{\mathbf{d}}$ is computed in the first clock cycle and that the first iteration is performed in the second clock cycle. The multiplexers are used for storing the computed value of $\hat{\mathbf{d}}$ and to initialize $\hat{\mathbf{x}}^{(0)}$ in the first clock cycle. Hence, one more clock cycle is required here compared to the Fast algorithm for the same number of iterations.
\begin{figure}
	\centering
	\includegraphics[scale=1.2]{multarchitecture.pdf}
	\caption{Proposed architecture of Mult algorithm partially implemented in the logarithmic domain with number of parallel signals annotated.}
	\label{fig:MULarchitecture}
\end{figure}

To compute the logarithm and exponentiation, either a table or using Mitchell's logarithm approximation \cite{Michell1962} is considered. There are many methods suggested to either perform approximations of the logarithm and exponent functions using smaller tables or to improve the accuracy of Mitchell's algorithm. However, as seen in the results section, the accuracy of the Mitchell approximation is enough.

The idea of Mitchell's logarithm approximation is to utilize the approximation
\begin{equation}\label{eq:logapprox}
\log (x) \approx x - 1
\end{equation}
in the range $1 \leq x \leq 2$.

To benefit from this, the position of the leading one, i.e., the most significant bit being one, of the input $x$ is determined\footnote{Here, the input is always non-negative, so the negative case where the leading zero should be detected does not have to be considered.} and then shifted to the range $1 \leq x < 2$. The position of the first one also determines the integer part of the logarithm. Then, the leading one is removed and the fractional part of the logarithm is obtained as per (\ref{eq:logapprox}).

For computing the exponent, the opposite procedure is performed. A one is put in front of the fractional part and the final result is obtained by shifting the fractional result with the leading one according to the number of integer bits.

\subsubsection{A Note on Pipelining}
It should be noted here that pipelining will not provide any advantage as we need the result of one iteration before being able to start the next one.
If required, one may interleave several detections by introducing pipelining, although that is not considered here since it will increase the detection latency and storage requirements.
Unfolding may provide a small increase of detection rate as unfolding $N$ times will most likely not decrease the clock frequency $N$ times since the synthesis tool probably can optimize some parts better.
Naturally, fully unfolding the architecture the required number of iterations will enable pipelining and streaming processing.
However, the obtained detection rate is most likely more than enough for most practical current and envisioned scenarios.

\section{Results}\label{sec:mainresults}
For the implementation example, we consider a system with parameters as shown in Table~\ref{tab:param}. These parameters are selected based on the number of users and user activity, where the number of antennas, orthogonal pilots, and transmitted pilots are found by simulation to obtain a reasonable detection performance. For more results on how these parameters relate, we refer to \cite{Becirovic2019}. Synthesis results are for 28~nm FD-SOI standard cells with a power supply voltage of 1.0~V, using Design Compiler. Power estimates for the final designs are obtained by simulating the synthesized netlist using SDF data to produce switching activities that are back annotated into the synthesis tool. Compared to the preliminary version \cite{Sarband2020}, here the low leakage library with a slow-slow characterization at 125$^\circ$C is used.

\begin{table}
	\centering
	\caption{Parameters for the Implemented Example.\label{tab:param}}
	\begin{tabular}{l|r}
		Parameter & Value \\
		\hline
		Number of users, $K$ & 1024 \\
		User activity & $2^{-6}$\\
		Number of antennas, $M$ & 96 \\
		Number of orthogonal pilots, $\tau_\text{p}$ & 16 \\
		Number of transmitted pilots, $T$ & 8 \\
	\end{tabular}
\end{table}

\subsection{Number of Iterations}
Determining the number of iterations required for robust activity detection is important since this relates the clock frequency of the implementation with the number of detections per second. Here, we will consider floating-point implementation of the algorithms, while fixed-point implementation is further discussed in Section~\ref{sec:wordlength}.

First, the behavior of the algorithms is illustrated. When iterating the algorithms, the values in $\hat{\mathbf{x}}$ corresponding to active users will be close to one, while the values corresponding to inactive users will be close to zero. A typical sequence of values for a successful decoding is shown in  Fig.~\ref{fig:compconverge}, where the values of  $\hat{\mathbf{x}}$ are shown after each iteration. Here, one can see that the Mult algorithm initially has a faster convergence for active users, but that around iteration 50, the Fast algorithm has caught up and reaches the final values slightly faster. One can see that with about 100 iterations, it is possible to clearly separate the active and inactive users.
\begin{figure}
	\centering
	\includegraphics[scale=.3]{compconvergenceannotated.pdf}
	\caption{Convergence of $\hat{\mathbf{x}}$ using  Fast (dashed) and Mult (solid) algorithms successfully performing an activity detection. Large values corresponds to active users and small values to inactive users. There is a clear separation between the active and inactive users.}
	\label{fig:compconverge}
\end{figure}

Typically, one does not use a fixed-threshold value, but studies the effect of the number of active users that are not detected (missed detections) and number of inactive users that are detected as active (false alarms) with a varying threshold.
The probability of missed detections, $p_{\text{m}}$, and false alarm, $p_{\text{fa}}$, are calculated as
\begin{IEEEeqnarray}{rCl}
	p_{\text{m}}& =\frac{\#\text{undetected active users}}{\#\text{active users}}    \label{eq:pm}\\
	p_{\text{fa}}& =\frac {\#\text{detected inactive users}} {\#\text{inactive users}}\label{eq:pfa}
\end{IEEEeqnarray}
and plotted against each other in a receiver operating characteristic (ROC) curve. The ROC curves for the two considered algorithms are shown in Fig.~\ref{fig:ROCconvergedfloating} when using 4000 iterations, clearly more than in Fig.~\ref{fig:compconverge}, and a large number of instances. Here, it can be seen that the Fast algorithm has a slightly better behavior compared to the Mult algorithm. Each point on the curve corresponds to a threshold and the position, the number of missed detections and false alarms if that threshold was used. Hence, the values in the top left corner corresponds to large thresholds and the values in the bottom right corner corresponds to small thresholds. For example, if a threshold of, say, 0.9 is used, there will almost never be any false alarms, but it is quite common that there are missed detections.
\begin{figure}
	\centering
	\includegraphics[scale=.3]{convergedROC.pdf}
	\caption{Missed detection, $p_{\text{m}}$, versus false alarm rates, $p_{\text{fa}}$, for 4000 iterations results of floating-point realizations of Fast and Mult algorithms.}
	\label{fig:ROCconvergedfloating}
\end{figure}

From Fig.~\ref{fig:compconverge}, one can realize that for good convergence cases, it is possible to get a good separation of the active and inactive users is a few hundred iterations.

Instead, consider a case with bad convergence, as shown in Fig.~\ref{fig:badconverge}. Here, the channel conditions and the selection of active users interacts in such a way that it is very hard to distinguish the active and inactive users. In these cases, one can expect that the curves of some active users are actually below the curves of some inactive users and that there are no distinct value that will separate the two groups.
\begin{figure}
	\centering
	\subfigure[]{\label{fig:badfast}\includegraphics[scale=.28]{Fastbadconvergence.pdf}}
	\hspace{-2mm}
	\subfigure[]{\label{fig:badmul}\includegraphics[scale=.28]{Multbadconvergence.pdf}}
	\caption{Convergence of $\hat{\textbf{x}}$ using  (a) Fast and (b) Mult algorithms where the active and inactive users are not clearly separated.

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

To study this further, the same case as in Fig.~\ref{fig:badconverge}, was solved using 10000 iterations, with the results shown in Fig.~\ref{fig:badconverge2}. Here, it can be seen that the Fast algorithm in Fig.~\ref{fig:badfast2} reaches a more or less stable value in about 1000 iterations, while for the Mult algorithm in Fig.~\ref{fig:badmul2} some variables have not reached a stable value even after 10000 iterations.
\begin{figure}
	\centering
	\subfigure[]{\label{fig:badfast2}\includegraphics[scale=.28]{Fastbadconvergence2.pdf}}\hspace{-2mm}
	\subfigure[]{\label{fig:badmul2}\includegraphics[scale=.28]{Multbadconvergence2.pdf}}
	\caption{Convergence of $\hat{\textbf{x}}$ using  (a) Fast and (b) Mult algorithms where the active and inactive users are not clearly separated. Extended version of Fig.~\ref{fig:badconverge} covering 10000 iterations.

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

Therefore, it can be argued that although Mult gives worse results in Fig.~\ref{fig:ROCconvergedfloating}, the degradation comes primarily from cases where the results will contain detection errors anyway. Using even more iterations will make the curves in Fig.~\ref{fig:ROCconvergedfloating} identical. Hence, we will use the Fast algorithm curve in Fig.~\ref{fig:ROCconvergedfloating} as a reference curve for a converged solution.

In Fig.~\ref{fig:ROCiteretedfloating}, the missed detection, $p_{\text{m}}$, and false alarm rates, $p_{\text{fa}}$, are shown for the floating-point implementations with 128, 192, and 256 iterations,  for the Fast and Mult algorithms. As seen, the results are better for the Fast algorithm, but we know from previous discussions that the Mult algorithm converges faster in the good cases, so this will also be further considered.
 \begin{figure}
	\centering
	\subfigure[]{\label{fig:ROCfastfloatIT}\includegraphics[scale=0.3]{floatingFast.pdf}}
	\hspace{-1mm}
	\subfigure[]{\label{fig:ROCMULfloatIT}\includegraphics[scale=.3]{floatingMult.pdf}}
	\caption{Missed detection, $p_{\text{m}}$, versus false alarm rates, $p_{\text{fa}}$, for floating-point implementations, 128 iterations, 192 iterations and 256 iterations for (a) Fast algorithm, and (b) Mult algorithm.}
	\label{fig:ROCiteretedfloating}
\end{figure}

\subsection{Word Length Optimization}\label{sec:wordlength}
A careful word length optimization has been performed to determine the minimum resolution of the signals while still being able to operate correctly.
In addition, the maximum values have been determined by a combination of interval arithmetic and simulations.
One example of the latter is the output of the $\mathbf{A}$ block that theoretically may require seven additional MSBs compared to the input as up to 87 $\mathbf{p}^{(k)}$\mbox{-}values are summed for the selected $\mathbf{A}$ since the maximum number of ones in a row of the selected $\mathbf{A}$ is 87.
However, simulation show that at most five additional MSBs are required in practice as not all the $\mathbf{p}^{(k)}$\mbox{-}values will be close to one.
Similarly, for the Fast algorithm, the subtraction of $\mathbf{x}^{(k)}$ from $\mathbf{x}^{(k+1)}$ results in a shorter output word length than the input word length. It can be seen from Figs.~\ref{fig:compconverge} and \ref{fig:badconverge} that the values are correlated and do not change much between two iterations.

In the following, $S$ and $U$ denote signed and unsigned representation, respectively.
The numbers following $S$ and $U$ are the number of integer and fractional bits, respectively.
A negative number of integer bits denotes that the corresponding number of fractional bits are not required. As examples, $S(2,2)$ denotes a four-bit signed number with values between $ -\frac{2^{4-1}}{2^2} = -2$ and $\frac{2^{4-1}-1}{2^2} = 1.75$ and $U(-2, 4)$ denotes a two-bit unsigned number with values between $0$ and $\frac{2^2-1}{2^4} = 0.1875$. In general, a $S(I,F)$ number goes from $-\frac{2^{I+F-1}}{2^F} = -{2^{I-1}}$ to $\frac{2^{I+F-1}-1}{2^F} = 2^{I-1} - {2^{-F}}$ and a $U(I,F)$ number goes from $0$ to $\frac{2^{I+F}-1}{2^F} = 2^{I} - {2^{-F}}$.

For the nodes where quantization is performed, i.e., the number of fractional bits are reduced, rounding is used.

The resulting ROC curves for the resulting word lengths, further discussed below, are shown in Fig.~\ref{fig:ROCfloatfixed}. Here, one can see that 256 and 192 iterations are required for the Fast and Mult algorithms to reach close to the floating-point converged performance, respectively. On the other hand, 128 iterations are enough for the Fast algorithm to reach a performance similar to the Mult algorithm.
\begin{figure}
	\centering
	\subfigure[]{\label{fig:ROCfast}\includegraphics[scale=0.3]{fixedFast.pdf}}
	\vfill
	\subfigure[]{\label{fig:ROCMULtable}\includegraphics[scale=0.3]{fixedtableROC.pdf}}
	\hspace{-1mm}
	\subfigure[]{\label{fig:ROCMULMitchell}\includegraphics[scale=.3]{fixedmitchellROC.pdf}}

	\caption{Missed detection, $p_{\text{m}}$, versus false alarm $p_{\text{fa}}$, rates for converged floating-point and fixed-point with 128, 192, and 256 iterations, for (a) Fast algorithm, (b) Mult algorithm using lookup tables for logarithm and exponent computations, and (c) Mult algorithm using Mitchell's approximated logarithm.}
	\label{fig:ROCfloatfixed}
\end{figure}

\subsubsection{Fast}
The value of $L$ for the considered $\mathbf{A}$ is 520.74.
Strictly, this should be rounded up to guarantee convergence under all possible inputs, but simulations show that selecting $L = 512$ works well, resulting in a shift instead of a multiplication.
The values of $s^{(k)}$ are pre-computed and stored in a lookup table (all values for $k \geq 58$ are $1-2^{-5}$).
The values for $\mathbf{p}^{k}$ closely follow the values of $\mathbf{x}^{k}$.
As seen from Fig.~\ref{fig:badfast}, sometimes the values for $\mathbf{x}^{k}$, and therefore $\mathbf{p}^{k}$, can take on values above one.
Hence, we saturate the values to at most one, using a single guard-bit. This operation is combined with the $\max$ operation in (\ref{eq:xupdate}) to the $\min(\max(.))$-operation in Fig.~\ref{fig:wl}.
It can also be noted that (\ref{eq:pupdate}) on rare occasions result in a small negative value, which is confirmed by simulations.
Hence, $\min(\max(.))$-operations are performed before storing $\mathbf{p}^{k}$ as well to keep the word length limited and enabling storage of unsigned values only.

The resulting architecture with annotated word lengths is shown in Fig.~\ref{fig:wl}.
Rounding is here also performed for the input to the matrix-vector multiplication by $\mathbf{A}$, which significantly reduces the complexity of the computation by halving the input word length.
\begin{figure}
	\centering
	\includegraphics[scale=1.2]{architecturewl.pdf}
	\caption{Modified proposed architecture for the Fast algorithm with annotated word lengths for the implemented example.}
	\label{fig:wl}
\end{figure}

The implemented Fast algorithm can therefore be described as
\begin{IEEEeqnarray}{rCl}
	\mathbf{x}^{(0)} & = &  \mathbf{0}, \ \mathbf{p}^{(0)} = \mathbf{0} \nonumber \\
	\mathbf{x}^{(k+1)} & = & \min\left\{\mathbf{m},  \max\left\{\mathbf{0},\ \mathbf{p}^{(k)}-\frac{1}{L}\mathbf{A}^T\left(\mathbf{A} \mathbf{p}^{(k)} -\mathbf{b}\right)\right\} \right\} \nonumber\\
	\mathbf{p}^{(k+1)} & = &  \min\left\{\mathbf{m},\right.  \nonumber \\
& & 	\left.\max\left\{\mathbf{0},\ \mathbf{x}^{(k+1)} + s^{(k)}\left(\mathbf{x}^{(k+1)}- \mathbf{x}^{(k)}\right)\right\} \right\} \nonumber
\end{IEEEeqnarray}
where $\mathbf{m}$ is a vector of the largest positive value representable, $1-2^{-12}$.
\subsubsection{Mult}
As discussed earlier, the Mult architecture is partially based on computation in the logarithmic domain. The results in Figs.~\ref{fig:ROCMULtable} and \ref{fig:ROCMULMitchell} show that the Mitchell logarithm approximation gives the same activity detection performance as using a lookup table.

As part of the word length optimization, the initial value of $\mathbf{x}^{(0)}$ was reduced from 1 to $2^{-7}$. This reduces the values of $\mathbf{e}^{(1)}$, allowing a reduced word length, without reducing the detection performance. In addition, saturation was introduced in order to handle potential overflows of $\hat{\mathbf{x}}^{(k+1)}$.

This leads to that the implemented Mult algorithm, partially in the logarithmic domain, is described as
\begin{IEEEeqnarray}{rCl}
	\hat{\mathbf{x}}^{(0)} & = & \mathbf{s} \\
	\hat{\mathbf{d}} & = & \log_2 \left(\mathbf{A}^T\mathbf{b}\right) \\
	\hat{\mathbf{e}}^{(k+1)} & = & \log_2 \left(\mathbf{A}^T\mathbf{A}\mathbf{x}^{(k)} \right) \\

	\hat{\mathbf{x}}^{(k+1)} & = & \text{sat}\left(\hat{\mathbf{d}} + \hat{\mathbf{x}}^{(k)}_i - {\hat{\mathbf{e}}^{(k+1)}} \right) \\
	\mathbf{x}^{(k+1)} & = & 2^{\hat{\mathbf{x}}^{(k+1)}}
\end{IEEEeqnarray}
where $\mathbf{s}$ is a vector with the initial value, $\log_2\left(2^{-7}\right) = -7$, and $\mathrm{sat}$ corresponds to saturating the result in case of overflow. The resulting architecture with annotated word lengths is shown in Fig.~\ref{fig:wlmul}.
\begin{figure}
	\centering
	\includegraphics[scale=1.2]{multarchitecturewl.pdf}
	\caption{Modified proposed architecture for the Mult algorithm partially implemented in the logarithmic domain with annotated word lengths for the implemented example.}
	\label{fig:wlmul}
\end{figure}

The structure of the used logarithm approximation with base 2 is shown in Fig.~\ref{fig:mitchelllog2}. It consists of a leading one detector (LOD) that outputs the position of the first one in the input word, with the MSB corresponding to 0. As there in general are $I$ integer bits, the value must be offset by $I-1$ to determine the correct integer part of the logarithmic value. The LOD value is also used to left-shift the input such that the leading one comes in the MSB position. This leading one is removed and the remaining part is used as the fractional part of the logarithm and combined with the integer part. Finally, as a zero input value will output the most negative output value, a simple saturation multiplexer is added. The programmable shift is realized using a 11-to-1-multiplexer, which provided the best synthesis results.
\begin{figure}
	\centering
	\includegraphics[scale=.7]{log2arch.pdf}
	\caption{Realization of the Mitchell logarithm approximation with a $U(I,F)$ input format and zero input saturating to the most negative output value.}
	\label{fig:mitchelllog2}
\end{figure}

The logarithm approximation in Fig.~\ref{fig:mitchelllog2} and the corresponding exponentiation require a significantly smaller area than using lookup tables and allow significantly higher clock frequencies. As the detection results are similar for both approaches, as seen from Figs.~\ref{fig:ROCMULtable} and \ref{fig:ROCMULMitchell}, only the Mitchell approximation version is considered in the following.

\subsection{Sub-Expression Sharing}
As mentioned in Section~\ref{sec:proposedarchitecture}, the number of adders for performing the matrix-vector multiplications by $\mathbf{A}$ and $\mathbf{A}^T$ can be reduced.
Here, we use an iterative two-term sub-expression sharing approach with the constraint that each addition should be performed at the minimum depth to keep the power consumption low.
The results of this are shown in Table~\ref{tab:ses}.
As seen, a significant reduction is obtained compared to implementing the additions straightforwardly.
For comparison, the other computations for the Fast approach are, as can be seen from Fig.~\ref{fig:architecture}, $3K + \tau_{\text{p}}T = 3200$ additions/subtractions and $2K$ multiplications, of which $K$ are reduced to a simple shift, so $K=1024$ multiplications. For the Mult approach, $2K = 2048$ additions/subtractions, $K=1024$ logarithm approximations, and $K=1024$ exponent approximations are used, as evident from Fig.~\ref{fig:MULarchitecture}.
\begin{table}
	\centering
	\caption{Number of Adders for the Matrix-Vector Multiplications without and with Sub-Expression Sharing in the Implemented Example.\label{tab:ses}}
	\begin{tabular}{c|cc|c}
		& \multicolumn{2}{c|}{Adders} & \\
		Matrix & Without & With & Reduction\\
		\hline
		$\mathbf{A}$ & 8064 & 5699 & 29.3\%\\
		$\mathbf{A}^T$ & 7168 & 4709 & 34.3\%\\
		\hline
		Total & 15232 & 10408 & 31.7\%
	\end{tabular}
\end{table}

To further clarify the consequences of sub-expression sharing, the $\mathbf{A}$ and $\mathbf{A}^T$ matrix-vector multiplications have been synthesized separately for using the word lengths of the Fast implementation, i.e., the word lengths in Fig.~\ref{fig:wl}. The area and power results are shown in Figs.~\ref{fig:apa}~and~\ref{fig:apat} for $\mathbf{A}$ and $\mathbf{A}^T$, respectively. It can be seen that although the reduction in area is not the $\approx 30\%$ expected, more for $\mathbf{A}$ and less for $\mathbf{A}^T$, sub-expression sharing provides a significant area and power reduction.

\begin{figure}
	\centering
\subfigure[]{\label{fig:aarea}\includegraphics[scale=0.3]{areaAses.pdf}}
	\hspace{-3mm}
\subfigure[]{\label{fig:apower}\includegraphics[scale=.3]{powerAses.pdf}}  \caption{ (a) Area  and (b) power results at different clock frequencies for the matrix-vector multiplication with $\mathbf{A}$ with and without using sub-expression sharing.}
\label{fig:apa}
\end{figure}

\begin{figure}
	\centering
	\subfigure[]{\label{fig:atarea}\includegraphics[scale=0.3]{areaATses.pdf}}
		\hspace{-3mm}
	\subfigure[]{\label{fig:atpower}\includegraphics[scale=.3]{powerATses.pdf}}  \caption{ (a) Area  and (b) power consumption at different clock frequencies for the matrix-vector multiplication with $\mathbf{A}^T$ with and without using sub-expression sharing.}
	\label{fig:apat}
\end{figure}

\subsubsection{Selecting $\mathbf{A}$}

The number of ones in $\mathbf{A}$ is always $KT$, as each user uses $T$ pilots. This leads to that the number of adders without sharing is the same independent of the pilot hopping sequences: $KT-\tau_{\text{p}}T$ for $\mathbf{A}$ and $KT-K$ for $\mathbf{A}^T$.
However, by selecting different pilot hopping sequences, different sub-expression sharing results can be obtained.
Therefore, one can try to optimize $\mathbf{A}$ to give fewer adders, while still having good detection performance.
One may note that qualitatively, detection performance increases the more different the pilot hopping sequences are, while the sharing increases with similarities in pilot hopping sequences.

Similarly, one can see that decreasing the length of the pilot hopping sequences $T$ will lead to fewer adders, assuming that the number of users is constant.
In \cite{Becirovic2019}, it was shown that detection performance scales with $\tau_{\text{p}}T$.
Hence, to decrease the implementation complexity with similar detection performance, one may decrease $T$ and increase $\tau_{\text{p}}$.
However, this will lead to that less time is available to transmit data.

\subsection{Implementation Results\label{sec:results}}
Based on the previous analyses, three different design cases are considered. A design using the Fast algorithm and 256 iterations provides detection results close to the converged floating-point results as shown in Fig.~\ref{fig:ROCfast}. For the Mult algorithm, 192 iterations are used since increasing the number of iterations beyond that only gives a limited detection gain. As discussed earlier, the Mult algorithm converges faster for the good cases, but slower for the bad cases, so depending on channel conditions etc, the actual detection performance may differ, and the different approaches may be beneficial in different scenarios. Finally, the Fast algorithm with 128 iterations gives about the same detection performance as the Mult algorithm with 192 iterations and is also included. The design is identical, and it is only the number of iterations in the power simulations that differ. The corresponding ROC curves are shown in Fig.~\ref{fig:implementedROC}.
\begin{figure}
	\centering
	\includegraphics[scale=.38]{implementedROC.pdf}
	\caption{Missed detection, $p_{\text{m}}$, versus false alarm rates, $p_{\text{fa}}$, of the three considered implementation examples and the converged floating-point results.}
	\label{fig:implementedROC}
\end{figure}

The designs were synthesized for increasing clock frequency and power simulated. The results in terms of area and power are shown in Fig.~\ref{fig:apFastMult}. It can be seen that the area and power is significantly smaller for the Mult algorithm, although the maximum clock frequency is slightly higher for the Fast algorithm, 320~MHz, as compared to 300~MHz for the Mult algorithm\footnote{Can be increased to 320~MHz by using a different programmable shift with a small area increase.}.
\begin{figure}
	\centering
	\subfigure[]{\label{fig:aFastMult}\includegraphics[scale=0.3]{areafastmul.pdf}}
	\hfill
	\subfigure[]{\label{fig:pFastMult}\includegraphics[scale=.3]{newpower.pdf}}  \caption{a) Area and (b) power consumption at different clock frequencies for the implementations of Fast and Mult algorithms. }

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

It may at first seem a bit unintuitive that the power consumption for the Fast approach using 128~iterations is higher than when using 256 iterations. However, it should be noted that these results are against clock frequency, so the approach with 128~iterations can perform twice as many detections. Also, as illustrated in the following, the switching activity and therefore power consumption is reduced when the results converge. Similarly, the Mult approach using 192~iterations can perform more detections compared to the Fast approach with 256~iterations, but fewer than the Fast approach with 128~iterations.

In the following, the 250~MHz clock frequency designs are used. Considering the area curves in Fig.~\ref{fig:apFastMult}, this should be close to the minimum area design, as synthesizing for a lower clock rate will not decrease the area significantly.

To provide further insights here, the average power consumption was analyzed in blocks of 16 iterations. The results are shown in Fig.~\ref{fig:blockpower}. As expected, the power consumption is high in the first iterations and decreases as the results converge.
\begin{figure}
	\centering
	\subfigure[]{\label{fig:blockpower}\includegraphics[scale=.3]{powerBlocksFastMUlFr250.pdf}}
	\hfill
	\subfigure[]{\label{fig:energy}\includegraphics[scale=.3]{energyBlocksFastMulnormalized.pdf}}
	\caption{(a) Average power and (b) accumulated energy consumption over blocks of 16 iterations for Fast and Mult approaches. }
	\label{fig:blockpowerenergy}
\end{figure}

Similarly, one can compute the energy required for one detection instance. This is shown in Fig.~\ref{fig:energy}, where the 16 iteration block averages are used as underlying data. Here, it is clear that the total energy required for one detection is significantly lower for the Mult algorithm compared to the Fast algorithm. Even if 128~iterations are used for the Fast algorithm, about 250~nJ is required, compared to about 100~nJ for the Mult algorithm. Still, the Fast algorithm has better detection performance when more iterations are performed.

Finally, it should be pointed out that the current approaches can perform about one million detections per second. This is most likely too much in the current settings. However, it is possible that the future brings scenarios with higher detection rate requirements. Hence, in the meanwhile one can reduce the energy consumption further by power supply voltage scaling, leading to even lower power and energy consumption.

\subsection{Comparison with Previous Results}
As discussed earlier, to the best of the authors' knowledge, the only other existing previous implementations of mMTC activity detection algorithms are \cite{Tran2019,Henriksson2020}. However, this is for the other type of scenario, non-orthogonal sequences, instead of the pilot-hopping scheme considered here. Also, \cite{Tran2019} considers the case where each user can send different sequences and therefore transmit a few bits of information (four bits). Hence, the comparison is not fully relevant as such.  The reported detection performance is better for both \cite{Tran2019} and \cite{Henriksson2020}. However, this will apart from the algorithm also depend on, e.g., SNR, sequence length, number of users, user activity etc.

The implementations in \cite{Tran2019} and \cite{Henriksson2020} can perform about 11600 and 289 detections per second, respectively, so significantly less than the proposed architectures. The implementation occupies about 5.1~mm$^2$ and 1.14~mm$^2$ chip area for \cite{Tran2019} and \cite{Henriksson2020}, respectively in the same standard cell process, memories not included. This should be compared to about 0.7~mm$^2$ and 0.4~mm$^2$ for the Fast and Mult algorithms, respectively. Hence, both the proposed architectures occupies significantly less chip area.

The energy consumption for one detection is about 0.13~mJ and 4.5~mJ for \cite{Tran2019} and \cite{Henriksson2020}, respectively. The Fast algorithm presented here requires about 250~nJ and 300~nJ for 128 and 256 iterations, respectively. The Mult algorithm requires about 100~nJ. Hence, the proposed approaches are several orders of magnitude more energy efficient.

Another advantage of the approaches in \cite{Tran2019,Henriksson2020} is that the pilot sequences are easy to change. In the approaches proposed in the current work, they are hard coded into the structure of the matrix-vector multiplications with $\mathbf{A}$ and $\mathbf{A}^T$. While this provides a very efficient implementation, it may have other drawbacks. However, one may imagine that a system comes with prepared user modules. This will for example guarantee that all users have unique pilot-hopping sequences. One can also consider to synthesize the proposed architectures to FPGAs, in which case the sequences are readily configurable through a re-synthesis of the $\textbf{A}$ and $\textbf{A}^T$ matrix-vector multiplications\footnote{The spectral radius is about the same for all matrices with the same parameters that have been tested.}.

\section{Conclusion}
\label{sec:conclusions}
Two different architectures for activity detection in grant-free random access massive machine type communication based on pilot-hopping sequences have been proposed. They both solve a non-negative least squares problem, but with different algorithms. Both algorithms are iterative and with deterministic computations, enabling a high-speed implementation. It is shown that the fixed-point implementation of the fast projected gradient algorithm, Fast, provides detection results close to the converged floating-point realization. The multiplicative updates algorithm, Mult, is shown to converge faster for good cases, but slower for bad cases. It is partially implemented in the logarithmic domain to replace multiplications and divisions with additions and subtractions. Here, 192 iterations are selected, which gives a detection performance similar to 128 iterations of the Fast algorithm. However, the area and power/energy consumptions of the multiplicative updates algorithm is significantly lower. The energy consumed for one detection is about 100, 250, and 300 nJ for Mult with 192 iterations, Fast with 128 iterations, and Fast with 256 iterations, respectively, in a scenario with 1024 users and eight pilots, each with a length of 16.

The detection rate is more than one million detections per second, which is most likely too much for contemporary scenarios. However, future scenarios will most likely require an increased detection speed. It is also possible to reduce the clock frequency and power supply voltage to benefit from a lower energy consumption. Currently, the pilot hopping matrices are hard coded in the structure, requiring that the terminals must send the correct sequences. When implementing the same architecture on an FPGA, it is possible to change the pilot hopping matrices by reprogramming the FPGA. This is readily obtained as the architectures have code generators developed.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}