Distributed Caching for Data Dissemination in the Downlink of Heterogeneous Networks
Jun Li Youjia Chen Zihuai Lin Wen Chen Branka Vucetic Lajos Hanzo
eess.SP cs.IT math.IT

theoremTheorem
propositionProposition
definitionDefinition
lemmaLemma
corollaryCorollary
remarkRemark
constructionConstruction
algorithmAlgorithm

















IEEE Transactions on Communications, Accepted for Publication
Li et al.: Distributed Caching for Data Dissemination in the Downlink of Heterogeneous Networks

Distributed Caching for Data Dissemination in the Downlink of Heterogeneous Networks
    Jun Li Member, IEEE, Youjia Chen, Zihuai Lin Senior Member, IEEE,
 Wen Chen Senior Member, IEEE, Branka Vucetic Fellow, IEEE, and Lajos Hanzo Fellow, IEEE
Jun Li is with the School of Electronic and Optical Engineering, Nanjing University of Science and Technology, Nanjing, 210094, CHINA. Email: jleesr80@gmail.com.
Youjia Chen, Zihuai Lin, and Branka Vucetic are with School of Electrical and Information Engineering, The University of Sydney, NSW, 2006, AUSTRALIA. E-mail: {youjia.chen,zihuai.lin,branka.vucetic}@sydney.edu.au.
Wen Chen is with the Department of Electronic Engineering, Shanghai Jiao Tong University, CHINA. E-mail: wenchen@sjtu.edu.cn.
Lajos Hanzo is with  Department of Electronics and Computer Science, University of Southampton, U.K. E-mail: lh@ecs.soton.ac.uk.
This work was supported by Australian Research Council Programs DP120100405, by National 973 Project 2012CB316106, by NSF China 61328101, 61271230 and 61472190, by STCSM Science and Technology Innovation Program 13510711200, and by SEU National Key Lab on Mobile Communications 2013D11 and 2013D02.
    version of December 30, 2023
================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



Heterogeneous cellular networks (HCN) with embedded small cells are considered, where multiple mobile users wish to download network content of different popularity. By caching data into the small-cell base stations (SBS), we will design distributed caching optimization algorithms via belief propagation (BP) for minimizing the downloading latency. First, we derive the delay-minimization objective function (OF) and formulate an optimization problem. Then we develop a framework for modeling the underlying HCN topology with the aid of a factor graph. Furthermore, distributed BP algorithm is proposed based on the network's factor graph. Next, we prove that a fixed point of convergence exists for our distributed BP algorithm. In order to reduce the complexity of the BP, we propose a heuristic BP algorithm. Furthermore, we evaluate the average downloading performance of our HCN for different numbers and locations of the base stations (BS) and mobile users (MU), with the aid of stochastic geometry theory. By modeling the nodes distributions using a Poisson point process, we develop the expressions of the average factor graph degree distribution, as well as an upper bound of the outage probability for random caching schemes. We also improve the performance of random caching. Our simulations show that (1) the proposed distributed BP algorithm has a near-optimal delay performance, approaching that of the high-complexity exhaustive search method, (2) the modified BP offers a good delay performance at a low communication complexity, (3) both the average degree distribution and the outage upper bound analysis relying on stochastic geometry match well with our Monte-Carlo simulations, and (4) the optimization based on the upper bound provides both a better outage and a better delay performance than the benchmarks.



Wireless caching, heterogeneous cellular networks, belief propagation, stochastic geometry





§ INTRODUCTION

Wireless data traffic is expected to increase by a factor of 40 over the next five years, from the current level of 93 Petabytes to 3600 Petabytes per month <cit.>, driven by a rapid increase in the number of mobile users (MU) and aggravated by their bandwidth-hungry mobile applications. A promising approach to enhancing the network capacity is to embed small cells relying on low-power base stations (BS) into the existing macro-cell based networks. These networks, which are referred to as heterogeneous cellular networks (HCN) <cit.>, typically contain regularly deployed macro-cells and embedded femto-cells as well as pico-cells <cit.> that are served by macro-cell BSs (MBS) and small-cell BSs (SBS), respectively. The aim of these flexibly deployed low-power SBSs is to eliminate the coverage holes and to increase the capacity in hot-spots.

There is evidence that the MUs' downloading of video on-demand files is the main reason for the growth of data traffic over cellular networks <cit.>. According to the prediction of Cisco on mobile data traffic, the mobile video streaming traffic will occupy 72% percentage of the overall mobile data traffic by 2019. Often, there are numerous repetitive downloading requests of popular contents, such as online blockbusters, leading to redundant data streaming. The redundancy of data transmissions can be reduced by locally storing popular data, known as caching, into the local SBSs, effectively forming a local cloud caching system (LCCS). The LCCS brings the content closer to the MUs and alleviates redundant data transmissions via redirecting the downloading requests to local SBSs. Also, the SBSs are willing to cache files into their buffers as long as they can, since caching is capable of significantly reducing the tele-traffic load on their back-haul channels, which are expensive.

In <cit.>, the authors study the caching strategies of delay-tolerant vehicular networks, where the data subscribers and “helpers” are always moving and the links between them are opportunistic. By proposing an efficient algorithm to carefully allocate the network resources to mobile data, the decision is made as to which content should use the erasure coding, as well as conceiving the coding policy for each mobile data. In <cit.>, optimal cache replacement policies are investigated. The cache replacement process takes place after the data caching process has been completed, and determines which particular data item should be deleted from the cache, when the available storage space is insufficient for accommodating an item to be cached.

Since the HCN structure has been widely adopted in current cellular networks and will prevail in near-future networks, we are interested in the SBS-based LCCS in the context of HCNs. In contrast to the vehicular networks discussed in <cit.>, where the mobility and the opportunistic communication contact are important issues, in the context of HCNs, the BSs are always fixed, and the MUs are assumed to be moving at a low speed. Thus, we ignore the mobility issues in the HCNs and assume that each MU is associated with a fixed BS during file-downloading. At the time of writing, there are already technical reports highlighting the advantages of caching in HCNs <cit.>. Based on these reports, the LCCS with SBS caching for HCNs is capable of efficiently 1) reducing the transmission latency due to short distance between the SBSs and the MUs, 2) offloading redundant data streams from MBSs, and 3) alleviating heavy burdens on the back-haul channels of the SBSs. Therefore, SBS-based caching will bring about significant breakthroughs for future HCNs.

The concept of caching is common in wireline networks and computer systems. However, research on efficient caching design for wireless cellular networks relying on small cells is still in its infancy <cit.>. Usually, data caching consists of two phases: data placement and data transmission. During the data placement phase, data is cached into local SBSs in order to form an LCCS. In the data transmission phase, MUs request data from the LCCS. The focus of wireless caching research is mainly on the optimization of data placement for ensuring that the downloading latency is minimized. The caching optimization is a non-trivial problem. This is due to the massive scale of video contents to be stored in the limited memory of the SBSs.

The survey papers <cit.> report on a range of attractive caching architectures conceived for future cellular networks. In <cit.>, a caching scheme is proposed for a device-to-device (D2D) based cellular network on the MUs' caching of popular data. In this scheme, the D2D cluster size was optimized for reducing the downloading delay. In <cit.>, the authors propose a caching scheme for wireless sensor networks, where the protocol model of <cit.> is adopted. In <cit.>, a femto-caching scheme is proposed for a cellular network combined with SBSs, where the data placement at the SBSs is optimized in a centralized manner for reducing the transmission delay imposed. However, <cit.> considers an idealized system, where neither the interference nor the impact of wireless channels is taken into account. The associations between the MUs and the SBSs are pre-determined without considering the specific channel conditions encountered. Furthermore, this centralized optimization method assumes that the MBS has perfect knowledge of all the channel state information (CSI) between the MUs and SBSs, which is impractical.

Against this background, in this paper, we consider distributed caching solutions for HCNs operating under more practical considerations. Our contributions consist of two parts.


In the first part, we propose distributed caching algorithms for enhancing the downloading performance via belief propagation (BP) <cit.>. The BP algorithm is capable of decomposing a global optimization problem into multiple sub-problems, thereby offering an efficient distributive approach of solving the global optimization problem <cit.>. As the BP method has been widely adopted for distributively solving resource allocation in cellular networks, we arrange file placement via BP algorithms by viewing files as a type of resource.


In the second part, we analyze the average caching performance based on stochastic geometry theory <cit.>. We are interested in optimizing the average performance of a set of HCNs, where the channels exhibit Rayleigh fading and the distributions of network nodes obey a Poisson point process (PPP) <cit.>.



