# An Effective Handover Analysis for the Randomly Distributed Heterogeneous
Cellular Networks

Bin Fang  Wuyang Zhou  Bin Fang -<EMAIL_ADDRESS>Wuyang Zhou-
<EMAIL_ADDRESS>Wireless Information Network Lab, University of Science and
Technology of China, Hefei, China, 230026

###### Abstract

Handover rate is one of the most import metrics to instruct mobility
management and resource management in wireless cellular networks. In the
literature, the mathematical expression of handover rate has been derived for
homogeneous cellular network by both regular hexagon coverage model and
stochastic geometry model, but there has not been any reliable result for
heterogeneous cellular networks (HCNs). Recently, stochastic geometry modeling
has been shown to model well the real deployment of HCNs and has been
extensively used to analyze HCNs. In this paper, we give an effective handover
analysis for HCNs by stochastic geometry modeling, derive the mathematical
expression of handover rate by employing an infinitesimal method for a
generalized multi-tier scenario, discuss the result by deriving some
meaningful corollaries, and validate the analysis by computer simulation with
multiple walking models. By our analysis, we find that in HCNs the handover
rate is related to many factors like the base stations’ densities and
transmitting powers, user’s velocity distribution, bias factor, pass loss
factor and etc. Although our analysis focuses on the scenario of multi-tier
HCNs, the analytical framework can be easily extended for more complex
scenarios, and may shed some light for future study.

## Keywords

Stochastic geometry modeling, handover rates, heterogeneous cellular network.

##  1 Introduction

With the dramatically increasing of wireless traffic as well as the population
of wireless terminals, the traditional homogeneous cellular network cannot
provide sufficient bandwidth for all the wireless terminals. In response to
the capacity challenges, smaller coverage base stations (BSs) are deployed in
hotspots to offload and have a range of tens of meters to several hundred
meters. This brings heterogeneity to traditional cellular network and gives
birth to the heterogeneous cellular network (HCN). Heterogeneity is expected
to be a key feature of the next generation of cellular networks, and an
essential means for providing higher network capacity as well as expanded
indoor and cell-edge coverage. In general, HCNs comprise a conventional
cellular network overlaid with a diverse set of lower-power BSs such as micro
cells, femtocells and perhaps relay BSs.

The BSs in different tiers of HCNs (the tiers of BSs are ordered by transmit
power) may share the same spectra and have different coverage. Handover
happens when a user leaves the coverage of its serving BS and handover rate is
defined as the number of handovers per unit time. The handover in HCNs can be
divided into two types: horizontal handover and vertical handover. Horizontal
handover is the handover between two BSs in the same tier and vertical
handover is the handover between two BSs in different tiers. Compared with
horizontal handovers, vertical handovers are more difficult to implement
because the HCNs may be deployed by different service providers. Thus, to
implement vertical handovers, extra communication overhand between the the
HCNs is essential. Moreover, vertical handovers would lead to additional
transfer delay, jitter and high risk of dropping, which degrade the service
quality. Therefore, handover rate especially vertical handover rate is one of
the most metrics to instruct the deployment of mobility management and
resource management.

The mathematical expression of handover rate in the homogeneous cellular
network has been derived by the regular hexagon coverage model [cite_last01]
and stochastic geometry model [cite_new03][cite_new04]. But for the
heterogeneous cellular network, there has not been any reliable and
generalized handover expression, due to the randomness of the BS positions of
HCNs and the different transmitting power of different tiers. In the real
deployment, the HCN BSs are distributed irregularly, sometimes in an anywhere
plug-and-play manner, which results in a high level of spatial position
randomness. BSs in different tiers have different transmitting power for
communication, leading to different cell size for different tiers. As a
consequence, it is difficult to characterize the cell boundaries and to track
boundary crossings by UEs (i.e. handovers) in the global networks. Few
previous works have resolved the above challenges and give a reliable handover
rate expression.

In the recent years, stochastic geometry modeling has shown its admirable
ability of modeling the position distribution of HCN BSs [cite_06][cite_07],
and has provided tractable accurate performance bounds for cellular wireless
networks. Stochastic geometry is a very powerful mathematical and statistical
tool for the modeling, analysis and design of HCNs with random topologies. For
instance, the modeling is employed for the capacity analysis of random channel
access schemes like ALOHA [cite_04] and carrier sensing multiple access (CSMA)
[cite_05], the capacity analysis of single and multiple tier cellular networks
[cite_last02] and the capacity analysis of cognitive-based networks
[cite_last03].

Hence, in this paper, we investigate the handover rates including horizontal
and vertical handover rates by the stochastic geometry modeling. By employing
an infinitesimal method, the mathematical expression of instantaneous handover
rate is derived for a typical moving UE. From the derivation, we find that the
instantaneous handover rate is related to the instantaneous moving speed, and
is independent of the moving direction. That means only the moving speed
distribution of the memoryless walking model contributes to the handover rates
and the handover rates can be derived through averaging the instantaneous
handover rates by the moving speed distribution. Thus, the derived handover
rate expressions are applicable for all the memoryless walking models. The
derived expressions are validated by computer simulation with multiple walking
models and the impacts of system parameters like BS density, transmit power,
moving velocity of UE, path loss factor are evaluated. Although our analysis
focuses on the scenario of multi-tier HCNs, the analytical framework can be
easily extended for more complex scenarios.

##  2 Downlink System Model

A fairly general model of HCNs considered in this paper contains $N$ tiers of
BSs that are distinguished by their spatial densities, transmit powers, path
loss exponents and biasing factors. For instance, as shown in Figure 1, high-
power macrocell BS networks are overlaid with successively denser and lower
power picocells and femtocells. Macrocell BSs and femtocell BSs can be well
modeled by spatial random processes [cite_06][cite_new02]. Under this model,
the positions of BSs in the _n_ th-tier are modeled according to a homogeneous
PPP (Poisson point process) $\Phi_{n}$ with intensity $\lambda_{n}$ in an
Euclidean plane.

Each BS in the _n_ th-tier has the same transmit power $P_{n}$, and has the
same path loss exponents $\alpha_{n}>2$, ${n=1,...,N}$. Assume that UEs are
uniformly distributed in the Euclidean plane with density of $f_{u}$, and the
movements of UEs are memoryless and are independent of the distributions of
BSs. Memoryless here means the current position of a UE is only related to its
latest position and is independent of its more earlier position, i.e.
$\mathbb{P}[S(t_{0})|S(t_{1}),S(t_{2}),\cdots]=\mathbb{P}[S(t_{0})|S(t_{1})]$,
where $\mathbb{P}[x]$ is the probability of x, $S(t_{i})$ is the position of a
UE at time $t_{i}$, and $t_{i}>t_{i+1}$, i = 0,1,2,$\cdots$.

We assume open access which means a user is allowed to access any tier’s BSs.
And consider a cell association based on maximum biased-received-power (BRP)
(termed biased association), where a mobile UE is associated with the
strongest BS in terms of long-term averaged BRP at the UE. The BRP from the
_j_ th BS in the _n_ th-tier is $P_{r,nj}$ that can be given by

$P_{r,nj}=P_{n}L_{0}(R_{nj}/r_{0})^{-\alpha_{n}}B_{n}$ (1)

where $R_{nj}$ is the distance of the _j_ th BSs in the _n_ th-tier from the
origin, $L_{0}$ is the path loss at the reference distance $r_{0}$ (typically
about $(4\pi/\nu)^{-2}$ for $r_{0}=1$, where $\nu$ denotes the wavelength).
And $B_{n}$ is the bias factor of admission [cite_new05], that could extends
the cell range (or coverage) of the _n_ th-tier by employing $B_{n}>1$. The
considered BRP is a long-term averaged value and fading is averaged out, and
so does not include fading. We assume that handover happens only when the UE
is going across the boundary of the current BS’s coverage, which is determined
by the long-term averaged BRP and is shown as Figure 1.

##  3 Problem Formulation

Consider a typical UE that is at the origin and is admitted to the _k_ th BS
of the _m_ th-tier initially. As the typical UE moves, it may immigrate to
other BSs in the same tier or other tiers. So its admission state can be
depicted as Figure 2, that at time _t_ , the typical UE is admitted to the _k_
th BS with probability $P_{a,k}(t)$, or is admitted to other BS with
probability $P_{a,k}(0)-P_{a,k}(t)$. Thus, the instantaneous transition rate
from k state to the $\overline{k}$ state at time _t_ is

$H_{k}^{m}(t)=-\frac{dP_{a,k}(t)}{dt}$ (2)

Since the movement of the typical UE is a memoryless process and is
independent of BSs’ distributions, then the instantaneous transition rate is
stable and can be given by

$H_{k}^{m}\triangleq-\lim\limits_{t\rightarrow 0}\frac{dP_{a,k}(t)}{dt}$ (3)

$H_{k}^{m}$ is the instantaneous handover rate of the typical UE indeed. Then,
the handover rate in a region with area $S$ can be given by

$\lambda_{h}=\mathbb{E}[\sum\limits_{m=1}^{N}H_{k}^{m}f_{u}S]=\sum\limits_{m=1}^{N}(\mathbb{E}[H_{k}^{m}])f_{u}S$
(4)

where $\mathbb{E}[x]$ is the expectation of variable $x$, and
$\mathbb{E}[H_{k}^{m}]$ is the average handover rate of the typical UE. As the
handover is assumed to happen at the boundary of BS coverage, not at the
boundary of a specified region with area S, the average handover rate is
independent of the shape of the specified region. In the following section, we
would derive the arithmetic expression of $\mathbb{E}[H_{k}^{m}]$.

##  4 Derivation of Handover Rate

In this section, we would derive $\mathbb{E}[H_{k}^{m}]$ in two steps:
firstly, we use an infinitesimal method to derive the instantaneous handover
rate at time 0 of the typical UE with instantaneous velocity $v$, then average
the instantaneous handover rate by the distribution of the velocity $v$. Thus,
the impacts of walking model on the handover rate can be reflected by the
distribution of the velocity $v$, i.e. the walking model decides the
distribution of velocity $v$, then further affects the average handover rate.
So the analysis is applicable for all the memoryless walking models with that
thought. Note that the moving direction of the typical UE would not affect the
instantaneous handover rate and we would give the explanation in the following
derivation.

In the following, we would give the position distribution of the associated BS
firstly, and then derive the handover probability of
$(P_{a,k}(0)-P_{a,k}(t))$, the instantaneous handover rate ($H_{k}^{m}$) and
the average handover rate ($\mathbb{E}[H_{k}^{m}]$) in turn.

###  4.1 Position Distribution of Associated BS

Denote $(R_{nj},\theta_{nj})$ as the polar coordinate of the _j_ th BS in the
_n_ th-tier. Assume that the typical UE is admitted to the _k_ th BS of the
_m_ th-tier initially, thus $P_{r,mk}>P_{r,nj}$ for all
$n\in\\{1,\cdots,N\\}$. According to the max-BRP based association and the BRP
definition in equation (1), the distance boundary condition of these
unassociated BSs can be derived as equation (5) based on $P_{r,mk}>P_{r,nj}$,

$R_{nj}>(\frac{P_{n}B_{n}}{P_{m}B_{m}})^{\frac{1}{\alpha_{n}}}(R_{mk})^{\frac{\alpha_{m}}{\alpha_{n}}}\triangleq
R_{n}^{lb}$ (5)

where $R_{n}^{lb}$ is defined as the distance lower bound of the _n_ th-tier
BSs for the clarity of expression.

According to the distance lower bound of each tier, the probability density
function (PDF) of the associated BS’s distance ($R_{mk}$) can be derived by
using the null probability of a 2-D Poisson process with density $\lambda$ in
an area A, which is $\exp(-\lambda A)$. By setting $\lambda=\lambda_{n}$ and
$A=\pi(R_{n}^{lb})^{2}$ for the _n_ th tier, (n=1,$\cdots$,N)[cite_new05], we
could give the probability density function (PDF) of $R_{mk}$ as

$f(R_{mk})=2\pi\lambda_{m}R_{mk}\exp\\{-\pi\sum\limits_{n=1}^{N}\lambda_{n}(R_{n}^{lb})^{2}\\}$
(6)

where $R_{m}^{lb}=R_{mk}$. Note that
$\int\limits_{0}^{\infty}f(R_{mk})dR_{mk}\triangleq\gamma_{m}$ is the
probability that the typical UE is admitted to a _m_ th-tier BS.

Since BSs are deployed as a PPP (Poisson point process), the $\theta_{mk}$ is
uniformly distributed in the range of $[0,2\pi]$, and its PDF is given by

$f(\theta_{mk})=\frac{1}{2\pi}$ (7)

Since the distributions of $\theta_{mk}$ and $R_{mk}$ are independent of the
coordinate axis and the moving direction of UE, we can assume that the X axis
is along the moving direction of the typical UE at time 0.

###  4.2 Derivation of Handover Probability

When the typical UE moves an infinitesimal distance of $r$, i.e. the typical
UE moves to the point $(r,0)$ with $r\rightarrow 0$, the BRP from the _j_ th
BS in the _n_ th-tier is

$P_{r,nj}^{new}\\!\\!=\frac{P_{n}B_{n}L_{0}}{(\sqrt{(R_{nj}\cos(\theta_{nj})\\!-\\!r)^{2}\\!+\\!(R_{nj}\sin(\theta_{nj}))^{2}})^{\alpha_{n}}}$
(8)

According to the max-BRP based association, after the typical UE moves to the
new point, it is still admitted to the primary BS _k_ only when the BRP from
BS _k_ is larger than the BRP from anyone else. Given the position of BS _k_ ,
the probability that the typical UE keeps the primary link to the BS k is
denoted by $P_{a,k}(r|R_{mk},\theta_{mk})$ and is given by equation (9).

