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

One of the key goals of the ultrarelativistic nuclear collisions is to gain information on the hadron-parton phase diagram, which is characterized by different phases and different types of the phase transitions \cite{banks1983deconfining}. Quantum Chromodynamics (QCD), the gauge field theory that describes the strong interactions of colored quarks and gluons and their colorless bound states, has two important intensive state parameters at equilibrium, namely temperature $T$ and baryon chemical potential $\mu_{\mathtt{b}}$. A remarkable world-wide theoretical and experimental effort has been dedicated to the study of strongly interacting matter under extreme condition of temperature and baryon chemical potential. The lattice QCD simulations provide an \textit{a priori} non-perturbative regularization of QCD that makes it compliant with analytic and computational methods with no model assumptions other than QCD itself being needed to formulate the theory. The temperature and density (chemical potential) dependence of the bulk thermodynamic quantities, commonly summarized as the equation of state (EoS), provides the most basic characterization of equilibrium properties of the strongly interacting matter. Its analysis within the framework of lattice QCD has been refined ever since the early calculations performed in pure $SU(N)$ gauge theories \cite{engels1981high}. The EoS at vanishing chemical potentials does already provide important input into the modeling of the hydrodynamic evolution of hot and dense matter created in heavy-ion collisions \cite{Karsch:2003zq,Tawfik:2004sw}. While this is appropriate for the thermal conditions met in these collisions at the LHC and the highest RHIC beam energies, knowledge of the EoS at non-vanishing baryon, strangeness and electric charge chemical potentials is indispensable for the hydrodynamic models of the conditions met in the beam energy scan (BES) at RHIC \cite{Gyulassy:2004zy} and in future experiments at facilities like FAIR at GSI and NICA at JINR \cite{Tawfik:2016cot,Tawfik:2016wfm}.

Bulk thermodynamic observables such as pressure, energy density, and entropy density as well as the second order quantities such as the specific heat and velocity of sound have now been obtained at vanishing chemical potentials for the three lightest quark flavors \cite{Tawfik:2004vv}. By the analysis of the chiral transition temperature, $T_{\mathtt{c}}\simeq 154 \pm 9$ MeV \cite{bazavov2012chiral} it has also been shown that the bulk thermodynamic observables change smoothly in the transition region probed. Due to the well-known sign problem encountered in the lattice QCD formulations at finite chemical potential, a direct calculation of the EoS at non-zero chemical potential is unfortunately not fully reliable. A lot of efforts has been made to circumvent the divergences at non-zero chemical potential, such as Taylor expansion of the thermodynamic potential on course lattices besides other sophisticated computational techniques \cite{gavai2001quark,allton2002qcd,ejiri2006isentropic,borsanyi2012qcd,gunther2017qcd,d2017higher} which made it possible to conduct calculations covering the range $0 \leq \mu_{\mathtt{b}}/T \lesssim 3$ that is expected to be explored with the BES at RHIC by varying the beam energies in the range $7.7\leq \sqrt{s_{\mathtt{NN}}} \leq 200~$GeV \cite{Bzdak:2019pkr}. A promising approach in this quest is the investigation of hadron production. The hadron resonance gas (HRG) is customarily used in the lattice QCD calculations as a reference for  the hadronic sector \cite{Huovinen:2009yb,Borsanyi:2010bp}. At low temperatures, they are found to be in quite good agreement with the HRG model calculations \cite{Tawfik:2014eba}, although some systematic deviations have been observed, which may be attributed to the existence of additional resonances which are not taken into account in HRG model calculations based on well established resonances listed by the particle data group \cite{Tanabashi:2018oca} and perhaps to need to extend the model to incorporate interactions.

In the HRG model, the thermodynamics of a strongly interacting system is conjectured to be approximated as an ideal gas composed of hadron resonances with masses $\lesssim2~$GeV \cite{Tawfik:2014eba,Karsch:2003zq} that are treated as a free gas, exclusively in the hadronic phase, i.e. below $T_c$. Therefore, the hadronic phase in the confined phase of QCD could be modeled as a non-interacting gas of the hadron resonances. It is reported in recent literature that the standard performance of the HRG model seems to be unable to describe all the available data that is predicted by recent lattice QCD simulations \cite{samanta2019exploring,vovchenko2017multicomponent}. The conjecture to incorporate various types of interactions has been worked out in various studies \cite{Tawfik:2004sw,Huovinen:2017ogf,Vovchenko:2017xad,Vovchenko:2017drx}. When comparing the thermodynamics calculated within the HRG framework with the corresponding data obtained using lattice QCD methods, one has to decide how to incorporate interactions among the hadrons.

Arguments based on the S-matrix approach \cite{dashen1969s,venugopalan1992thermal,lo2017s} suggest that the HRG model includes attractive interactions between hadrons which lead to the formation of resonances. More realistic hadronic models take into account the contribution of both attractive and repulsive interactions between the component hadrons. Repulsive interactions in the HRG model had previously been considered in the framework of the relativistic cluster and virial expansions \cite{venugopalan1992thermal}, via repulsive mean fields \cite{olive1981thermodynamics,olive1982quark}, and via excluded volume (EV) corrections \cite{hagedorn1980hot,gorenstein1981phase,hagedorn1983pressure,kapusta1982thermodynamics,
rischke1991excluded,anchishkin1995generalization}. In particular, the effects of EV interactions between hadrons on HRG thermodynamics \cite{satarov2009equation,andronic2012interacting,bhattacharyya2014fluctuations,
albright2014matching,vovchenko2015hadron,albright2015baryon,redlich2016thermodynamics,
alba2017excluded} and on observables in heavy-ion collisions \cite{braun1999chemical,cleymans2006comparison,begun2013hadron,fu2013higher,
vovchenko2017examination,alba2018flavor,satarov2017new,vovchenko2017van} has extensively been studied in the literature. Recently, repulsive interactions have received renewed interest in the
context of lattice QCD data on fluctuations of conserved charges. It was shown that large deviations of several fluctuation observables from the ideal HRG baseline could well be interpreted
in terms of repulsive baryon-baryon interactions \cite{huovinen2018hadron,vovchenko2017van,vovchenko2017equations}.

The present script is organized as follows: In Section \label{sec:Mod}, we review the detailed formalism of the conventional ideal (uncorrelated) HRG model, then we develop a non-ideal (correlated) statistical correction to the ideal HRG model inspired by Beth-Uhlenbeck (BU) quantum theory of non-ideal gases. The calculations of the HRG thermodynamics based on the proposed correction are discussed in Section \label{sec:Res}. Section \label{sec:Cncls} is devoted to the conclusions and outlook.

\section{Model Description}
\label{sec:Mod}

In this study, we use the particle interaction probability term originally implemented in the expression for the second virial coefficient worked out in ref. \cite{uhlenbeck1937quantum} in order to suggest a statistical correction to the uncorrelated HRG model. Beth and Uhlenbeck suggested a connection between the virial coefficients and the probabilities of finding pairs, triples and so on, of particles near each other \cite{uhlenbeck1936quantum}. In the classical limit, which is usually designated by sufficiently high temperatures and/or low particle densities, it was shown that these probabilities (explicit expressions are to follow next section) can be expressed by Boltzmann factors so long as the de Broglie wavelength, which is a common measure of the significance of the quantum non-localiy, is small enough compared with the particle spacial extent measured by the particle {''diameter''} \cite{uhlenbeck1936quantum}. Such a particle diameter can be considered as a measure of the spatial extent within which a particle can undergo hardcore (classical) interactions. Based on a comparison of model with experimental results, it was concluded  that at sufficiently low temperatures for which the thermal de Broglie wave length is comparable with the particle diameter, deviations from the classical excluded volume model due quantum effects will be significant \cite{uhlenbeck1936quantum}.

An extension has been made to the quantum mechanical model of the particle interactions as proposed in ref. \cite{uhlenbeck1936quantum} by considering the influence of Bose or Fermi statistics in addition to the effect of the inclusion of discrete quantum states for a general interaction potential that is not necessarily central \cite{uhlenbeck1937quantum}. The expression for the second virial expansion developed in ref. \cite{uhlenbeck1936quantum} and then extended in ref. \cite{uhlenbeck1937quantum} was later generalized using the cluster integral to describe particle interactions provided that those particles don't form bound states \cite{dashen1969s,venugopalan1992thermal,kostyuk2001second}.

The quantum mechanical Beth-Uhlenbeck (BU) approach were quite recently used to model the repulsive interactions between baryons in a hadron gas \cite{vovchenko2018beth}. The second virial coefficient or the “excluded volume” parameter was calculated within the BU approach and found to be temperature dependent, and found also to differ dramatically from the classical excluded volume (EV) model result. At temperatures $T=100-200$ MeV, the widely used classical EV model \cite{huang1963statistical,landau1980lifshitz,greiner2012thermodynamics} underestimates the EV parameter for nucleons at a given value of the nucleon hard-core radius (assumed $\simeq 0.3$ fm) by large factors of 3-4. It was thus concluded in \cite{vovchenko2018beth} that previous studies, which employed the hard-core radii of hadrons as an input into the classical EV model, have to be re-evaluated using the appropriately rescaled quantum mechanical EV parameters.

In this section, we first introduce the basic formulation of the ideal HRG model. Then we develop a statistical model inspired by Beth-Uhlenbeck (BU) quantum theory of non-ideal (correlated) gases \cite{uhlenbeck1937quantum} as a correction to the ideal (uncorrelated) HRG model. We thereby implement in our calculations a modified version of the partition function of a typical ideal gas. In the framework of bootstrap picture \cite{fast1963statistical,fast1963nuovo,eden1966polkinghorne}, an equilibrium thermal model for an interaction free gas has a partition function $Z(T, \mu, V)$ from which the thermodynamics of such a system can be deduced by taking the proper derivatives.

\subsection{Non-correlated ideal HRG}

In a grand canonical ensemble, the partition function reads \cite{Tawfik:2014eba,Karsch:2003vd,Karsch:2003zq,Redlich:2004gp,Tawfik:2004sw,Tawfik:2005qh}
\begin{eqnarray}
Z(T,V,\mu)=\mbox{Tr}\left[\exp\left(\frac{{\mu}N-H}{T}\right)\right], \label{eq:Z}
\end{eqnarray}
where $H$ is Hamiltonian combining all relevant degrees of freedom and $N$ is the number of constituents of the statistical ensemble. Eq. (\ref{eq:Z}) can be expressed as a sum over all hadron resonances taken from recent particle data group (PDG) \cite{Tanabashi:2018oca} with masses up to $2.5~$ GeV,
\begin{eqnarray}
\ln  Z(T,V,\mu)=\sum_i{{\ln Z}_i(T,V,\mu)} = V \sum_i \frac{g_i}{2{\pi}^2}\int^{\infty}_0{\pm p^2 dp {\ln} {\left[1\pm {\lambda}_i \exp\left(\frac{-{\varepsilon}_i(p)}{T} \right) \right]}}, \label{eq:lnZ}
\end{eqnarray}
where the pressure can be derived as $T\partial \ln  Z(T,V,\mu)/\partial V$, $\pm$ stands for fermions and bosons, respectively. $\varepsilon_{i}=\left(p^{2}+m_{i}^{2}\right)^{1/2}$ is the dispersion relation and $\lambda_i$ is the fugacity factor of the $i$-th particle \cite{Tawfik:2014eba},
\begin{eqnarray}
\lambda_{i} (T,\mu)=\exp\left(\frac{B_{i} \mu_{\mathtt{b}}+S_{i} \mu_{\mathtt{S}}+Q_{i} \mu_{\mathtt{Q}}}{T} \right), \label{eq:lmbd}
\end{eqnarray}
where $B_{i} (\mu_{\mathtt{b}})$, $S_{i} (\mu_{\mathtt{S}})$, and $Q_{i} (\mu_{\mathtt{Q}})$ are baryon, strangeness, and electric charge quantum numbers (their corresponding chemical potentials) of the $i$-th hadron, respectively. From phenomenological point of view, the baryon chemical potential $\mu_{\mathtt{b}}$ - along the chemical freezeout boundary, where the production of particles is conjectured to cease - can be related to the nucleon-nucleon center-of-mass energy $\sqrt{s_{\mathtt{NN}}}$ \cite{tawfik2015thermal}
\begin{eqnarray}
\mu_{\mathtt{b}} &=& \frac{a}{1+b \sqrt{s_{\mathtt{NN}}}}, \label{eq:mue}
\end{eqnarray}
where
$a=1.245\pm0.049~$GeV and $b=0.244\pm0.028~$GeV$^{-1}$. In addition to pressure, the number and energy density, respectively, and likewise the entropy density and other thermodynamics can straightforwardly be derived from the partition function by taking the proper derivatives
\begin{eqnarray}
n_i(T,\mu) &=&
\sum_{i}\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}{p^{2} dp \frac{1}{\exp\left[\frac{\mu_{i} - \varepsilon_{i}(p)}{T}\right] \pm 1}}, \\
\rho_i(T,\mu) &=&
\sum_{i}\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}{p^{2} dp \frac{-\varepsilon_{i(p)}\pm \mu_i}{\exp\left[\frac{\mu_{i}-\varepsilon_i(p)}{T}\right] \pm 1}}. \label{eq:e}
\end{eqnarray}

