# Optimal Deployment of Multistatic Radar System Using Multi-Objective
Particle Swarm Optimization

Yichuan Yang, Tianxian Zhang, Wei Yi, Lingjiang Kong, Xiaolong Li, Bing Wang,
Xiaobo Yang Y. Yang, T. Zhang, W. Yi, L. Kong (Corresponding author), X. Li,
B. Wang and X. Yang are with University of Electronic Science and Technology
of China, School of Electronic Engineering, Chengdu, Sichuan, P.R. China,
611731. Fax: +86-028-61830064, Tel: +86-028-61830768, E-mail:
<EMAIL_ADDRESS>
###### Abstract

We consider an optimization deployment problem of multistatic radar system
(MSRS). Through the antenna placing and the transmitted power allocating, we
optimally deploy the MSRS for two goals: 1) the first one is to improve the
coverage ratio of surveillance region; 2) the second goal is to get a even
distribution of signal energy in surveillance region. In two typical working
modes of MSRS, we formulate the optimization problem by introducing two
objective functions according to the two mentioned goals, respectively.
Addressing on two main challenges of applying multi-objective particle swarm
optimization (MOPSO) in solving the proposed optimization problem, we propose
a deployment algorithm based on multi-objective particle swarm optimization
with non-dominated relative crowding distance (MOPSO-NRCD). For the challenge
of value difference, we propose a novel selection method with a non-dominated
relative crowding distance. For the challenge of particle allocation, a multi-
swarm structure of MOPSO is also introduced. Finally, simulation results are
given out to prove the advantages and validity of the proposed deployment
algorithm. It is shown that with same number of employed particles, the
proposed MOPSO-NRCD algorithm can achieve better optimization performance than
that of traditional multi-objective particle swarm optimization with crowding
distance (MOPSO-CD).

###### Index Terms:

multistatic radar system; multi-objective particle swarm optimization; optimal
deployment.

##  I Introduction

Taking advantages of higher amount of information available, multistatic radar
system (MSRS) has been proved with ability to improve performance in many
aspects (e.g., detection and estimation) [1, 2, 3, 4, 5, 6, 7]. Meanwhile,
such radar system has potential capability to be effective measures to counter
the threats of stealth targets and anti-radar missiles in military
applications [8]. Comparing with the conventional radar with a single transmit
antenna, MSRS provides more degrees of freedom, which support flexible
management modes, lead to further improved on many aspect of performance [3].
In fact, the management of MSRS, which can be divides into many modes, not
only improves but also affects the performance of MSRS to a larger extend.
Therefore, the management of MSRS should be optimized in order to get desire
performances.

Within all the management modes, the deployment of MSRS, which involved with
massive parameter settings and adjustment (e.g., placing the dispersed
antennas and allocating transmitted power), is of great importance to be a
basis that affect the performance of MSRS to a large extent [9, 10, 11, 12,
13, 14]. Specifically, the antenna optimal placement [9, 11, 13, 14] and power
optimal allocation [10, 12, 13], which are two main aspects of MSRS
deployment, have been proved with the ability to improve the localization [9,
10, 11, 12] and detection [13, 14] performances of MSRS. For the purpose of
providing a proper basis for MSRS to fully explore its potential (e.g., better
performance and more flexibility), the optimal deployment of MSRS should be
ascertained. As we know, the establishment of the optimal deployment can be
translated to an optimization problem with two key elements such as the
optimization model and the optimization algorithm. Many investigations have
been done on these two elements.

###  I-A Overview

For the optimization model: The optimization model is important as it
determines the deployment purpose and variables that can be adjusted.
Specifically, by employing different performance metrics as the objectives,
different optimal deployments of MSRS can be obtained. Some published works
focus on the localization performance metrics [9, 10, 11, 12]. The derived
Cramér-Rao bound (CRB) on velocity estimation of target [9] and the target
location estimation accuracy of an ultrawideband multiple-input-multiple-
output (MIMO) noise radar system [11] are introduced as objectives to explore
the optimal antenna placement of MSRS. In [10] and [12], the optimal power
allocations are obtained by considering two objectives, i.e., the CRB on
target position estimation and the posterior CRB on target state estimation.
Besides, detection performance metrics are also been taken into account, where
a weighted average detection performance and a coverage ratio are respectively
adopted as objectives in [13, 14]. For the investigating of MSRS deployment
problem, the above mentioned works just focus on one objective. However, radar
systems usually need to simultaneously adapt to different situations, wherein
different performance metrics are mainly concerned. Therefore, the deployment
of MSRS need to simultaneously take more than one objective into account.

For the optimization algorithm: The optimization algorithm is another
important issue since it determines the optimization performance. In [9], the
optimal antenna placement is obtained from an derived analytical expression of
CRB on target velocity estimation. For the purpose of achieving performance
better than a predefined threshold by optimally allocating the power, a
relaxation and domain decomposition method is introduced in [10]. In addition,
a geometrical and numerical method is employed to improve the estimation
accuracy of target location in [11].

The above mentioned works only consider the situation that a target locates in
one specific resolution cell. In real application, instead of focusing on one
single resolution cell, it is more practical to consider the MSRS deployment
problem for a region which contains multiple resolution cells. However, by
considering a region with multiple resolution cells and taking some
performance metrics as objectives, the deployment problems of MSRS are
normally of great computational complexity (i.e., huge computational load).
Considering the detection performance of MSRS in a region, a sequentially
exhaustive enumeration (SEE) method is introduced in [13]. While the SEE is
still suffered the challenge of a huge computational load, especially when the
concerned region contains many resolution cells and there are many variables
to be adjusted (i.e., many antenna to be placed). To address this problem, we
propose a particle swarm optimization (PSO) based placement algorithm of MSRS
in our previous work [14]. The computational load can be significantly
reduced.

However, the above works pay no attention to a more general situation, wherein
more than one objective are simultaneously considered. Therefore, the
optimization problem becomes a multi-objective problem (MOP), which is
normally more complex than single-objective optimization problem [15, 16, 17].
Recently, to expand the PSO algorithm in solving MOP, i.e., multi-objective
particle swarm optimization (MOPSO), has got an increasing interest in both
theoretical research [29, 30, 31, 32, 33] and real application [35, 36, 37,
38]. For the researches on real applications, some works take multiple antenna
metrics (e.g., the side lobe levels, antenna gains, the reflection coefficient
over the operating frequency band, and so on) as objectives in [35, 36]. In
[38], the coverage and lifetime are taken into consideration to study the
deployment of sensor networks. Nevertheless, none of these works consider a
practical and common case, in which the concerned objectives are significant
different to each other (e.g., values of multiple objective functions).

###  I-B Motivation

The motivations of the current paper are respectively focus on the two above
aspects:

1) In real radar applications, radar systems usually need to adapt to
different situations in a working environment of continuous changes, wherein
different tasks and signal processings are required. Thus, to enhance the
utilization efficiency of resources and flexibility of system, the deployment
of MSRS should provide a basis contains improvement of objectives
characterized by different performance metrics, which can evaluate or affect
the performance of different radar tasks.

2) Considering the performances of MSRS in a region, the objectives of the
optimization deployment problem may be significantly different to each other.
Then, the selection method of global best solution, which plays a key role to
determines the optimization performance of MOPSO to a lager extend, becomes
challenging since it is easily to be effected by the concerned objectives.
Based on the concept of nearest neighbor density estimation [33], a typical
and widely used method with crowding distance is introduced to select the
global best solution in order to guide the search of MOPSO [32]. While
obtained by directly summing all components (i.e., sub-crowding-distances in
the dimensionality of all objectives), the crowding distance can be called as
absolute crowding distance, substantially. However, direct sum of all
components is mainly effected by the components whose values are much larger
than the rests (e.g., the differences between values of multiple objective
functions could be several orders of magnitude). Therefore, the absolute
crowding distance can not adequately show the true crowding information in
some practical situations, wherein the values of multiple objective functions
as well as the sub-crowding-distances of objectives are significantly
different. Such problem is very common in real application including the
deployment of MSRS, while there is no effective method to solve it.

###  I-C Original Contributions

The original contributions of this paper are as follows.

1) We propose a general optimization deployment model, in which different
adjustable variables and objectives can be taken into account for similar
systems (e.g., sensor networks) with different purposes. Then, for the purpose
of improving the surveillance performance of MSRS equipped with different
detection methods, we introduce a coverage ratio and a lowest ratio of total
signal energy to noise power as objective functions to formulate a specific
optimization problem, wherein the antenna positions and transmitted power are
considered as adjustable variables for deployment of MSRS.

2) Since the values of multiple objective functions are commonly of
significant difference in real applications, we propose a non-
dominated111Taking maximization problem as example, given two vectors
$\boldsymbol{\alpha},\boldsymbol{\beta}$, if
$\boldsymbol{\alpha}\geqslant\boldsymbol{\beta}$ for
$\alpha_{k}\geqslant\beta_{k},k=1,2,\dots,K$. Then, we say that
$\boldsymbol{\alpha}$ dominates $\boldsymbol{\beta}$ (denotes by
$\boldsymbol{\alpha}\prec\boldsymbol{\beta}$) if
$\boldsymbol{\alpha}\geqslant\boldsymbol{\beta}$ and
$\boldsymbol{\alpha}\neq\boldsymbol{\beta}$ [33]. relative crowding distance
based method to properly select global best solution. Specifically, we firstly
select candidates from the non-dominated solutions. The sub-crowding-distances
of these selected candidates in dimensionality of all concerned objectives are
non-dominated. Then, we calculate the relative crowding distance by
normalizing and summing all the sub-crowding-distances for each candidate
solution. Finally, global best solutions will be selected out from the
candidates according to their relative crowding distances.

###  I-D Organization

The rest of this paper is organized as follows. In Section II, we formulate
the optimization problem for MSRS deployment by choosing two typical working
modes and proposing two objective functions according to the concerned
objectives. Then, focusing on the challenge that the values of multiple
objective functions are of great difference, we propose a novel MOPSO
algorithm based on non-dominated relative crowding distance in Section III.
The presented method is demonstrated by using some numerical examples in
Section IV. Finally, the conclusion is given in Section V.

##  II Problem Formulation

Figure 1: Sketch of surveillance region, placement region and effective
coverage region.

Though with difficulty due to the setting and adjustments of massive variables
(e.g., placing antennas and allocating transmitting power), the deployment is
of a great importance to be a foundation of MSRS and a determination for the
performances of MSRS to a large extent. For the purpose of fully exploring the
potential and achieving the desirable performances, the optimal deployment
should be established. Such establishment could be usually formulated as an
optimization problem in a general form:

$\displaystyle\hat{\mathbf{\Theta}}=\text{arg}\max\limits_{\mathbf{\Theta}}\mathbf{G}(\mathbf{\Theta}):=[g_{1}(\mathbf{\Theta}),g_{2}(\mathbf{\Theta}),\dots,g_{K}(\mathbf{\Theta})],$
(1)

where $g_{k}(\mathbf{\Theta})$ ($k=1,2,\dots,K$) denotes the $k$th objective
characterised by some performance metrics, $K$ is the number of concerned
objectives. $\mathbf{\Theta}$ is a vector that contains the adjustable
variables for MSRS deployment. We refer to $\mathbf{\Theta}$ as a deployment
vector (DV). We also refer to $\mathbf{G}(\mathbf{\Theta})$ as a objective
vector (OV) which contains all the concerned objectives. In practical
applications, radar systems usually need to face the requirement of different
tasks whose performances are evaluated and effected by different metrics. To
fully consider all the performance metrics may make the problem too complex to
be solved. Moreover, in real radar application, instead of taking all missions
into account at the same period, radar systems usually just need to focus on
the major missions and ignore the insignificant missions. Thus, in this paper,
we focus on the radar detection performance to study the MSRS deployment
problem. Meanwhile, our optimization model and algorithms can be easily
expanded to other radar applications or missions.

###  II-A Selection of Objectives and Optimization Parameters

Generally, a common duty of radar systems in both civil and military
applications is monitoring a surveillance region with targets of interest or
potential threats. To fulfill such duty, radar systems should be capable to
reliably detect targets in arbitrary resolution cell of a surveillance region.
As we all know, the effectiveness of target detection can be usually regarded
as that the detection probability ($P_{d}$) of a target in a resolution cell
is higher than a preset threshold ($P_{dt}$) [18]. Then, the radar
surveillance performance is determined by the effectiveness of target
detection for all the resolution cells within the surveillance region.
Similarly, surveillance performance of MSRS is also determined by target
detection. Thus, we need to find proper objectives to optimally deploy the
MSRS. Similar to our previous work [14], we select the coverage ratio
($\mathcal{C_{R}}$) as one objective of the optimal MSRS deployment problem.