$\begin{array}[]{l}P_{a,k}(r|R_{mk},\theta_{mk})=\prod\limits_{n=1}^{N}\mathbb{P}[P_{r,mk}^{new}\geq
P_{r,nj}^{new}]=\prod\limits_{n=1}^{N}\mathbb{P}[\cos(\theta_{nj})\leq\frac{R_{nj}^{2}+x_{nj}}{2rR_{j}}]\end{array}$
(9)

where $\mathbb{P}[x]$ denotes the probability of event x, and $x_{nj}$ is
defined as

$x_{nj}\\!=\\!r^{2}\\!-\\!(\frac{B_{n}P_{n}}{B_{m}P_{m}})^{\frac{2}{\alpha_{n}}}(R_{mk}^{2}\\!-\\!2rR_{mk}\cos(\theta_{mk})\\!+\\!r^{2})^{\frac{\alpha_{m}}{\alpha_{n}}}$
(10)

And $\cos(\theta_{nj})\leq\frac{R_{nj}^{2}+x_{nj}}{2rR_{j}}$ is derived
according to $P_{r,mk}^{new}\geq P_{r,nj}^{new}$ in equation (9).

Hence, according to equation (9), the typical UE would keep its primary link
if all the BSs in the _n_ th-tier (n=1,$\cdots$,N) are in the region of
$\cos(\theta_{nj})\leq\frac{R_{nj}^{2}+x_{nj}}{2rR_{j}}$, or wound not other
wise.

We call the region of $\cos(\theta_{nj})>\frac{R_{nj}^{2}+x_{nj}}{2rR_{j}}$ as
the bad region of _n_ th-tier BSs, which is shown as Figure 3. It means that
if there are _n_ th-tier BSs in the _n_ th-tier bad region, the typical UE
would immigrate from the serving BS to one of those BSs. Thus,
$P_{a,k}(r|R_{k},\theta_{k})$ equals to the probability that no BS is in its
bad region for all tiers.

Denote the area of the _n_ th-tier bad region as
$A_{mn}(r,R_{mk},\theta_{mk})$, then we can give the null probability of the
the PPP $\Phi_{n}$ in the _n_ th-tier bad region as
$\exp(\lambda_{n}A_{mn}(r,R_{mk},\theta_{mk}))$. Since all the PPPs
$\\{\Phi_{n}\\}_{n=1,\cdots,N}$ are independent, $P_{a,k}(r|R_{k},\theta_{k})$
is the product of those null probabilities and can be given by the following
equation

$P_{a,k}(r|R_{k},\theta_{k})=\prod\limits_{n=1}^{N}\exp(-\lambda_{n}A_{mn}(r,R_{mk},\theta_{mk}))$
(11)

According to $\cos(\cdot)\leq 1$ and the definition of bad region, the
$\theta_{nj}$ boundary conditions of the _n_ th-tier bad region can be given
as

$1\geq\cos(\theta_{nj})>\frac{R_{nj}^{2}+x_{nj}}{2rR_{j}},$ (12)

Thus, based on the $R_{nj}$ boundary condition in equation (5) and
$\theta_{nj}$ boundary conditions in equation (12), we can further derive the
boundary conditions of the _n_ th-tier bad region in the Appendix 7.1, and
give the results as equation (13), where $\vartheta_{nj}$ is defined as
$\vartheta_{nj}=\theta_{nj}^{max}-\theta_{nj}^{min}$, $\theta_{nj}^{max}$ and
$\theta_{nj}^{min}$ are the upper bound and lower bound of $\theta_{nj}$,
respectively. The shapes of bad regions with different boundary conditions in
equation (13) can be depicted as Figure 3.

As Figure 3 (a) shows, when $R_{n}^{lb}<R_{mk}$ holds, according to the
derivation in Appendix 7.1, the range of $R_{nj}$ is
$[R_{n}^{lb},r+\sqrt{r^{2}-x_{nj}}]$ and the range of $\theta_{nj}$ is
$[-\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}),\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}})]$
for a certain $R_{nj}$.

As Figure 3 (b) shows, when both $R_{n}^{lb}>R_{mk}$ and
$|\theta_{mk}|<\arccos(\frac{-R_{mk}}{R_{n}^{lb}})$ hold, the ranges of
$R_{nj}$ and $\theta_{nj}$ are the same as the ranges of $R_{nj}$ and
$\theta_{nj}$ in Figure 3 (a). When both $R_{n}^{lb}>R_{mk}$ and
$|\theta_{mk}|>\arccos(\frac{-R_{mk}}{R_{n}^{lb}})$ hold, as shown by Figure 3
(c), according to the analysis in Appendix 7.1, the range of $R_{nj}$ is
$[R_{n}^{lb},r+\sqrt{r^{2}-x_{nj}}]$, and the range of $\theta_{nj}$ is
$[-\pi,\pi]$ if $R_{nj}\in[R_{n}^{lb},-r+\sqrt{r^{2}-x_{nj}}]$, or is
$[-\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}),\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}})]$
if $R_{nj}\in[-r+\sqrt{r^{2}-x_{nj}},r+\sqrt{r^{2}-x_{nj}}]$.

$\scriptsize\left\\{\begin{array}[]{lc}\\!\\!\\!R_{nj}\in\emptyset,\vartheta_{nj}=0,&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\cos(\theta_{mk})>\frac{R_{mk}}{R_{n}^{lb}}\\\
\\!\\!\\!R_{nj}\in[R_{n}^{lb},r+\sqrt{r^{2}-x_{nj}}],\vartheta_{nj}=\left\\{\begin{array}[]{lc}\\!\\!\\!2\pi,&\\!\\!\\!R_{nj}\in[R_{n}^{lb},-r+\sqrt{r^{2}-x_{nj}})\\\
\\!\\!\\!2\arccos\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}},&\\!\\!\\!R_{nj}\in[-r+\sqrt{r^{2}-x_{nj}},r+\sqrt{r^{2}-x_{nj}}]\end{array}\right.&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\cos(\theta_{mk})<-\frac{R_{mk}}{R_{n}^{lb}}\\\
\\!\\!\\!R_{nj}\in(R_{n}^{lb},r+\sqrt{r^{2}-x_{nj}}),\vartheta_{nj}=2\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}),&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-\frac{R_{mk}}{R_{n}^{lb}}\leq\cos(\theta_{mk})\leq\frac{R_{mk}}{R_{n}^{lb}}\end{array}\right.$
(13)

Based on the boundary conditions, the area of the _n_ th-tier bad region
($A_{mn}(r,R_{mk},\theta_{mk})$) can be derived as the equation (14), where
$I(c)$ is the index function that equals 1 if the condition c holds or 0
otherwise. The first term in equation (14) corresponds to the case of
$\frac{-R_{mk}}{R_{n}^{lb}}<\cos(\theta_{mk})<\frac{R_{mk}}{R_{n}^{lb}}$,
which is depicted by Figure 3 (a) and (b), and the second and third terms
correspond to the case of $\cos(\theta_{mk})<\frac{-R_{mk}}{R_{n}^{lb}}$ that
is depicted by Figure 3 (c).

$\begin{array}[]{l}A_{mn}(r,R_{mk},\theta_{mk})=[\int\limits_{R_{n}^{lb}}^{r+\sqrt{r^{2}-x_{nj}}}2\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}})R_{nj}dR_{nj}]I(-\frac{R_{mk}}{R_{n}^{lb}}\leq\cos(\theta_{mk})\leq\frac{R_{mk}}{R_{n}^{lb}})\\\
+[\int\limits_{-r+\sqrt{r^{2}-x_{nj}}}^{r+\sqrt{r^{2}-x_{nj}}}2\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}})R_{nj}dR_{nj}+\int\limits_{R_{n}^{lb}}^{-r+\sqrt{r^{2}-x_{nj}}}2\pi
R_{nj}dR_{nj}]I(\cos(\theta_{mk})<-\frac{R_{mk}}{R_{n}^{lb}})\end{array}$ (14)

Thus, the handover probability can be obtained as
$(P_{a,k}(0)-\mathbb{E}_{\\{R_{mk},\theta_{mk}\\}}P_{a,k}(r|R_{mk},\theta_{mk}))$.

###  4.3 Handover Rate Derivation

Since the instantaneous handover rate $H_{k}^{m}$ is the derivative of
handover probability, according to equations (3) and (11), the instantaneous
handover rate $H_{k}^{m}$ can be derived as equation (15).

$\begin{array}[]{l}H_{k}^{m}=-\lim\limits_{t\rightarrow
0}\frac{dP_{a,k}(t)}{dt}=-\lim\limits_{t\rightarrow
0}\frac{d\mathbb{E}_{\\{R_{mk},\theta_{mk}\\}}[P_{a,k}(r|R_{mk},\theta_{mk})]}{dt}\\\
\\!\\!=\\!-\\!\mathbb{E}_{\\{R_{mk},\theta_{mk}\\}}[\lim\limits_{t\rightarrow
0}\frac{dP_{a,k}(r|R_{mk},\theta_{mk})}{dt}]\\\
\\!\\!=\\!\\!-\mathbb{E}_{\\{R_{mk},\theta_{mk}\\}}[\lim\limits_{r\rightarrow
0}\\!(\\!\frac{dP_{a,k}(r|R_{mk},\theta_{mk})}{dr}\frac{dr}{dt})]\\\
\\!\\!\mathop{=}\limits^{(a)}\\!\\!-\mathbb{E}_{\\{R_{mk},\theta_{mk}\\}}[\lim\limits_{r\rightarrow
0}\frac{d}{dr}(-\sum\limits_{n=1}^{N}\lambda_{n}A_{mn}(r,R_{mk},\theta_{mk})\cdot\exp(-\sum\limits_{n=1}^{N}\lambda_{n}A_{mn}(r,R_{mk},\theta_{mk})))]\cdot
v\\\
\\!\\!\mathop{=}\limits^{(b)}\\!\mathbb{E}_{\\{R_{mk},\theta_{mk}\\}}[\sum\limits_{n=1}^{N}(\lambda_{n}\cdot\lim\limits_{r\rightarrow
0}\frac{dA_{mn}(r,R_{mk},\theta_{mk})}{dr})]\cdot v\end{array}$ (15)

where $v=\lim\limits_{r\rightarrow 0}\frac{dr}{dt}$ is the instantaneous
velocity of the typical UE at time $t=0$, (a) is obtained according to the
expression of $P_{a,k}(r|R_{mk},\theta_{mk})$ in equation (11), and (b) is
obtained due to $\lim\limits_{r\rightarrow 0}A_{mn}(r,R_{mk},\theta_{mk})=0$
and $\lim\limits_{r\rightarrow
0}\exp(-\sum\limits_{n=1}^{N}\lambda_{n}A_{mn}(r,R_{mk},\theta_{mk}))=1$.

Denote $H_{k}^{m-n}$ as the instantaneous handover rate from the _k_ th BS in
the _m_ th-tier to the BSs in the _n_ th-tier. Thus

$H_{k}^{m}=\sum\limits_{n=1}^{N}H_{k}^{m-n}$ (16)

According to the derivation of $H_{k}^{m}$ and the expression of $H_{k}^{m}$
in equation (15), we can give $H_{k}^{m-n}$ as

$H_{k}^{m-n}=\mathbb{E}_{\\{R_{mk},\theta_{mk}\\}}[\lambda_{n}\cdot\lim\limits_{r\rightarrow
0}\frac{dA_{mn}(r,R_{mk},\theta_{mk})}{dr}]\cdot v$ (17)

Based on the expression of $A_{mn}(r,R_{mk},\theta_{mk})$ in equation (14),
$\lim\limits_{r\rightarrow 0}\frac{dA_{mn}(r,R_{mk},\theta_{mk})}{dr}$ is
derived in Appendix 7.2, and the result is given as equation (18).

$\begin{array}[]{l}\lim\limits_{r\rightarrow
0}\frac{dA_{mn}(r,R_{mk},\theta_{mk})}{dr}=I(\cos(\theta_{mk})<-\frac{R_{mk}}{R_{n}^{lb}})[-2\pi\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})]\\\
+I(-\frac{R_{mk}}{R_{n}^{lb}}<\cos(\theta_{mk})<\frac{R_{mk}}{R_{n}^{lb}})[-2\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})\arccos(\frac{R_{n}^{lb}}{R_{mk}}\cos(\theta_{mk}))+2\sqrt{(R_{n}^{lb})^{2}-(\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk}))^{2}}]\end{array}$
(18)

######  Theorem 1.

The instantaneous handover rate from _m_ th-tier BSs to _n_ th-tier BSs for a
UE with instantaneous velocity $v$ is $H_{k}^{m-n}$ and is given by equation
(19).

$H_{k}^{m-n}=8\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}[\sqrt{\frac{1-z^{2}}{(\frac{R_{mk}}{R_{n}^{lb}})^{2}-z^{2}}}+\sqrt{\frac{(\frac{R_{mk}}{R_{n}^{lb}})^{2}-z^{2}}{1-z^{2}}}]dz(R_{n}^{lb})^{2}\exp\\{-\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}$
(19)

###### Proof.

Based on equations (17) and (18), $H_{k}^{m-n}$ is further derived in Appendix
7.3. ∎

Hence, the average handover arrival rate $\lambda_{h}$ in equation (4) can be
given by equation (20), where $f_{v}(v)$ is the probability density function
of the velocity $v$ in the specified region with area $S$ and is determined by
the walking model of UEs.

$\begin{array}[]{l}\lambda_{h}=\sum\limits_{m=1}^{N}(\mathbb{E}[H_{k}^{m}f_{u}S])=\sum\limits_{m=1}^{N}\sum\limits_{n=1}^{N}\mathbb{E}[H_{k}^{m-n}f_{u}S]\\\
=\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv\\!\sum\limits_{m=1}^{N}\\!\sum\limits_{n=1}^{N}\\!\\{\\!8\lambda_{n}\lambda_{m}\\!\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz(R_{n}^{lb})^{2}\\!\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\\}\end{array}$
(20)