Specifically, our contributions in the first part are follows.


  * We commence by deriving the delay as our optimization objective function (OF) and formulate the problem as optimizing the file placement.

  * We develop a framework for modeling the associated factor graph based on the topology of the network. A distributed BP algorithm is proposed based on the factor graph, which allows the file placement to be optimized in a distributed manner between the MUs and SBSs.

  * We prove that a fixed point exists in the proposed BP algorithm and show that the BP algorithm is capable of converging to this fixed point under certain conditions.

  * To reduce the communication complexity, we propose a heuristic BP algorithm.


Our contributions in the second part are follows.


  * By following the stochastic geometry framework, we model the MUs and SBSs in the HCN as different ties of a PPP. Furthermore, we develop the average degree distribution of the factor graph in the BP algorithm.

  * A random caching scheme is proposed, where each SBS will cache a file with a pre-determined probability. We can characterize the average downloading performance by outage probability (OP) and develop a tight upper bound of the OP expression with a closed form under the random caching scheme.

  * Based on the upper bound derived, we further improve the OP performance of random caching by optimizing the probabilities for caching different files.


In the simulations, we first investigate the average degree distribution of the factor graph, as well as the OP and the delay of the random caching schemes, in conjunction with various PPP parameters and power settings. It is shown that both the degree distribution and our upper bound analysis match well with the results of Monte-Carlo simulations. Furthermore, the optimization based on the upper bound provides both a better OP and a better delay than the benchmarks. Then we evaluate the distributed BP algorithm in our HCNs having a fixed number of BSs and MUs. It is shown that the proposed distributed BP algorithm has a near-optimal performance, approaching that of the exhaustive search method. The heuristic BP also offers a relatively good performance, despite its significantly reduced communication complexity.

The rest of this paper is organized as follows. We describe the system model in Section <ref> and present the distributed file downloading problem relying on caching in Section <ref>. We then propose a distributed BP algorithm in Section <ref>, where the proof of existence for a fixed point is also presented. In Section <ref>, a heuristic BP algorithm is proposed for reducing the associated communication complexity. Our stochastic geometry based analysis is detailed in Section <ref>, where the average degree distribution of the factor graph and the OP of the random caching scheme are developed. Our simulation results are summarized in Section <ref>, while our conclusions are provided in Section <ref>.



§ SYSTEM MODEL

Let us consider an HCN consisting of a single MBS and K SBSs illuminating both femto-cells and pico-cells, while supporting J MUs randomly located in the network. Let us denote by ℬ_0 the MBS and by ℬ={ℬ_1,ℬ_2,⋯,ℬ_K} the set of the SBSs, where ℬ_k, k∈𝒦={1,2,⋯,K}, represents the k-th SBS. Furthermore, denote by 𝒰={𝒰_1,𝒰_2,⋯ 𝒰_J} the set of the MUs, where 𝒰_j, j∈𝒥={1,2,⋯,J}, represents the j-th MU. The MBS ℬ_0 caches files into the memories of the SBSs during off-peak time via back-haul channels. Once the caching process is completed, the MBSs and SBSs are ready to act upon the downloading requests of the MUs.

We assume that a dedicated frequency band of bandwidth W is allocated to the downlink channels spanning from the SBSs to the MUs for file-dissemination. For reasons of careful load balancing, we consider the “SBS-first” constraint, where each MU will try to download data from its adjacent SBSs, unless the required files cannot be found in these SBSs. In this case, the MU will turn to the MBS for retrieving the required files. For the sake of simplicity, we assume that the MBS will support a fixed download rate, denoted by C_0, for the MUs in the channels which are orthogonal to those spanning from the SBSs to MUs.

In order to satisfy the “SBS-first” constraint for offloading data from the MBS, some incentives may be provided for the MUs. For example, downloading from the SBSs is much cheaper than from the MBS. Here, we assume that the download rate C_0 supported by the MBS is never higher than the lowest download rate supported by the SBSs. This limit imposed on the download rate from the MBS will not only encourage the MUs to download from the SBSs first, but also effectively control the data traffic of the MBS imposed by file downloading.

Denote by P_k the transmission power of the k-th SBS, and by σ^2 the noise power at each MU. The path-loss between ℬ_k and the MU 𝒰_j is modeled as d_k,j^-α, where d_k,j is the distance between ℬ_k and 𝒰_j, and α is the path-loss exponent. The random channel between ℬ_k and 𝒰_j is Rayleigh fading, whose coefficient h_k,j has the average power of one. We assume that all the downlink channels spanning from the SBSs to the MUs are independent and identically distributed (i.i.d.).

Suppose that each file is split into multiple chunks and each chunk can be downloaded by an MU in a short time slot. Due to the short downloading time of a chunk, we assume furthermore that the probability of having two MUs streaming a chunk at the same time (or within a relative delay of a few seconds) from the same SBS is basically zero <cit.>. Hence, neither direct multicasting by exploiting the broadcast nature of the wireless medium nor network coding is considered. Furthermore, we focus our attention on the saturated scenario, where the SBSs keep transmitting data to the MUs <cit.>. Hence, each MU is subject to the interference imposed by all the other SBSs in ℬ, when downloading files from its associated SBS. Given a channel realization h_j=[h_1,j,⋯,h_K,j], the channel capacity between ℬ_k and 𝒰_j can be calculated based on the signal-to-interference-plus-noise ratio (SINR) as

    C_k,j=Wlog(1+h_k,j^2d_k,j^-αP_k/q∈𝒦\{k}∑h_q,j^2d_q,j^-αP_q+σ^2).