To better understand, the introduction of $\mathcal{C_{R}}$ begins with the
definition of effective coverage region (ECR). If the $P_{d}$ of a target in
some resolution cells is higher than $P_{dt}$ (i.e., $P_{d}\geqslant P_{dt}$),
we call these cells as effectively covered resolution cell. The ECR is
constructed by all effectively covered resolution cells. Under a commonly
setting of MSRS, resolution cells of MSRS are usually approximated as square
grids [3]. Let $\mathbb{A}$ and $\mathbb{C}$ denote the surveillance region
and the effective coverage region, respectively (see Fig.1). Then, we define
the coverage ratio of surveillance region ($\mathcal{C}_{\mathcal{R}}$) as:

$\displaystyle\
\mathcal{C}_{\mathcal{R}}=\frac{\textrm{area}(\mathbb{C}\cap\mathbb{A})}{\textrm{area}(\mathbb{A})}.$
(2)

Obviously, a bigger $\mathcal{C}_{\mathcal{R}}$ represents a better reliable
detection performance and the radar system will effectively monitor the whole
surveillance region (i.e., $\mathbb{A}$) when the $\mathcal{C}_{\mathcal{R}}$
reaches $1$ (i.e., $\mathbb{C}\cap\mathbb{A}=\mathbb{A}$). The goal of
$\mathcal{C_{R}}=1$ may be not difficult to achieve for a strong target.
However, such goal is extremely hard for a small target. It is because that
the signal energy to noise power ratio of the small target is usually very
low, which decreases its detection probability as well as the
$\mathcal{C_{R}}$. To address this problem, some detection algorithms have
been proposed by increasing the complexity of signal processing, such as
track-before-detect (TBD) [41, 42] and adaptive waveform design technologies
[43, 44]. However, due to lacking real-time capabilities [42], small
application area [44] and huge computation complexity [41, 43, 44], these
detection algorithms is not suitable for real time monitoring of a
surveillance region. Recently, we have proposed a method, which can
significantly increase the $\mathcal{C_{R}}$ by optimally placing the antennas
of MSRS [14]. Nevertheless, as just only focusing on $\mathcal{C_{R}}$, the
signal energy will be unevenly distributed in the surveillance region. Thus,
in some parts of the surveillance region, the signal energy to noise power
ratio could be reduced to a excessively low level, which would degrade the
detection performance sharply. Therefore, some radar coverage gaps that can
not be effectively monitored will occur. Such phenomenon is usually
unacceptable to the surveillance requirement of a radar system.

To avoid this problem, we propose a solution contains two deployment goals to
improve the surveillance performance of radar system: 1) The first goal is to
enhance the $\mathcal{C_{R}}$ of radar system; 2) for the purpose of getting a
even distribution of signal energy in surveillance region, which can provide a
proper basis to facilitate the usage of the complex methods to achieve
desirable performance, the second goal is to roundly improve the ratio of
total signal energy to noise power (RTSN). To this end, based on worst-case
performance optimization, we choose to enhance the RTSN of the received
returning signal whose RTSN is the lowest among others. Assume there are $U$
resolution cells in surveillance region, the second objective can be expressed
as:

$\displaystyle\mathcal{L}_{\mathcal{R}}=\min\limits_{u=1,2,\dots,U}\chi_{u},$
(3)

where $\mathcal{L}_{\mathcal{R}}$ denotes the lowest RTSN of all resolution
cells in surveillance region, $\chi_{u}$ denotes the RTSN of the returning
signals from the $u$th resolution cell.

After the introduction of the goals, their corresponding objective functions
as well as the calculation methods are needed to build the problem, which will
be introduced in the following subsections.

###  II-B Introductions of Working Modes of MSRS

Due to the existing of various adjustable parameters, the MSRS is a system of
high degree of freedom and it has many different working modes, under which
the calculation methods of one objective could be different. According to the
reception mechanism of return signals, MSRS roughly have two typical working
modes and we briefly introduce them here.

####  II-B1 Cooperative Mode

In this mode, the signals transmitted by transmitters of different radar nodes
are orthogonal to each other, each receiver of radar nodes can receive and
separate the returning signals reflected by targets. All the received
returning signals will be jointly processed. Taking advantages of the
cooperation of dispersed radar nodes to improve performance in many aspects,
this working mode has received a significant attention in recent years[3, 4,
5, 7, 9, 10]. However, in real applications of such mode, the control and
synchronize of the widely dispersed radar nodes are difficult (e.g., the
problem of pulse chasing will occur when controlling all the radar nodes to
illuminate a region simultaneously [28]). Meanwhile, to satisfied the
requirement of cooperation (e.g., communication requirement), system load of
this working mode would be big.

####  II-B2 Non-cooperative Mode

In this mode, each radar node of MSRS works independently without cooperation
with others and it can only receive and process the signal transmitted by
itself. Though without performance improvement of cooperation, such working
mode does not have to settle the challenge of controlling and synchronization
of dispersed radar nodes. Meanwhile, the system resources (e.g., computation
and communication resources) consumption of this mode is much less than that
of cooperative mode, which facilitates its application in real application.
Thus, non-cooperative mode is an effective supplement to cooperative mode and
we need to take it into account when considering the optimal deployment of
MSRS.

###  II-C Introduction of Objective Functions

After the introduction of deployment goals ($\mathcal{C_{R}}$ and
$\mathcal{L_{R}}$), we still need their objective functions to solve the
problem. In this subsection, under the two typical working mode, we propose
the corresponding objective functions as well as their calculation methods for
the two deployment goals.

####  II-C1 Objective function of $\mathcal{C_{R}}$

Since the shape of ECR of MSRS is normally irregular, we have to analyse the
$P_{d}$ of a target for each resolution cell within a surveillance region in
order to get $\mathcal{C_{R}}$. Meanwhile, the surveillance region often
contains massive resolution cells and the calculation of $P_{d}$ of MSRS is
difficult in both theoretical research (e.g., complex integral) [6] and real
application (e.g., phase and amplitude information of different transmit-
receive (T-R) pairs are inconsistent) [3]. Therefore, the calculation of
$\mathcal{C_{R}}$ is usually of great computational complexity.

To reduce the computational complexity, the square-law detector (as a commonly
used non-coherent detector) is employed for MSRS. Specifically, square-law
detector firstly sums the square of samples of the received returning signals,
then compares it with a fixed threshold to detect the existence of targets
[18]. Such detector has two major advantages: 1) We can acquire a simple
$P_{d}$ calculation method of this detector under some reasonable assumptions.
2) Without considering the phase information, which avoid the difficulty of
phase aligning of signals of different T-R pairs. Therefore, the square-law
detector is easy to implement for real application.

Similar to our previous works [14], we consider the same signal model in [5]
and ignoring clutter. Under assumptions as follows: a MSRS constituted by $J$
monostatic radars, there are $U$ resolution cells in the surveillance region,
the radar cross section (RCS) of the target is unfluctuating, the background
noise is white Gaussian noise, and each radar only transmits one pulse in one
observation. For the cooperative mode, as all the received returning signals
are jointly processed, $P_{d}$ of target locates in the $u$th resolution cell
(whose position is $(x_{u},y_{u})$) under this mode can be calculated as [19,
20]:

$\displaystyle P_{d}(\mathbf{\Theta})$ $\displaystyle=Q_{J\times
J}(\sqrt{2\times\chi_{u}(\mathbf{\Theta})},\sqrt{2\times\gamma_{T}}),$ (4)
$\displaystyle P_{fa}$
$\displaystyle=e^{-\gamma_{T}}\sum\nolimits_{i=1}^{J\times
J-1}\gamma_{T}^{i}/i!,$

where $\chi_{u}(\mathbf{\Theta})$ denotes the RTSN of the returning signals
that determines the performance of MSRS in cooperative mode, $Q$ denotes the
Marcum Q function, $\gamma_{T}$ is the detection threshold which is set
according to false alarm probability ($P_{fa}$).

For the non-cooperative mode, all the radar nodes of MSRS work independently.
The receiver of each radar node only receives the returning signal that
transmitted by its own transmitter. Instead of jointly processing all the
received signals, each radar node only processes the returning signal received
by itself. All the radar nodes have neither cooperation nor communication to
each other. Specifically, in this mode, all radar nodes independently detect
the target and they all get their own $P_{d}$ of target, and thus the $P_{d}$
of MSRS is equal to that of the radar node whose $P_{d}$ of target is the
highest among all radar nodes. Therefore, $P_{d}$ of target locates in the
$u$th resolution cell can be calculated as:

$\displaystyle P_{d}(\mathbf{\Theta})$
$\displaystyle=Q(\sqrt{2\times\chi_{u}(\mathbf{\Theta})},\sqrt{2\times\gamma_{T}}),$
(5) $\displaystyle P_{fa}$ $\displaystyle=e^{-\gamma_{T}}\gamma_{T},$

where $\chi_{u}(\mathbf{\Theta})$ denotes the RTSN of the returning signals
that determines the performance of MSRS in non-cooperative mode. The
$\chi_{u}(\mathbf{\Theta})$ in different mode (i.e., cooperative and non-
cooperative mode) are different and we will introduce them as well as their
calculation methods in the second part of this subsection (i.e., Objective
function of $\mathcal{L_{R}}$).

Using the calculation method of $P_{d}$ for both the cooperative mode (i.e.,
use (4)) and the non-cooperative mode (i.e., use (5)), we can calculate the
$P_{d}$ of every resolution cell within the surveillance region. Assume that
we can get $U^{\prime}(\mathbf{\Theta})$ effectively covered resolution cells
(i.e., in which the $P_{d}(\mathbf{\Theta})$s are higher than $P_{dt}$) for a
given $\mathbf{\Theta}$ in the surveillance region. As mentioned before, there
are $U$ resolution cells in surveillance region in total. Then, we can obtain
the $\mathcal{C_{R}}$ as bellow:

$\displaystyle\mathcal{C}_{\mathcal{R}}(\mathbf{\Theta})$
$\displaystyle=\frac{\textrm{area}(\mathbb{C}(\mathbf{\Theta})\cap\mathbb{A})}{\textrm{area}(\mathbb{A})}=\frac{U^{\prime}(\mathbf{\Theta})}{U}.$
(6)

####  II-C2 Objective function of $\mathcal{L_{R}}$

Generally, a MSRS contains many T-R pairs [3]. For the target in the $u$th
resolution cell, the T-R pair that transmitted from the transmitter of the
$m$th radar node and received at the receiver of the $n$th radar node, the
RTSN of returning signal in this T-R pair is given as:

$\chi^{mn}_{u}=\frac{P_{tm}\tau
A_{tm}A_{rn}\sigma^{mn}\lambda^{2}F_{tm}^{2}F_{rn}^{2}}{{4\pi}^{3}KT_{s}(R_{tm}R_{rn})^{2}C_{B}L_{tm}L_{rn}},$
(7)

where $R_{tm}$, $R_{rn}$, $P_{tm}$, $\tau$, $A_{tm}$, $A_{rn}$, $\sigma^{mn}$,
$\lambda$, $F_{tm}$, $F_{rn}$, $C_{B}$, $T_{s}$, $K$, $L_{tm}$, $L_{rn}$ are
respectively transmitter-to-target range, receiver-to-target range,
transmitted power of the $m$th radar node, pulse width, transmit antenna power
gain, receive antenna power gain, bistatic radar target cross section,
wavelength, pattern propagation factor for transmitter-to-target path, pattern
propagation factor for target-to-receiver path, bandwidth correction factor,
receive system noise temperature, Boltzmann’s constant, transmit system
losses, receive system loss. $\chi^{mn}_{u}$ denotes the RTSN of returning
signal in this T-R pair.

Generally, most of these above radar parameters are constant and it is
unnecessary to take all of them into account. Besides, to concretely set up
all these involved parameters is beyond the main scope of this paper. We need
to find proper method to calculate the RTSN for each T-R pair conveniently.
Similar with our previous work [14], we use a method based on radar range
equation to approximatively calculate the RTSN.

Firstly, as we mentioned before, the concerned MSRS is assumed to be
constituted by the same kind of monostatic radars. It is reasonable assume
that the transmitter and receiver of each radar node share one antenna and all
the above radar parameters except transmitted power (we need to optimally
allocate them) are the same. Under such assumption, we can get that
$A_{tm}=A_{rn}=A$ and $L=L_{tm}L_{rn},\text{for}\hskip 4.2679pt1\leqslant
n,m\leqslant J$. Furthermore, we ignore the effect of pattern propagation
factors (i.e., $F_{tm}=F_{rn}=1,\text{for}\hskip 4.2679pt1\leqslant
n,m\leqslant J$). Therefore, the T-R pair observed at $n$th receiver due to
the transmission at $m$th transmitter is given by[18] (still use $u$th
resolution cell as example):

$\chi^{mn}_{u}=\frac{P_{tm}\tau
A^{2}\sigma^{mn}\lambda^{2}}{{4\pi}^{3}KT_{s}R_{tm}^{2}R_{rn}^{2}C_{B}L}.$ (8)

Then, it is well known that once the radar system parameters are established,
we can get a maximum detection range for each of the concerned radar node. The
maximum detection range can be given by radar range equation [18]:

$R_{\max}=[\frac{P_{t}\tau
A^{2}\sigma\lambda^{2}}{{4\pi}^{3}KT_{s}D_{0}C_{B}L}]^{\frac{1}{4}},$ (9)

where $\sigma$ is the radar target cross section and $R_{\max}$ is the maximal
detection range when radar node works independently, $D_{0}$ is the
detectability factor and $P_{t}$ is a standard transmitted power satisfied
with a fixed performance characterized by $R_{\max}$, $\sigma$, $P_{fa}$ and
$P_{dt}$.

Finally, the RTSN of returning signal in the T-R pair which observed at $n$th
receiver due to the transmission at $m$th transmitter can be obtained by
substituting (9) into (8) and with some simple derivation:

$\displaystyle\chi^{mn}_{u}$
$\displaystyle=D_{0}\frac{\rho_{m}\sigma_{m,n}(R_{\max})^{4}}{\sigma(R_{tm}R_{rn})^{2}},$
(10) $\displaystyle\rho_{m}$ $\displaystyle=\frac{P_{tm}}{P_{t}}.$ (11)

Instead of setting specific values, here we allocate the transmitted power by
adjusting a ratio ($\rho$) between the standard one ($P_{t}$) and the actual
transmit ones ($P_{tm}$s) of each radar node. For the cooperative mode, there
are $J\times J$ T-R pairs. This is because that each receiver of radar nodes
can receive and separate all the returning signals reflected by targets. As
jointly processing the received signals from all T-R pairs, the performance of
MSRS in this mode is effected by all the T-R pairs. Meanwhile, the noise power
of different T-R pairs are the same due to the same kinds of receiver in
different radar nodes. Then, the RTSN of the $u$th resolution cell can be
obtained as:

$\displaystyle\chi_{u}(\mathbf{\Theta})$
$\displaystyle=\frac{E_{S}^{t}}{\sigma_{0}}$
$\displaystyle=\frac{\sum\nolimits_{m=1}^{J}\sum\nolimits_{n=1}^{J}E_{S}^{mn}}{\sigma_{0}}$
$\displaystyle=\sum\nolimits_{m=1}^{J}\sum\nolimits_{n=1}^{J}\frac{E_{S}^{mn}}{\sigma_{0}}$
$\displaystyle=\sum\limits_{m=1}^{J}\sum\limits_{n=1}^{J}\chi^{mn}_{u},$ (12)

where $E_{S}^{t}$, $E_{S}^{mn}$ and $\sigma_{0}$ respectively denote the total
signal energy, the signal energy of a T-R pair (i.e., the pair transmitted by
$m$th node and received by $n$th node) and the noise power of each T-R pair.

For the non-cooperative mode, there are $J$ T-R pairs. This is because that
each radar node can only receive and process the signal transmitted by itself.
Since all radar nodes works independently without neither cooperation nor
communication to each other, the performance of the MSRS is determined by the
T-R pair whose own performance is better than that of other T-R pairs.
Therefore, we can obtain the $\chi_{r}$ as:

$\displaystyle\chi_{u}(\mathbf{\Theta})=\max\limits_{\scriptstyle{m=1,2,\dots,J}\atop\scriptstyle{n=m}}\chi^{mn}_{u},$
(13)

and then we can calculate $\mathcal{L_{R}}$ as:

$\displaystyle\mathcal{L_{R}}(\mathbf{\Theta})=\min\limits_{u=1,2,\dots,U}\chi_{u}(\mathbf{\Theta}).$
(14)

After the introduction of the objective functions, we give the specific form
of our problem as below:

$\displaystyle\hat{\mathbf{\Theta}}=\text{arg}\max\limits_{\mathbf{\Theta}}\mathbf{G}(\mathbf{\Theta}):=(g_{1}(\mathbf{\Theta}),g_{2}(\mathbf{\Theta})),$
(15) $\displaystyle
g_{1}(\mathbf{\Theta})=\mathcal{C}_{\mathcal{R}}(\mathbf{\Theta})=\frac{\textrm{area}(\mathbb{C}(\mathbf{\Theta})\cap\mathbb{A})}{\textrm{area}(\mathbb{A})},$
$\displaystyle
g_{2}(\mathbf{\Theta})=\mathcal{L}_{\mathcal{R}}(\mathbf{\Theta})=\min\limits_{r=1,2,\dots,R}\chi_{r}(\mathbf{\Theta}),$
$\displaystyle s.t.\quad\sum\limits_{j=1}^{J}P_{j}=P_{T}=JP_{t}.$

We remark that the definition of the deployment vector (DV) in this problem is
$\mathbf{\Theta}=[{\mbox{\boldmath{$\theta$}}}_{1},\dots,{\mbox{\boldmath{$\theta$}}}_{j},\dots,{\mbox{\boldmath{$\theta$}}}_{J},P_{1},\dots,P_{j},\dots,P_{J}]$,
where ${\mbox{\boldmath{$\theta$}}}_{j}=(x_{j},y_{j})\in\mathbb{B}$ and
$P_{j}$ denote the positions of antenna and the transmitted power of the $j$th
radar node, respectively. $P_{T}=JP_{t}$ denotes the fixed total transmitted
power of the whole radar system and $\mathbb{B}$ denotes the placement region
that antenna can be placed in. By optimally placing the antenna and allocating
the transmitted power of each radar node, the purpose of the optimization
problem is to simultaneously enhance the $\mathcal{C_{R}}$ and
$\mathcal{L_{R}}$, with limitation of number and total transmitted power of
all radar nodes.

Such problem is not easy to solve: Firstly, since we need to jointly consider
the positions of all antennas and the transmitted power of each radar, the
problem is of high dimensionality which would cause huge computational load
for the optimization processing. Secondly, with potential conflict and
differences between the two objectives, the optimization of this problem could
be more complicated than conventional single-objective ones. Focusing on these
issues, we propose an novel optimization algorithm based MOPSO to solve the
proposed problem in the next section.

##  III Introduction of MOPSO and A Novel MOPSO Algorithm

By utilizing a swarm of candidate solutions (named as particles) to search for
optimal solution, PSO has been proved its validity in many applications [21,
22, 23, 24, 25] since it is originally inspired by social behavior of bird
flocking [26]. Such successes have motivated researchers to extend the usage
of PSO into multi-objective optimization problems (MOP) [29, 30, 31, 32, 34,
33, 35, 36, 37, 38]. Nevertheless, to our best knowledge, the use of MOPSO for
the deployment of MSRS has barely been investigated. Meanwhile, in order to
adopt the MOPSO into real applications, some problems still remain to be
solved (e.g., proper nearest neighbor density estimation in situation that the
values of objective functions are significant different). In this section,
some brief introductions of basic theories of MOP and MOPSO will be given
first. Then, we will propose a novel MOPSO algorithm for the optimization
deployment of MSRS.

###  III-A Introduction of Multi-objective Particle Swarm Optimization

Obviously, the problem described in (15) is a multi-objective optimization
problem (MOP). Generally, the objectives (e.g., $\mathcal{C_{R}}$ and
$\mathcal{L_{R}}$) of MOP are usually in conflict with respect to each other,
which means the improvements in some objectives will probably result in
degradations of other objectives. Therefore, there is hardly any single ideal
optimal solutions for MOP [34]. Generally, based on the concept of Pareto-
optimal, a group of non-dominated solutions that represent the best possible
compromises among the objectives are needed to solve MOP [29, 30, 31, 32, 33,
34, 35, 36, 37, 38]. Specifically, we can get a OV (i.e.,
$\mathbf{G}(\mathbf{\Theta})$) once a DV (i.e., $\mathbf{\Theta}$) is
determined in the context of (15), and then the non-dominated solutions are
the DVs whose corresponding OVs are non-dominated. In the multi-objective
space, the non-dominated OVs constitute a front which is called as Pareto
front (denotes by $\mathcal{P_{F}}$). Many techniques (including PSO) are
adopted to solve MOPs in both theoretical research [29, 30, 31, 32, 33, 34]
and real applications [35, 36, 37, 38]. Under the background of the proposed
optimization problem, we will briefly introduce the principle of MOPSO and
then proposed a novel MOPSO algorithm.

To solve the proposed optimization problem by MOPSO, each particle stands for
a candidate DV and we suppose there are $S$ particles (means swarm size is
$S$). Then, the optimization process of MOPSO begins with a random
initialization of these $S$ DVs (i.e., the positions of antennas and
transmitted power of transmitters are randomly set up within fixed value
ranges). As commonly setting in the MOPSO, an external archive is set up to
store the non-dominated solutions (i.e., the non-dominated DVs) [31]. Then,
from the initial DVs, we select the satisfactory ones out and store them in
this external archive. Besides, the initial individual best solution of each
particle is set to be the initial position of itself, and then the global best
solution is selected from the initial external archive with some methods.
After the initialization, all particles will update themselves as follow
(without losing of generality, we using the $s$th particle as example, $1\leq
s\leq S$) in every iteration:

$\displaystyle\mbox{\boldmath{$v$}}^{s}(t+1)=$
$\displaystyle\omega(t)\times\mbox{\boldmath{$v$}}^{s}(t)+c_{1}\mbox{\boldmath{$r$}}_{1}(t)(\mbox{\boldmath{$p$}}^{s}(t)-\boldsymbol{\Theta}^{s}(t))$
$\displaystyle+c_{2}\mbox{\boldmath{$r$}}_{2}(t)(\mbox{\boldmath{$p$}}^{g}(t)-\boldsymbol{\Theta}^{s}(t)),$
(16) $\displaystyle\boldsymbol{\Theta}^{s}(t+1)=$
$\displaystyle\boldsymbol{\Theta}^{s}(t)+\mbox{\boldmath{$v$}}^{s}(t+1),$ (17)

where $\mbox{\boldmath{$v$}}^{s}$, $\boldsymbol{\Theta}^{s}$ and
$\mbox{\boldmath{$p$}}^{s}$ denote the velocity, the position (i.e., DV) and
the individual best solution of $s$th particle (i.e., the best DV proposed by
the $s$th particle so far), $\mbox{\boldmath{$p$}}^{g}$ denotes the global
best solution (i,e, the best DV proposed by all particles so far). $\omega(t)$
is inertia weight. $c_{1}$ and $c_{2}$ are acceleration constants.
$\mbox{\boldmath{$r$}}_{1}(t)$ and $\mbox{\boldmath{$r$}}_{2}(t)$ are vectors
whose elements are random independent variables in the range [0,1]. The
velocity of particle in each dimension is limited by a fixed $V_{max}$.

After the updating of position and velocity, we calculate the OV of each new
particle and then update the individual best solutions as (still using the
$s$th solution for instance):

