Tensor Completion for Radio Map Reconstruction using Low Rank and Smoothness

TL;DR

This paper introduces two novel algorithms for reconstructing radio maps from limited measurements by leveraging low tensor rank and smoothness, demonstrating superior performance over existing methods in synthetic tests.

Contribution

The paper proposes two new algorithms that combine low-rank tensor completion with smoothness constraints for radio map reconstruction, improving accuracy over prior techniques.

Findings

Algorithms outperform state-of-the-art methods in synthetic experiments.

Reconstruction accuracy improves with combined low-rank and smoothness constraints.

Synthetic ray tracing data validates the effectiveness of the proposed methods.

Abstract

Radio maps are important enablers for many applications in wireless networks, ranging from network planning and optimization to fingerprint based localization. Sampling the complete map is prohibitively expensive in practice, so methods for reconstructing the complete map from a subset of measurements are increasingly gaining attention in the literature. In this paper, we propose two algorithms for this purpose, which build on existing approaches that aim at minimizing the tensor rank while additionally enforcing smoothness of the radio map. Experimental results with synthetic measurements derived via ray tracing show that our algorithms outperform state of the art techniques.

Figures

Equations

$$x_{n+1}:=x_n+t_n\left({\mathrm{prox}}_{\gamma f}\left(2{\mathrm{prox}}_{\gamma g}{x}_n-x_n\right)-{\mathrm{prox}}_{\gamma g}{x}_n\right),$$

$${\mathrm{prox}}_{\gamma f}:\mathcal{H}\to \mathcal{H}:x\mapsto \underset{y\in \mathcal{H}}{\mathrm{argmin}}\left(f(y)+\frac{1}{2\gamma }{\Vert x-y\Vert }_{\mathcal{H}}^2\right)$$

$$\underset{\mathrm{X}\in \mathcal{T}}{\mathrm{minimize}}\sum _{i=1}^M{f}_i(\mathrm{X})$$

$$\underset{\mathrm{Z}\in {\mathcal{H}}_0}{\mathrm{minimize}}f(\mathrm{Z})+g(\mathrm{Z}),$$

$$f$$

$$g$$

