Scalable Spectrum Allocation for Large Networks Based on Sparse Optimization\footnote{This work was supported by a grant from Futurewei.
\IEEEauthorblockN{Binnan Zhuang
Joint allocation of spectrum and user association is considered for a large cellular network. The objective is to optimize a network utility function such as average delay given traffic statistics collected over a slow timescale. A key challenge is scalability: given $n$ Access Points (APs), there are $O(2^n)$ ways in which the APs can share the spectrum. The number of variables is reduced from $O(2^n)$ to $O(nk)$, where $k$ is the number of users, by optimizing over local overlapping neighborhoods, defined by interference conditions, and by exploiting the existence of sparse solutions in which the spectrum is divided into $k+1$ segments. We reformulate the problem by optimizing the assignment of subsets of active APs to those segments. An $\ell_0$ constraint enforces a one-to-one mapping of subsets to spectrum, and an iterative (reweighted $\ell_1$) algorithm is used to find an approximate solution. %based on a reweighted $\ell_1$ approximation. Numerical results for a network with 100 APs serving several hundred users show the proposed method achieves a substantial increase in total throughput relative to benchmark schemes. %Joint allocation of spectrum and user association %is considered over a large cellular network. %The objective is to optimize a network utility function %such as average delay given traffic statistics collected %over a slow timescale. A key challenge is scalability: %interference across cells causes dependencies across the entire %network, making the problem intractable for a large network. %We present an algorithm for approximating the solution %that can be applied to a network with 100 Access Points (APs) %serving several hundred mobiles. %This is achieved by optimizing over local overlapping neighborhoods, %defined by interference conditions, and by exploiting the sparse %nature of the solution. Namely, the spectrum is divided %into $k+1$ segments, where $k$ is the number of mobiles, and %the problem is reduced to optimizing the assignment of different subsets of active APs to the spectrum segments. %A convex relaxation is proposed for this, which relies on a reweighted %$\ell_1$ approximation of an $\ell_0$ constraint. %Numerical comparisons with benchmark schemes show that the proposed %method achieves a substantial increase in achievable throughput. %The algorithm exploits the fact that the solution divides the spectrum into $n$ segments, %where $n$ is the number of APs. Each segment is mapped to a particular set, or {\em pattern, %of active APs. The problem is then to find the mapping of segments to patterns, %and to optimize the width of the segments. Finding a scalable algorithm with $n$ %becomes challenging since the number of patterns grows as $O(2^n)$. % % % In a network with $n$ access points (APs), the spectrum is divided into $2^n$ segments, corresponding to all possible patterns of active APs. The objective is to optimize the widths of the segments jointly with user association. Hence, the resulting optimization problem has $O(2^n)$ variables and becomes computationally intractable for large networks. The number of variables is reduced to $O(nk^2)$ by reformulating the problem over local overlapping neighborhoods. Exploiting the sparse structure of the solution then leads to a convex formulation with $\ell_0$ (cardinality) constraints. An iterative algorithm for finding an approximate solution is developed based on a reweighted $\ell_1$ approximation. Numerical results suggest the solution obtained by the proposed algorithm achieves substantial average packet delay reduction and network throughput gain compared with full-spectrum-reuse even for a network with 100 APs.
isit2017_arxiv.tex
%\documentclass[10pt,conference]{style_file}
\documentclass[conference]{IEEEtran}

\usepackage{cite}
\usepackage{graphicx}
% \usepackage[cmex10]{amsmath}
\usepackage{amsmath}
\usepackage{bbm}
% \usepackage{amsfonts,amssymb,amsbsy}
\usepackage{algorithmicx}
\usepackage{algpseudocode}
\usepackage{algorithm}
\usepackage{dsfont}
\usepackage{color}
\usepackage{xcolor}
\usepackage{pdfcomment}
\usepackage[normalem]{ulem}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\usepackage{tikz}
\usetikzlibrary{positioning,shapes,shapes.geometric,calc,chains,shapes.multipart,arrows}

% \ifCLASSOPTIONcompsoc
%   \usepackage[caption=false,font=normalsize,labelfont=sf,textfont=sf]{subfig}
% \else
\usepackage[caption=false,font=footnotesize]{subfig}
% \fi

\tikzset{
	rubberduck/.style={
		%draw=blue!30,
		shape=isosceles triangle,
		fill=red!30,
		minimum height=18pt,
		minimum width=12pt,
		shape border rotate=#1,
		isosceles triangle stretches,
		inner sep=0pt,
	},
	ap/.style={rubberduck=+90},
	ducky/.style={rubberduck=-90},
	queue/.style={
		rectangle split,
		minimum width=2em,
		%minimum height=2pt,
		%rectangle split vertical,
		rectangle split parts=5,
		draw,
		anchor=south,
	}
}

\makeatletter
\newenvironment{varsubequations}[1]
 {%
  \addtocounter{equation}{-1}%
  \begin{subequations}
  \renewcommand{\theparentequation}{#1}%
  \def\@currentlabel{#1}%
 }
 {%
  \end{subequations}
 }
\makeatother

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}{Proposition}
\newtheorem{corollary}{Corollary}

\DeclareMathOperator*{\minimize}{minimize}
\DeclareMathOperator*{\maximize}{maximize}
\DeclareMathOperator{\subjectto}{subject~to}

\algnewcommand\algorithmicinput{\textbf{INPUT:}}
\algnewcommand\INPUT{\item[\algorithmicinput]}

\algnewcommand\algorithmicoutput{\textbf{OUTPUT:}}
\algnewcommand\OUTPUT{\item[\algorithmicoutput]}

\algnewcommand\algorithmicinit{\textbf{Initialization:}}
\algnewcommand\Init{\item[\algorithmicinit]}

\algnewcommand\algorithmicpost{\textbf{Post Processing:}}
\algnewcommand\Post{\item[\algorithmicpost]}

%\newcommand{\EW}[1]{{\color{magenta}#1}}
\newcommand{\EW}[1]{#1}
\renewcommand{\sout}[1]{}
\newcommand{\bs}{\boldsymbol{s}}
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bv}{\boldsymbol{v}}
\newcommand{\bw}{\boldsymbol{w}}
\newcommand{\bx}{\boldsymbol{x}}
\newcommand{\by}{\boldsymbol{y}}
\newcommand{\bz}{\boldsymbol{z}}
\newcommand{\bh}{\boldsymbol{h}}
\newcommand{\balpha}{\boldsymbol{\alpha}}
\newcommand{\bbeta}{\boldsymbol{\beta}}
\newcommand{\bgamma}{\boldsymbol{\gamma}}
\newcommand{\NN}{\mathcal{N}}
\newcommand{\KK}{\mathcal{K}}
\newcommand{\Ncal}{\mathcal{N}}
\newcommand{\Kcal}{\mathcal{K}}
\newcommand{\NU}{\mathcal{U}}
\newcommand{\NA}{\mathcal{A}}

\begin{document}
\title{Scalable Spectrum Allocation for Large Networks Based on Sparse Optimization\footnote{This
work was supported by a grant from Futurewei.}}

\author{\IEEEauthorblockN{Binnan Zhuang}
\IEEEauthorblockA{Modem R\&D Lab
\\Samsung Semiconductor, Inc.
\\San Diego, CA}
\and
\IEEEauthorblockN{Dongning Guo, Ermin Wei, and Michael L. Honig}
\IEEEauthorblockA{Department of Electrical Engineering and Computer Science
\\Northwestern University
\\Evanston, IL}}

\maketitle

\begin{abstract}
Joint allocation of spectrum and user association is considered for a large cellular network. The objective is to optimize a network utility function such as average delay given traffic statistics collected over a slow timescale. A key challenge is scalability: given $n$ Access Points (APs),
there are $O(2^n)$ ways in which the APs can
share the spectrum. The number of variables is reduced from $O(2^n)$ to $O(nk)$, where $k$ is the number of users, by optimizing over local overlapping neighborhoods,
defined by interference conditions, and by exploiting the existence of
sparse solutions in which the spectrum is divided into $k+1$ segments.
We reformulate the problem by optimizing the assignment of subsets
of active APs to those segments.
An $\ell_0$ constraint enforces a one-to-one mapping of subsets to spectrum,
and an iterative (reweighted $\ell_1$) algorithm is used to find an approximate solution.
%based on a reweighted $\ell_1$ approximation.
Numerical results for a network with 100 APs
serving several hundred users show the proposed method achieves
a substantial increase in total throughput relative to benchmark schemes.

%Joint allocation of spectrum and user association
%is considered over a large cellular network.
%The objective is to optimize a network utility function
%such as average delay given traffic statistics collected
%over a slow timescale. A key challenge is scalability:
%interference across cells causes dependencies across the entire
%network, making the problem intractable for a large network.
%We present an algorithm for approximating the solution
%that can be applied to a network with 100 Access Points (APs)
%serving several hundred mobiles.
%This is achieved by optimizing over local overlapping neighborhoods,
%defined by interference conditions, and by exploiting the sparse
%nature of the solution. Namely, the spectrum is divided
%into $k+1$ segments, where $k$ is the number of mobiles, and
%the problem is reduced to optimizing the assignment of different subsets of active APs to the spectrum segments.
%A convex relaxation is proposed for this, which relies on a reweighted
%$\ell_1$ approximation of an $\ell_0$ constraint.
%Numerical comparisons with benchmark schemes show that the proposed
%method achieves a substantial increase in achievable throughput.

%The algorithm exploits the fact that the solution divides the spectrum into $n$ segments,
%where $n$ is the number of APs. Each segment is mapped to a particular set, or {\em pattern},
%of active APs. The problem is then to find the mapping of segments to patterns,
%and to optimize the width of the segments. Finding a scalable algorithm with $n$
%becomes challenging since the number of patterns grows as $O(2^n)$.
%
%
% In a network with $n$ access points (APs), the spectrum is divided into $2^n$ segments, corresponding to all possible patterns of active APs. The objective is to optimize the widths of the segments jointly with user association. Hence, the resulting optimization problem has $O(2^n)$ variables and becomes computationally intractable for large networks. The number of variables is reduced to $O(nk^2)$ by reformulating the problem over local overlapping neighborhoods. Exploiting the sparse structure of the solution then leads to a convex formulation with $\ell_0$ (cardinality) constraints.  An iterative algorithm for finding an approximate solution is developed based on a reweighted $\ell_1$ approximation. Numerical results suggest the solution obtained by the proposed algorithm achieves substantial average packet delay reduction and network throughput gain compared with full-spectrum-reuse even for a network with 100 APs.
\end{abstract}


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

Heterogeneous cellular networks with dense deployment of access points (APs)
are anticipated to be a major component of 5G networks.
%has been proposed since the Long Term Evolution Advanced (LTE-A)~\cite{3GPP_TR-36.814,hwang2013holistic,nakamura2013trends}.
%In the coming 5G networks, hyper dense deployment is widely accepted in order to meet the increasing traffic demands.
%Two important problems arise with the densification, namely,
A challenge with such dense deployments is interference management.
%Existing studies suggest multi-cell coordination,
%and even network-wide global radio resource management is the key to solve such problems.
Coordinated radio resource allocation across multiple cells
is one approach for mitigating inter-cell interference.
That includes
%In time domain, it can be achieved
joint scheduling across multiple cells
over a fast timescale~\cite{yu2013multicell,wang2013joint,HonLuo13TransSP}
as well as the assignment of resources across cells
over a slower timescale~\cite{zhuang2014traffic,kuang2014joint}.
Whereas the former approach requires instantaneous knowledge of channel gains,
the latter approach relies on statistical knowledge of interference.
%however, that generally requires knowledge of instantaneous channels gains.
%In contrast, mitigation of average interference over a slow timescale
%can be achieved through
%frank2010cooperative,
%where interference management is often optimized together with power control and beamforming. In frequency \EW{domain}, the corresponding technique is dynamic spectrum allocation~\cite{
%stolyar2008self,
%chang2009multicell, %ChaTao09ICC,
%ali2009dynamic,
%madan2010cell,
%liao2014base
%},
%which \EW{is} usually studied for orthogonal frequency-division multiple access (OFDMA) systems.
%assignment of users to Access Points (APs) over a slower timescale
%can also
%take into account the interference environment (e.g., see
%is traditionally realized by assigning each UE to the AP with the strongest signal, which is known as maximum reference signal receive power (maxRSRP) association. The so called {\em range extension} technique~
%\cite{KhaBhuJi10EW,DamMonYon11TransWC,
%has been introduced to add a positive bias to the pico cell RSRP for load balancing between macro and pico tiers. More sophisticated user association can be achieved by solving non-convex mixed-integer programs~
%shen2014distributed,fooladivanda2013joint,).
%KuaSpeDro12VTC,Kua14arXiv,LinYu13GLOBE,YeRonChe13TransWC}.


%To achieve the expected capacity enhancement in densely deployed networks, the key is interference management. Coordinated radio resource allocation across multiple cells has been proved effective in mitigating inter-cell interference. In time domain, the multi-cell joint scheduling have been studied in~\cite{yu2013multicell,frank2010cooperative,wang2013joint}, where the interference aware scheduling is often optimized with power control and beamforming. In the frequency domain, the dynamic spectrum allocation problem has been studied in~\cite{
%stolyar2008self,
%chang2009multicell, %ChaTao09ICC,
%ali2009dynamic,
%madan2010cell,
%liao2014base
%}

We consider the joint allocation of spectrum and user association in a large network
with many APs. The objective is to optimize a network utility, such as average delay,
given traffic statistics over a geographic region that change slowly relative
to channel fading. Our approach builds on our prior work
%over a slow timescale framework, given  proposed in~
\cite{zhuang2016scalable,zhuang2014traffic} in which for a network of $n$ APs
and $k$ mobiles (or User Equipments (UEs)), the spectrum is partitioned
into $2^n$ {\em patterns}, corresponding to all possible subsets
of active APs. The problem is to optimize
the widths of spectrum segments, associated with the different patterns,
along with the user association under each pattern.
This has been shown to provide significant performance improvement
in throughput enhancement and delay reduction
\cite{zhuang2014traffic,zhuang2016scalable,kuang2014joint,kuang2016energy}.

The original convex problem formulation in \cite{zhuang2014traffic}
is not useful for large networks because the number of patterns
grows as $O(2^n)$.
%Although it is proved in~\cite{zhuang2014traffic} and~\cite{2015energy} that at most $k$ out of $2^n$ patterns will be \EW{active} in the delay minimizing solution, it is in general NP-hard to decide the active patterns.
In prior work~\cite{zhuang2016scalable}, we have reduced the number of variables
to $O(n)$ by recognizing that each link rate depends only on {\em local} patterns
of active APs. The problem can then be redefined over sets of overlapping
{\em interference neighborhoods}, associated with those local patterns.
Each AP has its own interference neighborhood, which
captures the interference from nearby APs. The challenge with
the approach in~\cite{zhuang2016scalable} is to ensure
that the spectrum assigned to each particular AP
is consistent across the neighborhoods to which it belongs.
To accomplish that, the spectrum is discretized and a coloring algorithm is
proposed to ensure that the local patterns of active APs are globally consistent.
%Even with those constraints, however, it may not be possible to find
%consistent {\em placements} of the spectrum assignments for each AP
%within the available band. A discrete
%i.e., patterns of the APs within range. A local cluster based formulation reduces the number of relevant patterns at each AP from $O(2^n)$ to $O(1)$. A discrete coloring algorithm is then presented to ensure that the local patterns are globally consistent.

Here we take a different approach to addressing the scalability problem.
This is based on the fact that the solution is {\em sparse},
meaning that at most $k+1$ out of the $2^n$ possible patterns
appear in the optimal allocation for a general network utility function.
Hence we reformulate the problem by dividing the spectrum into $k+1$ segments,
rather than $2^n$, and attempt to identify the pattern
that should be associated with each segment. This effectively
reverses the approach in~\cite{zhuang2016scalable}, which attempts
to assign a segment of spectrum to each pattern.
In this reformulation, we initially assume that any combination of patterns
can be assigned to each of the $k+1$ segments.
This problem is a convex relaxation of our original problem.
The one-to-one mapping of spectrum segments to patterns is then enforced
with an $\ell_0$ (cardinality) constraint.
An algorithm for finding an approximate solution to this problem is presented
based on a reweighted $\ell_1$ approximation of this constraint \cite{candes2008enhancing}.

The approach to scalability presented here has the advantage
of eliminating the combinatorial coloring problem that arises in~\cite{zhuang2016scalable}.
Instead, a new $\ell_0$ constraint is introduced. Although this does not
simplify the original problem, it helps in finding an approximate solution,
since the reweighted $\ell_1$ approximations for the $\ell_0$ constraints
are known to perform well.
%Second, the number of variables is reduced
%from $O(2^n)$ in~\cite{zhuang2016scalable} to $O(nk^2)$, facilitating scalability.
Numerical results indicate that this method generally gives better performance
for a fixed computational complexity than the method in~\cite{zhuang2016scalable}.
% We also show that how the centralized iterative algorithm presented here can also be decomposed into subproblems for parallel computing based on the Alternating Directions Method of Multipliers (ADMM)~\cite{boyd2011distributed}.

%The rest of the paper is organized as follows.
%The system model is presented in Section~\ref{sec:SysMod}. The original formulation with $2^n$ global patterns is shown in Section~\ref{sec:GlobalForm}. The scalable formulation based on sparsity constraint is shown in Section~\ref{sec:MainRes}. An efficient iterative algorithm
% and its distributed solution based on ADMM are
%is introduced in Section~\ref{sec:l1}.
% and Section~\ref{sec:ADMM}, respectively.
%The simulation results are presented in Section~\ref{sec:NumRes}. Section~\ref{sec:Con} concludes the paper.


\section{System Model}
\label{sec:SysMod}

% \begin{figure}
% 	\centering
%   \small
% %	\resizebox{\columnwidth}{!}{%
% 	\begin{tikzpicture}
% 	\begin{axis}[
% 	xbar stacked,
% 	ytick=data,
% 	y axis line style= { draw opacity=0 },
% 	axis y line*=none,
% 	axis x line*=bottom,
% 	xtick={0,1},
% 	width=3.75in, %textwidth,
% 	bar width=6mm,
% 	% xlabel={Time in Seconds},
% 	yticklabels={AP3,AP2,AP1},
% 	xmin=0,
% 	xmax=1,
% 	% area legend,
% 	y=8mm,
% 	enlarge y limits={abs=0.625},
% 	]
% 	\addplot[white,fill=white] coordinates
% 	{(0.1,0) (0.1,1) (0,2)};
% 	\addplot[white,fill=white] coordinates
% 	{(0.15,0) (0,1) (0,2)};
% 	\addplot[cyan,fill=cyan] coordinates
% 	{(0,0) (0,1) (0.1,2)};
% 	\addplot[cyan,fill=cyan] coordinates
% 	{(0,0) (.15,1) (.15,2)};
% 	\addplot[cyan,fill=cyan] coordinates
% 	{(0.2,0) (0.2,1) (.2,2)};
% 	\addplot[cyan,fill=cyan] coordinates
% 	{(.25,0) (0,1) (.25,2)};
% 	\addplot[white,fill=white] coordinates
% 	{(0,0) (0.25,1) (0,2)};
% 	\addplot[cyan,fill=cyan] coordinates
% 	{(0,0) (0.075,1) (0,2)};
% 	\addplot[white,fill=white] coordinates
% 	{(0.075,0) (0,1) (0.075,2)};
% 	\addplot[cyan,fill=cyan] coordinates
% 	{(0.125,0) (0.125,1) (0,2)};
% 	\addplot[cyan,fill=cyan] coordinates
% 	{(0.1,0) (0,1) (0,2)};
%   %\addplot +[mark=none] coordinates {(.1, -.1) (.1, 1)};
%   \addplot[black] coordinates {(0.1,0) (0.1,1)};
% 	\end{axis}%
% 	\node at (2.53,-.3) {$y^{\{1,2\}}$};
% 	\draw (1.45,-.2) -- (1.45,2.5);
% 	\node at (0.725,-.3) {$y^{\{1\}}$};
% 	\draw (3.61,-.2) -- (3.61,2.5);
% 	\node at (5.045,-.3) {$y^{\{1,2,3\}}$};
% 	\draw (6.48,-.2) -- (6.48,2.5);
% 	\node at (8.275,-.3) {$y^{\{1,3\}}$};
% 	\draw (10.07,-.2) -- (10.07,2.5);
% 	\node at (10.61,-.3) {$y^{\{2\}}$};
% 	\draw (11.15,-.2) -- (11.15,2.5);
% 	\node at (12.045,-.3) {$y^{\{2,3\}}$};
% 	\draw (12.94,-.2) -- (12.94,2.5);
% 	\node at (13.72,-.3) {$y^{\{3\}}$};
% 	\draw[dashed] (8.31,0.12) -- (8.31,.96);
% 	\draw[dashed] (7.75,1.75) -- (7.75,2.5);
% 	\node at (7.115,2.1) {UE1};
% 	\node at (7.115,2.82) {$w_{1\to 1}^{\{1,3\}}$};
% 	\node at (8.91,2.1) {UE2};
% 	\node at (8.91,2.82) {$w_{1\to 2}^{\{1,3\}}$};
% 	\node at (7.395,.5) {UE1};
% 	\node at (7.395,1.2) {$w_{3\to 1}^{\{1,3\}}$};
% 	\node at (9.19,.5) {UE2};
% 	\node at (9.19,1.2) {$w_{3\to 2}^{\{1,3\}}$};
% 	\draw[decorate, decoration={brace}, very thick, blue] (6.48,2.44) -- (7.75,2.44);
% 	\draw[decorate, decoration={brace}, very thick, blue] (7.75,2.44) -- (10.07,2.44);
% 	\draw[decorate, decoration={brace}, very thick, blue] (6.48,0.83) -- (8.31,.83);
% 	\draw[decorate, decoration={brace}, very thick, blue] (8.31,0.83) -- (10.07,.83);
% \end{tikzpicture}%}
% 	\caption{Illustration of all patterns of a 3-AP 2-UE network. The allocations to 2 UEs are revealed under pattern $\{1,3\}$.}
% 	\label{fig1}
% \end{figure}
\vspace{-5pt}
\begin{figure}%[!t]
\centering
\includegraphics[width=0.9\columnwidth]{specbar}\vspace{-5pt}
\caption{An example off spectrum allocation among 3 APs.}
\label{fig:Mod}
\end{figure}

As in~\cite{zhuang2016scalable}, we consider a network containing
the set of APs $\NN=\{1,\cdots,n\}$ and the set of UEs $\KK=\{1,\cdots,k\}$.
Each ``UE'' is actually associated with a particular location, and could refer
to a group of nearby mobiles.
The $n$ APs share $W$ Hz of spectrum. For convenience, we normalize $W=1$.
Each AP can transmit on any part(s) of the spectrum, which is homogeneous
(has the same distance attenuation). APs sharing the same spectrum
interfere at UEs within range of both transmitters. All transmissions are assumed
to be omni-directional, although a similar set of problems can be reformulated
with directional transmissions.

We define a subset of active (transmitting) APs
$A\subset\NN$ as a {\em pattern}.
There are $2^n$ patterns, and a particular allocation of spectrum maps
each pattern to a slice of spectrum.
We denote an allocation as $\{y_A\}_{A\subset \NN}$,
where $y_A$ denotes the amount of bandwidth assigned to pattern $A$.
An example with three APs is depicted in Fig.~\ref{fig:Mod}.
AP 1 \EW{owns ${\{1\}}$ exclusively}; shares ${\{1,2\}}$ with AP 2;
shares ${\{1,3\}}$ with AP 3; and shares ${\{1,2,3\}}$ with both AP 2 and AP 3.

User association is determined by how the spectrum assigned to a particular AP
is allocated to different UEs.  \EW{Specifically,} denote the fraction of total bandwidth
used by AP $i$ to serve UE $j$ under pattern $A$ as $x_A^{i\to j}$,
%which is only defined
for $i\in A$.
UE $j$ is then assigned to AP $i$ if $x_A^{i\to j}>0$ for some $A \subset \mathcal{N}$.
%meaning AP $i$ transmits to UE $j$ under some pattern(s).
Since the total bandwidth assigned to pattern $A$ is $y_A$, we have
\begin{align}
\label{eq:sum_x=y}
\sum_{j\in\KK}x^{i\to j}_A \leq y_A,~\forall A\subset\NN,\;i\in A.
\end{align}
%(We do not allow multicasting to multiple UEs using the same spectrum.)
The total bandwidth allocated to all patterns
%(including the empty set $\emptyset$)
is then
\begin{align}
\label{eq:bandwidth}
\sum_{A\subset \NN}y_A=1.
\end{align}
%On a slow timescale, we assume the average spectral efficiencies under different patterns are known constants or can be measured.

Let $s_A^{i\to j}$ denote the spectral efficiency of the link from AP $i$ to UE $j$ under pattern $A$.
This is measured over a {\em slow timescale}, and is therefore an average over short-term fading.
The value of $s_A^{i\to j}$ is therefore determined by the distance between AP $i$ and UE $j$,
shadowing, and similar long-term characteristics of the interference links from APs in $A$
to UE location $j$.
For concreteness
%in obtaining numerical results,
we assume
\begin{align} \label{eq:SE}
  s^{i\rightarrow j}_A
  =
  \frac{W}{\tau}\mathbbm{1}_{\{i\in A\}}
  \log_2\left(1+\frac{p_ig^{i\to j}}{\sum\limits_{i'\in A\setminus\{i\}}p_{i'}g^{i'\to j}+n_j}\right) % &~\text{packets/second}
\end{align}
%\noindent
%packets per second,
where $\mathbbm{1}_{\{i\in A\}}=1$ if $i\in A$ and $0$ otherwise,
$p_i$ is the transmit power spectral density (PSD) at AP $i$, $g^{i\to j}$ is the power gain of link
{$i\to j$}, and $n_j$ is the noise PSD
%\EW{Define PSD please}
at UE $j$. We assume fixed, flat transmit PSDs over the slow timescale considered.
The factor $W/\tau$, where $\tau$ is the average packet length (bits), gives
the units in packets/sec. The link gain $g^{i\to j}$ includes pathloss and shadowing effects.
%In a large network the link gain vanishes with distance.
%\EW{Why is this sentence here? Are we saying  ``For instance, in a large network, the link gain vanishes with distance, then $g^{i\to j}$ will be 0''?}
Clearly, $s_A^{i\to j}=0$ if $i\not\in A$, i.e., AP $i$ does not transmit on pattern $A$.
In practice, the spectral efficiencies can be measured as time-averaged channel gains.
The total rate received by UE $j$ is therefore
\begin{align}
\label{eq:rate}
r^j=\sum_{A\subset \NN}\sum_{i\in A}s_A^{i\to j}x^{i\to j}_A,~\forall j\in\KK.
\end{align}


\section{Problem Formulation with Global Patterns}
\label{sec:GlobalForm}
The problem is to maximize a network utility function over the spectrum allocation and user association
designated by
$\bx=\big(x^{i\to j}_A\big)_{A\subset \NN, i\in A, j\in\KK}$\EW{, $\by = (y_A)_{A\subset \NN}$:}
%\section{The Optimization Problem}
%\label{sec:OptProb}
%\subsection{Problem formulation}
%The optimization problem based on global patterns is formulated as~\ref{eq:Opt}:
\begin{varsubequations}{P0}
\label{eq:Opt}
\begin{align}
  % \underset{\br,\;\bx,\;\by}{\text{maximize}}~&~ u(r^1,\cdots,r^k)&\label{eq:Obj-Opt}\\
  \maximize_{\br,\;\bx,\;\by}~
  & u(r^1,\cdots,r^k) \label{eq:Obj-Opt}\\
%\text{subject to}~&~r_j=\sum_{i\in\NN}\sum_{A\subset N} s_{i\to j,A}x_{i\to j,A},&~\forall j\in\KK \label{eq:Con1-Opt}\\
%&~r_j>\lambda_j,~\forall j\in\KK\label{eq:Con2-Opt}\\
%&~y_A=\sum_{j\in\KK} x_{i\to j,A}&~\forall A\subset N,\;i\in A \label{eq:Con2-Opt}\\
  \subjectto~&
  %& r^j=\sum_{A\subset\NN}\sum_{i\in A} s_A^{i\to j} x^{i\to j}_A,
  %\quad \forall j\in\KK \label{eq:Con1-Opt}\\
  %& \sum_{j\in\KK}x^{i\to j}_A \le y_A,
  %\quad \forall A\subset\NN,\;i\in A\label{eq:Con2-Opt}\\
  %& \sum_{A\subset\NN}y_A=1 \label{eq:Con3-Opt}\\
   x^{i\to j}_A\geq 0,
  \quad \forall A\subset \NN,\;\forall i\in A,~\forall j\in \KK \label{eq:Con4-Opt}
\end{align}
\end{varsubequations}\noindent
and constraints \eqref{eq:rate}, \eqref{eq:sum_x=y}, and \eqref{eq:bandwidth},
where the network utility function $u$ depends on the service rates to all UEs.
%\EW{The constraint~\eqref{eq:Con2-Opt} is based on~\eqref{eq:sum_x=y} and the observation that the inequality will be tight at optimal allocation.}
The optimization problem is convex if $u$
is concave in $\br=[r^1,\cdots,r^k]$. \sout{and element-wise non-decreasing.}
%\EW{Monotonicity is not necessary for convexity}
%The class of utility functions that makes~\ref{eq:Opt} convex include: sum rate, minimum user service rate, sum log-rate.
As in~\cite{zhuang2016scalable}, we will take $u$ to be the average packet delay, given by
%minimizing the average packet delay, which is given by:
\begin{align}\label{eq:AvgDelay}
  u(r^1,\cdots,r^k)=-\sum_j\frac{\lambda^j}{(r^j-\lambda^j)^+}
\end{align}\noindent
where $\frac{1}{(x)^+}$ equals $\frac{1}{x}$ if $x>0$ and $+\infty$ otherwise, and $\lambda_j$ is the Poisson packet arrival rate for UE $j$.
This assumes exponential packet lengths and backlogged interference~\cite{zhuang2014traffic}.

Solving~\ref{eq:Opt} becomes prohibitively expensive % or infeasible, not sure which one reads better...}
as the network size grows, due to the inherit complexity from the $2^n$ global patterns.
%We next reformulate the problem with $O(nk^2)$ variables by defining local interference neighborhoods.
However, a key property of solution(s)
to \eqref{eq:Opt} is that at least one is {\em sparse}, i.e., contains at most $k+1$ patterns.
%\EW{This sparsity property will serve as the building block for development of scalable algorithms.}
\begin{proposition}{(\cite{zhuang2015energy})}
\label{thm:Spec}
\ref{eq:Opt} has a solution that divides
the spectrum into at most $k+1$ segments, i.e.,
\begin{align} \label{eq:ThmSpec}
  \left|\{A\subset \NN~|~y_A>0\}\right|\leq k+1.
\end{align}\nonumber
Furthermore, if $u$ is element-wise nondecreasing in $\br$,
$k+1$ in \eqref{eq:ThmSpec} is reduced to $k$.
\end{proposition}
%This property forms the basis of our approach to solving~\ref{eq:Opt}.
%The sparse optimal solution described in Proposition~\ref{thm:Spec} is suitable for practical implementation.
Determining which $k+1$ active patterns appear in a solution
is then the key challenge in solving~\ref{eq:Opt}.

\section{Reformulation With Sparsity Constraints}
\label{sec:MainRes}
\subsection{Local Neighborhoods}
\label{subsec:local}
We first reduce the number of variables in \eqref{eq:Opt} from $O(2^n)$ to $O(n)$
by reformulating \eqref{eq:Opt} over {\em interference neighborhoods}
based on local patterns~\cite{zhuang2016scalable}. Due to pathloss, we assume interference vanishes
beyond a certain distance. Let $L\subset \NN\times \KK$ denote the set of links with nonzero gains.
%\footnote{$\NN \times \KK$ is the Cartesian product of $\NN$ and $\KK$.},
Hence each UE only receives power from APs within its neighborhood:
\vspace{-0.2cm}\begin{align}
\label{eq:SetA}
\mathcal{A}_j=\{i|(i\to j)\in L\}.
\end{align}
%\vspace{-0.1cm}
From the AP side, each AP can only transmit to a collection of UEs in a AP neighborhood with positive rates:
\vspace{-0.1cm}\begin{align}
\label{eq:SetU}
\mathcal{U}_i=\{j|(i\to j)\in L\}.
\end{align}\vspace{-1pt}
We define the interference neighborhood for each AP $i$ as:
%\vspace{-5pt}
\begin{align}
\label{eq:SetN}
\mathcal{N}_i = \cup_{j\in\mathcal{U}_i}\mathcal{A}_j,
\end{align}
i.e., $\mathcal{N}_i$ includes AP $i$ and all APs that interfere with it.

An example with three APs (denoted by $\{1,2,3\}$) and 2 UEs (denoted by $\{a,b\}$) is shown in Fig~\ref{fig:LocalPattern}. The set of nonzero links are $L=\{1\to a,\;2\to a,\;2\to b,\;3\to b\}$.
The AP neighborhoods are $\mathcal{U}_1=\{a\}$, $\mathcal{U}_2=\{a,b\}$, and $\mathcal{U}_3=\{b\}$; and the UE neighborhoods are $\mathcal{A}_a=\{1,2\}$ and $\mathcal{A}_b=\{2,3\}$. AP 1's interference neighborhood is $\mathcal{N}_1=\{1,2\}$, as AP 2 interferes with it at UE $a$. AP 3's interference neighborhood is $\mathcal{N}_3=\{2,3\}$, as AP 2 interferes with it at UE $b$. AP 2's interference neighborhood is $\mathcal{N}_2=\{1,2,3\}$, since AP 1 and AP 3 interfere with AP 2 at UE $a$ and UE $b$, respectively.

\begin{figure}%[!t]
\centering
\includegraphics[width=.6\columnwidth]{venn} %LocalPattern}
\caption{Neighborhoods in the case of three APs and two UEs.}
\label{fig:LocalPattern}
\end{figure}

The spectral efficiency of link $i\to j$ only depends on the local patterns in $\mathcal{A}_j$, according to the definition of UE neighborhoods in~\eqref{eq:SetA}. Therefore, for any link $i\to j$, \EW{the spectral efficiency under }a global pattern $A$ is equivalent to its intersection with UE $j$'s neighborhood $\mathcal{A}_j$:
\vspace{-5pt}\begin{align}
\label{eq:LocalSE}
s^{i\to j}_A=s^{i\to j}_{A\cap \mathcal{A}_j} ~\forall j\in\KK,\; \forall i\in \mathcal{A}_j,\; A\subset\NN.
\end{align}

%\EW{Having related the spectral efficiency for a global pattern to a local pattern,
\vspace{-5pt}We next express each rate in terms of local patterns.
%From AP $i$'s point of view, \sout{all the patterns within its interference cluster must be differentiated.}
%\EW{patterns need to be differentiated if they overlap with different interference neighborhoods.}  For example, in Fig.~\ref{fig:LocalPattern}
%AP 2 can serve both UEs $a$ and $b$. Although $\{1\}$ and $\{1,3\}$ are equivalent from UE $a$'s perspective, they are different from UE $b$'s perspective as $s^{3\to b}_{\{1\}}=0$ and $s^{3\to b}_{\{1,3\}}>0$.
We define bandwidth allocation variables within a local interference neighborhood $\NN_i$ as:
% as the fraction of the bandwidth allocated to link $i\to j$ under the local pattern:
\vspace{-5pt}\begin{align}
\label{eq:LocalGlobal}
z^{i\to j}_B=\sum_{A\subset\NN:A\cap\NN_i=B}x^{i\to j}_A,~i\in\NN,\;j\in\mathcal{U}_i,\;\forall B\subset\NN_i.\vspace{-5pt}
\end{align}
\sout{If we treat an interference cluster $\NN_i$ as a partition of all the global patterns, many global patterns are equivalent.}
%\EW{%Not really a partition... Try this ``
%For any
Because all global patterns sharing the same overlap with $\NN_i$ contribute to the same local pattern of AP $i$, $z^{i\to j}_B$ represents the bandwidth allocated to link $i\to j$ under local pattern $B$.
% , then locally (for AP $i$) they are treated as equivalent and can be combined.''}
Hence from AP $i$'s perspective, the bandwidth assigned to a local pattern $B$ must be the sum of the bandwidths assigned to all global patterns containing $B$.
%Given~\eqref{eq:LocalGlobal} and~\eqref{eq:LocalSE},
The service rate of link $i\to j$, as defined in~\eqref{eq:rate}, can then be calculated
as \cite{zhuang2016scalable}:
%\begin{subequations}
\begin{equation}
r^j
%&=\sum_{A\subset \NN}\sum_{i\in A}s_A^{i\to j}x^{i\to j}_A\label{eq:eq1}\\
%&=\sum_{A\subset \NN}\sum_{i\in A}s^{i\to j}_{A\cap\mathcal{A}_j}x^{i\to j}_{A}\label{eq:eq2}\\
%&=\sum_{i\in \mathcal{A}_j}\sum_{B\subset\NN_i}\sum_{A\subset\NN:A\cap\NN_i=B}s^{i\to j}_{A\cap\mathcal{A}_j}x^{i\to j}_{A}\label{eq:eq3}\\
=\sum_{i\in \mathcal{A}_j}\sum_{B\subset\NN_i}s^{i\to j}_{B\cap\mathcal{A}_j}z^{i\to j}_B % \label{eq:eq4}
\label{eq:LocalRate}
\end{equation}%
%\end{subequations}%
% From~\eqref{eq:eq1} to~\eqref{eq:eq2}
where we have used~\eqref{eq:LocalSE} and~\eqref{eq:LocalGlobal},
%to arrive at~\eqref{eq:eq2},%. From~\eqref{eq:eq2} to
%~\eqref{eq:eq3} is because
%any $\NN_i$ can be used as a partition of all global patterns, and
and the fact that only APs in $\mathcal{A}_j$ transmit to UE $j$ with positive rate.
%. From~\eqref{eq:eq3} to~\eqref{eq:eq4}
%and we have used%~\eqref{eq:LocalSE} and
%~\eqref{eq:LocalGlobal} to obtain~\eqref{eq:LocalRate}.

\subsection{Sparse Optimization}%  with $O(nk^2)$ variables}
\EW{Motivated by Proposition \ref{thm:Spec},
%i.e., only $k+1$ out of $2^n$ patterns will be active in the optimal solution,
we reformulate \eqref{eq:Opt} by dividing the spectrum into $k+1$ segments,
and seek to assign a single pattern to each segment.} Proposition \ref{thm:Spec}
implies that by optimizing this assignment we obtain a solution to \ref{eq:Opt}. \sout{The number of local patterns in each local interference cluster is bounded despite the size of the network. The general idea is provided here.}
%\EW{To formulate this equivalent problem, we first} divide the spectrum into $k$ segments,
Let $h_l$ denote the bandwidth of the $l$th segment (to be optimized).
We now associate a set of variables $z$ (for local neighborhoods)
and $y$ {\em for each segment} $l$.
We can therefore rewrite the rate in~\eqref{eq:LocalRate} as
%\vspace{-0.2cm}
\begin{align}
\label{eq:Con1}
r^j=\sum_{l=1}^{k+1}\sum_{i\in\mathcal{A}_j}\sum_{B\subset\NN_i} s^{i\to j}_{B\cap\mathcal{A}_j}\;z^{i\to j}_{B,l},~\forall j\in \KK
\end{align}
where the first sum is over the $k+1$ spectrum segments.
Apart from the addition of segment index, the $z$ variable is the same local bandwidth allocation variable defined in~\eqref{eq:LocalGlobal}.
The amount of bandwidth assigned to local pattern $B$ in AP $i$'s interference neighborhood
$\NN_i$ within segment $l$ is
%\vspace{-0.2cm}
\begin{align}
\label{eq:Con2}
y^i_{B,l}=\sum_{j\in\mathcal{U}_i}z^{i\to j}_{B,l}, ~\forall i\in\NN,\;\forall B\subset\NN_i .
%l=1,\cdots,k+1
\end{align}
The total amount of spectrum assigned to $\NN_i$ satisfies
%\vspace{-0.2cm}
\begin{align}
\label{eq:Con5}
\sum_{B\subset \NN_i}y^i_{B,l}\leq h_l,~\forall i.
\end{align}

We also introduce the following {\em consistency constraint} to ensure that
the amount of spectrum allocated to an AP is consistent across  any
two neighborhoods $\NN_i$ and $\NN_m$ that contain it~\cite{zhuang2016scalable}.
That is, for any
%$l=1,\cdots,k+1$, $i\in\NN$, $m\in \NN_i$, and
nonempty $C\subset\NN_i\cap\NN_m$, %  $C\neq\emptyset$,
%\vspace{-0.2cm}
\begin{align}
\label{eq:Con3}
\sum_{B\subset \NN_i:B\cap\NN_m=C}y^i_{B,l}=\sum_{B\subset \NN_m:B\cap\NN_i=C}y^m_{B,l}.
% , &~\forall i\in\NN_m,\;\forall m\in \NN_i,\;\forall C\subset\NN_i\cap\NN_m,\;C\neq\emptyset,\;\forall l\in \KK
\end{align}
See, for example, $\NN_1$ and $\NN_3$ in Fig.~\ref{fig:LocalPattern}. In interference cluster $\NN_1$, AP 2 transmits under pattern $\{2\}$ and $\{1,2\}$; in interference cluster $\NN_3$, AP 2 transmits under pattern $\{2\}$ and $\{2,3\}$. The total bandwidth used by AP 2 must be consistent across
neighborhoods so that
$y^1_{\{2\},l}+y^1_{\{1,2\},l}=y^3_{\{2\},l}+y^3_{\{2,3\},l}$, where $\{2\}$ is the overlapping pattern.
This example can be extended to a {\em set} of APs, which are members of two
interference neighborhoods, giving~\eqref{eq:Con3}.

To ensure a one-to-one mapping of patterns to the $k+1$ segments,
% the optimal structure described in Proposition~\ref{thm:Spec},
we add the $\ell_0$-norm constraint
%\vspace{-0.2cm}
\begin{align}
\label{eq:Con6}
\sum_{B\subset\NN_i}|y^i_{B,l}|_0\leq1,&~\forall i\in\NN,\; l=1,\cdots,k+1
\end{align}
where $|x|_0=1$ if $x\neq 0$, and $|x|_0=0$ if $x=0$.
That is, we constrain each AP to use at most one active pattern in each segment.
%Combining~\eqref{eq:Con3} and~\eqref{eq:Con6}, there is only one active pattern on each segment.

We can now reformulate \ref{eq:Opt} in terms of the local interference variables $\bz$
across the $k+1$ spectrum segments:
%Constraint~\eqref{eq:Con6-Opt3} limits the number of active patterns on each segment to be 1, where $|x|_0=1$ if $x\neq 0$, and $|x|_0=0$ if $x=0$. Intuitively, we force each AP to use at most one active pattern over each segment. Combining~\eqref{eq:Con3-Opt3} and~\eqref{eq:Con5-Opt3}, the APs are forced to use a single unified pattern on each segment. Hence, a feasible solution to~\ref{eq:Opt3} always satisfies Proposition~\ref{thm:Spec}.
%Based on~\eqref{eq:Con1}--\eqref{eq:Con6},~\ref{eq:Opt} is reformulated as~\ref{eq:Opt3}:
\begin{varsubequations}{P1}
\label{eq:Opt3}
%\vspace{-0.3cm}
\begin{align}
&\maximize_{\br,\;\by,\;\bz} \;
u(r^1,\cdots,r^k) & \label{eq:Obj-Opt3}\\
& \subjectto \;
%  & r^j=\sum_{l=1}^k\sum_{i\in\mathcal{A}_j}\sum_{B\subset\NN_i} s^{i\to j}_{B\cap\mathcal{A}_j}\;z^{i\to j}_{B,l},
%  \; \forall j\in \KK \label{eq:Con1-Opt3}\\
%  & \sum_{j\in\mathcal{U}_i}z^{i\to j}_{B,l} \le y^i_{B,l},
%  \quad \forall i\in \NN,\;\forall B\subset\NN_i,\;l=1,\cdots,k+1 \label{eq:Con2-Opt3}\\
%  & \sum_{B\subset \NN_i:B\cap\NN_m=C}y^i_{B,l}
%  = \sum_{B\subset \NN_m:B\cap\NN_i=C}y^m_{B,l}, \nonumber\\
%  & \qquad\quad l=1,\cdots,k+1,\;\forall i\in\NN,\;\forall m\in \NN_i,\nonumber\\
%  & \qquad\quad \forall C\subset\NN_i\cap\NN_m:\, C\neq\emptyset \label{eq:Con3-Opt3}\\
%  & \sum_{B\subset \NN_i}y^i_{B,l}\leq h_l,
%  \quad \forall i\in\NN,\;l=1,\cdots,k+1 \label{eq:Con5-Opt3}\\
%  & \sum_{B\subset\NN_i}|y^i_{B,l}|_0\leq1,
%  \quad \forall i\in\NN,\;l=1,\cdots,k+1 \label{eq:Con6-Opt3}\\
 % y_{B,l}^i\geq 0, ~~
%  \quad \forall B\subset\NN_i,\;l=1,\cdots,k+1 \label{eq:Con8-Opt3}\\
  z^{i\to j}_{B,l}\geq 0,~~
  \forall (i\to j)\in L,\;
  \forall B\subset\NN_i,\;
  %l=1,\cdots,k+1
  \label{eq:Con9-Opt3}\vspace{-2cm}\\
  &
  %l=1,\cdots k+1,
  \qquad \sum_{l=1,\cdots,k+1}h_l=1, ~~~h_l\geq0,
%   ~~~ l=1,\cdots,k+1,
   \label{eq:Con4-Opt3}\vspace{-0.2cm}
\end{align}
\end{varsubequations}\noindent
and constraints \eqref{eq:Con1}-\eqref{eq:Con6} for each $l=1,\cdots,k+1$.
%, \eqref{eq:Con2}, \eqref{eq:Con5}, \eqref{eq:Con3},
%and \eqref{eq:Con6}.
%The equivalence between~\ref{eq:Opt3} and~\ref{eq:Opt} is given in the following theorem.
\begin{theorem}
\label{thm:main}
\ref{eq:Opt} and~\ref{eq:Opt3} are equivalent (have the same solutions)
given the local neighborhood definitions~\eqref{eq:SetA}-\eqref{eq:SetN}.
\end{theorem}

The proof is omitted due to limited space.
Hence a solution to~\ref{eq:Opt3} always satisfies Proposition~\ref{thm:Spec}.
This reformulation leads to a computationally efficient approximation algorithm.


\section{Iterative %$\ell_1$
Approximation}
\label{sec:l1}
The number of variables is reduced from $O(2^n)$ in~\ref{eq:Opt} to $O(nk)$ in~\ref{eq:Opt3}. However, the $\ell_0$ norm constraint makes the problem non-convex and difficult to solve.
We apply the reweighted $\ell_1$ approach, described in~\cite{candes2008enhancing},
to approximate the $\ell_0$ constraint.
Specifically, the $\ell_0$ norm is approximated by the weighted $\ell_1$ norm, $\sum_i w_i|x_i|$,
where $w_i$ is iteratively adapted, and taken to be inversely proportional
to $|x_i|$ computed at the preceding iteration.
This has the effect of suppressing small nonzero entries with large
weights.\footnote{This algorithm has also been used in~\cite{zhuang2015energy}
to solve an AP activation problem.}

The $\ell_1$ approximation cannot be directly applied to~\ref{eq:Opt3}
since the $\ell_0$ norm constraints are coupled through~\eqref{eq:Con6} and
\eqref{eq:Con4-Opt3}.
%\eqref{eq:Con7-Opt3}.
Hence we present an iterative algorithm based on the $\ell_1$ reweighted heuristic,
where the weights depend on both the $\ell_1$ norm and the bandwidth of each segment $h_l$.
In each iteration of the algorithm, shown in Algorithm~\ref{alg:l1},
we solve~\ref{eq:Opt3}, but with the $\ell_0$ norm constraint~\eqref{eq:Con6}
replaced by the weighted sum
%~\eqref{eq:Con5-Opt6}.
%\begin{varsubequations}{P2}
%\label{eq:Opt6}
%\begin{align}
 % \maximize_{\br,\;\by,\;\bz} \;
%  &u(r_1,\cdots,r_k) \label{eq:Obj-Opt6}\\
 % \subjectto \;
 \vspace{-5pt}
\begin{equation}
\sum_{B\subset\NN_i}w^i_{B,l}y^i_{B,l}\leq1 .
 % \quad \forall i\in\NN,\;l=1,\cdots,k+1 \label{eq:Con5-Opt6}\\
 % & \eqref{eq:Con1}\text{--}%,\; \eqref{eq:Con2-Opt3},\; \eqref{eq:Con3-Opt3},\;
 %   \eqref{eq:Con6},\;
 %   \eqref{eq:Con8-Opt3}\text{--}%,\; \eqref{eq:Con9-Opt3},\; \eqref{eq:Con7-Opt3},\;
  %  \eqref{eq:Con4-Opt3}.
 % \nonumber
%\end{align}
%\end{varsubequations}%\noindent
\label{eq:l1approx}\vspace{-3pt}
\end{equation}
In Algorithm~\ref{alg:l1}, we refer to this reweighted $\ell_1$ version of \ref{eq:Opt3}
as P2.


%\EW{Can we move the output to the end of the algorithm?}

%The iterative algorithm for solving~\ref{eq:Opt6} is shown in Algorithm~\ref{alg:l1}.
Algorithm \ref{alg:l1} shows a random initialization of the weights to introduce
asymmetry in the first iteration. Otherwise, e.g., if all $w^i_{B,l}$'s are initialized
to a constant,
%,~\forall i\in\NN,B\subset\NN_i,l=1,\cdots,k+1$,
the solution stays symmetric over all segments,
because the constraints and objectives are identical for each segment.
A symmetric solution, i.e., $y^i_{B,1}=y^i_{B,2}=,\cdots,=y^i_{B,k},~\forall i\in \NN,~B\subset\NN_i$, \sout{is obviously not optimal for an arbitrary network topology}\EW{is generally not a solution to the original problem}. In each iteration, we solve P2 with the current weights $\bw$ to obtain the current $\bx,~\by,~\bz$ and $\bh=[h_1,\cdots,h_k]$. Then the weights are updated
%using the current solution by $w^i_{B,l}=\frac{1}{y^i_{B,l}+\alpha h_l}$.
as shown.
%\begin{align} \label{eq:wiBl}
%  w^i_{B,l}=\frac{1}{y^i_{B,l}+\alpha h_l}.
%\end{align}
The iterations terminate when the variables converge
or the maximum number of iterations $t_{max}$ is reached.
%The post processing in Algorithm~\ref{alg:l1} is to further optimize the performance and guarantee feasibility of the obtained spectrum allocation, which will be introduced in the following section.

The weight update in Algorithm~\ref{alg:l1} is obtained by approximating the $\ell_0$ norm
%Another way to understand this is to consider
%% $\frac{1}{y^i_{B,l}+\alpha h_l}$
%the right hand side of~\eqref{eq:wiBl} as the derivative of $\log(y^i_{B,l}+\alpha h_l)$,
with $\log(x+\epsilon)$ for small $\epsilon$ (see~\cite{candes2008enhancing} and the references therein).
Adding $\epsilon$ in the denominator allows small components to be propagated to the next iteration.
%As explained in~\cite{candes2008enhancing}, the rule of thumb for selecting $\epsilon$ is to make it close and a little less than the optimal nonzero point.
Here the goal is to obtain a {\em single} nonzero $y^i_{B,l}=h_l,~B\in\NN_i$, so we take
$\epsilon=\alpha h_l$. This is different from the fixed $\epsilon$
proposed in~\cite{candes2008enhancing} in that $\alpha h_l$ changes
with $h_l$ in each iteration.
That is, Algorithm~\ref{alg:l1} simultaneously searches
for the optimal reuse pattern and its associated bandwidth.

It is possible that Algorithm~\ref{alg:l1} produces multiple reuse patterns for a segment at termination.
%An exact consensus on a single reuse pattern may not be reached for each segment at the end of the iterations in Algorithm~\ref{alg:l1}.
In those cases AP $i$ chooses the dominant pattern
$B^i_l= \arg\max_{B\in\NN_i} y^i_{B,l}$. The global reuse pattern assigned to segment $l$
is then given by $A^*_l=\cup_{i:i\in B^i_l}\{i\}$. %Given the $k$ global reuse patterns, we can decide
The spectral efficiencies for each segment are determined by
the corresponding global reuse pattern (set of interfering APs):
$\bar{s}^{i\to j}_l=s^{i\to j}_{A^*_l\cap \mathcal{A}_j}$.
%\begin{align}
%\label{eq:SpeEff1}
%\bar{s}^{i\to j}_l=s^{i\to j}_{A^*_l\cap \mathcal{A}_j}
%\bar{s}^{i\to j}_l=
%%\begin{cases}
%%s^{i\to j}_{A^*_l\cap \mathcal{A}_j},~\text{if} A^*_l\cap \mathcal{A}_j\neq\emptyset,\\
%%0,~\text{Otherwise}.
%%\end{cases}
%\end{align}\noindent
Given the assignment of patterns to each spectrum segment,
the widths of the segments, $\bh=(h_l)_{l=1,\cdots,k+1}$ along with
the assignment of spectrum to mobiles, $\bar{\bx}= (\bar{x}^{i\to j}_l)_{j\in \KK,i\in\mathcal{A}_j,l=1,\cdots,k+1}$,
can then be re-optimized by solving the convex problem:
\begin{varsubequations}{P3}\vspace{-3pt}
\label{eq:Opt7}
\begin{align}
  \maximize_{\bar{\bx},\;\bh} \;
  & u(r_1,\cdots,r_k) \label{eq:Obj-Opt7}\\\vspace{-3pt}
   \subjectto \;
  & r_j = \sum_{l=1}^k\sum_{i\in \mathcal{A}_j} \bar{s}^{i\to j}_l \bar{x}_l^{i\to j},
  \quad \forall j\in\KK \label{eq:Con1-Opt7}\\
  &  \sum_{j=1}^k \bar{x}^{i\to j}_l\leq h_l,~\forall j\in \KK, i\in \mathcal{A}_j
  \label{eq:Con2-Opt7}
%  & \sum_{l=1}^k h_l =1, ~~~h_l, \bar{x}^{i\to j}_l \geq 0, ~~l=1,\cdots,k+1  \label{eq:Con3-Opt7}\\
% . \label{eq:Con4-Opt7}
\end{align}\vspace{-5pt}
\end{varsubequations}
$\bar{x}^{i\to j}_l \geq 0$ and \eqref{eq:Con4-Opt3} for $l=1,\cdots,k+1$.





%
%\subsection{Post processing}
%An exact consensus on a single reuse pattern may not be reached for each segment at the end of the iterations in Algorithm~\ref{alg:l1}. Instead, a dominating pattern is usually picked by each AP. Let the dominating pattern in AP $i$'s interference cluster on segment $l$ be $B^i_l= \arg\max_{B\in\NN_i} y^i_{B,l}$. The reuse patterns on segment $l$ can be unified by letting all APs using $A^*_l=\cup_{i:i\in B^i_l}\{i\}$. For the dominating pattern $B^i_l$, we usually have $y^i_{B^i_l,l}\approx h^l$ and $y^i_{B,l}\approx 0,\forall B\neq B^i_l$. If $i\in B^i_l$ then AP $i$ transmits on segment $l$. Otherwise it leaves segment $l$ blank for other APs. The global reuse pattern on segment $l$ includes all active APs on the segment. Given the $k$ global reuse patterns, we can decide the spectral efficiency for each link $(i\to j)$ by:
%\begin{align}
%\label{eq:SpeEff1}
%\bar{s}^{i\to j}_l=
%\begin{cases}
%s^{i\to j}_{A^*_l\cap \mathcal{A}_j},~\text{if} A^*_l\cap \mathcal{A}_j\neq\emptyset,\\
%0,~\text{Otherwise}.
%\end{cases}
%\end{align}
%
%The user association and spectrum allocation can be refined assuming each segment uses the unified pattern $A^*_l$. The corresponding optimization problem is~\ref{eq:Opt7}:
%\begin{varsubequations}{P3}
%\label{eq:Opt7}
%\begin{align}
%  \maximize_{\bar{\bx},\;\bh} \;
%  & u(r_1,\cdots,r_k) \label{eq:Obj-Opt7}\\
%  \subjectto \;
%  & r_j = \sum_{l=1}^k\sum_{i\in \mathcal{A}_j} \bar{s}^{i\to j}_l \bar{x}_l^{i\to j},
%  \quad \forall j\in\KK \label{eq:Con1-Opt7}\\
%  & \sum_{j=1}^k \bar{x}^{i\to j}_l\leq h_l,
%  \quad \forall l\in\KK\label{eq:Con2-Opt7}\\
%  & \bar{x}^{i\to j}_l,~h_l\geq0,
%  \quad \forall j\in \KK, i\in \mathcal{A}_j, l\in\KK\label{eq:Con3-Opt7}\\
%  & \sum_{l=1}^k h_l =1. \label{eq:Con4-Opt7}
%\end{align}%
%\end{varsubequations}%
%The number of variables and the number of constraints in~\ref{eq:Opt7} are $O(nk^2)$ and $O(k)$, respectively. Hence, it can be easily solved even with large $n$ and $k$. The spectrum is divided into $k$ pieces\footnote{This does not mean $k$ reuse patterns are used, since the selected pattern $A^*_l$ on different segments can be the same.} in the optimal solution to~\ref{eq:Opt7}. The bandwidth of the $l$th segment is $h_l$. On each segment $l$, AP $i$ will allocate $\bar{x}^{i\to j}_l$ amount of the spectrum to UE $j$.\footnote{If AP $i$ is inactive on segment $l$, then $\bar{x}^{i\to j}_l=0,\;\forall j\in\KK$.}
\begin{algorithm}
\caption{Iterative re-weighted $\ell_1$ approximation.}
\label{alg:l1}
\begin{algorithmic}[]
\INPUT {$(s^{i\to j}_C)_{j\in \KK, i\in \mathcal{A}_j, C\subset \mathcal{A}_j}$, and $(\lambda_j)_{j\in \Kcal}$.}
\OUTPUT{The widths of $k+1$ segments  $(h_l)_{l=1,\cdots,k+1}$, the $k+1$ active patterns $(A^*_l)_{l=1,\cdots,k+1}$, and the spectrum allocated to link $i\to j$ on segment $l$, $(\bar{x}^{i\to j}_l)_{j\in \KK,i\in\mathcal{A}_j,l=1,\cdots,k+1}$}
\Init{Choose $w_{B,l}^i$'s randomly in (0,1).\\
Set iteration counter $t=0$.}
\While{
%$(y^i_{B,l})_{i\in \NN,B\subset\NN_i,l\in \KK}$
Variables have not converged and $t<t_{max}$}
    \State 1. Solve P2
    %(\eqref{eq:Opt3} with \eqref{eq:l1approx}),
    with the current weights $\bw$.
    \State 2. Update $w^i_{B,l}=\frac{1}{y^i_{B,l}+\alpha h_l}$.
    \State 3. $t=t+1$.
\EndWhile
\Post {Determine the reuse patterns $\{A^*_l \}$
and spectral efficiencies $\{ \bar{s}^{i\to j}_l\}$ across segments;
solve P3.}
%Then optimize user allocations $(\bar{x}^{i\to j}_l)_{j\in \KK,i\in \mathcal{A}_j,l\in \KK}$.}
\end{algorithmic}
\end{algorithm}\vspace{-0.1cm}



\section{Numerical Results}
\label{sec:NumRes}
%\begin{table}%[!t]
%\caption{Parameter configurations.}
%\label{tab:par}
%\centering
%\begin{tabular}{||c|c||}
%\hline\hline
%Parameter & Value/Function\\
%\hline
%pathloss exponent & 3\\
%standard deviation of shadow fading & 3\\
%macro transmit PSD & 5 $\mu$W/Hz\\
%pico transmit PSD & 1 $\mu$W/Hz\\
%noise PSD & $1\times 10^{-7}$ $\mu$W/Hz\\
%total bandwidth & 20 MHz\\
%average packet length & 1 Mb\\
%\hline\hline
%\end{tabular}
%\end{table}
%The solutions obtained by Algorithm~\ref{alg:l1} are evaluated using numerical simulations.
%Throughout this section, we consider networks with $n$ APs and $k$ UEs.
%Unless specified particularly, the general assumption are given as follows: Among the $n$ APs,
The simulation results assume one macro AP is located at the center of the area,
and the remaining $n-1$ pico APs are randomly dropped around it.
%To generate the following simulation results,
%macro APs is placed at specific locations and pico APs are randomly deployed.
% of the network.
The $k$ UEs are placed on a rectangular lattice.
Link gains include both pathloss and shadowing.
Additional simulation parameters are shown
in the footnote.\footnote{The pathloss exponent is 3, standard deviation
of shadow fading is 3, macro-transmit PSD is 5 $\mu$W/Hz,
pico transmit PSD is 1 $\mu$W/Hz,
noise PSD is $10^{-7}$ $\mu$W/Hz,
total bandwidth is 20 MHz, and
average packet length is 1 Mb.}
%Table~\ref{tab:par}.
%In addition, the local cluster $\mathcal{A}_j$ of each UE $j$ is formed by picking the $4$ strongest APs around it.


\subsection{Small Network}
\label{subsec:small}
We first compare solutions to~\ref{eq:Opt} and~\ref{eq:Opt3} for a small network with $n=10$ and $k=32$. %Among the 10 APs, one macro AP is placed at the network center.
Since the number of variables is relatively small, we solve both versions of~\ref{eq:Opt} with and without the local neighborhood approximation using a standard convex optimization solver. The local neighborhoods are constructed by including the strongest four APs for each UE.
The solution to~\ref{eq:Opt3} is obtained using Algorithm~\ref{alg:l1}.
%
%To solve~\ref{eq:Opt3}, we can either use Algorithm~\ref{alg:l1} with a standard convex optimization solver or use the distributed approach based on Algorithm~\ref{alg:ADMM}.
%Two other simple schemes are also compared here. One is
We compare those with full spectrum reuse where each user is assigned
to the AP with the strongest signal (maxRSRP association), and also
the optimal orthogonal allocation,\footnote{Both spectrum allocation
and user association are optimized assuming each AP exclusively occupies a fraction of the spectrum.} i.e., only $\{y_{\{i\}}\}_{i\in\NN}$ are active.

Fig. ~\ref{fig:delay-N10K23} shows delay versus traffic arrival rate for all schemes.
The end of each curve represents the maximum arrival rate the scheme can support.
%The optimal orthogonal allocation (marked by round marker) quickly becomes saturated, as the orthogonal spectrum allocation is quite inefficient even with $n=10$. The full spectrum allocation with maxRSRP association (without any marker) achieves much higher throughput. However, the delay is also much higher as all APs interfere with each other.
The curves obtained by solving~\ref{eq:Opt} with and without local neighborhood approximation are very close, which indicates considering the four strongest interferers is enough in such a small network.
%The solution to~\ref{eq:Opt3} obtained by Algorithm~\ref{alg:l1} and Algorithm~\ref{alg:ADMM} are almost on top of each other, which proves the validity of the ADMM based solution. Hence, in the subsequent results, we will only show solution using Algorithm~\ref{alg:l1}.
The solution to~\ref{eq:Opt3} incurs slightly larger delay than the solution to~\ref{eq:Opt}.
%The throughput regions obtained through solving~\ref{eq:Opt} and~\ref{eq:Opt3} are the same. For all the simulations in this section, we limit the maximum number of iterations in Algorithm~\ref{alg:l1} to eight.
The jointly optimized spectrum allocations and user associations achieve substantial delay reduction as well as eight times throughput compared to full frequency reuse with maxRSRP association.

The optimized spectrum allocation and user associations obtained
%by solving~\ref{eq:Opt} without the local neighborhood approximation, and
solving~\ref{eq:Opt3} are depicted in
%Fig.~\ref{fig:Sol1} and
Fig.~\ref{fig:Sol3}.
%In Fig.~\ref{fig:TopoN10},
The macro and pico APs are represented by the bigger and smaller towers;
each handset represents a UE. Each solid line shows a connection between
the corresponding AP and UE. The grid for each UE shows
the spectrum used to serve that UE.
The traffic arrival rate (in (0, 100)) for each UE is shown under each grid.
%The allocation achieved by the solution to~\ref{eq:Opt} with local neighborhood approximation is omitted, since it is almost identical to the solution to~\ref{eq:Opt} without local neighborhood approximation.
The user association in Fig.~\ref{fig:Sol3} is close but not identical to that
obtained by solving~\ref{eq:Opt} (not shown)
%in Fig.~\ref{fig:Sol1}
due to the approximation in Algorithm~\ref{alg:l1},
%Take the UE on the up left corner as an example, more APs transmit to it using a larger portion of the spectrum in Fig~\ref{fig:Sol3}. This also
which explains the performance difference
%between the solution to~\ref{eq:Opt} and the solution to~\ref{eq:Opt3}
shown in Fig.~\ref{fig:delay-N10K23}.

%\begin{figure}%[!t]
%\begin{minipage}{\columnwidth}
%	\centering
%	\includegraphics[width=\columnwidth]{N10K23}
%	\caption{Delay versus traffic arrival rate curves for a HetNet with $n=10$ and $k=23$.}
%	\label{fig:delay-N10K23}
%\end{minipage}\hfill
%\begin{minipage}{\columnwidth}
%		\centering
%		\includegraphics[width=\columnwidth]{N100K200}
%		\caption{Delay versus traffic arrival rate curves for a HetNet with $n=100$ and $k=195$.}
%		\label{fig:delay-N100K195}
%\end{minipage}
%\end{figure}

\begin{figure}%[!t]
\centering
\includegraphics[width=0.9\columnwidth, height=65mm]{N10K23}\vspace{-0.1cm}
\caption{Delay versus traffic arrival rate curves for a small network with $n=10$ and $k=23$.}
\label{fig:delay-N10K23}
\end{figure}
%\vspace{-3cm}
\begin{figure} %*}
\centering \vspace{-0.2cm}
%\subfloat[The solution to~\ref{eq:Opt}]{
%\includegraphics[width=0.45\columnwidth]{TopoN10K23P1}
%\label{fig:Sol1}}
%\subfloat[The solution to~\ref{eq:Opt3}]{
\includegraphics[width=0.65\columnwidth, height=60mm]{TopoN10K23P3}\vspace{-0.1cm}
\caption{Example spectrum allocations and user associations from \eqref{eq:Opt3}.}
\setlength{\belowcaptionskip}{-10pt}\vspace{-0.5cm}
\label{fig:Sol3}
%\label{fig:TopoN10}
\end{figure} %*}

%\subsection{Performance in Medium-Size Networks}
%\label{subsec:med}
%We present the performance comparison of different algorithms in a larger network with $n=30$ APs and $k=46$ UEs, deployed on a $600\times600$ square meter area. Both~\ref{eq:Opt} and~\ref{eq:Opt3} become computationally prohibitive due to the $2^{30}$ variables for global patterns. Hence we compare the solution to the simple full-spectrum-reuse with maxRSRP association, the optimal orthogonal allocation and the optimal user association under full spectrum allocation. A simplified version of~\ref{eq:Opt3} is also compared. Instead of using 46 segments, 5 segments are used, which limit the solution to no more than five active patterns. The delay versus average traffic arrival rate curves obtained by the five different schemes are shown in Fig.~\ref{fig:delay-N30K46}. The optimal orthogonal allocation becomes more inefficient. As the number of AP increase, each AP gets a smaller fraction of the entire spectrum on average. The solutions to~\ref{eq:Opt3} still achieves 4 times network throughput and substantial delay reduction compared with full-spectrum-reuse with maxRSRP association. Interestingly, using 5 segments in~\ref{eq:Opt3} achieves almost the same performance as using 46 segments. This is because there are only seven active patterns in the solution to~\eqref{eq:Opt3} with 46 segments. Many segments use the same active pattern in the solution. Optimizing user association under full-spectrum-allocation also has superior performance over full-spectrum-allocation with maxRSRP association. However, it can only support half of the maximum arrival rate achieved by the proposed solution.
%
%\begin{figure}
%\centering
%\includegraphics[width=\columnwidth]{N30K46}
%\caption{Delay versus traffic arrival rate curves for a HetNet with $n=30$ and $k=46$.}
%\label{fig:delay-N30K46}
%\end{figure}
\subsection{Large Network}
Fig. \ref{fig:delay-N100K195} compares the performance of different allocation schemes
in a large network with $n=100$ APs and $k=200$ UEs.
%The network is deployed on a $1250\times 1250$ square meter area.
To facilitate the simulations, we reduce the size of each UE neighborhood from four to three,
i.e., each UE can only be served by the three strongest APs.
%The optimal orthogonal allocation cannot support even very light load.
We compare full-spectrum-reuse with maxRSRP association, full-spectrum-reuse with optimized associations and the solution to~\ref{eq:Opt3} with 50 segments.
%Optimizing the user association under full-spectrum-reuse
%has much lower delay and marginal throughput gain compared
%to full-spectrum-reuse with maxRSRP association.
The solution to~\ref{eq:Opt3} achieves 1.5 times the throughput
of the full-spectrum-reuse with optimized association.
The solutions to~\ref{eq:Opt3} use no more than 23 active patterns ($<50$ available segments).
%which indicates it is enough to use 50 segments to formulate~\ref{eq:Opt3}.
The throughput gain achieved by the solution to~\ref{eq:Opt3} in this large network
is smaller than that in Section~\ref{subsec:small}.
This is mainly because we only consider the three strongest interferers,
which compromises the benefits from interference management.

\begin{figure}
\centering
\includegraphics[width=0.96\columnwidth,height=60mm]{N100K200}\vspace{-0.2cm}
\caption{Delay versus traffic arrival rate curves for a large network with $n=100$ and $k=195$.}\vspace{-0.4cm}
\label{fig:delay-N100K195}
\end{figure}

\section{Conclusion}
\label{sec:Con}
An approach to joint allocation of spectrum with user association has been
presented which exploits the sparsity of the optimal solution.
%The optimization occurs over a slow timescale, corresponding to the rate at which
%traffic statistics change. The computational
%complexity scales as $O(nk)$ instead of $O(2^n)$.  hence
%The approach can be applied
%to large  networks, e.g., on the order of 100 APs with several hundred mobile locations.
The proposed algorithm has been observed to achieve near-optimal
performance for small to medium-size networks for which the optimal solution can be computed,
and provides substantial gains relative to full frequency reuse.
Although not considered here, the formulation can be extended
to accommodate spatial selectivity and different power levels.
The performance-complexity tradeoff for such extensions is left for future work.

%The paper introduces a new look of future 5G networks with densely deployed APs through global radio resource management. By jointly optimizing spectrum allocation and user association across all APs on a slow timescale, substantial performance improvement can be achieved. The complexity issue of such an approach becomes critical in a large network with increased interactions among APs. With the help of advanced optimization techniques and the exponential decay nature of a wireless link, we explored the possibility of global radio resource management in large networks.
%The proposed framework and fully scalable solution suggests a way for potential city-wide global radio resource management. Power control and load dependent interference are not considered in the current problem formulation, which are future research directions.


%\section*{Acknowledgment}
%This work was supported in part by a gift from Futurewei Technologies.  The authors thank Dr.~Weimin Xiao and Dr.~Jialing Liu for useful discussions.


%\appendix
%%\section{Appendix}
%% \label{sec:App}
%
%We first present an equivalent formulation of~\ref{eq:Opt} using the local patterns and local variables introduced in Section~\ref{subsec:local}:
%\begin{varsubequations}{P8} \label{eq:Opt1}
%\begin{align}
%  \maximize_{\br,\;\bx,\;\by\;\bz} \;
%  & u(r^1,\cdots,r^k) \\
%  \subjectto \;
%  & r^j =\sum_{i\in \NA_j}\sum_{B\subset \NN_i}s^{i\to j}_{B\cap\NA_j}z^{i\to j}_B,
%  \quad \forall j\in\KK\label{eq:Con2-Opt1}\\
%  & z^{i\to j}_B=\sum_{A\subset\NN,\;A\cap\NN_i=B}x^{i\to j}_A, \nonumber\\
%  & \qquad \forall i\in\NN,\;\forall j\in \mathcal{U}_i,\;\forall B\subset\NN_i\label{eq:Con1-Opt1}\\
%  & \sum_{j\in\KK}x^{i\to j}_A \le y_A,
%  \quad \forall A\subset\NN \label{eq:Con4-Opt1}\\
%  & \sum_{A\subset \NN}y_A=1, \\
%  & x^{i\to j}_A\geq0,
%  \quad \forall i\in \NN,~\forall j\in \KK,~\forall A\subset \NN. \label{eq:Con3-Opt1}
%\end{align}%
%\end{varsubequations}%
%It is easy to see~\ref{eq:Opt1} is equivalent to~\ref{eq:Opt} under the local neighborhood assumption, based on our discussion about~\eqref{eq:LocalGlobal} and~\eqref{eq:LocalRate}. To prove Theorem~\ref{thm:main}, we next prove~\ref{eq:Opt3} is equivalent to~\ref{eq:Opt1}.
%
%To show the equivalence between~\ref{eq:Opt1} and~\ref{eq:Opt3}, we introduce two auxiliary optimization problems. First, we apply the $k$ constituents transformation to~\ref{eq:Opt1} to obtain~\ref{eq:Opt4}:
%\begin{varsubequations}{P9}
%\label{eq:Opt4}
%\begin{align}
%  \maximize_{\br,\;\bx,\;\by} \;
%  & u(r^1,\cdots,r^k) \\
%  & r^j=\sum_{l=1}^k\sum_{i\in\mathcal{A}_j}\sum_{B\subset\NN_i} s^{i\to j}_{B\cap\mathcal{A}_j}\;z^{i\to j}_{B,l},
%  \quad \forall j\in \KK \label{eq:Con1-Opt4}\\
%  & z^{i\to j}_{B,l}=\sum_{A\cap\NN_i=B}x^{i\to j}_{A,l}, \nonumber\\
%  &\qquad \forall i\in\NN,\;\forall j\in\mathcal{U}_i,\;\forall B\subset\NN_i,\;\forall l\in\KK \label{eq:Con2-Opt4}\\
%  & \sum_{j\in\KK}x^{i\to j}_{A,l}=y_{A,l},
%  \quad \forall A\subset \NN,\;\forall l\in\KK\label{eq:Con3-Opt4}\\
%  & \sum_{A\subset \NN}y_{A,l}\leq h_l,
%  \quad \forall l\in\KK \label{eq:Con5-Opt4}\\
%  & \sum_{A\in\NN}|y_{A,l}|_0\leq 1,
%  \quad \forall l\in\KK\label{eq:Con6-Opt4}\\
%  & \sum_{l\in\KK}h_l=1 \label{eq:Con7-Opt4}\\
%  %&y_{A,l}\geq0,&\forall A\subset\NN,\; l\in\KK\label{eq:Con4-Opt4}\\
%  & x^{i\to j}_{A,l}\geq 0,
%  \quad \forall i\in\NN,\;\forall j\in\KK,\; \forall A\subset \NN,\;\forall l\in\KK \label{eq:Con8-Opt4}\\
%  & h_l\geq0,
%  \quad \forall l\in\KK. \label{eq:Con9-Opt4}
%\end{align}
%\end{varsubequations}\noindent
%\ref{eq:Opt4} can be considered as first reformulating~\ref{eq:Opt1} by having $k$ constituents of the variables in~\ref{eq:Opt1} for $k$ spectrum segments, where each segment has bandwidth $h_l$. Again, cardinality constraint~\eqref{eq:Con6-Opt4} ensures only one active pattern in each segment. Given the formulation, we have the following property:
%
%\begin{proposition}\label{prop:Opt1=Opt4}
%  \ref{eq:Opt1} and~\ref{eq:Opt4} are equivalent.
%\end{proposition}
%
%\begin{IEEEproof}
%%\label{prf:Opt1=Opt4}
%We first show that without the cardinality constraint~\eqref{eq:Con6-Opt4},~\ref{eq:Opt4} is equivalent to~\ref{eq:Opt1}. To be precise, without~\eqref{eq:Con6-Opt4}, if all variables with subscript $l\neq1$ is set to 0 except for $l=1$,~\ref{eq:Opt4} reduces to~\ref{eq:Opt1}. Thus~\ref{eq:Opt4} without~\eqref{eq:Con6-Opt4} is a relaxation of~\ref{eq:Opt1}. On the other hand, from any solution to~\ref{eq:Opt4}, we can combine the variables for segments $l=1,\cdots,k$ to one feasible solution of~\ref{eq:Opt1}. Hence, without~\eqref{eq:Con6-Opt4}, the equivalence holds.
%It remains to show with the $\ell_0$ constraint~\eqref{eq:Con6-Opt4}, the optimal solution to~\ref{eq:Opt1} is still in the feasible set of~\ref{eq:Opt4}. This is because the optimal solution to~\eqref{eq:Opt1} activates at most $k$ global patterns $y_A$, according to Theorem~\ref{thm:Spec}. If we take the optimal solution to~\ref{eq:Opt1} and let the $k$ active reuse patterns each correspond to a distinct index $l$ in~\ref{eq:Opt4}, we can find all constraints in~\ref{eq:Opt4} are satisfied.
%To be precise, suppose the $k$ active reuse patterns found by solving~\ref{eq:Opt1} are $A_1,\cdots,A_k$. Then the variables in~\ref{eq:Opt4} can be constructed based on the optimal solution to~\ref{eq:Opt1}, $\left(y^*_A\right)_{A\subset\NN}$ and $\left(x^{*i\to j}_A\right)_{i\in\NN, j\in\mathcal{U}_i, A\subset\NN}$ as follows:
%\begin{align}\label{eq:Opt1->Opt4}
%h_l&=y^*_{A_l}&\\\nonumber
%y_{A,l}&=\left\{
%\begin{array}{lr}
%y^*_{A_l}  & \text{if} A=A_l\\
%0          & \text{otherwise}
%\end{array}
%\right.\\\nonumber
%x^{i\to j}_{A,l}&=\left\{
%\begin{array}{lr}
%x^{*i\to j}_{A_l} & \text{if} A=A_l\\
%0                         & \text{otherwise}
%\end{array}
%\right.
%\end{align}
%for $l=1,\cdots,k$. Then it is easy to see that $r_j=r_j^*$ for every $j\in\KK$ and all constraints in~\ref{eq:Opt4} are satisfied.
%\end{IEEEproof}
%
%The second auxiliary optimization problem is formulated by adding the global variables back to~\ref{eq:Opt3} to obtain~\ref{eq:Opt5}:
%\begin{varsubequations}{P10} \label{eq:Opt5}
%\begin{align}
%  \maximize_{\br,\;\bz,\;\by} \;
%  & u(r^1,\cdots,r^k) \label{eq:Obj-Opt5}\\
%  \subjectto \;
%  & r^j=\sum_{l=1}^k\sum_{i\in\mathcal{A}_j}\sum_{B\subset\NN_i} s^{i\to j}_{B\cap\mathcal{A}_j}\;z^{i\to j}_{B,l},
%  \quad \forall j\in \KK \label{eq:Con1-Opt5}\\
%  & \sum_{j\in\mathcal{U}_i}z^{i\to j}_{B,l}=\sum_{A\subset\NN: A\cap\NN_i =B}y_{A,l},
%  \nonumber\\
%  & \qquad \forall i\in\NN,\;\forall B\subset\NN_i,\;\forall l\in \KK\label{eq:Con2-Opt5}\\
%  & \sum_{A\subset \NN}y_{A,l}\leq h_l,
%  \quad \forall l\in \KK \label{eq:Con5-Opt5}\\
%  & \sum_{A\subset\NN}|y_{A,l}|_0\leq1,
%  \quad \forall l\in \KK \label{eq:Con6-Opt5}\\
%  & y_{A,l}\geq0,
%  \quad \forall A\subset\NN\label{eq:Con7-Opt5}\\
%  & z^{i\to j}_{B,l}\geq 0,
%  \quad \forall i\in\NN\; \forall j\in\mathcal{U}_i\;\forall B\subset\NN_i,\;l\in\KK\label{eq:Con9-Opt5}\\
%  & \sum_{l\in \KK}h_l=1 \label{eq:Con8-Opt5}\\
%  & h_l\geq0,
%  \quad \forall l\in \KK. \label{eq:Con4-Opt5}
%\end{align}%
%\end{varsubequations}%
%We will show the following equivalence property:
%\begin{proposition}
%\label{prop:Opt3=Opt5}
%\ref{eq:Opt3} and~\ref{eq:Opt5} are equivalent.
%\end{proposition}
%
%\begin{IEEEproof}
%\label{prf:Opt3=Opt5}
%The difference between~\ref{eq:Opt3} and~\ref{eq:Opt5} are only in the $y$ variables, where $y^i_{B,l}$ in~\ref{eq:Opt3} is associated with the local reuse pattern $B\subset\NN_i$, whereas $y_{A,l}$ in~\ref{eq:Opt5} is associated with global reuse pattern $A\subset\NN$. The objectives of~\ref{eq:Opt3} and~\ref{eq:Opt5} are identical. The other constraints are equivalent except for constraints~\eqref{eq:Con2-Opt3}--\eqref{eq:Con8-Opt3} and~\eqref{eq:Con2-Opt5}--\eqref{eq:Con7-Opt5}. It remains to show that the two optimization problems share the same optimal solution. Let the $y$ variables given by these two sets of constraints be related by:
%\begin{align}
%\label{eq:Opt3->Opt5}
%y^i_{B,l}=\sum_{A\subset\NN:A\cap\NN_i=B}y_{A,l},\;\forall l\in\KK,\;\forall i\in\NN,\;\forall B\subset\NN_i
%\end{align}
%Hence~\eqref{eq:Con2-Opt3} becomes~\eqref{eq:Con2-Opt5}.
%Suppose all constraints in~\ref{eq:Opt5} hold. Then~\eqref{eq:Con8-Opt3} is trivial by~\eqref{eq:Opt3->Opt5}, and we also have:
%\begin{align}
%\label{eq:y_eq}
%\sum_{B\subset\NN_i:B\cap\NN_m=C}y^i_{B,l}&=\sum_{B\subset\NN_i:B\cap\NN_m=C}\sum_{A\subset\NN:A\cap\NN_i=B}y_{A,l}\\\nonumber
%&=\sum_{A\subset\NN:A\cap\NN_i\cap\NN_m=C}y_{A,l}
%\end{align}
%which implies~\eqref{eq:Con3-Opt3}. Moreover,~\eqref{eq:Opt3->Opt5} also implies
%\begin{align}
%\label{eq:sumOpt3Opt5}
%\sum_{B\subset\NN_i}y^i_{B,l}&=\sum_{B\subset\NN_i}\sum_{A\subset\NN:A\cap\NN_i=B}y_{A,l}\\\nonumber
%&=\sum_{A\subset\NN}y_{A,l}\leq h_l
%\end{align}
%Therefore~\eqref{eq:Con5-Opt3} is satisfied. Similarly,~\eqref{eq:Con6-Opt3} is established as:
%\begin{align}
%\label{eq:L0Opt3Opt5}
%\sum_{B\subset\NN_i}|y^i_{B,l}|_0 &=\sum_{B\subset\NN_i}\left|\sum_{A\subset\NN:A\cap\NN_i=B}y_{A,l}\right|_0\\\nonumber
%&\leq\sum_{B\subset\NN_i}\sum_{A\subset\NN:A\cap\NN_i=B}|y_{A,l}|_0\\\nonumber
%&\leq\sum_{A\subset\NN}|y_{A,l}|_0\\\nonumber
%&\leq1
%\end{align}
%Therefore, any solution in the feasible set of~\ref{eq:Opt5} will be a feasible solution of~\ref{eq:Opt3}, which also holds for the optimal solution to~\ref{eq:Opt5}.
%
%It remains to show the optimal solution to~\ref{eq:Opt3} also belongs to the feasible set of~\ref{eq:Opt5}. The key is to reconstruct global variables $y_{A,l}$ from the local variables $y^i_{B,l}$. Let us focus on a specific $l\in\KK$. Constraint~\eqref{eq:Con6-Opt3} dictates that there is at most one active local reuse pattern in every neighborhood. Namely, we can identify one active pattern $B^i_l\subset\NN_i$ for every $i\in\NN$, such that $y^i_{B,l}=0,~\forall B\neq B^i_l$. In fact all the active patterns will be assigned the same bandwidth for the same $l$. If $B^i_l$ and $B^m_l$ overlaps, then by~\eqref{eq:Con3-Opt3} we have:
%\begin{align}
%\label{eq:Opt5->Opt3}
%y^i_{B^i_l,l}&=\sum_{B\subset\NN_i:B\cap\NN_m=B^i_l\cap B^m_l}y^i_{B,l}\\\nonumber
%&=\sum_{B\subset\NN_m:B\cap\NN_i=B^i_l\cap B^m_l}y^m_{B,l}\\\nonumber
%&=y^m_{B^m_l,l}
%\end{align}
%For those active patterns not overlapping, i.e., $B^i_l\cap B^m_l=\emptyset$, $y^i_{B^i_l,l}$ and $y^m_{B^m_l,l}$ are not necessarily equal. However, due to the total bandwidth constraint on segment $l$, all active pattern in the optimal solution should have:
%\begin{align}
%\label{eq:BwOpt3Opt5}
%y^{*i}_{B^i_l,l}=h_l, \forall i\in \NN.
%\end{align}
%This is because if any $y^{*i}_{B^i_l,l}<h_l$ the objective can always be reduced by increasing $y^{*i}_{B^i_l,l}$ and those active patterns overlapping with it. Hence, we can determine a global reuse pattern for segment $l$ as:
%\begin{align}
%\label{eq:GlobalPattern}
%A_l = \cup_{i\in\NN}B^i_l
%\end{align}
%The global variables are thus given by:
%\begin{align}
%\label{eq:GlobalY}
%y^{*}_{A,l}=\left\{
%\begin{array}{lr}
%h_l &\text{if} A=A_l\\
%0 &\text{otherwise}
%\end{array}
%\right.
%\end{align}
%It is easy to verify $\big(y^{*}_{A,l}\big)_{A\subset\NN,l\in\KK}$ satisfies~\eqref{eq:Con2-Opt5}--\eqref{eq:Con7-Opt5}. Hence, the optimal solution to~\ref{eq:Opt3} is in the feasible set of~\ref{eq:Opt5}. The equivalence is therefore established.
%\end{IEEEproof}
%
%The last step is to prove the following equivalence:
%\begin{proposition}
%\label{prop:Opt4=Opt5}
%\ref{eq:Opt4} is equivalent to~\ref{eq:Opt5}.
%\end{proposition}
%
%\begin{IEEEproof}
%\label{prf:Opt4=Opt5}
%The objective functions of~\ref{eq:Opt4} and~\ref{eq:Opt5} are identical. Except for the constraints~\eqref{eq:Con2-Opt4}, ~\eqref{eq:Con3-Opt4} and the constraint~\eqref{eq:Con2-Opt5}, the reset of the constraints in the two problems are also identical. It is easy to see $x^{i\to j}_{A,l}$ is auxiliary variable in~\ref{eq:Opt4}. Hence,~\eqref{eq:Opt4} and~\eqref{eq:Opt5} are equivalent.
%\end{IEEEproof}
%
%Proposition~\ref{prop:Opt1=Opt4}--\ref{prop:Opt4=Opt5} immediately implies~\ref{eq:Opt3} is equivalent to~\ref{eq:Opt1}, and therefore Theorem~\ref{thm:main}.

\bibliographystyle{ieeetr}
\begin{thebibliography}{1}

\bibitem{yu2013multicell}
W.~Yu, T.~Kwon, and C.~Shin, ``Multicell coordination via joint scheduling,
  beamforming, and power spectrum adaptation,'' {\em {IEEE} Trans. Wireless
  Commun.}, vol.~12, no.~7, pp.~1--14, 2013.

\bibitem{wang2013joint}
F.~Wang, L.~Song, Z.~Han, Q.~Zhao, and X.~Wang, ``Joint scheduling and resource
  allocation for device-to-device underlay communication,'' in {\em Proc.\
  Conf.\ Wireless \ Comm. and Networking}, pp.~134--139, IEEE, 2013.

\bibitem{HonLuo13TransSP}
M.~Hong and Z.-Q. Luo, ``Distributed linear precoder optimization and base
  station selection for an uplink heterogeneous network,'' {\em {IEEE} Trans.
  Signal Process.}, vol.~61, pp.~3214--3228, June 2013.

\bibitem{zhuang2014traffic}
B.~Zhuang, D.~Guo, and M.~L. Honig, ``Traffic-driven spectrum allocation in
  heterogeneous networks,'' {\em {IEEE} J. Sel. Areas Commun.}, vol.~PP,
  no.~99, pp.~1--1, 2015.

\bibitem{kuang2014joint}
Q.~Kuang, ``Joint user association and reuse pattern selection in heterogeneous
  networks,'' in {\em 2014 11th International Symposium on Wireless
  Communications Systems (ISWCS)}, pp.~401--405, IEEE, 2014.

\bibitem{zhuang2016scalable}
B.~Zhuang, D.~Guo, E.~Wei, and M.~L. Honig, ``Scalable spectrum allocation and
  user association in networks with many small cells,''
 {\em https://arxiv.org/abs/1701.03247,} submitted to {\em IEEE Trans. Commun.}.

\bibitem{kuang2016energy}
Q.~Kuang and W.~Utschick, ``Energy management in heterogeneous networks with
  cell activation, user association, and interference coordination,'' {\em IEEE
  Trans. Wireless Commun.}, vol.~15, June 2016.

\bibitem{candes2008enhancing}
E.~J. Candes, M.~B. Wakin, and S.~P. Boyd, ``Enhancing sparsity by reweighted
  $l_1$ minimization,'' {\em Journal of Fourier analysis and applications},
  vol.~14, no.~5-6, pp.~877--905, 2008.

\bibitem{zhuang2015energy}
B.~Zhuang, D.~Guo, and M.~L. Honig, ``Energy-efficient cell activation, user
  association, and spectrum allocation in heterogeneous networks,'' {\em IEEE
  J. Sel. Areas Commun.}, vol.~PP, no.~99, pp.~1--1,
  2016.

\end{thebibliography}

%\bibliography{IEEEabrv,def,Thesis}
\end{document}

1702.05679.extracted_bib
{yu2013multicell W.~Yu, T.~Kwon, and C.~Shin, ``Multicell coordination via joint scheduling, beamforming, and power spectrum adaptation,'' {\em {IEEE Trans. Wireless Commun., vol.~12, no.~7, pp.~1--14, 2013.
1702.05679.extracted_bib
{wang2013joint F.~Wang, L.~Song, Z.~Han, Q.~Zhao, and X.~Wang, ``Joint scheduling and resource allocation for device-to-device underlay communication,'' in {\em Proc.\ Conf.\ Wireless \ Comm. and Networking, pp.~134--139, IEEE, 2013.
1702.05679.extracted_bib
{HonLuo13TransSP M.~Hong and Z.-Q. Luo, ``Distributed linear precoder optimization and base station selection for an uplink heterogeneous network,'' {\em {IEEE Trans. Signal Process., vol.~61, pp.~3214--3228, June 2013.
1702.05679.extracted_bib
{zhuang2014traffic B.~Zhuang, D.~Guo, and M.~L. Honig, ``Traffic-driven spectrum allocation in heterogeneous networks,'' {\em {IEEE J. Sel. Areas Commun., vol.~PP, no.~99, pp.~1--1, 2015.
1702.05679.extracted_bib
{kuang2014joint Q.~Kuang, ``Joint user association and reuse pattern selection in heterogeneous networks,'' in {\em 2014 11th International Symposium on Wireless Communications Systems (ISWCS), pp.~401--405, IEEE, 2014.
1702.05679.extracted_bib
{zhuang2016scalable B.~Zhuang, D.~Guo, E.~Wei, and M.~L. Honig, ``Scalable spectrum allocation and user association in networks with many small cells,'' {\em https://arxiv.org/abs/1701.03247, submitted to {\em IEEE Trans. Commun..
1702.05679.extracted_bib
{kuang2016energy Q.~Kuang and W.~Utschick, ``Energy management in heterogeneous networks with cell activation, user association, and interference coordination,'' {\em IEEE Trans. Wireless Commun., vol.~15, June 2016.
1702.05679.extracted_bib
{candes2008enhancing E.~J. Candes, M.~B. Wakin, and S.~P. Boyd, ``Enhancing sparsity by reweighted $l_1$ minimization,'' {\em Journal of Fourier analysis and applications, vol.~14, no.~5-6, pp.~877--905, 2008.
1702.05679.extracted_bib
{zhuang2015energy B.~Zhuang, D.~Guo, and M.~L. Honig, ``Energy-efficient cell activation, user association, and spectrum allocation in heterogeneous networks,'' {\em IEEE J. Sel. Areas Commun., vol.~PP, no.~99, pp.~1--1, 2016.