$\displaystyle\mbox{\boldmath{$p$}}^{s}(t+1)=\left\\{\begin{array}[]{ll}\boldsymbol{\Theta}^{s}(t+1),&if\hskip
1.42262pt\mathbf{G}(\boldsymbol{\Theta}^{s}(t+1))\prec\mathbf{G}(\mbox{\boldmath{$p$}}^{s}(t)),\\\
\mbox{\boldmath{$p$}}^{s}(t),&others.\end{array}\right.$ (20)

Other particles will be updated in the same way. Afterwards, similar with
[31], we update the external archive, and then the global bests (i.e.,
$\mbox{\boldmath{$p$}}^{g}$) will be selected out from the updated external
archive with some selection methods. Such iteration process will continue
until $t$ reaches the $T_{\max}$, where the optimization process will stop and
the optimization results (i.e., the DVs that stored in the external archive)
will be output.

Generally, since with great importance to affect the search process as well as
the optimization performance, the selection of global best solution from the
set of non-dominated solutions becomes a key component of the MOPSO approach
designment [33]. Many selection methods are proposed in order to get better
optimization performance. Within all the published works, the selection method
with absolute crowding distance, which is based on the concept of neighbour
density estimation, has been proved with its validity in theoretical research
[32, 33]. Nevertheless, such method may not able to reasonably select out the
global best solution when the values of the objective functions are
significantly different (e.g., the difference could be several orders of
magnitude), which is a common phenomenon in real applications. Focusing on
this issue, we proposed a novel MOPSO algorithm with a non-dominated relative
crowding distance. Before the introduction of the proposed MOPSO, we would
briefly introduce the MOPSO with absolute crowding distance (MOPSO-CD)
firstly.

Figure 2: An example of the $i$th cuboid that belongs to
$\mathbf{G}(\mathbf{\Theta}_{i})$ in a two-objective space.

The absolute crowding distance of one solution (i.e., one DV) is the
estimation of the size of the largest cuboid (see Fig.2) in multi-objective
space. One cuboid can enclose one position point of a specific OV without
including any other position points. For instance, the absolute crowding
distance of the $i$th DV can be calculated as (suppose there are $I$ non-
dominated DVs stored in the external archive):

$\displaystyle\xi_{cd}^{i}$
$\displaystyle=\sum\nolimits_{k=1}^{K}\eta_{k}^{i},\quad 2\leqslant i\leqslant
I-1,$ (21) $\displaystyle\eta_{k}^{i}$
$\displaystyle=|g_{k}(\mathbf{\Theta}_{i+1})-g_{k}(\mathbf{\Theta}_{i-1})|,$
(22) $\displaystyle\xi_{cd}^{1}$
$\displaystyle=\xi_{cd}^{I}=\max\limits_{2\leqslant i\leqslant
I-1}\xi_{cd}^{i},$ (23)

where $\xi_{cd}^{i}$ denotes the absolute crowding distance of the $i$th DV
(i.e., $\mathbf{\Theta}_{i}$), $\eta_{k}^{i}$ denotes the sub-crowding-
distance in the dimension of $k$th objective. The DV with the largest absolute
crowding distance has the largest probability to be selected as the global
best [32]. As mentioned before, since obtained by directly summing all
components, the selection method with absolute crowding distance may face two
challenges in real applications:

The challenge of value difference: In real applications, the value of
different objective are usually of different orders of magnitude. In (15), the
value range of $\mathcal{C_{R}}$ (i.e., it is in the range of [0,1]) is much
smaller than that of $\mathcal{L_{R}}$ (e.g., the difference between its lower
bound and upper bound could be tens of dB). As obtained by directly summing of
all the sub-crowding-distances (i,e,. $\eta_{k}^{i}$), the absolute crowding
distance is mainly effected by the $\eta_{k}^{i}$s whose values are much
larger than the others. Therefore, the crowding distance information of the
$\eta_{k}^{i}$s whose values are small will be covered. Thus, the absolute
crowding distance is not adequately enough in such scenarios.

For better understanding, we use an example to illuminate this problem, in the
case of maximization problem which has two objectives ($g_{1}$ and $g_{2}$),
assume that the value range of $g_{1}$ is [0,1] and the the value range of
$g_{2}$ is [0,100]. We demonstrate an $\mathcal{P_{F}}$ (which is constituted
by 6 non-dominated solutions) and its projections at axis of $g_{1}$ and axis
of $g_{2}$ in Fig.3 (we scale down the axis of $g_{2}$ for convenient
demonstration).

Figure 3: An example of 6 position points characterized by 6 DVs in a two-
objective space, the two objectives are of different orders of magnitude: the
value range of $g_{1}$ is [0,1] and the the value range of $g_{2}$ is [0,100].

Considering about the value range of each objective and the intervals between
neighbouring solutions jointly, we can see that the distribution of the
projection on the axis of $g_{2}$ is of better uniformity than that of
$g_{1}$. Specifically, the $\mathbf{G(\Theta_{4})}$ in $\mathcal{P_{F}}$ is
obviously sparsest and we should search the space around $\mathbf{\Theta_{4}}$
to find new solutions. However, for the aspect of absolute crowding distance,
we can see that (in TABLE I) the $\mathbf{G(\Theta_{2})}$ has the biggest
absolute crowding distance (ignoring the $\mathbf{G(\Theta_{1})}$ and the
$\mathbf{G(\Theta_{6})}$, because they located on the edge of
$\mathcal{P_{F}}$ and their absolute crowding distance equal to
$\mathbf{G(\Theta_{2})}$ according to (23)), which means the space around
$\mathbf{\Theta_{2}}$ has the highest priority to be searched. Such result of
searching priority is unreasonable and it is due to the reason that the value
range of $g_{1}$ is much smaller than that of $g_{2}$ (i.e., at least two
orders of magnitude). Therefore, the absolute crowding distance as well as its
based selection method of global best solution are not proper in such kind of
real applications.

TABLE I: Absolute Crowding Distance of Solutions in Fig.3 |  $\mathbf{G(\Theta_{1})}$ |  $\mathbf{G(\Theta_{2})}$ |  $\mathbf{G(\Theta_{3})}$ |  $\mathbf{G(\Theta_{4})}$ |  $\mathbf{G(\Theta_{5})}$ |  $\mathbf{G(\Theta_{6})}$
---|---|---|---|---|---|---
$\xi^{i}_{cd}$ |  41.2 |  41.2 |  40.51 |  39.71 |  40.4 |  41.2

The challenge of particle allocation: Comparing with assigning a same mission
to all particles, the cooperation of particles based on labour division is
with potential to achieve better optimization performance. For instance, from
the principle of MOPSO-CD, since
$\xi_{cd}^{1}=\xi_{cd}^{I}=\max\limits_{2\leqslant i\leqslant
I-1}\xi_{cd}^{i}$ and in every iteration that the solution with the largest
absolute crowding distance value has the largest probability to be selected as
the global best, the particles of MOPSO-CD should search more solutions on the
explored $\mathcal{P_{F}}$ and extend the value range of each objective
function to get a bigger (wider) $\mathcal{P_{F}}$ jointly. Obviously, we can
divide the optimization goal into two parts (i.e., to find more solutions on
$\mathcal{P_{F}}$ and to get a bigger $\mathcal{P_{F}}$) and then allocate
different swarm of particles to take on these two parts of optimization task.

Addressing on both of these two challenges, we will propose a novel MOPSO
algorithm in the next subsection.

###  III-B MOPSO with Non-dominated Relative Crowding Distance

For the challenge of value difference: we proposed a selection method with
non-dominated relative crowding distance. The selection method contains three
steps. We will introduce the detail of them in the following.

Firstly: we calculate a crowding distance vector (CDV) and a relative crowding
distance for each solution stored in the external archive. For example, for
the $i$th solution (i.e., $\mathbf{\Theta}_{i}$), the calculation of its CDV
is given as follow:

$\mathbf{H}^{i}=[\eta_{1}^{i},\eta_{2}^{i},\dots,\eta_{K}^{i}]^{\text{T}},\\\
$ (24)

where $\mathbf{H}^{i}$ denotes the CDV of the $i$th solution,
$\eta_{k}^{1}=2\times|g_{k}(\mathbf{\Theta}_{2})-g_{k}(\mathbf{\Theta}_{1})|$
and
$\eta_{k}^{I}=2\times|g_{k}(\mathbf{\Theta}_{I})-g_{k}(\mathbf{\Theta}_{I-1})|$.
The calculation of relative crowding distance is given as follow (still using
the $i$th solution as example):

$\displaystyle\xi_{rcd}^{i}=\sum\nolimits_{k=1}^{K}\frac{\eta_{k}^{i}}{\max\limits_{1\leqslant
i\leqslant I}g_{k}(\mathbf{\Theta}_{i})-\min\limits_{1\leqslant i\leqslant
I}g_{k}(\mathbf{\Theta}_{i})},$ (25)

where $\xi_{rcd}^{i}$ denotes the relative crowding distance of the $i$th
solution. Then, we selected out the ones whose CDVs are non-dominated as
candidates of global best. Afterwards, within all these candidates, we choose
$Y$ global best solutions whose relative crowding distance are bigger than the
rest from these candidates. Finally, we divide the swarm of particles into $Y$
groups. We allocate one global best solution for each divided group of
particles in order to update their positions as well as their velocity. The
number of particles contains in each group will be fixed as:

$N_{y}=S\times\frac{\xi_{rcd}^{y}}{\sum\nolimits_{y=1}^{Y}\xi_{rcd}^{y}},$
(26)

where $\xi_{rcd}^{y}$ is the relative crowding distance of the $y$th global
best solution, $S$ is the number of particles in swarm and $N_{y}$ is the
number of particles in the $y$th group.

We still use the example in Fig.3 to illuminate this selection process.
Firstly, the CDVs (i.e., $\mathbf{H}$), the crowding distance (i.e.,
$\xi_{cd}$) and the relative crowding distance (i.e., $\xi_{rcd}$) of all the
6 solutions are given out in TABLE II. Then, with non-dominated CDV, the
$\mathbf{\Theta_{2}}$, $\mathbf{\Theta_{3}}$ and $\mathbf{\Theta_{4}}$ are
qualified to be the candidates. Thus, the global best solution will be
selected out according to their $\xi_{rcd}$ (i.e., the selected ones will be
marked with circle). For instance, since with the largest $\xi_{rcd}$, the
$\mathbf{\Theta_{4}}$ has the highest priority to be selected as global best
in any case, while the $\mathbf{\Theta_{2}}$ is with the lowest priority (only
when $Y=3$) for its smallest $\xi_{rcd}$. Finally, we assume there are $100$
particles in total and we give out the result of particle dividing when
$Y=1,2,3$ (i.e., the numbers of particles in the group assigned to each
selected solution are given out after each circle). We notice that the
$\mathbf{\Theta_{5}}$ is not be selected for any case though its relative
crowding distance higher than that of the $\mathbf{\Theta_{2}}$, it is because
that the CDV of the $\mathbf{\Theta_{5}}$ is dominated by the
$\mathbf{\Theta_{3}}$, which excluding its possibility to be selected as
global best. We also give out the pseudo code of this selection method in
Algorithm 1.

TABLE II: Absolute Crowding Distance, Crowding Distance Vector and Relative Crowding Distance of Solutions in Fig.3 |  $\xi^{i}_{cd}$ |  $\mathbf{H}^{i}$ |  $\xi^{i}_{rcd}$ |  $Y=1$ |  $Y=2$ |  $Y=3$
---|---|---|---|---|---|---
$\mathbf{G(\Theta_{1})}$ |  41.2 |  $[0.18,40]^{\text{T}}$ |  0.58 |  |  |
$\mathbf{G(\Theta_{2})}$ |  41.2 |  $\mathbf{[0.2,41]^{\text{T}}}$ |  0.61 |  |  |  $\bigcirc$(23)
$\mathbf{G(\Theta_{3})}$ |  41.51 |  $\mathbf{[0.51,40]^{\text{T}}}$ |  0.91 |  |  $\bigcirc$(45) |  $\bigcirc$(35)
$\mathbf{G(\Theta_{4})}$ |  39.71 |  $\mathbf{[0.71,39]^{\text{T}}}$ |  1.1 |  $\bigcirc$(100) |  $\bigcirc$(55) |  $\bigcirc$(42)
$\mathbf{G(\Theta_{5})}$ |  40.4 |  $[0.4,40]^{\text{T}}$ |  0.8 |  |  |
$\mathbf{G(\Theta_{6})}$ |  41.2 |  $[0.18,40]^{\text{T}}$ |  0.58 |  |  |

  * •

$Y_{\text{actual}}$ is the number of solutions whose $\mathbf{H}$s are non-
dominated after check all the solutions that stored in external archive

  * •

$Y_{\text{u}}$ is the fixed maximum number of global best solutions

  * •

$Y$ is the number of global best solutions

Calculation the $\mathbf{H}$ and $\xi_{ncd}$ for every solution that stored in
external archive;

Select out the solutions whose $\mathbf{H}$s are non-dominated;

if  _ $Y_{\text{actual}}\leqslant Y_{\text{u}}$ _ then

Choose all the $Y_{\text{actual}}$ solutions as global best solutions;

$Y=Y_{\text{actual}}$;

for  _ $y=1,2,\ldots,Y$ _ do

Use (26) to calculate $N_{y}$ for the $y$th selected global best solution;

end for

end if

if  _ $Y_{\text{actual}}>Y_{\text{u}}$ _ then

Sort the $Y_{\text{actual}}$ solutions according their relative crowding
distance (i.e., $\xi_{rcd}$) in descending order, pick the first $Y_{u}$ ones
(whose relative crowding distance are higher than that of the rest) out as
global best solutions;

$Y=Y_{\text{u}}$;

for  _ $y=1,2,\ldots,Y$ _ do

Use (26) to calculate $N_{y}$ for the $y$th selected global best solution;

end for

end if

Algorithm 1 Pseudo code of selection method based on non-dominated relative
crowding distance

For the challenge of particle allocation: we reference to a classic type of
MOPSO with multi-swarm structure [30]. $K+1$ swarms of particles are employed
for the searching and different swarms has different missions. We use $K$ of
these swarms, which are called as sub-swarms, to enlarge the shape of
$\mathcal{P_{F}}$. It is well known that the extremum of each objective
function constitutes the edge of the obtained $\mathcal{P_{F}}$. Then, we let
each sub-swarm do a single-objective PSO by employing one of the objectives as
fitness function. By assigning each of the concerned objective a swarm of
particles to find a better objective function value, the shape of
$\mathcal{P_{F}}$ could be enlarged effectively. For the rest one swarm, which
is called as main swarm, we assign it to find non-dominated solutions within
the edge of the obtained $\mathcal{P_{F}}$. Thus, with the particle
cooperation based on labour dividing, the search process could be of better
pertinence and we can get better optimization performance.

Figure 4: Flow diagram of MOPSO-NRCD.

With using a global best selection method, which is based on a non-dominated
relationship of CDV and a relative crowding distance, we call the proposed
algorithm as MOPSO with non-dominated relative crowding distance (MOPSO-NRCD)
and we give its flow diagram in Fig.4. We can see that the main swarm and all
the sub-swarms work in a parallel way. In every iteration, the main swarm is
dynamically divided into $Y$ groups to chase the $Y$ selected global best
solutions in order to find more non-dominated solutions. Meanwhile, each sub-
swarm aims to find better values for its assigned objective function and all
sub-swarms work together in order to get a bigger $\mathcal{P_{F}}$. The
pseudo code of MOPSO-NRCD is given at Algorithm 2 (the calculation method of
objective functions is given in Algorithm 3, whose detail will be given in the
next section).

  * •

$K$ is the number of the objective functions and the number of sub-swarms

  * •

$S_{m}$ is the number of particles in main swarm

  * •

$S_{k}$ is the number of particles in the $k$th sub-swarm

for  _ $s=1,2,\ldots,S_{m}$ _ do

Initialize particle and its velocity;

Use Algorithm 3 to calculate objective functions of particle;

Set the initial individual best solution of particle as its initial position;

end for

for  _ $k=1,2,\ldots,K$ _ do

for  _ $s=1,2,\dots,S_{k}$ _ do

Initialize particle and its velocity;

Use Algorithm 3 to calculate objective functions of particle;

Set the initial individual best solution of particle as its initial position;

end for

Use $g_{k}$ as fitness function to Initialize the global best solution of the
$k$th sub-swarm;

end for

Select non-dominated ones from all DV that proposed by particles in main swarm
and sub-swarms;

Store the non-dominated DVs in external archive;

for  _ $t=1,2,\ldots,T_{max}$_ do

Use Algorithm 1 to select $Y$ global bests solutions and divide the main swarm
into $Y$ groups (the $y$th group contains $N_{y}$ particles);

for  _ $y=1,2,\ldots,Y$ _ do

utilize the $y$th global best solution to update the $N_{y}$ particles in
$y$th group;

Use Algorithm 3 to calculate the objective functions of particle;

Update the individual best solution of each particle in $y$th group;

end for

for  _ $k=1,2,\ldots,K$ _ do

for  _ $s=1,2,\dots,S_{k}$ _ do

Update the velocity and position of particle;

Use Algorithm 3 to calculate the objective functions of particle;

Use $g_{k}$ as fitness function to update the individual best solution of
particle;

end for

Use $g_{k}$ as fitness function to update the global best solution of the
$k$th sub-swarm;

end for

Use the obtained non-dominated DVs in current iteration to update external
archive;

end for

Output the solutions that stored in external archive.

Algorithm 2 Pseudo code of MOPSO-NRCD

##  IV Numerical Results

TABLE III: Parameter Setting of Simulation Parameter Name | Value | Parameter Name | Value
---|---|---|---
$P_{dt}$ | $0.8$ | $P_{fa}$ | $10^{-6}$
$D_{0}$ | 12.5dB | $\sigma$ | $1$
$T_{max}$ | $2000$ | $R_{max}$ | $6\hskip 1.42262pt\text{km}$
$Y_{u}$ | 3 | Size of grid | $2.5\hskip 1.42262pt\text{km}^{2}$
$c_{1}$ | $2$ | $c_{2}$ | $2$
$V_{max}$ | $4$ |  $S$ of MOPSO-CD | $200$
$S_{m}$ of MOPSO-NRCD | $100$ |  $S_{k}$ of MOPSO-NRCD | $50$

Before the introduction of simulation results, we firstly give out some common
simulation settings in TABLE III. For fairly comparison, the number of the
employed particles of both MOPSO-CD and MOPSO-NRCD are set as the same (i.e.,
the number is set as $200$). According to the designed performance and the
assumption that every radar transmits one pulse during an observation time, we
can find out $D_{0}=12.5$ dB [27]. We also set that
$\sigma=1,\sigma_{m,n}=|\alpha|^{2}$, where $\alpha$ is assumed to be zero-
mean complex gaussian. The $\omega(t)$ is set to be decreased with iteration
time as $\omega(t)=0.9-0.5*(t/T_{\max})$. Let placement region and
surveillance region are coincident ($\mathbb{A}=\mathbb{B}=50\hskip
1.42262pt\text{km}\times 50\hskip 1.42262pt\text{km}$). The calculation method
of the objective functions of (15) is given in Algorithm 3.

  * •

$R$ is the number of resolution cells in surveillance region

  * •

$J$ is the number of radar nodes in MSRS

for  _ $r=1,2,\ldots,R$ _ do

Calculate the transmitter-to-target ranges (i.e., $R_{tm}$s, $m=1,2,\dots,J$)
and the receiver-to-target ranges (i.e., $R_{rn}$s, $n=1,2,\dots,J$) for each
radar nodes;

Calculate the ratio of signal energy to noise power for every T-R pair by
substituting the values of $D_{0}$, $\sigma$. $\sigma_{m,n}$, and $R_{max}$
into (10);

Calculate the $\chi_{r}$ for one resolution cell with (II-C2) and (13) for
cooperative and non-cooperative mode, respectively;

Calculate the $P_{d}$ by submitting the obtained $\chi_{r}$ and the fixed
$P_{fa}$ into (4) and (5) for cooperative and non-cooperative mode,
respectively;

end for

Calculate the $\mathcal{C_{R}}$ and $\mathcal{L_{R}}$ by (6) and (14),
respectively;

Algorithm 3 Pseudo code of calculation methods of $\mathcal{C_{R}}$ and
$\mathcal{L_{R}}$ in (15)

(a)

(b)

(c)

(d)

Figure 5: The performance of solutions obtained by random deployment, MOPSO-CD
and MOPSO-NRCD under condition of cooperative mode. (a). 5 antennas; (b). 6
antennas; (c). 7 antennas; (d). 8 antennas.

(a)

(b)

(c)

(d)

Figure 6: The performance of solutions obtained by random deployment, MOPSO-CD
and MOPSO-NRCD under condition of non-cooperative mode. (a). 5 antennas; (b).
6 antennas; (c). 7 antennas; (d). 8 antennas.

We execute 5 independent optimization tests for both MOPSO-CD and NOPSO-NRCD
in different simulation cases, which are characterized by different antenna
number and different working mode. Meanwhile, we also give the performance of
50 solutions obtained by random deployment in each case for comparison. From
Fig. 5 and Fig. 6, we can see that the MOPSO-CD can obtain better solutions
comparing with that of random deployment. Specifically, all the solutions of
random deployment are dominated by at least one solution obtained by the
MOPSO-CD (i.e., the $\mathcal{C_{R}}$ and the $\mathcal{L_{R}}$ of each
solution obtained by random deployment are both smaller than that of at least
one solutions obtained by MOPSO-CD). Meanwhile, we can see that each solution
obtained by MOPSO-CD are dominated by at least one solution of MOPSO-NRCD in
any cases. This proves that the MOPSO-NRCD is capable to achieve better
optimization performance than MOPSO-CD.

It is hard to directly quantify the improvement on each concerned objective
function (e.g., MOPSO-CD compares with random deployment and MOPSO-NRCD
compares with MOPSO-CD) from Fig. 5 and Fig. 6. Therefore, we propose a metric
to solve this problem. The introduction of the metric begins with some
definition. For comparison, we call the group of solutions with better
performance as improved group of solutions (IGS), while the group of solutions
with worse performance as control group of solutions (CGS). Assume that there
are $I^{z}$ solutions in IGS that dominate the $z$th solution in CGS (assume
that CGS contains $Z$ non-dominated solutions). Then, we can give the
calculation method of the proposed metric as:

$\mathcal{A_{I}}^{k}=\frac{1}{Z}\sum\nolimits_{z=1}^{Z}(\frac{1}{I^{z}}\sum\nolimits_{i^{\prime}=1}^{I^{z}}(g_{k}(\mathbf{\Theta}_{{i^{\prime}}})-g_{k}(\mathbf{\Theta}_{z}))).$
(27)

Obviously, the $\mathcal{A_{I}}^{k}$ represents the average improvement in
$k$th objective function of IGS comparing with CGS without worsening in other
objective functions. Using the obtained solutions in Fig. 5 and Fig. 6, we
analyse the $\mathcal{A_{I}}^{k}$s of MOPSO-CD comparing with random
deployment (i.e., the solutions of MOPSO-CD constitute the IGS and the
solutions of random deployment constitute the CGS) as well as that of MOPSO-
NRCD comparing with MOPSO-CD (i.e., the solutions of MOPSO-NRCD constitute the
IGS and the solutions of MOPSO-CD constitute the CGS). The analysing results
are shown in TABLE IV and TABLE V. We can see that the optimal solutions of
MOPSO-CD can bring significant improvement when comparing with the solutions
of random deployment, and thus the validity of MOPSO is proved. Besides, since
the solutions of MOPSO-NRCD bringing further improvement in both of these two
concerned objectives when comparing with that of MOPSO-CD, the advantages of
MOPSO-NRCD is also proved.

TABLE IV: $\mathcal{A_{I}}^{k}$ of MOPSO-CD compares with random deployment And MOPSO-NRCD compares with MOPSO-CD in Fig.5 $J$ | MOPSO-CD compares | MOPSO-NRCD compares
---|---|---
| with random deployment | with MOPSO-CD
|  $\mathcal{L_{R}}$ (dB) |  $\mathcal{C_{R}}$ (%) |  $\mathcal{L_{R}}$ (dB) |  $\mathcal{C_{R}}$ (%)
5 |  5.47 |  0.0616 |  0.71 |  0.0097
6 |  6.09 |  0.0698 |  0.69 |  0.0137
7 |  5.97 |  0.0849 |  0.70 |  0.0144
8 |  5.74 |  0.0922 |  0.81 |  0.0174
TABLE V: $\mathcal{A_{I}}^{k}$ of MOPSO-CD compares with random deployment And MOPSO-NRCD compares with MOPSO-CD in Fig.6 $J$ | MOPSO-CD compares | MOPSO-NRCD compares
---|---|---
| with random deployment | with MOPSO-CD
|  $\mathcal{L_{R}}$ (dB) |  $\mathcal{C_{R}}$ |  $\mathcal{L_{R}}$ (dB) |  $\mathcal{C_{R}}$
5 |  11.61 |  0.0722 |  1.03 |  0.0115
6 |  11.61 |  0.0823 |  1.07 |  0.0144
7 |  11.06 |  0.0948 |  0.86 |  0.0151
8 |  12.52 |  0.1083 |  0.59 |  0.0129

From another aspect, instead of analysing the concerned objective one by one,
we employ the size of dominated space to jointly evaluate the optimization
performance of MOPSO-CD and MOPSO-NRCD. The size of dominated space has been
introduced by Zitzler in [39]. It is a measure of the amount of objective
space weakly dominated by a given $\mathcal{P_{F}}$ which is constituted by a
set of solutions. Generally, a bigger size of dominated space means a better
optimization performance. To calculate the size of dominated space, we set the
reference point as $\mathcal{L_{R}}=-15\hskip 1.42262pt\text{dB}$ and
$\mathcal{C_{R}}=0.15$. In TABLE VI, we can see that the average size of
dominated space of MOPSO-NRCD is larger than that of MOPSO-CD in any cases,
which improves the advantages of MOPSO-NRCD on the other hand.

TABLE VI: Average Size of Dominated Space of MOPSO-CD and MOPSO-NRCD in Cooperative Mode and Non-cooperative Mode of 5 Independent Optimization Tests $J$ | Cooperative Mode | Non-cooperative Mode
---|---|---
MOPSO-CD | MOPSO-NRCD | MOPSO-CD | MOPSO-NRCD
5 | 2.25 | 2.62 | 1.21 | 1.4
6 | 3.51 | 4.13 | 1.93 | 2.23
7 | 5 | 5.74 | 2.73 | 3.25
8 | 6.51 | 7.78 | 3.65 | 4.31

##  V Conclusion

An optimization deployment problem of MSRS radar was studied in this paper.
The optimization goal is to improve the surveillance performance of a
surveillance region by optimally placing the antennas of MSRS and allocating
their transmitted power. Firstly, by introducing two objectives, we had
formulated the problem under consideration of two typical working mode of
MSRS. The corresponding objective functions as well as their calculation
methods were also proposed as well. Then, focusing on the challenge of value
difference and the challenge of particle allocation, we had proposed a novel
MOPSO with non-dominated relative crowding distance (MOPSO-NRCD). Simulation
results showed that the proposed algorithm can achieve the better performance
than the MOPSO-CD with same number of employed particles for different working
mode and different antenna number of MSRS.

## References

  * [1] H.D. Griffiths, “Multistatic, MIMO and networked radar: the future of radar sensors?”, invited keynote presentation, EuRAD conference, Paris, 30 September 2010.
  * [2] H.D. Griffiths and C. Baker, “Towards the intelligent adaptive radar network,” in Radar Conference (RADAR), 2013 IEEE, April 2013.
  * [3] A. M. Haimovich, R. S. Blum, and L. J. Cimini, “MIMO radar with widely separated antennas.” IEEE Signal Process. Mag., vol. 25, no. 1, pp. 116-129, Jan. 2008.
  * [4] Q. He, N. H. Lehmann, R. S. Blum, and A. M. Haimovich, “MIMO radar moving target detection in homogeneous clutter,” IEEE Trans. Aerosp. Electron. Syst., vol. 46, no. 3, pp. 1290-1301, Jul. 2010.
  * [5] P. Wang, H. Li, and B. Himed, “Moving target detection using distributed MIMO radar in clutter with nonhomogeneous power.”, IEEE Trans. Signal Process., vol. 59, no. 10, pp. 4809-4820, Oct. 2011.
  * [6] A. Maio and M. Lops, “Design principles of MIMO radar detectors, IEEE Trans. Aerosp. Electron. Syst., vol. 43, pp. 886-898, Jul. 2007.
  * [7] H. Godrich, A. M. Haimovich, and R. S. Blum, “Target localization accuracy gain in MIMO radar-based systems,” IEEE Trans. Inf. Theory, vol. 56, pp. 2783-2803, Jun. 2010.
  * [8] L. Neng-Jing, “Radar ECCMs new area: Anti-stealth and anti-ARM,” IEEE Trans. Aerosp. Electron. Syst., vol. 31, no. 3, pp. 1120-1127, Jul. 1995.
  * [9] Q. He, R. S. Blum, H. Godrich, and A. M. Haimovich, “Target velocity estimation and antenna placement for MIMO radar with widely separated antennas,” IEEE J. Sel. Topics Signal Process., vol. 4, no. 1, pp.79-100, 2010.
  * [10] H. Godrich, A. Petropulu, and H. V. Poor, “Power allocation strategies for target localization in distributed multiple-radar architectures,” IEEE Trans. Signal Process., vol. 59, no. 7, pp. 3226-3240, Jul. 2011.
  * [11] W. Chen and R. M. Narayanan, “Antenna placement for minimizing target localization error in UWB MIMO noise radar,” IEEE Antennas Wireless Propag. Lett., vol. 10, pp. 135-138, 2011.
  * [12] P. Chavali and A. Nehorai, “Scheduling and power allocation in a cognitive radar network for multiple-target tracking,” IEEE Trans. Signal Process., vol. 60, no. 2, pp. 715-729, Feb. 2012.
  * [13] M. Radmard, M. M. Chitgarha, M. Nazari-Majd, and M. M. Nayebi, “Antenna placement and power allocation optimization in MIMO detection.” IEEE Trans. Aerosp. Electron. Syst., vol. 50, no. 2, pp. 1468-1478, 2014.
  * [14] Y. Yang, W. Yi, T. Zhang, G. Cui, L. Kong, X. Yang, J. Yang, “Fast Optimal Antenna Placement for Distributed MIMO Radar with Surveillance Performance,” IEEE Signal Process. Lett., vol.22, pp.1955-1959, Nov. 2015
  * [15] L. Ke, Q. Zhang, and R. Battiti, “MOEA/D-ACO: A multiobjective evolutionary algorithm using decomposition and ant colony, IEEE Trans. Cybern., vol. 43, no. 6, pp. 1845-1859, Dec. 2013.
  * [16] S. W. Jiang, J. Zhang, Y. S. Ong, A. N. S. Zhang, and P. S. Tan, “A Simple and Fast Hypervolume Indicator-based Multiobjective Evolutionary Algorithm , IEEE Trans. Cybern., vol. 45, no. 10, pp. 2202-2213, 2014.
  * [17] V. A. Shim, K. C. Tan, and H. J. Tang, “Adaptive memetic computing for evolutionary multi-objective optimization , IEEE Trans. Cybern., vol. 45, no. 4, pp. 610-621, 2014.
  * [18] M. I. Skolnik, Radar Handbook, 2nd ed. New York: McGraw-Hill, Jan. 1990.
  * [19] D. A. Shnidman, “Radar detection probabilities and their calculation,” IEEE Trans. Aerosp. Electron. Syst., vol. 31, no. 3, pp. 928-950, July 1995.
  * [20] G. Cui, A. De Maio, and Piezzo, M. “Performance prediction of the incoherent radar detector for correlated generalized Swerling-chi fluctuating targets.” IEEE Trans. Aerosp. Electron. Syst., vol. 49, no. 2, pp. 356-368, 2013.
  * [21] N. Aziz, A. Mohemmed, and M. Alias, “A wireless sensor network coverage optimization algorithm based on particle swarm optimization and Voronoi diagram,” in Proc. Int. Conf. Networking, Sensing and Control (ICNSC), Mar. 2009, pp. 602-607.
  * [22] X. Wang and S. Wang, “Hierarchical Deployment Optimization for Wireless Sensor Networks,” IEEE Trans. Mobile Computing, vol. 10, on. 7, pp. 1028-1041, July. 2011.
  * [23] R. V. Kulkarni and G. K. Venayagamoorthy, Particle swarm optimization in wireless-sensor networks: A brief survey,” IEEE Trans. Syst., Man, Cybern. C, Appl. Rev., vol. 41, no. 2, pp. 262-267, 2011.
  * [24] S. Kiranyaz, T. Ince, A. Yildirim, and M. Gabbouj, Fractional particle swarm optimization in multidimensional search space, IEEE Trans. Syst., Man, Cybern., Part B, vol. 40, no. 2, pp. 298-319, Apr. 2010.
  * [25] Z. H. Zhan, J. Zhang, Y. Li, and H. S. H. Chung, “Adaptive particle swarm optimization, IEEE Trans. Syst., Man, Cybern. B, Cybern., vol. 39, no. 6, pp. 1362-1381, Dec. 2009.
  * [26] J. Kennedy and R. Eberhart, “Particle swarm optimization,” in Proc.IEEE Int. Conf. Neural Netw., 27 Nov. C1 Dec., 1995, vol. 4, pp. 1942-1948
  * [27] M. A. Richards, Fundementals of Radar Signal Processing. New York: McGraw-Hill, 2005.
  * [28] Chernyak. V, “Multisite Radar Systems Composed of MIMO Radars,” IEEE AES Syst. Mag., vol. 29, pp. 28-37, Dec. 2014.
  * [29] C. A. Coello Coello and M. S. Lechuga, “MOPSO: A proposal for multiple objective particle swarm optimization,” in Proc. Congr. Evol. Comput., 2002, pp. 1051-1056.
  * [30] K. E. Parsopoulos and M. N. Vrahatis, “Particle swarm optimization method in multiobjective problems,” in Proc. ACM Symp. Applied Computing 2002 (SAC 2002), 2002, pp. 603-607.
  * [31] C. A. Coello Coello, G. T. Pulido, and M. S. Lechuga, “Handling multiple objectives with particle swarm optimization,” IEEE Trans. Evol. Comput., vol. 8, no. 3, pp. 256-279, Jun. 2004.
  * [32] C. R. Raquel and P. C. Nava, “An effective use of crowding distance in multiobjective particle swarm optimization,” in Proc. Genetic Evol. Comput., 2005, pp. 257-264.
  * [33] M. Reyes-Sierra and C. A. Coello Coello, “Multi-objective particle swarm optimizers: A survey of the state-of-the-art,” Int. J. Comput. Intell. Res., vol. 2, no. 3, pp. 287-308, 2006.
  * [34] W. Hu and G. G. Yen, “Adaptive multiobjective particle swarm optimization based on parallel cell coordinate system, IEEE Trans. Evol. Comput., vol. 19, no. 1, pp. 1 C18, Feb. 2015.
  * [35] S. K. Goudos, Z. D. Zaharis, D. G. Kampitaki, I. T. Rekanos, and C. S. Hilas, “Pareto optimal design of dual-band base station antenna arrays using multiobjective particle swarm optimization with fitness sharing,” IEEE Trans. Magnetics., vol. 45, no. 3, pp. 1522-1525, Mar. 2009.
  * [36] S. Chamaani, S. A. Mirtaheri, and M. S. Abrishamian, “Improvement of time frequency-domain performance of antipodal Vivaldi antenna using multiobjective particle swarm optimization,” IEEE Trans. Antennas Propag., vol. 59, no. 5, pp. 1738-1742, May 2011.
  * [37] B. Xue, M. Zhang, W. Browne, “Particle swarm optimization for feature selec-tion in classification: a multi-objective approach,” IEEE Trans. Cybern., vol. 43, no. 6, pp. 1656-1671, Dec. 2012.
  * [38] P. M. Pradhan and G. Panda, “Connectivity constrained wireless sensor deployment using multiobjective evolutionary algorithms and fuzzy decision making,” Ad Hoc Netw., vol. 10, no. 6, pp. 1134-1145, 2012.
  * [39] E. Zitzler, “Evolutionary algorithms for multiobjective optimization: Methods and applications,” Doctoral dissertation ETH 13398, Swiss Federal Instit. Technol. (ETH), Zurich, Switzerland, 1999.
  * [40] S. P. Sira, D. Cochran, A. Papandreou-Suppappola, D. Morrell, W. Moran, S. D. Howard, and R. Calderbank, “Adaptive waveform design for improved detection of low-RCS targets in heavy sea clutter,” IEEE J. Sel. Topics Signal Process., vol. 1, no. 1, pp. 56-66, Jun. 2007.
  * [41] S. M. Tonissen and R. J. Evans, “Peformance of dynamic programming techniques for track-before-detect, IEEE Trans. Aerosp. Electron. Syst., vol. 32, no. 4, pp. 1440-1451, Oct. 1996.
  * [42] W. Yi, M. Morelande, L. Kong, and J. Yang. “A computationally efficient particle filter for multitarget tracking using an independence approximation.” IEEE Trans. on Signal Processing., 61(4):843-856, 2013.
  * [43] R. Calderbank, S. Howard, and B. Moran, “Waveform diversity in radar signal processing,” IEEE Signal Process. Mag., vol. 26, no. 1, pp. 32-41, Feb. 2009.
  * [44] T. Zhang, G. Cui, L. Kong, W. Yi, X. Yang, “Phase-Modulated Waveform Evaluation and Selection Strategy in Compound-Gaussian Clutter,” IEEE Trans. Signal Process., vol.61, no.5, pp.1143-1148, Mar, 2013


11institutetext:  IT University of Copenhagen
11email<EMAIL_ADDRESS>
http://thomasahle.com

# On the Problem of $p_{1}^{-1}$ in Locality-Sensitive Hashing

Thomas Dybdahl Ahle  0000-0001-9747-0479

###### Abstract

A Locality-Sensitive Hash (LSH) function is called
$(r,cr,p_{1},p_{2})$-sensitive, if two data-points with a distance less than
$r$ collide with probability at least $p_{1}$ while data points with a
distance greater than $cr$ collide with probability at most $p_{2}$. These
functions form the basis of the successful Indyk-Motwani algorithm (STOC 1998)
for nearest neighbour problems. In particular one may build a $c$-approximate
nearest neighbour data structure with query time $\tilde{O}(n^{\rho}/p_{1})$
where $\rho=\frac{\log 1/p_{1}}{\log 1/p_{2}}\in(0,1)$. That is, _sub-linear
time_ , as long as $p_{1}$ is not too small. This is significant since most
high dimensional nearest neighbour problems suffer from the curse of
dimensionality, and can’t be solved _exact_ , faster than a brute force
_linear-time_ scan of the database.

Unfortunately, the best LSH functions tend to have very low collision
probabilities, $p_{1}$ and $p_{2}$. Including the best functions for Cosine
and Jaccard Similarity. This means that the $n^{\rho}/p_{1}$ query time of
_LSH is often not sub-linear after all_ , even for approximate nearest
neighbours!

In this paper, we improve the general Indyk-Motwani algorithm to reduce the
query time of LSH to $\tilde{O}(n^{\rho}/p_{1}^{1-\rho})$ (and the space usage
correspondingly.) Since $n^{\rho}p_{1}^{\rho-1}<n\Leftrightarrow
p_{1}>n^{-1}$, our algorithm always obtains sublinear query time, for any
collision probabilities at least $1/n$. For $p_{1}$ and $p_{2}$ small enough,
our improvement over all previous methods can be _up to a factor $n$_ in both
query time and space.

The improvement comes from a simple change to the Indyk-Motwani algorithm,
which can easily be implemented in existing software packages.

###### Keywords:

locality-sensitive hashing nearest neighbor similarity search

##  1 Introduction

Locality Sensitive-Hashing (LSH) framework [18] is one of the most efficient
approaches to the nearest neighbour search problem in high dimensional spaces.
It comes with theoretical guarantees, and it has the advantage of easy
adaption to nearly any metric or similarity function one might want to search.

The $(r_{1},r_{2})$-near neighbour problem is defined as follows: Given a set
$X$ of points, we build a data-structure, such that given a query, $q$ we can
quickly find a point $x\in X$ with distance $<r_{2}$ to $q$, or determine that
$X$ has no points with distance $\leq r_{1}$ to $q$. Given a solution to this
“gap” problem, one can obtain a $r_{1}/r_{2}$-approximate nearest neighbour
data structure, or even an exact111In general we expect the exact problem to
be impossible to solve in sub-linear time, given the hardness results of [5,
1]. However for practical datasets it is often possible. solution using known
reductions [2, 14, 17].

For any measure of similarity, the gap problem can be solved by LSH: we find a
distribution of functions $H$, such that $p_{1}\geq\Pr_{h\sim H}[h(x)=h(y)]$
when $x$ and $y$ are similar, and $p_{2}\leq\Pr_{h\sim H}[h(x)=h(y)]$ when $x$
and $y$ are dissimilar. Such a distribution is called
$(r_{1},r_{2},p_{1},p_{2})$-sensitive. If $p_{1}>p_{2}$ the LSH framework
gives a data-structure with query time $\tilde{O}(n^{\rho}/p_{1})$ for
$\rho=\frac{\log 1/p_{1}}{\log 1/p_{2}}$, which is usually significantly
faster than the alternatives.

_At least when $p_{1}$ is not too small._

The two most common families of LSH is Cross-Polytope (or Spherical) LSH [6]
for Cosine similarity and MinHash [11, 10] for Jaccard Similarity.

Cross-Polytope is the basis of the Falconn software package [21], and solves
the $(r,cr)$-near neighbour problem on the sphere in time
$\tilde{O}(n^{1/c^{2}}/p_{1})$. Here
$p_{1}=d^{-\frac{\tau^{2}}{4-\tau^{2}}}(\log d)^{-\Omega(1)}$, where
$\tau=\|p-q\|_{2}\in[0,2]$ is the distance between two close points. We see
that already at $\tau\approx\sqrt{2}$ (which corresponds to near orthogonal
vectors) the $1/p_{1}$ factor results in a factor $d$ slow-down. For larger
$\tau\in(\sqrt{2},2]$ the slow-down can grow arbitrary large. Using
dimensionality reduction techniques, like the Johnson Lindenstrauss transform,
one may assume $d=\varepsilon^{-2}\log n$ at the cost of a factor
$1+\varepsilon$ distortion of the distances. However if $\varepsilon$ is just
$1/100$, the slow-down factor of $d$ is still worse than, say, $n^{1/2}$ for
datasets of size up to $10^{8}$, and so if $c\leq\sqrt{2}$ we get that
$n^{\rho}/p_{1}$ s larger than $n$. So _worse than a brute force scan of the
database_!

The MinHash algorithm was introduced by Broder et al. for the Alta Vista
search engine, but is used today for similarity search on sets in everything
from natural language processing to gene sequencing. MinHash solves the
$(j_{1},j_{2})$ gap similarity search problem, where $j_{1}\in(0,1)$ is the
Jaccard Similarity of similar sets, and $j_{2}$ is that of dissimilar sets, in
time $\tilde{O}(n^{\rho}/j_{1})$ where $\rho=\frac{\log 1/j_{1}}{\log
1/j_{2}}$. (In particular MinHash is $(j_{1},j_{2},j_{1},j_{2})$-sensitive in
the sense defined above.) Now consider the case $j_{1}=n^{-1/4}$ and
$j_{2}=n^{-3/10}$. This is fairly common as illustrated in fig. 1(a). In this
case $\rho=\frac{\log 1/j_{1}}{\log 1/j_{2}}=5/6$, so we end up with
$n^{\rho}/j_{1}=n^{13/12}$. Again _worse than a brute force scan of the
database_!

In this paper we reduce the query time of LSH to $n^{\rho}/p_{1}^{1-\rho}$,
which is less than $n$ for all $p_{1}>1/n$. In the MinHash example above, we
get $n^{\rho}/p_{1}^{1-\rho}=n^{5/6+1/4(1-5/6)}=n^{7/8}$. More than a factor
$n^{.208}$ improvement(!) In general the improvement of $p_{1}^{-\rho}$ may be
as large as a factor of $n$ when $p_{1}$ and $p_{2}$ are both close to $1/n$.
This is illustrated in fig. 1(b).

The improvements to LSH comes from a simple observation: During the algorithm
of Indyk and a certain “amplification” procedure has to be applied $\frac{\log
n}{\log 1/p_{2}}$ times. When $\log 1/p_{2}$ does not divide $n$, which is
extremely likely, the amount of amplification has to be approximated by the
next integer. We propose instead an ensemble of LSH tables with different
amounts of amplification, which when analysed sufficiently precisely yields
the improvements described above.

(a) Density plots of pairwise Jaccard Similarities in the datasets studied by
Mann et al. [20]. The similarities are normalized by the dataset sizes, so we
can compare the effect of $1/p_{1}$ with the effect of $n^{\rho}$. We see that
reasonable values for $j_{1}=p_{1}$ range between $n^{-1/3}$ and $n^{-1/6}$ on
those datasets.

(b) Saving possible, in query time and space, over classical LSH as a function
of $p_{1}$ and $p_{2}$. With $p_{1}=n^{-1/4}$ and $p_{2}=n^{-1/3}$ we save a
factor of $n^{3/16}=n^{.1875}$.

Figure 1: Overview over available savings

###  1.1 Related Work

We will review various directions in which LSH has been improved and
generalized, and how those results related to what is presented in the present
article.

In many cases, the time required to sample and evaluate the hash functions
dominate the time required by LSH. Recent papers [16, 12] have reduced the
number of distinct calls to the hash functions which is needed. The most
recent paper in the line of work is [12], which reduces the number of calls to
$(\frac{\log n}{\log 1/p_{2}})^{2}/p_{1}$. On top of that, however, they still
require $n^{\rho}/p_{1}$ work, so the issue with small $p_{1}$ isn’t touched
upon. In fact, some of the some of the algorithms in [12] _increase_ the
dependency from $n^{\rho}/p_{1}$ to $n^{\rho}/(p_{1}-p_{2})$.

Other work has sought to generalize the concept of Locality Sensitive Hashing
to so-called Locality Sensitive Filtering, LSF [9]. However, the best work for
set similarity search based on LSF [4, 13] still have factors similar in
spirit to $p_{1}^{-1}$. E.g., the Chosen Path algorithm in [13] uses query
time $\tilde{O}(n^{\rho}/b_{1})$, where $b_{1}$ is the similarity between
close sets.

A third line of work has sought to derandomize LSH. The result is so-called
Las Vegas LSH [3, 22]. Here the families $H$ are built combinatorially, rather
than at random, to guarantee the data structure always return a near
neighbour, when one exists. While these methods don’t have probabilities, they
still end up with similar factors for similar reasons.

As mentioned, the reason $p_{1}^{-1}$ shows up in all these different
approaches, is that they all rely on the same amplification procedure, which
has to be applied an integer number of times. One might wonder if tree based
methods, which do an adaptive amount of amplification, could get rid of the
$1/p_{1}$ dependency. However as evidenced by the classical and current work
[8, 7, 14, 15] these methods still have a factor $1/p_{1}$. We leave it open
whether this might be avoidable with better analysis, perhaps inspired by the
results in this paper.

##  2 Preliminaries

Before we give the new LSH algorithm, we will recap the traditional analysis.
For a more comprehensive introduction to LSH, see the Mining of Massive
Datasets book [19], Chapter 3. In the remainder of the article we will use the
notation $[n]=\\{1,\dots,n\\}$.

Assume we are given a $(r_{1},r_{2},p_{1},p_{2})$-sensitive LSH family, $H$,
as defined in the introduction. Let $k$ and $L$ be some integers defined
later, and let $[m]$ be the range of the hash functions, $h\in H$. Let $n$ be
an upper bound on the number of points to be inserted. 222If we don’t know how
many points will be inserted, several black box reductions allow transforming
LSH into a dynamic data structure. The Indyk-Motwani data-structure consists
of $L$ hash tables, each with $m^{k}$ hash buckets.

To insert a point, $x$, we draw $L\cdot k$ functions from $H$, denoted by
$(h_{i,j})_{i\in[L],j\in[k]}$. In each table $i\in[L]$ we insert $x$ into the
bucket keyed by $(h_{i,1}(x),h_{i,2}(x),\dots,h_{i,k}(x))$. Given a query
point $q$, the algorithm iterates over the $L$ tables and retrieves the data
points hashed into the same buckets as $q$. The process stops as soon as a
point is found within distance $r_{1}$ from $q$.

The algorithm as described has the performance characteristics listed below.
Here we assume the hash functions can be sampled and evaluated in constant
time. If this is not the case, one can use the improvements discussed in the
related work.

  * •

Query time: $O(L(k+np_{2}^{k}))=O(n^{\rho}p_{1}^{-1}\log n)$.

  * •

Space: $O(nL)=O(n^{1+\rho}p_{1}^{-1})$ plus the space to store the data
points.

  * •

Success probability $99\%$.

To get these bounds, we have defined $k=\lceil\frac{\log n}{\log
1/p_{2}}\rceil$ and

$\displaystyle L=\lceil p_{1}^{-k}\rceil\leq\exp\left(\log
1/p_{1}\cdot\lceil\tfrac{\log n}{\log 1/p_{2}}\rceil\right)+1\leq
n^{\rho}/p_{1}+1.$ (1)

It’s clear from this analysis that the $p_{1}^{-1}$ factor is only necessary
when $\frac{\log n}{\log 1/p_{2}}$ is not an integer. However in those cases
it is clearly necessary, since there is no obvious way to make a non-integer
number of function evaluations. We also cannot round $k$ down instead of up,
since the number of false positives would explode: rounding down would result
in a factor of $p_{2}^{-1}$ instead of $p_{1}^{-1}$ — much worse.

##  3 LSH with High-Low Tables

The idea of the main algorithm is to create some LSH tables with $k$ rounded
down, and some with $k$ rounded up. We call those respectively “high
probability” tables and “low probability” tables. In short “LSH with High-Low
Tables”.

The main theorem is the following:

###### Theorem 3.1

Let $H$ be a $(r_{1},r_{2},p_{1},p_{2})$-sensitive LSH family, and let
$\rho=\frac{\log 1/p_{1}}{\log 1/p_{2}}$. Assume $p_{1}>1/n$ and $p_{2}>1/n$.
Then there exists a solution to the $(r_{1},r_{2})$-near neighbour problem
with the following properties:

  * •

Query time: $O(n^{\rho}p_{1}^{\rho-1}\log n)$.

  * •

Space: $O(nL)=O(n^{1+\rho}p_{1}^{\rho-1})$ plus the space to store the data
points.

  * •

Success probability $99\%$.

###### Proof

Assume $r_{1},r_{2},p_{1},p_{2}$ are given. Define $\rho=\frac{\log
1/p_{1}}{\log 1/p_{2}}$, $\kappa=\frac{\log n}{\log 1/p_{2}}$, and
$\alpha=\lceil\kappa\rceil-\kappa\in[0,1)$. We build $\lfloor a\rfloor+\lceil
b\rceil$ tables (for $a,b\geq 0$ to be defined), where the first $\lfloor
a\rfloor$ use the hash function concatenated $\lfloor\kappa\rfloor$ times as
keys, and the remaining $\lceil b\rceil$ use it concatenated
$\lceil\kappa\rceil$ times.

The total number of tables to query is then $\lfloor a\rfloor+\lceil b\rceil$.
The expected total number of far points we have to retrieve is

$\displaystyle n(\lfloor a\rfloor p_{2}^{\lfloor\kappa\rfloor}+\lceil b\rceil
p_{2}^{\lceil\kappa\rceil})$ $\displaystyle=n(\lfloor a\rfloor
p_{2}^{\kappa-1+\alpha}+\lceil b\rceil p_{2}^{\kappa+\alpha})$ (2)
$\displaystyle=\lfloor a\rfloor p_{2}^{-1+\alpha}+\lceil b\rceil
p_{2}^{\alpha}$ (3) $\displaystyle\leq ap_{2}^{-1+\alpha}+(b+1)p_{2}^{\alpha}$
(4) $\displaystyle\leq ap_{2}^{-1+\alpha}+bp_{2}^{\alpha}+1.$ (5)

For the second equality, we used the definition of $\kappa$:
$p_{2}^{\kappa}=1/n$. We only count the expected number of points seen that
are at least $r_{2}$ away from the query. This is because the algorithm, like
classical LSH, terminates as soon as it sees a point with distance less than
$r_{2}$.

Given any point in the database within distance $r_{1}$ we must be able to
find it with high enough probability. This requires that the query and the
point shares a hash-bucket in one of the tables. The probability that this is
not the case is

$\displaystyle(1-p_{1}^{\lfloor\kappa\rfloor})^{\lfloor
a\rfloor}(1-p_{1}^{\lceil\kappa\rceil})^{\lceil b\rceil}$
$\displaystyle\leq(1-p_{1}^{\lfloor\kappa\rfloor})^{a-1}(1-p_{1}^{\lceil\kappa\rceil})^{b}$
(6)
$\displaystyle\leq\exp(-ap_{1}^{\lfloor\kappa\rfloor}-bp_{1}^{\lceil\kappa\rceil})(1-p_{1}^{\lfloor\kappa\rfloor})^{-1}$
(7)
$\displaystyle=\exp(-(ap_{1}^{-1+\alpha}+bp_{1}^{\alpha})n^{-\rho})(1-p_{1}^{\lfloor\kappa\rfloor})^{-1}$
(8)
$\displaystyle\leq\exp(-(ap_{1}^{-1+\alpha}+bp_{1}^{\alpha})n^{-\rho})\cdot
2.$ (9)

For the equality, we used the definition of $\kappa$ and $\rho$:
$p_{1}^{\kappa}=p_{2}^{\rho\kappa}=n^{-\rho}$. For the last inequality we have
assumed $p_{2}>1/n$ so $\lfloor\kappa\rfloor\geq 1$, and that $p_{1}<1/2$,
since otherwise we could just get the theorem from the classical LSH
algorithm.

We now define $a,b$ such that $ap_{2}^{-1+\alpha}+bp_{2}^{\alpha}=a+b$ and
$ap_{1}^{-1+\alpha}+bp_{1}^{\alpha}=n^{\rho}$. By the previous calculations
this will guarantee the number of false positives is not more than the number
of tables, and a constant success probability.

We can achieve this by taking

$\displaystyle\begin{bmatrix}a\\\
b\end{bmatrix}=\begin{bmatrix}p_{2}^{-1+\alpha}-1&p_{2}^{\alpha}-1\\\
p_{1}^{-1+\alpha}&p_{1}^{\alpha}\end{bmatrix}^{-1}\begin{bmatrix}0\\\
n^{\rho}\end{bmatrix}=\frac{n^{\rho}}{(p_{2}^{-1+\alpha}-1)p_{1}^{\alpha}+(1-p_{2}^{\alpha})p_{1}^{-1+\alpha}}\begin{bmatrix}1-p_{2}^{\alpha}\\\
p_{2}^{-1+\alpha}-1\end{bmatrix}.$ (10)

We note that both values are non-negative, since $\alpha\in[0,1]$.

When actually implementing the LSH algorithm width High-Low buckets, these are
the values you should use for the number of respectively the high and low
probability tables. That will ensure you take full advantage of when $\alpha$
is not worst case, and you may do even better than the theorem assumes.

To complete the theorem we need to prove $a+b\leq n^{\rho}p_{1}^{\rho-1}$. For
this we bound

$\displaystyle\frac{a+b}{n^{\rho}}$
$\displaystyle=\frac{p_{2}^{-1+\alpha}-p_{2}^{\alpha}}{(p_{2}^{-1+\alpha}-1)p_{1}^{\alpha}+(1-p_{2}^{\alpha})p_{1}^{-1+\alpha}}$
(11) $\displaystyle\leq\left(\frac{(p_{1}-p_{2})\log 1/p_{1}}{(1-p_{1})\log
p_{1}/p_{2}}\right)^{\rho}\left(\frac{(1-p_{2})\log
p_{1}/p_{2}}{(p_{1}-p_{2})\log 1/p_{2}}\right)$ (12)
$\displaystyle=\exp\left(D\left(\rho\,\bigg{\|}\,\frac{1/p_{1}-1}{1/p_{2}-1}\right)\right)$
(13) $\displaystyle\leq p_{1}^{\rho-1}.$ (14)

Here $D(r\|x)=r\log\tfrac{r}{x}+(1-r)\log\tfrac{1-r}{1-x}$ is the Kullback-
Leibler divergence. The two inequalities are proven in the Appendix as lemma 1
and lemma 3. The first bound comes from maximizing over $\alpha$, so in
principle we might be able to do better if $\kappa=\frac{\log n}{\log
1/p_{2}}$ is close to an integer. The second bound is harder, but the
realization that the left hand side can be written on the form of a divergence
helps a lot. The bound is tight up to a factor 2, so no significant
improvement is possible.

Finally we can boost the success probability from $1-2/e\approx 0.26$ to 99%
by repeating the entire data-structure 16 times.

## References

  * [1] Abboud, A., Rubinstein, A., Williams, R.: Distributed pcp theorems for hardness of approximation in p. In: 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS). pp. 25–36. IEEE (2017)
  * [2] Ahle, T.D., Aumüller, M., Pagh, R.: Parameter-free locality sensitive hashing for spherical range reporting. In: Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 239–256. SIAM (2017)
  * [3] Ahle, T.D.: Optimal las vegas locality sensitive data structures. In: 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS). pp. 938–949. IEEE (2017)
  * [4] Ahle, T.D., Knudsen, J.B.T.: Subsets and supermajorities: Optimal hashing-based set similarity search. arXiv preprint arXiv:1904.04045 (2020)
  * [5] Ahle, T.D., Pagh, R., Razenshteyn, I., Silvestri, F.: On the complexity of inner product similarity join. In: Proceedings of the 35th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems. pp. 151–164. ACM (2016)
  * [6] Andoni, A., Indyk, P., Laarhoven, T., Razenshteyn, I., Schmidt, L.: Practical and optimal lsh for angular distance. In: Advances in Neural Information Processing Systems. pp. 1225–1233 (2015)
  * [7] Andoni, A., Razenshteyn, I., Nosatzki, N.S.: Lsh forest: Practical algorithms made theoretical. In: Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 67–78. SIAM (2017)
  * [8] Bawa, M., Condie, T., Ganesan, P.: Lsh forest: self-tuning indexes for similarity search. In: Proceedings of the 14th international conference on World Wide Web. pp. 651–660 (2005)
  * [9] Becker, A., Ducas, L., Gama, N., Laarhoven, T.: New directions in nearest neighbor searching with applications to lattice sieving. In: Proceedings of the Twenty-Seventh Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 10–24. SIAM (2016)
  * [10] Broder, A.Z.: On the resemblance and containment of documents. In: Compression and Complexity of Sequences 1997. Proceedings. pp. 21–29. IEEE (1997)
  * [11] Broder, A.Z., Charikar, M., Frieze, A.M., Mitzenmacher, M.: Min-wise independent permutations. In: Proceedings of the thirtieth annual ACM symposium on Theory of computing. pp. 327–336. ACM (1998)
  * [12] Christiani, T.: Fast locality-sensitive hashing frameworks for approximate near neighbor search. In: International Conference on Similarity Search and Applications. pp. 3–17. Springer (2019)
  * [13] Christiani, T., Pagh, R.: Set similarity search beyond minhash. In: Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2017, Montreal, QC, Canada, June 19-23, 2017. pp. 1094–1107 (2017). https://doi.org/10.1145/3055399.3055443, https://doi.org/10.1145/3055399.3055443
  * [14] Christiani, T., Pagh, R., Thorup, M.: Confirmation sampling for exact nearest neighbor search. arXiv preprint arXiv:1812.02603 (2018)
  * [15] Christiani, T.L., Pagh, R., Aumüller, M., Vesterli, M.E.: Puffinn: Parameterless and universally fast finding of nearest neighbors. In: European Symposium on Algorithms. pp. 1–16 (2019)
  * [16] Dahlgaard, S., Knudsen, M.B.T., Thorup, M.: Fast similarity sketching. In: 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS). pp. 663–671. IEEE (2017)
  * [17] Datar, M., Immorlica, N., Indyk, P., Mirrokni, V.S.: Locality-sensitive hashing scheme based on p-stable distributions. In: Proceedings of the twentieth annual symposium on Computational geometry. pp. 253–262. ACM (2004)
  * [18] Indyk, P., Motwani, R.: Approximate nearest neighbors: towards removing the curse of dimensionality. In: Proceedings of the thirtieth annual ACM symposium on Theory of computing. pp. 604–613. ACM (1998)
  * [19] Leskovec, J., Rajaraman, A., Ullman, J.D.: Mining of massive data sets. Cambridge university press (2020)
  * [20] Mann, W., Augsten, N., Bouros, P.: An empirical evaluation of set similarity join techniques. Proceedings of the VLDB Endowment 9(9), 636–647 (2016)
  * [21] Razenshteyn, I., Schmidt, L.: Falconn-fast lookups of cosine and other nearest neighbors (2018)
  * [22] Wei, A.: Optimal las vegas approximate near neighbors in lp. In: Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 1794–1813. SIAM (2019)