Similarly, the average handover arrival rate from _m_ th-tier BSs to _n_ th-
tier BSs (denoted as $\lambda_{h}^{m-n}$) can be derived as equation (21).

$\begin{array}[]{l}\lambda_{h}^{m-n}=\mathbb{E}[H_{k}^{m-n}f_{u}S]\\\
=\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv8\lambda_{n}\lambda_{m}\\!\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz(R_{n}^{lb})^{2}\\!\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\end{array}$
(21)

Although the derived handover rate expressions are not closed-form, these
expressions can be efficiently computed numerically as opposed to the usual
Monte Carlo methods that rely on repeated random sampling to compute these
results.

###  4.4 Discussions of Handover Rates

In equations of (20) and (21), the general expressions of handover rates have
been derived. In this section, some corollaries and special cases of handover
rates in the stochastic modeling of the HCNs would be given.

######  Corollary 1.

$\lambda_{h}^{m-n}=\lambda_{h}^{n-m}$ holds for any $m,n\in\\{1,\cdots,N\\}$,
that is, the forward and reverse handover rates between any two tiers are the
same.

###### Proof.

See Appendix 7.4. ∎

Corollary 1 holds under the condition that the UE movements in different tiers
are homogeneous and UEs are uniformly distributed. Corollary 1 indicates that
the mobility between any two tiers would reach statical balance over time.

######  Corollary 2.

When $N=1$, all the BSs are homogeneous and $\\{\lambda_{n}\\}=\lambda$, the
expression of average handover rate $\lambda_{h}$ can be further simplified as

$\lambda_{h}(N=1)=\frac{4\sqrt{\lambda}}{\pi}f_{u}S\int\limits_{0}^{+\infty}vf_{v}(v)dv$
(22)

###### Proof.

When $N=1$, according to equation of (20), $\lambda_{h}$ can be further
derived as equation (23), where (a) is obtained due to $R_{m}^{lb}=R_{mk}$.

$\begin{array}[]{l}\lambda_{h}(N=1)=\sum\limits_{m=1}^{N}(\mathbb{E}[H_{k}^{m}f_{u}S])\mathop{=}\limits^{(a)}\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv8\lambda^{2}\\!\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{1}\\!\\!2dzR_{mk}^{2}\\!\exp\\{\\!-\\!\pi\lambda
R_{mk}^{2})\\}dR_{mk}\\\
=\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv16\lambda^{2}\int\limits_{0}^{+\infty}R_{mk}^{2}\exp\\{\\!-\\!\pi\lambda
R_{mk}^{2})\\}dR_{mk}\\\
=\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv\frac{8\lambda}{\pi\sqrt{\lambda}}\mathcal{Q}(0)\\\
=\frac{4\sqrt{\lambda}}{\pi}f_{u}S\int\limits_{0}^{+\infty}vf_{v}(v)dv\end{array}$
(23)

∎

Corollary 2 is consistent with the handover rate expression of homogeneous
cellular network given in [cite_last01][cite_new03].

######  Corollary 3.

When $\\{\alpha_{n}\\}=\alpha$, i.e. BSs in all the tiers have the same path
loss exponent, the handover rate of $\lambda_{h}^{m-n}$ can be simplified as
equation (24),

$\begin{array}[]{l}\lambda_{h}^{m-n}(\\{\alpha_{n}\\}=\alpha)\\!\\!=\\!\\!\frac{2\lambda_{n}\lambda_{m}\beta_{n}^{2}f_{u}S}{\pi(\sum\limits_{i=1}^{N}\lambda_{i}\beta_{i}^{2})^{1.5}}\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv\cdot\\!\\!\\!\\!\\!\\!\int\limits_{0}^{\min(1,\beta_{n})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{\beta_{n}^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{\beta_{n}^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz\end{array}$
(24)

where $\beta_{n}\triangleq(\frac{P_{n}B_{n}}{P_{m}B_{m}})^{\frac{1}{\alpha}}$.

###### Proof.

When $\\{\alpha_{n}\\}=\alpha$, according to equation (5), $R_{n}^{lb}$ can be
simplified as
$R_{n}^{lb}=(\frac{P_{n}B_{n}}{P_{m}B_{m}})^{\frac{1}{\alpha}}R_{mk}=\beta_{n}R_{mk}$.
Then $\lambda_{h}^{m-n}$ in equation (21) can be further derived as equation
(25).

$\begin{array}[]{l}\lambda_{h}^{m-n}(\\{\alpha_{n}\\}=\alpha)=\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv8\lambda_{n}\lambda_{m}\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{\min(1,\beta_{n})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{\beta_{n}^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{\beta_{n}^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz\beta_{n}^{2}R_{mk}^{2}\\!\exp\\{\\!-\\!\pi(\sum\limits_{i=1}^{N}\lambda_{i}\beta_{i}^{2})R_{mk}^{2}\\}dR_{mk}\\\
=\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv8\lambda_{n}\lambda_{m}\beta_{n}^{2}\int\limits_{0}^{\min(1,\beta_{n})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{\beta_{n}^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{\beta_{n}^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz\\!\int\limits_{0}^{+\infty}R_{mk}^{2}\\!\exp\\{\\!-\\!\pi(\sum\limits_{i=1}^{N}\lambda_{i}\beta_{i}^{2})R_{mk}^{2}\\}dR_{mk}\\\
=\\!f_{u}\\!S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv8\lambda_{n}\lambda_{m}\beta_{n}^{2}\int\limits_{0}^{\min(1,\beta_{n})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{\beta_{n}^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{\beta_{n}^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz\frac{1}{4\pi(\sum\limits_{i=1}^{N}\lambda_{i}\beta_{i}^{2})^{1.5}}\\\
=\frac{2\lambda_{n}\lambda_{m}\beta_{n}^{2}f_{u}S}{\pi(\sum\limits_{i=1}^{N}\lambda_{i}\beta_{i}^{2})^{1.5}}\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv\int\limits_{0}^{\min(1,\beta_{n})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{\beta_{n}^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{\beta_{n}^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz\end{array}$
(25)

∎

######  Corollary 4.

For a UE with constant velocity $v$, its residence time in a BS coverage of
_m_ th-tier (denoted by $T_{r}^{m}$) is exponential distributed with average
value $\frac{\gamma_{m}}{H_{k}^{m}}$, i.e. the PDF of $T_{r}^{m}$ can be given
by

$f(T_{r}^{m})=\frac{H_{k}^{m}}{\gamma_{m}}\exp(-\frac{H_{k}^{m}}{\gamma_{m}}T_{r}^{m})$
(26)

where $\gamma_{m}$ is the probability that a UE is associated with the _m_ th-
tier BS and is given by equation (27) referred to [cite_new05].

$\gamma_{m}=2\pi\lambda_{m}\int\limits_{0}^{\infty}R_{mk}\exp\\{-\pi\sum\limits_{n=1}^{N}\lambda_{n}(R_{n}^{lb})^{2}\\}dR_{mk}$
(27)

###### Proof.

When the velocity $v$ is constant, the instantaneous handover rate $H_{k}^{m}$
would keep constant and do not change with time. So the transition rate from k
state to $\overline{k}$ in Figure 2 is constant, i.e.
$\frac{1}{\gamma_{m}}\lim\limits_{t\rightarrow
t_{0}}\frac{dP_{a,k}(t-t_{0})}{dt}=-\frac{H_{k}^{m}}{\gamma_{m}}$ is constant
for all $t_{0}\geq 0$. Thus

$P_{a,k}(t)=\gamma_{m}\exp(-\frac{H_{k}^{m}}{\gamma_{m}}t)$ (28)

due to $P_{a,k}(0)=\gamma_{m}$. So,

$\begin{array}[]{l}f(T_{r}^{m})=-\frac{1}{\gamma_{m}}\frac{P_{a,k}(t)}{dt}=\frac{H_{k}^{m}}{\gamma_{m}}\exp(-\frac{H_{k}^{m}}{\gamma_{m}}T_{r}^{m})\end{array}$
(29)

∎

######  Corollary 5.

When the UEs are not uniformly distributed in the whole region and UEs in each
tier BS coverage are uniformly distributed, the handover rate from _m_ th-tier
BS to _n_ th-tier BS in a specified region with area S can be given by
equation (30).

$\begin{array}[]{l}\lambda_{h}^{m-n}(\\{f_{m,u}\\}\neq
f_{u})=\mathbb{E}[H_{k}^{m-n}f_{m,u}S]\\\
=f_{m,u}S\\!\int\limits_{0}^{+\infty}\\!vf_{v}(v)dv8\lambda_{n}\lambda_{m}\\!\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz(R_{n}^{lb})^{2}\\!\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\end{array}$
(30)

where $f_{m,u}$ is the UE density in the _m_ th tier BS coverage. And the
total handover rate in the specified region is $\lambda_{h}(\\{f_{m,u}\\}\neq
f_{u})=\sum\limits_{m=1}^{N}\sum\limits_{n=1}^{N}\lambda_{h}^{m-n}(\\{f_{m,u}\\}\neq
f_{u})$.

###### Proof.

For a UE admitted to a _m_ th-tier BS, its average handover rate to a _n_ th-
tier BS is $\frac{H_{k}^{m}}{\gamma_{m}}$. On the other hand, the average
coverage area of _m_ th-tier BSs in the specified region is $S\gamma_{m}$. So
the total handover rate from _m_ th-tier BSs to _n_ th-tier BSs in the
specified region is
$\lambda_{h}^{m-n}=\mathbb{E}[\frac{H_{k}^{m-n}}{\gamma_{m}}f_{m,u}S\gamma_{m}]=\mathbb{E}[H_{k}^{m-n}f_{m,u}S]$.
∎

######  Corollary 6.

When the walking models of UEs in different tier BS coverage are different,
the handover rate from the _m_ th-tier BS to _n_ th-tier BS in a specified
region with area S can be given by equation (31).

$\begin{array}[]{l}\lambda_{h}^{m-n}(\\{f_{m,v}(v)\\}\neq
f_{v}(v))=\mathbb{E}[H_{k}^{m-n}f_{u}S]\\\
=f_{u}S\\!\int\limits_{0}^{+\infty}\\!vf_{m,v}(v)dv8\lambda_{n}\lambda_{m}\\!\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}\\!\\![\sqrt{\frac{1\\!-\\!z^{2}}{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}}\\!\\!+\\!\\!\sqrt{\frac{(\frac{R_{mk}}{R_{n}^{lb}})^{2}\\!-\\!z^{2}}{1\\!-\\!z^{2}}}]\\!dz(R_{n}^{lb})^{2}\\!\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\end{array}$
(31)

where $f_{m,v}(v)$ is the velocity distribution in the _m_ th-tier BS
coverage. And the total handover rate in the specified region is
$\lambda_{h}(\\{f_{m,v}(v)\\}\neq
f_{v}(v))=\sum\limits_{m=1}^{N}\sum\limits_{n=1}^{N}\lambda_{h}^{m-n}(\\{f_{m,v}(v)\\}\neq
f_{v}(v))$.

###### Proof.

According to the derivation of $H_{k}^{m-n}$, $H_{k}^{m-n}$ is only related to
the velocity distribution in the _m_ th-tier BS coverage and is independent of
the velocity distribution in other tier BS coverage. Thus, Corollary 6 is can
be directly derived from the derivation of $H_{k}^{m-n}$ in the former
section. ∎

From the above analysis, we could give the general steps for handover rate
analysis by stochastic geometry modeling: firstly, obtain the PDF of the
associated BS’s position, secondly, by the infinitesimal method, get the area
of the bad region, based on which, then derive the instantaneous handover
probability, and derive the instantaneous handover rate through taking the
derivative, and at last, average the instantaneous handover rate by the
distribution of UE’s velocity.

##  5 Simulation Results

###  5.1 Validation of the Analysis

Now that we have developed the general expression of handover rate for HCNs by
stochastic geometry modeling, it is important to see how well the analytical
results match the computer simulation. Here, we consider two types of walking
models, i.e. straight-line walking model and RWP walking model. For the
straight-line walking model, UE would move without changing its moving
direction, while for the RWP model, UE would change its moving direction for a
randomly chosen direction of $[0,2\pi]$ and then keep the direction for a
randomly chosen duration of [0,100s]. For both models, the velocity is
uniformly distributed in $[0,2v]$.

In the simulation, we considered a circle region with radius of 10Km,
generated BSs of each tier and their positions by the PPP $\Phi_{n}$ with the
density of $\lambda_{n}$ for each try, and generated the UEs with the density
of $f_{u}$. For each try, all the UEs move with a specified walking model for
104 seconds. We counted the number of handover in a specified region with area
S = 1km2 and gave the handover rates by averaging the results.

A total of two tiers are modeled according to PPP in the simulation (N=2). For
the simulation, some parameters keep constant and these parameters are $f_{u}$
= 100/km2, S = 1km2, $P_{1}=1$, $P_{2}=0.2$, $B_{1}=1$, $\alpha_{1}=3.5$ and
$\lambda_{1}=1$/km2. Other parameters would vary in different tries and these
parameters are the velocity $v$, the 2nd-tier BS density $\lambda_{2}$ and,
2nd-tier BS’s path loss factor $\alpha_{2}$ and the 2nd-tier BS’s bias factor
$B_{2}$.

Figure 4 compares the analytical and experimental handover rates of
$\lambda_{h}^{1-1}$, $\lambda_{h}^{1-2}$ and $\lambda_{h}^{2-2}$. In Figure 4,
the dotted lines represent the analytical results, the circles and triangles
represent the simulation results obtained by straight-line walking model and
RWP walking model, respectively, And we set $\lambda_{2}=2$/km2, $B_{2}=1$ and
$\alpha_{2}=3.5$. In Figure 4, it is observed that the handover rates of
$\lambda_{h}^{1-1}$, $\lambda_{h}^{1-2}$ and $\lambda_{h}^{2-2}$ increase
linearly with the average velocity of UE. This matches the derived expression
of handover rate. In Figure 4, we can see that $\lambda_{h}^{1-1}$ is the
largest, $\lambda_{h}^{1-2}$ takes the second place and $\lambda_{h}^{2-2}$ is
the minimum for any value of average velocity. This can be explained as
follows. Since $P_{1}=1$ and $P_{2}=0.2$, the coverage of 1st-tier BS is much
larger than the coverage of 2nd-tier BS. So the 1st-tier BS has longer
boundary line than the 2nd-tier BS. Thus, the total length of the boundary
line between two 1st-tier BSs is the longest, the total length of the boundary
line between a 1st-tier BS and a 2nd-tier BS is shorter and the total length
of boundary line between two 2nd-tier BSs is the shortest. Hence, the handover
is most likely to happen at the boundary between two 1st-tier BSs, is less
likely to happen at the boundary between a 1st-tier BS and a 2nd-tier BS, and
is least likely to happen at the boundary between two 2nd-tier BSs. From the
figure, it can be seen that the relative error between the analytical results
and the simulation results are less than 3% for both walking models. The
relative error is most likely brought by the limited simulation time. The good
matching validates that the analysis is reliable with the variation of UE
velocity. In the simulation, we find that the straight-line walking model and
the RWP walking model without pause time almost have the same handover rates
when they have the same velocity distribution.

Figure 5 illustrates the analytical and experimental results of the total
handover rate ($\lambda_{h}$) and the handover rate between the two tiers
($\lambda_{h}^{1-2}$) with different 2nd-tier BSs density ($\lambda_{2}$) for
different average UE velocity ($v$). From the figure, we can observe that
handover rates increase linearly with the average UE velocity, and both
$\lambda_{h}$ and $\lambda_{h}^{1-2}$ increase with the 2nd-tier BS density.
This is because larger BS density leads to smaller cells and UEs are more
likely to move out the smaller cells. The analytical results match well with
experimental results of straight-line walking model and RWP walking model,
with relative error less than 3%. This validates that the analysis is reliable
with the variation of both BS density and moving velocity.

Figure 6 demonstrates the analytical and experimental results of handover
rates with variations of 2nd-tier BS density ($\lambda_{2}$) and path loss
factor ($\alpha_{2}$). The analytical results match well with the experimental
results for both walking models with relative error less than 3%. This
validates that the analysis is reliable with the variations of both BS density
and path loss factors. From this figure, we can observe that all the handover
rates increase with the 2nd-tier BS density. The handover rate of
$\lambda_{h}^{1-1}$ increases with 2nd-tier BS path loss factor $\alpha_{2}$,
$\lambda_{h}^{1-2}$ increases with $\alpha_{2}$ firstly, then decreases with
$\alpha_{2}$, and $\lambda_{h}^{2-2}$ decreases with $\alpha_{2}$. This is
because the coverage of the 2nd-tier BS decreases with $\alpha_{2}$ and the
coverage of the 1st-tier BS increases with $\alpha_{2}$ according to equation
(1). Thus, the boundary line of the 1st-tier BSs increases with the
$\alpha_{2}$ and the boundary line of the 2nd-tier BSs decreases with the
$\alpha_{2}$. So, the total length of the boundary line between two 1st-tier
BSs increases with $\alpha_{2}$, the total length of the boundary line between
two 2nd-tier BSs decreases with $\alpha_{2}$, and the total length of the
boundary line between a 1st-tier BS and a 2nd-tier BS increases firstly with
$\alpha_{2}$, then decreases with $\alpha_{2}$. It is the reason that leads to
the increasing of $\lambda_{h}^{1-1}$, the decreasing of $\lambda_{h}^{2-2}$,
and the changing of $\lambda_{h}^{1-2}$.

Figure 7 demonstrates the experimental results of the forward and reverse
handover rates between the two tiers for the RWP walking model with the
variations of 2nd-tier BS density, path loss factor and average velocity.
Figure 7 shows that the forward handover rates ($\lambda_{h}^{1-2}$) match
well with reverse handover rates ($\lambda_{h}^{2-1}$) between the two tiers
with relative error less than 3%. This validates the reliability of Corollary
1.

Figure 8 demonstrates the analytical and experimental results of residual time
distributions with constant velocity. The analysis results (denoted by ‘Theo.’
in the figure) is calculated by equation (26), and the experimental results
(denoted by ‘Simu.’ in the figure) is obtained by the RWP model with constant
velocity $v$=5m/s. From the figure, we can see that the experimental results
match well with the analysis results, which validates the effectiveness of
Corollary 4. As the figure illustrated, the residual time in the 1st-tier BS
is much larger than the residual time in the 2nd-tier BS due to the fact that
the 1st-tier BS coverage is much larger than the 2nd-tier BS coverage. As the
figure shown, both the residual time in the 1st-tier BS and in the 2nd-tier BS
decrease with the 2nd-tier BS density $\lambda_{2}$, but the decrement of
residual time in the 1st-tier BS is more significant. This can be explained as
follows. With the increase of 2nd-tier BS density, the average coverage of
each tier BS decreases, so the residual time decreases for each tier. On the
other hand, with the increase of 2nd-tier BS density, the average 1st-tier BS
coverage decreases more significantly because the coverage of each 1st-tier BS
is much larger than the 2nd-tier BS coverage and the increased 2nd-tier BSs
would be more likely to occupy the 1st-tier BSs coverage. Thus, the decrement
of residual time in the 1st-tier BS is more significant.

###  5.2 Effect of Bias Factor

Figure 9 shows the numerical results of handover rates with the variation of
2nd-tier BS’s bias factor ($B_{2}$). From this figure, we can observe that the
total handover rate ($\lambda_{h}$) decreases with the 2nd-tier BS’s bias
factor, $\lambda_{h}^{1-1}$ and $\lambda_{h}^{1-2}$ decrease with $B_{2}$ and
$\lambda_{h}^{2-2}$ increases with $B_{2}$ when $B_{2}\in[1,2]$, which is a
reasonable range of $B_{2}$. This can be explained as follows. According to
equation (20), when $R_{mk}$ equals $R_{n}^{lb}$, $\lambda_{h}^{m-n}$ could
reach the minimal value, and $\lambda_{h}^{m-n}$ decreases with $R_{mk}$ if
$R_{mk}<R_{n}^{lb}$ holds. Thus, when $B_{2}$ increases and
$R_{1k}<R_{2}^{lb}$ holds, $R_{1k}$ increases and $\lambda_{h}^{1-2}$
decreases. When $B_{2}$ increases, the coverage of 1st-tier BS decreases and
the coverage of 2nd-tier BS increases according to equation (1), so the
boundary line of 1st-tier BS decreases and the boundary line of the 2nd-tier
BS increases, which lead to the results that $\lambda_{h}^{1-1}$ decreases
with $B_{2}$ and $\lambda_{h}^{2-2}$ increases with $B_{2}$. The results
demonstrate that we can decrease the total handover rate by reasonably
adjusting the bias factors.

##  6 Conclusion

In the literature, there has not been any general handover rate derivation for
the heterogeneous cellular networks. Thus, in this paper, we give a
generalized handover analytical framework by employing the stochastic geometry
modeling for heterogeneous cellular networks, derive the arithmetic expression
of handover rate, give some meaningful corollories and validate the analysis
by computer simulation. The analysis may shed some light for future extension
and study.

## Acknowledgements

The authors wish to thank the anonymous reviewers for their constructive
comments. This paper is supported by National programs for High Technology
Research and Development (2012AA01A50603).

##  7 Appendix

###  7.1 Boundary conditions of bad region

According to the bad region boundary conditions in equation (12), we can
further derive $\vartheta_{nj}$ as follow:

$\left\\{\begin{array}[]{lc}\vartheta_{nj}=0,&\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}>1\\\
\vartheta_{nj}=2\arccos\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}},&-1\leq\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}\leq
1\\\
\vartheta_{nj}=2\pi,&\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}<-1\end{array}\right.$
(32)

For the case of $-1\leq\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}\leq 1$, we can
further derive the inequations of $\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}\leq 1$
and $\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}\geq-1$ as equations (33) and (34),
respectively,

$\begin{array}[]{l}R_{nj}^{2}-2rR_{nj}+x_{nj}\leq 0\\\ \Rightarrow
R_{nj}\\!\in\\![r-\sqrt{r^{2}-x_{nj}},r+\sqrt{r^{2}-x_{nj}}]\end{array}$ (33)

$\begin{array}[]{l}R_{nj}^{2}+2rR_{nj}+x_{nj}\geq 0\\\ \Rightarrow
R_{nj}\\!\in\\!(\\!-\\!\infty,\\!-\\!r\\!-\\!\sqrt{r^{2}\\!-\\!x_{nj}})\\!\bigcup(-\\!r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}},\\!+\\!\infty)\end{array}$
(34)

Hence, we can give the range of $R_{nj}$ as equation (35) according to
equation (5) and $R_{nj}>0$.

$\min(R_{n}^{lb},|\\!-\\!r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}}|)\\!\leq\\!R_{nj}\\!\leq\\!r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}}$
(35)

As $-r+\sqrt{r^{2}\\!-\\!x_{nj}}>0$ when $r\rightarrow 0^{+}$, we can further
simplify equation (35) as

$\min(R_{n}^{lb},\\!-\\!r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}})\\!\leq\\!R_{nj}\\!\leq\\!r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}}$
(36)