It should be noticed that both  $T$ and $\mu=B_{i} \mu_{\mathtt{b}}+S_{i} \mu_{\mathtt{S}}+\cdots$ are related to each other and to $\sqrt{s_{\mathtt{NN}}}$ \cite{Tawfik:2014eba}. As an overall thermal equilibrium is assumed, $\mu_{\mathtt{S}}$ is taken as a dependent variable to be estimated due to the strangeness conservation, i.e. at given $T$ and $\mu_{\mathtt{b}}$, the value assigned to $\mu_{S}$ is the one assuring $\langle n_{\mathtt{S}}\rangle-\langle n_{\bar{\mathtt{S}}}\rangle=0$. Only then, $\mu_{\mathtt{S}}$ is combined with $T$ and $\mu_{\mathtt{b}}$ in determining the thermodynamics, such as the particle number, energy, entropy, etc. The chemical potentials related to other quantum charges, such as the electric charge and the third-component isospin, etc. can also be determined as functions of $T$, $\mu_{\mathtt{b}}$, and $\mu_{\mathtt{S}}$ and each of them must fulfill the corresponding conservation laws.

\subsection{Quantum-statistically correlated HRG}

As introduced in the previous section, for a quantum gas of fermions and bosons with mass $m_{i}$ and correlation (interaction) distance $r$, at temperature $T$ and vanishing $\mu_{\mathtt{b}}$, a two-particle interaction probability of the form
\begin{eqnarray}
1\pm \exp\left(-4\pi^{2} m_{i} T r^{2}\right)
\end{eqnarray}
was first introduced by Beth and Uhlenbeck \cite{uhlenbeck1937quantum} in an attempt to model the interactions of a quantum gas of particles assuming a general potential and neglecting the possibility for bound states formation. The Boltzmann-like term $\exp\left(-4\pi^{2} {m}_{i} T r^{2}\right)$ remains in effect even for an ideal gas, which is a typical approximation at sufficiently high temperatures. The $\pm$ sign expresses the apparent attraction (repulsion) between bosons (fermions) due to change of statistics \cite{uhlenbeck1937quantum}. Inspired by such a correction, we introduce a correction for the probability term in the expression for the the ideal hadron gas partition function given in Eq. (\ref{eq:Z}). We propose a new probability term of the form
\begin{eqnarray}
1\pm \lambda_i \exp\left(\frac{-{\varepsilon}_i(p)}{T}\right) \left[1\pm \exp\left(-4\pi^{2} {m}_{i} T r^{2}\right)\right].
\end{eqnarray}

This corrected probability function obviously incorporates interactions in the the hadron resonance gas in the sense of Beth and Uhlenbeck quantum correlations \cite{uhlenbeck1937quantum} with $r$ being the  correlation (interaction) length between any two hadrons at equilibrium temperature $T$. Based on our proposed corrected probability function, we modify the non-correlated HRG partition function $Z(T, \mu, V)$ to have the following form
\begin{eqnarray}
\ln Z ^{\prime}  \textprime(T,V,\mu) =
\sum_i V \frac{g_i}{2 {\pi}^2}\int^{\infty}_0{\pm p^2 dp {\ln} \left[1 \pm \lambda_{i} \exp\left(\frac{-{\varepsilon}_{i(p)}}{T} \right)  \left[1\pm \exp\left(4 \pi^{2} m_{i} T r^{2}\right)\right]\right]},  \label{eq:lnZ/}
\end{eqnarray}
which apparently sums over all hadron resonances following the same recipe described in motivating Eq. (\ref{eq:lnZ}) for the case of non-correlated HRG. The thermodynamics of the correlated HRG can thus be calculated by taking the proper derivatives of $\ln Z \textprime$ as explicitly stated in the corresponding  non-correlated HRG case discussed above.

\section{Calculation Results}
\label{sec:Res}

We confront the data of the thermodynamics calculated using our statistically corrected HRG model based on Eq. (\ref{eq:lnZ/}) with the corresponding lattice thermodynamics data from \cite{bazavov2017qcd,bazavov2014equation} in the temperature range $T\in [130, 200~$MeV$]$. These temperatures are rather typical for the phenomenological applications in the context of heavy-ion collisions and lattice QCD equation-of-state. In refs. \cite{bazavov2017qcd,bazavov2014equation}, the authors calculated the QCD equation of state using Taylor expansions that include contributions from up to sixth order in the baryon, strangeness and electric charge chemical potentials. Calculations have been performed with a highly improved staggered quark action in the temperature range $T\in [130, 330~$MeV$]$ using up to four different sets of lattice cut-offs. The lattice data we are confronting to our model are shown on Figure 7. of ref. \cite{bazavov2017qcd}, (Left) the total pressure in (2+1)-flavor QCD and (Right) the total energy density in (2+1)-flavor QCD for several values of $\mu_{b}/T$.

  Figure \ref{fig:one} of this letter depicts the temperature dependence of the normalized pressure $P/T^{4}$, normalized energy density $\rho/T^{4}$, and trace anomaly $\left(\rho-3P\right)/T^{4}$ (dashed curves) calculated using our statistically corrected HRG model based on Eq. (\ref{eq:lnZ/}). Moreover, in Figure \ref{fig:one}, our model data are confronted to the corresponding lattice data taken from ref. \cite{bazavov2017qcd} (symbols with error bars) at vanishing baryon chemical potential $\mu_{b}=0 MeV$. Comparison is made for four different values of the correlation length $r$. Table \ref{tab1} lists $\chi^{2}$/dof statistic for the normalized pressure $P/T^{4}$, trace anomaly $\left(\rho-3P\right)/T^{4}$ and normalized energy density $\rho/T^{4}$ calculated in our statistically corrected hadron resonance gas (HRG) model and confronted to the corresponding lattice data from \cite{bazavov2017qcd} for four values of baryon chemical potential $\mu_{b}=0,170, 340$, and $425~$MeV.

At vanishing chemical potential, the best fit to the lattice data occurs for the case of zero correlation length which, in this case, corresponds to ideal HRG model. However, a slight exaggeration ($\chi^{2}/dof$ $\simeq$ 0.06 in the trace anomaly data) of our model's thermodynamics is observed in the vicinity of the critical phase transition temperature ($ T_{c}$ $\backsimeq$ 160 MeV). In the range $T\in [130, 200~$MeV$]$, the discrepancy between our model thermodynamics and the corresponding  lattice data is amplified. Generally, it is obvious that increasing the correlation length emphasizes the mismatch between our model and lattice data.

\begin{figure}[!htb]
\includegraphics[width=8.25cm]{HRGM_mubt=0Mev_r=0.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=0Mev_r=005.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=0Mev_r=01.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=0Mev_r=015.eps}
\caption{Normalized pressure $P/T^{4}$, normalized energy density $\rho/T^{4}$ and trace anomaly $\left(\rho-3P\right)/T^{4}$ (dashed curves) calculated using our statistically corrected HRG model and confronted to the corresponding lattice data taken from ref. \cite{bazavov2017qcd} (symbols with error bars), at $\mu_{b}=0~$MeV. Comparison is made for four different values of the correlation length $r$.}
\label{fig:one}
\end{figure}

For the case $\mu_{b}=170~$MeV and as it appears in Fig. \ref{fig:two}, the best fit generally occurs for $r=0$ and for $r=0.05~$fm. This good-fit temperature range extends from temperatures well below $T_{c}$ till $T \gtrsim T_{c} \backsimeq 160~$MeV. It is quite obvious here that the model fits the lattice data better compared to the corresponding vanishing chemical potential case(s). However, in the temperature range $T\in [170, 200~$MeV$]$, the mismatch of our model with the lattice data becomes more pronounced compared to the corresponding range of the vanishing chemical potential case(s).

\begin{figure}[!htb]
\includegraphics[width=8.25cm]{HRGM_mubt=170Mev_r=0.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=170Mev_r=005.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=170Mev_r=01.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=170Mev_r=015.eps}
\caption{The same as in Fig. \ref{fig:one} but at $\mu_{b}=170~$MeV.}
\label{fig:two}
\end{figure}

For the case of $\mu_{b}=340$ MeV, see Fig. \ref{fig:three}, the only interesting observation is that for $r=0.15~$fm, the model data well below and in the vicinity of $T_{c}$ and up to $T \backsimeq$ 170 MeV significantly approaches the corresponding lattice data. The data mismatch then diverges for higher temperatures.

\begin{figure}[!htb]
\includegraphics[width=8.25cm]{HRGM_mubt=340Mev_r=0.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=340Mev_r=005.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=340Mev_r=01.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=340Mev_r=015.eps}
\caption{The same as in Fig. \ref{fig:one} but at $\mu_{b}=340~$MeV.}
\label{fig:three}
\end{figure}

For the case of $\mu_{b}=425$ MeV, Fig. \ref{fig:four}, the data mismatch is generally too large to suggest any plausible correlation at any of values of $r$ and for all temperatures of interest.

\begin{figure}[!htb]
\includegraphics[width=8.25cm]{HRGM_mubt=425Mev_r=0.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=425Mev_r=005.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=425Mev_r=01.eps}
\includegraphics[width=8.25cm]{HRGM_mubt=425Mev_r=015.eps}
\caption{The same as in Fig. \ref{fig:one} but at $\mu_{b}=425~$MeV.}
\label{fig:four}
\end{figure}

\begin{table}[htb]
\begin{tabular}{|c|c|c|c|c|}
\hline
 $\mu_{b} $ MeV & $ r $ fm & $\chi^{2}$/dof for $P/T^{4}$ & $\chi^{2}$/dof for $\left(\rho-3P\right)/T^{4}$ & $\chi^{2}$/dof for  $\rho/T^{4}$ \\
\hline
  & 0 & $0.020 \pm 0.00065$ & $0.061 \pm 0.0025$ & $0.080 \pm 0.0034$ \\
\cline{2-5}
  & 0.05 & $0.013 \pm 0.00037$ & $0.097 \pm 0.0011$ &  $0.089 \pm 0.0012$ \\
\cline{2-5}
0  & 0.1 & $0.009 \pm 0.00051 $& $0.091 \pm 0.0021$ & $0.174 \pm 0.0033$ \\
\cline{2-5}
  &0.15  &$0.102 \pm 0.0010$  & $0.024 \pm 0.0018$ & $0.370 \pm 0.0052$ \\
\hline
  & 0 & $0.002 \pm 0.0002$ & $0.3698 \pm 0.0050$  & $0.036 \pm 0.0010$ \\
\cline{2-5}
  & 0.05 & $0.008 \pm 0.00041$ & $0.3190 \pm 0.0045$ & $0.062 \pm 0.002 $\\
\cline{2-5}
170  & 0.1 & $0.010 \pm 0.00062$ & $0.1452 \pm 0.0011$  & $0.075 \pm 0.0035$ \\
\cline{2-5}
  &0.15  & $0.028 \pm 0.0009 $& $0.069 \pm 0.0017 $& $0.145 \pm 0.0051 $\\
\hline
  & 0 & $0.020 \pm 0.0014 $& $0.853 \pm 0.0028 $& $0.241 \pm 0.0040$ \\
\cline{2-5}
  & 0.05 & $0.010 \pm 0.0016$ & $2.46 \pm 0.0062$ & $1.206 \pm 0.011$ \\
\cline{2-5}
340  & 0.1 & $0.005 \pm 0.0010$ & $1.31 \pm 0.0075$ &  $0.399 \pm 0.0093 $\\
\cline{2-5}
  & 0.15  &$0.015 \pm 0.0007$ & $0.335 \pm 0.0029$  & $0.101 \pm 0.005$ \\
\hline
  & 0 & $ 0.065 \pm 0.0012 $ & $ 2.12 \pm 0.15 $ & $ 4.28 \pm 0.0007 $  \\
\cline{2-5}
  & 0.05 & $0.351 \pm 0.0030$ & $6.0 \pm 0.0019$ &  $3.3 \pm 0.024$ \\
\cline{2-5}
425  & 0.1 & $0.035 \pm 0.0091$  & $2.77 \pm 0.0012 $&  $1.68 \pm 0.017$ \\
\cline{2-5}
  &0.15  & $0.058 \pm 0.0015$ & $1.44 \pm 0.0038$ & $0.376 \pm 0.0082$ \\
\hline
\end{tabular}
\caption{$\chi^{2}$/dof statistic for the normalized pressure $P/T^{4}$, trace anomaly $\left(\rho-3P\right)/T^{4}$, and normalized energy density $\rho/T^{4}$ calculated in our statistically corrected hadron resonance gas (HRG) model confronted to the corresponding lattice data taken from ref. \cite{bazavov2017qcd,bazavov2014equation}, at four values of baryon chemical potential $\mu_{b}=0,170,340$, and $425~$MeV.}
\label{tab1}
\end {table}

\section{Conclusions}
\label{sec:Cncls}