##  4 Appendix

###### Lemma 1

For all $\alpha\in[0,1]$ we have

$\displaystyle
f(\alpha)=\frac{p_{2}^{-1+\alpha}-p_{2}^{\alpha}}{(p_{2}^{-1+\alpha}-1)p_{1}^{\alpha}+(1-p_{2}^{\alpha})p_{1}^{-1+\alpha}}\leq\left(\frac{(p_{1}-p_{2})\log
1/p_{1}}{(1-p_{1})\log p_{1}/p_{2}}\right)^{\rho}\left(\frac{(1-p_{2})\log
p_{1}/p_{2}}{(p_{1}-p_{2})\log 1/p_{2}}\right),$ (15)

where $\rho=\frac{\log 1/p_{1}}{\log 1/p_{2}}$.

###### Proof

We first show that $f(\alpha)$ is log-concave, which implies it is maximized
at the unique $\alpha$ such that $f^{\prime}(\alpha)=0$. Log-concavity follows
easily by noting

$\displaystyle\frac{d^{2}\log
f(\alpha)}{d\alpha^{2}}=-\frac{(1-p_{1})(p_{1}-p_{2})p_{2}^{1+\alpha}(\log\tfrac{1}{p_{2}})^{2}}{((1-p_{1})p_{2}+p_{2}^{\alpha}(p_{1}-p_{2}))^{2}}\leq
0.$ (16)