For case of $\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}>1$, similar to equation (33),
we can further derive the inequation of $\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}>1$
as $R_{nj}\in(\max(R_{nj}^{lb},r+\sqrt{r^{2}-x_{nj}}),+\infty)$.

For the case of $\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}<-1$, similar to equation
(34), we can further derive the inequation of
$\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}}<-1$ as
$R_{nj}\in(R_{nj}^{lb},-r+\sqrt{r^{2}-x_{nj}})$

On the other hand, according to the expression of $R_{n}^{lb}$ in equation (5)
and $x_{nj}$ in equation (10), $R_{n}^{lb}=-r+\sqrt{r^{2}\\!-\\!x_{nj}}$ holds
when $r=0$. And the derivative of $-r+\sqrt{r^{2}\\!-\\!x_{nj}}$ at $r=0$ is
derived as

$\begin{array}[]{l}\lim\limits_{r\rightarrow
0}\frac{d(-r+\sqrt{r^{2}\\!-\\!x_{nj}})}{dr}=(\frac{B_{n}P_{n}}{B_{m}P_{m}})^{\frac{1}{\alpha_{n}}}(R_{mk})^{\frac{\alpha_{m}}{\alpha_{n}}}\frac{-\cos(\theta_{mk})}{R_{mk}}-1=-\frac{R_{n}^{lb}}{R_{mk}}\cos(\theta_{mk})-1\end{array}$
(37)

Thus, $-r+\sqrt{r^{2}\\!-\\!x_{nj}}$ increases at $r=0$ if
$\cos(\theta_{mk})<-\frac{R_{mk}}{R_{n}^{lb}}$ or decreases otherwise. So,

$\left\\{\begin{array}[]{lc}\\!\\!R_{n}^{lb}\\!>\\!-r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}},&\\!\\!\cos(\theta_{mk})\\!>\\!-\frac{R_{mk}}{R_{n}^{lb}}\\\
\\!\\!R_{n}^{lb}\\!<\\!-r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}},&\\!\\!\cos(\theta_{mk})\\!<\\!-\frac{R_{mk}}{R_{n}^{lb}}\end{array}\right.$
(38)

Similarly, the following relationships can be derived,

$\left\\{\begin{array}[]{lc}\\!\\!R_{n}^{lb}\\!>\\!r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}},&\\!\\!\cos(\theta_{mk})\\!>\\!\frac{R_{mk}}{R_{n}^{lb}}\\\
\\!\\!R_{n}^{lb}\\!<\\!r\\!+\\!\sqrt{r^{2}\\!-\\!x_{nj}},&\\!\\!\frac{-R_{mk}}{R_{n}^{lb}}\\!\\!<\cos(\theta_{mk})\\!<\\!\frac{R_{mk}}{R_{n}^{lb}}\end{array}\right.$
(39)

Based on the above relationships, the boundary conditions of equation (32) can
be derived as equation (13).

###  7.2 Derivation of $\lim\limits_{r\rightarrow
0}\frac{dA_{mn}(r,R_{mk},\theta_{mk})}{dr}$

For the clarity of expression, we define the expressions of
$A_{mn}^{(1)}(r,R_{mk},\theta_{mk})$, $A_{mn}^{(2)}(r,R_{mk},\theta_{mk})$ and
$A_{mn}^{(3)}(r,R_{mk},\theta_{mk})$ in the equation of (40).

$\begin{array}[]{c}A_{mn}^{(1)}(r,R_{mk},\theta_{mk})=[\int\limits_{R_{n}^{lb}}^{r+\sqrt{r^{2}-x_{nj}}}2\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}})R_{nj}dR_{nj}]I(-\frac{R_{mk}}{R_{n}^{lb}}\leq\cos(\theta_{mk})\leq\frac{R_{mk}}{R_{n}^{lb}})\\\
A_{mn}^{(2)}(r,R_{mk},\theta_{mk})=[\int\limits_{-r+\sqrt{r^{2}-x_{nj}}}^{r+\sqrt{r^{2}-x_{nj}}}2\arccos(\frac{R_{nj}^{2}+x_{nj}}{2rR_{nj}})R_{nj}dR_{nj}]I(\cos(\theta_{mk})<-\frac{R_{mk}}{R_{n}^{lb}})\\\
A_{mn}^{(3)}(r,R_{mk},\theta_{mk})=[\int\limits_{R_{n}^{lb}}^{-r+\sqrt{r^{2}-x_{nj}}}2\pi
R_{nj}dR_{nj}]I(\cos(\theta_{mk})<-\frac{R_{mk}}{R_{n}^{lb}})\end{array}$ (40)

