\section{Introduction}

Massive machine-type  communications (mMTC) is a class of services planned for 5G and 6G systems to provide access for enormous numbers of connected devices, possibly on the order of tens of billions~\cite{Saad2020,Carvalho2017,Bockelmann2016}. Massive multiple-input multiple-output (mMIMO) technology has been proposed as a way to improve the spectral and energy efficiency of such systems by exploiting the spatial multiplexing provided by the large numbers of antennas~\cite{Larsson2014,Bjornson2016}. To take advantage of the boost in spectral and energy efficiency provided by large-scale antenna systems, channel state information (CSI) is required at the base station (BS). Uplink CSI is obtained through the transmission of pilot sequences from the user equipments (UEs) to the BS, which estimates the channel using the received signal by employing any estimation technique such as least squares or minimum mean square error. A time-division duplexing (TDD) protocol is likely to be implemented in mMIMO networks in order to minimize signaling overhead and to avoid downlink channel estimation by exploiting the reciprocity between the uplink and downlink channels~\cite{Al-hubaishi2019}.

One of the main challenges of mMTC is the design of access protocols that support the huge number of MTC devices. Due to the large number of UEs, resources have to be share
, which precludes the allocation of orthogonal pilot sequences to all UEs. \textit{Pilot reuse} can be used to overcome the lack of resources, at the cost of leading to interference between UEs sharing the same pilot sequence, also known as \textit{pilot contamination}~\cite{Bjornson2017}.
While legacy systems rely on code and frequency multiplexing, mMTC can exploit spatial multiplexing through mMIMO to avoid/alleviate pilot contamination~\cite{Carvalho2017,Sanguinetti2020}.

Different approaches to mitigate pilot contamination have been proposed in the literature. In~\cite{Akbar2018}, the authors proposed a location-aware pilot assignment for multi-cell systems with Rician fading channels which requires knowledge about the large-scale fading, angle of arrival (AoA), and the channel Rician factor.
In~\cite{Al-hubaishi2019}, a pilot assignment scheme was devised to deal with pilot interference in multi-cell mMIMO based on the UEs' large-scale coefficients, which are assumed to be known in all cells. Their objective is to assign orthogonal pilot sequences to UEs with poor channel quality.
You {\it et al.}~\cite{You2015} proposed a pilot reuse strategy based on the channel covariance matrices for single-cell mMIMO systems with correlated channels. They proposed a pilot allocation algorithm that assigns orthogonal pilot sequences to UEs with similar channel covariance matrices, measured by the covariance matrix distance (CMD) metric.

Recently, we proposed a pilot reuse scheme in \cite{Ribeiro2020} that utilizes \textit{channel charting} (CC) to exploit the spatial information existing in CSI, aiming to maximize the AoA distances between the UEs using the same pilot sequence. CC is a framework proposed in~\cite{Studer2018} to estimate the relative position of devices in an unsupervised manner which maps the information obtained from the measured long-term CSI at the BS into a low-dimensional chart, in which the relative positions of UEs are preserved.
Our results~\cite{Ribeiro2020} showed that CSI can be utilized to allocate the pilot sequences, and improve the channel estimation accuracy and symbol error rate (SER), for single-cell correlated channels systems.

In this paper, we propose pilot reuse strategies to cope with the pilot contamination in multi-sector single-cell mMIMO networks under spatially correlated channels. To exploit the degree of orthogonality between the UEs' channels, granted by the spatial correlation, the pilot reuse is performed based on two features encompassing the channel correlation, especially in the angular domain: 1) the CMD metric and 2) the CC created based on the CMD metric. Utilizing the proposed features along with the pilot allocation algorithm devised in~\cite{Ribeiro2020}, the numerical results show that the angular information can be exploited also for multi-sector systems with spatially correlated channels, resulting in enhanced achievable rate and SER compared to baselines.

\section{System model}

\begin{figure}[t]
    \centering

    \includegraphics[scale=0.3]{Channel-ring-modelV6.pdf}
    \caption{Uplink mMTC scenario with $K$ active (green) and $N-K$ inactive users (gray). The BS is equipped with $S=3$ $M$-element ULAs.}\vspace{-3mm}
    \label{fig:channel}
\end{figure}

We consider an uplink communication scheme with a set ${\mathcal{N}=\{1,\ldots,N\}}$ of $N$ single-antenna UEs uniformly distributed within a cell, from which only $K<N$ are randomly active at any given time. In order to focus on the development of pilot assignment strategies to improve the channel estimation accuracy, and, subsequently, enhance the system performance in mMTC, we assume that the BS knows the set of active UEs at each transmission instant\footnote{This assumption can be invoked by the fact that there exists a multitude of compressed sensing based approaches, e.g., \cite{chen2018sparse,Senel2018,Djelouat2020Joint}, to solve the UE activity detection problem. However, this is outside of the scope of this paper and is left for future work.}.

The cell is divided into $S$ sectors\footnote{In practical scenarios, the coverage area is usually divided into three or six sectors of $120^\circ$ or $60^\circ$~\cite{Morales2019}. Therefore, we consider $S=3$ in the examples.} and the
transmitted signal is received
through $S$ uniform linear arrays (ULAs), each having $M$ antenna elements, as shown in Fig.~\ref{fig:channel}. The array response vector for a ULA is given by
\begin{equation}
\vec{a}_{\mathrm{r}}(\theta) = \left[{1, e^{-j2\pi \Delta_{\mathrm{r}}\cos(\theta)}, \ldots, e^{-j2\pi (M-1)\Delta_{\mathrm{r}} \cos(\theta)}}\right]^{\mbox{\scriptsize T}}
\end{equation}
where $\Delta_{\mathrm{r}}$ is the normalized spacing between the antenna elements in units of wavelengths and $\theta$ is the AoA, i.e., the incident angle of the received signal on the antenna array\cite[Sec. 7.2.1]{TseBook}.

We adopt the one-ring channel model~\cite[p. 236]{Bjornson2017} which assumes that the multi-path components are concentrated around the UEs while the BS is located in an elevated position lacking scatterers close to it. Thus, the uplink channel vector for user $n\in\mathcal{N}$ associated with the ULA of sector $s\in \mathcal{S}=\{1,\ldots,S\}$ is modelled as a superposition of $L$ propagation paths as
\begin{equation}
\label{eq::hk-sector}
\vec{h}_{n,s} = \frac{1}{\sqrt{L}}\sum_{l=1}^{L}\sqrt{\beta_{n,s,l}}\alpha_{n,s,l}\vec{a}_{\mathrm{r}}(\theta_{n,s,l}),
\end{equation}
where $\alpha_{n,s,l}$ is the complex gain of the $l$th path assumed to be an independent and identically distributed (i.i.d.) complex Gaussian random variable with zero mean and $\mathbb{E}\{|\alpha_{n,s,l}|^2\}=1$. In \eqref{eq::hk-sector}, the large-scale propagation effects and the BS antenna gain in the channel are captured in $\beta_{n,s,l} \in \mathbb{R}$, which follows the free space path loss model as described in~\cite[Eq.~(2.7)]{Goldsmith2005}, which is defined as
 \begin{equation}
 \beta_{n,s,l} \mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}} 10^{\frac{G_{\mathrm{A}}(\theta_{n,s,l})}{10}}\left(\frac{\lambda}{4\pi d_{n}}\right)^2,
 \label{eq::beta}
 \end{equation}
where $\lambda$ is the wavelength and $d_{n}$ is the distance between UE $n$ and the BS. In~\eqref{eq::beta}, $G_{\mathrm{A}}(\theta_{n,s,l})\in \mathbb{R}$ represents the antenna gain through the $l$th path for user $n$ at the ULA $s$, and is given in dB by~\cite[Table 7.1-1]{TR36873,Morales2019}
\begin{equation}
    G_{\mathrm{A}}(\theta_{n,s,l}) = G_{\mathrm{A_{max}}} {-\textrm{min}\left[12\left(\frac{\theta_{n,s,l}}{\theta_{3\textrm{dB}}}\right)^2,\;A_\textrm{max}\right]},
    \label{eq::AntennaGain}
\end{equation}
where $G_{\mathrm{A_{max}}}$ is the maximum antenna gain, $\theta_{3\textrm{dB}}$ is the half power beamwidth, $A_\textrm{max}$ is the maximum attenuation in dB at the ULAs, and $\theta_{n,s,l}$ is the AoA of the $l$th path which is modelled as an i.i.d. random variable with uniform distribution $U(\theta_{n,s}^{\mathrm{min}},\theta_{n,s}^{\mathrm{max}})$, with $\theta_{n,s}^{\mathrm{min}}=\Bar{\theta}_{n,s}-\sqrt{3}\sigma_{\theta}$ and $\theta_{n,s}^{\mathrm{max}}=\Bar{\theta}_{n,s}+\sqrt{3}\sigma_{\theta}$. Here, $\Bar{\theta}_{n,s}\in [0,2\pi]$ is the incident angle between user $n$ and the ULA of sector $s$, $\sigma_{\theta}$ is the angular standard deviation, which specifies the AoA interval $\mathcal{A}_{n,s}=[\theta_{n,s}^{\mathrm{min}},\theta_{n,s}^{\mathrm{max}}]$ for the incoming multi-path components arriving from user $n$ at the ULA of sector $s$.

\section{Pilot Transmission and Channel Estimation}
To exploit the symmetry between downlink and uplink channels, we consider the TDD protocol, as depicted in Fig.~\ref{fig:TDD}. We assume that during one coherence block, the channels are time-invariant and flat-fading. At each coherence block, the set of active users transmit $\tau$ known symbols, of equal power $p_{\mathrm{u}}$, to the BS for channel estimation. Right after transmitting the pilot symbols, the $K$ active UEs transmit their data to the BS. Then, the remaining time, within the coherence block, is used for downlink communication.
\begin{figure}
    \centering
    \includegraphics[scale=0.45]{TDD.pdf}
    \caption{The channel response is time-invariant and frequency-flat within each coherence block. At the beginning of each coherence block, $\tau$ symbols are sent to estimate the channel.}\vspace{-2mm}
    \label{fig:TDD}
\end{figure}

We consider that the UEs are assigned with pilot sequences taken from a pool of $\tau$ \textit{orthogonal} sequences. However, due to a vast number of UEs in mMIMO networks, we consider that $\tau\ll N$, i.e., global \textit{pilot reuse} is employed in the cell area. Consequently, in the channel estimation phase, the same pilot is shared by $K/\tau$ UEs on the average. Let $\mathcal{K}=\{1,\ldots,K\}\subseteq\mathcal{N}$ represent the set of active UEs at a given transmission interval. Let ${\mathcal{T}=\{1,\ldots,\tau\}}$ be the set of indices of available pilot sequences. User $k\in\mathcal{K}$ transmits a pilot signal ${\vec{\boldsymbol{\psi}}_k=\sqrt{p_{\mathrm{u}}}\boldsymbol{\phi}_{\pi_k}}$,
where $\pi_k\in\mathcal{T}$ is the index of the pilot sequence assigned to UE $k$ and $\boldsymbol{\phi}_{\pi_k}\,\in\,\mathbb{C}^{\tau}$ is the corresponding pilot sequence from the orthogonal pilot book $\vec{\Phi}=[\vec{\boldsymbol{\phi}}_1,\ldots,\vec{\boldsymbol{\phi}}_{\tau}] \in \mathbb{C}^{\tau\times \tau}$. We define the set of UEs sharing the same pilot sequence as UE $k$, including UE $k$ itself, as ${\mathcal{G}_k=\{j\mid j\in\mathcal{N}, \, \pi_j=\pi_k\}}$.