Meanwhile

$\displaystyle\frac{df(\alpha)}{d\alpha}=\frac{p_{1}(1-p_{2})(p_{2}/p_{1})^{\alpha}}{((1-p_{1})p_{2}+p_{2}^{\alpha}(p_{1}-p_{2}))^{2}}\left[(p_{1}-p_{2})p_{2}^{\alpha}\log\tfrac{1}{p_{1}}-(1-p_{1})p_{2}\log\tfrac{p_{2}}{p_{1}}\right],$
(17)

which implies $f(\alpha)$ is maximized in

$\displaystyle\alpha=\log\frac{(1-p_{1})p_{2}\log\tfrac{p_{2}}{p_{1}}}{(p_{1}-p_{2})\log\tfrac{1}{p_{1}}}\big{/}\log
p_{2}.$ (18)

Plugging into $f$ yields the lemma.

Note that $f(\alpha)$ is not regularly concave as $p_{1}$ and $p_{2}$ gets
small enough. Hence the use of log-concavity is necessary.

Next, we state a useful inequality, which is needed for the last proof.

###### Lemma 2

Let $p,r\in(0,1)$, then

$\displaystyle 1-\frac{1-p}{r}\leq p^{1/r}\leq\frac{pr}{1-p(1-r)}.$ (19)