We confronted a novel statistical correction of the HRG model with recent lattice data  \cite{bazavov2017qcd,bazavov2014equation}. All our model calculations considered in this study doesn't seem to satisfactorily mimic the corresponding lattice data in the full temperature range under investigation, $T \in [130, 200~$MeV$]$. However, the best matching occurs locally in the vicinity of $T_{c} $ in the range $T \in [140, 170~$MeV$]$ for the case of $\mu_{b}=170~$MeV at zero and $0.05~fm$ correlation radii, $r$, respectively. Another remarkable matching between our model data with the corresponding lattice data occurs for the the case of $\mu_{b}=340~$MeV, at $0.1$ and $0.15~$fm correlation radii, respectively for temperatures $T \lesssim T_{c}$ and up to $T \backsimeq 170~$MeV. In the lower temperature range, $T \in [130, 160~$MeV$]$, most of the cases investigated in this research show reasonable match with the corresponding lattice data for different correlation lengths except for the case in which $\mu_{b}=425~$MeV where no good fitting is observed for any correlation length.

\bibliographystyle{aip}


\section{Introduction}

Depth completion has been widely studied in the field of robot navigation, computer vision, and autonomous driving. Its goal is to convert a sparse depth map from active depth sensors such as LiDAR or RGB-D cameras to a dense depth map.
In practice, robots with accurate sensing and dense depth maps have a lower chance of collision \cite{dinesh1} and are better able to compute safe navigation routes through challenging environments \cite{dinesh2,dinesh3}. Currently, most recent methods for depth completion rely on supervised learning \cite{CSPN++,2020Eldesokey,GuideNet}, requiring large amounts of high-quality and dense depth maps as ground truth.
Unfortunately, generating dense depth maps is expensive and challenging due to the sparse and noisy depth values from other active sensors. Even existing expensive 3D LiDAR sensors only have a limited number of scan lines and provide sparse depth measurements.

Self-supervised methods with monocular videos that are able to train depth completion neural networks without the need for dense depth maps. These methods \cite{2017Godard,SfMLearner,2018DDVO,Monodepth2} are based on an unsupervised learning framework that use depth and camera ego-motion from monocular videos and use photometric loss as supervision. This loss measures the difference between a reference image and a synthesized image obtained by the depth-guided reprojection of other views into the reference view. However, the previous self-supervised training process for depth completion \cite{2019ICRAMa, 2019CVPRYang} may not work well on public indoor datasets such as NYUv2 \cite{NYUv2}. This is due to large areas of textureless regions and scenes captured by complex ego-motions of a handheld camera \cite{MovingIndoor,weakrect}.
Large indoor environments are more challenging not only because of  the above issues but also because of large regions of glossy and transparent surface,
non-Lambertian surfaces, longer and diverse ranges, and moving people. To overcome these issues, our proposed training process requires effective ways to compute complex camera ego-motions attached to a robot, mask out regions where photometric loss is unreliable, and output depth maps robust to diverse depth range. In this paper, we address the problem of developing self-supervised depth completion methods that work robustly on general indoor environments.

 \noindent\textbf{Main Results:} We improve the training framework by leveraging learning-based local features to obtain accurate relative poses, making robust matching with auto-masking and minimum reprojection loss \cite{Monodepth2}, and predicting inverse depth using a sigmoid function. However, we observe that merely masking out all regions still leads to depth artifacts including blurred boundaries and irregular patterns copied from sparse depth input. Therefore, we also propose a novel architecture with two unique components: sparsity-aware convolutions that extract features properly from sparse input and pixel-adaptive convolutions which propagate this information to the whole depth prediction with the guidance of image features. We highlight the effectiveness of our method by evaluating our approach on challenging datasets such as the Metro Station or the Department Store (Fig. \ref{intro_figure}). We also demonstrate the applicability of our method on datasets for which prior self-supervised depth estimation and completion fail to provide reliable depth maps. \textbf{Our novel contributions} can be summarized as follows:

\begin{figure}[t]

    \centering

    \includegraphics[width=0.87\linewidth]{ICRA_figures/introv2.png}
    \vspace{-2mm}
    \caption{Depth completion results in challenging scenes: Metro Station (top) and the Department Store (bottom) from NAVERLABS indoor dataset. Our method takes an RGB image and a sparse depth map (projection of LiDAR pointcloud) as input, and predicts a dense depth map.

    (low \protect\includegraphics[width=0.5cm,height=0.2cm]{ICRA_figures/colorbar.jpeg} high; grey means empty depth values.)}
    \label{intro_figure}
    \vspace{-4mm}
\end{figure}

\begin{itemize}
    \item We develop a robust and stable training framework for self-supervised depth completion in challenging indoor environments.
    \item We propose a novel architecture that can extract sparse depth features guided by high-level image features.
    \item Our method outperforms prior self- and fully- supervised methods on challenging indoor datasets like NAVERLABS. We also achieve the state-of-the-art performance among self-supervised methods on public benchmark datasets, including KITTI and the NYUv2 dataset based on various evaluation metrics.
\end{itemize}

\section{Related Work}
Recently, depth completion methods using a deep neural network have shown promising results. Uhrig et al. \cite{SparsityCNN} propose a sparsity-invariance convolution layer to handle sparse depth or features by using mask normalizations. Eldesokey et al. \cite{2018Eldesokey,2019Eldesokey} propose the normalized convolution to produce both a confidence map and dense depth output. In \cite{2020Eldesokey}, they further extend their work into uncertainty estimation with depth completion methods. Teixeira et al. \cite{Aerial} use this confidence propagation from \cite{2018Eldesokey} to compute dense depth from LiDAR.
Cheng et al. \cite{2018ECCVcspn,CSPN++} utilize the neural network to learn the affinity among neighboring pixels in order to perform the propagation process. In addition to using RGB images, some methods take advantage of information from other modalities, e.g., surface normals \cite{2018Funkhouser,deeplidar,depthnormal} or segmentation \cite{2018segmentation}. Yang et al. \cite{2019CVPRYang} predict the posterior over depth maps with a conditional prior network. Some works \cite{2019ICCVWIndoor,2020BMVCDecoder} only focus on an indoor depth completion task. However, all aforementioned works on indoor scenarios are limited to a small-scale indoor environment
and handle non-LiDAR indoor data (e.g., on the NYU-Depth-v2 \cite{NYUv2} and Matterport3D \cite{2018Funkhouser} datasets). In this paper, we introduce challenging indoor scenarios such as the department store or metro station containing LiDAR sensor data.

The lack of dense depth maps is a significant limitation of supervised depth completion. As an alternative, several current works study self-supervised depth completion methods.
Ma et al. \cite{2019ICRAMa} present self-supervised methods based on a photometric loss with pose estimation using the PnP method. Yang et al. \cite{2019CVPRYang} present a conditional prior network (CPN) that predicts a probability of depths at each pixel. Zhang et al. \cite{Dfinenet} propose an end-to-end framework to jointly estimate pose and depth without the PnP method. Yoon et al. \cite{2020IROSYoon} combine monocular depth estimation and Gaussian process-based depth regression for depth completion.

\section{Proposed Method}
\subsection{Training Framework}
In Fig. \ref{figure1}, our proposed depth completion model takes a target view image \(I_{t}\) and a sparse depth \(S_{t}\) as input, and outputs a dense inverse-depth \(\hat{d}_{t}\). We invert \(\hat{d}_{t}\) to compute a dense depth map \(\hat{D}_{t}\). In our method, the source view images contains its two adjacent temporal frames, i.e., \(I_{s} \in \{I_{t-1}, I_{t+1}\}\), although including a larger temporal context is possible. Also, we estimate the relative pose \(T_{t\rightarrow s}\) between the target image \(I_{t}\) and the set of source images \(I_{s}\) by using a Perspective-n-Point (PnP) \cite{PnP} with Random Sample Consensus (RANSAC) \cite{RANSAC}.
\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{ICRA_figures/figure1v3.pdf}
\vspace{-6mm}
\caption{Overview of our proposed framework.
The black block with a dash-dotted line shows the depth completion network (details in Fig. \ref{figure2}).
}
\vspace{-4mm}
\label{figure1}
\end{figure}

\subsubsection{Relative Pose Estimation}
Although most prior works \cite{Monodepth2,Dfinenet} adopt learning-based pose estimation, we observe that learning-based pose estimation reveals degradation in performance when applied to challenging environments such as indoor scenes \cite{MovingIndoor,towardsbetter,weakrect}. Additionally, according to recent findings \cite{sattlerlearningpose}, learning-based pose estimation methods suffer from scale ambiguity problems and make it difficult to provide absolute-scale depth supervision. Unlike self-supervised depth estimation, pose estimation via the PnP method is aligned with the depth completion task because our model needs to predict absolute scale depth.

Our method first detects learning-based features and computes their corresponding descriptors via learning-based algorithms. We select R2D2 \cite{r2d2}, which shows better performance than traditional hand-crafted features on the visual localization benchmarks.
Given the predicted feature in source and target images, we can compute high-quality correspondences by using a FLANN-based search algorithm \cite{FLANN} and then output the relative camera pose by solving PnP with corresponding sparse depths.
Compared to hand-crafted features used by \cite{2019ICRAMa}, dense correspondences obtained by R2D2 enable us to fully utilize sparse depth measurements with corresponding feature points for PnP pose estimation.

\subsubsection{Photometric Loss}
Given the camera intrinsic matrix \(K\) and the relative pose \(T_{t\rightarrow s}\), we compute the warped pixel coordinates and synthesize target image \(I^\prime_{t}\) from source images \(I_{s}\) via a bilinear sampling function. Following the self-supervised depth estimation \cite{2017Godard,SfMLearner}, we use the combination of the L1 and SSIM \cite{SSIM} as the photometric loss \(L_{ph}\):

\begin{equation}
\begin{split}
    L_{ph}(I_{t}, I^\prime_{t}) &= \alpha \frac{1-\text{SSIM}(I_{t}, I^\prime_{t})}{2}
    + (1-\alpha)\parallel I_{t}-I^\prime_{t} \parallel, \\
    I^\prime_{t}(x) &= I_{s}\langle \pi( KT_{t\rightarrow s}\hat{D}_{t}(x)K^{-1}x_{h}) \rangle\:.
\end{split}
\label{photometric_loss}
\end{equation}

where \(x_{h} = [x^{T}, 1]^{T}\) are the homogeneous coordinates from the target image, \(\pi\) gives pixel coordinates,  and \(\langle \cdot \rangle\) denotes the bilinear sampling function. We adopt the minimum reprojection loss and auto-masking from \cite{Monodepth2}. Our methods often face the occlusion issue, meaning that some pixels from frame \(I_{t}\) do not appear in either \(I_{t-1}\) or \(I_{t+1}\) frame due to moving people. To mitigate this occlusion issue, our per-pixel minimum photometric loss is defined as follows,
\begin{equation}
\begin{split}
    L_{photo} = \min_{I_{s}} L_{ph}(I_{t}, I^\prime_{t}).
\end{split}
\label{minimum_loss}
\end{equation}
Additionally, auto-masking encourages our method to mask out static pixels for which the appearance does not change between the current frame and temporally adjacent frames. This auto-masking results in removing pixels belonging to low-texture scenes (e.g. ceiling or fluorescent light) that have undesirable effects on the minimum photometric loss. Furthermore, these masks often remove non-Lambertian surfaces (e.g. mirrors or digital bilborad screens), which are common in the department store or the metro station.

\subsubsection{Sparse Depth Loss}
Our proposed model can learn to encode absolute scale from sparse depth loss by aligning the depth prediction \(\hat{D}_{t}\) with sparse depth \(S_{t}\). We use the L1-norm to penalize the difference between \(\hat{D}_{t}\) and \(S_{t}\) over \(\Omega\), where sparse depth is available,
\begin{equation}
\begin{split}
    L_{depth} = \sum_{x \in \Omega}|\hat{D}_{t}(x) - S_{t}(x)| \:.
\end{split}
\label{depth_loss}
\end{equation}

\subsubsection{Total Loss Function}
Following \cite{2017Godard,2019ICRAMa,Monodepth2}, we also incorporate an edge-aware smoothness loss function \(L_{smooth}\) over the predicted depth map to regularize the depth in texture-less low-image gradient regions.
\begin{equation}
    L_{smooth} = {\mid\partial_{x}d^{*}\mid}e^{-{\mid\partial_{x}I_{t}\mid}} + {\mid\partial_{y}d^{*}\mid}e^{-{\mid\partial_{y}I_{t}\mid}},
\label{smooth_loss}
\end{equation}
where \(\partial_{x}, \partial_{y}\) denotes the gradients along the either the x or y direction. \(d^{*} = \hat{d}_{t} / {\overline{\hat{d}_{t}}}\) is the normalized inverse depth \cite{2018DDVO}.
Our overall loss function is a weighted sum over the aforementioned losses,
\begin{equation}
    L_{tot} = L_{photo} + \lambda_{d}L_{depth} + \lambda_{s}L_{smooth}
\label{total_loss}
\end{equation}
where \(\lambda_{d}\) and \(\lambda_{s}\) denote weighting terms selected through a grid search.

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{ICRA_figures/Mask.jpg}
\vspace{-3mm}
\caption{Examples of auto-masks computed after one epoch. Auto masks (black pixels) usually contain low-texture scenes or non-Lambertain surfaces where photometric loss is unreliable.}
\label{figure:auto-masks}
\vspace{-4mm}
\end{figure}