Based on those expressions and equation (14), the following relationship can
be obtained,

$A_{mn}=A_{mn}^{(1)}+A_{mn}^{(2)}+A_{mn}^{(3)}$ (41)

Then $\lim\limits_{r\rightarrow 0}\frac{dA_{mn}(r,R_{mk},\theta_{mk})}{dr}$
can be derived by deriving $\lim\limits_{r\rightarrow
0}\frac{dA_{mn}^{(i)}(r,R_{mk},\theta_{mk})}{dr}$, i = 1,2,3. For simplicity,
we define $\varphi^{u}(r)$, $\varphi^{d}(r)$ and $f(r,R_{nj})$ as
$\varphi^{u}(r)=r+\sqrt{r^{2}-x_{nj}}$,
$\varphi^{d}(r)=-r+\sqrt{r^{2}-x_{nj}}$, and
$f(r,R_{nj})=2R_{nj}\arccos(\frac{x_{nj}+R_{nj}^{2}}{2rR_{nj}})$,
respectively.

$\begin{array}[]{l}\lim\limits_{r\rightarrow
0}\frac{dA_{mn}^{(1)}}{dr}=\lim\limits_{r\rightarrow
0}\frac{d}{dr}\int\limits_{R_{n}^{lb}}^{\varphi^{u}(r)}f(r,R_{nj})dR_{nj}I^{(1)}\\\
=\\!\\!\lim\limits_{r\rightarrow
0}[f(r,\varphi^{u}(r))\frac{d}{dr}\varphi^{u}(r)\\!\\!+\\!\\!\\!\\!\int\limits_{R_{n}^{lb}}^{\varphi^{u}(r)}\\!\frac{d}{dr}f(r,R_{nj})dR_{nj}]I^{(1)}\\\
\mathop{=}\limits^{(a)}\lim\limits_{r\rightarrow
0}[\int\limits_{R_{n}^{lb}}^{\varphi^{u}(r)}\\!\frac{d}{dr}f(r,R_{nj})dR_{nj}]I^{(1)}\\\
=\lim\limits_{r\rightarrow
0}[\int\limits_{R_{n}^{lb}}^{\varphi^{u}(r)}\frac{-2R_{nj}}{\sqrt{1-(\frac{R_{nj}^{2}+x_{nj}}{2R_{nj}r})^{2}}}\frac{\frac{dx_{nj}}{dr}r-(R_{nj}^{2}+x_{nj})}{2R_{nj}r^{2}}dR_{nj}]I^{(1)}\\\
\mathop{=}\limits^{y=R_{nj}^{2}}\lim\limits_{r\rightarrow
0}I^{(1)}[\int\limits_{(R_{n}^{lb})^{2}}^{(\varphi^{u}(r))^{2}}\frac{-1}{\sqrt{1-\frac{(y+x_{nj})^{2}}{4yr^{2}}}}\frac{\frac{dx_{nj}}{dr}r-(y+x_{nj})}{r^{2}}\frac{1}{2\sqrt{y}}dy]\\\
=\lim\limits_{r\rightarrow
0}I^{(1)}[\int\limits_{(R_{n}^{lb})^{2}}^{(\varphi^{u}(r))^{2}}\frac{-1}{\sqrt{4r^{4}-4r^{2}x_{nj}-(y+x_{nj}-2r^{2})^{2}}}\cdot\frac{(\frac{dx_{nj}}{dr}r-2r^{2})-(y+x_{nj}-2r^{2})}{r}dy]\\\
=\lim\limits_{r\rightarrow
0}-I^{(1)}[-\frac{\frac{dx_{nj}}{dr}r-2r^{2}}{r}\arccos\frac{y+x_{nj}-2r^{2}}{\sqrt{4r^{4}-4r^{2}x_{nj}}}+\frac{1}{r}\sqrt{4r^{4}-4r^{2}x_{nj}-(y+x_{nj}-2r^{2})^{2}}]|_{(R_{n}^{lb})^{2}}^{(\varphi^{u}(r))^{2}}\\\
=\lim\limits_{r\rightarrow
0}-I^{(1)}[\frac{\frac{dx_{nj}}{dr}r-2r^{2}}{r}\arccos(\frac{(R_{n}^{lb})^{2}+x_{nj}-2r^{2}}{\sqrt{4r^{4}-4r^{2}x_{nj}}})-\frac{1}{r}\sqrt{4r^{4}-4r^{2}x_{nj}-((R_{n}^{lb})^{2}+x_{nj}-2r^{2})^{2}}]\\\
\mathop{=}\limits^{(b)}I^{(1)}[-\frac{2(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})\arccos(\frac{R_{n}^{lb}}{R_{mk}}\cos(\theta_{mk}))+2\sqrt{(R_{n}^{lb})^{2}-(\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk}))^{2}}]\end{array}$
(42)

where
$I^{(1)}=I(-\frac{R_{mk}}{R_{n}^{lb}}\leq\cos(\theta_{mk})\leq\frac{R_{mk}}{R_{n}^{lb}})$,
and $(a)$ is given by $f(r,\varphi^{u}(r))=0$, $(b)$ is given due to
$\lim\limits_{r\rightarrow 0}x_{nj}=-(R_{n}^{lb})^{2}$,
$\lim\limits_{r\rightarrow
0}\frac{d}{dr}x_{nj}=\frac{2(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})$ and
the L’Hopital rule.

$\begin{array}[]{l}\lim\limits_{r\rightarrow
0}\frac{dA_{mn}^{(2)}}{dr}=\lim\limits_{r\rightarrow
0}\int\limits_{\varphi^{d}(r)}^{\varphi^{u}(r)}f(r,R_{nj})dR_{nj}I^{(2)}\\\
=\lim\limits_{r\rightarrow
0}[f(r,\varphi^{u}(r))\frac{d}{dr}\varphi^{u}(r)-f(r,\varphi^{d}(r))\frac{d}{dr}\varphi^{d}(r)+\int\limits_{\varphi^{d}(r)}^{\varphi^{u}(r)}\frac{d}{dr}f(r,R_{nj})dR_{nj}]I^{(2)}\\\
\mathop{=}\limits^{(c)}[2\pi
R_{n}^{lb}(\frac{R_{n}^{lb}}{R_{mk}}\cos(\theta_{mk})\\!+1\\!)\\!-\\!2\pi\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})]I^{(2)}\end{array}$
(43)

where $I^{(2)}=I(\cos(\theta_{mk})<-\frac{R_{mk}}{R_{n}^{lb}})$, and $(c)$ is
obtained by $f(r,\varphi^{u}(r))=0$, $f(r,\varphi^{d}(r))=2\pi$,
$\lim\limits_{r\rightarrow
0}\frac{d}{dr}\varphi^{d}(r)=-\frac{R_{n}^{lb}}{R_{mk}}\cos(\theta_{mk})-1$
and $\lim\limits_{r\rightarrow
0}\int\limits_{\varphi^{d}(r)}^{\varphi^{u}(r)}\frac{d}{dr}f(r,R_{nj})dR_{nj}=-\\!2\pi\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})$,
which can be derived similarly as equation (42).

$\begin{array}[]{l}\lim\limits_{r\rightarrow
0}\frac{dA_{mn}^{(3)}}{dr}=\lim\limits_{r\rightarrow
0}\frac{d}{dr}\int\limits_{R_{n}^{lb}}^{\varphi^{d}(r)}2\pi
R_{nj}dR_{nj}I^{(2)}\\\ =\lim\limits_{r\rightarrow
0}[2\pi\varphi^{d}(r)\frac{d\varphi^{d}(r)}{dr}]I^{(2)}\\\ =-2\pi
R_{n}^{lb}(\frac{R_{n}^{lb}}{R_{mk}}\cos(\theta_{mk})+1)I^{(2)}\end{array}$
(44)

which is obtained due to $\lim\limits_{r\rightarrow
0}\varphi^{d}(r)=R_{n}^{lb}$.

Thus, add equations of (42) (43) and (44) together, we can give
$\lim\limits_{r\rightarrow 0}\frac{dA_{mn}(r,R_{mk},\theta_{mk})}{dr}$ as
equation (18).

###  7.3 Derivation of $H_{k}^{m-n}$

According equation (17), $H_{k}^{m-n}$ can be further expressed as equation
(45),

$\begin{array}[]{l}H_{k}^{m-n}\\!\\!=\\!\\!\lambda_{n}v\\!\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{2\pi}\lim\limits_{r\rightarrow
0}\\!\frac{dA_{mn}}{dr}\\!f(\theta_{mk})\\!f(R_{mk})d\theta_{mk}dR_{mk}\\\
=\\!\\!\lambda_{n}v\\!\\!\int\limits_{0}^{+\infty}\\!\int\limits_{0}^{2\pi}\lim\limits_{r\rightarrow
0}\\!\frac{dA_{mn}}{dr}\\!\lambda_{m}\\!R_{mk}\\!\exp(\\!-\\!\pi\\!\lambda_{m}\\!R_{mk}^{2}\\!)d\theta_{mk}dR_{mk}\end{array}$
(45)

For simplicity, we define $h_{k,1}^{m-n}$, $h_{k,2}^{m-n}$ and $h_{k,3}^{m-n}$
in equation of (46).

$\begin{array}[]{l}h_{k,1}^{m-n}=I^{(2)}[-2\pi\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})]\\\
h_{k,2}^{m-n}\\!\\!=\\!\\!I^{(1)}[\\!-\\!2\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})\\!\arccos(\frac{R_{n}^{lb}}{R_{mk}}\cos(\theta_{mk}))]\\\
h_{k,3}^{m-n}=I^{(1)}[2\sqrt{(R_{n}^{lb})^{2}-(\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk}))^{2}}]\end{array}$
(46)

Thus

$H_{k}^{m-n}=\mathbb{E}[h_{k,1}^{m-n}+h_{k,2}^{m-n}+h_{k,3}^{m-n}]$ (47)

And $\mathbb{E}[h_{k,1}^{m-n}]$, $\mathbb{E}[h_{k,2}^{m-n}]$ and
$\mathbb{E}[h_{k,3}^{m-n}]$ can be derived as equations of (48), (49) and
(50), respectively.

$\begin{array}[]{l}\mathbb{E}[h_{k,1}^{m-n}]=2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{\pi-\arccos(\min(1,\frac{R_{mk}}{R_{n}^{lb}}))}^{\pi}-2\pi\frac{(R_{n}^{lb})^{2}}{R_{mk}}\cos(\theta_{mk})R_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}d\theta_{mk}dR_{mk}\\\
=2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}[2\pi\frac{(R_{n}^{lb})^{2}}{R_{mk}}(-\sin(\theta_{mk}))|_{\pi-\arccos(\min(1,\frac{R_{mk}}{R_{n}^{lb}}))}^{\pi}]R_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\\\
=2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}[2\pi\frac{(R_{n}^{lb})^{2}}{R_{mk}}\sqrt{1-(\min(1,\frac{R_{mk}}{R_{n}^{lb}}))^{2}}]R_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\end{array}$
(48)

$\small\begin{array}[]{l}\mathbb{E}[h_{k,2}^{m-n}]\mathop{=}\limits^{z=\cos(\theta_{mk})}2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}^{-\min(1,\frac{R_{mk}}{R_{n}^{lb}})}[-2\frac{(R_{n}^{lb})^{2}}{R_{mk}}z\arccos(\frac{R_{n}^{lb}}{R_{mk}}z)]\frac{-1}{\sqrt{1-z^{2}}}dzR_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\\\
=2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}2\frac{(R_{n}^{lb})^{2}}{R_{mk}}[(-\sqrt{1-z^{2}}\arccos(\frac{R_{n}^{lb}}{R_{mk}}z))|_{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}^{-\min(1,\frac{R_{mk}}{R_{n}^{lb}})}+2\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}\frac{\sqrt{1-z^{2}}}{\sqrt{(\frac{R_{mk}}{R_{n}^{lb}})^{2}-z^{2}}}]R_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\\\
=2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}2\frac{(R_{n}^{lb})^{2}}{R_{mk}}[-\pi\sqrt{1-(\min(1,\frac{R_{mk}}{R_{n}^{lb}}))^{2}}+2\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}\frac{\sqrt{1-z^{2}}}{\sqrt{(\frac{R_{mk}}{R_{n}^{lb}})^{2}-z^{2}}}]R_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\end{array}$
(49)