\iffalse
In the channel estimation phase, user $k\green{\in\mathcal{K}}$ transmits a pilot signal ${\vec{\boldsymbol{\psi}}_k=\sqrt{p_{\mathrm{u}}}\boldsymbol{\phi}_{\pi_k}}$,
where $\pi_k$ is the index of the pilot sequence assigned to UE $k$ and $\boldsymbol{\phi}_{\pi_k}\,\in\,\mathbb{C}^{\tau}$ is the corresponding pilot sequence from the orthogonal pilot book $\vec{\Phi}=[\vec{\boldsymbol{\phi}}_1,\ldots,\vec{\boldsymbol{\phi}}_{\tau}] \in \mathbb{C}^{\tau\times \tau}$. Let ${\mathcal{N}=\{1,\ldots,N\}}$ represent the set of all UEs connected to the BS, $\mathcal{K}=\{1,\ldots,K\}\green{\subseteq\mathcal{N}}$ the set of active UEs at a given transmission interval, and ${\mathcal{T}=\{1,\ldots,\tau\}}$ the set of indices of available pilot sequences. Since a large number of UEs are expected in mMIMO networks, it is assumed that $\tau\ll N$. Therefore, we define the set of UEs sharing the same pilot sequence as UE $k$, including UE $k$ itself, as
${\mathcal{G}_k=\{j\mid j\in\mathcal{N}, \, \pi_j=\pi_k
\}}$.
\fi

The received signal for channel estimation at the ULA of sector $s$, $\vec{Y}_{s}=\left[\vec{y}_{s}^1,\ldots,\vec{y}_{s}^\tau\right]\in\mathbb{C}^{M \times \tau}$, can be written as
\begin{equation}
\vec{Y}_{s} = \vec{H}_{s}\vec{\Psi} + \vec{N}_{s},
\label{eq::receivedSignal}
\end{equation}
where $\vec{H}_{s}=\left[\vec{h}_{1,s},\ldots,\vec{h}_{K,s}\right]\in\mathbb{C}^{M \times K}$ is the channel matrix for the active UEs, $\vec{\Psi}=\left[\vec{\boldsymbol{\psi}}_1,\ldots,\vec{\boldsymbol{\psi}}_K\right]^{\mbox{\scriptsize T}}\in\mathbb{C}^{K\times \tau}$ is the pilot signal matrix, and $\vec{N}_{s}=\left[\vec{n}_{s}^{1},\ldots,\vec{n}_{s}^\tau\right]\in\mathbb{C}^{M \times \tau}$ is the noise matrix. We model the noise as an i.i.d. complex Gaussian random variable ${\vec{n}_{s} \sim \mathcal{CN}(0,\sigma_\mathrm{n}^2)}$, where $\sigma_\mathrm{n}^2$ is the noise power at each antenna element.

We consider that the linear minimum mean square error (LMMSE) receiver is deployed at the BS to jointly estimate the active UEs' channel vectors across the $S$ sectors. Therefore, the compound channel $\vec{h}_{k}\in\mathbb{C}^{MS}$
between UE $k\in\mathcal{K}$ and the BS is given by
 \begin{equation}
\label{eq::hk}
\vec{h}_{k} =
\begin{bmatrix}
\vec{h}_{k,1} \\
\vdots\\
\vec{h}_{k,S}
\end{bmatrix},
\end{equation}
with covariance matrix $\vec{R}_{k}=\mathbb{E}[\vec{h}_{k}\vec{h}_{k}^{\mbox{\scriptsize H}}]\in\mathbb{C}^{MS\times MS}$. The received pilot signal,
$
 \vec{Y}=\left[\vec{Y}_1^{\mbox{\scriptsize T}},\ldots,\vec{Y}_S^{\mbox{\scriptsize T}}\right]^{\mbox{\scriptsize T}}\in\mathbb{C}^{MS\times\tau}.
$

We assume that the channel covariance matrices of all UEs (active and inactive ones) are known at the BS. In order to be able to
retrieve the UEs' individual signals from the received compound signal, the BS needs to keep track of $\vec{R}_n,\,\forall n \in \mathcal{N}$. In practice, an initial training phase is required to obtain the first estimate for the covariance matrices. After this initial phase, the BS can keep updating each $\vec{R}_{n}$ based on the estimated channel.

The LMMSE estimate of the communications channel between user $k$ and the BS, $\vec{h}_{k}$ in \eqref{eq::hk}, is given as~\cite{You2015},
\begin{equation}
\label{eq::hk_hat}
    \vec{\hat{h}}_k = \vec{R}_{k}\vec{Q}_{k}^{-1}\vec{y}_{k}^{\mathrm{p}}.
\end{equation}
Here, $\vec{y}_{k}^{\mathrm{p}}$ represents the processed received signal for UE $k$ after correlating the received signal with the pilot sequence assigned to user $k$, i.e.,
\begin{equation}
\begin{split}
    \vec{y}_{k}^{\mathrm{p}} &= \frac{1}{p_{\mathrm{u}} \tau}\vec{Y}\vec{\boldsymbol{\psi}}_{k}^{*} \\
    &=\vec{h}_k+\underbrace{\sum_{j\,\in\, \mathcal{I}_k}{\vec{h}_j}}_\text{Pilot Interference}+\frac{1}{p_{\mathrm{u}}\tau}\vec{N}\vec{\boldsymbol{\psi}}_{k}^{*} ,
\end{split}
\label{eq::LS}
\end{equation}
where $\mathcal{I}_k=\{j|j\neq k,\,j\in\mathcal{G}_k\cap\mathcal{K}\}$ is the set of interfering users to user $k$,
and $\vec{Q}_{k}\in\mathbb{C}^{MS\times MS}$ in~\eqref{eq::hk_hat} is the covariance matrix of the received signal, given as
\begin{equation}
    \vec{Q}_{k}=\vec{R}_k+\sum_{j\in \mathcal{I}_k}\vec{R}_j + \frac{\sigma_\mathrm{n}^2}{p_u\tau}\vec{I}_{MS}.
    \label{eq::Rerror}
\end{equation}

Due to the orthogonality principle of the MMSE estimator~\cite[Sec. 12.4]{Kay1993}, the channel estimation error, $\vec{\tilde{h}}_{k}\sim \mathcal{CN}(\vec{0},\vec{R}_{\tilde{\vec{h}}_{k}})$, is independent of $\vec{\hat{h}}_{k}$. Therefore, we can decompose the channel $\vec{h}_{k}$ as $\vec{h}_{k}=\vec{\hat{h}}_{k}+\vec{\tilde{h}}_{k}$. Thus, the error covariance matrix for user $k$ is~\cite{Studer2018}
\begin{equation}
    \vec{R}_{\tilde{\vec{h}}_{k}}=\vec{R}_{k}-\vec{R}_{k}\vec{Q}_{k}^{-1}\vec{R}_{k}.
    \label{eq::err-covariance}
\end{equation}

\section{Uplink Data Transmission}

Let $\vec{x}=\left[x_1,\ldots,x_K\right]^{\mbox{\scriptsize T}}\in\mathbb{C}^{K}$ be the transmitted symbol vector at a given time instant. The corresponding received signal at the BS, $\vec{y}_{}=[\vec{y}_{1}^{\mbox{\scriptsize T}},\ldots,\vec{y}_{s}^{\mbox{\scriptsize T}}]^{\mbox{\scriptsize T}} \in \mathbb{C}^{MS}$, is given by
\begin{equation}
\vec{y}_{} = \vec{H}_{}\vec{x} + \vec{n}_{},
\end{equation}
where $\vec{H}_{}=\left[\vec{h}_{1},\ldots,\vec{h}_{K}\right]\in\mathbb{C}^{MS\times K}$ is the channel matrix, and $\vec{n}_{}=[\vec{n}_{1}^{\mbox{\scriptsize T}},\ldots,\vec{n}_{s}^{\mbox{\scriptsize T}}]^{\mbox{\scriptsize T}} \in \mathbb{C}^{MS}$ is the noise vector at the $M\times S$ receiver antenna elements. We assume that all UEs transmit with the same power, i.e., $p_{\mathrm{u}}=|x_k|^2$ is the transmit symbol power for user $k$.

Given the estimated channel $\vec{\hat{H}}=[\vec{\hat{h}}_{1},\ldots,\vec{\hat{h}}_{K}]\in\mathbb{C}^{MS\times K}$, we use the LMMSE receiver, $\vec{w}_k\in\mathbb{C}^{MS}$, derived in~\cite{You2015},
\begin{equation}
    \vec{w}_k=\left(\vec{\hat{H}}\vec{\hat{H}}^{\mbox{\scriptsize H}}+\sum_{k=1}^{K}{\vec{R}_{\tilde{\vec{h}}_k}}+\frac{\sigma_\mathrm{n}^2}{p_u}\vec{I}_{MS}\right)^{-1}\vec{\hat{h}_k},
    \label{eq::LMMSECombinervector}
\end{equation}
which takes into account the error covariance matrix $\vec{R}_{\tilde{\vec{h}}_k}$.

Therefore, the received symbol vector $\hat{\vec{r}}\in\mathbb{C}^{MS}$, after employing the LMMSE receiver $\vec{W}\in\mathbb{C}^{MS\times K}$, is given by
\begin{equation}
    \hat{\vec{r}}=\vec{W}^{\mbox{\scriptsize H}} \vec{y}.
\end{equation}

Given the communication model depicted in Fig.~\ref{fig:TDD}, the uplink spectral efficiency is lower bounded by~\cite[Th. 4.1]{Bjornson2017}
\begin{equation}\label{eq::net-se}
    R^{\mathrm{up}}_k=\left(1-\frac{\tau}{T_c}\right)R^{\mathrm{ach,up}}_k,
\end{equation}
where $T_c$ is the coherence time and
the corresponding achievable uplink rate for user $k$, $R^{\mathrm{ach,up}}_k$, is expressed as
\begin{equation}
    R^{\mathrm{ach,up}}_k=\mathbb{E}\left\{\log_2{\left(1+\gamma_k^{\textrm{up}}\right)}\right\},