Due to the `SBS-first' constraint, we have C_0≤ C_k,j, ∀ k∈𝒦,j∈𝒥.

Denote by ℱ the library or set of files, which consists of Q popular files to be requested frequently by the MUs. The popularity distribution among the set ℱ is represented by 𝒫={p_1,p_2,⋯,p_Q}, where the MUs make independent requests of the f-th file, f=1,⋯,Q, with the probability of p_f. Without any loss of generality, all these files have the same size of M bits. We assume that ℬ_0 has a sufficiently large memory and hence accommodates the entire library of files, while the storage of each SBS is limited to G files, where we have G<Q.

Without a loss of generality, we assume that Q/G is an integer. The Q files in ℱ are divided into N=Q/G file groups (FG), with each FG containing G files. The f-th file, ∀ f∈{(n-1)G+1,⋯,nG}, is included in the n-th FG, n∈𝒩={1,⋯,N}. We denote by ℱ_n the n-th FG, and by P_ℱ_n the probability that the MUs request a file in ℱ_n. Based on 𝒫, we have

    P_ℱ_n=∑_f=(n-1)G+1^nGp_f.

File caching is then carried out on the basis of FG, i.e., each SBS caches one of the N FGs.



§ DISTRIBUTED FILE DOWNLOADING RELYING ON CACHING

The caching-based distributed file downloading protocol consists of two stages. The first stage, or file placement stage, includes file content broadcasting and caching. In this stage, ℬ_0 broadcasts the FGs to the SBSs via the back-haul during off-peak periods. At the same time, the SBSs listen to the broadcasting from ℬ_0, and cache the FGs needed. The second stage, or file downloading stage, includes MU-SBS associations and file content transmissions. In this stage, each MU makes decisions as to which SBSs it should be associated with, and then starts to download files from the associated SBSs. When the requested files are not found in the adjacent SBSs, the MUs will turn to the MBS for these files.



 §.§ File Placement Matrix

For assigning the N FGs to the K SBSs, we set up a file placement matrix Λ of size K× N. The entry λ_k,n∈{0,1} in Λ indicates whether ℱ_n is cached by ℬ_k or not. We have λ_k,n=1 if ℱ_n is cached by ℬ_k, while λ_k,n=0 otherwise. The k-th row of Λ indicates which FG is cached by ℬ_k, and the n-th column indicates which BS caches ℱ_n. The number of the SBSs which cache ℱ_n can be calculated as ∑_k∈𝒦λ_k,n. Since each SBS caches one FG, we have ∑_n∈𝒩λ_k,n=1.



 §.§ MU-SBS Association

Denote by ℋ(j) the subscript set of the specific SBSs, which are capable of providing a sufficiently high SINR for the MU 𝒰_j. The SBSs in ℋ(j) are the candidates for 𝒰_j to be potentially associated with. By setting an SINR threshold δ, ℬ_k will be included in ℋ(j) if and only if

    h_k,j^2d_k,j^-αP_k/q∈𝒦\{k}∑h_q,j^2d_q,j^-αP_q+σ^2≥δ.


When requesting a file in ℱ_n, 𝒰_j first communicates with one of the SBSs in ℋ(j) which caches ℱ_n. It is possible that more than one SBS in ℋ(j) caches ℱ_n. In this case, 𝒰_j will associates with the optimal SBS, which imposes the minimum downloading delay.

It is clear that the downloading delay is inversely proportional to the downlink transmission rate. According to the file request assumption stipulated in the previous section, there is only a single MU connected to an SBS at each time. Thus, the maximum transmission rate from ℬ_h to 𝒰_j, ∀ h∈ℋ(j), is the channel capacity between them, i.e., C_h,j. When 𝒰_j tries to download a file in ℱ_n, it follows the maximum-capacity association criterion. Hence, 𝒰_j associates with B_ĥ such that

    ĥ=h∈ℋ(j)max{λ_h,nC_h,j}.


When none of the SBSs in ℋ(j) caches ℱ_n, i.e., we have λ_h,n=0, ∀ h∈ℋ(j), 𝒰_j will associate with the MBS for the requested file.



 §.§ Optimization Problem Formulation

We now optimize the matrix Λ for minimizing the average delay of downloading a file. Only when the optimal Λ has been determined will the file-placement stage commence, where the files are placed according this optimal matrix. Once the MU-SBS associations have been determined, we can optimize the matrix Λ for minimizing the average delay of downloading a file. First, given the channel coefficients and the specific location of 𝒰_j, the delay of downloading a file in ℱ_n by 𝒰_j can be calculated as

    D_j,n=M/max_h∈ℋ(j){λ_h,nC_h,j},   ∃λ_h,n≠ 0, ∀ h∈ℋ(j)
    M/C_0,   otherwise.

Based on the request probability of each FG, the delay for 𝒰_j to download a file from ℱ can be written as D_j=∑_n∈𝒩P_ℱ_nD_j,n. Thus, the average delay for each MU can be calculated as

    D=1/J∑_j∈𝒥D_j.


By setting D as the OF, let us hence formulate the delay optimization problem as follows:

    minimize      D
    s.t.     ∑_n∈𝒩λ_k,n=1, ∀  k∈𝒦,
         Λ∈{0,1}^K× N.

The optimization problem in (<ref>) is an integer programming problem, which is NP-complete. In <cit.>, similar optimization problems have been solved by sub-optimal solutions, such as the classic greedy algorithm (GA). However, the existing solutions are typically based on centralized optimization. As we can see from (<ref>), a centralized minimization of D at ℬ_0 requires the global CSI between ℬ and 𝒰, which is impractical. Hence, we will dispense with this assumption and optimize Λ in a distributed manner at a low complexity.




§ DISTRIBUTED BELIEF PROPAGATION ALGORITHM

In this section, we propose a distributed algorithm based on BP for solving the optimization problem of (<ref>) as follows: 1) We first develop a factor graph for describing the message passing in the BP algorithm. 2) Then we map the resultant factor graph to the network for the sake of facilitating the distributed BP optimization. 3) This solved by solving our optimization problem by proposing a distributed BP algorithm. 4) Finally, the proof of existence for a fixed point of convergence in the BP algorithm is presented.



 §.§ Factor Graph Model

In our BP algorithm, the factor graph has to be first established based on the underlying network as a standard bipartite graphical representation of a mathematical relationship between the local delay functions and file allocation variables. Then the BP algorithm is implemented by iteratively passing messages between the local functions and their related variables. Our optimization problem is thus solved by the proposed BP algorithm based on the factor graph.

Based on the topology of the HCN, we develop a factor graph model 𝒢=(𝒱,ℰ), where 𝒱 is the vertex set, and ℰ is the edge set. The vertex set 𝒱 consists of factor nodes and variable nodes. Each factor node is related to an MU and each variable node is related to an SBS. To simplify the notations, we denote by j∈𝒥 the j-th factor node and denote by k∈𝒦 the k-th variable node. Hence, the vertex set 𝒱 is composed of 𝒥 and 𝒦, i.e., 𝒱={𝒥,𝒦}.

As mentioned in the previous section, ℬ_k will be a candidate for 𝒰_j to potentially associate with, but only if the received SINR at 𝒰_j from ℬ_k is no less than the threshold δ. Correspondingly, in our factor graph, an edge in the edge set ℰ connecting 𝒰_j and ℬ_k, denoted by (j,k), exists if the received SINR at 𝒰_j from ℬ_k is no less than δ. The node k is named as a neighboring node of j, if there is an edge (j,k). Actually, ℋ(j) defined previously represents the set of the neighboring nodes of the factor node j. Furthermore, denote by ℋ(k) the set of neighboring node for the variable node k. Fig. <ref> illustrates a factor graph extracted from an HCN with 5 SBSs and 10 MUs. Take ℬ_1 in the factor graph for example. The edges exist between ℬ_1 and 𝒰_2 as well as 𝒰_4, which means that ℬ_1 can provide a sufficient large SINR for both 𝒰_2 and 𝒰_4.

The distributed BP algorithm is based on the factor graph 𝒢. The factor nodes in 𝒥 represent the local utility functions generated from the decomposition results of the global utility function, which will be discussed later in this subsection. The variable nodes in 𝒦 represent the variables to be optimized, i.e., the entries of Λ. The factor nodes and variable nodes are connected by edges in ℰ, indicating the message flows in the BP algorithm. That is, messages are only passing between a node and its neighbors. We now illustrate the optimization problem on the factor graph.



  §.§.§ Factor Nodes

According to Eq. (<ref>), the OF can be decomposed into J local contributions as D_1,⋯, D_J. These local contributions are calculated based on Eq. (<ref>). Since the BP algorithm solves maximization problems, we define a series of utility functions as F≜ -D and F_j≜ -D_j. Then our optimization problem can be rewritten as

    max_ΛF(Λ),    F= 1/J∑_j∈𝒥F_j.

We use the j-th factor node to represent the j-th local utility function F_j, which is related to 𝒰_j. Hence, the maximization of F can be achieved by maximizing F_j at 𝒰_j, ∀ j∈𝒥.


  §.§.§ Variable Nodes

Each variable node is related to an SBS. Here, we use the k-th variable node to represent the k-th row of Λ, denoted by λ_k, which is related to ℬ_k. The location of `1' in λ_k indicates which specific FG is stored by ℬ_k. Note that the first constraint in (<ref>) means that each SBS only stores a single FG. Given this constraint, λ_k has N possible values according to N different locations of `1'. We denote by λ_k^[1], ⋯, λ_k^[N] the N values of λ_k. When we have λ_k=λ_k^[n], this implies that the FG ℱ_n is stored by ℬ_k. Take N=2 for example, where λ_k=λ_k^[1]=[1 0] indicates that the FG ℱ_1 is stored in the SBS ℬ_k, while λ_k=λ_k^[2]=[0 1] indicates that ℱ_2 is stored in ℬ_k. The variables λ_k, k=1,⋯,K, are the parameters to be optimized for maximizing F in (<ref>). For simplicity, we use the matrix Λ to represent the set of the variables λ_k in the factor graph.



 §.§ Distributed Belief Propagation

In standard BP, the variables are optimized by estimating their marginal probability distributions <cit.>. Note that the utility function F is a function of the file placement matrix Λ. We define the probability mass function (PMF) p(Λ) of Λ based on the utility function F(Λ) as

    p(Λ)≜1/Zexp(μ F(Λ)),

where μ is a positive number and Z is the normalization factor. According to <cit.>, the result of large deviations shows that when μ→∞, p(Λ) concentrates around the maxima of F(Λ), i.e., lim_μ→∞𝔼(Λ)=ΛmaxF(Λ), where 𝔼(Λ) is the expectation of Λ. Once we obtain 𝔼(Λ), we can have a good estimate of the specific Λ which maximizes F(Λ).

In our distributed BP, the maximization of F can be decomposed into J maximization operations on F_j at 𝒰_j, j=1,⋯, J. Correspondingly, the estimation of Λ is decomposed into J estimations of its subsets Λ_j at 𝒰_j, where Λ_j={λ_h,∀ h∈ℋ(j)}. The PMF of Λ_j is written as p_j(Λ_j)=1/Z_jexp(μ F_j(Λ_j)), where Z_j is the normalization factor. Since all the variables are independent, the estimation of Λ_j at 𝒰_j can be further decomposed into the estimation of each individual λ_h via calculating its PMF p_j(λ_h), which is the marginal PMF of p_j(Λ_j) with respect to the variable λ_h. Hence we have p_j(λ_h)=𝔼_∼λ_h(p_j(Λ_j)), where 𝔼_∼λ_h(·) represents the expectation over the elements in Λ_j, except for λ_h. The PMF p_j(λ_h) is viewed as the message, which is iteratively updated between 𝒰_j and ℬ_h, ∀ h∈ℋ(j). The PMF p_j(λ_h) consists of N probabilities estimated by 𝒰_j, i.e., (λ_h=λ_h^[1]),⋯,(λ_h=λ_h^[N]), where (λ_h=λ_h^[n]) represents the probability that ℱ_n is stored by ℬ_h.

Without a loss of generality, we assume that the edge (j,k) does exist in the factor graph. We represent the iteration index by t and denote by p_k→ j^(t)(λ_k) and p_j→ k^(t)(λ_k) the belief messages emanated from ℬ_k to 𝒰_j and from 𝒰_j to ℬ_k during the t-th iteration, respectively. The steps describing the distributed BP are as follows.



  §.§.§ Initialization