$\begin{array}[]{l}\mathbb{E}[h_{k,3}^{m-n}]\mathop{=}\limits^{z=\cos(\theta_{mk})}2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}^{-\min(1,\frac{R_{mk}}{R_{n}^{lb}})}[2\sqrt{1-(\frac{R_{n}^{lb}}{R_{mk}}z)^{2}}\frac{-1}{\sqrt{1-z^{2}}}R_{n}^{lb}]R_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\\\
=2\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}2\frac{(R_{n}^{lb})^{2}}{R_{mk}}[2\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}\frac{\sqrt{(\frac{R_{mk}}{R_{n}^{lb}})^{2}-z^{2}}}{\sqrt{1-z^{2}}}]R_{mk}\exp\\{\\!-\\!\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\end{array}$
(50)

Thus, $H_{k}^{m-n}$ can be derived as equation of (19) by adding the results
of equations of (48), (49) and (50) together.

###  7.4 Proof of Corollary 1

According to the equation (21), $\lambda_{h}^{m-n}=\lambda_{h}^{n-m}$ can be
proved by proving that the equation of $H_{k}^{m-n}=H_{k}^{n-m}$ holds. Based
on the expression of $H_{k}^{m-n}$ in equation (19), the proof can be given by
the equation (51), where (d) follows from plugging
$R_{mk}=x^{\frac{\alpha_{n}}{\alpha_{m}}}(\frac{P_{m}B_{m}}{P_{n}B_{n}})^{\frac{1}{\alpha_{m}}}$,
(e) follows from plugging $z_{1}=\frac{R_{m,n}^{lb}}{x}z$, and
$R_{i,n}^{lb}=(\frac{P_{i}B_{i}}{P_{n}B_{n}})^{\frac{1}{\alpha_{i}}}x^{\frac{\alpha_{n}}{\alpha_{i}}}$,
(f) follows from plugging $x=R_{nj}$, where $R_{nj}$ is the nearest distance
of the _n_ th-tier BSs to the origin.

$\begin{array}[]{l}H_{k}^{m-n}=8\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{0}^{\min(1,\frac{R_{mk}}{R_{n}^{lb}})}[\sqrt{\frac{1-z^{2}}{(\frac{R_{mk}}{R_{n}^{lb}})^{2}-z^{2}}}+\sqrt{\frac{(\frac{R_{mk}}{R_{n}^{lb}})^{2}-z^{2}}{1-z^{2}}}]dz(R_{n}^{lb})^{2}\exp\\{-\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i}^{lb})^{2})\\}dR_{mk}\\\
\mathop{=}\limits^{(d)}8\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{0}^{\min(1,\frac{R_{m,n}^{lb}}{x})}[\sqrt{\frac{1-z^{2}}{(\frac{R_{m,n}^{lb}}{x})^{2}-z^{2}}}+\sqrt{\frac{(\frac{R_{m,n}^{lb}}{x})^{2}-z^{2}}{1-z^{2}}}]dzx^{\frac{\alpha_{n}}{\alpha_{m}}+1}(\frac{P_{m}B_{m}}{P_{n}B_{n}})^{\frac{1}{\alpha_{m}}}\exp\\{-\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i,n}^{lb})^{2})\\}dx\\\
\mathop{=}\limits^{(e)}8\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{0}^{\min(1,\frac{x}{R_{m,n}^{lb}})}[\sqrt{\frac{(\frac{x}{R_{m,n}^{lb}})^{2}-z_{1}^{2}}{1-z_{1}^{2}}}+\sqrt{\frac{1-z_{1}^{2}}{(\frac{x}{R_{m,n}^{lb}})^{2}-z_{1}^{2}}}]dz_{1}(R_{m,n}^{lb})^{2}\exp\\{-\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i,n}^{lb})^{2})\\}dx\\\
\mathop{=}\limits^{z=z_{1}}8\lambda_{n}\lambda_{m}v\int\limits_{0}^{+\infty}\int\limits_{0}^{\min(1,\frac{x}{R_{m,n}^{lb}})}[\sqrt{\frac{1-z^{2}}{(\frac{x}{R_{m,n}^{lb}})^{2}-z^{2}}}+\sqrt{\frac{(\frac{x}{R_{m,n}^{lb}})^{2}-z^{2}}{1-z^{2}}}]dz(R_{m,n}^{lb})^{2}\exp\\{-\pi\sum\limits_{i=1}^{N}(\lambda_{i}(R_{i,n}^{lb})^{2})\\}dx\\\
\mathop{=}\limits^{(f)}H_{k}^{n-m}\end{array}$ (51)

## Figures

### Figure 1 - Example of donwlink HCNs with three tiers of BSs: high-power
macrocell BSs (red square) are overlaid with successively denser and lower
power picocells (red triangle) and femtocells (red circle).

### Figure 2 - The admission state and its transition of the typical UE.

### Figure 3 - The bad region when the typical UE moves to (r,0).

### Figure 4 - The average handover rates between tiers in 1$km^{2}$ region
with different $\overline{v}$.

### Figure 5 - The total handover rate ($\lambda_{h}$) and handover rate
between different tiers ($\lambda_{h}^{1-2}$) in 1$km^{2}$ region with
different average velocity $\overline{v}$ and 2nd-tier BSs density
$\lambda_{2}$.

### Figure 6 - The handover rates between tiers in 1$km^{2}$ region with
different average velocity $\overline{v}$ and 2nd-tier path loss factor
$\alpha_{2}$.

### Figure 7 - The forward and reverse handover rates between the two tiers in
1$km^{2}$ region with different average velocity $\overline{v}$, 2nd-tier BSs
density $\lambda_{2}$ and 2nd-tier path loss factor $\alpha_{2}$.

### Figure 8 -The CDF of residual time in the 2-tier BSs with different
$\lambda_{2}$.

### Figure 9 -The average handover arrival rates in 1$km^{2}$ region with
different 2nd-tier bias factor ($B_{2}$).


# TEQ: Trainable Equivalent Transformation for Quantization of LLMs