\label{eq::Rachievable}
\end{equation}
where the expectation is taken over the channel realizations,
and $\gamma_k^{\textrm{up}}$ is the \textit{instantaneous} uplink signal-to-interference-and-noise ratio (SINR), given as
\begin{equation}
    \gamma_k^{\textrm{up}}=\frac{|\vec{w}_k^{\mbox{\scriptsize H}}\vec{\hat{h}}_k|^2}{\vec{w}_k^{\mbox{\scriptsize H}}\left(\sum_{j\neq k}{\vec{\hat{h}}_j\vec{\hat{h}}_j^{\mbox{\scriptsize H}}+\sum_{n=1}^{K}{\vec{R}_{\vec{\tilde{h}}_n}}+\frac{\sigma_\mathrm{n}^2}{p_u}\vec{I}_{MS}}\right)\vec{w}_k}.
    \label{eq::SINR_instantaneous}
\end{equation}

\section{Pilot Reuse Algorithm}
The coordination of pilot assignment is crucial for mMTC. Due to the massive amount of UEs in such systems, the reuse of the orthogonal pilot sequences in the same cell becomes inevitable. As highlighted by Björnson in~\cite[p. 246]{Bjornson2017}, the strongest interference usually originates from UEs in the same cell. Therefore, a proper assignment of the pilots is essential to mitigate the pilot contamination.

Next, we tackle the pilot contamination problem by proposing pilot reuse strategies that exploit the spatial information present in CSI through the second-order statistics of the radio environment. In particular, we propose two features encompassing the channel correlation, especially in the angular domain, for the pilot reuse: 1) the CMD metric and 2) the CC created based on the CMD metric. These features are fed to a modified version of our greedy pilot allocation algorithm proposed in~\cite{Ribeiro2020} (Algorithm~\ref{alg:NNPA} detailed in Sect.\ \ref{subsec:NNPA}), which completes the pilot assignment, aiming at minimizing the pilot interference between the users sharing the same pilot.

The goal is to assign orthogonal pilot sequences to users with overlapping AoA intervals to avoid pilot contamination. To capture the spatial information embedded in the second-order statistics and reveal the angular domain relationship among the UEs, we use the CMD~\cite{Herdin2005} metric. Thus, we compute a dissimilarity matrix $\vec{D}\,\in\mathbb{R}^{N\times N}$ and define the dissimilarity measure associated with UE $n$ and $j$ as
\begin{equation}\label{eq::CMD}
    d_{n,j}=1-\frac{\trace{\vec{R}_n^{\mbox{\scriptsize H}}\vec{R}_j}}{\fro{\vec{R}_n}\fro{\vec{R}_j}}.
\end{equation}

Next, we elaborate on two different features that utilize the CMD metric in \eqref{eq::CMD} as a basis. The features will serve as input to our \textit{Nearest Neighbor Pilot Assignment} Algorithm 1.

\subsection{CMD-aided Pilot Assignment}
Since the dissimilarity matrix $\vec{D} =[\vec{d}_1,\ldots,\vec{d}_N]$ carries the spatial information needed to allocate the pilot sequences, we propose to directly utilize the distance vectors, $\vec{d}_n = [d_{n,1},\ldots,d_{n,N}]^{\mbox{\scriptsize T}}\,\in\mathbb{R}^{N}$, as input features, ${\vec{f}_n}\,\in\mathbb{R}^{N}$, to the proposed pilot allocation algorithm as detailed in~\ref{subsec:NNPA}. Thus, we form $\vec{f}_n$ by setting
$\vec{f}_n = \vec{d}_n$. After getting $\vec{f}_n$ we apply the Algorithm~\ref{alg:NNPA} to allocate the pilots.

\subsection{CC-aided Pilot Assignment}

CC is as a framework developed in~\cite{Studer2018} to generate unsupervised radio environment mappings utilizing the CSI. The idea behind CC is to find a suitable feature from CSI and then apply a dimensionality reduction (DR) technique to get a lower dimensional embedding, which preserves the relative position of the UEs.
Given the CMD distance in~\eqref{eq::CMD}, we propose to utilize it as an input to CC framework.
The target is to obtain a feature that well characterizes the angular domain structures in the multi-sector multi-user scenario, in order to feed the pilot allocation Algorithm~\ref{alg:NNPA}.

As shown in~\cite{Ribeiro2020}, for single-cell scenarios, CC can be used to coordinate the pilot assignment under pilot reuse regimes, to achieve better network performance than other methods from the literature. In~\cite{Ribeiro2020}, we have utilized CC to retrieve the angular relationship between the UEs by applying the discrete Fourier transform on the channel covariance matrices and then taking the absolute value.
Here, we differently propose to use the CMD dissimilarity matrix
$\vec{D}$ in~\eqref{eq::CMD} as an input to the DR technique, to construct the CC. Thus, we apply a function $\mathcal{C}$ that generates CC by mapping these high-dimensional features $\vec{d}_n\in\mathbb{R}^{N}$ to a low-dimensional domain, i.e.,
\begin{equation}
  \mathcal{C}:\vec{d}_n\mapsto\vec{f}_n,\quad n\in\mathcal{N},
  \label{eq::cc-mapping}
\end{equation}
where ${\vec{f}_n\in\mathbb{R}^{C}}$ is the point in the $C$-dimensional
CC corresponding to $\vec{d}_n$, where typically $C=2$ or $C=3$.

Several unsupervised DR techniques have been proposed to map the extracted features into a lower dimension embedding\cite{Studer2018,Lei2019,Huang2019,Agostini2020}. We propose the use of Laplacian Eigenmaps (LE) as a DR technique as it aims to preserve the local structure of the high-dimensional embedding by minimizing the distances between data points and its $\nu$ nearest neighbors, where $\nu$ is a design parameter. Note that the CC method is not limited to utilize LE, and different DR techniques could be used as well. Clearly, as the low-dimensional embedding is not unique, they may lead to different pilot assignments, thus different network performance.

\subsection{Nearest Neighbor Pilot Assignment Algorithm}
\label{subsec:NNPA}

\begin{algorithm}[t]
\caption{Nearest neighbor pilot assignment}
\label{alg:NNPA}
\SetAlgoNoLine
\DontPrintSemicolon
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
\SetKwFor{For}{for}{do}{end for}

\Input{1) The set of UEs $\mathcal{N}=\{1,\ldots,N\}$, 2) the index set of orthogonal pilots $\mathcal{T}=\{1,\ldots,\tau\}$, and 3) the suitable features $\vec{f}_n,\, n \in\,\mathcal{N}$.}

\Output{A pilot assignment $\vec{\Psi}=\left[\vec{\boldsymbol{\psi}}_1,\ldots,\vec{\boldsymbol{\psi}}_K\right]^{\mbox{\scriptsize T}}$.}

Initialize the set of unassigned UEs $\mathcal{N}^{\mathrm{un}}=\mathcal{N}$, and the set of unassigned pilots $\mathcal{T}^{\mathrm{un}}=\mathcal{T}$.\;

Select a random UE $n$ and initialize the auxiliary variable $n'$ with it, i.e., $n'=n$.\;

Assign $\boldsymbol{\phi}_1$ to user $n$ and update the set of unassigned UEs and pilots, i.e., $\mathcal{N}^{\mathrm{un}}\gets \mathcal{N}^{\mathrm{un}}\setminus \{n\}$ and $\mathcal{T}^{\mathrm{un}}\gets \mathcal{T}^{\mathrm{un}}\setminus \{1\}$.\;

Initialize the auxiliary variable: $p=2$.\;
\While{$\mathcal{N}^{\mathrm{un}}\neq \emptyset$}{
  \If{$\mathcal{T}^{\mathrm{un}}=\emptyset$}{
  Reinitialize: $\mathcal{T}^{\mathrm{un}}=\mathcal{T}$ and $p=1$.\;
   }

  Assign pilot $\boldsymbol{\phi}_{p}$ to user $n$, i.e., $\vec{\boldsymbol{\psi}}_k=\boldsymbol{\phi}_p$, that satisfies $n=\underset{n\in \mathcal{N}^{\mathrm{un}}}{\argmin}\Vert \vec{f}_n-\vec{f}_{n'}\Vert^2$.\;

  Update the set of unassigned UEs, $\mathcal{N}^{\mathrm{un}}\gets \mathcal{N}^{\mathrm{un}}\setminus \{n\}$, and the set of unassigned pilots, $\mathcal{T}^{\mathrm{un}}\gets \mathcal{T}^{\mathrm{un}}\setminus \{p\}$.\;

  Update $n'=n$ and $p=p+1$.\;
 }
\end{algorithm}

After extracting the CMD and CMD + CC features, we apply a slightly modified version of the low-complexity pilot algorithm proposed in~\cite{Ribeiro2020} to account for the different input features and incorporate the inactive users as well, as shown in Algorithm~\ref{alg:NNPA}.
The main premise of the algorithm is that the features preserve the angular relationship among UEs. Thus, the greedy algorithm allocates the orthogonal pilot sequences in an \textit{ordered} way $\boldsymbol{\phi}_1,\boldsymbol{\phi}_2,\ldots,\boldsymbol{\phi}_\tau,\boldsymbol{\phi}_1,\boldsymbol{\phi}_2,\ldots$ by finding the unassigned UE with the small Euclidean distance between the features, aiming at maximizing
the distances between the same pilot sequences.

The first step of Algorithm~\ref{alg:NNPA} is to allocate the first pilot sequence, $\boldsymbol{\phi}_1$, to one of the users. Then, it greedily allocates the next orthogonal pilot sequence, $\boldsymbol{\phi}_2$, to the unassigned UE with the smallest distance to the previously allocated user $n'$, i.e., it finds $n$ that minimizes $\norm{\vec{f}_n-\vec{f}_{n'}}^2$. Then, it repeats this process until all orthogonal pilot sequences have been allocated. Once the last orthogonal pilot sequence $\boldsymbol{\phi}_\tau$ has been allocated, after $\tau$ iterations, it starts to reuse pilot sequences by allocating the first sequence $\boldsymbol{\phi}_1$ to the closest unassigned UE from the precedent allocated user. It repeats this process until all UEs have been assigned pilot sequences.

\section{Numerical and Simulation results}
\label{sec:results}

We consider $N=512$ UEs uniformly distributed in a cell with $K=64$ active UEs at any given time, and a BS equipped with $S=3$ $M$-element ULAs, with $M=64$ critically spaced, $\Delta_r = 0.5$, elements. The propagation channel between each user and the BS consists of $L = 200$ paths with a fixed angular standard deviation $\sigma_\theta=15^{\circ}$. For the antenna parameters in~\eqref{eq::AntennaGain}, we set the maximum antenna gain as $G_{\mathrm{A_{max}}}=0$~dB, the maximum attenuation as $A_{\mathrm{max}}=30$~dB, and the half-power beamwidth as $\theta_{\mathrm{3dB}}=65^\circ$~\cite{TR36873}. We use binary phase shift keying (BPSK) for the channel estimation and quadrature phase shift keying (QPSK) for the data transmission.