At the variable nodes, set t=1 and let p_k→ j^(1)(λ_k) to be the initial distribution of λ_k, e.g., the a priori popularity distribution 𝒫.



  §.§.§ Variable Node Update

During the t-th iteration, each SBS ℬ_k updates the message p_k→ j^(t)(λ_k) to be sent to 𝒰_j based on the messages gleaned from ℬ_k's neighboring MUs other than 𝒰_j in the previous iteration. This includes the calculations of N probabilities. Given λ_k=λ_k^[n], ∀ n∈𝒩, we have

    p_k→ j^(t)(λ_k^[n])=1/Z_k∏_ħ∈ℋ(k)\{j}p_ħ→ k^(t-1)(λ_k^[n]),

where Z_k is the normalization factor so that we have ∑_n∈𝒩p_k→ j^(t)(λ_k^[n])=1.



  §.§.§ Factor Node Update

In the t-th iteration, 𝒰_j updates the N probabilities of the message p_j→ k^(t)(λ_k) to be sent to ℬ_k, which is based on the messages received from 𝒰_j's neighboring SBSs, except for ℬ_k. The messages updated at the factor nodes are calculated according to the marginal PMF. Given λ_k=λ_k^[n], ∀ n∈𝒩, we have

    p_j→ k^(t)(λ_k^[n])
       =𝔼_∼λ_k(exp(μ F_j(λ_k^[n],{λ_h,∀ h∈ℋ(j)\{k}})))
       
    =∑_h∈ℋ(j)\{k}∑_λ_h=λ_h^[1]^λ_h^[N](∏_q∈ℋ(j)\{k}
    p_q→ j^(t)(λ_q)·
           exp(μ F_j(λ_k^[n],{λ_h,∀ h∈ℋ(j)\{k}}))).




  §.§.§ Final Solution

Let us assume that there are t=T iterations in the distributed BP algorithm. After T iterations, the probability that ℱ_n is stored by ℬ_k can be obtained by

    (λ_k=λ_k^[n])=1/Z_k∏_ħ∈ℋ(k)p_ħ→ k^(T)(λ_k^[n]).

Based on (<ref>), the decision as to which file should be stored by ℬ_k can be made by choosing the specific file that has the maximum a posteriori probability (λ_k=λ_k^[n]), ∀ n∈𝒩.



 §.§ Convergence to a Fixed Point

Let us now investigate the existence of a fixed point of convergence in our distributed BP algorithm. The essence of the distributed BP algorithm is to keep updating the PMF p_j(λ_k) before reaching its final estimate. Based on (<ref>) and (<ref>), the evolution of p_j(λ_k) during the t-th iteration can be obtained from the PMFs in the (t-1)-th iteration as

    p_k→ j^(t)(λ_k)=   1/Z_k∏_ħ∈ℋ(k)\{j}∑_h∈ℋ(ħ)\{k}∑_λ_h=λ_h^[1]^λ_h^[N]
       (exp(μ F_ħ(Λ_ħ))·∏_q∈ℋ(ħ)\{k}
    p_q→ħ^(t-1)(λ_q)).

We view the PMF p_k→ j^(t)(λ_k) as a probability vector of length N. We define the probability vector set ℳ^(t)≜{p_k→ j^(t)(λ_k)} for all k∈𝒦 as well as j∈𝒥, and define the message mapping function Γ:ℝ^N× KJ→ℝ^N× KJ based on (<ref>) so that ℳ^(t)=Γ(ℳ^(t-1)). Then we have the following lemma.

The message mapping function Γ is a continuous mapping.

Proof: Please refer to Appendix <ref>.


Given Lemma <ref>, we have the following theorem.

A fixed point of convergence exists for the proposed distributed BP algorithm.

Proof: Please refer to Appendix <ref>.


The question of convergence to the fixed point is, unfortunately, not well understood in general <cit.>. Generally, if the factor graph contains no cycles, the belief propagation can be shown to converge to a fixed solution point in a finite number of iterations. The performance, including the optimality and the convergence rate, of the BP crucially depends on the choice of the objective function, as well as the scale, the sparsity and the number of cycles in the underlying factor graph. As such, the theoretical analysis of the BP algorithm's optimality and convergence rate remains an open challenge.



§ A HEURISTIC BP WITH REDUCED COMPLEXITY

In the context of the BP algorithm, the message p_j(λ_k) exchanged between 𝒰_j and ℬ_k in each iteration, includes N probability values, which are real numbers. Hence, the communication overhead of the message passing is relatively high. Hence, we propose a heuristic BP (HBP) algorithm for reducing the communication overhead imposed. The rationale behind the term “heuristic BP” is that we still follow the classic concept of belief propagation, but use a different format of the beliefs from the conventional one.

Assuming that the edge (j,k) exists, in the t-th iteration of the HBP, instead of forwarding the N probabilities stored in p^(t)_j→ k(λ_k) to ℬ_k, 𝒰_j randomly selects an FG according to these N probabilities. Then the integer index n^(t)_j→ k of the FG selected will be forwarded to the SBS ℬ_k.

At the SBS side, the SBS ℬ_k receives |ℋ(k)| integers, i.e., n^(t)_ħ→ k, ∀ħ∈ℋ(k), from its neighboring MUs, where |·| denotes the cardinality of a set. Based on n^(t)_ħ→ k, the SBS ℬ_k infers the number of those MUs, which indicate that ℱ_n should be stored in the SBS ℬ_k, for n=1,⋯, N. Let us assume now that in the t-th iteration, there are J_k,n^(t) MUs specifically indicating that ℱ_n should be stored in ℬ_k, where we have ∑_n∈𝒩 J_k,n^(t)=|ℋ(k)|. We can view J_k,n^(t)/|ℋ(k)| as the probability that the specific FG ℱ_n is stored by the SBS ℬ_k.

In this case, the probability p_k→ j^(t)(λ_k^[n]) in (<ref>) will be recalculated as

    p_k→ j^(t)(λ_k^[n])=J_k,n^(t-1)-1/|ℋ(k)|-1,   if   n=n^(t-1)_j→ k,
    J_k,n^(t-1)/|ℋ(k)|-1,   if   n≠ n^(t-1)_j→ k.

Note that in (<ref>), the information n^(t-1)_j→ k transmitted from the MU 𝒰_j to the SBS ℬ_k is excluded when calculating p_k→ j^(t)(λ_k^[n]), for the sake of ensuring that only uncorrelated information is exchanged throughout the HBP.

At the MU side, it is clear that the MU 𝒰_j has to obtain p_k→ j^(t)(λ_k^[n]) for the sake of updating the output information. However, there is no need for the SBS ℬ_k to transmit the N probabilities p_k→ j^(t)(λ_k^[n]) to each of its neighboring MUs. Alternatively, ℬ_k broadcasts the N integers, J_k,1^(t),⋯,J_k,N^(t) to the neighboring MUs for reducing the transmission overhead. After receiving the N integers from the SBS ℬ_k, the MU 𝒰_j calculates p_k→ j^(t)(λ_k^[n]) in (<ref>).

Based on the above discussions, the HBP algorithm can be summarized as follows.



  §.§.§ Initialization

At the variable nodes, we set t=1. The SBS ℬ_k randomly generates |ℋ(k)| independent integers, n_1,⋯, n_|ℋ(k)|, according to the popularity distribution 𝒫. These integers are viewed as the indexes of the FGs. We then set J_n,k^(1) to be the number of the integers that are equal to n.



  §.§.§ Variable Node Update

In the t-th iteration, ℬ_k updates and broadcasts the N integers J_n,k^(t), for n=1,⋯, N, to the neighboring MUs. The resulting calculations performed on these N integers J_n,k^(t) are based on the integers n^(t-1)_ħ→ k, ∀ħ∈ℋ(k), received from the neighboring MUs during the last iteration. Specifically, the n-th integer J_n,k^(t) is obtained by counting the number of n^(t-1)_ħ→ k that are equal to n.



  §.§.§ Factor Node Update

The MU 𝒰_j first calculates the probabilities p_h→ j^(t)(λ_k^[n]), ∀ h∈ℋ(j) according to Eq. (<ref>) based on the integers gleaned from the SBS ℬ_h. Then based on p_h→ j^(t)(λ_k^[n]), ∀ h∈ℋ(j)\{k}, 𝒰_j calculates p_j→ k^(t)(λ_k^[n]) according to Eq. (<ref>). After obtaining the N probabilities p_j→ k^(t)(λ_k^[n]), n=1,⋯, N, 𝒰_j randomly chooses an FG according to these N probabilities and sends the index n_j→ k^(t) of the FG to the SBS ℬ_k.



  §.§.§ Final Solution