\subsection{Network Architecture}
In our training framework, the objective function defines a trade-off between a sparse depth loss, which enables our network to predict absolute scale depth on sparse regions, and a photometric loss, which encourages the network to predict dense depth where sparse information does not exist. In challenging indoor environments, self-supervised methods with photometric loss have shown that finding the optimal depth value is very difficult because the loss only comes from the appearance difference between target and synthesized images; measuring these differences is difficult due to the large non-texture regions, repeating structures, and occlusions.

The sparse depth loss might easily overwhelm the whole self-supervised training, and the depth predictions often result in overfitting on sparse depth input. Both feature extraction from the sparse depth and feature fusion schemes of multi-modal data are critical to avoid such negative effects. Our network adopts a late fusion scheme that uses each encoder for different sources because image and sparse depth features are heterogeneous. Throughout several experiments in Sec. \ref{subsubsection:naverlabs}, a late fusion scheme is shown to be particularly suitable for self-supervised depth completion. We focus on designing an effective architecture that leverages sparse depth features to fill the empty regions rather than getting stuck producing sparse depth prediction.

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{ICRA_figures/figure2.pdf}
\caption{ Our model consists of three main components: (1) an RGB encoder that extracts RGB features, reducing the spatial resolution of the input images to consider the global context; (2) a depth encoder that extracts sparse depth features from sparse input via sparse convolutional blocks; (3) a decoder that takes both RGB and depth features, and transforms them into dense depth maps. Data flow contains convolution operations.
}
\label{figure2}
\vspace{-6mm}
\end{figure}

\subsubsection{Sparse Convolution}
In order to extract features from sparse depth, recent depth completion methods fall into two groups. One group \cite{SparsityCNN,2018Eldesokey} designs shallow networks with sparse convolution layers due to the poor performance of standard convolution layers on sparse input. In contrast, the other group \cite{2018ECCVcspn,2019ICRAMa,CSPN++} demonstrates that employing standard convolution layers yields better performance with their training schemes; these methods also perform well on different levels of sparsity. Comparing both approaches, we conclude that deep stacks of standard convolution layers
are key to handling sparse depth input. Therefore, we choose to combine these two approaches by stacking sparse convolution layers \cite{SparsityCNN} deeper. Sparse convolution blocks take a sparse depth feature and a binary observation mask as inputs, which both have the same spatial size. Since an observation mask indicates the reliability of each location of the feature maps, the convolution operations only attend on reliable input features which are involved with sparse depth input. As a result, the depth encoder is able to propagate depth information from reliable sparse depth input to its surrounding depth features.
In contrast to sparsity invariant convolutions \cite{SparsityCNN}, we do not apply mask normalization, which causes a degradation problem due to the small values extracted from the depth encoder when stacking deep layers.

\subsubsection{Pixel-Adaptive Convolution (PAC)}
Recent works utilize operations like concatenation or element-wise addition to fuse the information from both the RGB and the depth encoder. However, this simple method is not the optimal way to fuse the information from two modalities because extracted depth features are not guided by RGB images in a late-fusion manner. To achieve this, we require spatially-variant convolution weights that different convolution kernels depending on image features are applied to different spatial positions of the depth features. Thus, we adopt pixel-adaptive convolutions (PACs) introduced by Su et al. \cite{2019CVPRpacnet}. Assume that the PACs take a sparse depth feature map \(\textbf{v}\) and an image feature \(\textbf{f}\) as input, and produce the output sparse depth feature \(\textbf{v}^\prime\). We denote the convolution weight as \(\textbf{W}\) and the bias term as \(\textbf{b}\). PACs augment the spatially invariant convolution weight \(\textbf{W}\) with a local-adaptive kernel function \(\textbf{K}\). Following the notation in \cite{2019CVPRpacnet}, the PAC convolution is defined as
\begin{equation}
\begin{split}
    \textbf{v}^\prime_{i} = \sum_{j\in \Omega(i)} K(\textbf{f}_i, \textbf{f}_j) \textbf{W}[\textbf{p}_i - \textbf{p}_j]\textbf{v}_j + \textbf{b}
\end{split}
\label{photometric_loss}
\end{equation}
where \(\textbf{f}\) are image features from the RGB encoder and act as guidance for the PACs in the depth encoder. High-level image features are preferable because they contain image context information. The kernel function \(\textbf{K}\) has a standard Gaussian form: \(K(\textbf{f}_i, \textbf{f}_j) = exp(-\frac{1}{2}(\textbf{f}_i - \textbf{f}_j)^T(\textbf{f}_i - \textbf{f}_j))\). This kernel computes the correlation between image features and guides the standard convolutional weights to generate depth features depending on the content learned from the RGB encoder. This spatially-variant kernel is helpful because, for instance, depth features on a person should not be applied to generate depth maps on the wall. Thus, PACs enable our proposed architecture to generate intermediate depth features consistent with image content.

\subsubsection{Decoder}
The decoder is similar to Ma et al \cite{2019ICRAMa} with one major difference. Unlike their method, our proposed network predicts inverse depth instead of depth. Thus, we replace the ReLU nonlinearity function with sigmoid activation. We observe that applying a sigmoid function is critical for training stability and is more robust to diverse depth ranges on challenging indoor scenes. Without inverse depth prediction, we fail to train the network proposed by \cite{2019ICRAMa} on our challenging indoor datasets.

\section{Implementation and Performance}
\subsection{Datasets} \label{subsection:Datasets}
To demonstrate that our proposed method performs robustly in the challenging environment, we experiment with NAVERLABS indoor dataset. Additionally, we evaluate our method on two benchmarks: the NYUv2 dataset and KITTI. For indoor scenes, we report five evaluation metrics: root mean squared error (RMSE), mean absolute relative error (ABS Rel), and $\delta_{i}$ ,which means the percentage of predicted pixels for which the relative error is less than a threshold \(i\). The RMSE is represented in \(mm\) scale. For the KITTI benchmark, we use four official metrics for evaluation.

\subsubsection{NAVERLABS indoor dataset} \label{subsubsection:naverlabs}

NAVERLABS dataset
\cite{spoxelnet}
comprises data collected from three different places: Department Store B1 (Dept. B1), Department Store 1F (Dept. 1F), and Metro Station (MS). All data were collected using a mapping robot equipped with two 16-channel LiDAR sensors and six RGB cameras with 2592\(\times\)2048 resolution. LiDAR SLAM is used to estimate sensor poses for our dataset, and this pose information enables us to generate depth maps. We captured 17K images for Dept. 1F, 32K images for Dept. B1, and 61K images for MS.
The MS dataset was captured at a highly crowded metro station.

We accumulate consecutive LiDAR sweeps and project the point clouds onto each of six cameras in order to generate depth maps. These depth maps can be used as sparse input and ground truth for evaluation. In order to create sparse input, we accumulate 0.3 seconds of LiDAR scans, which correspond to one LiDAR scan since the LiDAR is sampled at 10 Hz (see the second column in Fig. \ref{figure4}). In the Dept. 1F and Dept. B1 datasets, we collect 1 second of LiDAR scans to generate ground truth depth for training supervised depth completion algorithms (see the third column in Fig. \ref{figure4}). However, the depth maps generated from the MS dataset are too noisy to serve as ground truth due to noise caused by moving objects (see SS-(c) in Fig. \ref{figure4}). For evaluation, we select the high-quality depth maps generated from projecting a LiDAR 3D model onto the camera frames in the Dept. 1F and Dept. B1 datasets (see 1F-(d) and B1-(d) in Fig. \ref{figure4}). In contrast, we use the SfM method \cite{SfM_revisited} to reconstruct a sparse 3D model of the MS dataset and MVS system \cite{MVS} to produce dense depth maps instead of noisy depth maps from LiDAR point clouds.
Given the pose from LiDAR SLAM, we can obtain a 3D model aligned with the LiDAR data. SS-(d) in Fig. \ref{figure4} generated by the SfM is the example of depth maps for evaluation.

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{ICRA_figures/dataset_v4.pdf}
\vspace{-3mm}
\caption{Samples from the NAVERLABS indoor dataset. We include more samples and demonstrate many challenging features from our datasets in supplementary video.
}
\vspace{-6mm}
\label{figure4}
\end{figure}

\subsubsection{NYU-Depth-v2 dataset}
The NYUv2 dataset \cite{NYUv2} comprises RGB images and depth maps collected from Microsoft Kinect in 464 indoor scenes. The raw training dataset contains 268K images. We downsampled 10 times to remove redundant frames, resulting in 47k frames sampled for the training set. Our method is tested on the 654 official labeled test set for evaluation.

\subsubsection{KITTI depth completion dataset}
The KITTI depth completion dataset \cite{SparsityCNN} provides RGB images, sparse lidar points, and dense ground truth depth. Following the official split, it contains 86k frames for the training set, 1k frames for the validation set, and 1k frames for the evaluation.

\subsection{Implementation Details}

We adopt ResNet34 \cite{ResNet} as the RGB encoder, and it was pretrained on ImageNet \cite{ImageNet}.
The proposed deep networks were implemented with PyTorch \cite{PyTorch} and trained on 4 Tesla V100 GPUs.
We trained our network from scratch for 30 epochs, with a batch size of 16 using Adam \cite{Adam}, where \(\beta_{1}=0.9\), \(\beta_{2}=0.999\); we used an input resolution of 1024\(\times\)760 for NAVERLABS, 640\(\times\)480 for NYUv2, and 1024\(\times\)320 for KITTI. We used an initial learning rate of \(10^{-4}\) for the first 10 epochs and halved it every 10 epochs. In loss functions, we set the \(\alpha\) to 0.85, \(\lambda_{d}\) to 0.001 and \(\lambda_{s}\) to 0.1.

\subsection{Experiments on NAVERLABS indoor datasets}

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{ICRA_figures/figure3_v2.pdf}
\vspace{-3.2mm}
\caption{
SD represents a sparse depth. We design 4 architectures: (a) 1e-1d: an encoder and a decoder; (b) 2e-1d: an RGB encoder, a depth encoder, and a decoder; (c) 2es-1d: an RGB encoder, a depth encoder with sparse convolutions, and a decoder; (d) 2es-1dp: the same architecture as Fig. \ref{figure2}.
}
\vspace{-2mm}
\label{fig:Ablation}
\end{figure}

\begin{table}[t]
\begin{center}
\caption{\label{table:NAVERLABS} Depth completion results on NAVERLABS indoor dataset.}
\resizebox{0.48\textwidth}{!}{
\begin{tabular}{lcccccc}
\toprule
\multicolumn{7}{c}{(a)  Dept. 1F Dataset}\\
\midrule
Method & T & RMSE \(\downarrow\) & ABS Rel \(\downarrow\) & $\delta_{1.25}$ \(\uparrow\) & $\delta_{1.25^2}$ \(\uparrow\) & $\delta_{1.25^3}$ \(\uparrow\) \\
\midrule
UNET \cite{2019ICRAMa}& S & 5043.5 & 0.571 & 0.273 & 0.463 & 0.553 \\
nUNET \cite{2019Eldesokey}& S & 4993.8 & 0.171 & 0.898 & 0.933 & 0.951 \\
CSPN \cite{2018ECCVcspn}& S  & 2857.1 & 0.266 & 0.808 & 0.886 & 0.924 \\
MD2 \cite{Monodepth2}& SS  & 3862.4 & 0.493 & 0.488 & 0.716 & 0.818 \\
\midrule
1e-1d & SS  & 2899.6 & 0.256 & 0.784 & 0.869 & 0.911 \\
2e-1d & SS  & 2801.4 & 0.265 & 0.783 & 0.868 & 0.909 \\
2es-1d & SS  & 2700.3 & 0.242 & 0.790 & 0.875 & 0.915 \\
Ours & SS & \textbf{2692.9} & \textbf{0.234} & \textbf{0.815} & \textbf{0.887} & \textbf{0.926} \\
\bottomrule
\\
\toprule
\multicolumn{7}{c}{(b) Dept. B1 Dataset}\\
\midrule
Method & T & RMSE \(\downarrow\) & ABS Rel \(\downarrow\) & $\delta_{1.25}$ \(\uparrow\) & $\delta_{1.25^2}$ \(\uparrow\) & $\delta_{1.25^3}$ \(\uparrow\) \\
\midrule
UNET \cite{2019ICRAMa}& S & 3048.5 & 0.346 & 0.599 & 0.711 & 0.764 \\
nUNET \cite{2019Eldesokey}& S & 2963.1 & 0.324 & 0.635 & 0.732 & 0.779  \\
CSPN \cite{2018ECCVcspn}& S  & 2162.2 & 0.199 & 0.774 & 0.890 & 0.938  \\
MD2 \cite{Monodepth2}& SS  & 3592.7 & 0.338 & 0.47 & 0.729 & 0.862  \\
\midrule
1e-1d & SS  & 2210.5 & 0.186 & 0.823 & 0.914 & 0.947 \\
2e-1d & SS  & 2157.5 & 0.170 & 0.832 & 0.916 & 0.949 \\
2es-1d & SS  & 2116.8 & 0.183 & 0.829 & 0.916 & 0.949\\
Ours & SS & \textbf{2024.1} & \textbf{0.178} & \textbf{0.841} & \textbf{0.922} & \textbf{0.953} \\
\bottomrule
\\
\toprule
\multicolumn{7}{c}{(c) Metro Station Dataset}\\
\midrule
Method & T & RMSE \(\downarrow\) & ABS Rel \(\downarrow\) & $\delta_{1.25}$ \(\uparrow\) & $\delta_{1.25^2}$ \(\uparrow\) & $\delta_{1.25^3}$ \(\uparrow\) \\
\midrule
UNET \cite{2019ICRAMa}& S & 4029.7 & 0.749 & 0.201 & 0.219 & 0.23  \\
nUNET \cite{2019Eldesokey}& S & 2705.1 & 0.472 & 0.444 & 0.528 & 0.582 \\
MD2 \cite{Monodepth2}& SS  & 1321.2 & 0.166 & 0.8 & 0.916 & 0.939 \\
Ours & SS & \textbf{1248.4} & \textbf{0.138} & \textbf{0.866} & \textbf{0.923} & \textbf{0.942} \\
\bottomrule
\end{tabular}
}
\end{center}
``S'' represents \textit{supervised} training methods and ``SS'' denotes \textit{self-supervised} training methods. Different from other depth completion methods, MD2 is based on self-supervised depth estimation with monocular video sequences. Compared to prior self- and fully- supervised methods, our proposed method demonstrates improvement in RMSE metric by 5.8\% in Dept. 1F, by 6.4\% in Dept. B1 and by 5.5\% in Metro Station.
\vspace{-6mm}
\end{table}

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{ICRA_figures/Analysisv2.jpg}
\vspace{-2mm}
\caption{Examples of depth artifacts generated by 1e-1d model. The top (NYUv2) and bottom (Dept. 1F) figures in third column shows irregular sparse patterns copied from sparse inputs and blurred depth boundaries}
\label{figure:Analysis}
\vspace{-6mm}
\end{figure}