\begin{figure}[t]
\centering
\begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth,trim=1mm 0mm 10mm 2mm,clip]{SER-AS15TAU64K64N512DUCE.pdf}
    \caption{SER as a function of SNR.}
    \label{}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth,trim=1mm 0mm 10mm 2mm,clip]{SER-M.pdf}
    \caption{SER as a function of the number of antennas per ULA $M$.}
    \label{}
\end{subfigure}
\caption{Symbol error rate (SER) for a fixed pilot length of $\tau=64$. (a) shows the effect of the SNR and (b) the effect of the number of antennas for a fixed SNR=0~dB.}\vspace{-2mm}
\label{fig:SER}
\end{figure}

\begin{figure}[t]
\centering
\begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth,trim=1mm 0mm 10mm 2mm,clip]{CDF-AR-ASD15TAU64K64N512DUCE-SNR10-linear.pdf
    \caption{CDF.}
    \label{}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth,trim=1mm 0mm 10mm 2mm,clip]{SUMAR-TAU.pdf
    \caption{Achievable sum rate.}
    \label{}
\end{subfigure}
\caption{Achievable sum rate at 10~dB SNR. (a) CDF of users' uplink achievable rate for $\tau=64$.(b) Achievable sum rate vs. pilot length.}\vspace{-2mm}
\label{fig:AR}
\end{figure}

We consider three baseline pilot assignment methods: 1) A random pilot assignment scheme; 2) A real position method, which relies on the exact UEs' positions to compute the angular separation between the UEs, and then deploys Algorithm~\ref{alg:NNPA}. Note that this baseline acts as a lower bound to the CC-based method, since it uses the real position instead of the relative position estimated via CC; 3) The statistical greedy pilot scheduling (SGPS) developed in~\cite{You2015} that, similarly to our method, relies on the knowledge of the channel covariance matrices. Differently from~\cite{You2015}, we extract designed features from the channel covariance matrices and then apply Algorithm~\ref{alg:NNPA}; this also obviates the need to store the channel covariance matrices at the BS.

Fig.~\ref{fig:SER} shows how SER of the different pilot assignment strategies changes with the SNR and the number of antennas $M$. A fixed pilot length of $\tau=64$ is considered, which gives a pilot reuse factor of $N/\tau=8$, meaning that each pilot sequence is reused 8 times across all the UEs connected to the BS. For benchmark, we also show the lower bound for SER, which is achieved with perfect CSI knowledge. From Fig.~\ref{fig:SER}(a) we see that the performance of the CC-based pilot allocation approaches the real position method, beating all the other baseline methods. The CMD-based approach, although implementing the same pilot allocation algorithm, shows poor performance as compared to CC-based, which highlights the importance of choosing a good feature to the pilot allocation Algorithm~\ref{subsec:NNPA}. Fig.~\ref{fig:SER}(b) presents the effect of changing the number of ULA elements $M$ for SNR=0~dB.
By increasing the number of antennas, the performance increases because the UEs' channels become more orthogonal.

We also evaluate the performance of the proposed features in terms of the uplink achievable rate. The cumulative distribution function (CDF) for the average uplink rate per UE is depicted in Fig.~\ref{fig:AR}(a) for a fixed pilot length of $\tau=64$, at 10~dB SNR. It can be seen that, on average, at least 90~\% of the UEs using the CMD-based pilot allocation strategy have a higher uplink achievable rate than the UE with the highest rate employing the random pilot allocation. On the other hand, this percentage grows to 100~\% when the UEs are employing CC-based or SGPS schemes.

Fig.~\ref{fig:AR}(b) plots the uplink sum rate against the pilot length.
One can notice that the CC-based method performs very well, approaching the real position based pilot allocation. However, it slightly degrades its performance for
$\tau>100$ symbols, i.e., when the pilot reuse factor $N/\tau$ decreases below 5. On the other hand, in this regime, the SGPS method performs as well as the real position one. Although the CMD-based method greatly improves the performance as compared to the random pilot allocation, it underperforms when compared to the other methods.

\section{Conclusions}
This paper investigated the reuse of orthogonal pilot sequences in a single-cell multi-sector network for mMIMO systems. In these systems, the transmission rate performance is interference limited due to pilot contamination. We proposed to apply new features capturing the angular information present in CSI in the algorithm developed in~\cite{Ribeiro2020} to accommodate the multi-sector framework and to improve the average uplink rate for the UEs. The CC-based feature was shown to approach the real position method, outperforming the other baselines in terms of SER and achievable uplink rate.
The CMD-based feature performs close
to SGPS and CC-based pilot allocation, whilst having lower complexity. The proposed features showed to fit well into the multi-sector framework by presenting a competitive performance for SER and achievable rate as compared to existing pilot allocation schemes.

\section*{Acknowledgements}

This work has been financially supported in part by the Academy of Finland project ROHM (grants 319484 and 319485). The work of L. Ribeiro, M. Leinonen, and M. Juntti has been financially supported in part by the Academy of Finland 6Genesis Flagship (grant 318927). The work of M. Leinonen has also been financially supported in part by Infotech Oulu and the Academy of Finland (grant 323698).

\bibliographystyle{IEEEtran}
\bibliography{jour_short,conf_short,References.bib}

\end{document}


\section*{Introduction}
\emph{Multi-valued Decision Diagrams} (MDD) are a generalization of
\emph{Binary Decision Diagrams} (BDD) which have long been used in the
verification, e.g., for model checking purposes
\cite{mcmillan:92:symbolic-model-checking}. Recently, these graphical
models have drawn the attention of researchers from the CP and OR communities.
The popularity of these decision diagrams (DD) stems from their ability to
provide a compact representation of large solution spaces as in the case of the
table constraint \cite{regin:15:mdd,helene:18:compact-mdd}.
One of the
research streams which emerged from this increased interest about MDDs is
\emph{decision-diagram-based optimization} (DDO) \cite{bergman:16:theoretical}.
Its purpose is to efficiently solve combinatorial optimization problems by
exploiting problem structure through DDs.
This paper belongs to the DDO sub-field and intends
to further improve the efficiency of DDO solvers through the introduction of two
bounding techniques: local-bounds pruning (LocB) and rough-upper-bound pruning
(RUB).