After T iterations, the SBS ℬ_k makes the decision that the FG ℱ_n̂ should be stored for ensuring that

    n̂=n∈𝒩maxJ_k,n^(T).


The overhead of the HBP is significantly lower than that of the original BP introduced in the previous section. From a communication complexity perspective, in each iteration of the HBP, an SBS ℬ_k broadcasts N integers, while an MU 𝒰_j transmits |ℋ(j)| integers. On the other hand, in the original BP, ℬ_k transmits N|ℋ(k)| real numbers, while 𝒰_j transmits N|ℋ(j)| real numbers for each iteration. From a computational complexity perspective, in a single iteration of the HBP, the computational complexity is on the order of O(N) at the SBS ℬ_k, and O(|ℋ(j)|N^|ℋ(j)|) at the MU 𝒰_j. On the other hand, in the original BP, the computational complexity is O(N|ℋ(k)|^2) at ℬ_k, and O(|ℋ(j)|N^|ℋ(j)|) at 𝒰_j for each iteration.



§ PERFORMANCE ANALYSIS BASED ON STOCHASTIC GEOMETRY

In this section, we analyze both the average degree distribution of the factor graph and the average downloading performance based on stochastic geometry theory. We model the distribution of the MUs as a PPP Φ_U having the intensity of λ_U, and that of the SBSs as an independent PPP Φ_B with the intensity λ_B <cit.>. For simplicity, we assume that all the SBSs have the same transmission power P. In the following, both the degree distribution and the downloading performance are averaged over both the channels' fading coefficients and over the PPP distributions of the nodes.



 §.§ Average Degree Distributions of the Factor Graph

Let us now investigate the degree distribution of the factor graph averaged over PPP. Note that the degree of a factor node j is defined as the number of its neighboring variable nodes, given by the cardinality |ℋ(j)|, while the degree of a variable node k is defined as the number of its neighboring factor nodes, i.e., |ℋ(k)|. Then we have the following theorem.

The factor nodes in the factor graph have the average degree

    ζ_U=2πλ_BZ(λ_B, P, α, δ),

and the variable nodes have the average degree

    ζ_B=2πλ_UZ(λ_B, P, α, δ),

where we have

    Z(λ_B, P, α, δ)=
    ∫^∞_0exp{-2λ_Bπ/αδ ^2/αB(2/α,1-2/α)r^2-δσ^2/Pr^α}rdr

and the Beta function B(x,y)=∫^1_0t^x-1(1-t)^y-1dt.

Proof: Please refer to Appendix <ref>.


When neglecting the noise, we have the following corollary based on Theorem <ref>.

When neglecting the noise, Z(λ_B, P, α, δ) in (<ref>) can be rewritten as

    Z(λ_B, P, α, δ)=α/4πλ_B B(2/α,1-2/α) δ^2/α.

Then we can simplify the average degree of the factor nodes in Eq. (<ref>) to

    ζ_U=α/2 δ^2/α B(2/α, 1-2/α),

and the average degree of the variable nodes in Eq. (<ref>) to

    ζ_B=λ_Uα/2 λ_Bδ^2/α B(2/α, 1-2/α).


Proof: Please refer to Appendix <ref>.


Equations (<ref>) and (<ref>) can be seen as approximations of (<ref>) and (<ref>), respectively, when the effects of the noise are neglected. These approximations are significantly accurate for the HCN, since the interference effects are dominant due to the dense deployments of the SBSs.

From (<ref>), we can see that ζ_U is only related to δ and α, but is independent of λ_U, P and λ_B. In other words, the factor node degree has no relation with the intensities of the MUs and SBSs or with the power of the SBSs. The intuitive reason is that although increasing both the PPP intensities and the power of the SBSs can increase the total signal power, the interference also increases at the same time, which keeps the degree ζ_U of the factor nodes constant. Similarly, observe from (<ref>) that ζ_B is independent of the power P, i.e., increasing the transmission power of the SBSs will not influence the average degree distribution of the factor graph.


We observe that B(2/α, 1-2/α)=π when α=4. Thus, we have closed-form expressions for ζ_U and ζ_B in (<ref>) and (<ref>), respectively, when α=4.




 §.§ Downloading Performance of Random Caching

Since the performance of BP based caching remains difficult for mathematical analysis in closed form, we propose a random caching scheme and analyze its performance based on stochastic geometry theory. The random caching is realized by randomly picking out Ω_ℱ_n· K (0≤Ω_ℱ_n≤1) SBSs from the entire set of K SBSs for caching the FG ℱ_n.

To evaluate the downloading performance, we first define an outage 𝒬_n as the event of an MU's failing to find the FG ℱ_n in its neighboring SBSs. The following theorem states an upper bound of the OP of 𝒬_n. As mentioned before, since the interference is the dominant factor predetermining the network performance, we ignore the noise effects in the following performance analysis to simplify our derivations.

The OP for downloading a file in ℱ_n can be upper-bounded by

    (𝒬_n)≤C(δ, α)(1-Ω_ℱ_n)+A(δ,α)Ω_ℱ_n/C(δ, α)(1-Ω_ℱ_n)+A(δ,α)Ω_ℱ_n+Ω_ℱ_n,

where we have C(δ,α)≜2/αδ^2/αB(2/α,1-2/α), A(δ,α)≜2δ/α-2 _2F_1(1,1-2/α;2-2/α;-δ), and _2F_1 represents the hypergeometric function.

Proof: Please refer to Appendix <ref>.


When the path-loss exponent α=4, we have C(δ,4)=√(δ)/2π and A(δ,4)=δ _2F_1(1,1/2;3/2,-δ). It becomes clear from (<ref>) that (𝒬_n) is only related to δ and Ω_ℱ_n, where a higher δ leads to a higher (𝒬_n). This is because a larger δ will reduce the number of possibly eligible serving SBSs, resulting in an increase of OP. We can see that a higher Ω_ℱ_n leads to a lower (𝒬_n).

Let us define the averaged OP 𝒬 over all the files. Based on the file popularity, the OP of 𝒬 can be upper-bounded by

    (𝒬)=∑_n∈𝒩P_ℱ_n(𝒬_n)
    ≤∑_n∈𝒩P_ℱ_n(C(δ, α)(1-Ω_ℱ_n)+A(δ,α)Ω_ℱ_n)/C(δ, α)(1-Ω_ℱ_n)+A(δ,α)Ω_ℱ_n+Ω_ℱ_n.

The average delay D̅ of each MU can be obtained based on the average OP, i.e.,

    D̅=(1-(𝒬))D̅_s+(𝒬)M/C_0,

where D̅_s is the average delay of downloading from the SBSs. The delay D̅ can be seen as the average value of D in Eq. (<ref>) over both the PPP and the channel fading. Note that D̅_s is usually challenging to calculate and does not have a closed form in the PPP analysis.

Next, we optimize Ω_ℱ_n for improving the downloading performance. Since we do not have a closed-form expression for D̅, we minimize the upper bound of (𝒬) in (<ref>), i.e.,

    max_{Ω_ℱ_n}   ∑_n∈𝒩P_ℱ_nΩ_ℱ_n/Ω_ℱ_n(A(δ,α)-C(δ,α)+1)+C(δ,α),
    s.t.   ∑_n∈𝒩Ω_ℱ_n=1, 
       Ω_ℱ_n≥ 0.

By relying on the classic Lagrangian multiplier, we arrive at the optimal solution as

    Ω^⋆_ℱ_n=max{√(P_ℱ_n/ξ)-C(δ,α)/A(δ,α)-C(δ,α)+1,0},

where ξ=(∑^n^*_q=1√(P_ℱ_q))^2/(n^*C(δ,α_s)+A(δ,α_s)-C(δ,α_s)+1)^2, and n^* satisfies the constraint that Ω_ℱ_n≥ 0.








§ SIMULATION RESULTS

In this section, we first focus on the HCNs associated with PPP distributed nodes, where we investigate the average degree distribution of the factor graph and the performance of the random caching scheme. Then we consider an HCN supporting a fixed number of nodes. We investigate the delay optimized by the BP algorithm and compare it to other benchmarks, including both the random caching and the optimal scheme using exhaustive search.

Note that the physical layer parameters in our simulations, such as the path-loss exponent, noise power, transmit power of the SBSs, and the intensity of the SBSs, are chosen to be practical and in line with the values set by 3GPP standards. For instance, the transmit power of an SBS is typically 2 Watt in 3GPP. The unit of power, such as noise power and transmit power, is the classic Watt. The intensities of the SBSs and MUs are expressed in terms of the numbers of the nodes per square kilometer. Unless specified otherwise, we set the path loss to α=4, the number of files to Q=100, transmit power to P=2, and the noise power to σ^2=10^-10. All the simulations are executed with Matlab. Also, we consider the performance averaged over a thousand network cases, where the locations of network nodes are uniformly distributed in each case, and randomly changed from case to case.



 §.§ Average Degree Distributions of Factor Graph