\begin{figure*}[t]
    \centering

    \includegraphics[width=0.85\linewidth]{ICRA_figures/comparisonv2.pdf}
    \vspace{-2mm}
    \caption{Qualitative comparison with ``UNET''  \cite{2019ICRAMa}, ``nUNET''  \cite{2019Eldesokey}, and ``CSPN'' \cite{2018ECCVcspn} on the NAVERLABS indoor dataset.

    As mentioned in Sec. \ref{subsubsection:naverlabs}, the groundtruth of the Metro Station dataset is the results of MVS.

    (low \protect\includegraphics[width=1cm,height=0.2cm]{ICRA_figures/colorbar.jpeg} high; grey means empty depth values.)
    }
    \label{fig:comparison}
    \vspace{-4mm}
\end{figure*}

We compare our proposed method with 4 baseline approaches with public codes: 1) UNET is a U-shaped ResNet \cite{2019ICRAMa}; 2) nUNET is an early-fusion encoder-decoder network combined with multiple normalization convolutions \cite{2019Eldesokey}; 3) CSPN learns the affinity to perform sparse depth propagation \cite{2018ECCVcspn}; 4) Monodepth2 (MD2) is based on ResNet34 \cite{Monodepth2}. The same relative pose generated from our method is applied to training MD2 for a fair comparison. MD2 benefits from the use of ground truth median scaling at inference, while other depth completion methods utilize unmodified network prediction. Except for 4), the aforementioned methods are based on supervised training with ground truth depth. As mentioned in Section \ref{subsection:Datasets}, it is difficult to create ground truth depth for the Metro Station dataset (MS). In the case of the MS dataset, we only use the sparse depth maps as weak supervisions for supervised depth completion. Quantitative results are shown in Table \ref{table:NAVERLABS}, and the associated qualitative comparisons are shown in Fig. \ref{fig:comparison}. We train the model from Ma et al. \cite{2019ICRAMa} with its network architecture and the relative poses from our method, which are more accurate. Note that Ma et al. failed to produce meaningful results. We observe that the network just copies the sparse inputs, meaning that the sparse depth loss only serves as supervision for training while the photometric loss is unreliable.

In contrast, our method successfully performs depth completion in the NAVLERLABS dataset. This indicates that our training framework is more robust to variant supervision compared to other self-supervised depth completion methods. The performance of our method is superior to existing supervised completion methods. From columns 4 through 6 of Fig. \ref{fig:comparison}, we can observe that supervised depth completion methods rely heavily on the density of ground truth depth. We notice that these algorithms fail to fill the holes where LiDAR information does not exist (e.g., the upper part of the image) due to the sensor configuration. In contrast, thanks to the effectiveness of photometric loss, the proposed method succeed in filling empty regions such as the ceiling, which LiDAR sensors are not able to measure, with accurate depth values. In addition, our method even predicts the correct depth map on shiny or transparent surface which cannot be scanned accurately by sensors (see First row in Fig. \ref{fig:comparison}). The experimental results of CSPN bolster this observation because our method achieves performance comparable to the supervised setting in Table \ref{table:NAVERLABS}-(a), (b), while its training with weak supervision collapsed on the MS dataset. If we collect the high-quality dense depth maps (see the last column in Fig. \ref{fig:comparison}), the supervised depth completion methods can be superior to our method. However, collecting dense depth maps requires a high cost and can be impractical in certain environments such as the MS dataset.

\begin{table}[t]
\begin{center}
\caption{\label{table:NYUv2} depth completion results on the NYUv2 dataset}
\resizebox{0.48\textwidth}{!}{
\begin{tabular}{lcccccc}
\toprule
Method & T & RMSE \(\downarrow\) & ABS Rel \(\downarrow\) & $\delta_{1.25}$ \(\uparrow\) & $\delta_{1.25^2}$ \(\uparrow\) & $\delta_{1.25^3}$ \(\uparrow\) \\
\midrule
\multicolumn{7}{l}{(a) 500 Samples}\\
Ma et al. \cite{2018ICRAMa} & S & 0.204 & 0.043 & 97.8 & 99.6 & 99.9 \\
CSPN \cite{2018ECCVcspn} & S & 0.117 & 0.016 & 99.2 & 99.9 & 100 \\
Ma et al. \cite{2019ICRAMa} & SS & 0.271 & 0.068 & - & - & - \\
Ours & SS & 0.178 & 0.033  & 98.1 & 99.7 & 100 \\
\midrule
\multicolumn{7}{l}{(b) 200 Samples}\\
Ma et al. \cite{2018ICRAMa} & S & 0.230 & 0.044 & 97.1 & 99.4 & 99.8 \\
Yang et al. \cite{2019CVPRYang} & SS & 0.569 & 0.171  & - & - & - \\
Yoon et al. \cite{2020IROSYoon} & SS & 0.309 & -  & - & - & - \\
Ours  & SS & 0.240 & 0.048 & 96.6 & 99.4 & 99.9 \\
\bottomrule
\end{tabular}
\vspace{-3mm}
}

\caption{\label{table:KITTI} depth completion results on the KITTI validation set}
\resizebox{0.48\textwidth}{!}{
\begin{tabular}{lccccc}
\toprule
Method & T & RMSE [mm] & MAE [mm] & iRMSE [1/km] & iMAE [1/km] \\
\midrule
Uhrig et al. \cite{SparsityCNN} & S & 1601.33 & 481.27 & 4.94 & 1.78  \\
Ma et al. \cite{2019ICRAMa} & S & 814.73 & 249.95 & 2.80 & 1.21  \\
\midrule
Yoon et al. \cite{2020IROSYoon} & SS & 1593.37 & 547.00 & 27.98 & 2.36 \\
Ma et al. \cite{2019ICRAMa} & SS & 1343.33 & 358.66 & 4.28 & 1.64  \\
Yang et al. \cite{2019CVPRYang} & SS & 1310.03 & 347.17 & - & - \\
Ours & SS & 1212.89 & 346.12 & 3.54 & 1.29 \\
\bottomrule
\end{tabular}
}

\end{center}
``S'' represents \textit{supervised} training methods and ``SS'' denotes \textit{self-supervised} training methods.
Compared to prior self-supervised method, our method reduces the RMSE by 34.3\% in NYUv2 and by 7.4\% in KITTI.
\vspace{-6mm}
\end{table}

To investigate the impact of our network architecture, we conduct ablation studies for different network configurations on the NAVERLABS dataset. Figure \ref{fig:Ablation} illustrates the network variants, and their quantitative results are shown in Table \ref{table:NAVERLABS}. Comparing the results of 1e-1d with those of 2e-1d, we observe that fusing the RGB and depth feature at the earliest stage worsens the results. With the early fusion, we often observe depth with blurrier depth discontinuities, as seen in Fig. \ref{figure:Analysis}. With PAC, depth maps can eliminate either irregular lattice patterns from LiDAR or sparse points from sampled depth maps. We conclude that the late fusion is more effective for self-supervised learning because image and depth features are heterogeneous data, and early fusion is more likely to produce noisy initial depth values.

\subsection{Experiments on Public Benchmark Datasets}

To verify the performance of our method on public datasets, we train and evaluate our method on the NYUv2 dataset. Following previous works, we uniformly sample 200 and 500 sparse depth points separately for sparse input. Except for supervised depth completion methods, our proposed method achieves the best performance among self-supervised depth completion methods in Table \ref{table:NYUv2}.

In order to prove the generalization capability of our method, we train and evaluate our method on an outdoor dataset, the KITTI depth completion dataset. The comparison results are shown in Table \ref{table:KITTI}. We observe that our method outperforms other self-supervised depth completion methods on every metric. Although Yang \cite{2019CVPRYang} pretrains their network using a virtual dataset such as the Virtual KITTI dataset \cite{VKITTI}, which is similar to the benchmark dataset, our method performs better than their method. We include qualitative results in the supplementary video.

\section{Conclusions, Limitations, and Future Work}
In this paper, we investigated the self-supervised depth completion framework for a challenging indoor environment. We provide two main research contributions. First, we introduced a training framework that leads to stable self-supervised training. Second, we proposed a novel architecture to combine image features and sparse depth features effectively. In indoor environments, our method outperforms other self-supervised methods and is competitive with supervised methods. Furthermore, it generalizes better to all public benchmark datasets. As part of future work, we would like to apply this method to complement depth sensors for mobile robots in indoor environments and and evaluate its benefits for robot navigation. A limitation of the current method is the high computational cost for embedded platforms. For robotic tasks, we would like to explore network compression algorithms to reduce computational complexity and latency.

\addtolength{\textheight}{0cm}

\section*{ACKNOWLEDGMENT}
This work was supported by the Institute of Information \& Communications Technology Planning \& Evaluation(IITP) grant funded by the Korea government(MSIT) (No. 2019-0-01309, Development of AI Technology for Guidance of a Mobile Robot to its Goal with Uncertain Maps in Indoor/Outdoor Environments). This work was supported in part by ARO Grants W911NF1910069, W911NF2110026  and W911NF1910315 and NSF grant 2031901.

\newpage


\section{Introduction}

\acresetall \bstctlcite{IEEEexample:BSTcontrol}Cell-free massive \ac{MIMO}
can be considered as the most recent development of distributed massive
\ac{MIMO} \cite{Ngo2017cfmm}, whereby a large number of \acp{AP}
serve a group of users in a large area of service. Cell-free massive
\ac{MIMO} has been receiving increasing attention as a potential
backbone technology for beyond 5G networks due to the massive performance
gains compared to the colocated massive \ac{MIMO} \cite{Ngo2017cfmm}.
While the research on cell-free massive \ac{MIMO} has been focusing on
the downlink, that is not the case for the uplink. For uplink cell-free
massive \ac{MIMO}, linear receivers are commonly employed at \acp{AP}
to combine users' signals. As a result, an important problem is to
jointly design linear combining and users' powers to maximize a performance
measure. In the following, we attempt to provide a comprehensive survey,
but by no means inclusive, on the previous studies on this line of
research.

In \cite{Ngo2018rician}, the \ac{MRC} technique was employed at the
\acp{AP} based on local channel estimates to detect signals from single-antenna
users. Closed-form expression for achievable \ac{SE} is derived,
which is then used to formulate a max-min fairness problem. To solve
this problem, the work of \cite{Ngo2018rician} applies an \ac{AO} method by solving two subproblems: the linear combining
coefficient design which is in fact a generalized eigenvalue problem,
and the power control problem which is quasi-linear, and thus, can be solved by a bisection method with linear programming. A similar problem
is considered in \cite{Mai2019uplinkMU}, except the zero-forcing scheme
was used at \acp{AP}, and a similar solution to \cite{Ngo2018rician} is
proposed. In \cite{chamalee2019alternating}, the \ac{L-MMSE} combining
is adopted at each \ac{AP}. A similar \ac{AO} method was presented to
solve the max-min fairness problem. In particular, the power control
subproblem (when the \ac{L-MMSE} coefficients are fixed) is approximated
as a \ac{GP} problem. In \cite{Bashar2018mixQoS}, the problem of
max-min fairness for a group of users is considered, while imposing
\ac{QoS} constraints on the remaining users. This problem is also
solved using a \ac{GP} approach, which avoids the need of a bisection search as done in \cite{Ngo2018rician}. Other studies in the existing literature
such as \cite{Bashar2018maxmin,Bashar2019maxmin} also exploit \ac{GP}
to solve the power control problem. This is also the case in \cite{Bashar2019EEuplink}
where the problem of the total \ac{EE} maximization, subject to the per-user power and per-user \ac{QoS} constraints, is considered.
The problem of \ac{SE} maximization for cell-free massive \ac{MIMO}
with \ac{NOMA} is in studied \cite{Zhang2019noma}. Again, the considered
problem is solved using the sequential successive convex approximation method in combination
with \ac{GP}.