This paper starts by covering the necessary background on DDO. Then, it presents
the  local-bound and rough-upper-bound pruning techniques in Sections
\ref{sec:local-bounds} and \ref{sec:rub}. After that, it presents an
experimental study which we conducted using `ddo'
\cite{gillard:20:ddo}\footnote{\url{https://github.com/xgillard/ddo}}, our open
source fast and generic MDD-based optimization library. This experimental study
investigates the relevance of RUB and LocB through four disinct NP-hard
problems: the Weighted
Maximum Independent Set Problem (MISP), Maximum Cut Problem (MCP), Maximum 2
Satisfiability Problem (MAX2SAT) and the Traveling Salesman Problem with Time
Windows (TSPTW). Finally, section \ref{sec:related} discusses previous related
work before drawing conclusions.

\section{Background}\label{sec:background}
The coming paragraphs give an overview of discrete optimization with decision
diagrams. Most of the formalism presented here originates from
\cite{bergman:16:branch-and-bound}. Still, we reproduce it here for the sake of
self-containedness.
\paragraph{Discrete optimization.}\label{sec:discrete-optimization}
A discrete optimization problem is a constraint \emph{satisfaction} problem with
an associated objective function to be maximized. The discrete optimization
problem $\pb{P}$ is defined as $\max \set{f(x) \mid x \in D \wedge C(x)}$ where
$C$ is a set of constraints, $x = \tuple{x_0, \dots, x_{n-1}}$ is an assignment
of
values to variables, each of which has an associated finite domain $D_i$ s.t. $D
= D_0 \times \dots \times D_{n-1}$ from where the values are drawn. In that
setup,
the function $f: D \to \mathbb{R}$ is the objective to be maximized.

Among the set of feasible solutions $Sol(\pb{P}) \subseteq D$ (i.e. satisfying
all constraints in $C$), we denote the optimal solution by $x^*$. That is, $x^*
\in Sol(\pb{P})$ and $\forall x \in Sol(\pb{P}) : f(x^*) \ge f(x)$.

\paragraph{Dynamic programming.}
Dynamic programming (DP) was introduced in the mid 50's by Bellman
\cite{bellman:54}. This strategy is significantly popular and is at the heart of
many classical algorithms (e.g., Dijkstra's
algorithm \cite[p.658]{cormen:09:algorithms} or Bellman-Ford's
\cite[p.651]{cormen:09:algorithms}).

Even though a dynamic program is often thought of in terms of recursion, it is
also natural to consider it as a labeled transition system. In that case, the
\emph{DP model} of a given discrete optimization problem $\pb{P}$ consists
of:
\begin{itemize}
    \item a set of state-spaces $S_0, \dots, S_{n}$ among which one
    distinguishes the \emph{initial state} $r$, the \emph{terminal
        state} $t$ and the \emph{infeasible state} $\bot$.
    \item a set of transition functions $t_i : S_i \times D_i \to S_{i+1}$ for
    $i = 0, \dots, n-1$ taking the system from one state $s^{i}$ to the next
    state $s^{i+1}$ based on the value $d$ assigned to variable $x_i$ (or
    to $\bot$ if assigning $x_i = d$ is infeasible). These functions
    should never allow one to recover from infeasibility ($t_i(\bot, d)
    = \bot$ for any $d \in D_i$).
    \item a set of transition cost functions $h_i: S_i \times D_i \to \mathbb{R}$
    representing the immediate reward of assigning some value $d \in D_i$ to the
    variable $x_i$ for $i = 0, \dots, n-1$.
    \item an initial value $v_r$.
\end{itemize}

On that basis, the objective function $f(x)$ of $\pb{P}$ can be formulated as
follows:

\begin{align*}
& \text{maximize~} f(x) = v_r + \sum_{i = 0}^{n-1} h_i(s^i, x_i) \\
& \text{subject to}\\
& s^{i+1} = t_i(s^i, x_i) \text{~for~} i = 0, \dots, n-1 ; x_i \in D_i \wedge
C(x_i) \\
& s^i \in S_i \text{~for~} i = 0, \dots, n
\end{align*}
where $C(x_i)$ is a predicate that evaluates to $true$ when the partial
assignment $\tuple{x_0, \dots ,x_i}$ does not violate any constraint in $C$.

The appeal of such a formulation stems from its simplicity and its
expressiveness which allows it to effectively capture the problem structure.
Moreover, this formulation naturally lends itself to a DD representation; in
which case it represents an exact DD encoding the complete set $Sol(\pb{P})$.

\subsection{Decision diagrams} \label{sec:dd}
Because DDO aims at solving constraint \emph{optimization} problems and not just
constraint \emph{satisfaction} problems, it uses a particular DD flavor known as
reduced weighted DD -- DD as of now. As initially posed by
Hooker\cite{hooker:13:dd-dp}, DDs can be perceived as a compact representation
of the search trees. This is achieved, in this context, by superimposing
isomorphic subtrees.

To define our DD more formally, we will slightly adapt the notation from
\cite{bergman:16:theoretical}.
A DD $\dd{B}$ is a layered directed acyclic graph $\dd{B} = \tuple{n, U, A, l,
    d, v, \sigma}$ where $n$ is the number of variables from the encoded
    problem,
$U$ is a set of nodes; each of which is associated to some state $\sigma(u)$.
The mapping $l: U \to \set{0 \dots n}$ partitions the nodes from $U$ in
disjoint layers $L_0 \dots L_{n}$ s.t. $L_i = \set{u \in U: l(u) = i}$ and the
states of all the nodes belonging to the same layer pertain to the same
DP-state-space ($\forall u \in L_i : \sigma(u) \in S_i$ for $i = 0, \dots,
n$). Also, it should be the case that no two distinct nodes of one same layer
have the same state ($\forall u_1, u_2 \in
L_i: u_1 \ne u_2 \implies \sigma(u_1) \ne \sigma(u_2)$, for $i = 0, \dots,
n$).

The set $A \subseteq U \times U$ from our formal model is a set of directed
arcs connecting the nodes from $U$. Each such arc $a = \arc{u_1, u_2}$ connects
nodes from subsequent layers ($l(u_1) = l(u_2) - 1$) and should be regarded as
the materialization of a branching decision about variable $x_{l(u_1)}$. This is
why all arcs are annotated via the mappings $d: A \to D$ and $v: A \to \mathbb{R}$
which  respectively associate a decision and value (weight) with the given arc.

\begin{example}
    An arc $a$ connecting nodes $u_1 \in L_3$ to $u_2 \in L_4$, annotated with
    $d(a) = 6$ and $v(a) = 42$ should be understood as the assignment $x_3 = 6$
    performed from state $\sigma(u_1)$. It should also be understood that
    $t_3(\sigma(u_1), 6) = \sigma(u_2)$ and the benefit of that assignment is
    $v(a) = h_3(\sigma(u_1), 6) = 42$.
\end{example}

Because each $r$-$t$ path describes an
assignment that satisfies $\pb{P}$, we will use $Sol(\dd{B})$ to denote the set
of all the solutions encoded in the r-t paths of DD $\dd{B}$. Also, because
unsatisfiability is irrecoverable, r-$\bot$ paths are typically omitted from
DDs. It follows that a nice property from using a DD representation
$\dd{B}$ for the DP formulation of a problem $\pb{P}$, is that finding $x^*$ is
as simple as finding the longest r-t path in $\dd{B}$ (according to the relation
$v$ on arcs).

\paragraph{Exact-MDD.} \label{sec:exact}
For a given problem $\pb{P}$, an exact MDD $\dd{B}$ is an MDD that exactly
encodes the solution set $Sol(\dd{B}) = Sol(\pb{P})$ of the problem $\pb{P}$. In
other words, not only do all r-t paths encode valid solutions of $\pb{P}$, but
no feasible solution is present in $Sol(\pb{P})$ and not in $\dd{B}$.
An exact MDD for $\pb{P}$ can be compiled in a top-down fashion\footnote{
    An incremental refinement \emph{a.k.a. construction by separation}
    procedure is detailed in \cite[pp. 51--52]{cire:14:decision-diagrams} but we
    will not cover it here for the sake of conciseness.}.
This naturally follows from the above definition. To that end, one simply
proceeds by a repeated unrolling of the transition relations until all
variables are assigned.

\subsection{Bounded-Size Approximations} \label{sec:approx}
In spite of the compactness of their encoding, the construction of DD suffers
from a potentially exponential memory requirement in the worst
case\footnote{
    Consequently, it also suffers from a potentially exponential time
    requirement in the worst case. Indeed, time is constant in the final number
    of nodes (unless the transition functions themselves are exponential in the
    input).}.
Thus, using DDs to exactly encode the solution space of a problem is often
intractable. Therefore, one must resort to the use of \emph{bounded-size}
approximation of the exact MDD. These are compiled generically by inserting a
call to a width-bounding procedure to ensure that the width (the number
$|L_i|$ of distinct nodes belonging to the $L_i$) of the current layer $L_i$
does not exceed a given bound $W$. Depending on the behavior of that procedure,
one can either compile a  restricted-MDD (= an under-approximation) or a
relaxed-MDD (= an over-approximation).

\paragraph{Restricted-MDD: Under-approximation.} \label{sec:restricted}
A restricted-MDD provides an under-approximation of some exact-MDD. As such, all
paths of a restricted-MDD encode valid solutions, but some solutions might be
missing from the MDD. This is formally expressed as follows: given the DP
formulation of a problem $\pb{P}$, $\dd{B}$ is a restricted-MDD iff $Sol(\dd{B})
\subseteq Sol(\pb{P})$.

To compile a restricted-MDD, it is sufficient to simply delete certain nodes
from the current layer until its width fits within the specified bound $W$. To
that end, the width-bounding procedure simply selects a subset of the nodes from
$L_i$ which are heuristically assumed to have the less impact on the tightness
of the bound. Various heuristics have been studied in the literature
\cite{bergman:14:optimization-bounds}, and \emph{minLP} was shown to be the
heuristic that works best in practice. This heuristic decides to select (hence
remove) the nodes having the shortest longest path from the root.

\paragraph{Relaxed-MDD: Over-approximation.} \label{sec:relaxed}
A relaxed-MDD $\dd{B}$ provides a bounded-width over-approximation of some
exact-MDD. As such, it may hold paths that are no solution to $\pb{P}$, the
problem being solved. We have thus formally that $Sol(\dd{B}) \supseteq
Sol(\pb{P})$.

Compiling a relaxed-MDD requires one to be able to \emph{merge} several nodes
into an inexact one. To that end, we use two operators:
\begin{itemize}
    \item $\oplus$ which yields a new node combining the states of a selection
    of nodes so as to over-approximate the states reachable in the selection.
    \item $\Gamma$ which is used to possibly relax the weight of arcs incident to
    the selected nodes.
\end{itemize}

These operators are used as follows. Similar to the restricted-MDDs case, the
width-bounding procedure starts by heuristically selecting the least promising
nodes and removing them from layer $L_i$. Then the states of these selected
nodes are combined with one another so as to create a merged node $\mathcal{M} =
\oplus(selection)$. After that, the inbound arcs incident to all selected nodes
are $\Gamma$-relaxed and redirected towards $\mathcal{M}$. Finally, the
result of the merger ($\mathcal{M}$) is added to the layer in place of the
initial selection of nodes.

\paragraph{Summary.}
Fig.~\ref{fig:mdd:exact-restricted-relaxed} summarizes the information from
sections \ref{sec:dd} and \ref{sec:relaxed}. It displays the three MDDs
corresponding to one same example problem having four variables. The exact MDD
(a) encodes the complete solution set and, equivalently, the state space of the
underlying DP encoding. One easily notices that the
restricted DD (b) is an under approximation of (a) since it achieves its width
boundedness by removing nodes d and e and their children (i, j). Among others,
it follows that the solution $[ x_0= 0, x_1 = 0, x_2 = 0, x_3 = 0 ]$ is not
represented in (b) even though it exists in (a). Conversely, the relaxed diagram
(c) achieves a maximum layer with of 3 by merging nodes d, e and h into a new
inexact node $\mathcal{M}$ and by relaxing all arcs entering one of the merged
nodes. Because of this, (c) introduces solutions that do not exist in (a) as is
for instance the case of the assignment $[ x_0= 0, x_1 = 0, x_2 = 3, x_3 = 1 ]$.
Moreover, because the operators $\oplus$ and $\Gamma$ are correct\footnote{The
    very definition of these operators is problem-specific. However,
    \cite{hooker:17:job-sequencing} formally defines the conditions that are
    necessary to correctness.}, the length of the longest path in (c) is an
    upper
bound on the optimal value of the objective function. Indeed, one can see that
the length of the longest path in (a) (= the exact optimal solution) has a value
of 25 while it amounts to 26 in (c).

\begin{figure*}[h]
    \centering
    \scalebox{.7}{
        \scriptsize
        \tikzfig{graphs/exact_restricted_relaxed}
    }
    \caption{
        The exact (a), restricted (b) and relaxed (c) versions of an
        MDD with four variables. The width of MDDs (b) and (c) have been
        bounded to a maximum layer width of three. The decision labels of the
        arcs are shown above the layers separation lines (dashed). The arc
        weights are shown below the layer separation lines. The longest path of
        each MDD is boldfaced. In (c), the node $\mathcal{M}$ is the result of
        merging nodes d, e and h with the $\oplus$ operator. Arcs that
        have been relaxed with the $\Gamma$ operator are pictured with a double
        stroke. Note, because these arcs have been $\Gamma$-relaxed, their value
        might be greater than that of corresponding arcs in (a), (b).
        Similarly, all ``inexact'' nodes feature a double border.
    }
    \label{fig:mdd:exact-restricted-relaxed}
\end{figure*}

\subsection{The Dynamics of Branch-and-Bound with DDs}\label{sec:bab-dd}
Being able to derive good lower and upper  bounds for some optimization problem
$\pb{P}$ is useful when the goal is to use these bounds to strengthen algorithms
\cite{davarnia:18:approx-non-linear,tjandraatmadja:18:decision-diagrams,tjandraatmadja:19:cuts-from-relaxed}.
But it is not the only way these approximations can be used. A complete
and efficient branch-and-bound algorithm relying on those approximations was
proposed in \cite{bergman:16:branch-and-bound} which we hereby reproduce
(\alg{alg:branch-and-bound}).

This algorithm works as follows: at start, the node $r$ is created for the
initial state of the problem and placed onto the \emph{fringe} -- a global
priority queue that tracks all nodes remaining to explore and orders them from
the most to least promising. Then, a loop consumes the nodes from that fringe
(line \ref{alg:branch-and-bound:1a}), one at a time and explores it until the
complete state space has been exhausted. The \emph{exploration} of a node $u$
inside that loop proceeds as follows: first, one compiles a restricted DD
$\rst{B}$ for the sub-problem rooted in $u$ (line \ref{alg:branch-and-bound:1b}).
Because all paths in a restricted DD are feasible solutions, when the lower
bound $v^*(\rst{B})$ derived from the restricted DD $\rst{B}$ improves over the
current best known solution $\lb{v}$; then the longest path of $\rst{B}$ (best
sol. found in $\rst{B}$) and its length $v^*(\rst{B})$ are memorized (lines
\ref{alg:branch-and-bound:improve-lb-start}-\ref{alg:branch-and-bound:improve-lb-end}).

In the event where $\rst{B}$ is exact (no restriction occurred during the
compilation of $\rst{B}$), it covers the complete state space of the sub-problem
rooted in $u$. Which means the processing of $u$ is complete and we may safely
move to the next node. When this condition is not met, however, some additional
effort is required
In that case, a \emph{relaxed} DD $\rlx{B}$ is compiled from $u$(line
\ref{alg:branch-and-bound:relax}). That relaxed DD serves two purposes: first,
it is used to derive an upper bound $v^*(\rlx{B})$ which is compared to the
current best known solution (line \ref{alg:branch-and-bound:2b}). This gives us
a chance to prune the unexplored state space under $u$ when $v^*(\rlx{B})$
guarantees it does not contain any better solution than the current best.
The second use of $\rlx{B}$ happens when $v^*(\rlx{B})$ cannot provide such a
guarantee. In that case, the exact cutset of $\rlx{B}$ is used to enumerate
residual sub-problems which are enqueued onto the fringe (lines
\ref{alg:bab-cutset-start}-\ref{alg:bab-cutset-end}).

A cutset for some relaxed DD $\rlx{B}$ is a subset $\mathcal{C}$ of
the nodes from $\rlx{B}$ such that any $r-t$ path of $\rlx{B}$ goes through at
least one node $\in \mathcal{C}$. Also, a node $u$  is said to be exact iff all
its incoming paths lead to the same state $\sigma(u)$. From there, an exact
cutset of $\rlx{B}$ is simply a cutset whose nodes are all exact.
Based on this definition, it is easy to convince
oneself that an exact cutset constitutes a frontier up to which the relaxed DD
$\rlx{B}$ and its exact counterpart $\dd{B}$ have not diverged. And, because it
is a
cutset, the nodes composing that frontier cover all paths from both $\dd{B}$ and
$\rlx{B}$; which guarantees the completeness of
\alg{alg:branch-and-bound} \cite{bergman:16:branch-and-bound}.

Any relaxed-MDD admits at least one exact cutset -- e.g. the trivial $\set{r}$
case. Often though, it is not unique and different options exist as to what
cutset to use. It was experimentally shown by \cite{bergman:16:branch-and-bound}
that most of  the time, the Last Exact Layer (LEL) is superior to all other
exact cutsets in practice. LEL consists of the \emph{deepest} layer of the
relaxed-MDD having all its nodes exact.
\begin{example}
    In Fig.-\ref{fig:mdd:exact-restricted-relaxed} (c), the first
    inexact node $\mathcal{M}$ occurs in layer $L_2$. Hence, the LEL cutset
    comprises all nodes (a, b, c) from the layer $L_1$. Because $\mathcal{M}$
    is inexact, and because it is a parent of nodes i, j and k, these three
    nodes are considered inexact too.
\end{example}
\vspace{-2em}
\noindent
\begin{center}
\scalebox{0.8}{
\begin{minipage}[t]{.6\textwidth}
\begin{algorithm}[H]
    \begin{algorithmic}[1]
        \STATE Create node $r$ and add it to $Fringe$
        \label{alg:branch-and-bound:1a}
        \STATE $\lb{x} \gets \bot$
        \STATE $\lb{v} \gets -\infty$
        \WHILE{$Fringe$ is not empty}
        \STATE $u \gets Fringe.pop()$ \label{alg:branch-and-bound:1b}
        \STATE $\rst{B} \gets Restricted(u)$ \label{alg:branch-and-bound:2a}
        %
        \IF{$v^*(\rst{B}) > \lb{v}$}\label{alg:branch-and-bound:improve-lb-start}
        \STATE $\lb{v} \gets v^*(\rst{B})$
        \STATE $\lb{x} \gets x^*(\rst{B})$
        \ENDIF\label{alg:branch-and-bound:improve-lb-end}
        %
        \IF{$\rst{B}$ is not exact} \label{alg:branch-and-bound:prune}
        \STATE $\rlx{B} \gets Relaxed(u)$\label{alg:branch-and-bound:relax}
        \IF{$v^*(\rlx{B}) > \lb{v}$} \label{alg:branch-and-bound:2b}
        \FORALL{$u' \in \rlx{B}.exact\_cutset()$} \label{alg:bab-cutset-start}
        \label{alg:branch-and-bound:cutset}
        \STATE $Fringe.add(u')$
        \ENDFOR \label{alg:bab-cutset-end}
        \ENDIF \label{alg:branch-and-bound:3}
        \ENDIF
        \ENDWHILE
        \RETURN $\lst{\lb{x}, \lb{v}}$ \label{alg:branch-and-bound:4}
    \end{algorithmic}
    \caption{Branch-And-Bound with DD}
    \label{alg:branch-and-bound}
\end{algorithm}
\end{minipage}
}
\hfill
\scalebox{0.8}{
    \begin{minipage}[t]{.6\textwidth}
        \begin{algorithm}[H]
            \begin{algorithmic}[1]
                \color{gray}
                \STATE Create node $r$ and add it to $Fringe$
                \STATE $\lb{x} \gets \bot$
                \STATE $\lb{v} \gets -\infty$
                \WHILE{$Fringe$ is not empty}
                    \STATE $u \gets Fringe.pop()$
                    \color{black}
                    \IF{$v|_u^* \le \lb{v}$}
                        \label{alg:reasoning:local-bounds:pop-skip}
                        \STATE \textbf{continue}
                    \ENDIF
                    \color{gray}
                    \STATE $\rst{B} \gets Restricted(u)$
                    %
                    \IF{$v^*(\rst{B}) > \lb{v}$}
                        \STATE $\lb{v} \gets v^*(\rst{B})$
                        \STATE $\lb{x} \gets x^*(\rst{B})$
                    \ENDIF
                    %
                    \IF{$\rst{B}$ is not exact}
                        \STATE $\rlx{B} \gets Relaxed(u)$
                        \IF{$v^*(\rlx{B}) > \lb{v}$}
                            \color{black}
                            \FORALL{$u' \in \rlx{B}.exact\_cutset()$}
                                \IF{$v|_{u'}^* > \lb{v}$}
                                    \label{alg:reasoning:local-bounds:skip-enqueue}
                                    \STATE $Fringe.add(u')$
                                    \label{alg:reasoning:local-bounds:do-enqueue}
                                \ENDIF
                            \ENDFOR
                            \color{gray}
                        \ENDIF
                    \ENDIF
                \ENDWHILE
                \RETURN $\lst{\lb{x}, \lb{v}}$
            \end{algorithmic}
            \caption{Local bound pruning}
            \label{alg:reasoning:local-bounds}
        \end{algorithm}
    \end{minipage}
}
\end{center}

\section{Improving the filtering of branch-and-bound MDD}
In the forthcoming paragraphs, we introduce the local bound and present the
rough upper bound: two reasoning techniques to reinforce the pruning strength
of \alg{alg:branch-and-bound}.

\subsection{Local bounds (LocB)}
\label{sec:local-bounds}
Conceptually, pruning with local bounds is rather simple: a relaxed MDD
$\rlx{B}$ provides us with \emph{one} upper bound $v^*(\rlx{B})$
on the optimal value of the objective function for some given sub-problem.
However, in the event where $v^*(\rlx{B})$ is greater than the best known lower
bound $\lb{v}$ (best current solution) nothing guarantees that all nodes from
the exact cutset of $\rlx{B}$ admit a longest path to $t$ with a length of
$v^*(\rlx{B})$. Actually, this is quite unlikely. This is why we propose to
attach a \emph{``local'' upper bound} to each node of the cutset. This local
upper bound -- denoted $v|_u^*$ for some cutset node $u$ -- simply records
the length of the longest r-t path passing through $u$ in the relaxed MDD
$\rlx{B}$.

In other words, LocB  allows us to refine the information provided by a relaxed
DD $\rlx{B}$. On one hand, $\rlx{B}$ provides us with $v^*(\rlx{B})$ which is
the length of the longest r-t path in $\rlx{B}$. As such, it provides an upper
bound on the optimal value that can be reached from the root node of $\rlx{B}$.
With the addition of LocB, the relaxed DD provides us with an additional piece
of information. For each individual node $u$ in the exact cutset of $\rlx{B}$,
it defines the value $v|_u^*$ which is an upper bound on the value attainable
from that node.

As shown in \alg{alg:reasoning:local-bounds}, the value $v|_u^*$ can prove
useful at two different moments. First, in the event where $v|_u^* \le \lb{v}$,
this value can serve as a justification to not enqueue the subproblem $u$ (line
\ref{alg:reasoning:local-bounds:skip-enqueue}) since exhausting this subproblem
will yield no better solution than $\lb{v}$.
More formally, by definition of a cutset and of LocB, it must be the case that
the longest r-t path of $\rlx{B}$ traverses one of the cutset nodes $u$ and thus
that $v^*(\rlx{B}) = v|_u^*$ (where $v|_u^*$ is the local bound of $u$). Hence
we have: $\exists u \in \text{~cutset of~} \rlx{B} : v^*(\rlx{B}) = v|^*_u$.
However, because $v^*(\rlx{B})$ is the length of the \emph{longest} r-t path of
$\rlx{B}$, there may exist cutset nodes that only belong to r-t paths shorter
than $v^*(\rlx{B})$. That is: $
\forall u’ \in \text{~cutset of~} \rlx{B} : v^*(\rlx{B}) \ge v|_{u’}^*$.
Which is why $v|_{u’}^*$ can be stricter than $v^*(\rlx{B})$ and hence let
LocB be stronger at pruning nodes from the frontier.

The second time when $v|_u^*$  might come in handy occurs when the node $u$ is
popped out of the fringe (line \ref{alg:reasoning:local-bounds:pop-skip}).
Indeed, because the fringe is a global priority queue, any node that has been
pushed on the fringe can remain there for a long period of time. Thus, chances
are that the value $\lb{v}$ has increased between the moment when the node was
pushed onto the fringe (line \ref{alg:reasoning:local-bounds:do-enqueue}) and
the moment when it is popped out of it. Hence, this gives us an additional
chance to completely skip the exploration of the sub-problem rooted in $u$.

Let us illustrate that with the relaxed MDD shown on
Fig.\ref{fig:reasoning:local-bounds}, for which the exact cutset comprises the
highlighted nodes $a$ and $b$. Please note that because this scenario may occur
at any time during the problem resolution, we will assume that the fringe is not
empty when it starts.
Assuming that the current best solution $\lb{v}$ is $20$ when one explores the
pictured subproblem, we are certain that exploring the subproblem rooted in $a$
is a waste of time, because the local bound $v|_a^*$ is only $16$.
Also, because the fringe was not empty, it might be the case that $b$ was left
on the
fringe for a long period of time. And because of this, it might be the case that
the best known value $\lb{v}$ was improved between the moment when $b$ was
pushed on the fringe and the moment when it was popped out of it.
Assuming that $\lb{v}$ has improved to $110$ when
$b$ is popped out of the fringe, it may safely be skipped because $v|_b^*$
guarantees that an exploration of $b$ will not yield a better solution than
$102$.

\alg{alg:computing-local-bounds} describes the procedure to compute the local
bound $v|^*_u$ of each node $u$ belonging to the exact cutset of a relaxed MDD
$\rlx{B}$. Intuitively, this is achieved by doing a bottom-up traversal of
$\rlx{B}$, starting at $t$ and stopping when the traversal crosses the last
exact layer (line \ref{alg:computing-local-bounds:end-cond}). During that
bottom-up traversal, the algorithm marks the nodes that are reachable from $t$.
This way, it can avoid the traversal of dead-end nodes. Also,
\alg{alg:computing-local-bounds} maintains a value $v^*_{\uparrow_t}(u)$ for
each node $u$ it encounters. This value represents the length of the longest u-t
path. Afterwards (line \ref{alg:computing-local-bounds:save-lb}), it is summed
with the length of the longest r-u path $v^*_{r-u}$ to derive the exact value of
the local bound $v|^*_u$.

\begin{algorithm}
    \begin{algorithmic}[1]
        \STATE $lel \gets$ \text{Index of the last exact layer}
        \STATE $v^*_{\uparrow_t}(u) \gets -\infty ~\text{\textbf{for each
        node}}~u \in \rlx{B}$ \hfill\texttt{// init. longest u-t path}
        \STATE $mark(t) \gets $true
        \STATE $v^*_{\uparrow_t}(t) \gets 0$ \hfill\texttt{// longest t-t path}
        \FORALL{$i = n$ to $lel$}\label{alg:computing-local-bounds:end-cond}
           \FORALL{node $u \in L_{i}$}
                \IF{$mark(u)$}
                    \FORALL{arc $a = (u', u)$ incident to $u$}
                        \STATE $mark(u') \gets$ true
                        \STATE $v^*_{\uparrow_t}(u') \gets
                        \max ( v^*_{\uparrow_t}(u'), v^*_{\uparrow_t}(u) + v(a)
                        )
                        $  \hfill\texttt{// longest u'-t path}
                        \label{alg:computing-local-bounds:lp_from_bot}
                    \ENDFOR
                \ENDIF
           \ENDFOR
        \ENDFOR
        \FORALL{node $u \in \rlx{B}.exact\_cutset()$}
            \IF{$mark(u)$}
                \STATE $v|^*_u \gets v^*_{r-u} + v^*_{\uparrow_t}(u)$
                \hfill\texttt{// longest r-u path + longest u-t path}
                \label{alg:computing-local-bounds:save-lb}
            \ELSE
                \STATE $v|^*_u \gets -\infty$
                \label{alg:computing-local-bounds:ditch-not-marked}
            \ENDIF
        \ENDFOR
    \end{algorithmic}
    \caption{Computing the local bounds}
    \label{alg:computing-local-bounds}
\end{algorithm}

\paragraph{Caveat for the First Exact Layer Cutset (FEL)}
\alg{alg:computing-local-bounds} is sufficient to cope with the cases where
$\rlx{B}$ uses either a LEL cutset or an FC cutset. However, in the event where
$\rlx{B}$ would implement another strategy -- such as ie the First Exact Layer
(FEL), then \alg{alg:computing-local-bounds} would need to be adapted and
perform a complete backward traversal of $\rlx{B}$. Concretely, this means that
line \ref{alg:computing-local-bounds:end-cond} would need to read $\text{for
    all } i \in n \text{ down to } 0$ instead of $\text{down to } lel$. This
change is mandatory. Failing to implement it could result in too many nodes
being pruned away because of the line
\ref{alg:computing-local-bounds:ditch-not-marked}.

\subsubsection{Complexity matters...}
As shown in \alg{alg:computing-local-bounds}, the computation of the local
bound $v|^*_u$ of some node $u$ from the exact cutset of a relaxed MDD
$\rlx{B}$ only requires the ability to remember if $u$ is reachable from $t$
($mark(u)$) and to sum the lengths of the longest r-u path $v^*_{r-u}$ and the
longest u-t path $v^*_{\uparrow_t}(u)$. To that end, it is sufficient to only
store one flag and two integers in each node. Thus, the spatial complexity of
implementing local bounds is $\Theta(1)$ per node.

Similarly, the backwards traversal of $\rlx{B}$ as performed in
\alg{alg:computing-local-bounds} only visits a subset of the nodes and arcs
that have been created during the compilation of $\rlx{B}$. The
time complexity of \alg{alg:computing-local-bounds} is thus bounded by
$O(|U|\times|A|)$. Actually, we even know that this bound is lesser or equal to
that of the top-down compilation because the width-bounding procedure
reduces the number of nodes in $\rlx{B}$. Hence it reduces
the number of nodes potentially visited by \alg{alg:computing-local-bounds}. It
follows that the use of local bounds has no impact on the time complexity of
the derivation of a relaxed MDD.

\subsection{Rough upper bound (RUB)}\label{sec:rub}
Rough upper bound pruning departs from the following observation: assuming the
knowledge of a lower bound $\underline{v}$ on the value of $v^*$, and assuming
that one is able to swiftly compute a rough upper bound $\overline{v_s}$ on the
optimal value $v^*_s$ of the subproblem rooted in state $s$; any node $u$ of a
MDD having a rough upper bound
$\overline{v_{\sigma(u)}} \le \underline{v}$ may be discarded as it is
guaranteed not to improve the best known solution. This is pretty much the same
reasoning that underlies the whole branch-and-bound idea. But here, it is used
to prune portions of the search space explored \emph{while compiling}
approximate MDDs.

To implement RUB, it suffices to adapt the MDD compilation procedure (top-down ,
iterative refinement, ...) and introduce a check that avoids creating a node
$u'$ with state $next$ when $\overline{v_{next}} \le \underline{v}$.

The key to RUB effectiveness is that RUB is used while compiling the restricted
and relaxed DDs. As such, its computation does not directly appear in Alg. 1,
but rather is accounted within the compilations of $Restricted(u)$ and
$Relaxed(u)$
from Alg. 1. Thus, it really is not used as yet-an-other-bound competing with
that of line 12, but instead to speed up the computation of restricted and
relaxed DDs. More precisely, this speedup occurs because the compilation of the
DDs discards some nodes that would otherwise be added to the next layer of the
DD and then further expanded, which are ruled out by RUB.
A second benefit of using RUBs is that it helps tightening the bound derived
from a relaxed DD (Alg.1 line 12). Because the layers that are generated in a
relaxed DD are narrower when applying RUB, there are fewer nodes exceeding the
maximum layer width. The operator $\oplus$ hence needs to merge a smaller set of
nodes in order to produce the relaxation.

The dynamics of RUB is graphically illustrated by
Fig.-\ref{fig:reasoning:rough-ub} where the set of highlighted nodes can be
safely elided since the (rough) upper bound computed in node $s$ is lesser than
the best lower bound.

\begin{figure}[ht]
    %
    \begin{minipage}{.48\textwidth}
        \centering
        \scriptsize
        \ctikzfig{graphs/local_bounds}
        \caption{
            An example relaxed-MDD having an exact cutset $\set{a, b}$ with
            local bounds $v|_a^*$ and $v|_b^*$. The nodes with a simple border
            represent exact nodes and those with a double border represent
            ``inexact'' nodes. The edges along the longest path are displayed in
            bold.}
        \label{fig:reasoning:local-bounds}
    \end{minipage}
    \hfill
    \begin{minipage}{.48\textwidth}
        \centering
        \scriptsize
        \ctikzfig{graphs/rough_ub}
        \caption{
            Assuming a lower bound $\lb{v}$ of 100 and a rough upper bound
            $\ub{v_s}$ of 42 for the node $s$, all the highlighted nodes (in
            red, with a dashed border) may be pruned from the MDD.}
        \label{fig:reasoning:rough-ub}
    \end{minipage}
    %
\end{figure}

\subsubsection{Important Note}
It is important to understand that because the RUB is computed at each node of
each restricted and relaxed MDD compiled during the instance resolution, it must
be extremely inexpensive to compute. This is why RUB is best obtained from a
fast and simple problem specific procedure.

\section{Experimental Study}
In order to evaluate the impact of the pruning techniques
proposed above, we conducted a series of experiments on four problems.
In particular, we conducted experiments on the Maximum Independent Set Problem
(MISP), the Maximum Cut Problem (MCP), the Maximum Weighted 2-Satisfiablility
Problem (MAX2SAT) and the Traveling Salesman Problem with Time Windows (TSPTW).
For the first three problems, we generated sets of random instances which we
attempted to solve with different configurations of our own open source solver
written in Rust
\cite{gillard:20:ddo}\footnote{https://github.com/xgillard/ddo}. For TSPTW, we
reused openly available sets of benchmarks which are usually used to assess the
efficiency of new solvers for TSPTW\cite{lopez:20:tsptw}.
Thanks to the generic nature of our framework, the model and all heuristics used
to solve the instances were the same for all experiments. This allowed us to
isolate the impact of RUB and LocB on the solving performance and neutralize
unrelated factors such as variable ordering. Indeed, the only variations between
the different solver flavors relate to the presence (or absence) of RUB and LocB.
All experiments were run on the same physical machine equipped with an AMD6176
processor and 48GB of RAM. A maximum time limit of 1800 seconds was allotted to
each configuration to solve each instance.

The models we used for MISP, MCP and MAX2SAT are direct translations of
the ones described in \cite{bergman:16:branch-and-bound}. The details of the
RUBs we formulated for these problems are given in the appendixes
\ref{appendix:rub:misp}, \ref{appendix:rub:mcp} and \ref{appendix:rub:max2sat}.
The DP formulation we used for TSPTW, its relaxation and RUB are detailed in
appendix \ref{appendix:rub:tsptw}.

\paragraph{MISP.}\label{sec:xp:misp}
To assess the impact of RUB and LocB on MISP, we generated random graphs based
on the Erdos-Renyi model G(n, p) \cite{erdos-renyi:59} with the number of
vertices n = 250, 500, 750, 1000, 1250, 1500, 1750 and the probability of having
an edge connecting any two vertices p = 0.1, 0.2, ... , 0.9. The weight of the
edges in the generated graphs were drawn uniformly from the set $\{-5, -4, -3,
-2, -1, 1, 2, 3, 4, 5\}$. We generated 10 instances for each combination of size
and density (n, p).

\paragraph{MCP.}\label{sec:xp:mcp}
In line with the strategy used for MISP, we generated random MCP instances as
random graphs based on the Erdos-Renyi model G(n, p). These graphs were generated
with the number of vertices n = 30, 40, 50 and the probability p
of connecting any two vertices = 0.1, 0.2, 0.3, .., 0.9. The
weights of the edges in the generated graphs were drawn uniformly among $\{-1,
1\}$. Again, we generated 10 instances per combination n, p.

\paragraph{MAX2SAT.}\label{sec:xp:max2sat}
Similar to the above, we used random graphs based the Erdos-Renyi model G(n, p)
to derive MAX2SAT instances. To this end, we produced graphs with n =
60, 80, 100, 200, 400, 1000 (hence instances with 30, 40, 50, 100, 200 and 500
variables) and p = 0.1, 0.2, 0.3, .. , 0.9. For each combination of size
(n) and density (p), we generated 10 instances. The weights of the clauses
in the generated instances were drawn uniformly from the set $\{1, 2, 3, 5, 6,
7, 8, 9, 10\}$.

\paragraph{TSPTW.}\label{sec:xp:tsptw}
To evaluate the effectiveness of our rules on TSPTW, we used the 467 instances
from the following suites of benchmarks, which are usually used to assess the
efficiency of new TSPTW solvers. AFG \cite{bench:afg}, Dumas \cite{bench:dumas},
Gendreau-Dumas \cite{bench:gendreau-dumas}, Langevin \cite{bench:langevin},
Ohlmann-Thomas \cite{bench:ohlmann-thomas},
Solomon-Pesant \cite{bench:solomon-pesant} and
Solomon-Potvin-Bengio \cite{bench:potvin-bengio}.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{graphs/combo_new.pdf}
    \caption{Number of solved instances over time for each considered problem}
    \label{fig:xp:solved}
\end{figure}

Figure \ref{fig:xp:solved} gives an overview of the results from our
experimental study. It respectively depicts the evolution over time
of the number of instances solved by each technique for MISP (a), MCP (b) and
MAX2SAT (c) and TSPTW (d).

As a first step, our observation of the graphs will focus on the differences
that arise between the single threaded configurations of our \emph{ddo} solvers.
Then, in a second phase, we will incorporate an existing state-of-the-art ILP
solver (Gurobi 9.0.3) in the comparison. Also, because both Gurobi and our
\emph{ddo} library come with built-in parallel computation capabilities, we will
consider both the single threaded and parallel (24 threads) cases. This second
phase, however, only bears on MISP, MCP and MAX2SAT by lack of a Gurobi TSPTW
model.

\paragraph{DDO configurations}
The first observation to be made about the four graphs in
Fig.\ref{fig:xp:solved}, is that for all considered problems, both RUB and
LocB outperformed the 'do-nothing' strategy; thereby showing the relevance of
the rules we propose. It is not clear however which of the two rules brings the
most improvement to the problem resolution. Indeed, RUB seems to be the driving
improvement factor for MISP (a) and TSPTW (d) and the impact of LocB appears to
be moderate or weak on these problems. However, it has a much higher impact
for MCP (b) and MAX2SAT (c). In particular, LocB appears to be the driving
improvement factor for MCP (b). This is quite remarkable given that LocB
operates in a purely black box fashion, without any problem-specific knowledge.
Finally, it should also be noted that the use of RUB and LocB are not mutually
exclusive. Moreover, it turns out that for all considered problems, the
combination RUB+LocB improved the situation over the use of any single rule.

Furthermore, Fig.\ref{fig:xp:gap} confirms the benefit of
using both RUB \emph{and} LocB together rather than using any single technique.
For each problem, it measures the “performance” of using RUB+LocB vs the best
single technique through the end gap. The end gap is defined as
$\left(100 * \frac{|UB|-|LB|}{|UB|}\right)$. This metric allows us to account
for all instances, including the ones that could not be solved to optimality.
Basically, a small end gap means that the solver was able to confirm a tight
confidence interval of the optimum. Hence, a smaller gap is better. On each
subgraphs of Fig.\ref{fig:xp:gap}, the distance along the x-axis represents the
end gap for reach instance when using both RUB and LocB whereas the distance
along y-axis represents the end gap when using the best single technique for the
problem at hand. Any mark above the diagonal shows an instance for which using
both RUB and LocB helped reduce the end gap and any mark below that line
indicates an instance where it was detrimental.

From graphs \ref{fig:xp:gap}-a, \ref{fig:xp:gap}-c and \ref{fig:xp:gap}-d it
appears that the combination RUB+LocB supersedes the use of RUB only.
Indeed the vast majority of the marks sit above the diagonal and the rest on it.
This indicates a beneficial impact of using both techniques even for the hardest
(unsolved) instances. The case of MCP (graph \ref{fig:xp:gap}-b) is less clear
as most of the marks sit on the diagonal. Still, we can only observe three marks
below the diagonal and a bit more above it. Which means that even though the use
of RUB in addition to LocB is of little help in the case of MCP, its use does
not degrade the performance for that considered problem.

\begin{figure}
    \centering
    \includegraphics[scale=.55]{graphs/combo_gaps.pdf}
    \caption{End gap: The benefit of using both techniques vs the best single
    one}
    \label{fig:xp:gap}
\end{figure}

\paragraph{Comparison with Gurobi 9.0.3}
The first observation to be made when comparing the performance of Gurobi vs the
DDO configurations, is that when running on a single thread, ILP outperforms the
basic DDO approach (without RUB and LocB). Furthermore, Gurobi turns out to be
the best single threaded solver for MCP by a fair margin. However, in the MISP
and MAX2SAT cases, Fig.~\ref{fig:xp:solved} shows that the DDO solvers
benefitting
from RUB and LocB were able to solve more instances and to solve them faster
than Gurobi. Which underlines the importance of RUB and LocB.

When lifting the one thread limit, one can see that the DD-based approach
outperform ILP on each of the considered problems. In particular, in the case of
MCP for which Gurobi is the best single threaded option; our DDO solver was
able to find and prove the optimality of all tested instances in a little less
than 800 seconds. The ILP solver, on the other end, was not able to prove the
optimality of the 9 hardest instances within 30 minutes. Additionally, we also
observe that in spite of the performance gains of MIP when running in parallel,
Gurobi fails to solve as many MISP and MAX2SAT instances and to solve them as
fast as the single threaded DDO solvers with RUB and LocB. This emphasizes once
more the relevance of our techniques.
It also shows that the observation from \cite{bergman:14:parallel} still hold
today: despite the many advances of MIP the DDO approach still scales better
than MIP on the considered problems when invoked in parallel.

\section{Previous work}\label{sec:related}
DDO emerged in the mid' 2000's when \cite{hooker:06:discrete-optimization}
proposed to use decision diagrams as a way to solve discrete optimization
problems to optimality. More or less concomitantly, \cite{hooker:07:relaxed-mdd}
devised relaxed-MDD even though the authors envisioned its use as a CP
constraint store rather than a means to derive tight upper bounds for
optimization problems. Then, the relationship between decision diagrams and
dynamic programming was clarified by \cite{hooker:13:dd-dp}.

Recently, Bergman, Ciré and van Hoeve investigated the various ways to
compile decision diagrams for optimization (top-down, construction by
separation) \cite{cire:14:decision-diagrams}. They also investigated the
heuristics used to parameterize these DD compilations. In particular, they
analyzed the impact of variable ordering in
\cite{cire:14:decision-diagrams,bergman:14:optimization-bounds} and node
selection heuristics (for merge and deletion) in
\cite{bergman:14:optimization-bounds}. Doing so, they empirically demonstrated
the crucial impact of variable ordering on the tightness of the derived bounds
and highlighted the efficiency of minLP as a node selection heuristic. Later
on, the same authors proposed a complete branch-and-bound algorithm based on DDs
\cite{bergman:16:branch-and-bound}. This is the algorithm which we propose to
adapt with extra reasoning mechanisms and for which we provide a generic
open-source implementation in Rust \cite{gillard:20:ddo}. The impressive
performance of DDO triggered some theoretical research to analyze the quality of
approximate MDDs \cite{bergman:16:theoretical} and the correctness of the
relaxation operators \cite{hooker:17:job-sequencing}.

This gave rise to new lines of work. The first one focuses on the
resolution of a larger class of optimization problems; chief of which
multi-objective problems \cite{bergman:16:multiobjective} and
problems with a non-linear objective function. These are either solved by
decomposition \cite{bergman:16:multiobjective} or by using DDO to strengthen
other IP techniques \cite{davarnia:18:approx-non-linear}. A second trend aims
at hybridizing DDO with other IP techniques. For instance, by using Lagrangian
relaxation \cite{hooker:19:improved-job-sequencing} or by solving a MIP
\cite{bergman:17:optimal-bdd-relax} to derive with very tight bounds. But the
other direction is also under active investigation: for example,
\cite{tjandraatmadja:18:decision-diagrams,tjandraatmadja:19:cuts-from-relaxed}
use DD to derive tight bounds which are used to replace LP relaxation in a
cutting planes solver. Very recently, a third hybridization approach has been
proposed by Gonzàlez et al.\cite{gonzalez:20:mip-dp-ml}. It adopts the
branch-and-bound MDD perspective, but whenever an upper bound is to be derived,
it uses a trained classifier to decide whether the upper bound is to be computed
with ILP or by developing a fixed-width relaxed MDD.

The techniques (ILP-cutoff pruning and ILP-cutoff heuristic) proposed by
Gonzalez et al.\cite{gonzalez:20:mip-dp-ml} are related to RUB and LocB in the
sense that all techniques aim at reducing the search space of the problem.
However, they fundamentally differ as ILP-cutoff pruning acts as a  replacement
for the compilation of a relaxed MDD whereas the goal of RUB is to speed up the
development of that relaxed MDD by removing nodes \emph{while} the MDD is being
generated. The difference is even bigger in the case of ILP-cutoff heuristic vs
LocB: the former is used as a primal heuristic while LocB is used to filter out
sub-problems that can bear no better solution. In that sense, LocB belongs more
to the line of work started by
\cite{hooker:07:relaxed-mdd,hooker:08:propagting-separable-equalities,hooker:10:mdd-cp}:
it enforces the constraint $lb \le f(x) \le ub$ and therefore provokes the
deletion of nodes and arcs that cannot lead to the optimal solution.

More recently, Horn et al explored an idea in \cite{horn:2021} which closely
relates to RUB. They use ``fast-to-compute dual bounds'' as an admissible
heuristic to guide the compilation of MDDs in an A* fashion for the
prize-collecting TSP. It prunes portions of the state space during the MDD
construction, similarly to when RUB is applied. Our approach differs from that
of \cite{horn:2021} in that we attempt to incorporate problem specific knowledge
in a framework that is otherwise fully generic. More precisely, it is perceived
here as a problem-specific pruning that exploits the combinatorial structure
implied by the state variables. It is independent of other MDD compilation
techniques, e.g., our techniques are compatible with node merge ($\oplus$)
operators and other methodologies defined in the DDO literature. We also
emphasize that, as opposed to more complex LP-based heuristics that are now
typical in A* search, we investigate quick methodologies that are also easy to
incorporate in a MDD branch and bound.

\section*{Conclusion and future work}\label{sec:conclusion}
This paper presented and evaluated the impact of the local bound and rough upper
bound techniques to strengthen the pruning of the branch-and-bound MDD
algorithm. Our experimental study on MISP, MCP, MAX2SAT and TSPTW confirmed the
relevance of these techniques. In particular, our experiments have shown that
devising a fast and simple rough upper bound is worth the effort as it can
significantly boost the efficiency of a solver. Similarly, our experiments
showed that the use of local bound can significantly improve the efficiency of
DDO solver despite its problem agnosticism. Furthermore, it revealed that a
combination of RUB and LocB supersedes the benefit of any single reasoning
technique.
These results are very promising and we believe that the public
availability of an open source DDO framework implementing RUB and LocB might
serve as a basis for novel DP formulation for classic problems.

\newpage
\bibliographystyle{splncs04}