$${\iota }_D(\mathrm{Z})=\{\begin{array}{ll}0,& \text{if }{\mathrm{Z}}_1=\cdots ={\mathrm{Z}}_M\\ +\infty & \text{otherwise}.\end{array}$$

$${\mathrm{prox}}_{\gamma f}\mathrm{Z}$$

$${\mathrm{prox}}_{\gamma g}\mathrm{Z}$$

$$\underset{\mathrm{X}\in \mathcal{T}}{\mathrm{minimize}}$$

$$f_i(\mathrm{X})$$

$${\mathrm{prox}}_{\gamma f_i}\mathrm{X}=\mathrm{refold}(\mathrm{shrink}(X_{(i)},\gamma ))\text{for }i=1,\dots ,N$$

$$({\mathrm{prox}}_{\gamma f_{N+1}}\mathrm{X}{)}_j=\{\begin{array}{ll}\frac{\gamma }{\lambda \gamma +1}(\lambda {\mathcal{A}}^{\ast }(\mathrm{b})+\frac{1}{\gamma }\mathrm{X}{)}_j& \text{if }j\in \Omega \\ {\mathrm{X}}_j& \text{otherwise.}\end{array}$$

$$\begin{array}{rl}\underset{\mathrm{X}\in \mathcal{T}}{\mathrm{minimize}}& \sum _{i=1}^N{\alpha }_i{\mathrm{Tv}}_2(X_{(i)})+\sum _{i=1}^N{\Vert X_{(i)}\Vert }_{\ast }\\ & +\frac{\lambda }{2}{\Vert \mathcal{A}(\mathrm{X})-\mathrm{b}\Vert }_2^2,\end{array}$$

$${\mathrm{Tv}}_2:{\mathbb{R}}^{J\times K}\to \mathbb{R}:X\mapsto \sum _{j=1}^J\sum _{k=1}^{K-1}{\left(x_{j,k+1}-x_{j,k}\right)}^2,$$

$$\underset{\mathrm{Z}\in {\mathcal{H}}_0}{\mathrm{minimize}}f(\mathrm{Z})+g(\mathrm{Z}),$$

$$f_i(\mathrm{X})$$

$${\mathrm{prox}}_{\gamma f_i}:\mathrm{X}\mapsto \underset{\mathrm{Y}\in \mathcal{T}}{\mathrm{argmin}}\left({\alpha }_i{\mathrm{Tv}}_2(Y_{(i)})+\frac{1}{2\gamma }{\Vert \mathrm{X}-\mathrm{Y}\Vert }^2\right).$$

$${\left({\mathrm{prox}}_{\gamma f_i}\mathrm{X}\right)}_{j_1\dots j_{i-1}:j_{i+1}\dots j_N}={\mathrm{prox}}_{\gamma {\tilde{f}}_i}\left({\mathrm{x}}_{j_1\dots j_{i-1}:j_{i+1}\dots j_N}\right),$$

$${\mathrm{prox}}_{\gamma {\tilde{f}}_i}$$

$$h_{\mathrm{x}}(\mathrm{y})$$

$$\frac{\partial h_{\mathrm{x}}(\mathrm{y})}{\partial y_j}$$

$$\frac{\partial h_{\mathrm{x}}(\mathrm{y})}{\partial y_1}$$

$$\frac{\partial h_{\mathrm{x}}(\mathrm{y})}{\partial y_K}$$

$$A\mathrm{y}=\frac{1}{\gamma }\mathrm{x},$$

$$A=\left[\begin{array}{ccccc}(4{\alpha }_i+\frac{1}{\gamma })& -2{\alpha }_i& 0& \dots & 0\\ -2{\alpha }_i& (4{\alpha }_i+\frac{1}{\gamma })& -2{\alpha }_i& \dots & 0\\ 0& -2{\alpha }_i& (4{\alpha }_i+\frac{1}{\gamma })& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & (4{\alpha }_i+\frac{1}{\gamma })\end{array}\right]$$

$$\mathrm{y}$$

$${\mathrm{prox}}_{\gamma f_i}:\mathrm{X}\mapsto \mathrm{X}{\times }_i\frac{1}{\gamma }{A}^{-1}.$$

$$\begin{array}{rl}\underset{\mathrm{X}\in \mathcal{T}}{\mathrm{minimize}}& \sum _{i=1}^N{\alpha }_i{\mathrm{Tv}}_1(X_{(i)})+\sum _{i=1}^N{\Vert X_{(i)}\Vert }_{\ast }\\ & +\frac{\lambda }{2}{\Vert \mathcal{A}(\mathrm{X})-\mathrm{b}\Vert }_2^2,\end{array}$$

$${\mathrm{Tv}}_1(X):=\sum _{j=1}^J\sum _{k=1}^{K-1}\left|x_{j,k+1}-x_{j,k}\right|.$$

$$f_i(\mathrm{X})$$

$${\text{NMSE}}_{\text{dB}}=10{\log }_{10}\left(\frac{\sum _{j\in {\Omega }_{\mathrm{m}}}{\left({\hat{x}}_j-x_j\right)}^2}{\sum _{j\in {\Omega }_{\mathrm{m}}}{x}_j^2}\right),$$

Full text

©2012 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

I Introduction

Accurate and reliable characterization of radio frequency environments is a key ingredient in the design and operation of modern radio networks. Knowledge of the spatial distribution of the path loss is essential for many applications in wireless networks, particularly those expected in future 5G systems. One common application in which accurate estimation of path loss maps is required is proactive resource allocation, where the knowledge of future propagation conditions along a user’s trajectory is utilized to allocate resources[1]. Path loss maps can also be used for fingerprint-based localization schemes, where the geographical locations of users are estimated based on their measured channel states [2, 3].

For the above reasons, path loss estimation has gained a great deal of attention from both the academia and the industry, but accurate path loss estimation in wireless networks remains a challenging task to date [4]. A plethora of theoretical and empirical path loss models have been proposed by researchers over the years, and, although some have been proved to be useful in selected applications, many existing models fail to deliver good predictive performance in real world settings [5]. As a result, data-driven methods based on learning tools are increasingly gaining attention, and they have been shown to produce results with higher reliability and accuracy than purely model-driven methods[6].

Many of these data-driven methods for radio map reconstruction utilize matrix completion algorithms[2, 3, 7, 8], and, in the field of image and video processing, extensions of these matrix reconstruction methods to higher-order tensors have been shown to provide superior performance[9, 10]. However, tensor methods have rarely been applied to radio map reconstruction problems.

Against this background, we propose novel approaches for radio map reconstruction based on tensor reconstruction algorithms. One of the main advantages of the proposed algorithms is that they can easily take into account information such as path loss in multiple heights, frequency bands, and towards multiple base stations, among others. In addition, the proposed algorithms exploit the fact that radio maps are typically smooth to increase the quality of the reconstruction. We evaluate the performance of the proposed algorithms by using synthetic measurements that were obtained from ray tracing simulations performed during the METIS2020 project[11].

The remainder of this paper is organized as follows. In Section II we introduce standard results in mathematics that are crucial to the algorithms we propose. In Section III we derive novel algorithms to solve the aforementioned problems. Numerical experiments are given in Section IV.

Notation: Scalars, column vectors, matrices, and tensors are denoted by italic lower case letters $x$, bold lower case letters ${\bm{\mathrm{x}}}$, italic upper case letters $X$, and bold upper case letters ${\bm{\mathrm{X}}}$, respectively. The element $(i_{1},i_{2},\dots,i_{N})$ of a tensor ${\bm{\mathrm{X}}}\in\mathbb{R}^{n_{1}\times\dots\times n_{N}}$ is denoted by $x_{i_{1}i_{2}\dots i_{N}}$. The vector space containing $N$th order tensors is denoted by $\mathcal{T}:=\mathbb{R}^{n_{1}\times\dots\times n_{N}}$. By equipping this vector space with the inner product $\left\langle\mkern 1.0mu{\bm{\mathrm{X}}},{\bm{\mathrm{Y}}}\mkern 1.0mu\right\rangle:=\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\dots\sum_{i_{N}=1}^{n_{N}}x_{i_{1}i_{2}\dots i_{N}}y_{i_{1}i_{2}\dots i_{N}}$ and its induced norm $\left\|{\bm{\mathrm{X}}}\right\|:=\sqrt{\left\langle\mkern 1.0mu{\bm{\mathrm{X}}},{\bm{\mathrm{X}}}\mkern 1.0mu\right\rangle}$ we obtain a Hilbert space of tensors. The high-order generalizations of columns and rows of a matrix are called fibers, which are extracted by fixing all coordinates except for one. In particular, mode-$m$ fibers are denoted by ${\bm{\mathrm{x}}}_{i_{1}\dots i_{m-1}:i_{m+1}\dots i_{N}}$, where a colon in the index is used to indicate all elements of a mode. To compute the mode-$m$ unfolding of a tensor, we rearrange all mode-$m$ fibers of the tensor as the columns of a matrix. More formally, $X_{(m)}\in\mathbb{R}^{n_{m}\times I_{m}}$ denotes the mode-$m$ unfolding of the tensor ${\bm{\mathrm{X}}}$, where $I_{m}=\prod_{\begin{subarray}{c}k=1\\ k\neq m\end{subarray}}^{N}n_{k}$ and the tensor element $(i_{1},\dots,i_{N})$ is mapped to the matrix element $(i_{m},j)$, with $j=1+\sum_{\begin{subarray}{c}k=1\\ k\neq m\end{subarray}}^{N}(i_{k}-1)J_{k}$ and $J_{k}=\prod_{\begin{subarray}{c}\ell=1\\ \ell\neq m\end{subarray}}^{k-1}n_{\ell}$. The $m$-mode multiplication of a tensor ${\bm{\mathrm{X}}}\in\mathbb{R}^{n_{1}\times\dots\times n_{N}}$ and a matrix $Y\in\mathbb{R}^{J\times n_{m}}$ is denoted by ${\bm{\mathrm{X}}}\times_{m}Y,$ and it is defined by the following relation based on the unfoldings: ${\bm{\mathrm{Z}}}={\bm{\mathrm{X}}}\times_{m}Y\ \Leftrightarrow\ Z_{(m)}=YX_{(m)}$. The nuclear norm of a matrix $X$ is defined as $\left\|X\right\|_{*}:=\sum_{i}\sigma_{i}(X)$, where $\sigma_{i}(X)$ denotes the $i$th largest singular value of $X$. The class of all lower semi-continuous convex functions from an arbitrary real Hilbert space $\mathcal{H}$ to $(-\infty,\infty]$ that are not identically equal to $+\infty$ are denoted by $\Gamma_{0}(\mathcal{H})$.

II Preliminaries

The proposed algorithms are heavily based on the Douglas-Rachford splitting method, so we briefly review this iterative method in Section II-A. Then, in Section II-B, we present the algorithm introduced in [12], which is the algorithm we improve in Section III by considering specific features of radio maps.

II-A Douglas-Rachford splitting

Given an arbitrary Hilbert space $\mathcal{H}$, the objective of the Douglas-Rachford splitting algorithm is to minimize the sum of two functions $f,g\in\Gamma_{0}(\mathcal{H})$ satisfying some additional very mild assumptions [13, Corollary 28.3] that are valid in the applications described in the next sections. This algorithm generates a sequence $(x_{n})_{n\in\mathbb{N}}$ given by

[TABLE]

where $t_{n}\in[0,2]$ is the step size and

[TABLE]

denotes the proximal map of a function $f\in\Gamma_{0}(\mathcal{H})$ with regularization parameter $\gamma\in(0,\infty)$[13, p. 211]. For convergence properties of the recursion in (1), we refer the readers to [13, Section 28.3].

Although the Douglas-Rachford method as shown above is used to minimize the sum of two functions in $\Gamma_{0}(\mathcal{H})$, we can straightforwardly modify the algorithm to incorporate more than two functions[12]. Consider for concreteness the following optimization problem in the tensor space $\mathcal{T}$:

[TABLE]

where $f_{i}\in\Gamma_{0}(\mathcal{T})$ for each $i\in\left\{1,\dots,M\right\}$. By defining the Hilbert space $\mathcal{H}_{0}:=\mathcal{T}\times\mathcal{T}\times\dots\times\mathcal{T}$ ($M$-fold Cartesian product) with elements ${\bm{\mathrm{X}}}=({\bm{\mathrm{X}}}_{1},\dots,{\bm{\mathrm{X}}}_{M})$ and inner product $\left\langle\mkern 1.0mu{\bm{\mathrm{X}}},{\bm{\mathrm{Y}}}\mkern 1.0mu\right\rangle_{\mathcal{H}_{0}}:=\frac{1}{M}\sum_{i=1}^{M}\left\langle\mkern 1.0mu{\bm{\mathrm{X}}}_{i},{\bm{\mathrm{Y}}}_{i}\mkern 1.0mu\right\rangle$, we can obtain a minimizer of (3) by solving

[TABLE]

where $f$ and $g$ are given by

[TABLE]

where

[TABLE]

We assume that the conditions for the convergence of the Douglas-Rachford splitting method are satisfied and refer the readers to [13, Corollary 28.3] for more details. Note that ${\bm{\mathrm{X}}}^{\star}$ is a solution to (3) if and only if $({\bm{\mathrm{Z}}}_{1}^{\star},\dots,{\bm{\mathrm{Z}}}_{M}^{\star})$ is a solution to (4) and ${\bm{\mathrm{X}}}^{\star}={\bm{\mathrm{Z}}}_{1}^{\star}=\dots={\bm{\mathrm{Z}}}_{M}^{\star}$.

To solve (4) using the iteration in (1), we can use the proximal maps of $f$ and $g$ given by [12]

[TABLE]

where $\operatorname{mean}({\bm{\mathrm{Z}}}):=\frac{1}{M}\sum_{i=1}^{M}{\bm{\mathrm{Z}}}_{i}$.

II-B Convex approximation of the rank minimization problem

Suppose we have the problem of completing an $N$th-order tensor ${\bm{\mathrm{X}}}\in\mathcal{T}$ of low rank from a subset of $p$ given entries, which we express by $\mathcal{A}({\bm{\mathrm{X}}})={\bm{\mathrm{b}}}$, where $\mathcal{A}:\mathcal{T}\to\mathbb{R}^{p}$ is a linear map that maps the values at positions of the tensor for which samples are available (indicated by the set $\Omega\subset\left\{1,\dots,n_{1}\right\}\times\dots\times\left\{1,\dots,n_{N}\right\}$) into a vector and ${\bm{\mathrm{b}}}\in\mathbb{R}^{p}$ is the vector of known samples with $p<\prod_{i=1}^{N}n_{i}$.

To estimate the unknown samples, the authors of [12] propose to solve the following minimization problem, which considers the sums of the nuclear norms of the unfoldings as a convex approximation for the tensor rank:

[TABLE]

where $\lambda$ is a regularization parameter (we will give more detail on choosing $\lambda$ in Section IV).

In order to apply the techniques outlined in Section II-A to this problem, we first have to compute the proximal maps for the following family of functions:

[TABLE]

The proximal maps of these functions have already been derived in [12]:

[TABLE]

Here, $\operatorname{refold}\left(\hskip 1.19995pt\cdot\hskip 1.19995pt\right)$ denotes the folding of the matrix into a tensor; i.e. the inverse operation of the unfolding. The shrinkage operator $\operatorname{shrink}\left(\hskip 1.19995pt\cdot\hskip 1.19995pt\right)$ performs singular value soft-thresholding on a matrix. First, the singular value decomposition $X=U\Sigma V^{*}$ is computed, where $\Sigma=\operatorname{diag}(\sigma_{1}(X),\dots,\sigma_{r}(X))$ is the diagonal matrix containing the singular values of $X$. Then, $\tilde{\Sigma}=\operatorname{diag}(\max(\sigma_{1}(X)-\tau,0),\dots,\max(\sigma_{r}(X)-\tau,0))$ is the diagonal matrix containing the singular values shrunk by $\tau$. The result of the shrinkage operator is then given by $\operatorname{shrink}(X,\tau):=U\tilde{\Sigma}V^{*}$. The operator $\mathcal{A}^{*}$ is the adjoint operator of $\mathcal{A}$, which maps the known entries in the vector ${\bm{\mathrm{b}}}$ to their corresponding position in the tensor, while setting all unknown entries to zero.

III Proposed Algorithms

When representing a radio map as a matrix, we can expect regular structures in the matrix inherited from buildings, roads, and many other items. As a result, it is reasonable to assume that rows and columns are in general linearly dependent, making such a matrix rank-deficient. Therefore, when stacking many matrices representing correlated physical phenomena (e.g. path loss in different heights), we can also expect the resulting tensor to have low rank. This low-rank property can be exploited for the tensor completion problem by applying the Douglas-Rachford algorithm given by (1) to the problem presented in Section II-B.

This approach, however, has severe limitations, since the tensor can be exactly reconstructed only if the number of available samples is large. This is consistent with [12], where the authors only provide results for a relatively large number of samples. For the radio map reconstruction problem at hand, the given number of samples is typically much lower than required by the algorithm presented in Section II-B. In such cases, the algorithm in Section II-B produces a (low-rank) solution that is in general no longer consistent with the tensor being estimated. Severe artifacts might appear because fibers with a very small number of samples, or even with no samples at all, can be replaced by arbitrary copies of other fibers without increasing the rank. These solutions, however, exhibit a very low spatial coherence, which is not a typical characteristic of radio maps[14], as can be seen in the measurements that are shown later in Fig. 1 in Section IV.

Therefore, in the next sections, we propose two novel algorithms to address these limitations. First, we add the $\ell_{2}$ norm of the total variation as regularizer to (10) in order to ensure that neighboring tensor elements have similar values. The second algorithm instead is based on the $\ell_{1}$ regularization of the total variation in order to reduce the number of discontinuities in the path loss map.

III-A L2 Norm Total Variation Regularization

To take into account that path loss maps have high spatial correlation, we add additional terms to the cost function in (10) that penalize abrupt changes in neighboring components of the tensor. In particular, in this section we use the $\ell_{2}$ norm of the total variation (Tv) as the penalty term in the first proposed approach. More precisely, the first approach aims at solving the following problem:

[TABLE]

where $\alpha_{i}$ and $\lambda$ are regularization parameters, and the $\operatorname{Tv}_{2}$ operator computes the sum of the squared differences between neighboring elements of each mode-$i$ fiber:

[TABLE]

and where $K=n_{i}$ and $J=\prod_{\begin{subarray}{c}k=1\\ k\neq i\end{subarray}}^{N}n_{k}$. A scheme for choosing good regularization parameters will be presented in Section IV.

To apply the algorithm in Section II-A to solve (14), we first rewrite the problem in (14) as

[TABLE]

where the functions $f$ and $g$ follow the definitions in (5) and (6), and where the $M=2N+1$ functions are given by

[TABLE]

Using the iterations in (1) to solve (16) requires the proximal maps of $f$ and $g$, which are given in (8) and (9). This in turn requires the proximal maps of $f_{i}$ for each $i\in\{1,\dots,2N+1\}$. The proximal maps of $f_{N+1},\dots,f_{2N+1}$ are given in (12) and (13), so it remains to derive the proximal maps of $f_{1},\dots,f_{N}$, which, by definition, are given by

[TABLE]

The $\operatorname{Tv}_{2}$ operator is an additively separable function w.r.t. the fibers of the tensor. Therefore, the proximal map can be computed for each fiber independently [15]:

[TABLE]

where $\tilde{f}_{i}({\bm{\mathrm{x}}})=\alpha_{i}\sum_{k=1}^{K-1}\left(x_{k+1}-x_{k}\right)^{2}$ is the penalty term for a single fiber. In turn, the proximal map for $\tilde{f}_{i}({\bm{\mathrm{x}}})$ is given by

[TABLE]

where

[TABLE]

We now proceed to derive the proximal map in (20) in closed form. To this end, we first derive the partial derivatives of $h_{\bm{\mathrm{x}}}({\bm{\mathrm{y}}})$, which are given in (22), with special cases given in (23) and (24):

[TABLE]

Since (21) is convex and the problem is unconstrained, we can find a global minimizer by setting all partial derivates to zero to obtain the following system of linear equations:

[TABLE]

where

[TABLE]

The above system of linear equations can be solved efficiently because $A$ is a tridiagonal Toeplitz matrix, which can be efficiently inverted[16], so the following vector is easy to obtain in practice:

[TABLE]

The result of applying this equation to all fibers of the tensor can be written more compactly as the $i$-mode multiplication of the tensor with $A^{-1}$:

[TABLE]

We have now derived the proximal map of $f_{i}$ for each $i\in\{1,\dots,2N+1\}$ in (12), (13), and (27), which, together with (8) and (9), enables us to apply the Douglas-Rachford splitting method to solve (16).

As we show later by simulations, since this algorithm uses additional information, it is able to reconstruct the tensor accurately with fewer samples compared to the algorithm presented in Section II-B. One potential limitation of the algorithm presented in this section is that transitions provoked by, for example, building edges, might be smoothed heavily due to the quadratic penalty of the $\ell_{2}$ norm.

Remark 1**.**

In numerical experiments, we have noticed that the magnitude of the edges of the tensor decreases in each iteration. To alleviate this effect, we propose, as a heuristic, to scale each row of the matrix $A^{-1}$ to have a sum of 1.

III-B L1 Norm Total Variation Regularization

To overcome the potential limitation of the algorithms introduced in the previous subsection, we now present an algorithm that uses an $\ell_{1}$ regularizer in place of the $\ell_{2}$ regularizer. Intuitively, by using the $\ell_{1}$ norm, we minimize the number of sharp transitions between neighboring elements of the tensor unfoldings. Formally, the $\ell_{1}$ norm total variation regularization problem is given as follows:

[TABLE]

where $\alpha_{i}$ and $\lambda$ are regularization parameters, and

[TABLE]

To solve (28), we can again use the algorithm presented in Section II-A with the following definition:

[TABLE]

The challenge in solving (28) with the Douglas-Rachford splitting method is the computation of the proximal maps for $f_{1},\dots,f_{N}$. The method presented in Section III-A is not applicable here since the absolute value function is non-differentiable. Typically, in these cases the proximal maps are computed with iterative methods such as those described in [17], which is the approach we use in the simulations. However, some authors have also proposed to use neural networks as fast approximation of proximal maps[18].

IV Numerical Experiments and Conclusion

To evaluate the performance of the algorithms, we used synthetic path loss measurements obtained by ray tracing simulations during the METIS2020 project[11]. The chosen scenario (Madrid scenario) represents an urban environment with buildings of different heights, parks and streets, all arranged on a regular grid. The path loss measurements for a height of $1.5\text{\,}\mathrm{m}$ are shown in Fig. 1, but measurements for $26\text{\,}\mathrm{m}$ and $43.5\text{\,}\mathrm{m}$ were also available.

The available measurements111available at https://metis2020.com/wp-content/uploads/simulations/Ray-Tracing-files-for-TC2-Part-I.zip contain three repetitions of the environment model in $x$ and $y$ direction to avoid border effects. For the simulations, however, the repetitions were removed, and only the central block was used. The resolution of the measurements was $3\text{\,}\mathrm{m}$, and the area had a size of 387\text{\,}\mathrm{m}$\times$552\text{\,}\mathrm{m}, resulting in a tensor of dimensions $129\times 184\times 3$.

For the simulations, we randomly selected a given percentage of tensor elements, and we used these as input to the algorithms. As the performance metric, we used the normalized mean squared error (NMSE), which is defined as

[TABLE]

where $\hat{{\bm{\mathrm{X}}}}$ and ${\bm{\mathrm{X}}}$ are reconstructed and correct tensors given in $\mathrm{d}\mathrm{B}$, respectively, and $\Omega_{\mathrm{m}}\subseteq\left\{1,\dots,n_{1}\right\}\times\dots\times\left\{1,\dots,n_{N}\right\}$ is the set of indizes that were not used as input to the algorithms.

The regularization parameters $(\alpha_{i})_{i\in\left\{1,\dots,N\right\}}$ were obtained by using cross-validation with the holdout method. More specifically, $25\text{\,}\%$ of the given tensor elements were removed from the measurements at random and then used as test set. We then selected a set of reasonable regularization parameter values and applied cross-validation to determine the best parameters out of all possible combinations. For given regularization parameters $\alpha_{i}$ we applied the continuation technique described in [12] to $\lambda$ to achieve a sufficiently small value for $\left\|\mathcal{A}({\bm{\mathrm{X}}})-{\bm{\mathrm{b}}}\right\|_{2}$. To this end, we ran the algorithm for a small value for $\lambda$ until no significant progress was achieved, then increased $\lambda$ and ran the algorithm again. We repeated this process until either a further increase of $\lambda$ had no influence on the reconstruction result or a given number of iterations was achieved.

For comparison, we implemented the radial basis function (RBF) algorithm, as described in [19, Sect. 5.1]. We achieved the best results with the multiquadric function $\phi(r)=\sqrt{1+\left(\frac{r}{\epsilon}\right)^{2}}$, where the hyperparameter $\epsilon$ was chosen by line search, using the same cross-validation procedure as described for the other algorithms.

As can be seen in Fig. 2, the algorithm presented in Section III-B (denoted as L1-TV + Rank) yields a constant improvement over the existing algorithm presented in Section II-B (denoted as Rank). On the other hand, the algorithm presented in Section III-A with the heuristic mentioned in Remark 1 (denoted as L2-TV + Rank (w heur)), although not particularly good for a large number of samples, provides a significant improvement for a smaller number of samples. Only when less than $2\text{\,}\mathrm{\char 37\relax}$ of all tensor elements are given for reconstruction is it slightly outperformed by the RBF algorithm.

The fact that the proposed algorithms greatly outperform state-of-the-art techniques indicates that the proposed algorithms are successfully exploiting the assumed structure of radio maps. The comparison to the RBF algorithm shows that the assumption of smoothness is not sufficient and exploiting the low-rank property is necessary to achieve good performance in a wide range of conditions. The proposed algorithms form a framework which is very flexible and can be applied to a wide range of problems in different areas.