The main feature of all the above-mentioned studies is that they need to solve a single
geometric program or a series of linear programs to solve the incurred power control problems. It is well
known that these second-order methods requires very high complexity, and thus, are not
preferred for large-scale problems. This is why the existing literature
for the uplink of cell-free massive \ac{MIMO} has been limited to small-scale
scenarios whereby the number of \acp{AP} is around a few hundreds
(less than 250 in all aforementioned papers, to be precise). Thus,
the potential of cell-free massive \ac{MIMO} in large-scale settings has
remained unknown. We aim to fill this gap of the current literature
in this paper.

We assume \ac{MRC} at \acp{AP} using local channel estimates
and consider the problem of max-min fairness subject to power constraint
at each user. Similar to the known approaches, the considered problem
is split into two optimization subproblems. While the linear receiver
coefficient design is a generalized eigenvalue problem which admits
closed-form solution as in previous studies, we propose a more scalable
numerical method that can reduce the run-time to solve the power control
problem significantly. In this regard, our contributions are summarized
as follows:
\begin{itemize}
\item We \emph{reformulate} the power control problem, which is \emph{nonconvex}
in the original form, into a convex program. Due to the nature of
the max-min optimization, the objective of the equivalent convex program
is \emph{nonsmooth}, which can pose some numerical challenges.
\item We propose the use of Nesterov\textquoteright s smoothing technique
\cite{Nesterov2005a} to approximate the objective of the equivalent
convex function. This allows us to apply an \ac{APG} method that solves the approximated problem very fast. In
particular, our proposed method only requires the computation of the
gradient of the objective and closed-form expressions for projections.
\item We avail of the proposed method to explore the performance gains of
cell-free massive \ac{MIMO} in large-scale settings. In particular, we find
that having more \acp{AP} may be more beneficial than more antennas per \ac{AP}, given the number of the total antennas is kept the same.
\end{itemize}
\emph{Notations}: Bold lower and upper case letters represent vectors
and matrices. $\mathcal{CN}(0,a)$ denotes a complex Gaussian random
variable with zero mean and variance $a$. $\mathbf{X}^{\mbox{\scriptsize T}}$ and
$\mathbf{X}^{\mbox{\scriptsize H}}$ stand for the transpose and Hermitian of $\mathbf{X}$,
respectively. $x_{i}$ is the $i$-th entry of vector $\mathbf{x}$;
$[\mathbf{X}]_{i,j}$ is the $(i,j)$-th entry of $\mathbf{X}$. Notation
$\mathbf{e}_{i}$ denotes the $i$-th unit vector, i.e., the
vector such that $e_{i}=1$ and $e_{j}=0,\forall j\neq i$. $\nabla f(\mathbf{x})$
represents the gradient of $f(\mathbf{x})$. $\Vert\cdot\Vert$ denotes
the Euclidean or $l_{2}$ norm; $|\cdot|$ is the absolute value of the argument.

\section{System Model and Problem Formulation}

We consider an uplink cell-free massive \ac{MIMO} scenario where $K$ users
are jointly served by $M$ \acp{AP} and each \ac{AP} is equipped
with $L$ antennas. The users and \acp{AP} are randomly distributed
in a coverage area. We model the channel
coefficients between the $m$-th \ac{AP} and the $k$-th user as
$\mathbf{h}_{mk}=\beta_{mk}^{1/2}\mathbf{g}_{mk}$, where $\beta_{mk}$
is large-scale fading coefficient and $\mathbf{g}_{mk}\in\mathbb{C}^{L\times1}\sim\mathcal{CN}(0,1)$
is the vector of small-scale fading coefficients for all antennas
at the $m$-th \ac{AP}.

\subsection{Uplink Training and Channel Estimation}

We denote the length of the coherence interval and the uplink training
phase in data symbols as $T_{c}$ and  $T_{p}$, respectively. By $\sqrt{T_{p}}\boldsymbol{\psi}_{k}\in\mathbb{C}^{T_{p}\times1}$,
we denote the pilot sequence transmitted from the $k$-th user to
the \ac{AP}, where $\Vert\boldsymbol{\psi}_{k}\Vert^{2}=1$. For $\zeta_{p}$
being the transmit power of each pilot, the $m$-th \ac{AP} receives the
$L\times T_{p}$ matrix given as
\begin{equation}
\mathbf{Y}_{m}^{p}=\sqrt{\zeta_{p}T_{p}}\sum_{k=1}^{K}\mathbf{h}_{mk}\boldsymbol{\psi}_{k}^{\mbox{\scriptsize H}}+\mathbf{N}_{m}^{p},\label{eq:pilotMatrix}
\end{equation}
where $\mathbf{N}_{m}^{p}\in\mathbb{C}^{N\times T_{p}}\sim\mathcal{CN}(0,1)$
is the noise matrix during pilot transmission. Note that we have normalized
the channel coefficients $\mathbf{h}_{mk}$ by the square root of
the true noise power, $\sigma_{N}$, and thus the noise power in \eqref{eq:pilotMatrix}
is unity. To estimate the channel estimate at each \ac{AP}, the matrix
in \eqref{eq:pilotMatrix} is first multiplied with $\boldsymbol{\psi}_{k}$
to get $\mathbf{y}_{mk}^{p}=\mathbf{Y}_{m}^{p}\boldsymbol{\psi}_{k}$,
and then finding the minimum mean-square error (MMSE) estimate as
\cite{Ngo2017cfmm}
\begin{align*}
\hat{\mathbf{h}}_{mk} & =\mathbb{E}\{\mathbf{h}_{mk}\bigl(\mathbf{y}_{mk}^{p}\bigr)^{\mbox{\scriptsize H}}\}\Bigl(\mathbb{E}\{\mathbf{y}_{mk}^{p}\bigl(\mathbf{y}_{mk}^{p}\bigr)^{\mbox{\scriptsize H}}\}\Bigr)^{-1}\mathbf{y}_{mk}^{p}\\
 & =c_{mk}\mathbf{y}_{mk}^{p},
\end{align*}
where
\begin{equation}
c_{mk}=\frac{\sqrt{\zeta_{p}T_{p}}\beta_{mk}}{\zeta_{p}T_{p}\sum_{i=1}^{K}\beta_{mi}\left|\boldsymbol{\psi}_{k}^{\mbox{\scriptsize H}}\boldsymbol{\psi}_{i}\right|^{2}+1}.
\end{equation}
The mean square of any element of $\hat{\mathbf{h}}_{mk}$ is given
by
\begin{equation}
\nu_{mk}=\mathbb{E}\{\left|[\hat{\mathbf{h}}_{mk}]_{n}\right|^{2}\}=\frac{\zeta_{p}T_{p}\beta_{mk}^{2}}{\zeta_{p}T_{p}\sum_{i=1}^{K}\beta_{mi}\left|\boldsymbol{\psi}_{i}^{\mbox{\scriptsize H}}\boldsymbol{\psi}_{k}\right|^{2}+1}.
\end{equation}

\subsection{Uplink Payload Data Transmission}

All the users simultaneously transmit payload data to the $m$-th
\ac{AP} during the payload interval $T_{c}-T_{p}$. (Here the downlink transmission is ignored.) We denote by
$s_{k}$ the symbol to be transmitted from the $k$-th user where
$\mathbb{E}\{|s_{k}|^{2}\}=1$. The transmitted signal from the $k$-th
user is $x_{k}=\sqrt{\eta_{k}}s_{k}$, where $\eta_{k}=\mathbb{E}\{|x_{k}|^{2}\}$
is the power at the $k$-th user. The received signal vector from
all users at the $m$-th \ac{AP} is expressed as
\begin{equation}
\mathbf{y}_{m}^{u}=\sum_{k=1}^{K}\mathbf{h}_{mk}\sqrt{\eta_{k}}s_{k}+\mathbf{n}_{m}^{u},
\end{equation}
where $\mathbf{n}_{m}^{u}\in\mathbb{C}^{L\times1}\sim\mathcal{CN}(0,1)$ is the noise at the $m$-th \ac{AP} during the uplink transmission.
The power at user $k$ must satisfy the following power constraint
\begin{equation}
\eta_{k}\leq\eta^{(\max)},\,\forall k,
\end{equation}
where $\eta^{(\max)}$ is the maximum transmit power. In this paper,
similar to \cite{Bashar2019EEuplink,Bashar2019maxmin}, \ac{MRC} is used at each \ac{AP}. For further improve the
performance of the signal detection at the \ac{CPU}, each \ac{AP} $m$ multiples
the resulting signal with a weighting coefficient $u_{mk}$ for the
$k$-th user. As a result, the aggregated received signal at the \ac{CPU}
is given by
\begin{equation}
r_{k}=\sum_{m=1}^{M}u_{mk}\hat{\mathbf{h}}_{mk}^{\mbox{\scriptsize H}}\mathbf{y}_{m}^{u}.\label{eq:combinedsignal}
\end{equation}

\subsection{Achievable Uplink Spectral Efficiency}
To detect $s_{k}$, the combined signal at the \ac{CPU} in \eqref{eq:combinedsignal}
is split into desired signal, beamforming uncertainty,
inter-user interference, and total noise
for the $k$-th user. From this step, it was shown in \cite{Ngo2017cfmm}
that the following spectral efficiency for the $k$-th user is achievable
\begin{equation}
\mathcal{S}_{ek}(\boldsymbol{\eta},\mathbf{u})=\left(1-\frac{T_{p}}{T_{c}}\right)\log_{2}\bigl(1+\gamma_{k}(\boldsymbol{\eta},\mathbf{u})\bigr),\label{eq:SEk}
\end{equation}
where $\gamma_{k}(\boldsymbol{\eta},\mathbf{u})$ is the corresponding \ac{SINR}
for the $k$-th user defined as
\begin{equation}
\gamma_{k}(\boldsymbol{\eta},\mathbf{u})=\frac{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\bigl(\boldsymbol{\nu}_{kk}\boldsymbol{\nu}_{kk}^{\mbox{\scriptsize H}}\eta_{k}\bigr)\mathbf{u}_{k}}{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\bigl(\sum_{i\neq k}^{K}\boldsymbol{\nu}_{ki}\boldsymbol{\nu}_{ki}^{\mbox{\scriptsize H}}\eta_{i}+\frac{1}{L}\sum_{i=1}^{K}\mathbf{D}_{ki}\eta_{i}+\frac{1}{L}\mathbf{R}_{k}\bigr)\mathbf{u}_{k}},
\end{equation}
where $\boldsymbol{\eta}=[\eta_{1};\eta_{2};\ldots;\eta_{K}]$, $\mathbf{u}_{k}=[u_{1k};u_{2k};\ldots;u_{Mk}]$,
$\mathbf{u}=[\mathbf{u}_{1};\mathbf{u}_{2};\ldots;\mathbf{u}_{K}$],
$\mathbf{D}_{ki}\in\mathbb{R}_{+}^{M\times M}$ is a diagonal matrix
with $[\mathbf{D}_{ki}]_{m,m}=\nu_{mk}\beta_{mi}$, $\mathbf{R}_{k}\in\mathbb{R}_{+}^{M\times M}$
is a diagonal matrix with $[\mathbf{R}_{k}]_{m,m}=\nu_{mk}$, and
\begin{equation}
\boldsymbol{\nu}_{ki}\triangleq\left|\boldsymbol{\psi}_{k}^{\mbox{\scriptsize H}}\boldsymbol{\psi}_{i}\right|\left[\nu_{1k}\frac{\beta_{1i}}{\beta_{1k}};\nu_{2k}\frac{\beta_{2i}}{\beta_{2k}};\ldots;\nu_{Mk}\frac{\beta_{Mi}}{\beta_{Mk}}\right].
\end{equation}

\subsection{Problem Formulation}

To deliver fairness of all users in the system, we consider the maximization
of the minimum \ac{SE} problem, which is mathematically stated as
\begin{equation}
(\mathcal{P}_{1}):\begin{cases}
\underset{\mathbf{u}_{k},\boldsymbol{\eta}}{\maximize} & \underset{1\leq k\leq K}{\min}\mathcal{S}_{ek}(\boldsymbol{\eta},\mathbf{u})\\
\st & \Vert\mathbf{u}_{k}\Vert=1,\,\forall k,\\
 & 0\leq\boldsymbol{\eta}\leq\eta^{(\max)}/L.
\end{cases}\label{eq:maxminProb}
\end{equation}
It is easy to see that the above problem is equivalent to maximizing
the minimum of the \ac{SINR} of individual users, which is given
by

\begin{equation}
(\mathcal{P}_{2}):\begin{cases}
\underset{\mathbf{u}_{k},\boldsymbol{\eta}}{\maximize} & \underset{1\leq k\leq K}{\min}\gamma_{k}(\boldsymbol{\eta},\mathbf{u})\\
\st & \Vert\mathbf{u}_{k}\Vert=1,\,\forall k,\\
 & 0\leq\boldsymbol{\eta}\leq\eta^{(\max)}/L.