We compare our Monte-Carlo simulations and analytical results in the HCNs at various transmission powers and node densities. Fig. <ref> shows the average degree of the factor nodes with different transmission power P, SBSs' intensity λ_B, and MUs' intensity λ_U. We can see that for a given δ, the degree ζ_U remains unaffected by the specific choice of P, λ_B, and λ_U. Observe that our analytical results are consistent with the simulations. Similarly, Fig. <ref> shows the average degree of the variable nodes of different powers and node intensities, demonstrating that the results are independent of the power P, but depend on the densities λ_B and λ_U. We can also see that the analytical results match well with the simulation results.






 §.§ Average Downloading Performance of Random Caching

Let us now evaluate the average downloading performance of the random caching scheme supporting PPP distributed nodes. The file distribution 𝒫={p_1,⋯,p_Q} is modeled by the Zipf distribution <cit.>, which can be expressed as

    p_f=1/f^s/∑_q=1^Q1/q^s,  for f=1,⋯,Q,

where the exponent 0<s≤1 is a real number, and it characterizes the popularity of files. Explicitly, a larger s corresponds to a higher content reuse, i.e., the most popular files account for the majority of requests. Note that P_ℱ_n can be obtained based on p_f via Eq. (<ref>).




For the simulation results of this subsection, we assume that each SBS caches G=5 files, hence there are N=Q/G=20 FGs. We commence by considering the OP. In our optimized random caching (ORC), we set Ω_ℱ_n as in (<ref>). For comparison, we also consider another random caching scheme from <cit.> as our the benchmark, namely, the file popularity based random caching (FPRC). In the FPRC, Ω_ℱ_n is chosen to be consistent with the file popularity, i.e., we have Ω_ℱ_n=P_ℱ_n.

Fig. <ref> shows the OPs (𝒬_n)· P_ℱ_n for individual FGs under both the ORC and the FPRC schemes, where we have δ=0.03 and s=0.5. The conditional OP (𝒬_n) (given a file in ℱ_n is requested) is calculated from Eq. (<ref>), while the request probability P_ℱ_n of ℱ_n is calculated from Eq. (<ref>). The FGs are arranged in descending order of popularity, i.e., the first FG has the highest popularity, while the last one has the lowest popularity. We can see from the figure that compared to the FPRC, FGs having a higher popularity have a lower OP, while the ones with lower popularity have higher OPs in the ORC. For example, the OP for the most popular FG is around 0.054 in the ORC in contrast to 0.099 in the FPRC, while the probability of the least popular FG is 0.27 in the ORC in contrast to 0.25 in the FPRC. This is because the ORC is reminiscent of the classic water-filling, allocating more SBSs for caching the higher popular FGs for ensuring the minimization of the average OP.

Let us now investigate the average OP (𝒬). Fig. <ref> and Fig. <ref> show (𝒬) for different δ and s values, respectively. In Fig. <ref>, we fix s=0.5, while in Fig. <ref>, we fix δ=0.03. The dashed lines with different marks are based on the simulations associated with various power and densities, while the solid lines represent the analytical upper bounds of Eq. (<ref>). We can see that the average OP is independent of both the power P and densities λ_B and λ_U. The ORC scheme has a lower average OP than the FPRC. Furthermore, as expected, a higher SINR threshold δ leads to a higher OP, as shown in Fig. <ref>. At the same time, it is interesting to observe from Fig. <ref> that a larger s, representing more imbalanced downloading requests on the different files, can dramatically reduce the OP. We can see that the upper bounds evaluated from Eq. (<ref>) match the simulations quite accurately.

Next, we consider the average delay D̅ in Eq. (<ref>), where we assume an SINR threshold of δ=0.03, a bandwidth of W=10^7 Hz, and a file size of M=10^9 bits. Since C_0 should be always less than the maximum possible downloading rate provided by the SBSs, we assume C_0=Wlog(1+δ). For δ=0.03, C_0 becomes 4.26× 10^5 bits/sec. Fig. <ref> illustrates the average downloading delay associated with different s values. We can see that the ORC scheme always outperforms the FPRC scheme, and that their performance gap becomes larger upon increasing s. Again, the observed performance does not depend on the powers and intensities of the nodes.





 §.§ Delay Performance of Distributed BP algorithms

Let us now study the delay performance of distributed BP-based optimizations. We consider HCNs having fixed numbers of SBSs and MUs, where the locations of these nodes are time-variant. We first consider a small network, in which the optimal solution is found with the aid of an exhaustive search. This will allow us to characterize the performance disparity between the proposed BP algorithm and the optimal search-based solution. Then we focus our attention on a larger network to show the robustness of our BP algorithms. In both scenarios, we set the SINR threshold to δ=0.1, the transmission power to P=2, the bandwidth to W=10^7 Hz, and the file size to M=10^9 bits. Similar to the previous subsection, we assume that the rate provided by the MBS as C_0=Wlog(1+δ). For δ=0.1, we have C_0 as 1.3× 10^6 bits/sec.

In the first scenario, the nodes are arranged in a 0.6×0.6km^2 area using 8 SBSs and 4 MUs. We assume that each SBS caches G=25 files, and there are N=Q/G=4 FGs. Fig. <ref> shows the average delay performance under various schemes, where `HBP' is the heuristic BP algorithm proposed in Section <ref>, `BP' is the original BP algorithm proposed in Section <ref>, and `Optimal' is the optimal scheme relying on an exhaustive search. We can see from Fig. <ref> that the original BP approaches the optimal scheme within a small delay margin. The proposed HBP performs slightly worse than the original BP, with a relatively modest delay degradation of around 5% or 20 seconds, while it outperforms the ORC scheme by about 10% or 40 seconds gain. The FPRC performs the worst among all the caching schemes, exhibiting a substantial delay gap between the FPRC scheme and the ORC scheme.

In the second scenario, the nodes are arranged in a 1.5×1.5km^2 area with 50 SBSs and 25 MUs. We consider two cases, namely Case1 and Case2. In Case1, we assume that each SBS caches G=20 files and there are N=Q/G=5 FGs, while in Case2, we assume that each SBS caches G=10 files and that we have N=Q/G=10. Fig. <ref> shows the average delay performance under various schemes. It is clear from Fig. <ref> that in both cases the BP algorithm performs the best, while the FPRC performs the worst. The HBP exhibits a tiny delay increase of around 3% performance loss compared to the original BP, although it dramatically reduces the communication complexity during the optimization process.

Note also in Fig. <ref> that the ORC suffers from a 5% performance loss compared to the HBP, but it is much less complex than the HBP and BP. The optimization in ORC is based on the statistical information available about both of channels and the locations of the nodes, while both the BP and the HBP exploit the relevant instantaneous information at a relatively high communication complexity. In this sense, the ORC constitutes an efficient caching scheme. Furthermore, we can see from Fig. <ref> that there is a tradeoff between the storage and delay , i.e., a larger storage at each SBS in Case1 leads to a lower downloading delays compared to Case2.

In the above BP simulations, we set the maximum number of iterations to T=15. Table <ref> shows the average number of iterations under different s values for the two scenarios. We can see that the HBP relies on more iterations than the BP. Nevertheless, the overall communication complexity of the HBP is still lower than that of the BP, as we have discussed in Section <ref>. Explicitly, for each iteration of the HBP, ℬ_k broadcasts N integers and 𝒰_j transmits |ℋ(j)| integers. By contrast, in the original BP, ℬ_k transmits N|ℋ(k)| real numbers and 𝒰_j transmits N|ℋ(j)| real numbers.



 §.§ Delay Performance in a Large Scale Network

Finally, we consider a large-scale network associated with Q=1000 files, 50 SBSs, and 100 MUs within an area of 5× 5 km^2. Furthermore, we consider a lower connection probability to the SBSs by setting δ=0.2. By assuming that each SBS is capable of caching 20 files, we have overall 50 file groups. Fig. <ref> shows the average delay performance. We can see from the figure that both BP algorithms perform better than the random caching schemes. Particularly, the HBP has a roughly 1% performance loss compared to the original BP, which imposes however a much reduced communication complexity. This implies that our BP algorithms are robust in large-scale networks associated with a large number of files and network nodes.

Further comparing Figs. <ref>, <ref>, and <ref>, it is interesting to observe that the gap between our BP and HBP algorithms becomes smaller when the network scale becomes larger. More particularly in Fig. <ref>, the performance of these two schemes almost overlaps. This indicate that in large scale networks, we may consider to use the HBP rather than BP to obtain a good performance at a much reduced complexity.



§ CONCLUSIONS