Wenhua Cheng  Yiyang Cai  Kaokao Lv  Haihao Shen
Intel
{wenhua.cheng, yiyang.cai, kaokao.lv<EMAIL_ADDRESS>
###### Abstract

As large language models (LLMs) become more prevalent, there is a growing need
for new and improved quantization methods that can meet the computationalast
layer demands of these modern architectures while maintaining the accuracy. In
this paper, we present TEQ, a trainable equivalent transformation that
preserves the FP32 precision of the model output while taking advantage of
low-precision quantization, especially 3 and 4 bits weight-only quantization.
The training process is lightweight, requiring only 1K steps and less than
$1\text{\textperthousand}$ of the original model’s trainable parameters.
Furthermore, the transformation does not add any computational overhead during
inference. Our results are on-par with the state-of-the-art (SOTA) methods on
typical LLMs. Our approach can be combined with other methods to achieve even
better performance. The code is available at https://github.com/intel/neural-
compressor.

##  1 Introduction

Large language models (LLMs) have not only shown breakthrough performance in a
wide range of benchmarks and tasks but played an increasingly important role
in daily life, e.g., ChatGPT (OpenAI, ) in information retrieval and Copilot
(Github, ) in programming. However, as LLMs’ model size keeps growing
dramatically, their significant memory footprint and heavy computation
requirements become a major bottleneck of their usage.

One of the most promising ways to alleviate this challenge is quantization,
which can reduce storage and computational overhead. Quantization converts
high-bit floating-point data to lower-bit representations, and it has become
an effective model compression technique.

Quantization methods can generally be divided into two categories:
quantization aware training (QAT) (Shen et al., 2021; Zhuang et al., 2021;
Gong et al., 2019; Esser et al., 2019; Louizos et al., 2018) and post-training
quantization (PTQ) (Frantar et al., 2022; Li et al., 2022; Xiao et al., 2022;
Wei et al., 2022; Frantar and Alistarh, 2022; Hubara et al., 2021; Nagel et
al., 2020; Hassibi et al., 1993; LeCun et al., 1989). Their effectiveness has
been validated for a wide range of models. However, several issues still need
to be addressed, especially for LLMs. QAT simulates the quantization behavior
in the training/finetuning phase, but such a process is very costly for LLMs
due to their unprecedented parameter scale. In contrast, PTQ requires no
training and thus has drawn rising attention. However, PTQ is prone to large
accuracy drops, especially for extreme low-bit quantization. This provides
LLMs’ PTQ methods with great opportunities for improvement.

Lower-bit quantization (e.g., Int4, W4) has recently been widely discussed
since memory bandwidth is becoming the main bottleneck of LLMs. However, most
existing works focus on computer vision models (He et al., 2016; Howard et
al., 2017) that are much smaller than current popular LLMs such as
BLOOM-176B(Scao et al., 2022), OPT-175B(Zhang et al., 2022). Other extreme
quantization methods (Bai et al., 2020; Zhang et al., 2020) rely on the
knowledge distillation technique, introducing extra overhead. GPTQ(Frantar et
al., 2022) tunes the weights based on optimal brain surgeon(Hassibi et al.,
1993) and successfully achieves low-bit quantization on LLMs with low
computation overhead.

Our proposed method reduces the compression error by introducing a trainable
equivalent transformation (Fig. 1), which keeps the mathematical equivalency
of model output at FP32 precision. Moreover, the training cost is
significantly low, only 1k steps of batch size 1 with around less than one-
thousandth trainable parameters of the original models. Also, our method is
orthogonal to current popular LLMs quantization methods, and better accuracy
results could be achieved by combining ours with them.

In summary, the contribution of this paper is threefold:

  * •

We introduce a trainable equivalent transformation for the quantization of
LLMs, which keeps the model output unchanged at FP32 precision. Besides, the
training is quite lightweight.

  * •

Experimental results show our method could achieve results on par with or
better than the SOTA methods.

  * •

We also show that our method could be combined to get the new SOTA
performance.

In the following, we first briefly introduce the work related to ours in
Section 2. We then present the trainable equivalent transformation in Section
3. Experiments and conclusion are described in Sections 4 and 5 respectively.

##  2 Related Work

#### Quantization-aware Training.

QAT methods are widely used in model compression. By enabling finetuning
process, quantized models’ accuracy can often be on par with or even better
than those of original models. (Louizos et al., 2018) introduce a
differentiable quantization procedure by converting original weights and
activations’ distribution to categorical distributions. OQAT (Shen et al.,
2021) proposes a combined training scheme of architecture and quantization to
acquire many quantized models. Afterward, they are converted to lower-bit
models and optimized. (Zhuang et al., 2021) propose a progressive quantization
scheme by quantizing activations after weights. Indeed, QAT methods are
popular in relatively small-scale models, but their application in LLMs is
limited due to the expensive training or even fine-tuning costs as mentioned
in Section 1.

#### Post-training Quantization.

A large number of post-training methods quantize weights step by step and
modify unquantized weights to compensate for errors produced by previously
quantized weights. Optimal Brain Damage (OBD) (LeCun et al., 1989) uses
second-derivative information (Hessian-based estimation) to predict the effect
of weights’ perturbation analytically. Optimal Brain Surgeon (OBS) (Hassibi et
al., 1993) applies such an idea by devising a second-order framework for
weight pruning. Afterward, Optimal Brain Quantization (OBQ) migrate OBS’s
pruning framework to quantization since pruning and quantization share the
common idea of introducing perturbation in original models. Finally, GPTQ
(Frantar et al., 2022) improves the original framework’s efficiency by fixing
the quantization order within the layer and calculating the Hessian matrix’s
Cholesky decomposition before quantization. Other PTQ methods use a better
rounding scheme than commonly used rounding-to-nearest (RTN). AdaRound (Nagel
et al., 2020) learns a rounding scheme using mean squared error (MSE) for
layer-wise activation. AQuant (Li et al., 2022) adds a learnable border
function for activation quantization.

#### Large Language Model Quantization.

Researchers are devoting efforts to compression methods particularly designed
for LLMs as more open-source releases are available. LLM.int8() (Dettmers et
al., 2022) discovers peak values in activation outliers’ particular channels.
It proposes methods to ensure that these channels are kept in higher
precision. SmoothQuant (Xiao et al., 2022) addresses the issues mentioned
above by migrating difficulties from activation to weights with a handcrafted
equivalent transformation. ZeroQuant (Yao et al., 2022) devises an end-to-end
quantization and inference pipeline with a novel layer-wise knowledge
distillation algorithm. However, the largest model it has quantized has only
1.3B parameters. GPTQ (Frantar et al., 2022) tunes the weights based on
optimal brain surgeon (Hassibi et al., 1993) and successfully achieves low-bit
quantization on LLMs with low computation overhead. More recent, AWQ (Lin et
al., 2023) propose to search the optimal scales to protect parts of weights,
since they can significantly reduce the error caused by quantization.

##  3 Methodology

Figure 1 presents a schematic illustration of equivalent transformation. In
the following, we introduce the quantization process first. Consider a feed-
forward neural network comprised of $L$ layers, which perform matmul or
convolution operations. We only consider the matmul layer for simplicity since
our method could be easily extended to convolution layers. The $l^{th}$ matmul
operation can be denoted by $y_{l}=w_{l}\cdot x_{l}$. In which $w_{l}$ and
$x_{l}$ are the weights and activation(input), and $y_{l}$ is the
corresponding output. To quantize a tensor, a quantization op presented below
could be applied.

$Q(v)=clip(\left[\frac{v}{s}\right],-n,n),n\in\mathbb{N}$ (1)

where $s$ denotes the quantization scale parameter and $\left[\cdot\right]$
denotes the round-to-nearest (RTN) operation, while $-n$ and $n$ denote the
integer thresholds for clipping. We ignore the zero point for simplicity. For
a normal int8 quantization, i.e., W8A8, we need to quantize activation and
weights both. And for weight-only quantization, only the weights need to be
quantized. Finally, a de-quantization operation will be appended to
reconstruct the float output, normally not equal to $y_{l}$. In summary, the
$L_{l}$’s output after normal quantization is converted to:

$\hat{y_{l}}=Q^{-1}(Q(w_{l})\cdot Q(x_{l}))$ (2)

where $\hat{y_{l}}$ denotes the $L_{l}$’s reconstructed output after
quantization. The value of $(y_{l}-\hat{y_{l}})^{2}$ is usually named as
quantization loss.

Figure 1: A schematic illustration of TEQ, where $s_{w1}$ and $s_{w2}$ are
trainable parameters. A per-channel scale is multiplied at activations while
an inverse scale is multiplied at weights, which could keep the output
equivalent.

###  3.1 Trainable Equivalent Transformation

PTQ tends to cause a noticeable accuracy drop as mentioned before. SmoothQuant
(Xiao et al., 2022) and AWQ (Lin et al., 2023) rely on handcrafted rules to
migrating quantization difficulties of weights and activations. However, these
rules often fall into sub-optimal solutions, which cannot minimize error
caused by quantization. To alleviate this issue, we introduce a trainable
equivalent transformation that enforces the Fp32 output as the same but
greatly improves the quantization robustness. To be more specific, suppose the
shape of $w_{l}$ is ${c^{in}_{l}\times c^{out}_{l}}$, which stands for their
respective input and output channel numbers. For each layer $L_{l}$, we can
multiply a per-channel scaling vector $s_{l}\in\mathbb{R}^{c^{in}_{l}}$ for
weights and append a corresponding inverse scale vector for activation.
Mathematically, this can be restated as

$y_{l}=w_{l}\cdot diag(s_{l})\cdot diag(s_{l})^{-1}\cdot x_{l}$ (3)

operator $diag(\cdot)$ denotes converting a column/row vector to a diagonal
matrix whose eigenvalues are identical to the original vector’s elements.

$diag\left(\begin{bmatrix}s_{1}\\\ s_{2}\\\ \vdots\\\
s_{n}\end{bmatrix}\right)=\begin{bmatrix}s_{1}&&&\\\ &s_{2}&&\\\ &&\ddots&\\\
&&&s_{n}\end{bmatrix}$ (4)

Our observation shows the optimal $s_{w}$ is useful to reduce the quantization
loss. Therefore, we quantize the transformed model rather than the original
one.

n_bits | Methods | OPT-6.7B | OPT-13B | BLOOM-3B | BLOOM-7B1 | LLAMA-7B | LLAMA-13B
---|---|---|---|---|---|---|---
32 | FP32 | 64.97 | 65.54 | 55.65 | 60.29 | 68.87 | 71.06
4 | RTN | 62.99 | 64.17 | 53.17 | 57.80 | 67.41 | 68.86
GPTQ | 63.09 | 64.83 | 54.65 | 58.26 | 64.70 | 70.00
Ours | 63.30 | 64.91 | 53.83 | 58.93 | 67.71 | 69.55
Ours+GPTQ | 63.94 | 65.03 | 54.42 | 59.62 | 65.27 | 69.73
4_g128 | RTN | 64.04 | 64.88 | 54.91 | 59.32 | 67.87 | 70.88
GPTQ | 64.76 | 65.37 | 55.68 | 59.59 | 66.33 | 70.92
Ours | 64.11 | 64.87 | 54.98 | 59.35 | 68.10 | 71.00
Ours+GPTQ | 64.77 | 65.20 | 55.49 | 59.60 | 66.56 | 70.96
Table 1: The w4 average accuracy($\uparrow$) of four tasks, e.g., HellaSwag, WinoGrande, PIQA, and LAMBADA, in LM-eval. g denotes group size. "Ours+GPTQ" means we apply TEQ first and then apply GPTQ afterward. For LLAMA-7B, the result of GPTQ is w/o act-order. Results of act-order are shown in Appendix A.2. n_bits | Methods | OPT-6.7B | OPT-13B | BLOOM-3B | BLOOM-7B1 | LLAMA-7B | LLAMA-13B
---|---|---|---|---|---|---|---
32 | FP32 | 10.86 | 10.12 | 13.48 | 11.36 | 5.68 | 5.09
4 | RTN | 12.10 | 11.32 | 14.75 | 12.09 | 6.29 | 5.53
GPTQ | 11.59 | 10.33 | 14.10 | 11.73 | 6.59 | 5.33
Ours | 11.68 | 10.59 | 14.72 | 12.21 | 6.30 | 5.50
Ours+GPTQ | 11.29 | 10.36 | 14.03 | 11.74 | 6.76 | 5.35
4_g128 | RTN | 11.16 | 10.32 | 13.85 | 11.60 | 5.97 | 5.26
GPTQ | 10.98 | 10.20 | 13.69 | 11.48 | 6.29 | 5.21
Ours | 11.11 | 10.28 | 13.82 | 11.58 | 5.97 | 5.26
Ours+GPTQ | 11.02 | 10.21 | 13.69 | 11.48 | 6.28 | 5.21
Table 2: The w4 perplexity($\downarrow$) on WikiText-2. For LLAMA-7B, the
result of GPTQ is w/o act-order. Results of act-order are shown in Appendix
A.2.

The transformation has two per-channel scale operations, which will introduce
computation overhead. We fuse the weight scale to the weight itself. For the
activation scale, following (Xiao et al., 2022), we fuse it to the previous
layers, such as layernorm(Ba et al., 2016), batchnorm(Ioffe and Szegedy, 2015)
and etc. In all our experiments, we only apply the transformation to the layer
whose scales could be fused, which introduces no extra overhead at deployment.

###  3.2 Training Details

We train the scales $s_{l}$ because there is little knowledge of the best
equivalent transformation due to various models and quantization
configurations. It’s worth mentioning that the count of trainable scales is
much less than the model’s parameters, and the model weights are frozen.

To train the transformation scales, we follow the basic QAT to simulate the
quantization behavior, which could be denoted as

$y_{l_{q}}=(Q^{-1}Q(w_{l}))(Q^{-1}Q(x_{l}))$ (5)

For weight-only quantization, activation quantization will be ignored. We
adopt straight-through estimator (STE) (Bengio et al., 2013) to backward the
gradients.

We use Adam(Kingma and Ba, 2014) optimizer, betas [0.9, 0.9], and weight decay
0. The learning rate is 1e-3 unless explicitly stated and the decay type is
linear. We only train 1000 steps. We use the same loss function as the
original one in the training phase. For example, CrossEntorpy loss is adopted
for LLMs. The $s_{l}$ is usually initialized with 1. However, sometimes
$1.0/sqrt(w_{cin})$ leads to better results, so we pick the better one in our
experiments.

##  4 Experiments

In this section, we evaluate our proposed TEQ’s in different aspects.
Initially, we briefly introduce LLM architectures and tasks included in our
evaluation. Secondly, we illustrate a detailed comparison of our method and
other state-of-the-art (SOTA) methods, and both quantization accuracy and time
are considered.

###  4.1 Experimental Settings

#### Large Language Models.

We conduct our experiments on the most popular LLM architectures, including
LLaMAs (Touvron et al., 2023), BLOOMs (Scao et al., 2022), and OPTs (Zhang et
al., 2022). Parameter scalings ranging from million to billion are all
included.

n_bits | Methods | OPT-6.7B | OPT-13B | BLOOM-3B | BLOOM-7B1 | LLAMA-7B | LLAMA-13B
---|---|---|---|---|---|---|---
32 | FP32 | 64.97 | 65.54 | 55.65 | 60.29 | 68.87 | 71.06
3_g128 | RTN | 56.03 | 49.59 | 52.54 | 57.53 | 64.92 | 67.68
GPTQ | 62.98 | 64.68 | 53.41 | 58.12 | 58.29 | 68.73
Ours | 61.41 | 63.27 | 52.69 | 57.79 | 65.25 | 68.32
Ours+GPTQ | 63.16 | 64.60 | 53.71 | 58.00 | 59.27 | 69.15
Table 3: The 3 bits with group size 128 average accuracy($\uparrow$) of four tasks,e.g., HellaSwag, WinoGrande, PIQA, and LAMBADA, in LM-eval. g denotes group size. For LLAMA-7B, the result of GPTQ is w/o act-order. Results of act-order are shown in Appendix A.2. n_bits | Methods | OPT-6.7B | OPT-13B | BLOOM-3B | BLOOM-7B1 | LLAMA-7B | LLAMA-13B
---|---|---|---|---|---|---|---
32 | FP32 | 10.86 | 10.12 | 13.48 | 11.36 | 5.68 | 5.09
3_g128 | RTN | 22.37 | 40.50 | 15.68 | 12.47 | 7.01 | 5.88
GPTQ | 11.42 | 10.51 | 14.67 | 11.99 | 8.28 | 5.64
Ours | 12.03 | 11.83 | 15.48 | 12.40 | 6.89 | 5.81
Ours+GPTQ | 11.40 | 10.52 | 14.64 | 11.98 | 7.71 | 5.64
Table 4: WikiText-2 perplexity($\downarrow$) of 3 bits with group size 128.
For LLAMA-7B, the result of GPTQ is w/o act-order. Results of act-order are
shown in Appendix A.2.

#### Evaluation and Datasets.

We make assessments on several language tasks to satisfy the task-agnostic
setting. Specifically, we report average accuracy result on four common sense
reasoning tasks by leveraging lm-eval-harness(Gao et al., 2021), including
HellaSwag (Zellers et al., 2019), WinoGrande (Sakaguchi et al., 2021), PIQA
(Bisk et al., 2020) and LAMBADA (Paperno et al., 2016). Furthermore, we
complement our evaluation with perplexity (PPL) analysis on WikiText2 (Merity
et al., 2016), PTB (Marcus et al., 1994) as well as C4 (Raffel et al., 2020).

#### Implementation Details.

Following GPTQ (Frantar et al., 2022), we focus on weight-only quantization
and exclude the last layer When quantifying. We used a single HW accelerator
to quantize models with a scale of around ten billion parameters. We use the
same calibration dataset
pile-10k111https://huggingface.co/datasets/NeelNanda/pile-10k for a fair
comparison.

#### Baseline.

Our primary baseline is vanilla round-to-nearest quantization (RTN) which has
a remarkable result at 4bits using a small group size of 128. We also compare
with a state-of-the-art method GPTQ (Frantar et al., 2022).

###  4.2 Results

As mentioned above, we compare our results with RTN and the SOTA GTPQ(Frantar
et al., 2022). Also, since our method is orthogonal to GPTQ, we report
Ours+GPTQ as well, which applies TEQ first and then runs GPTQ official
code222https://github.com/IST-DASLab/gptq afterward. We mainly focus on the
models around 10B which is commonly used.

#### W4 Quantization.

We first evaluate TEQ on popular 4 bits quantization. Table 1 shows the lm-
eval results of different LLM model architectures and parameter sizes. TEQ
outperforms RTN in all cases except one. Comparing with GPTQ, TEQ shows better
results in 6 out of 12 scenarios. After combining GPTQ, new state-of-the-art
results could be achieved in 5 scenarios. In summary, TEQ could be helpful in
8 out of 12 scenarios. Table 8 shows the hyper-parameters that we used in the
experiements.

We also evaluate WikiText2 ppl in table 2 w/o group size and group size 128.
TEQ is better or on par with RTN. Similarly, the combined approach (Ours and
GPTQ) shows comparable or better results than standalone GPTQ.

#### W3 Quantization.

We also evaluate TEQ at weight with 3 bits. We only consider group size 128,
because the performance drops a lot without group size and usually could not
be deployed in practice. Similar to 4 bits evaluation, we report the lm-eval
result and wikitext2 ppl result in table 3 and 4 respectively. TEQ outperforms
RTN in all scenarios and is inferior to GPTQ on certain models. However, TEQ
could bring improvement for 8 out of 12 scenarios if taking Ours+GPTQ into
account.

#### Quantization Time.

We report the quantization time in Table 5. We adopt
Deepspeed333https://github.com/microsoft/DeepSpeed for 10B+ models due to the
potential out-of-memory (OOM) issue. As TEQ needs training, our time cost is
reasonably higher than GPTQ, especially when the model does not fit into the
device memory. It’s possible to reduce the time further by using more
resources or optimizing the code, while it’s out of scope.

Models | GPTQ | Ours
---|---|---
OPT-6.7B | 841 | 1239
OPT-13B | 1523 | 8737*
BLOOM-3B | 345 | 506
BLOOM-7B1 | 661 | 1148
LLAMA-7B | 712 | 1249
LLAMA-13B | 1240 | 9501*
Table 5: Quantization time in seconds for 4-bit weight quantization. * denotes
DeepSpeed is adopted in training for 10B+ models.

#### Analysis of Scales in TEQ.

We visualize the magnitude distribution histograms of $s_{l}$ for the layers
to which TEQ can be applied. Figure 2 displays the results of models with
$s_{l}$ initialized as scalar ones. Several conclusions can be drawn from
these results. Most notably, the majority of trained scales remain close to
their initial values (e.g., 1), typically within the range of [0.75, 1.25].
This suggests that even minor changes to the model can significantly reduce
quantization loss. Additionally, some scales deviate considerably from 1,
indicating the presence of “outlier” channels. Furthermore, scales in middle
layers tend to remain closer to their initial values compared to other layers,
suggesting that the first and last layers are more sensitive to the
quantization loss. We also attach results of scales initialized with
$1.0/sqrt(w_{cin})$ in Appendix A.5.

Figure 2:  The magnitude distributions of scales in TEQ for BLOOM-3B,
BLOOM-7.1B, OPT-6.7B, LLAMA-7B. The quantization configurations are w3_g128,
w4_g128, w4, and w4 respectively. Different colors refer to layer indices in
models (blue stands for shallow layers which are close to the data layer,
while red stands for deeper layers).

##  5 Conclusion

In this paper, we propose TEQ, a trainable equivalent transformation that
preserves the FP32 precision of the model output while also taking advantage
of low-precision quantization, and its training process is lightweight. Plus,
TEQ is regarded as orthogonal support for other quantization methods to
improve their performance. Our task-agnostic experiments and comparison with
other methods show that TEQ or its combination with other methods can obtain
comparable or better results.

###  5.1 Limitations

We find that the required memory during training is still high, though the
number of training parameters remains low. Moreover, since we enforce the
transformation to be equivalent and keep the architecture and FP32 output
unchanged, our results in some scenarios are inferior to the SOTA methods,
which could be fixed by combining the SOTA methods.

###  5.2 Ethics Statement

We propose TEQ for LLMs quantization. The method can be either used
individually or combined with other quantization methods. Since TEQ only
requires a few steps of finetuning on original models. Thus, it is safe to say
that TEQ’s technical details have no significant ethical implications. Our
work provides an exploration of large language model quantization through
simple finetuning, making their application easier. We believe increasingly
more work like this will emerge, making LLMs’ quantization more powerful.

## References

  * Ba et al. (2016) Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. 2016.  Layer normalization.  _arXiv preprint arXiv:1607.06450_.
  * Bai et al. (2020) Haoli Bai, Wei Zhang, Lu Hou, Lifeng Shang, Jing Jin, Xin Jiang, Qun Liu, Michael Lyu, and Irwin King. 2020.  Binarybert: Pushing the limit of bert quantization.  _arXiv preprint arXiv:2012.15701_.
  * Bengio et al. (2013) Yoshua Bengio, Nicholas Léonard, and Aaron Courville. 2013.  Estimating or propagating gradients through stochastic neurons for conditional computation.  _arXiv preprint arXiv:1308.3432_.
  * Bisk et al. (2020) Yonatan Bisk, Rowan Zellers, Jianfeng Gao, Yejin Choi, et al. 2020.  Piqa: Reasoning about physical commonsense in natural language.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 34, pages 7432–7439.
  * Dettmers et al. (2022) Tim Dettmers, Mike Lewis, Younes Belkada, and Luke Zettlemoyer. 2022.  Llm. int8 (): 8-bit matrix multiplication for transformers at scale.  _arXiv preprint arXiv:2208.07339_.
  * Esser et al. (2019) Steven K Esser, Jeffrey L McKinstry, Deepika Bablani, Rathinakumar Appuswamy, and Dharmendra S Modha. 2019.  Learned step size quantization.  _arXiv preprint arXiv:1902.08153_.
  * Frantar and Alistarh (2022) Elias Frantar and Dan Alistarh. 2022.  Optimal brain compression: A framework for accurate post-training quantization and pruning.  _arXiv preprint arXiv:2208.11580_.
  * Frantar et al. (2022) Elias Frantar, Saleh Ashkboos, Torsten Hoefler, and Dan Alistarh. 2022.  Gptq: Accurate post-training quantization for generative pre-trained transformers.  _arXiv preprint arXiv:2210.17323_.
  * Gao et al. (2021) Leo Gao, Jonathan Tow, Stella Biderman, Sid Black, Anthony DiPofi, Charles Foster, Laurence Golding, Jeffrey Hsu, Kyle McDonell, Niklas Muennighoff, Jason Phang, Laria Reynolds, Eric Tang, Anish Thite, Ben Wang, Kevin Wang, and Andy Zou. 2021.  A framework for few-shot language model evaluation.
  * (10) Github.  Github: Copilot.
  * Gong et al. (2019) Ruihao Gong, Xianglong Liu, Shenghu Jiang, Tianxiang Li, Peng Hu, Jiazhen Lin, Fengwei Yu, and Junjie Yan. 2019.  Differentiable soft quantization: Bridging full-precision and low-bit neural networks.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 4852–4861.
  * Hassibi et al. (1993) Babak Hassibi, David G Stork, and Gregory J Wolff. 1993.  Optimal brain surgeon and general network pruning.  In _IEEE international conference on neural networks_ , pages 293–299. IEEE.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2016.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778.
  * Howard et al. (2017) Andrew G Howard, Menglong Zhu, Bo Chen, Dmitry Kalenichenko, Weijun Wang, Tobias Weyand, Marco Andreetto, and Hartwig Adam. 2017.  Mobilenets: Efficient convolutional neural networks for mobile vision applications.  _arXiv preprint arXiv:1704.04861_.
  * Hubara et al. (2021) Itay Hubara, Yury Nahshan, Yair Hanani, Ron Banner, and Daniel Soudry. 2021.  Accurate post training quantization with small calibration sets.  In _International Conference on Machine Learning_ , pages 4466–4475. PMLR.
  * Ioffe and Szegedy (2015) Sergey Ioffe and Christian Szegedy. 2015.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International conference on machine learning_ , pages 448–456. pmlr.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_.
  * LeCun et al. (1989) Yann LeCun, John Denker, and Sara Solla. 1989.  Optimal brain damage.  _Advances in neural information processing systems_ , 2.
  * Li et al. (2022) Zhengyi Li, Cong Guo, Zhanda Zhu, Yangjie Zhou, Yuxian Qiu, Xiaotian Gao, Jingwen Leng, and Minyi Guo. 2022.  Efficient activation quantization via adaptive rounding border for post-training quantization.  _arXiv preprint arXiv:2208.11945_.
  * Lin et al. (2023) Ji Lin, Jiaming Tang, Haotian Tang, Shang Yang, Xingyu Dang, and Song Han. 2023.  Awq: Activation-aware weight quantization for llm compression and acceleration.  _arXiv_.
  * Louizos et al. (2018) Christos Louizos, Matthias Reisser, Tijmen Blankevoort, Efstratios Gavves, and Max Welling. 2018.  Relaxed quantization for discretized neural networks.  _arXiv preprint arXiv:1810.01875_.
  * Marcus et al. (1994) Mitch Marcus, Grace Kim, Mary Ann Marcinkiewicz, Robert MacIntyre, Ann Bies, Mark Ferguson, Karen Katz, and Britta Schasberger. 1994.  The penn treebank: Annotating predicate argument structure.  In _Human Language Technology: Proceedings of a Workshop held at Plainsboro, New Jersey, March 8-11, 1994_.
  * Merity et al. (2016) Stephen Merity, Caiming Xiong, James Bradbury, and Richard Socher. 2016.  Pointer sentinel mixture models.  _arXiv preprint arXiv:1609.07843_.
  * Nagel et al. (2020) Markus Nagel, Rana Ali Amjad, Mart Van Baalen, Christos Louizos, and Tijmen Blankevoort. 2020.  Up or down? adaptive rounding for post-training quantization.  In _International Conference on Machine Learning_ , pages 7197–7206. PMLR.
  * (25) OpenAI.  Openai: Chatgpt.
  * Paperno et al. (2016) Denis Paperno, Germán Kruszewski, Angeliki Lazaridou, Quan Ngoc Pham, Raffaella Bernardi, Sandro Pezzelle, Marco Baroni, Gemma Boleda, and Raquel Fernández. 2016.  The lambada dataset: Word prediction requiring a broad discourse context.  _arXiv preprint arXiv:1606.06031_.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J Liu. 2020.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _The Journal of Machine Learning Research_ , 21(1):5485–5551.
  * Sakaguchi et al. (2021) Keisuke Sakaguchi, Ronan Le Bras, Chandra Bhagavatula, and Yejin Choi. 2021.  Winogrande: An adversarial winograd schema challenge at scale.  _Communications of the ACM_ , 64(9):99–106.
  * Scao et al. (2022) Teven Le Scao, Angela Fan, Christopher Akiki, Ellie Pavlick, Suzana Ilić, Daniel Hesslow, Roman Castagné, Alexandra Sasha Luccioni, François Yvon, Matthias Gallé, et al. 2022.  Bloom: A 176b-parameter open-access multilingual language model.  _arXiv preprint arXiv:2211.05100_.
  * Shen et al. (2021) Mingzhu Shen, Feng Liang, Ruihao Gong, Yuhang Li, Chuming Li, Chen Lin, Fengwei Yu, Junjie Yan, and Wanli Ouyang. 2021.  Once quantization-aware training: High performance extremely low-bit architecture search.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 5340–5349.
  * Touvron et al. (2023) Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, et al. 2023.  Llama: Open and efficient foundation language models.  _arXiv preprint arXiv:2302.13971_.
  * Wei et al. (2022) Xiuying Wei, Ruihao Gong, Yuhang Li, Xianglong Liu, and Fengwei Yu. 2022.  Qdrop: randomly dropping quantization for extremely low-bit post-training quantization.  _arXiv preprint arXiv:2203.05740_.
  * Xiao et al. (2022) Guangxuan Xiao, Ji Lin, Mickael Seznec, Julien Demouth, and Song Han. 2022.  Smoothquant: Accurate and efficient post-training quantization for large language models.  _arXiv preprint arXiv:2211.10438_.
  * Yao et al. (2022) Zhewei Yao, Reza Yazdani Aminabadi, Minjia Zhang, Xiaoxia Wu, Conglong Li, and Yuxiong He. 2022.  Zeroquant: Efficient and affordable post-training quantization for large-scale transformers.  _Advances in Neural Information Processing Systems_ , 35:27168–27183.
  * Zellers et al. (2019) Rowan Zellers, Ari Holtzman, Yonatan Bisk, Ali Farhadi, and Yejin Choi. 2019.  Hellaswag: Can a machine really finish your sentence?  _arXiv preprint arXiv:1905.07830_.
  * Zhang et al. (2022) Susan Zhang, Stephen Roller, Naman Goyal, Mikel Artetxe, Moya Chen, Shuohui Chen, Christopher Dewan, Mona Diab, Xian Li, Xi Victoria Lin, et al. 2022.  Opt: Open pre-trained transformer language models.  _arXiv preprint arXiv:2205.01068_.
  * Zhang et al. (2020) Wei Zhang, Lu Hou, Yichun Yin, Lifeng Shang, Xiao Chen, Xin Jiang, and Qun Liu. 2020.  Ternarybert: Distillation-aware ultra-low bit bert.  _arXiv preprint arXiv:2009.12812_.
  * Zhuang et al. (2021) Bohan Zhuang, Mingkui Tan, Jing Liu, Lingqiao Liu, Ian Reid, and Chunhua Shen. 2021.  Effective training of convolutional neural networks with low-bitwidth weights and activations.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 44(10):6140–6152.

##  Appendix A Appendix

###  A.1 Additional comparison with AWQ

Although both AWQ and TEQ use a small calibration set from Pile, TEQ’s
evaluation methodology closely follows that of GPTQ and only shares a few
common tasks with AWQ. It is important to acknowledge that this comparison
inherently lacks rigor due to our reliance on referencing AWQ’s data alone.
Consequently, this approach introduces the potential unfairness in the
evaluation process, primarily stemming from the utilization of different
calibration datasets.

LLaMA-7B | AWQ | Ours
---|---|---
nbits | Method | PIQA | Hella. | Wino. | PIQA | Hella. | Wino.
16 | FP16 | 78.35 | 56.44 | 67.09 | 78.35 | 56.42 | 66.85
W3G128 | RTN | 75.84 | 53.10 | 63.22 | 75.68 | 53.18 | 63.06
GPTQ | 70.89 | 46.77 | 60.93 | 72.58 | 47.10 | 59.91
Proposed | 76.66 | 53.63 | 66.14 | 76.01 | 53.30 | 63.06
W4G128 | RTN | 77.86 | 55.81 | 65.59 | 77.58 | 55.91 | 65.59
GPTQ | 77.20 | 53.98 | 65.67 | 77.58 | 55.83 | 66.54
Proposed | 78.07 | 55.76 | 65.82 | 78.02 | 55.76 | 66.54
Table 6: Reported results of AWQ and Ours

Table 6 presents the LLaMA-7B’s results of our common tasks alongside AWQ in
table below and all the results of AWQ are from their paper.

###  A.2 Additional comparison with GPTQ act-order

We show the results in Table 7. TEQ still outperforms GPTQ in most cases.

nbits / gs | Methods | lm-eval ($\uparrow$) |  wikitext2 ppl ($\downarrow$)
---|---|---|---
4 / -1 | GPTQ-AO | 0.6713 | 6.06
Ours | 0.6771 | 6.30
Ours+GPTQ-AO | 0.6736 | 6.03
4 / 128 | GPTQ-AO | 0.6809 | 5.82
Ours | 0.6813 | 5.97
Ours+GPTQ-AO | 0.6811 | 5.82
3 / 128 | GPTQ-AO | 0.6042 | 8.29
Ours | 0.6521 | 6.89
Ours+GPTQ-AO | 0.6647 | 6.61
Table 7: Comparing results of Llama-7B for GPTQ with act-order. AO denotes
act-order. TEQ still outperforms GPTQ in most cases.

###  A.3 Special hyperparameters and settings

Usually, we adopt the same hyperparameters mentioned in section 3.2. So we
only list all the particular settings in Table 8.

lr | initialization | Models
---|---|---
default | $1.0/sqrt(w_{cin})$ |  | opt13b_w4; bloom3b_w4;
---
bloom7b_w4;
opt6.7b_w4_g128;
opt13b_w3_g128
4e-4 | default | llama7b_w4;
2e-4 | default | llama13b_w4_g128;
Table 8: Special hyperparameters and settings. g denotes group size Figure 3:
TEQ’s trained transformation parameters’ magnitude distributions, using
maximum’s square root value for initialization. From top to down are BLOOM-3B,
BLOOM-7.1B, OPT-6.7B and LLAMA-7B respectively.

###  A.4 More visualization results for TEQ’s trained parameters.

Figure 3 shows the magnitude distribution of scales initialized with
$1.0/sqrt(w_{cin})$. Since the initial value is related to channel-wise
maximum values, it’s more challenging to analyze. However, some outliers could
be still observed.

###  A.5 Counts of trainable parameters introduced by TEQ

We provide more details about counts of trainable parameters introduced by TEQ
in Table 9. table presented below offers details regarding the applicable
layers of TEQ in several models. We handle linear layers that possess
transformation scales that can be assimilated by their preceding layers, such
as Layer Normalization, among others.

As an illustration, within a single transformer block of OPT-6.7B, the QKV
layers have the same preceding layers and therefore utilize the same set of
trainable parameters. Based on the statistics, we have observed that TEQ’s
training only requires a minimal number of parameters (around the order from
1e-5 to 1e-4), thereby making our approach light-weighted enough.

Models |  Blocks |  TEQ Applicable Linear Layers |  Total Linear Layers |  TEQ Parameter Groups |  TEQ Parameters Counts |  Total Parameters Counts |  Ratio TEQ params and Total Params
---|---|---|---|---|---|---|---
Bloom 3B | 30 | 60 | 121 | 60 | 153600 | 3644810240 | 0.00421%
Bloom 7B1 | 30 | 60 | 121 | 60 | 245760 | 8096620544 | 0.00304%
OPT 6.7B | 32 | 160 | 193 | 72 | 786432 | 6864388096 | 0.01146%
OPT 13B | 40 | 200 | 241 | 96 | 1228800 | 13110865920 | 0.00937%
Llama 7B | 32 | 160 | 225 | 64 | 262144 | 6738415616 | 0.00389%
Llama 13B | 40 | 200 | 281 | 80 | 409600 | 13015864320 | 0.00315%
Table 9: Analysis of TEQ Parameters. TEQ only require a minimal ratio of
original models’ parameters (around the order from 1e-5 to 1e-4).