\end{cases}\label{eq:equivProb}
\end{equation}
Due to the coupling of $\mathbf{u}_{k}$ and $\boldsymbol{\eta}$
in problem \eqref{eq:maxminProb}, it is quite natural to decouple
problem $(\mathcal{P}_{1})$ into two subproblems: the receiver coefficient
design and the power control problem. We remark that this method has
been widely used in the previous studies. We detail these two problems
in the following section and, in particular, present our novel solution
for the power control problem.

\section{Proposed Solution}

\subsection{Receiver Coefficient Design}

For the receiver coefficient design, we fix the power allocation variable
$\boldsymbol{\eta}$ and solve the following problem to find $\mathbf{u}_{k},\,\forall k$:
\begin{equation}
(\mathcal{P}_{3}):\begin{cases}
\underset{\mathbf{u}_{k}}{\maximize} & \underset{1\leq k\leq K}{\min}\gamma_{k}(\boldsymbol{\eta},\mathbf{u})\\
\st & \Vert\mathbf{u}_{k}\Vert=1,\,\forall k.
\end{cases}\label{eq:RxCoeff}
\end{equation}
This problem can be independently solved for each user to maximize
the \ac{SINR} at each user \cite{Bashar2018mixQoS}. Let us define
$\mathbf{A}_{k}=\boldsymbol{\nu}_{kk}\boldsymbol{\nu}_{kk}^{\mbox{\scriptsize H}}\eta_{k}$
and $\mathbf{B}_{k}=\sum_{i\neq k}^{K}\boldsymbol{\nu}_{ki}\boldsymbol{\nu}_{ki}^{\mbox{\scriptsize H}}\eta_{i}+\frac{1}{L}\sum_{i=1}^{K}\mathbf{D}_{ki}\eta_{i}+\frac{1}{L}\mathbf{R}_{k}$.
Then $(\mathcal{P}_{3})$ is reduced to the following problem for
each $k$:
\begin{equation}
(\mathcal{P}_{4}):\begin{cases}
\underset{\mathbf{u}_{k}}{\maximize} & \frac{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\mathbf{A}_{k}\mathbf{u}_{k}}{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\mathbf{B}_{k}\mathbf{u}_{k}}\\
\st & \Vert\mathbf{u}_{k}\Vert=1.
\end{cases}\label{eq:RxCoeff2}
\end{equation}
We note that problem $(\mathcal{P}_{4})$ is a generalized eigenvalue
problem. According to \cite[Lemma B.10]{Bjornson2017book}, the solution
to $(\mathcal{P}_{4})$ is given by

\begin{equation}
\mathbf{u}_{k}^{*}=\frac{\tilde{\mathbf{u}}_{k}}{\bigl\Vert\tilde{\mathbf{u}}_{k}\bigr\Vert}\label{eq:optimalU}
\end{equation}
 where $\tilde{\mathbf{u}}_{k}=\sqrt{\eta_{k}}\mathbf{B}_{k}^{-1}\boldsymbol{\nu}_{kk}$.
We remark that the receiver coefficient design has been presented
in the previous studies known to us. Our contributions in the paper
are those in handling the power control problem when the receiver
coefficients are held fixed. We provide the details in the next subsection.

\subsection{Power Control Problem}

After updating the receiver coefficients $\mathbf{u}_{k},\,\forall k$,
we need to find the power coefficients to further improve the objective
of $(\mathcal{P}_{3})$, which leads to the following power control
problem:
\begin{equation}
(\mathcal{P}_{5}):\begin{cases}
\underset{\boldsymbol{\eta}}{\maximize} & \underset{1\leq k\leq K}{\min}\gamma_{k}(\boldsymbol{\eta},\mathbf{u})\\
\st & 0\leq\boldsymbol{\eta}\leq\eta^{(\max)}/L.
\end{cases}\label{eq:powerAlloc}
\end{equation}
To solve the above problem, existing methods reformulate it using
an epigraph form, written as
\begin{equation}
(\mathcal{P}_{6}):\begin{cases}
\underset{\boldsymbol{\eta},t}{\maximize} & t\\
\st & 0\leq\boldsymbol{\eta}\leq\eta^{(\max)}/L,\\
 & t\leq\gamma_{k}(\boldsymbol{\eta},\mathbf{u}),\,\forall k.
\end{cases}\label{eq:GPprob}
\end{equation}
There are \emph{two main drawbacks} of such a method. First, the number
of constraints in $(\mathcal{P}_{6})$ has been increased by $K$.
This is not preferred for large-scale problems, i.e. when $K$ is
relatively large. Second, problem $(\mathcal{P}_{6})$ is in fact
a GP which can be solved by off-the-shelf convex solvers. The major
issue of GP is that it requires very high complexity, compared to
other standard convex problems, and thus is not suitable for large
scale problems. To overcome this complexity issue in the current literature,
we propose a more numerical scalable method in what follows.

\subsubsection{Convex Reformulation}

To derive the proposed method, we first reformulate the nonconvex
objective function in $(\mathcal{P}_{5})$ as
\begin{equation}
\maximize\ \underset{1\leq k\leq K}{\min}\frac{1}{\gamma_{k}^{-1}(\boldsymbol{\eta})}\Longleftrightarrow\minimize\ \underset{1\leq k\leq K}{\max}\gamma_{k}^{-1}(\boldsymbol{\eta}),
\end{equation}
where
\begin{align}
\gamma_{k}^{-1}(\boldsymbol{\eta}) & =\frac{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\bigl(\sum_{i\neq k}^{K}\boldsymbol{\nu}_{ki}\boldsymbol{\nu}_{ki}^{\mbox{\scriptsize H}}\eta_{i}+\frac{1}{L}\sum_{i=1}^{K}\mathbf{D}_{ki}\eta_{i}+\frac{1}{L}\mathbf{R}_{k}\bigr)\mathbf{u}_{k}}{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\bigl(\boldsymbol{\nu}_{kk}\boldsymbol{\nu}_{kk}^{\mbox{\scriptsize H}}\eta_{k}\bigr)\mathbf{u}_{k}}\nonumber\\
 & =\eta_{k}^{-1}\left(\sum_{i\neq k}^{K}a_{ki}\eta_{i}+\sum_{i=1}^{K}b_{ki}\eta_{i}+c_{k}\right).\label{eq:posynomialObj}
\end{align}
In the above equation, $a_{ki}$, $b_{ki}$, and $c_{k}$
are defined as
\begin{align}
a_{ki} & =\frac{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\boldsymbol{\nu}_{ki}\boldsymbol{\nu}_{ki}^{\mbox{\scriptsize H}}\mathbf{u}_{k}}{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\boldsymbol{\nu}_{kk}\boldsymbol{\nu}_{kk}^{\mbox{\scriptsize H}}\mathbf{u}_{k}},\
b_{ki} =\frac{(1/L)\mathbf{u}_{k}^{\mbox{\scriptsize H}}\mathbf{D}_{ki}\mathbf{u}_{k}}{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\boldsymbol{\nu}_{kk}\boldsymbol{\nu}_{kk}^{\mbox{\scriptsize H}}\mathbf{u}_{k}},\nonumber\\
c_{k} & =\frac{(1/L)\mathbf{u}_{k}^{\mbox{\scriptsize H}}\mathbf{R}_{k}\mathbf{u}_{k}}{\mathbf{u}_{k}^{\mbox{\scriptsize H}}\boldsymbol{\nu}_{kk}\boldsymbol{\nu}_{kk}^{\mbox{\scriptsize H}}\mathbf{u}_{k}}.
\end{align}
 Note that we have simply written $\gamma_{k}^{-1}(\boldsymbol{\eta})$
instead of $\gamma_{k}^{-1}(\boldsymbol{\eta},\mathbf{u})$ since
$\mathbf{u}$ is fixed. Now we use a change of variables to convert
the posynominal in \eqref{eq:posynomialObj} into a convex function.
Specifically, we define $\theta_{i}=\log\eta_{i}$ or $\eta_{i}=e^{\theta_{i}}$
and reformulate \eqref{eq:posynomialObj} as
\begin{align}
f_{k}(\boldsymbol{\theta}) & =\gamma_{k}^{-1}(\boldsymbol{\theta})=e^{-\theta_{k}}\Bigl(\sum_{i\neq k}^{K}a_{ki}e^{\theta_{i}}+\sum_{i=1}^{K}b_{ki}e^{\theta_{i}}+c_{k}\Bigr)\nonumber \\
  & =\sum_{i\neq k}^{K}a_{ki}e^{(\mathbf{e}_{i}-\mathbf{e}_{k})^{\mbox{\scriptsize T}}\boldsymbol{\theta}}+\sum_{i=1}^{K}b_{ki}e^{(\mathbf{e}_{i}-\mathbf{e}_{k})^{\mbox{\scriptsize T}}\boldsymbol{\theta}}+c_{k}e^{-\mathbf{e}_{k}^{\mbox{\scriptsize T}}\boldsymbol{\theta}}.\label{eq:convexapprox}
\end{align}
Let $f(\boldsymbol{\theta})=\underset{1\leq k\leq K}{\max}f_{k}(\boldsymbol{\theta})$.
Then $(\mathcal{P}_{5})$ can be reformulated as
\begin{equation}
(\mathcal{P}_{7}):\begin{cases}
\underset{\boldsymbol{\theta}}{\minimize} & f(\boldsymbol{\theta})=\underset{1\leq k\leq K}{\max}f_{k}(\boldsymbol{\theta})\\
\st & \boldsymbol{\theta}\leq\theta_{\max},
\end{cases}\tag{\ensuremath{\mathcal{P}_{7}}}\label{eq:convexpowercontrol}
\end{equation}
where $\theta^{(\max)}=\log(\eta^{\max}/L)$. As mentioned earlier,
we can consider an epigraph form of $(\mathcal{P}_{7})$, which can
be solved by convex solvers. But this will increase the complexity.
A numerical difficulty in solving $(\mathcal{P}_{7})$ is due to the
fact that $f(\boldsymbol{\theta})$ is $nonsmooth$. To overcome
this issue, we apply Nesterov's smoothing technique to approximate
the function $f(\boldsymbol{\theta})$ by the following log-sum-exp
function \cite{Nesterov2005a}
\begin{align}
f(\boldsymbol{\theta};\tau) & =\frac{1}{\tau}\log\frac{1}{K}\sum\nolimits _{k=1}^{K}e^{\tau f_{k}(\boldsymbol{\theta})},\label{eq:smoothappx}
\end{align}
where $\tau>0$ is the positive \emph{smoothness} parameter. It was
shown in \cite{Nesterov2005a} that $f(\boldsymbol{\theta})+\frac{\log K}{\tau}\geq f(\boldsymbol{\theta};\tau)\geq f(\boldsymbol{\theta})$.
In other words, $f(\boldsymbol{\theta};\tau)$ is a differentiable
approximation of $f(\boldsymbol{\theta})$ with a numerical accuracy
of $\frac{\log K}{\tau}$. Thus, with a sufficiently high $\tau$,
$f(\boldsymbol{\theta})$ can be replaced by $f(\boldsymbol{\theta};\tau)$
and we can consider the following problem
\begin{equation}
\boxed{(\mathcal{P}_{8}):\begin{cases}
\underset{\boldsymbol{\theta}}{\minimize} & f(\boldsymbol{\theta};\tau)\\
\st & \boldsymbol{\theta}\in\Theta\triangleq\{\boldsymbol{\theta}|\boldsymbol{\theta}\leq\theta_{\max}\}.
\end{cases}}\tag{\ensuremath{\mathcal{P}_{8}}}\label{eq:convexpowercontrol-1}
\end{equation}

\subsubsection{Accelerated Projected Gradient Method}

 We are now in a position to propose an efficient algorithm to solve
$(\mathcal{P}_{8})$, which is essentially an \ac{APG} \cite{Beck2009LIP}. The description of the proposed
algorithm is provided in Algorithm \ref{alg:mAPG}. Note that $\alpha>0$
is called the step size which should be sufficiently small to guarantee
its convergence. Also, the notation $P_{\Theta}(\mathbf{u})$ denotes
the projection onto $\Theta$, i.e., $P_{\Theta}(\mathbf{u})=\arg\min\bigl\{\Vert\mathbf{x}-\mathbf{u}\Vert\ |\ \mathbf{x}\in\Theta\bigr\}$.
\begin{algorithm}[th]
\caption{Proposed \ac{APG} Algorithm for Solving $(\mathcal{P}_{8})$}
\label{alg:mAPG}

\begin{algorithmic}[1]

\STATE Input: $\boldsymbol{\theta}^{0}>0,\alpha>0$

\STATE $\boldsymbol{\theta}^{1}=\boldsymbol{\theta}^{0}$; $t_{0}=t_{1}=1$;
$n=1$

\REPEAT

\STATE $\mathbf{y}^{n}=\boldsymbol{\theta}^{n}+\frac{t_{n-1}-1}{t_{n}}(\boldsymbol{\theta}^{n}-\boldsymbol{\theta}^{n-1})$

\STATE $\boldsymbol{\theta}^{n+1}=P_{\Theta}(\mathbf{y}^{n}-\alpha\nabla f(\mathbf{y}^{n}))$
\label{grady}

\STATE $t_{n+1}=0.5\left(1+\sqrt{4t_{n}^{2}+1}\right)$;

\STATE $n=n+1$

\UNTIL{convergence }

\STATE Output: $\boldsymbol{\theta}^{*}$