In this paper, we designed distributed caching optimization algorithms with the aid of BP for minimizing the downloading latency in HCNs. Specifically, a distributed BP algorithm was proposed based on the factor graph according to the network structure. We demonstrated that a fixed point of convergence exists for the distributed BP algorithm. Furthermore, we proposed a modified heuristic BP algorithm for further reducing the complexity. To have a better understanding of the average network performance under varying numbers and locations of the network nodes, we involved stochastic geometry theory in our performance analysis. Specifically, we developed the average degree distribution of the factor graph, as well as an upper bound of the OP for random caching schemes. The performance of the random caching was also optimized based on the upper bound derived. Simulations showed that the proposed distributed BP algorithm approaches the optimal performance of the exhaustive search within a small margin, while the modified BP offers a good performance at a very low complexity. Additionally, the average performance obtained by stochastic geometry analysis matches well with our Monte-Carlo simulations, and the optimization based on the upper bound derived provides a better performance than the benchmark of <cit.>.



§ PROOF OF LEMMA <REF>

To simplify the notation in the proof, we assume that ℋ(j)=𝒦, ∀ j∈𝒥 and ℋ(k)=𝒥, ∀ k∈𝒦. Consider a pair of probability vector sets ℳ^(t-1)={p_k→ j^(t-1)(λ_k)} and ℳ^(t-1)={p̃_k→ j^(t-1)(λ_k)}. Then we have the supremum norm

    ||Γ(ℳ^(t-1))-Γ(ℳ^(t-1))||_sup
       
    =max_k,j,n|p_k→ j^(t)(λ_k^[n])-p̃_k→ j^(t)(λ_k^[n])|
    
       =max_k,j,n|∏_i∈𝒥\{j}∑_h∈𝒦\{k}∑_λ_h=λ_h^[1]^λ_h^[N](exp(μ F_i(Λ_i))
    (∏_q∈𝒦\{k}...
         ...p_q→ i^(t-1)(λ_q)
    -∏_q∈𝒦\{k}p̃_q→ i^(t-1)(λ_q)))|
    
       (a)≤max_j∏_i∈𝒥\{j}∑_h∈𝒦\{k}∑_λ_h=λ_h^[1]^λ_h^[N]
         |
    ∏_q∈𝒦\{k}p_q→ i^(t-1)(λ_q)
    -∏_q∈𝒦\{k}p̃_q→ i^(t-1)(λ_q)|
    
       (b)≤(K-1)N^K-1max_j
         ∏_i∈𝒥\{j}max_q∈𝒦\{k},n|
    p_q→ i^(t-1)(λ_q^[n])-p̃_q→ i^(t-1)(λ_q^[n])|
    
       ≤(K-1)N^K-1max_j,q∈𝒦\{k},n|
    p_q→ i^(t-1)(λ_q^[n])-p̃_q→ i^(t-1)(λ_q^[n])|^J-1
       ≤(K-1)N^K-1max_j,k,n|
    p_k→ i^(t-1)(λ_k^[n])-p̃_k→ i^(t-1)(λ_k^[n])|
    
       =(K-1)N^K-1||ℳ^(t-1)-ℳ^(t-1)||_sup.


The inequality (a) in (<ref>) is derived by exploiting the following two facts: 1) 0<exp(μ F_i(Λ))≤ 1, since F_i(Λ) is non-positive and μ is positive, and 2) ∑_s|x_s|≤|∑_s(x_s)| for arbitrary x_s. The inequality (b) in (<ref>) can be obtained from: 1) the following lemma, and 2) the fact that ∑_h∈𝒦\{k}∑_λ_h=λ_h^[1]^λ_h^[N] has to carry out the additions of N^K-1 items.

Given 0≤ a_1,⋯,a_K≤1 and 0≤ã_1,⋯,ã_K≤1, we have

    max_k∈𝒦|∏_q∈𝒦\{k} a_q-∏_q∈𝒦\{k}ã_q|≤ (K-1)max_q∈𝒦\{k}|a_q-ã_q|.


Proof: Please refer to Appendix <ref>.


From (<ref>), we can infer that Γ is a continuous mapping, since the coefficient (K-1)N^K-1 is a constant, and this completes the proof.
▪



§ PROOF OF THEOREM <REF>

Let 𝒮 be the collection of the message set ℳ^(t). The mapping function Θ maps 𝒮 to 𝒮 with the aid of the function Γ. According to Lemma <ref>, Θ is continuous since Γ is continuous. Furthermore, it is clear that the set 𝒮 is convex, closed and bounded. Based on Schauder's fixed point theorem, Θ has a fixed point. This completes the proof.▪



§ PROOF OF THEOREM <REF>



 §.§ The Average Degree of Factor Nodes

Without a loss of generality, we carry out the analysis for a typical MU located at the origin and assume that the potential serving SBSs are located at the point x_B. The fading (power) is denoted by h_x_B, which is assumed to be exponentially distributed, i.e., we have h_x_B∼exp(1). The path-loss function is given by ‖ x_B‖^-α, where ‖·‖ denotes the Euclidian distance.

The average degree of a factor node in the factor graph is equivalent to the number of SBSs that can provide a high enough SINR (≥δ) for the typical MU, which can be formulated as

    N_B=∫_ℝ^2λ_B(ρ(x_B)≥δ)dx_B,

where ρ(x_B) represents the SINR at the typical MU received from the SBSs located at x_B.

We first focus on the probability (ρ(x_B)≥δ) in (<ref>) as follows.

    (ρ(x_B)≥δ)   =(Ph_x_B‖ x_B‖ ^-α/x_k∈Φ_B∑Ph_x_k‖ x_k‖ ^-α+σ^2≥δ)  
       =(h_x_B≥δ(I+σ^2)/P‖ x_B‖ ^-α) 
       =𝔼_I(exp(-sI))exp(-sσ^2),

where x_k denotes the location of an interfering SBS, I=x_k∈Φ_B∑Ph_x_k‖ x_k‖ ^-α represents the aggregate interference, and s=δ‖ x_B‖ ^α/P. The last step is due to the exponential distribution of h_x_B. Then, we derive 𝔼_I(exp(-sI)) in (<ref>) as

    𝔼_I(exp(-sI))(a)=
       𝔼_Φ_B(x_k∈Φ_B∏∫_0^∞exp(-sPh_x_k‖ x_k‖ ^-α)exp(-h_x_k)dh_x_k)   
       (b)=exp(-λ_B∫_ℝ^2(1-1/1+sP‖ x_k‖ ^-α)dx_k)  
       
    =exp( -2πλ_B1/α(sP)^2/αB(2/α,1-2/α)),

where (a) is based on the independence of channel fading, and (b) follows from 𝔼(x∏u(x))=exp(-λ∫_ℝ^2(1-u(x))dx), where x∈Φ and Φ is an PPP in ℝ^2 with the intensity λ <cit.>.

Based on the derivation above, the average degree of the typical MU can be calculated as

    N_B=λ_B∫_ℝ^2
       exp(-2πλ_B/αδ^2/αB(2/α,1-2/α)‖ x_B‖ ^2-δσ^2/P‖ x_B‖ ^α)dx_B
        =2πλ_B∫^∞_0
       exp(-2πλ_B/αδ^2/αB(2/α,1-2/α)r^2-δσ^2/Pr^α)rdr.




 §.§ The Average Degree of Variable Nodes

In this subsection, we consider a typical SBS which is located at the origin, and assume that an MU is located at the point x_U. The average degree of a variable node in the factor graph is equivalent to the number of MUs that can receive at a high enough SINR (≥δ) from the typical SBS, which can be formulated as

    N_U=∫_ℝ^2λ_U(ρ(x_U)≥δ)dx_U,

where ρ(x_U) represents the received SINR at the MU located at x_U from the typical SBS, i.e.,

    (ρ(x_U)≥δ)
    =(Ph_x_U‖ x_U‖ ^-α/x_k∈Φ_B∑Ph_x_k‖ x_k-x_U‖ ^-α+σ^2≥δ),

where x_k denotes the location of an interfering SBS.

Since the PPP is a stationary process, the distribution of ‖ x_k-x_U‖ is independent of the value of x_U, i.e., we have p(‖ x_k-x_U‖)=p(‖ x_k‖), where p(·) represents the probability density function. Then, we have similar results to Eq. (<ref>). That is, we have

    (ρ(x_U)>δ)=𝔼_I(exp(-sI))exp(-sσ^2),

where s=δ‖ x_U‖ ^α/P. Then we arrive at

    N_U=2πλ_U
    ∫^∞_0exp(-2πλ_B/αδ ^2/αB(2/α,1-2/α)r^2-δσ^2/Pr^α)rdr.

By combining Eqs. (<ref>) and (<ref>), we complete the proof.▪



§ PROOF OF COROLLARY <REF>