###### Proof

We have $\frac{d^{2}}{dp^{2}}p^{1/r}=p^{1/r}(pr)^{-2}(1-r)$, so $p^{1/r}$ is
convex as a function of $p$. Since $1-\frac{1-p}{r}$ is it’s tangent (at
$p=1$) we get the first inequality.

For the second inequality, define $g(p)=p^{1/r}/\frac{pr}{1-p(1-r)}$. Then
$g(1)=1$ and $g(p)$ is non-decreasing, since

$\displaystyle g^{\prime}(p)=p^{1/r}(pr)^{-2}(1-p)(1-r)\geq 0.$ (20)

This shows that for $p\leq 1$ we have $p^{1/r}/\frac{pr}{1-p(1-r)}\leq 1$,
which is what we wanted to prove.

###### Lemma 3

Let $p,r\in(0,1]$ and let $x=\frac{1/p\,-\,1}{1/p^{1/r}-1}$, then

$\displaystyle D(r\|x)\leq r\log\tfrac{r}{x}\leq(1-r)\log\tfrac{1}{p},$ (21)

where $D(r\|x)=r\log\tfrac{r}{x}+(1-r)\log\tfrac{1-r}{1-x}$.

###### Proof

Using the upper bound of lemma 2 it follows directly that $x\in(0,r)$. This
suffices to show the first inequality of (21), since for $x\leq r$ we have
$\frac{1-r}{1-x}\leq 1$ and so the second term of $D(r\|x)$ is non-positive.

For the second inequality, we note that it is equivalent after manipulations
to $x\geq rp^{1/r-1}$. Plugging in $x$, and after more simple manipulations,
that’s equivalent in the range to $p^{1/r}\geq 1-\frac{1-p}{r}$, which is
lemma 2.

This finishes the proof of lemma 3.

It’s somewhat surprising that the last argument in the proof of lemma 3 works,
since if we had plugged the lower bound from lemma 2 directly into the problem
we would have had

$\displaystyle r\log\tfrac{r}{x}\leq r\log\tfrac{pr}{p+r-1},$ (22)

which is much weaker than what we prove, and not even defined for $p+r<1$.