\end{algorithmic}
\end{algorithm}

There are two main operations of Algorithm \ref{alg:mAPG}, namely:
finding the gradient of $f(\boldsymbol{\theta};\tau)$, and the projection
onto the feasible set $\Theta$. The details of these two operations
are given in the following.

\subsubsection*{Gradient of $f(\boldsymbol{\theta};\tau)$}

It is easy to see that the gradient of $f(\boldsymbol{\theta};\tau)$
is given by
\begin{equation}
\nabla f(\boldsymbol{\theta};\tau)=\frac{1}{\sum\nolimits _{k=1}^{K}\exp\bigl(\tau f_{k}(\boldsymbol{\theta})\bigr)}\sum\nolimits _{k=1}^{K}e^{\tau f_{k}(\boldsymbol{\theta})}\nabla f_{k}(\boldsymbol{\theta}),\label{eq:gradconvexfunction}
\end{equation}
where $\nabla f_{k}(\boldsymbol{\theta})$ is found as
\begin{align}
\nabla f_{k}(\boldsymbol{\theta}) & =\sum_{i\neq k}^{K}a_{ki}e^{(\mathbf{e}_{i}-\mathbf{e}_{k})^{\mbox{\scriptsize T}}\boldsymbol{\theta}}\bigl(\mathbf{e}_{i}-\mathbf{e}_{k}\bigr)\nonumber \\
 & +\sum_{i=1}^{K}b_{ki}e^{(\mathbf{e}_{i}-\mathbf{e}_{k})^{\mbox{\scriptsize T}}\boldsymbol{\theta}}\bigl(\mathbf{e}_{i}-\mathbf{e}_{k}\bigr)  -c_{k}e^{-\mathbf{e}_{k}^{\mbox{\scriptsize T}}\boldsymbol{\theta}}\mathbf{e}_{k}.\label{eq:gradInvSINR}
\end{align}

\subsubsection*{Projection onto $\Theta$}

For a given $\mathbf{x}\in\mathbb{R}^{K\times1}$, $P_{\Theta}(\mathbf{x})$
is the solution to the following problem:
\begin{multline}
\underset{\boldsymbol{\theta}\in\mathbb{R}^{K\times1}}{\minimize}\ \Bigl\{\Vert\boldsymbol{\theta}-\mathbf{x}\Vert^{2}\ \Bigl|\ \boldsymbol{\theta}\leq\theta^{(\max)}\Bigr\}.
\end{multline}
It is straightforward to check that the solution to the above problem
is $\boldsymbol{\theta}=[\theta_{1};\theta_{2};\ldots;\theta_{K}]$,
where
\begin{equation}
\theta_{k}=\begin{cases}
x_{k}, & x_{k}\leq\theta_{k}\\
\theta^{(\max)}, & x_{k}>\theta_{k}
\end{cases},k=1,\ldots,K.\label{eq:projPolyhedral}
\end{equation}

In summary, combining the receiver coefficient design and the power
control problem, the proposed algorithm for the considered problem
is outlined in Algorithm \ref{alg:minmax-APG}.
\begin{algorithm}[th]
\caption{Proposed \ac{APG}-based Algorithm for Solving $(\mathcal{P}_{1})$}
\label{alg:minmax-APG}

\begin{algorithmic}[1]

\STATE Input: $\boldsymbol{\theta}^{0}>0$; $m=1$

\REPEAT

\STATE For fixed $\boldsymbol{\theta}^{m-1}$, find optimal $\mathbf{u}^{m}=[\mathbf{u}_{1}^{m};\mathbf{u}_{2}^{m};\ldots;\mathbf{u}_{K}^{m}]$
using \eqref{eq:optimalU}\label{line:rxcoeff}

\STATE For fixed $\mathbf{u}^{m}$, find $\boldsymbol{\theta}^{m}=\boldsymbol{\theta}^{*}$
using Algorithm \ref{alg:mAPG}

\STATE $m=m+1$

\UNTIL{convergence}

\STATE Output: $\boldsymbol{\theta}^{m},\mathbf{u}^{m}$

\end{algorithmic}
\end{algorithm}

\subsection{Complexity Analysis}

We now provide the complexity analysis of the proposed algorithm for
one iteration using the big-O notation. It is clear that the complexity
of Algorithm \ref{alg:minmax-APG} is dominated by the computation
of four quantities: the receiver coefficients $\mathbf{u}^{m}$, the
smooth approximate objective function $f(\boldsymbol{\theta};\tau)$,
the gradient of $f(\boldsymbol{\theta};\tau)$, and the projection
onto $\Theta$. The computation of each $\mathbf{u}_{k}$
requires the calculation of the inverse of $\mathbf{B}_{k}$. Normally
this step requires the complexity of $\mathcal{O}(M^{3})$. However,
by exploiting the specific structure of $\mathbf{B}_{k},$ in the Appendix
we show that the complexity of this step is reduced to $\mathcal{O}\bigl((K-1)M^{2}\bigr).$
This \emph{massive computation saving when $M$ is large} has not
been reported in the existing literature. Accordingly, the complexity
for updating the receiver coefficient $\mathbf{u}^{m}$ is $\mathcal{O}(K^{2}M^{2})$.
It is easy to see that $KM$ multiplications are required to compute
$\mathcal{S}_{ek}(\boldsymbol{\eta})$. Therefore, the complexity
of finding $f(\boldsymbol{\eta})$ is $\mathcal{O}(K^{2}M).$ Similarly,
we can find that the complexity of $\nabla f_{k}(\boldsymbol{\theta})$
is $\mathcal{O}(K^{2}M)$. The projection of $\boldsymbol{\eta}$
onto $\Theta$ is given in \eqref{eq:projPolyhedral}, and thus the
complexity of the projection is $\mathcal{O}(K)$. In summary, the
per-iteration complexity of the proposed algorithm for solving $(\mathcal{P}_{4})$
is $\mathcal{O}(K^{2}M^{2})$.

\section{Numerical Results}

We evaluate the performance of our proposed method in terms of the
achievable rate and the time complexity. In the considered cell-free
massive \ac{MIMO} system, we randomly distribute the \acp{AP} and
the users over a $D\times D$ \si{\km\squared}. The large-scale fading
coefficient between the $m$-th AP and the $k$-th user is generated
as $\bar{\beta}_{mk}=\mathrm{PL}_{mk}.z_{mk},$ where $\mathrm{PL}_{mk}$
is the path loss between the $m$-th \ac{AP} and the $k$-th user,
and $z_{mk}$ represents the log-normal shadowing between the $m$-th
\ac{AP} and the $k$-th user with mean zero and standard deviation
$\sigma_{\textrm{sh}}$, respectively. In this paper, we adopt the
three-slope path loss model and model parameters as in \cite{Farooq2020}.
The noise power is given by $B\times k_{B}\times T_{0}\times W$,
where $B=20$ MHz denotes the bandwidth, $k_{B}=1.381\times10^{-23}$ (Joule/Kelvin)
represents the Boltzmann constant, $T_{0}=290$ (Kelvin) denotes the
noise temperature, and $W=9$ dB shows the noise figure. The length
of the coherence time and the uplink training phase are set to $T_{c}=200$, $T_{p}=20$, respectively. If not otherwise mentioned, we set the power transmit
power for downlink data transmission and uplink training phase (before
normalization) as $\zeta_{u}=0.2$ W and $\zeta_{p}=0.2$ W.

In the first numerical experiment, we compare Algorithm \ref{alg:minmax-APG}
with the method in \cite{Bashar2018mixQoS,Bashar2019maxmin} which
solves \eqref{eq:maxminProb} using \ac{GP} for the power control
problem. We refer to this method as the GP-based method in this section.
For this method, we use the convex solver MOSEK \cite{MOSEKApS15}
through the modeling tool CVX \cite{cvx} in the GP mode to
solve the power control problem in $(\mathcal{P}_{1})$. Figure \ref{Fig:GPconvergence}
shows the convergence of Algorithm \ref{alg:minmax-APG} and the \ac{GP}-based
method to solve $(\mathcal{P}_{1})$ for three randomly generated
channel realizations.
\begin{figure}[tbh]
\centering\includegraphics[bb=62bp 552bp 291bp 738bp]{Fig2_3.eps} \caption{Convergence of Algorithm \ref{alg:minmax-APG} and the \ac{GP}-based
algorithm for three random channel realizations. The relevant parameters
are taken as $M=150,K=20,L=1,D=1$.}
\label{Fig:GPconvergence}
\end{figure}

We can see that both methods reach almost the same objective for all
three considered channel realizations. The marginal difference is
due to a moderate value of the smoothness parameter $\tau$ in \eqref{eq:smoothappx}.
The main advantage of our proposed method over the \ac{GP}-based
method is that each iteration of the proposed method is very memory
efficient and computationally cheap, and hence, can be executed very
fast. As a result, the total run-time of the proposed method is far
less than that of the \ac{GP}-based method as shown in Table~\ref{table: Table 1}.
In Table~\ref{table: Table 1}, we report the actual run-time of
both methods to solve the max-min problem. Here, we run our codes
on a 64-bit Windows operating system with 16 GB RAM and Intel CORE
i7, 3.7 GHz. Both iterative methods are terminated when the difference
of the objective for the latest two iterations is less than $10^{-5}$.
\begin{table}[tbh]
\caption{Comparison of run-time (in seconds) between Algorithm \ref{alg:minmax-APG}
and the \ac{GP}-based method. Here, $K=20$ and $D=1$.}
\label{table: Table 1} \centering{}%
\begin{tabular}{c|c|c}
\hline
APs & \ac{GP}-based Method & Proposed Method\tabularnewline
\hline
120 & 3.67 & \textbf{0.80}\tabularnewline
\hline
160 & 3.79 & \textbf{1.21}\tabularnewline
\hline
200 & 4.08 & \textbf{1.86}\tabularnewline
\hline
240 & 4.42 & \textbf{2.45}\tabularnewline
\hline
\end{tabular}
\end{table}

In the next experiment, we plot the \ac{CDF} of the per-user \ac{SE} (bit/s/Hz) for $300$ channel realizations shown in Fig. \ref{Fig:CDFplot}.
We consider two large-scale scenarios: (i) $K=25,M=500$ and (ii)
$K=50,M=1000$. Note that the ratio between the number of APs to the
number of users is the same for two cases. It can be seen in Fig.
\ref{Fig:CDFplot} that as the number of the users increases, the
per-use \ac{SE} slightly decreases. We also observe that the per-user \ac{SE} differs in the range of 0.1 bit/s/Hz which means that
the fairness is indeed achieved among the users.

\begin{figure}[H]
\centering\includegraphics[bb=62bp 552bp 291bp 738bp]{Fig5.eps}\caption{\ac{CDF} of per-user achievable \ac{SE} for two scenarios $K=25,M=500$
and $K=50,M=1000$. The number of antennas at each AP is $L=1$.}
\label{Fig:CDFplot}
\end{figure}

Finally, we investigate the effect of increasing the number of antennas
per \ac{AP} on the minimum achievable \ac{SE}. Specifically, we plot
the minimum achievable \ac{SE} with respect to the number of antennas
for both $M=250$ and $M=500$ \acp{AP}. The number of users is fixed to
$K=50$. As expected, the minimum achievable \ac{SE} increases with the
number of antennas per \ac{AP} but the increase tends to be small when
the number of antennas is sufficiently large. The reason is that for
a large number of \acp{AP} channel harderning and favorable propagation
can be achieved by a few antennas per \ac{AP}. Specifically, we can see
that the \ac{SE} for the case of $500$ \acp{AP} and $5$ antennas per \ac{AP}
is larger than the \ac{SE} for the case of 250 \acp{AP} and 10 antennas per \ac{AP}. Therefore, for large-scale cell-free massive MIMO having more
APs with a few antennas each tends to be more beneficial than having
fewer APs with more antennas each.
\begin{figure}[H]
\includegraphics[bb=62bp 550bp 294bp 738bp]{Fig4.eps}\caption{Maximized minimum achievable \ac{SE} with respect to the number of antenna
at each AP for $K=50,D=1.$}
\label{Fig:RatevsAntennas}
\end{figure}

\section{Conclusion}

In this paper, we have proposed a low-complexity method for maximizing
the minimum achievable \ac{SE} in the uplink of the cell free massive
\ac{MIMO} subject to per-user power constraint, assuming the \ac{MRC}
technique at \acp{AP}. As in previous studies, we have divided the min-max
problem into two subproblems: the receiver coefficient design and
the power control problem. For the power control problem, the existing
solutions rely on \ac{GP} which requires high complexity. To overcome
this issue, we have reformulated the power control problem into a
convex form and then proposed an efficient numerical algorithm based
on Nesterov's smoothing technique and \ac{APG}.
Our proposed solution only requires the first order information of
the objective and the projection, both of which are given in the closed-form.
We have numerically shown that our proposed solution practically achieves
the same objective as the \ac{GP}-based method but in much lesser
time. We have used the proposed method to study the performance of
the large-scale systems up to 1000 \acp{AP} for which the known methods
are not suitable. We have also shown that a large-scale cell-free
massive \ac{MIMO} having more \acp{AP} with a few antennas has better
performance than a similar system having a fewer \acp{AP} with more
antennas per \ac{AP}.