When ignoring the noise, we have

    Z(λ_B, P, α, δ)=∫^∞_0exp(-2πλ_B/αδ^2/αB(2/α,1-2/α)r^2)rdr 
       =1/2∫_0^∞exp(-λ_B2π/αδ^2/αB(2/α,1-2/α)t)dt 
       =1/2λ_B2π/αδ^2/αB(2/α,1-2/α)=α/4πλ_B B(2/α,1-2/α) δ^2/α.

By substituting the above expression into (<ref>) and (<ref>), we obtain (<ref>) and (<ref>) respectively. This completes the proof.▪



§ PROOF OF THEOREM <REF>

We conduct the analysis for a typical MU that is located at the origin. We assume that when downloading a file in ℱ_n, the MU will always associate with its nearest SBS, which caches ℱ_n. Note that the OP derived under this assumption is an upper bound for the exact OP. This is because the MU will associate with the second-nearest SBS if it can provide a higher received SINR than that provided by the nearest SBS. Therefore, in some cases, the nearest SBS cannot provide a higher enough SINR (≥δ), while the second-nearest SBS can. According to our assumption, we will neglect these cases, which leads to a higher OP.

Let us denote by z the distance between the typical MU and the nearest SBS that caches ℱ_n. The location of the nearest SBS caching ℱ_n is denoted by x_Z. The fading (power) for an SBS located at x_B, ∀ x_B∈Φ_B, is denoted by h_x_B, which is assumed to be exponentially distributed, i.e., h_x_B∼exp(1). The path-loss function for a given point x_B is ‖ x_B‖^-α.

When random caching is adopted, the distribution of the SBSs that cache ℱ_n can be modeled as an PPP with the intensity of Ω_ℱ_nλ_B. The event that the typical MU can download a file in ℱ_n from an SBS means that the received SINR from the nearest SBS which caches ℱ_n is no less than the threshold δ. Let us denote by ρ(x_Z) the received SINR at the typical MU from the nearest SBS. Then the average probability that the MU can download the file from an SBS is

    (ρ(x_Z)≥δ)
       =∫_0^∞(.h_x_Zz^-α/x_k∈Φ_B\{x_Z}∑h_x_k‖ x_k‖ ^-α≥δ| z)f_Z(z)dz 
       =∫_0^∞(.h_x_Z≥δ(x_k∈Φ_B\{x_Z}∑h_x_k‖ x_k‖ ^-α)/z^-α|z)·
            2πΩ_ℱ_nλ_B z exp(-πΩ_ℱ_nλ_B z^2) dz 
       =∫_0^∞𝔼_I(exp(-z^αδ I))2πΩ_ℱ_nλ_B zexp(-πΩ_ℱ_nλ_B z^2) dz,

where we have I≜x_k∈Φ_B\{x_Z}∑h_x_k‖ x_k‖ ^-α, and the PDF of z, i.e., f_Z(z), is derived by the null probability of a Poisson process with the intensity of Ω_ℱ_nλ_B. Note that the interference I consists of I_1 and I_2, where I_1 is emanating from the SBSs caching the FGs ℱ_q, ∀ q∈𝒩,q≠ n, while I_2 is from the SBSs caching ℱ_n excluding x_Z. The SBSs contributing to I_1, denoted by Φ_n̅, have the intensity (1-Ω_ℱ_n)λ_B, while those contributing to I_2, denoted by Φ_n, have the intensity Ω_ℱ_nλ_B. Correspondingly, the calculation of 𝔼_I(exp(-z^αδ I)) will be split into the product of two expectations over I_1 and I_2. The expectation over I_1 directly follows (<ref>), i.e., we have

    𝔼_I_1(exp(-z^αδ I_1))=exp(-π(1-Ω_ℱ_n)λ_B C(δ,α)z^2),

where C(δ,α) has been defined as 2/αδ^2/αB(2/α,1-2/α). The expectation over I_2 has to take into account z as the distance from the nearest interfering SBS, i.e., we obtain

    𝔼_I_2(exp(-z^αδ I_2))
       =exp(-Ω_ℱ_nλ_B 2π∫_z^∞(1-1/1+ z^αδ r^-α)rdr)
       (a)=exp(-Ω_ℱ_nλ_B πδ^2/αz^22/α∫_δ^-1^∞x^2/α-1/1+x dx) 
       (b)=exp(-Ω_ℱ_nλ_Bπδ z^22/α-2 _2F_1(1,1-2/α;2-2/α;-δ)),

where (a) defines x≜δ^-1z^-αr^α, and _2F_1(·) in (b) is the hypergeometric function. Since we have defined A(δ,α)=2δ/α-2 _2F_1(1,1-2/α;2-2/α;-δ), by substituting (<ref>) and (<ref>) into (<ref>), we have

    (ρ(x_Z)≥δ)=∫_0^∞exp(-π(1-Ω_ℱ_n)λ_B C(δ,α)z^2)
       exp(-πΩ_ℱ_nλ_B z^2A(δ,α))2πΩ_ℱ_nλ_B zexp(-πΩ_ℱ_nλ_B  z^2)dz 
       =Ω_ℱ_n/C(δ, α)(1-Ω_ℱ_n)+A(δ,α)Ω_ℱ_n+Ω_ℱ_n.

It is clear that (𝒬_n)=1-(ρ(z)≥δ). This completes the proof.
▪



§ PROOF OF LEMMA <REF>

Without loss of generality, we assume k=1. Then (<ref>) becomes

    |∏_q=2^K a_q-∏_q=2^Kã_q|≤(K-1) max_q∈{2,⋯,K}|a_q-ã_q|.

Again, without loss of generality, we assume

    |a_2-ã_2|≥⋯≥|a_K-ã_K|.


First, we prove that |a_K-1a_K-ã_K-1ã_K|≤ 2|a_K-1-ã_K-1|, under the condition of |a_K-1-ã_K-1|≥|a_K-ã_K|. To prove this, we discuss the following possible cases.


  §.§.§ When a_K-1≥ã_K-1 and a_K≥ã_K

We have a_K≤ a_K-1-ã_K-1+ã_K. Then

    |a_K-1 a_K-ã_K-1ã_K|
           ≤|a_K-1(a_K-1-ã_K-1+ã_K)-ã_K-1ã_K|
           =|(a_K-1+ã_K)(a_K-1-ã_K-1)|
           ≤ 2|a_K-1-ã_K-1|.



  §.§.§ When a_K-1≥ã_K-1, a_K≤ã_K, and a_K-1a_K≥ã_K-1ã_K

We have

    |a_K-1a_K-ã_K-1ã_K|≤|a_K-1ã_K-ã_K-1ã_K|
               =|a_K-1-ã_K-1|ã_K≤ |a_K-1-ã_K-1|.



  §.§.§ When a_K-1≥ã_K-1, a_K≤ã_K, and a_K-1a_K≤ã_K-1ã_K

We have

    |ã_K-1ã_K-a_K-1a_K|≤|a_K-1ã_K-a_K-1a_K|
               =|a_K-ã_K|a_K-1≤ |a_K-1-ã_K-1|.



  §.§.§ When a_K-1≤ã_K-1, a_K≥ã_K, and a_K-1a_K≥ã_K-1ã_K

We have

    |a_K-1a_K-ã_K-1ã_K|≤|ã_K-1 a_K-ã_K-1ã_K|
               =|a_K-ã_K|ã_K-1≤ |a_K-1-ã_K-1|.



  §.§.§ When a_K-1≤ã_K-1, a_K≥ã_K, and a_K-1a_K≤ã_K-1ã_K

We have

    |ã_K-1ã_K-a_K-1a_K|≤|ã_K-1 a_K-a_K-1a_K|
               =|a_K-1-ã_K-1|a_K≤ |a_K-1-ã_K-1|.



  §.§.§ When a_K-1≤ã_K-1, a_K≤ã_K
 We have a_K≥ã_K+a_K-1-ã_K-1. Then

    |ã_K-1ã_K-a_K-1a_K|
           ≤|ã_K-1ã_K-a_K-1(ã_K+a_K-1-ã_K-1)|
           =|(a_K-1+ã_K)(ã_K-1-a_K-1)|
           ≤ 2|a_K-1-ã_K-1|.

From the above discussions, we can see that |a_K-1a_K-ã_K-1ã_K|≤ 2|a_K-1-ã_K-1|.

Second, as there is |a_K-1a_K-ã_K-1ã_K|≤ 2|a_K-1-ã_K-1|, we have |a_K-1a_K-ã_K-1ã_K|≤ 2|a_K-2-ã_K-2|. With this condition, we can prove that |a_K-2a_K-1a_K-ã_K-2ã_K-1ã_K|≤ 3|a_K-2-ã_K-2| by following the similar steps above. By doing this iteratively, we have

    |∏_q=2^K a_q-∏_q=2^Kã_q|≤(K-1) |a_2-ã_2|.

This completes the proof.▪



IEEEtran



